utilities.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// 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.
//
// ---------------------------------------------------------------------

#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
{
  namespace Utilities
  {
    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());

    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());

    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 = ComponentMask())
    {
      if (particle_handler.n_locally_owned_particles() == 0)
        {
          interpolated_field.compress(VectorOperation::add);
          return; // nothing else to do here
        }

      const auto &tria     = field_dh.get_triangulation();
      const auto &fe       = field_dh.get_fe();
      auto        particle = particle_handler.begin();

      // 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.dofs_per_cell);

      while (particle != particle_handler.end())
        {
          const auto &cell = particle->get_surrounding_cell(tria);
          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 auto &reference_location =
                particle->get_reference_location();

              const auto id = particle->get_id();

              for (unsigned int j = 0; j < fe.dofs_per_cell; ++j)
                {
                  const auto comp_j =
                    space_gtl[fe.system_to_component_index(j).first];
                  if (comp_j != numbers::invalid_unsigned_int)
                    interpolated_field[id * n_comps + comp_j] +=
                      fe.shape_value(j, reference_location) *
                      field_vector(dof_indices[j]);
                }
            }
        }
      interpolated_field.compress(VectorOperation::add);
    }

  } // namespace Utilities
} // namespace Particles
DEAL_II_NAMESPACE_CLOSE

#endif