utilities.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2020 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_particles_utilities_h
#define dealii_particles_utilities_h
 
#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/grid/grid_tools_cache.h>
 
#include <deal.II/lac/affine_constraints.h>
#include <deal.II/lac/sparsity_pattern_base.h>
 
#include <deal.II/matrix_free/fe_point_evaluation.h>
 
#include <deal.II/particles/particle_handler.h>
 
DEAL_II_NAMESPACE_OPEN
 
 
 
namespace Particles
{
  namespace Utilities
  {
    template <int dim, int spacedim, typename number = double>
    void
    create_interpolation_sparsity_pattern(
      const DoFHandler<dim, spacedim>                 &space_dh,
      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
      SparsityPatternBase                             &sparsity,
      const AffineConstraints<number>                 &constraints =
        AffineConstraints<number>(),
      const ComponentMask &space_comps = {});
 
    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 = {});
 
    template <int dim,
              int spacedim,
              typename InputVectorType,
              typename OutputVectorType>
    void
    interpolate_field_on_particles(
      const DoFHandler<dim, spacedim>                 &field_dh,
      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
      const InputVectorType                           &field_vector,
      OutputVectorType                                &interpolated_field,
      const ComponentMask                             &field_comps = {})
    {
      if (particle_handler.n_locally_owned_particles() == 0)
        {
          interpolated_field.compress(VectorOperation::add);
          return; // nothing else to do here
        }
 
      const auto &fe       = field_dh.get_fe();
      auto        particle = particle_handler.begin();
 
      // Take care of components
      const ComponentMask comps =
        (field_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
                                   field_comps);
      AssertDimension(comps.size(), fe.n_components());
      const auto n_comps = comps.n_selected_components();
 
      AssertDimension(field_vector.size(), field_dh.n_dofs());
      AssertDimension(interpolated_field.size(),
                      particle_handler.get_next_free_particle_index() *
                        n_comps);
 
      // 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++;
 
      std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
 
      while (particle != particle_handler.end())
        {
          const auto &cell = particle->get_surrounding_cell();
          const auto &dh_cell =
            typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &field_dh);
          dh_cell->get_dof_indices(dof_indices);
          const auto pic = particle_handler.particles_in_cell(cell);
 
          Assert(pic.begin() == particle, ExcInternalError());
          for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
            {
              const Point<dim> reference_location =
                particle->get_reference_location();
 
              const auto id = particle->get_id();
 
              for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
                {
                  const auto comp_j =
                    space_gtl[fe.system_to_component_index(j).first];
                  if (comp_j != numbers::invalid_unsigned_int)
                    interpolated_field[id * n_comps + comp_j] +=
                      fe.shape_value(j, reference_location) *
                      field_vector(dof_indices[j]);
                }
            }
        }
      interpolated_field.compress(VectorOperation::add);
    }
 
 
    template <int n_components,
              int dim,
              int spacedim,
              typename InputVectorType,
              typename OutputVectorType>
    void
    interpolate_field_on_particles(
      const DoFHandler<dim, spacedim>                 &field_dh,
      const Particles::ParticleHandler<dim, spacedim> &particle_handler,
      const InputVectorType                           &field_vector,
      OutputVectorType                                &interpolated_field,
      const ComponentMask                             &field_comps = {},
      const Mapping<dim, spacedim>                    &mapping =
        (ReferenceCells::get_hypercube<dim>()
#ifndef _MSC_VER
           .template get_default_linear_mapping<spacedim>()
#else
           .ReferenceCell<dim>::get_default_linear_mapping<spacedim>()
#endif
           ))
    {
      if (particle_handler.n_locally_owned_particles() == 0)
        {
          // This is a collective operation that must be matched by the
          // compress() at the end of the function on all other processes. The
          // values are written with operator=, hence we use insert here.
          interpolated_field.compress(VectorOperation::insert);
          return; // nothing else to do here
        }
 
      const auto &fe       = field_dh.get_fe();
      auto        particle = particle_handler.begin();
 
      // Take care of components
      const ComponentMask comps =
        (field_comps.size() == 0 ? ComponentMask(n_components, true) :
                                   field_comps);
      AssertDimension(comps.size(), n_components);
      const auto n_selected_comps = comps.n_selected_components();
 
      AssertDimension(field_vector.size(), field_dh.n_dofs());
      AssertDimension(interpolated_field.size(),
                      particle_handler.get_next_free_particle_index() *
                        n_selected_comps);
 
      std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
 
      // Generate an evaluator that will be used to interpolate the fields at
      // the particle location.
      FEPointEvaluation<n_components, dim, spacedim> evaluator(mapping,
                                                               fe,
                                                               update_values);
      std::vector<Point<dim>>            particle_reference_locations;
      std::vector<types::particle_index> particle_indices;
      Vector<double>                     local_dof_values(fe.dofs_per_cell);
 
      while (particle != particle_handler.end())
        {
          const auto &cell = particle->get_surrounding_cell();
          const auto &dh_cell =
            typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &field_dh);
          dh_cell->get_dof_indices(dof_indices);
          dh_cell->get_dof_values(field_vector, local_dof_values);
 
          const auto pic = particle_handler.particles_in_cell(cell);
          Assert(pic.begin() == particle, ExcInternalError());
 
          // Gather the reference location and ids of all particles
          particle_reference_locations.clear();
          particle_indices.clear();
          for (const auto &p : pic)
            {
              particle_reference_locations.emplace_back(
                p.get_reference_location());
              particle_indices.emplace_back(p.get_id());
            }
 
          evaluator.reinit(cell, particle_reference_locations);
          evaluator.evaluate(make_array_view(local_dof_values),
                             EvaluationFlags::values);
          for (unsigned int particle_index = 0; particle != pic.end();
               ++particle, ++particle_index)
            {
              const types::particle_index global_particle_id =
                particle_indices[particle_index];
 
              if constexpr (n_components == 1)
                {
                  // The single component may be deselected by the mask, in
                  // which case there is nothing to write (and the output
                  // vector has no entry for this particle).
                  if (comps[0])
                    interpolated_field[global_particle_id] =
                      evaluator.get_value(particle_index);
                }
              else
                {
                  unsigned int j_comp = 0;
                  for (unsigned int j = 0; j < n_components; ++j)
                    {
                      if (comps[j])
                        {
                          interpolated_field[global_particle_id *
                                               n_selected_comps +
                                             j_comp] =
                            evaluator.get_value(particle_index)[j];
                          j_comp++;
                        }
                    }
                }
            }
        }
 
      // The entries are written using operator=, possibly for particles whose
      // global id is owned by another process. Ship those values to their
      // owner and finalize the vector. This is collective and matches the
      // compress() in the early-return branch above.
      interpolated_field.compress(VectorOperation::insert);
    }
 
  } // namespace Utilities
} // namespace Particles
DEAL_II_NAMESPACE_CLOSE
 
#endif