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+unsubscr...@googlegroups.com.
To view this discussion on the web visit
https://groups.google.com/d/msgid/dealii/e7ff5cd5-6866-414c-8836-9890ddfc2fd7o%40googlegroups.com.
<<attachment: step-3.zip>>
