]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Particle interpolation sparsity.
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 30 Jan 2020 17:00:14 +0000 (18:00 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 16 Apr 2020 12:18:00 +0000 (14:18 +0200)
18 files changed:
include/deal.II/particles/utilities.h [new file with mode: 0644]
source/particles/CMakeLists.txt
source/particles/utilities.cc [new file with mode: 0644]
source/particles/utilities.inst.in [new file with mode: 0644]
tests/particles/particle_interpolation_01.cc [new file with mode: 0644]
tests/particles/particle_interpolation_01.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/particles/particle_interpolation_01.with_p4est=true.mpirun=2.output [new file with mode: 0644]
tests/particles/particle_interpolation_02.cc [new file with mode: 0644]
tests/particles/particle_interpolation_02.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/particles/particle_interpolation_02.with_p4est=true.mpirun=2.output [new file with mode: 0644]
tests/particles/particle_interpolation_03.cc [new file with mode: 0644]
tests/particles/particle_interpolation_03.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/particles/particle_interpolation_03.with_p4est=true.mpirun=2.output [new file with mode: 0644]
tests/particles/particle_interpolation_04.cc [new file with mode: 0644]
tests/particles/particle_interpolation_04.with_p4est=true.mpirun=1.output [new file with mode: 0644]
tests/particles/particle_interpolation_05.cc [new file with mode: 0644]
tests/particles/particle_interpolation_05.with_p4est=true.mpirun=1.with_trilinos=true.output [new file with mode: 0644]
tests/particles/particle_interpolation_05.with_p4est=true.mpirun=2.output [new file with mode: 0644]

diff --git a/include/deal.II/particles/utilities.h b/include/deal.II/particles/utilities.h
new file mode 100644 (file)
index 0000000..da9acdf
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/component_mask.h>
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_tools_cache.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include <deal.II/particles/particle_handler.h>
+
+
+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 <int dim,
+              int spacedim,
+              typename SparsityType,
+              typename number = double>
+    void
+    create_interpolation_sparsity_pattern(
+      const DoFHandler<dim, spacedim> &                space_dh,
+      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
+      SparsityType &                                   sparsity,
+      const AffineConstraints<number> &                constraints =
+        AffineConstraints<number>(),
+      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 <int dim, int spacedim, typename MatrixType>
+    void
+    create_interpolation_matrix(
+      const DoFHandler<dim, spacedim> &                space_dh,
+      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
+      MatrixType &                                     matrix,
+      const AffineConstraints<typename MatrixType::value_type> &constraints =
+        AffineConstraints<typename MatrixType::value_type>(),
+      const ComponentMask &space_comps = ComponentMask());
+  } // namespace Utilities
+} // namespace Particles
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 4bfae710b08902936bcfffdc78f98882f1e9f0d2..71dbc91ec8a42a75877c70df76551c9af3c85b40 100644 (file)
 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 (file)
