[Divyaprakash]

How do i access the cells selected by the operation sphereToCell?

vector rr(0.0475,0.0475,0.005); // Define the point

Foam::sphereToCell sph(mesh,rr,0.005,0);
Info << "SPH: " << sph << endl;

I have added the above lines, but I am not sure how to access the cells that are selected. The reason I am doing this within the solver and not using a dict is that there are multiple points which vary with time in the simulation and I need to select cells around each of them to do some operation.



Similarly I would also like to use boxToCell. But it requires treeBoundBoxList as one of the arguments and I am not sure what that is

[Tobermory]

Take a look at the coding in sphereToCell (and then topoSetSource, and then topoSet and then cellSet) and you'll see that it is generating a topoSet containing the cell numbers that lie within the sphere, i.e. a hash list of the cell numbers. You can access the cells by interrogating the hash list ...

... or you could bypass it entirely and code up the cell selection yourself since you say you want to do it in the solver. Look at the following from sphereToCell.C:

 void Foam::sphereToCell::combine(topoSet& set, const bool add) const
 {
     const pointField& ctrs = mesh_.cellCentres();
     const scalar radSquared = radius_*radius_;
 
     forAll(ctrs, celli)
     {
         scalar offset = magSqr(centre_ - ctrs[celli]);
         if (offset <= radSquared)
         {
             addOrDelete(set, celli, add);
         }
     }
 }

You can see how it is identifying the cells in the sphere. Just replace the addOrDelete() command with whatever it is that you want to do inside the sphere. Hopefully that's enough to get you started.

[Divyaprakash]

Thank you Tobermory!


I understand it now. However I see that it loops over all the cells to find the relevant cells. I think it would be time consuming if I have multiple such points for which I need cells contained in a sphere centered at those points. I was under the impression that some mesh connectivity info might be used to generate such sets.



I have also been working on the problem for a while now and I am able to generate a list of such cells using neighbors of the cell containing my desired point and neighbors of those neighbors 2 times over. Do you think this is better or should I stick to looping over all the cells?

[olesen]

Originally Posted by Divyaprakash

Thank you Tobermory!


I understand it now. However I see that it loops over all the cells to find the relevant cells. I think it would be time consuming if I have multiple such points for which I need cells contained in a sphere centered at those points. I was under the impression that some mesh connectivity info might be used to generate such sets.



I have also been working on the problem for a while now and I am able to generate a list of such cells using neighbors of the cell containing my desired point and neighbors of those neighbors 2 times over. Do you think this is better or should I stick to looping over all the cells?

I think you'll need to benchmark it, but consider that the sphereToCell is using the magSqr of the distance - ie, it's a bounding sphere - it should be fairly quick. If you had something with mesh connectivity, I'm not sure how that would actually help.


For a micro optimization, could write the distance as:

// const scalar d2 = magSqr(ctrs[elemi] - origin_);
const scalar d2 = origin_.distSqr(ctrs[elemi]);

This would save one intermediate. Would need to verify with godbolt it that actually does something.

[olesen]

Originally Posted by Divyaprakash

Thank you Tobermory!


I understand it now. However I see that it loops over all the cells to find the relevant cells. I think it would be time consuming if I have multiple such points for which I need cells contained in a sphere centered at those points. I was under the impression that some mesh connectivity info might be used to generate such sets.


I have also been working on the problem for a while now and I am able to generate a list of such cells using neighbors of the cell containing my desired point and neighbors of those neighbors 2 times over. Do you think this is better or should I stick to looping over all the cells?

Could try with findCell() for your mesh (uses an octree search), from there walk out the neighbours and neighbour-neighbours until you've exceeded the search radius. Will need to pay really, really close attention at processor boundaries. You might have to walk out the sphere locally and detect any possible processor crossing. If there is something, a returnReduceOr(...) so that all ranks know there is still something to do. Next you'd need a syncTools boundary swap to communicate which faces (on the other side) should continue walking things out. Continue on this walking out until you have nothing more locally. For the first implementation, you'll probably want to forget about processors that are attached via edges or points - until you get the first part working.


Take a look at some of the cell/face/point wave algorithms (or regionSplit) for more ideas.

[Divyaprakash]

Originally Posted by olesen

