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 [email protected] 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 <
> [email protected]> 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 [email protected].
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/dealii/e7ff5cd5-6866-414c-8836-9890ddfc2fd7o%40googlegroups.com
>
> <https://groups.google.com/d/msgid/dealii/e7ff5cd5-6866-414c-8836-9890ddfc2fd7o%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
> <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 [email protected].
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/7c4738a4-aeae-4710-9a9d-35011aaedd85n%40googlegroups.com.