index 0000000..0f75650
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#include <deal.II/lac/generic_linear_algebra.h>
+
+#include <deal.II/particles/utilities.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace Particles
+{
+  namespace Utilities
+  {
+    template <int dim, int spacedim, typename SparsityType, typename number>
+    void
+    create_interpolation_sparsity_pattern(
+      const DoFHandler<dim, spacedim> &                space_dh,
+      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
+      SparsityType &                                   sparsity,
+      const AffineConstraints<number> &                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<unsigned int> 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<types::global_dof_index> dof_indices(fe.dofs_per_cell);
+      std::vector<types::particle_index>   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<dim, spacedim>::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 <int dim, int spacedim, typename MatrixType>
+    void
+    create_interpolation_matrix(
+      const DoFHandler<dim, spacedim> &                space_dh,
+      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
+      MatrixType &                                     matrix,
+      const AffineConstraints<typename MatrixType::value_type> &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<unsigned int> 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<types::global_dof_index> dof_indices(fe.dofs_per_cell);
+      std::vector<types::particle_index>   particle_indices(
+        max_particles_per_cell * n_comps);
+
+      FullMatrix<typename MatrixType::value_type> 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<dim, spacedim>::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 (file)
index 0000000..4132f14
--- /dev/null
@@ -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<deal_II_dimension,
+                                                        deal_II_space_dimension,
+                                                        Sparsity,
+                                                        scalar>(
+      const DoFHandler<deal_II_dimension, deal_II_space_dimension> &space_dh,
+      const Particles::ParticleHandler<deal_II_dimension,
+                                       deal_II_space_dimension>
+        &                              particle_handler,
+      Sparsity &                       sparsity,
+      const AffineConstraints<scalar> &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<deal_II_dimension,
+                                              deal_II_space_dimension,
+                                              Matrix>(
+      const DoFHandler<deal_II_dimension, deal_II_space_dimension> &space_dh,
+      const Particles::ParticleHandler<deal_II_dimension,
+                                       deal_II_space_dimension>
+        &                                          particle_handler,
+      Matrix &                                     matrix,
+      const AffineConstraints<Matrix::value_type> &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 (file)
index 0000000..8ac53cc
--- /dev/null
@@ -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 <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/utilities.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  const unsigned int my_mpi_id =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim, spacedim> space_tria(
+    MPI_COMM_WORLD);
+  MappingQ1<dim, spacedim> 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<dim, spacedim> particle_handler(space_tria,
+                                                             mapping);
+
+  const unsigned int estimated_n_particles_per_processor = 1;
+
+  Particles::Generators::probabilistic_locations<dim, spacedim>(
+    space_tria,
+    ConstantFunction<spacedim, double>(1.),
+    true,
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) *
+      estimated_n_particles_per_processor,
+    particle_handler);
+
+
+  FE_Q<dim, spacedim> space_fe(1);
+
+  ComponentMask space_mask(space_fe.n_components(), true);
+
+  const auto n_comps = space_mask.n_selected_components();
+
+  DoFHandler<dim, spacedim> 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<double> 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 (file)
index 0000000..89b8dea
--- /dev/null
@@ -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 (file)
index 0000000..7c1dc68
--- /dev/null
@@ -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 (file)
index 0000000..b991e8e
--- /dev/null
@@ -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 <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/utilities.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  const unsigned int my_mpi_id =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim, spacedim> space_tria(
+    MPI_COMM_WORLD);
+  MappingQ1<dim, spacedim> 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<dim, spacedim> particle_handler(space_tria,
+                                                             mapping);
+
+  const unsigned int estimated_n_particles_per_processor = 1;
+
+  Particles::Generators::probabilistic_locations<dim, spacedim>(
+    space_tria,
+    ConstantFunction<spacedim, double>(1.),
+    true,
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) *
+      estimated_n_particles_per_processor,
+    particle_handler);
+
+
+  FESystem<dim, spacedim> space_fe(FE_Q<dim, spacedim>(1), 3);
+
+  ComponentMask space_mask({false, false, true});
+
+  const auto n_comps = space_mask.n_selected_components();
+
+  DoFHandler<dim, spacedim> 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<double> 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 (file)
index 0000000..d7bfc3f
--- /dev/null
@@ -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 (file)
index 0000000..fe622ac
--- /dev/null
@@ -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 (file)
index 0000000..41b3308
--- /dev/null
@@ -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 <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/utilities.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  const unsigned int my_mpi_id =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim, spacedim> space_tria(
+    MPI_COMM_WORLD);
+  MappingQ1<dim, spacedim> 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<dim, spacedim> particle_handler(space_tria,
+                                                             mapping);
+
+  const unsigned int estimated_n_particles_per_processor = 1;
+
+  Particles::Generators::probabilistic_locations<dim, spacedim>(
+    space_tria,
+    ConstantFunction<spacedim, double>(1.),
+    true,
+    Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD) *
+      estimated_n_particles_per_processor,
+    particle_handler);
+
+
+  FESystem<dim, spacedim> space_fe(FE_Q<dim, spacedim>(1), 3);
+
+  ComponentMask space_mask(space_fe.n_components(), true);
+
+  const auto n_comps = space_mask.n_selected_components();
+
+  DoFHandler<dim, spacedim> 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<double> 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 (file)
index 0000000..5a611f7
--- /dev/null
@@ -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 (file)
index 0000000..f16787f
--- /dev/null
@@ -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 (file)
index 0000000..48159a8
--- /dev/null
@@ -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 <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/utilities.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  const unsigned int my_mpi_id =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim, spacedim> space_tria(
+    MPI_COMM_WORLD);
+  MappingQ1<dim, spacedim> 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<dim, spacedim> particle_handler(space_tria,
+                                                             mapping);
+
+  const unsigned int estimated_n_particles_per_processor = 1;
+
+
+  FE_Q<dim, spacedim> space_fe(1);
+
+  ComponentMask space_mask(space_fe.n_components(), true);
+
+  Particles::Generators::regular_reference_locations<dim, spacedim>(
+    space_tria, space_fe.get_unit_support_points(), particle_handler);
+
+  const auto n_comps = space_mask.n_selected_components();
+
+  DoFHandler<dim, spacedim> 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<double> 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<double> 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 (file)
index 0000000..fceb50e
--- /dev/null
@@ -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 (file)
index 0000000..ab8ecaf
--- /dev/null
@@ -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 <deal.II/base/function_lib.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_pattern.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/particles/generators.h>
+#include <deal.II/particles/particle_handler.h>
+#include <deal.II/particles/utilities.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim, int spacedim>
+void
+test()
+{
+  const unsigned int my_mpi_id =
+    Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+
+  parallel::distributed::Triangulation<dim, spacedim> space_tria(
+    MPI_COMM_WORLD);
+  MappingQ1<dim, spacedim> 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<dim, spacedim> particle_handler(space_tria,
+                                                             mapping);
+
+  FE_Q<dim, spacedim> space_fe(1);
+
+  ComponentMask space_mask(space_fe.n_components(), true);
+
+  Particles::Generators::regular_reference_locations<dim, spacedim>(
+    space_tria, QGauss<dim>(2).get_points(), particle_handler);
+
+  const auto n_comps = space_mask.n_selected_components();
+
+  DoFHandler<dim, spacedim> 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<double> 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<spacedim> linear(exponents, space_mask.size());
+
+  VectorTools::interpolate(space_dh, linear, field);
+  matrix.vmult(interpolation, field);
+
+  Vector<double> 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 (file)
index 0000000..6078ef1
--- /dev/null
@@ -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 (file)
index 0000000..0bc4b94
--- /dev/null
@@ -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
+
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.