Could try with findCell() for your mesh (uses an octree search), from there walk out the neighbours and neighbour-neighbours until you've exceeded the search radius. Will need to pay really, really close attention at processor boundaries. You might have to walk out the sphere locally and detect any possible processor crossing. If there is something, a returnReduceOr(...) so that all ranks know there is still something to do. Next you'd need a syncTools boundary swap to communicate which faces (on the other side) should continue walking things out. Continue on this walking out until you have nothing more locally. For the first implementation, you'll probably want to forget about processors that are attached via edges or points - until you get the first part working.


Take a look at some of the cell/face/point wave algorithms (or regionSplit) for more ideas.

Thank your reply. I will need to start looking more into the source code to understand about how the parallelization happens. Thanks for giving me pointers.
For now I have implemented the following to select cells around point contained within a cube of a size about 5cells^3.

label cellIndex = mesh.findCell(rr);

labelList appended;
appended.append(cellIndex);

// Get the neighbouring cells 4 times
for (int l = 0; l < 4; ++l) {
    int nn = appended.size();
    for (int m = 0; m < nn; ++m) {
        labelList neighbors = mesh.cellCells()[appended[m]];
        appended.append(neighbors);
    }
    // Get rid of repeating values and also sort the array
    Foam::inplaceUniqueSort(appended);
}

Can this be parallelized easily? Basically the list "appended" that I obtain at the end will be looped through to perform some calculations on those cells.

[Tobermory]

On the speed issue, if the cells that you are searching for do not change throughout the simulation, then you only need to do the search once (at the start of the run) and then either save the cell numbers in a list or populate a (constant) volScalarField with a window function etc etc. In that case, it doesn't matter if the searching process is rather basic.

Even if the cells you are searching for are changing each time step, then you need to compare the search time against the time spent solving the pressure equation ... again, I'd be suprised if it is significant.

[olesen]

Originally Posted by Divyaprakash

Thank your reply. I will need to start looking more into the source code to understand about how the parallelization happens. Thanks for giving me pointers.
For now I have implemented the following to select cells around point contained within a cube of a size about 5cells^3.

label cellIndex = mesh.findCell(rr);

labelList appended;
appended.append(cellIndex);

// Get the neighbouring cells 4 times
for (int l = 0; l < 4; ++l) {
    int nn = appended.size();
    for (int m = 0; m < nn; ++m) {
        labelList neighbors = mesh.cellCells()[appended[m]];
        appended.append(neighbors);
    }
    // Get rid of repeating values and also sort the array
    Foam::inplaceUniqueSort(appended);
}

Can this be parallelized easily? Basically the list "appended" that I obtain at the end will be looped through to perform some calculations on those cells.

Still local-only, but please use DynamicList instead of List, this will reduce the number of allocations. More modern to use push_back() instead of append() and then also with push_uniq() to avoid adding duplicates in the first place.
Next "safety issue", do not assume that findCell() returns a positive value. It could also return -1 (ie, not found so skip that case. It will probably be okay in serial, but
in parallel most ranks will not contain the point.

Can't give more details about parallel (typing from my phone...)

[olesen]

Some of this walking of neighbours reminds me of the queuing/stack in the Cuthill-McKee algorithm (meshTools bandCompression) - can steal some ideas from there.

[olesen]

Should note that indexedOctree also has a findSphere() method. The octree itself can be (demand driven) returned from the polyMesh.

[Divyaprakash]

Thank you so much olesen for your suggestions for improvement. It is very difficult for a beginner with no formal training in programming to even know things such as push_uniq() exist. I'll definelty read more about them. May I ask how can one look for these? Since the listOps.H doesn't have it in their functions list.

Also regarding parallelization, can you suggest some source to start reading about it's implementation.

[olesen]

There is no single easy entry point to parallel programming and structures since it will always really depend on what you are trying to do. Som basic MPI resources are always a good start.

[Divyaprakash]

Originally Posted by Divyaprakash

Thank you so much olesen for your suggestions for improvement. It is very difficult for a beginner with no formal training in programming to even know things such as push_uniq() exist. I'll definelty read more about them. May I ask how can one look for these? Since the listOps.H doesn't have it in their functions list.

Also regarding parallelization, can you suggest some source to start reading about it's implementation.

I went through the documentation and in ESI openFoam it exists under a different name in the source code appendUniq (Line no. 288).

[olesen]

No idea why the "latest" API links to v2212 instead of v2312, but push_back() and push_uniq() are definitely current.
The cool thing about naming the method push_back() instead of the old append() name is that lets you wrap the DynamicList with a std back_inserter and then pass into std copy_if and other algorithms.