CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Accessing cells selected by sphereToCell

Register Blogs Community New Posts Updated Threads Search

Like Tree5Likes
  • 2 Post By Tobermory
  • 1 Post By olesen
  • 1 Post By Tobermory
  • 1 Post By olesen

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   March 9, 2024, 10:38
Default Accessing cells selected by sphereToCell
  #1
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
How do i access the cells selected by the operation sphereToCell?



Code:
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
Divyaprakash is offline   Reply With Quote

Old   March 10, 2024, 08:25
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 746
Rep Power: 14
Tobermory will become famous soon enough
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:
Code:
 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 and Severus like this.
Tobermory is offline   Reply With Quote

Old   March 10, 2024, 10:25
Default
  #3
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
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?
Divyaprakash is offline   Reply With Quote

Old   March 10, 2024, 15:51
Default
  #4
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
Quote:
Originally Posted by Divyaprakash View Post
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:

Code:

// 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.
Divyaprakash likes this.
olesen is offline   Reply With Quote

Old   March 10, 2024, 16:00
Default
  #5
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
Quote:
Originally Posted by Divyaprakash View Post
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.
olesen is offline   Reply With Quote

Old   March 11, 2024, 03:47
Default
  #6
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
Quote:
Originally Posted by olesen View Post
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.

Code:
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.
Divyaprakash is offline   Reply With Quote

Old   March 11, 2024, 06:59
Default
  #7
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 746
Rep Power: 14
Tobermory will become famous soon enough
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 likes this.
Tobermory is offline   Reply With Quote

Old   March 11, 2024, 08:47
Default
  #8
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
Quote:
Originally Posted by Divyaprakash View Post
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.

Code:
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 is offline   Reply With Quote

Old   March 11, 2024, 08:51
Default
  #9
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
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 is offline   Reply With Quote

Old   March 11, 2024, 10:29
Default
  #10
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
Should note that indexedOctree also has a findSphere() method. The octree itself can be (demand driven) returned from the polyMesh.
olesen is offline   Reply With Quote

Old   March 12, 2024, 01:47
Default
  #11
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
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.
Divyaprakash is offline   Reply With Quote

Old   March 14, 2024, 16:29
Default
  #12
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
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.
olesen is offline   Reply With Quote

Old   March 16, 2024, 02:28
Default
  #13
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
Quote:
Originally Posted by Divyaprakash View Post
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).
Divyaprakash is offline   Reply With Quote

Old   March 16, 2024, 07:26
Default
  #14
Senior Member
 
Mark Olesen
Join Date: Mar 2009
Location: https://olesenm.github.io/
Posts: 1,715
Rep Power: 40
olesen has a spectacular aura aboutolesen has a spectacular aura about
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.
Divyaprakash likes this.
olesen is offline   Reply With Quote

Old   December 10, 2024, 03:07
Default selecting cells
  #15
Member
 
Divyaprakash
Join Date: Jun 2014
Posts: 74
Rep Power: 12
Divyaprakash is on a distinguished road
I figured out another way to select cells programmatically. It basically consists of using the octree object for the entire mesh and then using the findbox method.
Code:
     // Use OpenFOAM's octree for efficient spatial searching     
const indexedOctree<treeDataCell>& cellTree = mesh.cellTree();      
// Find all cells that intersect with our search box    
 labelHashSet foundCells;    
 label count = cellTree.findBox(searchBox, foundCells);
The details are in my blog post here.
https://dpcfd.com/posts/2024/12/blog-post-4/

Last edited by Divyaprakash; December 10, 2024 at 03:07. Reason: typo
Divyaprakash is offline   Reply With Quote

Reply

Tags
spheretocell, toposet


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[snappyHexMesh] snappyHexMesh does not work Siassei OpenFOAM Meshing & Mesh Conversion 8 September 28, 2023 05:01
[snappyHexMesh] SnappyHexMesh Problem 96faizizzuddin OpenFOAM Meshing & Mesh Conversion 10 October 10, 2022 07:45
Averaging over iterations for steady-state simulation CFD student Fluent UDF and Scheme Programming 8 September 22, 2022 04:39
[snappyHexMesh] SnappyHexMesh running killed! Mark JIN OpenFOAM Meshing & Mesh Conversion 7 June 14, 2022 02:37
problem with parallel calculations OF8 mpirun otaolafr OpenFOAM Programming & Development 0 January 5, 2021 08:33


All times are GMT -4. The time now is 14:13.