[ParticleMechanicsApplication] Accessing number of particles in background element

[alecontri]

Hi all, I was wondering if there was a way to access the number of particles falling inside a certain Background Grid element in the ParticleMechanicsApplication from the "updated_lagrangian.cpp" element. As I understand, every MPM element created has a particle with an associated background geometry in which it falls, but accessing the number of integration points inside the element geometry by construction it always gives 1 if I am not wrong. I am grateful to anyone who can help me.

[philbucher]

@KratosMultiphysics/mpm

[VeronikaSinger]

Hi @alecontri you can define the number of particles which are created/initialized within one element of the body mesh. The variable particles_per_element can be changed in the ParticleMaterials.json or in the GiD Interface in the definition of the body (it is not a variable of the background mesh.) Best Veronika

[alecontri]

Thank you very much @VeronikaSinger for replying. Ok I understand, but that holds only in the initial step, right? What about accessing the effective number of particles inside a background element during a dynamic simulation (I need it at every time step since I am modifying the updated lagrangian element for my purposes) where particles can cross element borders and move wildly?

[VeronikaSinger]

Yes true - this only holds for the initial step (sorry I probably missunderstood your question). Up to now, there is no possibility to check the number of particles within one element during the calculation process. To resolve this quickly one can probably add a counter variable on the element level and check the background element IDs on the python level. But if you need this feature not only for testing we should discuss possible implementations. @tteschemacher do you have a suggestion?

[alecontri]

I need it for experimenting a stabilization on the mixed element. I though about ways to do it, but since I am not really a very experienced programmer in C++, I do not know what was the correct way to handle it. If it can be useful I put the ideas here just in case:

• adding a counter to every background element (since the number of integration points cannot be used for it) which is reset before every searching done at the beginning of the time stepping and is augmented during the search every time the same ID is identified as containing a particle. Since it is made in parallel maybe a problem is that an atomic access can slow down too much the parallel phase. Another problem is that a change to the element in the Kratos core is needed.
• adding a counter to every particle (mMP in the ParticleApplication) which is modified after the search with the number of particles in the same element where itself is located. At this point I still need to memorize a vector which is incremented during the search so again atomic access or a vector for every particle that is then summed (omp reduction)

Maybe you can tell me what's the best way to proceed.

[VeronikaSinger]

I think we should keep it in the ParticleMechanicsApplication as this is very specific and add a counter to the particles. If you just need the number of particles within one background element you possibly could check the GeometryParentId in a python process.

[rubenzorrilla]

For your first option (which is my preferred option) you don't need to change the base element. You can achieve the same behavior as follows:

1. Create a new variable in the MPM application to specify the number of material points per background element, namely TOTAL_MATERIAL_POINTS. This would be the "counter".
2. Then, each time you move the material points points you update this value of this variable, that is to say, initialize it to zero and use a search algorithm (the BinBasedFastPointLocator would most probably be the best option for this), to locate each material point. Then you can do p_found_element->GetValue(TOTAL_MATERIAL_POINTS) += 1. This is something that needs an atomic clause as you pointed up.

Concerning the second option, this is in my opinion no go as makes the MPI parallelization way complex, aside of the fact that you'll need to save and manipulate a dense large vector.

[rubenzorrilla]

Closing as it seems to be solved. @alecontri feel free to reopen if needed.