How to check a meshed medium by multiple spheres

[alisheikholeslam]

I have found this library today and it will help me a lot in my project. I have read the related article and thanks for the Python library. I have question that how could I get each mesh cells occupied volume by spheres (where is the documentation? Is it existed? Is it just the readme? If so, IMO, it needs to be improved) if we have a medium full of spheres.

1. Do we have to loop on each cell or sphere and find their corresponding spheres or cells, or we can do it with NumPy at once (not looping)? If we need loop, so we must order each cell nodes as it was explained (right hand rule outward)?
   I would be grateful if put a sample here which be representative to handle such medium, for hexahedron mesh cells (at least 2) and some spheres.
2. Also, I couldn't understand if the proposed method can handle smaller spheres, which are placed completely inside the cells. Thanks for clarifications.

[severinstrobl]

Thanks for your interest in using the library. You are correct, the documentation is somewhat minimal at this point, yet so is the interface. The original idea was not to obtain a solution being able to process a complete mesh and set of spheres, but to provide a building block that can easily be integrated into existing codes. Of course you can write a small wrapper which processes a number of cells/spheres in a format suitable for your application.

Regarding your questions:

1. The overlap library does not interface with any mesh format, as this is highly dependent on the application. So yes, you would have to loop over all the cells/spheres and perform a overlap calculation for each pairing and sum up the resulting overlap volumes accordingly. In case of larger numbers of cells or spheres, you might want to use some additional coarse collision detection scheme to reduce the computational cost. For a single hexahedron and sphere you can have a look at this test.
2. There is no limit on the size ratio between the cells and the spheres, at least to some degree, eventually numerical errors creep in (see Fig. 11 in the paper for details). I would say for a size ratio up to 100 it should be fine. The case of a (smaller) sphere being fully inside of a cell is handled as a special case to avoid most of the costly computations.

[alisheikholeslam]

Thanks for the answers.

I have doubt that two loops be efficient in terms of performance, because I need to do this in another loop. Also it needs to order the mesh in each loop :). I would be grateful if you would take time to help clarify my ambiguities:

1. Is the python library compatible with Numba library, which can accelerate the looping?
2. You have wrote a benchmark about 45000 spheres in the article. What does this sentence point to when say it done under 0.3 seconds. Is this showing the loop total time??

Assuming a system size of 45 000 particles, a figure common for coupled CFD–DEM simulations [8], this would lead to an average runtime of the exact volume calculation of less than 0.3 seconds on a single core of the CPU used for our benchmarks.

3. If we need to loop on all spheres, so what you mean by writing non-overlapping particles in the following sentence of the abstract:

The method can be used to … a system of spherical, non-overlapping particles

4. Is the model coordinate system need to be in a specified direction (I point it by ? on the image)?? I ask this because it was important in PFC software to be toward y direction.

5. Considering the shown coordinate (if have any effect) on the image, is the ordering true in this case (to be sure)?

6. My example will be something like the left pic, please let me know if there is any limitation to work on such a model (any recommendation).

7. Do you have any plan for developing this python library to handle the NumPy form for near future.

Many thanks for your time.
Yours sincerely

[severinstrobl]

@alisheikholeslam: For the moment, I converted the issue to a discussion, at least until we get a clearer idea of what is actually missing in the overlap library, I hope that's fine with you.

Let me try to answer the questions you raised in order:

1. I've never tried to use the library together with Numba. The Python bindings are a relatively new addition to the library, all the runs for the original publication used solely the C++ version. The Python layer is a very slim interface to the C++ implementation, but I have to admit I did not investigate the performance of the Python version so far. Could even be that the C++ code is not compiled with full optimizations for the Python package, I would have to look into that.
2. This ties into the previous answer: I was using the C++ version of the code and the numbers are for the case when the spheres and the elements are approximately the same size. For a tiny setup (1 sphere, 27 hexahedra), the average computation time was approximately 6 µs seconds. As the computation is not memory bound, this is extrapolated to a system with 45 000 spheres which would lead to the mentioned 0.3 seconds.
3. The non-overlapping refers to the spheres, e.g., a packing of hard spheres. There is no way to calculate the correct volume fraction using this library for each individual cell if the spheres themselves overlap. Based on the picture attached to your comment this does not seem to be the case, so you should be fine.
4. The only requirement is for the node ordering of the elements to match the one given by CGNS.
5. Yes, the node connectivity looks good to me.
6. As long as the spheres don't intersect, that looks like a perfect setup for using the library. I assume you are interested in the local volume fraction in each cell? Depending on the size of the system you might need to add an additional step to perform the coarse collision detection. This would limit the calls to the overlap library to the cells that actually touch/intersect the sphere. The specifics of this heavily depend on your application, which is why the overlap library does not address this part. For example, it could be that you already have for each cell a list of spheres intersecting this cell, then, to calculate the volume fraction, you would loop over all these spheres and subtract the sum from the volume of the cell. But if you don't have such a list, you would probably want to use a different library using octrees or similar techniques to build this list.
7. If I understand your question correctly, you are asking for some functionality to process entire lists of cells and spheres at once. I guess this is possible, but would require an additional library for the collision detection as discussed above. If you have experience in this area and can recommend to me a reliable Python package to perform this task, I'd be very interested. Implementing the functionality you are asking for should be relatively simple then.

[alisheikholeslam]

Thanks dear @severinstrobl for clarification.

I've never tried to use the library together with Numba.

just a little explanation for Numba library, this library make the loops very fast, and bring the python performance near the C one.

The only requirement is for the node ordering of the elements to match the one given by CGNS.

If the system coordinate be not important (as you confirmed), it will be enough for me.

Depending on the size of the system you might need to add an additional step to perform the coarse collision detection. If you have experience in this area and can recommend to me a reliable Python package to perform this task, I'd be very interested. Implementing the functionality you are asking for should be relatively simple then.

In this regard, to just devote a number of spheres in each loop, the fastest way that I know is to use cKDTree from SciPy library, but have some limitation based on the case. I have briefly explained the reasons in 1 p4-5, 2. I don't know if this module can search for nearest neighbors cubical, and because I have meshed the system by hexahedrons, so I wrote an algorithm the last day that find just overlapping spheres with each cubical elements and built such aforementioned lists. But, now I must check how the overlap library is written to try making this vectorized (of course just for hexahedrons that I need, in other cases, I think, cKDTree.query_ball_point from SciPy library and Balltree.query_radius from Sci-kit learn library can be very helpful). It will be a little time consuming for me due to my little knowledge in C++.

I am really enthusiastic about vectorizing this library.
Thanks for the answers
Best regards