From: Luca Heltai Date: Thu, 30 Jan 2020 17:00:14 +0000 (+0100) Subject: Particle interpolation sparsity. X-Git-Tag: v9.2.0-rc1~126^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0248119cf196c49e44433e537a42fd0eb31793e0;p=dealii.git Particle interpolation sparsity. --- diff --git a/include/deal.II/particles/utilities.h b/include/deal.II/particles/utilities.h new file mode 100644 index 0000000000..da9acdfb28 --- /dev/null +++ b/include/deal.II/particles/utilities.h @@ -0,0 +1,164 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii_particles_utilities +#define dealii_particles_utilities + +#include + +#include +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace Particles +{ + /** + * A namespace for functions offering tools to handle ParticleHandlers and + * their coupling with DoFHandlers + */ + namespace Utilities + { + /** + * Create an interpolation sparsity pattern for particles. + * + * Given a triangulation representing the domains $\Omega$, a particle + * handler of particles in $\Omega$, and a scalar finite element space + * $V(\Omega) = \text{span}\{v_j\}_{j=0}^n$, compute the sparsity pattern + * that would be necessary to assemble the matrix + * \f[ + * M_{ij} \dealcoloneq v_j(x_i) , + * \f] + * where $V(\Omega)$ is the finite element space associated with the + * `space_dh`, and the index `i` is given by the particle id whose position + * is `x_i`. + * + * In the case of vector valued finite element spaces, the components on + * which interpolation must be performed can be selected using a component + * mask. Only primitive finite element spaces are supported. + * + * When selecting more than one component, the resulting sparsity will have + * dimension equal to `particle_handler.n_global_particles() * + * mask.n_selected_components()` times `space_dh.n_dofs()`, and the + * corresponding matrix entries are given by + * \f[ + * M_{(i*n_comps+k),j} \dealcoloneq v_j(x_i) \cdot e_{comp_j}, + * \f] + * where `comp_j` is the only non zero component of the vector valued basis + * function `v_j` (equal to `fe.system_to_component_index(j).first`), and + * `k` corresponds to its index within the selected components of the mask. + * + * The `sparsity` is filled by locating the position of the particle with + * index `i` within the particle handler with respect to the embedding + * triangulation $\Omega$, and coupling it with all the local degrees of + * freedom specified in the component mask @p space_comps, following the + * ordering in which they are selected in the mask @p space_comps. + * + * If a particle does not fall within $\Omega$, it is ignored, and the + * corresponding rows of the sparsity will be empty. + * + * Constraints of the form supported by the AffineConstraints class may be + * supplied with the @p constraints argument. The method + * AffineConstraints::add_entries_local_to_global() is used to fill the + * final sparsity pattern. + * + * @author Bruno Blais, Luca Heltai, 2019 + */ + template + void + create_interpolation_sparsity_pattern( + const DoFHandler & space_dh, + const Particles::ParticleHandler &particle_handler, + SparsityType & sparsity, + const AffineConstraints & constraints = + AffineConstraints(), + const ComponentMask &space_comps = ComponentMask()); + + /** + * Create an interpolation matrix for particles. + * + * Given a triangulation representing the domains $\Omega$, a particle + * handler of particles in $\Omega$, and a scalar finite element space + * $V(\Omega) = \text{span}\{v_j\}_{j=0}^n$, compute the matrix + * \f[ + * M_{ij} \dealcoloneq v_j(x_i) , + * \f] + * where $V(\Omega)$ is the finite element space associated with the + * `space_dh`, and the index `i` is given by the particle id whose position + * is `x_i`. + * + * In the case of vector valued finite element spaces, the components on + * which interpolation must be performed can be selected using a component + * mask. Only primitive finite element spaces are supported. + * + * When selecting more than one component, the resulting sparsity will have + * dimension equal to `particle_handler.n_global_particles() * + * mask.n_selected_components()` times `space_dh.n_dofs()`, and the + * corresponding matrix entries are given by + * \f[ + * M_{(i*n_comps+k),j} \dealcoloneq v_j(x_i) \cdot e_{comp_j}, + * \f] + * where `comp_j` is the only non zero component of the vector valued basis + * function `v_j` (equal to `fe.system_to_component_index(j).first`), and + * `k` corresponds to its index within the selected components of the mask. + * + * The matrix is filled by locating the position of the particle with + * index `i` within the particle handler with respect to the embedding + * triangulation $\Omega$, and coupling it with all the local degrees of + * freedom specified in the component mask @p space_comps, following the + * ordering in which they are selected in the mask @p space_comps. + * + * If a particle does not fall within $\Omega$, it is ignored, and the + * corresponding rows of the matrix will be zero. + * + * Constraints of the form supported by the AffineConstraints class may be + * supplied with the @p constraints argument. The method + * AffineConstraints::distribute_local_to_global() is used to distribute + * the entries of the matrix to respect the given constraints. + * + * @author Bruno Blais, Luca Heltai, 2019 + */ + template + void + create_interpolation_matrix( + const DoFHandler & space_dh, + const Particles::ParticleHandler &particle_handler, + MatrixType & matrix, + const AffineConstraints &constraints = + AffineConstraints(), + const ComponentMask &space_comps = ComponentMask()); + } // namespace Utilities +} // namespace Particles +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/particles/CMakeLists.txt b/source/particles/CMakeLists.txt index 4bfae710b0..71dbc91ec8 100644 --- a/source/particles/CMakeLists.txt +++ b/source/particles/CMakeLists.txt @@ -17,22 +17,24 @@ INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) SET(_src + data_out.cc particle.cc particle_accessor.cc particle_iterator.cc particle_handler.cc generators.cc property_pool.cc - data_out.cc + utilities.cc ) SET(_inst + data_out.inst.in particle.inst.in particle_accessor.inst.in particle_iterator.inst.in particle_handler.inst.in generators.inst.in - data_out.inst.in + utilities.inst.in ) FILE(GLOB _header diff --git a/source/particles/utilities.cc b/source/particles/utilities.cc new file mode 100644 index 0000000000..0f75650b12 --- /dev/null +++ b/source/particles/utilities.cc @@ -0,0 +1,217 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + + +namespace Particles +{ + namespace Utilities + { + template + void + create_interpolation_sparsity_pattern( + const DoFHandler & space_dh, + const Particles::ParticleHandler &particle_handler, + SparsityType & sparsity, + const AffineConstraints & constraints, + const ComponentMask & space_comps) + { + if (particle_handler.n_locally_owned_particles() == 0) + return; // nothing to do here + + const auto &tria = space_dh.get_triangulation(); + const auto &fe = space_dh.get_fe(); + const auto max_particles_per_cell = + particle_handler.n_global_max_particles_per_cell(); + + + // Take care of components + const ComponentMask comps = + (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) : + space_comps); + AssertDimension(comps.size(), fe.n_components()); + + const auto n_comps = comps.n_selected_components(); + + // 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++; + + // [TODO]: when the add_entries_local_to_global below will implement + // the version with the dof_mask, this should be uncommented. + // // Construct a dof_mask, used to distribute entries to the sparsity + // Table<2, bool> dof_mask(max_particles_per_cell * n_comps, + // fe.dofs_per_cell); + // dof_mask.fill(false); + // for (unsigned int i = 0; i < space_fe.dofs_per_cell; ++i) + // { + // const auto comp_i = space_fe.system_to_component_index(i).first; + // if (space_gtl[comp_i] != numbers::invalid_unsigned_int) + // for (unsigned int j = 0; j < max_particles_per_cell; ++j) + // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true; + // } + + std::vector dof_indices(fe.dofs_per_cell); + std::vector particle_indices( + max_particles_per_cell * n_comps); + + auto particle = particle_handler.begin(); + while (particle != particle_handler.end()) + { + const auto &cell = particle->get_surrounding_cell(tria); + const auto dh_cell = + typename DoFHandler::cell_iterator(*cell, &space_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); + particle_indices.resize(n_particles * n_comps); + Assert(pic.begin() == particle, ExcInternalError()); + for (unsigned int i = 0; particle != pic.end(); ++particle, ++i) + { + const auto p_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) + constraints.add_entries_local_to_global( + {p_id * n_comps + comp_j}, {dof_indices[j]}, sparsity); + } + } + // [TODO]: when this works, use this: + // constraints.add_entries_local_to_global(particle_indices, + // dof_indices, + // sparsity, + // dof_mask); + } + } + + + + template + void + create_interpolation_matrix( + const DoFHandler & space_dh, + const Particles::ParticleHandler &particle_handler, + MatrixType & matrix, + const AffineConstraints &constraints, + const ComponentMask & space_comps) + { + if (particle_handler.n_locally_owned_particles() == 0) + { + matrix.compress(VectorOperation::add); + return; // nothing else to do here + } + + AssertDimension(matrix.n(), space_dh.n_dofs()); + + const auto &tria = space_dh.get_triangulation(); + const auto &fe = space_dh.get_fe(); + const auto max_particles_per_cell = + particle_handler.n_global_max_particles_per_cell(); + + // Take care of components + const ComponentMask comps = + (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) : + space_comps); + AssertDimension(comps.size(), fe.n_components()); + const auto n_comps = comps.n_selected_components(); + + AssertDimension(matrix.m(), + particle_handler.n_global_particles() * n_comps); + + + // 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++; + + // [TODO]: when the add_entries_local_to_global below will implement + // the version with the dof_mask, this should be uncommented. + // // Construct a dof_mask, used to distribute entries to the sparsity + // Table<2, bool> dof_mask(max_particles_per_cell * n_comps, + // fe.dofs_per_cell); + // dof_mask.fill(false); + // for (unsigned int i = 0; i < space_fe.dofs_per_cell; ++i) + // { + // const auto comp_i = space_fe.system_to_component_index(i).first; + // if (space_gtl[comp_i] != numbers::invalid_unsigned_int) + // for (unsigned int j = 0; j < max_particles_per_cell; ++j) + // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true; + // } + + std::vector dof_indices(fe.dofs_per_cell); + std::vector particle_indices( + max_particles_per_cell * n_comps); + + FullMatrix local_matrix( + max_particles_per_cell * n_comps, fe.dofs_per_cell); + + auto particle = particle_handler.begin(); + while (particle != particle_handler.end()) + { + const auto &cell = particle->get_surrounding_cell(tria); + const auto &dh_cell = + typename DoFHandler::cell_iterator(*cell, &space_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); + particle_indices.resize(n_particles * n_comps); + local_matrix.reinit({n_particles * n_comps, fe.dofs_per_cell}); + Assert(pic.begin() == particle, ExcInternalError()); + for (unsigned int i = 0; particle != pic.end(); ++particle, ++i) + { + const auto &reference_location = + particle->get_reference_location(); + + for (unsigned int d = 0; d < n_comps; ++d) + particle_indices[i * n_comps + d] = + particle->get_id() * n_comps + d; + + 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) + local_matrix(i * n_comps + comp_j, j) = + fe.shape_value(j, reference_location); + } + } + constraints.distribute_local_to_global(local_matrix, + particle_indices, + dof_indices, + matrix); + } + matrix.compress(VectorOperation::add); + } + +#include "utilities.inst" + + } // namespace Utilities +} // namespace Particles +DEAL_II_NAMESPACE_CLOSE diff --git a/source/particles/utilities.inst.in b/source/particles/utilities.inst.in new file mode 100644 index 0000000000..4132f14c1a --- /dev/null +++ b/source/particles/utilities.inst.in @@ -0,0 +1,51 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 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. +// +// --------------------------------------------------------------------- + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; + Sparsity : SPARSITY_PATTERNS; + scalar : REAL_SCALARS) + { +#if deal_II_dimension <= deal_II_space_dimension && deal_II_dimension > 1 + template void create_interpolation_sparsity_pattern( + const DoFHandler &space_dh, + const Particles::ParticleHandler + & particle_handler, + Sparsity & sparsity, + const AffineConstraints &constraints, + const ComponentMask & space_comps); +#endif + } + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS; + Matrix : SPARSE_MATRICES) + { +#if deal_II_dimension <= deal_II_space_dimension && deal_II_dimension > 1 + template void create_interpolation_matrix( + const DoFHandler &space_dh, + const Particles::ParticleHandler + & particle_handler, + Matrix & matrix, + const AffineConstraints &constraints, + const ComponentMask & space_comps); +#endif + } diff --git a/tests/particles/particle_interpolation_01.cc b/tests/particles/particle_interpolation_01.cc new file mode 100644 index 0000000000..8ac53ccb70 --- /dev/null +++ b/tests/particles/particle_interpolation_01.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// 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 an interpolation sparsity can be constructed for a single +// particle per processor for every valid dimension pair that exist + +#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); + + parallel::distributed::Triangulation space_tria( + MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(space_tria, -1, 1); + space_tria.refine_global(2); + + // Start by creating a particle handler that contains one particle per + // processor + Particles::ParticleHandler particle_handler(space_tria, + mapping); + + const unsigned int estimated_n_particles_per_processor = 1; + + Particles::Generators::probabilistic_locations( + space_tria, + ConstantFunction(1.), + true, + Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) * + estimated_n_particles_per_processor, + particle_handler); + + + FE_Q space_fe(1); + + ComponentMask space_mask(space_fe.n_components(), true); + + const auto n_comps = space_mask.n_selected_components(); + + DoFHandler space_dh(space_tria); + space_dh.distribute_dofs(space_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 + + << "Space FE: " << space_fe.get_name() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + DynamicSparsityPattern dsp(particle_handler.n_global_particles() * n_comps, + space_dh.n_dofs()); + + AffineConstraints constraints; + + // Build the interpolation sparsity + Particles::Utilities::create_interpolation_sparsity_pattern( + space_dh, particle_handler, dsp, constraints, space_mask); + + 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); + + SparsityTools::distribute_sparsity_pattern(dsp, + global_particles_index_set, + MPI_COMM_WORLD, + local_particle_index_set); + // Actually log the sparsity + { + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from(dsp); + sparsity_pattern.print_gnuplot(deallog.get_file_stream()); + } +} + +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/particle_interpolation_01.with_p4est=true.mpirun=1.output b/tests/particles/particle_interpolation_01.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..89b8dea85a --- /dev/null +++ b/tests/particles/particle_interpolation_01.with_p4est=true.mpirun=1.output @@ -0,0 +1,29 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 1 +DEAL:0:2d/2d::Space FE: FE_Q<2>(1) +DEAL:0:2d/2d::Space dofs: 25 +8 0 +13 0 +17 0 +21 0 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 1 +DEAL:0:2d/3d::Space FE: FE_Q<2,3>(1) +DEAL:0:2d/3d::Space dofs: 25 +8 0 +13 0 +17 0 +21 0 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 1 +DEAL:0:3d/3d::Space FE: FE_Q<3>(1) +DEAL:0:3d/3d::Space dofs: 125 +58 0 +59 0 +61 0 +62 0 +106 0 +107 0 +109 0 +110 0 diff --git a/tests/particles/particle_interpolation_01.with_p4est=true.mpirun=2.output b/tests/particles/particle_interpolation_01.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..7c1dc68330 --- /dev/null +++ b/tests/particles/particle_interpolation_01.with_p4est=true.mpirun=2.output @@ -0,0 +1,47 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 2 +DEAL:0:2d/2d::Space FE: FE_Q<2>(1) +DEAL:0:2d/2d::Space dofs: 25 +5 0 +8 0 +10 0 +13 0 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 2 +DEAL:0:2d/3d::Space FE: FE_Q<2,3>(1) +DEAL:0:2d/3d::Space dofs: 25 +5 0 +8 0 +10 0 +13 0 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 2 +DEAL:0:3d/3d::Space FE: FE_Q<3>(1) +DEAL:0:3d/3d::Space dofs: 125 +35 0 +36 0 +37 0 +38 0 +63 0 +64 0 +65 0 +66 0 + +7 -1 +8 -1 +16 -1 +17 -1 +7 -1 +8 -1 +16 -1 +17 -1 +78 -1 +80 -1 +82 -1 +83 -1 +87 -1 +89 -1 +91 -1 +92 -1 + diff --git a/tests/particles/particle_interpolation_02.cc b/tests/particles/particle_interpolation_02.cc new file mode 100644 index 0000000000..b991e8e804 --- /dev/null +++ b/tests/particles/particle_interpolation_02.cc @@ -0,0 +1,152 @@ +// --------------------------------------------------------------------- +// +// 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 an interpolation sparsity can be constructed for a single +// particle per processor for every valid dimension pair that exist +// when there are more than one components on the DoFHandler +// that is associated with the triangulation, and we select only one +// component + +#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); + + parallel::distributed::Triangulation space_tria( + MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(space_tria, -1, 1); + space_tria.refine_global(2); + + // Start by creating a particle handler that contains one particle per + // processor + Particles::ParticleHandler particle_handler(space_tria, + mapping); + + const unsigned int estimated_n_particles_per_processor = 1; + + Particles::Generators::probabilistic_locations( + space_tria, + ConstantFunction(1.), + true, + Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) * + estimated_n_particles_per_processor, + particle_handler); + + + FESystem space_fe(FE_Q(1), 3); + + ComponentMask space_mask({false, false, true}); + + const auto n_comps = space_mask.n_selected_components(); + + DoFHandler space_dh(space_tria); + space_dh.distribute_dofs(space_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 + + << "Space FE: " << space_fe.get_name() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + DynamicSparsityPattern dsp(particle_handler.n_global_particles() * n_comps, + space_dh.n_dofs()); + + AffineConstraints constraints; + + // Build the interpolation sparsity + Particles::Utilities::create_interpolation_sparsity_pattern( + space_dh, particle_handler, dsp, constraints, space_mask); + + 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); + + SparsityTools::distribute_sparsity_pattern(dsp, + global_particles_index_set, + MPI_COMM_WORLD, + local_particle_index_set); + // Actually log the sparsity + { + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from(dsp); + sparsity_pattern.print_gnuplot(deallog.get_file_stream()); + } +} + +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/particle_interpolation_02.with_p4est=true.mpirun=1.output b/tests/particles/particle_interpolation_02.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..d7bfc3f981 --- /dev/null +++ b/tests/particles/particle_interpolation_02.with_p4est=true.mpirun=1.output @@ -0,0 +1,29 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 1 +DEAL:0:2d/2d::Space FE: FESystem<2>[FE_Q<2>(1)^3] +DEAL:0:2d/2d::Space dofs: 75 +26 0 +41 0 +53 0 +65 0 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 1 +DEAL:0:2d/3d::Space FE: FESystem<2,3>[FE_Q<2,3>(1)^3] +DEAL:0:2d/3d::Space dofs: 75 +26 0 +41 0 +53 0 +65 0 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 1 +DEAL:0:3d/3d::Space FE: FESystem<3>[FE_Q<3>(1)^3] +DEAL:0:3d/3d::Space dofs: 375 +176 0 +179 0 +185 0 +188 0 +320 0 +323 0 +329 0 +332 0 diff --git a/tests/particles/particle_interpolation_02.with_p4est=true.mpirun=2.output b/tests/particles/particle_interpolation_02.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..fe622ac178 --- /dev/null +++ b/tests/particles/particle_interpolation_02.with_p4est=true.mpirun=2.output @@ -0,0 +1,47 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 2 +DEAL:0:2d/2d::Space FE: FESystem<2>[FE_Q<2>(1)^3] +DEAL:0:2d/2d::Space dofs: 75 +17 0 +26 0 +32 0 +41 0 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 2 +DEAL:0:2d/3d::Space FE: FESystem<2,3>[FE_Q<2,3>(1)^3] +DEAL:0:2d/3d::Space dofs: 75 +17 0 +26 0 +32 0 +41 0 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 2 +DEAL:0:3d/3d::Space FE: FESystem<3>[FE_Q<3>(1)^3] +DEAL:0:3d/3d::Space dofs: 375 +107 0 +110 0 +113 0 +116 0 +191 0 +194 0 +197 0 +200 0 + +23 -1 +26 -1 +50 -1 +53 -1 +23 -1 +26 -1 +50 -1 +53 -1 +236 -1 +242 -1 +248 -1 +251 -1 +263 -1 +269 -1 +275 -1 +278 -1 + diff --git a/tests/particles/particle_interpolation_03.cc b/tests/particles/particle_interpolation_03.cc new file mode 100644 index 0000000000..41b3308af9 --- /dev/null +++ b/tests/particles/particle_interpolation_03.cc @@ -0,0 +1,151 @@ +// --------------------------------------------------------------------- +// +// 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 an interpolation sparsity can be constructed for a single +// particle per processor for every valid dimension pair that exist +// when there are more than one components on the DoFHandler +// that is associated with the triangulation + +#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); + + parallel::distributed::Triangulation space_tria( + MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(space_tria, -1, 1); + space_tria.refine_global(2); + + // Start by creating a particle handler that contains one particle per + // processor + Particles::ParticleHandler particle_handler(space_tria, + mapping); + + const unsigned int estimated_n_particles_per_processor = 1; + + Particles::Generators::probabilistic_locations( + space_tria, + ConstantFunction(1.), + true, + Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) * + estimated_n_particles_per_processor, + particle_handler); + + + FESystem space_fe(FE_Q(1), 3); + + ComponentMask space_mask(space_fe.n_components(), true); + + const auto n_comps = space_mask.n_selected_components(); + + DoFHandler space_dh(space_tria); + space_dh.distribute_dofs(space_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 + + << "Space FE: " << space_fe.get_name() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + DynamicSparsityPattern dsp(particle_handler.n_global_particles() * n_comps, + space_dh.n_dofs()); + + AffineConstraints constraints; + + // Build the interpolation sparsity + Particles::Utilities::create_interpolation_sparsity_pattern( + space_dh, particle_handler, dsp, constraints, space_mask); + + 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); + + SparsityTools::distribute_sparsity_pattern(dsp, + global_particles_index_set, + MPI_COMM_WORLD, + local_particle_index_set); + // Actually log the sparsity + { + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from(dsp); + sparsity_pattern.print_gnuplot(deallog.get_file_stream()); + } +} + +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/particle_interpolation_03.with_p4est=true.mpirun=1.output b/tests/particles/particle_interpolation_03.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..5a611f756e --- /dev/null +++ b/tests/particles/particle_interpolation_03.with_p4est=true.mpirun=1.output @@ -0,0 +1,61 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 1 +DEAL:0:2d/2d::Space FE: FESystem<2>[FE_Q<2>(1)^3] +DEAL:0:2d/2d::Space dofs: 75 +24 0 +39 0 +51 0 +63 0 +25 -1 +40 -1 +52 -1 +64 -1 +26 -2 +41 -2 +53 -2 +65 -2 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 1 +DEAL:0:2d/3d::Space FE: FESystem<2,3>[FE_Q<2,3>(1)^3] +DEAL:0:2d/3d::Space dofs: 75 +24 0 +39 0 +51 0 +63 0 +25 -1 +40 -1 +52 -1 +64 -1 +26 -2 +41 -2 +53 -2 +65 -2 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 1 +DEAL:0:3d/3d::Space FE: FESystem<3>[FE_Q<3>(1)^3] +DEAL:0:3d/3d::Space dofs: 375 +174 0 +177 0 +183 0 +186 0 +318 0 +321 0 +327 0 +330 0 +175 -1 +178 -1 +184 -1 +187 -1 +319 -1 +322 -1 +328 -1 +331 -1 +176 -2 +179 -2 +185 -2 +188 -2 +320 -2 +323 -2 +329 -2 +332 -2 diff --git a/tests/particles/particle_interpolation_03.with_p4est=true.mpirun=2.output b/tests/particles/particle_interpolation_03.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..f16787fab7 --- /dev/null +++ b/tests/particles/particle_interpolation_03.with_p4est=true.mpirun=2.output @@ -0,0 +1,111 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 2 +DEAL:0:2d/2d::Space FE: FESystem<2>[FE_Q<2>(1)^3] +DEAL:0:2d/2d::Space dofs: 75 +15 0 +24 0 +30 0 +39 0 +16 -1 +25 -1 +31 -1 +40 -1 +17 -2 +26 -2 +32 -2 +41 -2 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 2 +DEAL:0:2d/3d::Space FE: FESystem<2,3>[FE_Q<2,3>(1)^3] +DEAL:0:2d/3d::Space dofs: 75 +15 0 +24 0 +30 0 +39 0 +16 -1 +25 -1 +31 -1 +40 -1 +17 -2 +26 -2 +32 -2 +41 -2 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 2 +DEAL:0:3d/3d::Space FE: FESystem<3>[FE_Q<3>(1)^3] +DEAL:0:3d/3d::Space dofs: 375 +105 0 +108 0 +111 0 +114 0 +189 0 +192 0 +195 0 +198 0 +106 -1 +109 -1 +112 -1 +115 -1 +190 -1 +193 -1 +196 -1 +199 -1 +107 -2 +110 -2 +113 -2 +116 -2 +191 -2 +194 -2 +197 -2 +200 -2 + +21 -3 +24 -3 +48 -3 +51 -3 +22 -4 +25 -4 +49 -4 +52 -4 +23 -5 +26 -5 +50 -5 +53 -5 +21 -3 +24 -3 +48 -3 +51 -3 +22 -4 +25 -4 +49 -4 +52 -4 +23 -5 +26 -5 +50 -5 +53 -5 +234 -3 +240 -3 +246 -3 +249 -3 +261 -3 +267 -3 +273 -3 +276 -3 +235 -4 +241 -4 +247 -4 +250 -4 +262 -4 +268 -4 +274 -4 +277 -4 +236 -5 +242 -5 +248 -5 +251 -5 +263 -5 +269 -5 +275 -5 +278 -5 + diff --git a/tests/particles/particle_interpolation_04.cc b/tests/particles/particle_interpolation_04.cc new file mode 100644 index 0000000000..48159a89f5 --- /dev/null +++ b/tests/particles/particle_interpolation_04.cc @@ -0,0 +1,153 @@ +// --------------------------------------------------------------------- +// +// 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 an interpolation matrix can be constructed for a single +// particle per vertex on a square matrix, and check we get an identity +// matrix + +#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); + + parallel::distributed::Triangulation space_tria( + MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(space_tria, -1, 1); + // space_tria.refine_global(1); + + // Start by creating a particle handler that contains one particle per + // processor + Particles::ParticleHandler particle_handler(space_tria, + mapping); + + const unsigned int estimated_n_particles_per_processor = 1; + + + FE_Q space_fe(1); + + ComponentMask space_mask(space_fe.n_components(), true); + + Particles::Generators::regular_reference_locations( + space_tria, space_fe.get_unit_support_points(), particle_handler); + + const auto n_comps = space_mask.n_selected_components(); + + DoFHandler space_dh(space_tria); + space_dh.distribute_dofs(space_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 + + << "Space FE: " << space_fe.get_name() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + DynamicSparsityPattern dsp(particle_handler.n_global_particles() * n_comps, + space_dh.n_dofs()); + + AffineConstraints constraints; + + // Build the interpolation sparsity + Particles::Utilities::create_interpolation_sparsity_pattern( + space_dh, particle_handler, dsp, constraints, space_mask); + + 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); + + SparsityTools::distribute_sparsity_pattern(dsp, + global_particles_index_set, + MPI_COMM_WORLD, + local_particle_index_set); + // Actually log the sparsity + + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from(dsp); + // Create a SparseMatrix + SparseMatrix matrix; + matrix.reinit(sparsity_pattern); + + // Build the interpolation matrix + Particles::Utilities::create_interpolation_matrix( + space_dh, particle_handler, matrix, constraints, space_mask); + + matrix.print(deallog.get_file_stream()); +} + +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/particle_interpolation_04.with_p4est=true.mpirun=1.output b/tests/particles/particle_interpolation_04.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..fceb50e5e8 --- /dev/null +++ b/tests/particles/particle_interpolation_04.with_p4est=true.mpirun=1.output @@ -0,0 +1,109 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 4 +DEAL:0:2d/2d::Space FE: FE_Q<2>(1) +DEAL:0:2d/2d::Space dofs: 4 +(0,0) 1.00000 +(0,1) 0.00000 +(0,2) 0.00000 +(0,3) 0.00000 +(1,1) 1.00000 +(1,0) 0.00000 +(1,2) 0.00000 +(1,3) 0.00000 +(2,2) 1.00000 +(2,0) 0.00000 +(2,1) 0.00000 +(2,3) 0.00000 +(3,3) 1.00000 +(3,0) 0.00000 +(3,1) 0.00000 +(3,2) 0.00000 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 4 +DEAL:0:2d/3d::Space FE: FE_Q<2,3>(1) +DEAL:0:2d/3d::Space dofs: 4 +(0,0) 1.00000 +(0,1) 0.00000 +(0,2) 0.00000 +(0,3) 0.00000 +(1,1) 1.00000 +(1,0) 0.00000 +(1,2) 0.00000 +(1,3) 0.00000 +(2,2) 1.00000 +(2,0) 0.00000 +(2,1) 0.00000 +(2,3) 0.00000 +(3,3) 1.00000 +(3,0) 0.00000 +(3,1) 0.00000 +(3,2) 0.00000 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 8 +DEAL:0:3d/3d::Space FE: FE_Q<3>(1) +DEAL:0:3d/3d::Space dofs: 8 +(0,0) 1.00000 +(0,1) 0.00000 +(0,2) 0.00000 +(0,3) 0.00000 +(0,4) 0.00000 +(0,5) 0.00000 +(0,6) 0.00000 +(0,7) 0.00000 +(1,1) 1.00000 +(1,0) 0.00000 +(1,2) 0.00000 +(1,3) 0.00000 +(1,4) 0.00000 +(1,5) 0.00000 +(1,6) 0.00000 +(1,7) 0.00000 +(2,2) 1.00000 +(2,0) 0.00000 +(2,1) 0.00000 +(2,3) 0.00000 +(2,4) 0.00000 +(2,5) 0.00000 +(2,6) 0.00000 +(2,7) 0.00000 +(3,3) 1.00000 +(3,0) 0.00000 +(3,1) 0.00000 +(3,2) 0.00000 +(3,4) 0.00000 +(3,5) 0.00000 +(3,6) 0.00000 +(3,7) 0.00000 +(4,4) 1.00000 +(4,0) 0.00000 +(4,1) 0.00000 +(4,2) 0.00000 +(4,3) 0.00000 +(4,5) 0.00000 +(4,6) 0.00000 +(4,7) 0.00000 +(5,5) 1.00000 +(5,0) 0.00000 +(5,1) 0.00000 +(5,2) 0.00000 +(5,3) 0.00000 +(5,4) 0.00000 +(5,6) 0.00000 +(5,7) 0.00000 +(6,6) 1.00000 +(6,0) 0.00000 +(6,1) 0.00000 +(6,2) 0.00000 +(6,3) 0.00000 +(6,4) 0.00000 +(6,5) 0.00000 +(6,7) 0.00000 +(7,7) 1.00000 +(7,0) 0.00000 +(7,1) 0.00000 +(7,2) 0.00000 +(7,3) 0.00000 +(7,4) 0.00000 +(7,5) 0.00000 +(7,6) 0.00000 diff --git a/tests/particles/particle_interpolation_05.cc b/tests/particles/particle_interpolation_05.cc new file mode 100644 index 0000000000..ab8ecaf802 --- /dev/null +++ b/tests/particles/particle_interpolation_05.cc @@ -0,0 +1,178 @@ +// --------------------------------------------------------------------- +// +// 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 an interpolation matrix can be constructed to interpolate +// Gauss quadrature points, and check the interpolation on a linear +// problem with Trilinos + +#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); + + parallel::distributed::Triangulation space_tria( + MPI_COMM_WORLD); + MappingQ1 mapping; + GridGenerator::hyper_cube(space_tria, -1, 1); + space_tria.refine_global(2); + + // Start by creating a particle handler that contains one particle per + // processor + Particles::ParticleHandler particle_handler(space_tria, + mapping); + + FE_Q space_fe(1); + + ComponentMask space_mask(space_fe.n_components(), true); + + Particles::Generators::regular_reference_locations( + space_tria, QGauss(2).get_points(), particle_handler); + + const auto n_comps = space_mask.n_selected_components(); + + DoFHandler space_dh(space_tria); + space_dh.distribute_dofs(space_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 + + << "Space FE: " << space_fe.get_name() << std::endl + << "Space dofs: " << space_dh.n_dofs() << std::endl; + + DynamicSparsityPattern dsp(particle_handler.n_global_particles() * n_comps, + space_dh.n_dofs()); + + AffineConstraints constraints; + + // Build the interpolation sparsity + Particles::Utilities::create_interpolation_sparsity_pattern( + space_dh, particle_handler, dsp, constraints, space_mask); + + 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); + + SparsityTools::distribute_sparsity_pattern(dsp, + global_particles_index_set, + MPI_COMM_WORLD, + local_particle_index_set); + // Create a SparseMatrix + TrilinosWrappers::SparseMatrix matrix; + matrix.reinit(local_particle_index_set, + locally_owned_dofs, + dsp, + MPI_COMM_WORLD); + + // Build the interpolation matrix + Particles::Utilities::create_interpolation_matrix( + space_dh, particle_handler, matrix, constraints, space_mask); + + // Create dofs and particle vectors + TrilinosWrappers::MPI::Vector field(locally_owned_dofs, MPI_COMM_WORLD); + TrilinosWrappers::MPI::Vector interpolation(local_particle_index_set, + MPI_COMM_WORLD); + + Tensor<1, spacedim> exponents; + for (unsigned int i = 0; i < spacedim; ++i) + exponents[i] = 1; + Functions::Monomial linear(exponents, space_mask.size()); + + VectorTools::interpolate(space_dh, linear, field); + matrix.vmult(interpolation, field); + + Vector values(space_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 (space_mask[i]) + if (std::abs(values[i] - interpolation(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/particle_interpolation_05.with_p4est=true.mpirun=1.with_trilinos=true.output b/tests/particles/particle_interpolation_05.with_p4est=true.mpirun=1.with_trilinos=true.output new file mode 100644 index 0000000000..6078ef1109 --- /dev/null +++ b/tests/particles/particle_interpolation_05.with_p4est=true.mpirun=1.with_trilinos=true.output @@ -0,0 +1,13 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 64 +DEAL:0:2d/2d::Space FE: FE_Q<2>(1) +DEAL:0:2d/2d::Space dofs: 25 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 64 +DEAL:0:2d/3d::Space FE: FE_Q<2,3>(1) +DEAL:0:2d/3d::Space dofs: 25 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 512 +DEAL:0:3d/3d::Space FE: FE_Q<3>(1) +DEAL:0:3d/3d::Space dofs: 125 diff --git a/tests/particles/particle_interpolation_05.with_p4est=true.mpirun=2.output b/tests/particles/particle_interpolation_05.with_p4est=true.mpirun=2.output new file mode 100644 index 0000000000..0bc4b941b3 --- /dev/null +++ b/tests/particles/particle_interpolation_05.with_p4est=true.mpirun=2.output @@ -0,0 +1,15 @@ + +DEAL:0:2d/2d::dim: 2, spacedim: 2 +DEAL:0:2d/2d::Total number of particles: 64 +DEAL:0:2d/2d::Space FE: FE_Q<2>(1) +DEAL:0:2d/2d::Space dofs: 25 +DEAL:0:2d/3d::dim: 2, spacedim: 3 +DEAL:0:2d/3d::Total number of particles: 64 +DEAL:0:2d/3d::Space FE: FE_Q<2,3>(1) +DEAL:0:2d/3d::Space dofs: 25 +DEAL:0:3d/3d::dim: 3, spacedim: 3 +DEAL:0:3d/3d::Total number of particles: 512 +DEAL:0:3d/3d::Space FE: FE_Q<3>(1) +DEAL:0:3d/3d::Space dofs: 125 + +