From: blaisb Date: Tue, 14 Apr 2020 03:26:04 +0000 (-0400) Subject: - Added documentation for the function X-Git-Tag: v9.2.0-rc1~53^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6ed02de21f20e180d449648ba45c672393aec839;p=dealii.git - Added documentation for the function - Added a test for the interpolate on particle The only missing thing is instantiation --- diff --git a/include/deal.II/particles/utilities.h b/include/deal.II/particles/utilities.h index dd22963894..71ca8cfb4f 100644 --- a/include/deal.II/particles/utilities.h +++ b/include/deal.II/particles/utilities.h @@ -38,6 +38,7 @@ DEAL_II_NAMESPACE_OPEN + namespace Particles { /** @@ -161,6 +162,33 @@ namespace Particles const ComponentMask &space_comps = ComponentMask()); + + /** + * Interpolate a vector field to the particles of a particle_handler + * + * Given a DoFHandler and a particle handler, interpolates a vector field + * at the position of the particles. The result is stored in an output + * vector whose size corresponds to the number of locally owned particles * + * number of components + * + * @param[in] dof_handler The DOF Handler which was used to generate the + * field vector that is to be interpolated. + * + * @param[in] particle_handler The particle handler whose particle serve as + * the interpolation points. + * + * @param[in] field_vector The vector of the field to be interpolated. This + * vector must be coherent with the dof_handler provided + * + * @param[in,out] interpolated_field The interpolated value of the field at + * the position of the particles. The size of the vector must be + * n_locally_owned_particles * n_components + * + * @param[in] field_comps An optional component mask that decides which + * subset of the vector fields are interpolated + * + * @author Bruno Blais, Luca Heltai, 2019 + */ template &particle_handler, const InputVectorType & field_vector, OutputVectorType & interpolated_field, - const ComponentMask &field_comps = ComponentMask()); + const ComponentMask &field_comps = ComponentMask()) + { + 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); + } + } // namespace Utilities } // namespace Particles DEAL_II_NAMESPACE_CLOSE diff --git a/source/particles/utilities.cc b/source/particles/utilities.cc index b416d8b6bc..39d67b9559 100644 --- a/source/particles/utilities.cc +++ b/source/particles/utilities.cc @@ -211,82 +211,86 @@ namespace Particles - 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); + // 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); + // } - 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" diff --git a/tests/particles/interpolate_on_particle_01.cc b/tests/particles/interpolate_on_particle_01.cc new file mode 100644 index 0000000000..c5dfd36d81 --- /dev/null +++ b/tests/particles/interpolate_on_particle_01.cc @@ -0,0 +1,155 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Test that a Monomial function interpolated to a field +// can be interpolated at the particle position correctly + +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +using namespace dealii; + +template +void +test() +{ + const unsigned int my_mpi_id = + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + // Create a triangulation with a non-trivial number of degree of freedom + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(tria, -1, 1); + tria.refine_global(2); + + // Generate particles at the gauss points where the field will be interpolated + Particles::ParticleHandler particle_handler(tria, mapping); + + Particles::Generators::regular_reference_locations( + tria, QGauss(2).get_points(), particle_handler); + + + FE_Q fe(1); + + ComponentMask mask(fe.n_components(), true); + const auto n_comps = mask.n_selected_components(); + + DoFHandler space_dh(tria); + space_dh.distribute_dofs(fe); + + IndexSet locally_owned_dofs = space_dh.locally_owned_dofs(); + IndexSet locally_relevant_dofs; + DoFTools::extract_locally_relevant_dofs(space_dh, locally_relevant_dofs); + + if (my_mpi_id == 0) + deallog << "dim: " << dim << ", spacedim: " << spacedim << std::endl + << "Total number of particles: " + << particle_handler.n_global_particles() << std::endl; + + const auto n_local_particles_dofs = + particle_handler.n_locally_owned_particles() * n_comps; + + auto particle_sizes = + Utilities::MPI::all_gather(MPI_COMM_WORLD, n_local_particles_dofs); + + const auto my_start = std::accumulate(particle_sizes.begin(), + particle_sizes.begin() + my_mpi_id, + 0u); + + IndexSet local_particle_index_set(particle_handler.n_global_particles() * + n_comps); + + local_particle_index_set.add_range(my_start, + my_start + n_local_particles_dofs); + + auto global_particles_index_set = + Utilities::MPI::all_gather(MPI_COMM_WORLD, n_local_particles_dofs); + + // Create dofs and particle vectors + TrilinosWrappers::MPI::Vector field_owned(locally_owned_dofs, MPI_COMM_WORLD); + TrilinosWrappers::MPI::Vector field_relevant(locally_relevant_dofs, + MPI_COMM_WORLD); + TrilinosWrappers::MPI::Vector interpolation_on_particles( + local_particle_index_set, MPI_COMM_WORLD); + + // Create interpolation function + Tensor<1, spacedim> exponents; + for (unsigned int i = 0; i < spacedim; ++i) + exponents[i] = 1; + Functions::Monomial linear(exponents, mask.size()); + + // Interpolate function to vector than interpolate field to particles + VectorTools::interpolate(space_dh, linear, field_owned); + field_relevant = field_owned; + Particles::Utilities::interpolate_field_on_particles( + space_dh, particle_handler, field_relevant, interpolation_on_particles); + + Vector values(mask.size()); + + for (auto particle : particle_handler) + { + const auto &location = particle.get_location(); + const auto &id = particle.get_id(); + linear.vector_value(location, values); + for (unsigned int i = 0, j = 0; i < values.size(); ++i) + if (mask[i]) + if (std::abs(values[i] - interpolation_on_particles(id * n_comps + + (j++))) > 1e-10) + deallog << "NOT OK" << std::endl; + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll init; + + deallog.push("2d/2d"); + test<2, 2>(); + deallog.pop(); + deallog.push("2d/3d"); + test<2, 3>(); + deallog.pop(); + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/interpolate_on_particle_01.with_p4est=true.mpirun=2.with_trilinos=true.output b/tests/particles/interpolate_on_particle_01.with_p4est=true.mpirun=2.with_trilinos=true.output new file mode 100644 index 0000000000..0798a93813 --- /dev/null +++ b/tests/particles/interpolate_on_particle_01.with_p4est=true.mpirun=2.with_trilinos=true.output @@ -0,0 +1,9 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 64 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 64 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 512 + +