Ciao Franco,
ParticleHandler was built with the intention of “losing” particles whenever
they fall out of the domain. This is the default behaviour of the
ParticleHandler, and unless you store a pointer to an existing particle
somewhere, you should not be in trouble.
If you iterate over the particles using the ParticleHandler, you should be
safe. The Utilities::interpolate_on_field cannot know how many particles are
actually there, so it asks the ParticleHandler for the
next_available_particle_index(), which should not change, even if you drop some
particles.
The way the interpolator works, is by looping *through* the ParticleHandler
(again, on the safe side), and deciding what to store where based on the
current particle id. Unless you renumber and reset all particle ids between the
removal and the interpolation, this will remain consistent.
Luca.
> On 16 Jul 2020, at 13:47, franco.m...@gmail.com <franco.milicc...@gmail.com>
> wrote:
>
> Thanks, Luca.
>
> I have read the documentation, that's why I was skeptic of my own code: I
> still think I am doing it wrong.
>
> As far as I understand from the documentation, particles should be removed
> from the handler after updating the cache, and thus the interpolated field
> should not contain the old removed one. Unfortunately I don't see any mention
> of "removal" in Particles::Utilities, but I could overlooking something here.
>
> Another question relative to the removal part is more of a worry. If I remove
> a particle, will any iterator be invalidated? If so, and say I need to remove
> a lot of particles, what is the best way to do this without too much
> overhead? As I understand, I have to remove (just one?) a particle and call a
> cache update at the end of the removal process.
>
> Thanks!
> Franco
> On Wednesday, July 15, 2020 at 7:47:58 PM UTC+2 luca....@gmail.com wrote:
> Franco,
>
> The interpolated field says that the field value is zero (on the line above
> your arrow). This is how it is documented: if a particle is removed, its
> interpolated value is left unchanged in the target vector. Zero in your case.
>
> Luca
>
>> Il giorno 15 lug 2020, alle ore 18:18, Franco Milicchio
>> <franco.m...@gmail.com> ha scritto:
>>
>>
>> Dear all,
>>
>> I am playing with 9.2 and particles, but I've encountered a weird problem.
>> If I remove a particle, the interpolated field on particle locations yields
>> a wrong answer.
>>
>> My code does (now) everything by hand:
>>
>> // Here be drangons...
>>
>> Point<2> p_0(0.1,0.5);
>> auto ref_cell_0 = GridTools::find_active_cell_around_point(dof_handler,
>> p_0);
>> Point<2> ref_p_0 =
>> StaticMappingQ1<2>::mapping.transform_real_to_unit_cell(ref_cell_0, p_0);
>>
>> Point<2> p_1(0.5,0.5);
>> auto ref_cell_1 = GridTools::find_active_cell_around_point(dof_handler,
>> p_1);
>> Point<2> ref_p_1 =
>> StaticMappingQ1<2>::mapping.transform_real_to_unit_cell(ref_cell_1, p_1);
>>
>> Point<2> p_2(0.9,0.5);
>> auto ref_cell_2 = GridTools::find_active_cell_around_point(dof_handler,
>> p_2);
>> Point<2> ref_p_2 =
>> StaticMappingQ1<2>::mapping.transform_real_to_unit_cell(ref_cell_2, p_2);
>>
>> particle_handler.insert_particle(Particles::Particle<2>(p_0, ref_p_0, 0),
>> ref_cell_0);
>> particle_handler.insert_particle(Particles::Particle<2>(p_1, ref_p_1, 1),
>> ref_cell_1);
>> particle_handler.insert_particle(Particles::Particle<2>(p_2, ref_p_2, 2),
>> ref_cell_2);
>>
>> particle_handler.update_cached_numbers();
>>
>> std::cout << "Printing particle id, location and reference location" <<
>> std::endl;
>> for(auto &p:particle_handler)
>> {
>> std::cout << "id: " << p.get_id() << " loc: " << p.get_location() << "
>> ref_loc: " << p.get_reference_location() << std::endl;
>> }
>>
>> field_on_particles.reinit(particle_handler.n_global_particles());
>> Particles::Utilities::interpolate_field_on_particles(dof_handler,
>> particle_handler, field, field_on_particles);
>>
>> std::cout << "Printing field value on particle location" << std::endl;
>> for(int pp = 0; pp< particle_handler.n_global_particles(); pp++)
>> {
>> std::cout << "fluid id: " << pp << " value: " <<
>> field_on_particles[pp] << std::endl;
>> }
>>
>> // ************************
>> // NOW REMOVE A PARTICLE...
>> // ************************
>>
>> for(auto itr=particle_handler.begin(); itr != particle_handler.end();
>> itr++)
>> {
>> if(itr->get_id()==1) particle_handler.remove_particle(itr);
>> }
>>
>> particle_handler.update_cached_numbers();
>>
>>
>> // ***********************
>> // ...AND VALUES ARE WRONG
>> // ***********************
>>
>> std::cout << "Printing particle id, location and reference location after
>> removing particle" << std::endl;
>>
>> for(auto itr_2=particle_handler.begin(); itr_2 != particle_handler.end();
>> itr_2++)
>> {
>> std::cout << "id: " << itr_2->get_id() << " loc: " <<
>> itr_2->get_location() << " ref_loc: " << itr_2->get_reference_location() <<
>> std::endl;
>> }
>>
>> field_on_particles.reinit(particle_handler.n_global_particles());
>> Particles::Utilities::interpolate_field_on_particles(dof_handler,
>> particle_handler, field, field_on_particles);
>>
>>
>> std::cout << "Printing field value on particles location after removing
>> particle" << std::endl;
>>
>> for(int ppp = 0; ppp < particle_handler.n_global_particles(); ppp++)
>> {
>> std::cout << "fluid id: " << ppp << " value: " <<
>> field_on_particles[ppp] << std::endl;
>> }
>>
>> The scalar field is generated with the X component, so its range is [0-1],
>> and as you can see before the removal the field is perfect. After, though, I
>> have weird numbers (forgive the ugly code and the useless reinit):
>>
>>
>> Number of active cells: 1024
>> Number of degrees of freedom: 1089
>> Printing particle id, location and reference location
>> id: 0 loc: 0.1 0.5 ref_loc: 0.2 1
>> id: 1 loc: 0.5 0.5 ref_loc: 1 1
>> id: 2 loc: 0.9 0.5 ref_loc: 0.8 1
>> Printing field value on particles location
>> fluid id: 0 value: 0.1 <--------------- OK!
>> fluid id: 1 value: 0.5 <--------------- OK!
>> fluid id: 2 value: 0.9 <--------------- OK!
>> Printing particle id, location and reference location after removing particle
>> id: 0 loc: 0.1 0.5 ref_loc: 0.2 1
>> id: 2 loc: 0.9 0.5 ref_loc: 0.8 1
>> Printing field value on particles location after removing particle
>> fluid id: 0 value: 0.1 <---------------------- STILL OK, I HOPE IT IS
>> fluid id: 1 value: 0. <---------------------- ????
>>
>> Am I using the particles the way they're not intended? Any hints are really
>> welcome!
>>
>> Thanks!
>> Franco
>>
>> PS. I am attaching a ZIP file with a minimal working code adapting the step
>> 3 of the tutorial.
>>
>>
>>
>> --
>> The deal.II project is located at http://www.dealii.org/
>> For mailing list/forum options, see
>> https://groups.google.com/d/forum/dealii?hl=en
>> ---
>> You received this message because you are subscribed to the Google Groups
>> "deal.II User Group" group.
>> To unsubscribe from this group and stop receiving emails from it, send an
>> email to dealii+un...@googlegroups.com.
>> To view this discussion on the web visit
>> https://groups.google.com/d/msgid/dealii/e7ff5cd5-6866-414c-8836-9890ddfc2fd7o%40googlegroups.com.
>> <step-3.zip>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see
> https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups
> "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to dealii+unsubscr...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/7c4738a4-aeae-4710-9a9d-35011aaedd85n%40googlegroups.com.
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to dealii+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/82DCB4F5-D1E9-4C4A-AF72-5D889F6032C2%40gmail.com.
