From c9c26c6086962c222a14bd51f3aab55a238a670d Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 31 Mar 2020 16:52:11 +0200 Subject: [PATCH] Added interpolate field on particles. --- include/deal.II/particles/utilities.h | 13 +++++ source/particles/utilities.cc | 79 +++++++++++++++++++++++++++ 2 files changed, 92 insertions(+) diff --git a/include/deal.II/particles/utilities.h b/include/deal.II/particles/utilities.h index e2e150f6c6..dd22963894 100644 --- a/include/deal.II/particles/utilities.h +++ b/include/deal.II/particles/utilities.h @@ -159,6 +159,19 @@ namespace Particles const AffineConstraints &constraints = AffineConstraints(), const ComponentMask &space_comps = ComponentMask()); + + + template + void + interpolate_field_on_particles( + const DoFHandler & field_dh, + const Particles::ParticleHandler &particle_handler, + const InputVectorType & field_vector, + OutputVectorType & interpolated_field, + const ComponentMask &field_comps = ComponentMask()); } // namespace Utilities } // namespace Particles DEAL_II_NAMESPACE_CLOSE diff --git a/source/particles/utilities.cc b/source/particles/utilities.cc index 2493aa9c91..b416d8b6bc 100644 --- a/source/particles/utilities.cc +++ b/source/particles/utilities.cc @@ -209,6 +209,85 @@ namespace Particles matrix.compress(VectorOperation::add); } + + + template + void + interpolate_field_on_particles( + const DoFHandler & field_dh, + const Particles::ParticleHandler &particle_handler, + const InputVectorType & field_vector, + OutputVectorType & interpolated_field, + const ComponentMask & field_comps) + { + if (particle_handler.n_locally_owned_particles() == 0) + { + interpolated_field.compress(VectorOperation::add); + return; // nothing else to do here + } + + const auto &tria = field_dh.get_triangulation(); + const auto &fe = field_dh.get_fe(); + auto particle = particle_handler.begin(); + const auto max_particles_per_cell = + particle_handler.n_global_max_particles_per_cell(); + + // Take care of components + const ComponentMask comps = + (field_comps.size() == 0 ? ComponentMask(fe.n_components(), true) : + field_comps); + AssertDimension(comps.size(), fe.n_components()); + const auto n_comps = comps.n_selected_components(); + + AssertDimension(field_vector.size(), field_dh.n_dofs()); + AssertDimension(interpolated_field.size(), + particle_handler.get_next_free_particle_index() * + n_comps); + // Add check on locally owned indices + + // Global to local indices + std::vector space_gtl(fe.n_components(), + numbers::invalid_unsigned_int); + for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i) + if (comps[i]) + space_gtl[i] = j++; + + std::vector dof_indices(fe.dofs_per_cell); + + while (particle != particle_handler.end()) + { + const auto &cell = particle->get_surrounding_cell(tria); + const auto &dh_cell = + typename DoFHandler::cell_iterator(*cell, &field_dh); + dh_cell->get_dof_indices(dof_indices); + const auto pic = particle_handler.particles_in_cell(cell); + const auto n_particles = particle_handler.n_particles_in_cell(cell); + + Assert(pic.begin() == particle, ExcInternalError()); + for (unsigned int i = 0; particle != pic.end(); ++particle, ++i) + { + const auto &reference_location = + particle->get_reference_location(); + + const auto id = particle->get_id(); + + for (unsigned int j = 0; j < fe.dofs_per_cell; ++j) + { + const auto comp_j = + space_gtl[fe.system_to_component_index(j).first]; + if (comp_j != numbers::invalid_unsigned_int) + interpolated_field[id * n_comps + comp_j] += + fe.shape_value(j, reference_location) * + field_vector(dof_indices[j]); + } + } + } + interpolated_field.compress(VectorOperation::add); + } + #include "utilities.inst" } // namespace Utilities -- 2.39.5