Reference documentation for deal.II version GIT relicensing-515-g7835781687 2024-05-02 18:40:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The step-89 tutorial program

This tutorial depends on step-67, step-87.

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

This program was contributed by Johannes Heinz, Maximilian Bergbauer, Marco Feder, and Peter Munch. Many ideas presented here are the result of common code development with Niklas Fehn, Luca Heltai, Martin Kronbichler, and Magdalena Schreter-Fleischhacker.

This tutorial is loosely based on the publication "High-order non-conforming discontinuous Galerkin methods for the acoustic conservation equations" by Johannes Heinz, Peter Munch, and Manfred Kaltenbacher [106].

Johannes Heinz was supported by the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowská–Curie Grant Agreement No. [812719].

Note
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation: 10.5281/zenodo.10033975

Introduction

This tutorial presents one way how to apply non-matching and/or Chimera methods within matrix-free loops in deal.II. We are following [106] to show that in some cases a simple point-to-point interpolation is not sufficient. As a remedy, Nitsche-type mortaring is used to suppress artificial modes observed for the acoustic conservation equations [106].

Acoustic conservation equations

Acoustic conservation equations are used to describe linear wave propagation. The set of equations consists of the conservation of mass and momentum

\[ \frac{\partial \, p}{\partial \, t} + \rho c^2 \nabla\cdot \mathbf{u} = 0,\\ \frac{\partial \, \mathbf{u}}{\partial \, t} + \frac{1}{\rho}\nabla p = \mathbf{0}. \]

Here, \(p\) is the acoustic pressure, \(\mathbf{u}\) the acoustic particle velocity, \(c\) the speed of sound, and \(\rho\) the mean density of the fluid in which waves are propagating. As stated above, the two equations are simply a different way of writing the wave equation: If you take the time derivative of the first equation, and the divergence of the second, i.e., compute

\[ \frac{\partial^2 \, p}{\partial \, t^2} + \rho c^2 \nabla\cdot \frac{\partial \mathbf{u}}{\partial t} = 0,\\ \frac{\partial \, \nabla \cdot \mathbf{u}}{\partial \, t} + \nabla \cdot \frac{1}{\rho}\nabla p = \mathbf{0}, \]

then you can substitute the second equation into the first one to obtain

\[ \frac{\partial^2 \, p}{\partial \, t^2} - \rho c^2 \nabla \cdot \frac{1}{\rho}\nabla p = \mathbf{0}, \]

which in the case of constant density \(\rho\) results in the more familiar form of the wave equation:

\[ \frac{\partial^2 \, p}{\partial \, t^2} - c^2 \Delta p = \mathbf{0}. \]

The reason one may want to consider the original form is that it has the form of a hyperbolic conservation law in which only first time and spatial derivatives appear. Whereas both the more familiar, second order form of the wave equation and the formulation as a first-order system conserve energy, it is often easier to devise numerical schemes that have the right amount of dissipation (necessary for numerical stability) using the well-developed machinery available for first-order systems.

For the discretization of this form, we make use of discontinuous Galerkin (DG) methods. DG methods are especially attractive for the acoustic conservation equations due to their low numerical dispersion errors. More importantly for this tutorial, DG methods natively extend to non-matching Nitsche-type methods [6]. I.e., numerical fluxes are not only used on inner element faces but also as non-matching coupling conditions.

The discretized equations read

\[ \int_{\Omega} q_h\frac{\partial \, p_h}{\partial \, t} +\int_{\Omega} q_h \rho c^2 \nabla\cdot\mathbf{u}_h +\int_{\partial\Omega} q_h\mathbf{n}\cdot\rho c^2(\mathbf{u}^*_h-\mathbf{u}_h)=0,\\ \int_{\Omega} \mathbf{w}_h\cdot\frac{\partial \,\mathbf{u}_h}{\partial \, t} +\int_{\Omega} \mathbf{w}_h\cdot \frac{1}{\rho} \nabla p_h +\int_{\partial\Omega} \mathbf{w}_h \cdot\mathbf{n} \frac{1}{\rho}(p^*_h-p_h)=\mathbf{0}, \]

where \(\mathbf{w}_h\) and \(q_h\) are test functions. The numerical fluxes are defined as follows [111] :

\[ p_h^*=p_h-\frac{\tau^-}{\tau^-+\tau^+}[p_h]+\frac{\tau^-\tau^+}{\tau^-+\tau^+}\jump{\mathbf{u}_h},\\ \mathbf{u}_h^*=\mathbf{u}_h-\frac{\gamma^-}{\gamma^-+\gamma^+}[\mathbf{u}_h]+\frac{\gamma^-\gamma^+}{\gamma^-+\gamma^+}\jump{p_h}, \]

with the penalty parameters \(\tau=\frac{\rho c}{2}\) and \(\gamma=\frac{1}{2\rho c}\). In these formulas, \([a] = a^- - a^+ \) denotes the jump of an arbitrary quantity \(a\) over element faces (face between elements \(K^-\) and \(K^+\)) and \(\jump{a} = a^- \mathbf{n}^- + a^+ \mathbf{n}^+\). For homogeneous materials, the fluxes reduce to standard Lax–Friedrichs fluxes ( \(\gamma^-=\gamma^+\) and \(\tau^-=\tau^+\))

\[ p_h^*=\average{p_h}+\tau\jump{\mathbf{u}_h},\\ \mathbf{u}_h^*=\average{\mathbf{u}_h}+\gamma\jump{p_h}. \]

The expression \(\average{a}=\frac{a^- + a^+}{2}\) denotes the averaging operator.

Non-matching discretizations

Non-matching discretizations can be used to connect mesh regions with different element sizes without the need for a transition region. Therefore, they are highly desirable in multiphysics applications. One example is a plate that radiates sound. The plate needs a much finer discretization than the surrounding air. In purely acoustic simulations, different materials require different element sizes to resolve the same wave because the speed of sound is directly proportional to the wavelength (we will simulate this example later on).

A different example of the usefulness of non-matching grids is where one wants to move the mesh in parts of the domain, but not others. A typical example is the simulation of windmills: One might want to enclose the rotating wings into a co-rotating mesh (to avoid having to remesh with every time step) but of course the mesh that describes the air above the surrounding landscape and around the tower on which the windmill is located should not rotate. In a case like this, one considers sliding rotating interfaces [72] between the co-rotating part of the mesh and the stationary part of the mesh, but this also requires the ability to handle non-matching discretizations.

Besides this, non-matching methods can be extended to Chimera methods (elements overlap). Chimera methods can help reduce the pressure on mesh generation tools since different regions of a mesh (that may overlap) can be meshed independently.

Different methods exist to treat non-matching interfaces. Within this tutorial, we concentrate on two methods: Point-to-point interpolation and Nitsche-type mortaring.

Point-to-point interpolation

Point-to-point interpolation is a naive approach. Whenever you need to compute integrals over the boundary of the cell at the left, for a coupled problem you then need to evaluate the solution or shape functions on the right at quadrature points of the face on the left, i.e., of the face of element \(K^-\). You can just evaluate these be interpolating the information on the right at these points, but this is in general expensive (read, for example, the documentation of VectorTools::point_value(), which implements this kind of functionality). As it can be seen from the picture this approach might be subject to aliasing in some cases.

Nitsche-type mortaring

Mortaring is the process of computing intersections and is not related to the Mortar method which enforces the coupling via Lagrange multipliers. Obtained intersections are also referred to as mortars. On each mortar a new integration rule is defined. The integral of the face of element \(K^-\) is computed on the intersections. The idea is that if we want to integrate something over a face \(f\subset \partial K^-\), that we break that integral into pieces:

\[ \int_f \cdots dx = \sum_i \int_{f_i} \cdots dx \]

where each of the \(f_i\) corresponds to the intersection of the original face \(f\) on the left with exactly one of the faces on the right; or, if we had multiple faces on the left, then the \(f_i\) would be the intersections of exactly one face on the left and one face on the right.

The point of this approach is first, that splitting the integral this way is exact. Secondly, and maybe more importantly, the terms we are integrating (the dots in the formula above) are now defined on one cell on each side, and consequently are smooth (whereas a finite element solution considered across multiple cells is, in general, not smooth). As a consequence, if we replace the integrals by the numerical integration (quadrature), then the result is exact as long as a sufficient number of integration points is used (at least for affine element shapes; for general curved elements, the integrand will contain rational expressions that are difficult to integrate exactly).

In this tutorial, the intersections are computed using CGAL, the Computational Geometry Algorithms Library. Therefore, deal.II has to be configured with DEAL_II_WITH_CGAL for the Nitsche-type mortaring implementation. See the deal.II Readme file for information about installation.

FERemoteEvaluation

In practice, for integrals as those mentioned above, we need to evaluate solutions (and shape functions) from cells across the non-matching interface. This is awkward enough if the other side is on the same processor, but outright difficult if the cells on the other side of the interface are owned by a different process in a parallel computation.

On regular meshes (say, doing things as we do in step-40), this is addressed by making sure that we are only computing integrals on locally owned cells and keeping around one layer of ghost cells for which we can query information. Ghost cells are the neighbors of locally owned cells, but in cases like the picture above, where the meshes are not matching, the cells on the other side of the interface are not neighbors in the logical sense – though they happen to be geometrically located adjacently. As a consequence, we need to find a way to efficiently query information on cells that are perhaps located on a different process.

FERemoteEvaluation is a wrapper class which provides a similar interface to, e.g., the FEEvaluation and FEFaceEvaluations classes to access values over non-matching interfaces in matrix-free loops. A detailed description on how to setup the class and how to use it in actual code is given below using hands-on examples. Within this tutorial we only show the usage for non-matching discretizations. Note however, that FERemoteEvaluation can also be used in other settings such as volume coupling. Under the hood, Utilities::MPI::RemotePointEvaluation is used to query the solution or gradients at quadrature points. A detailed description how this is done can be found in step-87. The main difference in the usage of FERemoteEvaluation compared to FEEvaluation is that the interpolated values/gradients of the finite element solution at all the quadrature points are precomputed globally before the loop over the cells/faces of the mesh (i.e., near the place where the communication takes place) instead of performing the interpolation out of the vector on a cell-by-cell basis.

The standard code to evaluate fluxes via FEEvaluation reads:

const auto face_function =
[&](const auto &data, auto &dst, const auto &src, const auto face_range) {
FEFaceEvaluation phi_m(data, true);
FEFaceEvaluation phi_p(data, false);
for (unsigned int f = face_range.first; f < face_range.second; ++f)
{
phi_m.reinit(f);
phi_p.reinit(f);
phi_p.gather_evaluate(src, EvaluationFlags::values); //compute values on face f
for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
phi_m.submit_value(phi_p.get_value(q), q); //access values with phi_p
phi_m.integrate_scatter(EvaluationFlags::values, dst);
}
};
matrix_free.template loop<VectorType, VectorType>(/* cell_operation= */{},
/* interior_face_operation= */ face_function,
/* boundary_face_operation= */{},
dst, src);

The code to evaluate fluxes via FERemoteEvaluation is shown below. For brevity, we assume all boundary faces are somehow connected via non-conforming interfaces for FERemoteEvaluation.

// Initialize FERemoteEvaluation: Note, that FERemoteEvaluation internally manages
// the memory to store precomputed values. Therefore, FERemoteEvaluation
// should be initialized only once to avoid continuous memory
// allocation/deallocation. At this point, remote_communicator is assumed
// to be initialized.
FERemoteEvaluation<dim,Number> phi_p_evaluator(remote_communicator);
// Precompute the interpolated values of the finite element solution at all
// the quadrature points outside the loop, invoking the vector entries and
// respective basis function at possibly remote MPI processes before communication.
phi_p_evaluator.gather_evaluate(src, EvaluationFlags::values);
const auto boundary_function =
[&](const auto &data, auto &dst, const auto &src, const auto face_range) {
FEFaceEvaluation phi_m(data, true);
// To access the values in a thread safe way each thread has
// to create a own accessor object. A small helper function
// provides the accessor.
auto phi_p = phi_p_evaluator.get_data_accessor();
for (unsigned int f = face_range.first; f < face_range.second; ++f)
{
phi_m.reinit(f);
phi_p.reinit(f);
for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
phi_m.submit_value(phi_p.get_value(q), q); //access values with phi_p
phi_m.integrate_scatter(EvaluationFlags::values, dst);
}
};
matrix_free.template loop<VectorType, VectorType>({}, {}, boundary_function, dst, src);

The object remote_communicator is of type FERemoteCommunicator and assumed to be correctly initialized in above code snippet. FERemoteCommunicator internally manages the update of ghost values over non-matching interfaces and keeps track of the mapping between quadrature point index and corresponding values/gradients. The update of the values/gradients happens before the actual matrix-free loop. FERemoteCommunicator, as well as FERemoteEvaluation behaves differently for the given template parameter value_type. If we want to access values at arbitrary points (e.g. in combination with FEPointEvaluation) value_type=Number. If the values are defined at quadrature points of a FEEvaluation object it is possible to get the values at the quadrature points of batches and value_type=VectorizedArray<Number>.

Overview

In the following, point-to-point interpolation and Nitsche-type mortaring is implemented.

At first we are considering the test case of a vibrating membrane, see e.g. [155]. Standing waves of length \(\lambda=2/M\) are oscillating with a period duration of \(T=2 / (M \sqrt{dim} c)\). \(M\) is the number of modes per meter, i.e. the number of half-waves per meter. The corresponding analytical solution reads as

\[ p =\cos(M \sqrt{d} \pi c t)\prod_{i=1}^{d} \sin(M \pi x_i),\\ u_i=-\frac{\sin(M \sqrt{d} \pi c t)}{\sqrt{d}\rho c} \cos(M \pi x_i)\prod_{j=1,j\neq i}^{d} \sin(M \pi x_j), \]

For simplicity, we are using homogeneous pressure Dirichlet boundary conditions within this tutorial. To be able to do so we have to tailor the domain size as well as the number of modes to conform with the homogeneous pressure Dirichlet boundary conditions. Within this tutorial we are using \(M=10\) and a domain that spans from \((0,0)\) to \((1,1)\).

For the point-to-point interpolation we observe aliasing which can be resolved using Nitsche-type mortaring.

In a more realistic example, we effectively apply the implementations to a test case in which a wave is propagating from one fluid into another fluid. The speed of sound in the left part of the domain the speed of sound is \(c=1\) and in the right part it is \(c=3\). Since the wavelength is directly proportional to the speed of sound, three times larger elements can be used in the right part of the domain to resolve waves up to the same frequency. The test case has been simulated with a different domain and different initial conditions, e.g. in [12].

The commented program

Include files

The program starts with including all the relevant header files.

  #include <deal.II/base/conditional_ostream.h>
  #include <deal.II/base/mpi.h>
 
  #include <deal.II/distributed/tria.h>
 
  #include <deal.II/fe/fe_dgq.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/mapping_q1.h>
 
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/grid_tools.h>
  #include <deal.II/grid/grid_tools_cache.h>
 
  #include <deal.II/matrix_free/fe_evaluation.h>
  #include <deal.II/matrix_free/matrix_free.h>
  #include <deal.II/matrix_free/operators.h>
 
  #include <deal.II/non_matching/mapping_info.h>
 
  #include <deal.II/numerics/data_out.h>
  #include <deal.II/numerics/vector_tools.h>

The following header file provides the class FERemoteEvaluation, which allows to access values and/or gradients at remote triangulations similar to FEEvaluation.

  #include <deal.II/matrix_free/fe_remote_evaluation.h>
 

We pack everything that is specific for this program into a namespace of its own.

  namespace Step89
  {
  using namespace dealii;
 

Initial conditions for vibrating membrane

Function that provides the initial condition for the vibrating membrane test case.

  template <int dim>
  class InitialConditionVibratingMembrane : public Function<dim>
  {
  public:
  InitialConditionVibratingMembrane(const double modes);
 

Function that the gives the initial pressure (comp 0) and velocity (comp 1 to 1 + dim).

  double value(const Point<dim> &p, const unsigned int comp) const final;
 
Definition point.h:111

Function that calculates the duration of one oscillation.

  double get_period_duration(const double speed_of_sound) const;
 
  private:
  const double M;
  };
 
  template <int dim>
  InitialConditionVibratingMembrane<dim>::InitialConditionVibratingMembrane(
  const double modes)
  : Function<dim>(dim + 1, 0.0)
  , M(modes)
  {
  static_assert(dim == 2, "Only implemented for dim==2");
  }
 
  template <int dim>
  double
  InitialConditionVibratingMembrane<dim>::value(const Point<dim> &p,
  const unsigned int comp) const
  {
  if (comp == 0)
  return std::sin(M * numbers::PI * p[0]) *
  std::sin(M * numbers::PI * p[1]);
 
  return 0.0;
  }
 
  template <int dim>
  double InitialConditionVibratingMembrane<dim>::get_period_duration(
  const double speed_of_sound) const
  {
  return 2.0 / (M * std::sqrt(dim) * speed_of_sound);
  }
 
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)

Gauss pulse

Function that provides the values of a pressure Gauss pulse.

  template <int dim>
  class GaussPulse : public Function<dim>
  {
  public:
  GaussPulse(const double shift_x, const double shift_y);
 

Function that the gives the initial pressure (comp 0) and velocity (comp 1 to 1 + dim).

  double value(const Point<dim> &p, const unsigned int comp) const final;
 
  private:
  const double shift_x;
  const double shift_y;
  };
 
  template <int dim>
  GaussPulse<dim>::GaussPulse(const double shift_x, const double shift_y)
  : Function<dim>(dim + 1, 0.0)
  , shift_x(shift_x)
  , shift_y(shift_y)
  {
  static_assert(dim == 2, "Only implemented for dim==2");
  }
 

Function that the gives the initial pressure (comp 0) and velocity (comp 1 to 1 + dim).

  template <int dim>
  double GaussPulse<dim>::value(const Point<dim> &p,
  const unsigned int comp) const
  {
  if (comp == 0)
  return std::exp(-1000.0 * ((std::pow(p[0] - shift_x, 2)) +
  (std::pow(p[1] - shift_y, 2))));
 
  return 0.0;
  }
 
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)

Helper functions

The following namespace contains free helper functions that are used in the tutorial.

  namespace HelperFunctions
  {

Helper function to check if a boundary ID is related to a non-matching face. A std::set that contains all non-matching boundary IDs is handed over additionally to the face ID under question. This function could certainly also be defined inline but this way the code is more easy to read.

  bool is_non_matching_face(
  const std::set<types::boundary_id> &non_matching_face_ids,
  const types::boundary_id face_id)
  {
  return non_matching_face_ids.find(face_id) != non_matching_face_ids.end();
  }
 

Helper function to set the initial conditions for the vibrating membrane test case.

  template <int dim, typename Number, typename VectorType>
  void set_initial_condition(MatrixFree<dim, Number> matrix_free,
  const Function<dim> &initial_solution,
  VectorType &dst)
  {
  VectorTools::interpolate(*matrix_free.get_mapping_info().mapping,
  matrix_free.get_dof_handler(),
  initial_solution,
  dst);
  }
 
const DoFHandler< dim > & get_dof_handler(const unsigned int dof_handler_index=0) const
const internal::MatrixFreeFunctions::MappingInfo< dim, Number, VectorizedArrayType > & get_mapping_info() const
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask={})

Helper function to compute the time step size according to the CFL condition.

  double
  compute_dt_cfl(const double hmin, const unsigned int degree, const double c)
  {
  return hmin / (std::pow(degree, 1.5) * c);
  }
 

Helper function that writes vtu output.

  template <typename VectorType, int dim>
  void write_vtu(const VectorType &solution,
  const DoFHandler<dim> &dof_handler,
  const Mapping<dim> &mapping,
  const unsigned int degree,
  const std::string &name_prefix)
  {
  DataOut<dim> data_out;
  flags.write_higher_order_cells = true;
  data_out.set_flags(flags);
 
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
  interpretation(
  std::vector<std::string> names(dim + 1, "U");
 
  names[0] = "P";
 
  data_out.add_data_vector(dof_handler, solution, names, interpretation);
 
  data_out.build_patches(mapping, degree, DataOut<dim>::curved_inner_cells);
  data_out.write_vtu_in_parallel(name_prefix + ".vtu",
  dof_handler.get_communicator());
  }
  } // namespace HelperFunctions
 
Abstract base class for mapping classes.
Definition mapping.h:318

Material access

This class stores the information if the fluid is homogeneous as well as the material properties at every cell. This class helps to access the correct values without accessing a large vector of materials in the homogeneous case.

  template <typename Number>
  class CellwiseMaterialData
  {
  public:
  template <int dim>
  CellwiseMaterialData(
  const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
  const std::map<types::material_id, std::pair<double, double>>
  &material_id_map)

If the map is of size 1, the material is constant in every cell.

  : homogeneous(material_id_map.size() == 1)
  {
  Assert(material_id_map.size() > 0,
  ExcMessage("No materials given to CellwiseMaterialData"));
 
  if (homogeneous)
  {
#define Assert(cond, exc)

In the homogeneous case we know the materials in the whole domain.

  speed_of_sound_homogeneous = material_id_map.begin()->second.first;
  density_homogeneous = material_id_map.begin()->second.second;
  }
  else
  {

In the in-homogeneous case materials vary between cells. We are filling a vector with the correct materials, that can be processed via read_cell_data().

  const auto n_cell_batches =
  matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
 
  speed_of_sound.resize(n_cell_batches);
  density.resize(n_cell_batches);
 
  for (unsigned int cell = 0; cell < n_cell_batches; ++cell)
  {
  speed_of_sound[cell] = 1.;
  density[cell] = 1.;
  for (unsigned int v = 0;
  v < matrix_free.n_active_entries_per_cell_batch(cell);
  ++v)
  {
  const auto material_id =
  matrix_free.get_cell_iterator(cell, v)->material_id();
 
  speed_of_sound[cell][v] =
  material_id_map.at(material_id).first;
  density[cell][v] = material_id_map.at(material_id).second;
  }
  }
  }
  }
 
  bool is_homogeneous() const
  {
  return homogeneous;
  }
 
  const AlignedVector<VectorizedArray<Number>> &get_speed_of_sound() const
  {
  Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()"));
  return speed_of_sound;
  }
 
  const AlignedVector<VectorizedArray<Number>> &get_density() const
  {
  Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()"));
  return density;
  }
 
  VectorizedArray<Number> get_homogeneous_speed_of_sound() const
  {
  Assert(homogeneous, ExcMessage("Use get_speed_of_sound()"));
  return speed_of_sound_homogeneous;
  }
 
  VectorizedArray<Number> get_homogeneous_density() const
  {
  Assert(homogeneous, ExcMessage("Use get_density()"));
  return density_homogeneous;
  }
 
  private:
  const bool homogeneous;
 

Materials in the in-homogeneous case.

Materials in the homogeneous case.

  VectorizedArray<Number> speed_of_sound_homogeneous;
  VectorizedArray<Number> density_homogeneous;
  };
 

To be able to access the material data in every cell in a thread safe way MaterialEvaluation is used. Similar to FEEvaluation, every thread creates its own instance and thus, there are no race conditions. For in-homogeneous materials, a reinit_cell() or reinit_face() function is used to set the correct material at the current cell batch. In the homogeneous case the _reinit() functions don't have to reset the materials.

  template <int dim, typename Number>
  class MaterialEvaluation
  {
  public:
  MaterialEvaluation(
  const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
  const CellwiseMaterialData<Number> &material_data)
  : phi(matrix_free)
  , phi_face(matrix_free, true)
  , material_data(material_data)
  {
  if (material_data.is_homogeneous())
  {

Set the material that is used in every cell.

  speed_of_sound = material_data.get_homogeneous_speed_of_sound();
  density = material_data.get_homogeneous_density();
  }
  }
 
  bool is_homogeneous() const
  {
  return material_data.is_homogeneous();
  }
 

Update the cell data, given a cell batch index.

  void reinit_cell(const unsigned int cell)
  {

In the homogeneous case we do not have to reset the cell data.

  if (!material_data.is_homogeneous())
  {

Reinit the FEEvaluation object and set the cell data.

  phi.reinit(cell);
  speed_of_sound =
  phi.read_cell_data(material_data.get_speed_of_sound());
  density = phi.read_cell_data(material_data.get_density());
  }
  }
 

Update the cell data, given a face batch index.

  void reinit_face(const unsigned int face)
  {

In the homogeneous case we do not have to reset the cell data.

  if (!material_data.is_homogeneous())
  {

Reinit the FEFaceEvaluation object and set the cell data.

  phi_face.reinit(face);
  speed_of_sound =
  phi_face.read_cell_data(material_data.get_speed_of_sound());
  density = phi_face.read_cell_data(material_data.get_density());
  }
  }
 

Return the speed of sound at the current cell batch.

  VectorizedArray<Number> get_speed_of_sound() const
  {
  return speed_of_sound;
  }
 

Return the density at the current cell batch.

  VectorizedArray<Number> get_density() const
  {
  return density;
  }
 
  private:

Members needed for the in-homogeneous case.

  FEEvaluation<dim, -1, 0, 1, Number> phi;
  FEFaceEvaluation<dim, -1, 0, 1, Number> phi_face;
 

Material defined at every cell.

  const CellwiseMaterialData<Number> &material_data;
 

Materials at current cell.

  VectorizedArray<Number> speed_of_sound;
  };
 
 

Boundary conditions

To be able to use the same kernel, for all face integrals we define a class that returns the needed values at boundaries. In this tutorial homogeneous pressure Dirichlet boundary conditions are applied via the mirror principle, i.e. \(p_h^+=-p_h^- + 2g\) with \(g=0\).

  template <int dim, typename Number>
  class BCEvaluationP
  {
  public:
  BCEvaluationP(const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m)
  : pressure_m(pressure_m)
  {}
 
  typename FEFaceEvaluation<dim, -1, 0, 1, Number>::value_type
  get_value(const unsigned int q) const
  {
  return -pressure_m.get_value(q);
  }
 
  private:
  const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m;
  };
 

We don't have to apply boundary conditions for the velocity, i.e. \(\mathbf{u}_h^+=\mathbf{u}_h^-\).

  template <int dim, typename Number>
  class BCEvaluationU
  {
  public:
  BCEvaluationU(const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m)
  : velocity_m(velocity_m)
  {}
 
  typename FEFaceEvaluation<dim, -1, 0, dim, Number>::value_type
  get_value(const unsigned int q) const
  {
  return velocity_m.get_value(q);
  }
 
  private:
  const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m;
  };
 

Acoustic operator

Class that defines the acoustic operator. The class is heavily based on matrix-free methods. For a better understanding in matrix-free methods please refer to step-67.

  template <int dim, typename Number, typename remote_value_type>
  class AcousticOperator
  {

If the remote evaluators are set up with a VectorizedArray we are using point-to-point interpolation. Otherwise we make use of Nitsche-type mortaring.

  static constexpr bool use_mortaring =
  std::is_floating_point_v<remote_value_type>;
 
  public:

In case of Nitsche-type mortaring, NonMatching::MappingInfo has to be provided in the constructor.

  AcousticOperator(
  const MatrixFree<dim, Number> &matrix_free,
  std::shared_ptr<CellwiseMaterialData<Number>> material_data,
  const std::set<types::boundary_id> &remote_face_ids,
  pressure_r_eval,
  velocity_r_eval,
  std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> c_r_eval,
  std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> rho_r_eval,
  std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> nm_info =
  nullptr)
  : matrix_free(matrix_free)
  , material_data(material_data)
  , remote_face_ids(remote_face_ids)
  , pressure_r_eval(pressure_r_eval)
  , velocity_r_eval(velocity_r_eval)
  , c_r_eval(c_r_eval)
  , rho_r_eval(rho_r_eval)
  , nm_mapping_info(nm_info)
  {
  if (use_mortaring)
  Assert(nm_info,
  ExcMessage(
  "In case of Nitsche-type mortaring NonMatching::MappingInfo \
  has to be provided."));
  }
 

Function to evaluate the acoustic operator.

  template <typename VectorType>
  void evaluate(VectorType &dst, const VectorType &src) const
  {

Update the precomputed values in corresponding the FERemoteEvaluation objects. The material parameters do not change and thus, we do not have to update precomputed values in c_r_eval and rho_r_eval.

  pressure_r_eval->gather_evaluate(src, EvaluationFlags::values);
  velocity_r_eval->gather_evaluate(src, EvaluationFlags::values);
 
  if constexpr (use_mortaring)
  {

Perform matrix free loop with Nitsche-type mortaring at non-matching faces.

  matrix_free.loop(
  &AcousticOperator::local_apply_cell<VectorType>,
  &AcousticOperator::local_apply_face<VectorType>,
  &AcousticOperator::local_apply_boundary_face_mortaring<VectorType>,
  this,
  dst,
  src,
  true,
  }
  else
  {

Perform matrix free loop with point-to-point interpolation at non-matching faces.

  matrix_free.loop(
  &AcousticOperator::local_apply_cell<VectorType>,
  &AcousticOperator::local_apply_face<VectorType>,
  &AcousticOperator::local_apply_boundary_face_point_to_point<
  VectorType>,
  this,
  dst,
  src,
  true,
  }
  }
 
  private:

This function evaluates the volume integrals.

  template <typename VectorType>
  void local_apply_cell(
  const MatrixFree<dim, Number> &matrix_free,
  VectorType &dst,
  const VectorType &src,
  const std::pair<unsigned int, unsigned int> &cell_range) const
  {
  FEEvaluation<dim, -1, 0, 1, Number> pressure(matrix_free, 0, 0, 0);
  FEEvaluation<dim, -1, 0, dim, Number> velocity(matrix_free, 0, 0, 1);
 

Class that gives access to the material at each cell

  MaterialEvaluation material(matrix_free, *material_data);
 
  for (unsigned int cell = cell_range.first; cell < cell_range.second;
  ++cell)
  {
  velocity.reinit(cell);
  pressure.reinit(cell);
 
  pressure.gather_evaluate(src, EvaluationFlags::gradients);
  velocity.gather_evaluate(src, EvaluationFlags::gradients);
 

Get the materials at the corresponding cell. Since we introduced MaterialEvaluation we can write the code independent if the material is homogeneous or in-homogeneous.

  material.reinit_cell(cell);
  const auto c = material.get_speed_of_sound();
  const auto rho = material.get_density();
  for (unsigned int q : pressure.quadrature_point_indices())
  {
  pressure.submit_value(rho * c * c * velocity.get_divergence(q),
  q);
  velocity.submit_value(1.0 / rho * pressure.get_gradient(q), q);
  }
 
  pressure.integrate_scatter(EvaluationFlags::values, dst);
  velocity.integrate_scatter(EvaluationFlags::values, dst);
  }
  }
 

This function evaluates the fluxes at faces between cells with the same material. If boundary faces are under consideration fluxes into neighboring faces do not have to be considered which is enforced via weight_neighbor=false. For non-matching faces the fluxes into neighboring faces are not considered as well. This is because we iterate over each side of the non-matching face separately (similar to a cell centric loop).

  template <bool weight_neighbor,
  typename InternalFaceIntegratorPressure,
  typename InternalFaceIntegratorVelocity,
  typename ExternalFaceIntegratorPressure,
  typename ExternalFaceIntegratorVelocity>
  inline DEAL_II_ALWAYS_INLINE void evaluate_face_kernel(
  InternalFaceIntegratorPressure &pressure_m,
  InternalFaceIntegratorVelocity &velocity_m,
  ExternalFaceIntegratorPressure &pressure_p,
  ExternalFaceIntegratorVelocity &velocity_p,
  const typename InternalFaceIntegratorPressure::value_type c,
  const typename InternalFaceIntegratorPressure::value_type rho) const
  {
#define DEAL_II_ALWAYS_INLINE
Definition config.h:109

Compute penalty parameters from material parameters.

  const auto tau = 0.5 * rho * c;
  const auto gamma = 0.5 / (rho * c);
 
  for (unsigned int q : pressure_m.quadrature_point_indices())
  {
  const auto n = pressure_m.normal_vector(q);
  const auto pm = pressure_m.get_value(q);
  const auto um = velocity_m.get_value(q);
 
  const auto pp = pressure_p.get_value(q);
  const auto up = velocity_p.get_value(q);
 

Compute homogeneous local Lax-Friedrichs fluxes and submit the corrsponding values to the integrators.

  const auto momentum_flux =
  0.5 * (pm + pp) + 0.5 * tau * (um - up) * n;
  velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
  if constexpr (weight_neighbor)
  velocity_p.submit_value(1.0 / rho * (momentum_flux - pp) * (-n), q);
 
  const auto mass_flux = 0.5 * (um + up) + 0.5 * gamma * (pm - pp) * n;
  pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
  if constexpr (weight_neighbor)
  pressure_p.submit_value(rho * c * c * (mass_flux - up) * (-n), q);
  }
  }
 

This function evaluates the fluxes at faces between cells with different materials. This can only happen over non-matching interfaces. Therefore, it is implicitly known that weight_neighbor=false and we can omit the parameter.

  template <typename InternalFaceIntegratorPressure,
  typename InternalFaceIntegratorVelocity,
  typename ExternalFaceIntegratorPressure,
  typename ExternalFaceIntegratorVelocity,
  typename MaterialIntegrator>
  void evaluate_face_kernel_inhomogeneous(
  InternalFaceIntegratorPressure &pressure_m,
  InternalFaceIntegratorVelocity &velocity_m,
  const ExternalFaceIntegratorPressure &pressure_p,
  const ExternalFaceIntegratorVelocity &velocity_p,
  const typename InternalFaceIntegratorPressure::value_type c,
  const typename InternalFaceIntegratorPressure::value_type rho,
  const MaterialIntegrator &c_r,
  const MaterialIntegrator &rho_r) const
  {

Interior material information is constant over quadrature points

  const auto tau_m = 0.5 * rho * c;
  const auto gamma_m = 0.5 / (rho * c);
 
  for (unsigned int q : pressure_m.quadrature_point_indices())
  {

The material at the neighboring face might vary in every quadrature point.

  const auto c_p = c_r.get_value(q);
  const auto rho_p = rho_r.get_value(q);
  const auto tau_p = 0.5 * rho_p * c_p;
  const auto gamma_p = 0.5 / (rho_p * c_p);
  const auto tau_sum_inv = 1.0 / (tau_m + tau_p);
  const auto gamma_sum_inv = 1.0 / (gamma_m + gamma_p);
 
  const auto n = pressure_m.normal_vector(q);
  const auto pm = pressure_m.get_value(q);
  const auto um = velocity_m.get_value(q);
 
  const auto pp = pressure_p.get_value(q);
  const auto up = velocity_p.get_value(q);
 
 

Compute inhomogeneous fluxes and submit the corresponding values to the integrators.

  const auto momentum_flux =
  pm - tau_m * tau_sum_inv * (pm - pp) +
  tau_m * tau_p * tau_sum_inv * (um - up) * n;
  velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
 
 
  const auto mass_flux =
  um - gamma_m * gamma_sum_inv * (um - up) +
  gamma_m * gamma_p * gamma_sum_inv * (pm - pp) * n;
 
  pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
  }
  }
 

This function evaluates the inner face integrals.

  template <typename VectorType>
  void local_apply_face(
  const MatrixFree<dim, Number> &matrix_free,
  VectorType &dst,
  const VectorType &src,
  const std::pair<unsigned int, unsigned int> &face_range) const
  {
  FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
  matrix_free, true, 0, 0, 0);
  FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_p(
  matrix_free, false, 0, 0, 0);
  FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
  matrix_free, true, 0, 0, 1);
  FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_p(
  matrix_free, false, 0, 0, 1);
 

Class that gives access to the material at each cell

  MaterialEvaluation material(matrix_free, *material_data);
 
  for (unsigned int face = face_range.first; face < face_range.second;
  face++)
  {
  velocity_m.reinit(face);
  velocity_p.reinit(face);
 
  pressure_m.reinit(face);
  pressure_p.reinit(face);
 
  pressure_m.gather_evaluate(src, EvaluationFlags::values);
  pressure_p.gather_evaluate(src, EvaluationFlags::values);
 
  velocity_m.gather_evaluate(src, EvaluationFlags::values);
  velocity_p.gather_evaluate(src, EvaluationFlags::values);
 
  material.reinit_face(face);
  evaluate_face_kernel<true>(pressure_m,
  velocity_m,
  pressure_p,
  velocity_p,
  material.get_speed_of_sound(),
  material.get_density());
 
  pressure_m.integrate_scatter(EvaluationFlags::values, dst);
  pressure_p.integrate_scatter(EvaluationFlags::values, dst);
  velocity_m.integrate_scatter(EvaluationFlags::values, dst);
  velocity_p.integrate_scatter(EvaluationFlags::values, dst);
  }
  }
 
 

Matrix-free boundary function for point-to-point interpolation

This function evaluates the boundary face integrals and the non-matching face integrals using point-to-point interpolation.

  template <typename VectorType>
  void local_apply_boundary_face_point_to_point(
  const MatrixFree<dim, Number> &matrix_free,
  VectorType &dst,
  const VectorType &src,
  const std::pair<unsigned int, unsigned int> &face_range) const
  {

Standard face evaluators.

  FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
  matrix_free, true, 0, 0, 0);
  FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
  matrix_free, true, 0, 0, 1);
 

Classes that return the correct BC values.

  BCEvaluationP pressure_bc(pressure_m);
  BCEvaluationU velocity_bc(velocity_m);
 

Class that gives access to the material at each cell

  MaterialEvaluation material(matrix_free, *material_data);
 

Remote evaluators.

  auto pressure_r = pressure_r_eval->get_data_accessor();
  auto velocity_r = velocity_r_eval->get_data_accessor();
  auto c_r = c_r_eval->get_data_accessor();
  auto rho_r = rho_r_eval->get_data_accessor();
 
  for (unsigned int face = face_range.first; face < face_range.second;
  face++)
  {
  velocity_m.reinit(face);
  pressure_m.reinit(face);
 
  pressure_m.gather_evaluate(src, EvaluationFlags::values);
  velocity_m.gather_evaluate(src, EvaluationFlags::values);
 
  if (HelperFunctions::is_non_matching_face(
  remote_face_ids, matrix_free.get_boundary_id(face)))
  {

If face is non-matching we have to query values via the FERemoteEvaluaton objects. This is done by passing the corresponding FERemoteEvaluaton objects to the function that evaluates the kernel. As mentioned above, each side of the non-matching interface is traversed separately and we do not have to consider the neighbor in the kernel. Note, that the values in the FERemoteEvaluaton objects are already updated at this point.

For point-to-point interpolation we simply use the corresponding FERemoteEvaluaton objects in combination with the standard FEFaceEvaluation objects.

  velocity_r.reinit(face);
  pressure_r.reinit(face);
 
  material.reinit_face(face);
 
  if (material.is_homogeneous())
  {

If homogeneous material is considered do not use the inhomogeneous fluxes. While it would be possible to use the inhomogeneous fluxes they are more expensive to compute.

  evaluate_face_kernel<false>(pressure_m,
  velocity_m,
  pressure_r,
  velocity_r,
  material.get_speed_of_sound(),
  material.get_density());
  }
  else
  {

If inhomogeneous material is considered use the in-homogeneous fluxes.

  c_r.reinit(face);
  rho_r.reinit(face);
  evaluate_face_kernel_inhomogeneous(
  pressure_m,
  velocity_m,
  pressure_r,
  velocity_r,
  material.get_speed_of_sound(),
  material.get_density(),
  c_r,
  rho_r);
  }
  }
  else
  {

If face is a standard boundary face, evaluate the integral as usual in the matrix free context. To be able to use the same kernel as for inner faces we pass the boundary condition objects to the function that evaluates the kernel. As detailed above weight_neighbor=false.

  material.reinit_face(face);
  evaluate_face_kernel<false>(pressure_m,
  velocity_m,
  pressure_bc,
  velocity_bc,
  material.get_speed_of_sound(),
  material.get_density());
  }
 
  pressure_m.integrate_scatter(EvaluationFlags::values, dst);
  velocity_m.integrate_scatter(EvaluationFlags::values, dst);
  }
  }
 

Matrix-free boundary function for Nitsche-type mortaring

This function evaluates the boundary face integrals and the non-matching face integrals using Nitsche-type mortaring.

  template <typename VectorType>
  void local_apply_boundary_face_mortaring(
  const MatrixFree<dim, Number> &matrix_free,
  VectorType &dst,
  const VectorType &src,
  const std::pair<unsigned int, unsigned int> &face_range) const
  {

Standard face evaluators for BCs.

  FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
  matrix_free, true, 0, 0, 0);
  FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
  matrix_free, true, 0, 0, 1);
 

For Nitsche-type mortaring we are evaluating the integrals over intersections. This is why, quadrature points are arbitrarily distributed on every face. Thus, we can not make use of face batches and FEFaceEvaluation but have to consider each face individually and make use of FEFacePointEvaluation to evaluate the integrals in the arbitrarily distributed quadrature points. Since the setup of FEFacePointEvaluation is more expensive than that of FEEvaluation we do the setup only once. For this we are using the helper function get_thread_safe_fe_face_point_evaluation_object().

  get_thread_safe_fe_face_point_evaluation_object<1>(
  thread_local_pressure_m_mortar, 0);
  get_thread_safe_fe_face_point_evaluation_object<dim>(
  thread_local_velocity_m_mortar, 1);
 
  BCEvaluationP pressure_bc(pressure_m);
  BCEvaluationU velocity_bc(velocity_m);
 
  MaterialEvaluation material(matrix_free, *material_data);
 
  auto pressure_r_mortar = pressure_r_eval->get_data_accessor();
  auto velocity_r_mortar = velocity_r_eval->get_data_accessor();
  auto c_r = c_r_eval->get_data_accessor();
  auto rho_r = rho_r_eval->get_data_accessor();
 
  for (unsigned int face = face_range.first; face < face_range.second;
  ++face)
  {
  if (HelperFunctions::is_non_matching_face(
  remote_face_ids, matrix_free.get_boundary_id(face)))
  {
  material.reinit_face(face);
 

First fetch the DoF values with standard FEFaceEvaluation objects.

  pressure_m.reinit(face);
  velocity_m.reinit(face);
 
  pressure_m.read_dof_values(src);
  velocity_m.read_dof_values(src);
 

Project the internally stored values into the face DoFs of the current face.

  pressure_m.project_to_face(EvaluationFlags::values);
  velocity_m.project_to_face(EvaluationFlags::values);
 

For mortaring, we have to consider every face from the face batches separately and have to use the FEFacePointEvaluation objects to be able to evaluate the integrals with the arbitrarily distributed quadrature points.

  for (unsigned int v = 0;
  v < matrix_free.n_active_entries_per_face_batch(face);
  ++v)
  {
  constexpr unsigned int n_lanes =
  velocity_m_mortar.reinit(face * n_lanes + v);
  pressure_m_mortar.reinit(face * n_lanes + v);
 

Evaluate using FEFacePointEvaluation. As buffer, simply use the internal buffers from the FEFaceEvaluation objects.

  velocity_m_mortar.evaluate_in_face(
  &velocity_m.get_scratch_data().begin()[0][v],
 
  pressure_m_mortar.evaluate_in_face(
  &pressure_m.get_scratch_data().begin()[0][v],
 
  velocity_r_mortar.reinit(face * n_lanes + v);
  pressure_r_mortar.reinit(face * n_lanes + v);
 
  if (material.is_homogeneous())
  {

If homogeneous material is considered do not use the inhomogeneous fluxes. While it would be possible to use the inhomogeneous fluxes they are more expensive to compute. Since we are operating on face v we call material.get_density()[v].

  evaluate_face_kernel<false>(
  pressure_m_mortar,
  velocity_m_mortar,
  pressure_r_mortar,
  velocity_r_mortar,
  material.get_speed_of_sound()[v],
  material.get_density()[v]);
  }
  else
  {
  c_r.reinit(face * n_lanes + v);
  rho_r.reinit(face * n_lanes + v);
 
  evaluate_face_kernel_inhomogeneous(
  pressure_m_mortar,
  velocity_m_mortar,
  pressure_r_mortar,
  velocity_r_mortar,
  material.get_speed_of_sound()[v],
  material.get_density()[v],
  c_r,
  rho_r);
  }
 

Integrate using FEFacePointEvaluation. As buffer, simply use the internal buffers from the FEFaceEvaluation objects.

  velocity_m_mortar.integrate_in_face(
  &velocity_m.get_scratch_data().begin()[0][v],
 
  pressure_m_mortar.integrate_in_face(
  &pressure_m.get_scratch_data().begin()[0][v],
  }
 

Collect the contributions from the face DoFs to the internal cell DoFs to be able to use the member function distribute_local_to_global().

  pressure_m.collect_from_face(EvaluationFlags::values);
  velocity_m.collect_from_face(EvaluationFlags::values);
 
  pressure_m.distribute_local_to_global(dst);
  velocity_m.distribute_local_to_global(dst);
  }
  else
  {

Same as in local_apply_boundary_face_point_to_point().

  velocity_m.reinit(face);
  pressure_m.reinit(face);
 
  pressure_m.gather_evaluate(src, EvaluationFlags::values);
  velocity_m.gather_evaluate(src, EvaluationFlags::values);
 
  material.reinit_face(face);
  evaluate_face_kernel<false>(pressure_m,
  velocity_m,
  pressure_bc,
  velocity_bc,
  material.get_speed_of_sound(),
  material.get_density());
 
  pressure_m.integrate_scatter(EvaluationFlags::values, dst);
  velocity_m.integrate_scatter(EvaluationFlags::values, dst);
  }
  }
  }
 
  const MatrixFree<dim, Number> &matrix_free;
 

CellwiseMaterialData is stored as shared pointer with the same argumentation.

  const std::shared_ptr<CellwiseMaterialData<Number>> material_data;
 
  const std::set<types::boundary_id> remote_face_ids;
 

FERemoteEvaluation objects are strored as shared pointers. This way, they can also be used for other operators without precomputing the values multiple times.

  const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
  pressure_r_eval;
  const std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
  velocity_r_eval;
 
  const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
  c_r_eval;
  const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
  rho_r_eval;
 
  const std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>>
  nm_mapping_info;
 

We store FEFacePointEvaluation objects as members in a thread local way, since its creation is more expensive compared to FEEvaluation objects.

  std::unique_ptr<FEFacePointEvaluation<1, dim, dim, Number>>>
  thread_local_pressure_m_mortar;
 
  std::unique_ptr<FEFacePointEvaluation<dim, dim, dim, Number>>>
  thread_local_velocity_m_mortar;
 
A class that provides a separate storage location on each thread that accesses the object.

Helper function to create and get FEFacePointEvaluation objects in a thread safe way. On each thread, FEFacePointEvaluation is created if it has not been created by now. After that, simply return the object corresponding to the thread under consideration.

  template <int n_components>
  get_thread_safe_fe_face_point_evaluation_object(
  &fe_face_point_eval_thread_local,
  unsigned int fist_selected_comp) const
  {
  if (fe_face_point_eval_thread_local.get() == nullptr)
  {
  fe_face_point_eval_thread_local = std::make_unique<
  *nm_mapping_info,
  matrix_free.get_dof_handler().get_fe(),
  true,
  fist_selected_comp);
  }
  return *fe_face_point_eval_thread_local.get();
  }
  };
 

Inverse mass operator

Class to apply the inverse mass operator.

  template <int dim, typename Number>
  class InverseMassOperator
  {
  public:
  InverseMassOperator(const MatrixFree<dim, Number> &matrix_free)
  : matrix_free(matrix_free)
  {}
 

Function to apply the inverse mass operator.

  template <typename VectorType>
  void apply(VectorType &dst, const VectorType &src) const
  {
  dst.zero_out_ghost_values();
  matrix_free.cell_loop(&InverseMassOperator::local_apply_cell<VectorType>,
  this,
  dst,
  src);
  }
 
  private:
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const

Apply the inverse mass operator onto every cell batch.

  template <typename VectorType>
  void local_apply_cell(
  VectorType &dst,
  const VectorType &src,
  const std::pair<unsigned int, unsigned int> &cell_range) const
  {
  FEEvaluation<dim, -1, 0, dim + 1, Number> phi(mf);
  minv(phi);
 
  for (unsigned int cell = cell_range.first; cell < cell_range.second;
  ++cell)
  {
  phi.reinit(cell);
  phi.read_dof_values(src);
  minv.apply(phi.begin_dof_values(), phi.begin_dof_values());
  phi.set_dof_values(dst);
  }
  }
 
  const MatrixFree<dim, Number> &matrix_free;
  };
 

Runge-Kutta time-stepping

This class implements a Runge-Kutta scheme of order 2.

  template <int dim, typename Number, typename remote_value_type>
  class RungeKutta2
  {
 
  public:
  RungeKutta2(
  const std::shared_ptr<InverseMassOperator<dim, Number>>
  inverse_mass_operator,
  const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
  acoustic_operator)
  : inverse_mass_operator(inverse_mass_operator)
  , acoustic_operator(acoustic_operator)
  {}
 

Set up and run time loop.

  void run(const MatrixFree<dim, Number> &matrix_free,
  const double cr,
  const double end_time,
  const double speed_of_sound,
  const Function<dim> &initial_condition,
  const std::string &vtk_prefix)
  {

Get needed members of matrix free.

  const auto &dof_handler = matrix_free.get_dof_handler();
  const auto &mapping = *matrix_free.get_mapping_info().mapping;
  const auto degree = dof_handler.get_fe().degree;
 

Initialize needed Vectors.

  VectorType solution;
  matrix_free.initialize_dof_vector(solution);
  VectorType solution_temp;
  matrix_free.initialize_dof_vector(solution_temp);
 
void initialize_dof_vector(VectorType &vec, const unsigned int dof_handler_index=0) const

Set the initial condition.

  HelperFunctions::set_initial_condition(matrix_free,
  initial_condition,
  solution);
 

Compute time step size: Compute minimum element edge length. We assume non-distorted elements, therefore we only compute the distance between two vertices

  double h_local_min = std::numeric_limits<double>::max();
  for (const auto &cell : dof_handler.active_cell_iterators())
  h_local_min =
  std::min(h_local_min,
  (cell->vertex(1) - cell->vertex(0)).norm_square());
  h_local_min = std::sqrt(h_local_min);
  const double h_min =
  Utilities::MPI::min(h_local_min, dof_handler.get_communicator());
 
T min(const T &t, const MPI_Comm mpi_communicator)
STL namespace.

Compute constant time step size via the CFL condition.

  const double dt =
  cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound);
 

Perform time integration loop.

  double time = 0.0;
  unsigned int timestep = 0;
  while (time < end_time)
  {

Write output.

  HelperFunctions::write_vtu(solution,
  matrix_free.get_dof_handler(),
  mapping,
  degree,
  "step_89-" + vtk_prefix +
  std::to_string(timestep));
 

Perform a single time step.

  std::swap(solution, solution_temp);
  time += dt;
  timestep++;
  perform_time_step(dt, solution, solution_temp);
  }
  }
 
  private:

Perform one Runge-Kutta 2 time step.

  void
  perform_time_step(const double dt, VectorType &dst, const VectorType &src)
  {
  VectorType k1 = src;
 

First stage.

  evaluate_stage(k1, src);
 

Second stage.

  k1.sadd(0.5 * dt, 1.0, src);
  evaluate_stage(dst, k1);
  dst.sadd(dt, 1.0, src);
  }
 

Evaluate a single Runge-Kutta stage.

  void evaluate_stage(VectorType &dst, const VectorType &src)
  {

Evaluate the stage

  acoustic_operator->evaluate(dst, src);
  dst *= -1.0;
  inverse_mass_operator->apply(dst, dst);
  }
 

Needed operators.

  const std::shared_ptr<InverseMassOperator<dim, Number>>
  inverse_mass_operator;
  const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
  acoustic_operator;
  };
 
 

Construction of non-matching triangulations

This function creates a two dimensional squared triangulation that spans from (0,0) to (1,1). It consists of two sub-domains. The left sub-domain spans from (0,0) to (0.525,1). The right sub-domain spans from (0.525,0) to (1,1). The left sub-domain has three times smaller elements compared to the right sub-domain.

  template <int dim>
  void build_non_matching_triangulation(
  std::set<types::boundary_id> &non_matching_faces,
  const unsigned int refinements)
  {
  const double length = 1.0;
 

At non-matching interfaces, we provide different boundary IDs. These boundary IDs have to differ because later on RemotePointEvaluation has to search for remote points for each face, that are defined in the same mesh (since we merge the mesh) but not on the same side of the non-matching interface.

  const types::boundary_id non_matching_id_left = 98;
  const types::boundary_id non_matching_id_right = 99;
 

Provide this information to the caller.

  non_matching_faces.insert(non_matching_id_left);
  non_matching_faces.insert(non_matching_id_right);
 

Construct left part of mesh.

  Triangulation<dim> tria_left;
  const unsigned int subdiv_left = 11;
  {subdiv_left, 2 * subdiv_left},
  {0.0, 0.0},
  {0.525 * length, length});
 
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)

The left part of the mesh has the material ID 0.

  for (const auto &cell : tria_left.active_cell_iterators())
  cell->set_material_id(0);
 

The right face is non-matching. All other boundary IDs are set to 0.

  for (const auto &face : tria_left.active_face_iterators())
  if (face->at_boundary())
  {
  face->set_boundary_id(0);
  if (face->center()[0] > 0.525 * length - 1e-6)
  face->set_boundary_id(non_matching_id_left);
  }
 

Construct right part of mesh.

  Triangulation<dim> tria_right;
  const unsigned int subdiv_right = 4;
  {subdiv_right, 2 * subdiv_right},
  {0.525 * length, 0.0},
  {length, length});
 

The right part of the mesh has the material ID 1.

  for (const auto &cell : tria_right.active_cell_iterators())
  cell->set_material_id(1);
 

The left face is non-matching. All other boundary IDs are set to 0.

  for (const auto &face : tria_right.active_face_iterators())
  if (face->at_boundary())
  {
  face->set_boundary_id(0);
  if (face->center()[0] < 0.525 * length + 1e-6)
  face->set_boundary_id(non_matching_id_right);
  }
 

Merge triangulations with tolerance 0 to ensure no vertices are merged, see the documentation of the function merge_triangulations().

  tria_right,
  tria,
  /*tolerance*/ 0.,
  /*copy_manifold_ids*/ false,
  /*copy_boundary_ids*/ true);
  tria.refine_global(refinements);
  }
 
void refine_global(const unsigned int times=1)
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)

Set up and run point-to-point interpolation

The main purpose of this function is to fill a FERemoteEvaluationCommunicator object that is needed for point-to-point interpolation. Additionally, the corresponding remote evaluators are set up using this remote communicator. Eventually, the operators are handed to the time integrator that runs the simulation.

  template <int dim, typename Number>
  void run_with_point_to_point_interpolation(
  const MatrixFree<dim, Number> &matrix_free,
  const std::set<types::boundary_id> &non_matching_faces,
  const std::map<types::material_id, std::pair<double, double>> &materials,
  const double end_time,
  const Function<dim> &initial_condition,
  const std::string &vtk_prefix)
  {
  const auto &dof_handler = matrix_free.get_dof_handler();
  const auto &tria = dof_handler.get_triangulation();
 
Triangulation< dim, spacedim > & get_triangulation()

Communication objects know about the communication pattern. I.e., they know about the cells and quadrature points that have to be evaluated at remote faces. This information is given via RemotePointEvaluation. Additionally, the communication objects have to be able to match the quadrature points of the remote points (that provide exterior information) to the quadrature points defined at the interior cell. In case of point-to-point interpolation a vector of pairs with face batch Ids and the number of faces in the batch is needed. FERemoteCommunicationObjectEntityBatches is a container to store this information.

The information is filled outside of the actual class since in some cases the information is available from some heuristic and it is possible to skip some expensive operations. This is for example the case for sliding rotating interfaces with equally spaced elements on both sides of the non-matching interface [72].

For the standard case of point to point-to-point interpolation without any heuristic we make use of the utility function compute_remote_communicator_faces_point_to_point_interpolation(). Please refer to this function to see how to manually set up the remote communicator from outside.

  std::vector<
  std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
  non_matching_faces_marked_vertices;
 
  for (const auto &nm_face : non_matching_faces)
  {

Sufficient lambda, that rules out all cells connected to the current side of the non-matching interface to avoid self intersections.

  auto marked_vertices = [&]() {

only search points at cells that are not connected to nm_face

  std::vector<bool> mask(tria.n_vertices(), true);
 
  for (const auto &cell : tria.active_cell_iterators())
  for (auto const &f : cell->face_indices())
  if (cell->face(f)->at_boundary() &&
  cell->face(f)->boundary_id() == nm_face)
  for (const auto v : cell->vertex_indices())
  mask[cell->vertex_index(v)] = false;
 
  return mask;
  };
 
  non_matching_faces_marked_vertices.emplace_back(
  std::make_pair(nm_face, marked_vertices));
  }
 
  auto remote_communicator =
  matrix_free, non_matching_faces_marked_vertices);
 
unsigned int n_vertices() const
unsigned int vertex_indices[2]
FERemoteEvaluationCommunicator< dim > compute_remote_communicator_faces_point_to_point_interpolation(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::vector< std::pair< types::boundary_id, std::function< std::vector< bool >()> > > &non_matching_faces_marked_vertices, const unsigned int quad_no=0, const unsigned int dof_no=0, const double tolerance=1e-9)

We are using point-to-point interpolation and can therefore easily access the corresponding data at face batches. This is why we use a VectorizedArray as remote_value_type

  using remote_value_type = VectorizedArray<Number>;
 

Set up FERemoteEvaluation object that accesses the pressure at remote faces.

  const auto pressure_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator, dof_handler, /*first_selected_component*/ 0);
 

Set up FERemoteEvaluation object that accesses the velocity at remote faces.

  const auto velocity_r =
  std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
  remote_communicator, dof_handler, /*first_selected_component*/ 1);
 

Set up cell-wise material data.

  const auto material_data =
  std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
 

If we have an inhomogeneous problem, we have to set up the material handler that accesses the materials at remote faces.

  const auto c_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator,
  matrix_free.get_dof_handler().get_triangulation(),
  /*first_selected_component*/ 0);
  const auto rho_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator,
  matrix_free.get_dof_handler().get_triangulation(),
  /*first_selected_component*/ 0);
 
  if (!material_data->is_homogeneous())
  {

Initialize and fill DoF vectors that contain the materials.

  matrix_free.get_dof_handler().get_triangulation().n_active_cells());
  matrix_free.get_dof_handler().get_triangulation().n_active_cells());
 
  for (const auto &cell : matrix_free.get_dof_handler()
  .get_triangulation()
  .active_cell_iterators())
  {
  c[cell->active_cell_index()] =
  materials.at(cell->material_id()).first;
  rho[cell->active_cell_index()] =
  materials.at(cell->material_id()).second;
  }
 
Point< 2 > first
Definition grid_out.cc:4623

Materials do not change during the simulation, therefore there is no need to precompute the values after the first gather_evaluate() again.

  c_r->gather_evaluate(c, EvaluationFlags::values);
  rho_r->gather_evaluate(rho, EvaluationFlags::values);
  }
 
 

Set up inverse mass operator.

  const auto inverse_mass_operator =
  std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
 

Set up the acoustic operator. Using remote_value_type=VectorizedArray<Number> makes the operator use point-to-point interpolation.

  const auto acoustic_operator =
  std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
  matrix_free,
  material_data,
  non_matching_faces,
  pressure_r,
  velocity_r,
  c_r,
  rho_r);
 

Compute the the maximum speed of sound, needed for the computation of the time-step size.

  double speed_of_sound_max = 0.0;
  for (const auto &mat : materials)
  speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
 
Point< 2 > second
Definition grid_out.cc:4624

Set up time integrator.

  RungeKutta2<dim, Number, remote_value_type> time_integrator(
  inverse_mass_operator, acoustic_operator);
 

For considered examples, we found a limiting Courant number of \(\mathrm{Cr}\approx 0.36\) to maintain stability. To ensure, the error of the temporal discretization is small, we use a considerably smaller Courant number of \(0.2\).

  time_integrator.run(matrix_free,
  /*Cr*/ 0.2,
  end_time,
  speed_of_sound_max,
  initial_condition,
  vtk_prefix);
  }
 

Set up and run Nitsche-type mortaring

The main purpose of this function is to fill a FERemoteEvaluationCommunicator object that is needed for Nitsche-type mortaring. Additionally, the corresponding remote evaluators are set up using this remote communicator. Eventually, the operators are handed to the time integrator that runs the simulation.

  template <int dim, typename Number>
  void run_with_nitsche_type_mortaring(
  const MatrixFree<dim, Number> &matrix_free,
  const std::set<types::boundary_id> &non_matching_faces,
  const std::map<types::material_id, std::pair<double, double>> &materials,
  const double end_time,
  const Function<dim> &initial_condition,
  const std::string &vtk_prefix)
  {
  #ifndef DEAL_II_WITH_CGAL
  (void)matrix_free;
  (void)non_matching_faces;
  (void)materials;
  (void)end_time;
  (void)initial_condition;
  (void)vtk_prefix;
 
  std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
 
  pcout << "In this function, mortars are computed using CGAL. "
  "Configure deal.II with DEAL_II_WITH_CGAL to run this function.\n";
 
  return;
  #else
 
  const auto &dof_handler = matrix_free.get_dof_handler();
  const auto &tria = dof_handler.get_triangulation();
  const auto &mapping = *matrix_free.get_mapping_info().mapping;
  const auto n_quadrature_pnts = matrix_free.get_quadrature().size();
 
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
Definition mpi.cc:107

In case of Nitsche-type mortaring a vector of pairs with cell iterator and face number is needed as communication object. FERemoteCommunicationObjectFaces is a container to store this information.

For the standard case of Nitsche-type mortaring without any heuristic we make use of the utility function compute_remote_communicator_faces_nitsche_type_mortaring(). Please refer to this function to see how to manually set up the remote communicator from outside and how to reinit NonMatching::MappingInfo.

  std::vector<
  std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
  non_matching_faces_marked_vertices;
 
  for (const auto &nm_face : non_matching_faces)
  {

Sufficient lambda, that rules out all cells connected to the current side of the non-matching interface to avoid self intersections.

  auto marked_vertices = [&]() {

only search points at cells that are not connected to nm_face

  std::vector<bool> mask(tria.n_vertices(), true);
 
  for (const auto &cell : tria.active_cell_iterators())
  for (auto const &f : cell->face_indices())
  if (cell->face(f)->at_boundary() &&
  cell->face(f)->boundary_id() == nm_face)
  for (const auto v : cell->vertex_indices())
  mask[cell->vertex_index(v)] = false;
 
  return mask;
  };
 
  non_matching_faces_marked_vertices.emplace_back(
  std::make_pair(nm_face, marked_vertices));
  }
 

Quadrature points are arbitrarily distributed on each non-matching face. Therefore, we have to make use of FEFacePointEvaluation. FEFacePointEvaluation needs NonMatching::MappingInfo to work at the correct quadrature points that are in sync with used FERemoteEvaluation object. Using compute_remote_communicator_faces_nitsche_type_mortaring() to reinit NonMatching::MappingInfo ensures this. In the case of mortaring, we have to use the weights provided by the quadrature rules that are used to set up NonMatching::MappingInfo. Therefore we set the flag use_global_weights.

Set up NonMatching::MappingInfo with needed update flags and additional_data.

  auto nm_mapping_info =
  std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
  mapping,
  additional_data);
 
  auto remote_communicator =
  matrix_free,
  non_matching_faces_marked_vertices,
  n_quadrature_pnts,
  0,
  nm_mapping_info.get());
 
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
FERemoteEvaluationCommunicator< dim > compute_remote_communicator_faces_nitsche_type_mortaring(const MatrixFree< dim, Number, VectorizedArrayType > &matrix_free, const std::vector< std::pair< types::boundary_id, std::function< std::vector< bool >()> > > &non_matching_faces_marked_vertices, const unsigned int n_q_pnts_1D, const unsigned int dof_no=0, NonMatching::MappingInfo< dim, dim, Number > *nm_mapping_info=nullptr, const double tolerance=1e-9)

Same as above but since quadrature points are aribtrarily distributed we have to consider each face in a batch separately and can not make use of VecorizedArray.

  using remote_value_type = Number;
 
  const auto pressure_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator, dof_handler, /*first_selected_component*/ 0);
 
  const auto velocity_r =
  std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
  remote_communicator, dof_handler, /*first_selected_component*/ 1);
 
  const auto material_data =
  std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
 
  const auto c_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator,
  matrix_free.get_dof_handler().get_triangulation(),
  /*first_selected_component*/ 0);
  const auto rho_r =
  std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
  remote_communicator,
  matrix_free.get_dof_handler().get_triangulation(),
  /*first_selected_component*/ 0);
 
  if (!material_data->is_homogeneous())
  {
  matrix_free.get_dof_handler().get_triangulation().n_active_cells());
  matrix_free.get_dof_handler().get_triangulation().n_active_cells());
 
  for (const auto &cell : matrix_free.get_dof_handler()
  .get_triangulation()
  .active_cell_iterators())
  {
  c[cell->active_cell_index()] =
  materials.at(cell->material_id()).first;
  rho[cell->active_cell_index()] =
  materials.at(cell->material_id()).second;
  }
 
  c_r->gather_evaluate(c, EvaluationFlags::values);
  rho_r->gather_evaluate(rho, EvaluationFlags::values);
  }
 

Set up inverse mass operator.

  const auto inverse_mass_operator =
  std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
 

Set up the acoustic operator. Using remote_value_type=Number makes the operator use Nitsche-type mortaring.

  const auto acoustic_operator =
  std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
  matrix_free,
  material_data,
  non_matching_faces,
  pressure_r,
  velocity_r,
  c_r,
  rho_r,
  nm_mapping_info);
 

Compute the the maximum speed of sound, needed for the computation of the time-step size.

  double speed_of_sound_max = 0.0;
  for (const auto &mat : materials)
  speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
 
 

Set up time integrator.

  RungeKutta2<dim, Number, remote_value_type> time_integrator(
  inverse_mass_operator, acoustic_operator);
 

Run time loop with Courant number \(0.2\).

  time_integrator.run(matrix_free,
  /*Cr*/ 0.2,
  end_time,
  speed_of_sound_max,
  initial_condition,
  vtk_prefix);
  #endif
  }
  } // namespace Step89
 
 

main()

Finally, the main() function executes the different versions of handling non-matching interfaces.

  int main(int argc, char *argv[])
  {
  using namespace dealii;
  constexpr int dim = 2;
  using Number = double;
 
  std::cout.precision(5);
  ConditionalOStream pcout(std::cout,
  (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
  0));
 
  const unsigned int refinements = 1;
  const unsigned int degree = 3;
 

Construct non-matching triangulation and fill non-matching boundary IDs.

Similar to step-87, the minimum requirement of this tutorial is MPI. parallel::distributed::Triangulation is used if deal.II is configured with p4est. Otherwise parallel::shared::Triangulation is used.

  #ifdef DEAL_II_WITH_P4EST
  #else
  parallel::shared::Triangulation<dim> tria(MPI_COMM_WORLD);
  #endif
 
  pcout << "Create non-matching grid..." << std::endl;
 
  std::set<types::boundary_id> non_matching_faces;
  Step89::build_non_matching_triangulation(tria,
  non_matching_faces,
  refinements);
 
  pcout << " - Refinement level: " << refinements << std::endl;
  pcout << " - Number of cells: " << tria.n_cells() << std::endl;
 
unsigned int n_cells() const

Set up MatrixFree.

  pcout << "Create DoFHandler..." << std::endl;
  DoFHandler<dim> dof_handler(tria);
  dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree), dim + 1));
  pcout << " - Number of DoFs: " << dof_handler.n_dofs() << std::endl;
 
  constraints.close();
 
  pcout << "Set up MatrixFree..." << std::endl;
  data.mapping_update_flags_inner_faces = update_values;
  data.mapping_update_flags_boundary_faces =
 
  MatrixFree<dim, Number> matrix_free;
  matrix_free.reinit(
  MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
 
 
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
@ update_gradients
Shape function gradients.

Run vibrating membrane test case

  pcout << "Run vibrating membrane test case..." << std::endl;

Vibrating membrane test case:

Homogeneous pressure DBCs are applied for simplicity. Therefore, modes can not be chosen arbitrarily.

  const double modes = 10.0;
  std::map<types::material_id, std::pair<double, double>> homogeneous_material;
  homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0);
  const auto initial_solution_membrane =
  Step89::InitialConditionVibratingMembrane<dim>(modes);
 
  pcout << " - Point-to-point interpolation: " << std::endl;
const types::material_id invalid_material_id
Definition types.h:277

Run vibrating membrane test case using point-to-point interpolation:

  Step89::run_with_point_to_point_interpolation(
  matrix_free,
  non_matching_faces,
  homogeneous_material,
  8.0 * initial_solution_membrane.get_period_duration(
  homogeneous_material.begin()->second.first),
  initial_solution_membrane,
  "vm-p2p");
 
  pcout << " - Nitsche-type mortaring: " << std::endl;

Run vibrating membrane test case using Nitsche-type mortaring:

  Step89::run_with_nitsche_type_mortaring(
  matrix_free,
  non_matching_faces,
  homogeneous_material,
  8.0 * initial_solution_membrane.get_period_duration(
  homogeneous_material.begin()->second.first),
  initial_solution_membrane,
  "vm-nitsche");
 

Run test case with in-homogeneous material

  pcout << "Run test case with in-homogeneous material..." << std::endl;

In-homogeneous material test case:

Run simple test case with in-homogeneous material and Nitsche-type mortaring:

  std::map<types::material_id, std::pair<double, double>>
  inhomogeneous_material;
  inhomogeneous_material[0] = std::make_pair(1.0, 1.0);
  inhomogeneous_material[1] = std::make_pair(3.0, 1.0);
  Step89::run_with_nitsche_type_mortaring(matrix_free,
  non_matching_faces,
  inhomogeneous_material,
  /*runtime*/ 0.3,
  Step89::GaussPulse<dim>(0.3, 0.5),
  "inhomogeneous");
 
 
  return 0;
  }

Results

Vibrating membrane: Point-to-point interpolation vs. Nitsche-type mortaring

We compare the results of the simulations after the last time step, i.e. at \(t=8T\). The \(y\)-component of the velocity field using Nitsche-type mortaring is depicted on the left. The same field using point-to-point interpolation is depicted on the right.

Besides this, the results for the pressure and the velocity in \(y\) direction are plotted along the horizontal line that spans from (0,0.69) to (1,0.69).

While the results of the pressure is similar, \(u_y\) differs clearly. At certain positions we can see aliasing errors for the point-to-point interpolation. For different mesh configurations and/or longer run times, the aliasing effects of the point-to-point simulation accumulate and the simulation gets instable. To keep the tutorial short we have chosen one mesh that can be used for all examples. For a configuration that yields instable results for a wide range of polynomial degrees, see [106].

Wave propagation through in-homogeneous fluid

This is just one example in which non-matching discretizations can be efficiently used to reduce the amount of DoFs. The example is nice, since results for a similar test case are shown in multiple publications. As before, we slightly adapted the test case to be able to use the same mesh for all simulations. The pressure field at \(t=0.3\) is depicted below.

As expected, we can easily see that the wave length in the right domain is roughly three times times the wave length in the left domain. Hence, the wave can be resolved with a coarser discretization.

Using the same element size in the whole domain, we can compute a reference result. The displayed reference result is obtained by choosing the same subdivision level for both sub-domains, i.e. subdiv_right = 11. In this particular example the reference result uses \(92928\) DoFs, while the non-matching result uses \(52608\) DoFs. The pressure result is plotted along the horizontal line that spans from (0,0.5) to (1,0.5).

The results, obtained with the non-matching discretization is in good agreement with the reference result.

Possibilities for extensions

All the implementations are done with overlapping triangulations in mind. In particular the intersections in the mortaring case are constructed such that they are computed correctly for overlapping triangulations. For this the intersection requests are of dimension \(dim-1\). The cells which are intersected with the intersection requests are of dimension \(dim\). For the simple case of non-conforming interfaces it would be sufficient to compute the intersections between \(dim-1\) and \(dim-1\) entities. Furthermore, the lambda could be adapted, such that cells are only marked if they are connected to a certain boundary ID (in this case, e.g. 99) instead of marking every cell that is not connected to the opposite boundary ID (in this case, e.g. 98). Marking less cells can reduce the setup costs significantly.

Note that for in-homogeneous material in this procedure is questionable, since it is not clear which material is present in the overlapping part of the mesh.

The plain program

/* ---------------------------------------------------------------------
*
* Copyright (C) 2023 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.
*
* ---------------------------------------------------------------------
*
*
* Authors: Johannes Heinz, TU Wien, 2023
* Maximilian Bergbauer, TUM, 2023
* Marco Feder, SISSA, 2023
* Peter Munch, University of Augsburg/Uppsala University, 2023
*/
namespace Step89
{
using namespace dealii;
template <int dim>
class InitialConditionVibratingMembrane : public Function<dim>
{
public:
InitialConditionVibratingMembrane(const double modes);
double value(const Point<dim> &p, const unsigned int comp) const final;
double get_period_duration(const double speed_of_sound) const;
private:
const double M;
};
template <int dim>
InitialConditionVibratingMembrane<dim>::InitialConditionVibratingMembrane(
const double modes)
: Function<dim>(dim + 1, 0.0)
, M(modes)
{
static_assert(dim == 2, "Only implemented for dim==2");
}
template <int dim>
double
InitialConditionVibratingMembrane<dim>::value(const Point<dim> &p,
const unsigned int comp) const
{
if (comp == 0)
return std::sin(M * numbers::PI * p[0]) *
std::sin(M * numbers::PI * p[1]);
return 0.0;
}
template <int dim>
double InitialConditionVibratingMembrane<dim>::get_period_duration(
const double speed_of_sound) const
{
return 2.0 / (M * std::sqrt(dim) * speed_of_sound);
}
template <int dim>
class GaussPulse : public Function<dim>
{
public:
GaussPulse(const double shift_x, const double shift_y);
double value(const Point<dim> &p, const unsigned int comp) const final;
private:
const double shift_x;
const double shift_y;
};
template <int dim>
GaussPulse<dim>::GaussPulse(const double shift_x, const double shift_y)
: Function<dim>(dim + 1, 0.0)
, shift_x(shift_x)
, shift_y(shift_y)
{
static_assert(dim == 2, "Only implemented for dim==2");
}
template <int dim>
double GaussPulse<dim>::value(const Point<dim> &p,
const unsigned int comp) const
{
if (comp == 0)
return std::exp(-1000.0 * ((std::pow(p[0] - shift_x, 2)) +
(std::pow(p[1] - shift_y, 2))));
return 0.0;
}
namespace HelperFunctions
{
bool is_non_matching_face(
const std::set<types::boundary_id> &non_matching_face_ids,
const types::boundary_id face_id)
{
return non_matching_face_ids.find(face_id) != non_matching_face_ids.end();
}
template <int dim, typename Number, typename VectorType>
void set_initial_condition(MatrixFree<dim, Number> matrix_free,
const Function<dim> &initial_solution,
VectorType &dst)
{
matrix_free.get_dof_handler(),
initial_solution,
dst);
}
double
compute_dt_cfl(const double hmin, const unsigned int degree, const double c)
{
return hmin / (std::pow(degree, 1.5) * c);
}
template <typename VectorType, int dim>
void write_vtu(const VectorType &solution,
const DoFHandler<dim> &dof_handler,
const Mapping<dim> &mapping,
const unsigned int degree,
const std::string &name_prefix)
{
DataOut<dim> data_out;
data_out.set_flags(flags);
std::vector<DataComponentInterpretation::DataComponentInterpretation>
interpretation(
std::vector<std::string> names(dim + 1, "U");
names[0] = "P";
data_out.add_data_vector(dof_handler, solution, names, interpretation);
data_out.build_patches(mapping, degree, DataOut<dim>::curved_inner_cells);
data_out.write_vtu_in_parallel(name_prefix + ".vtu",
dof_handler.get_communicator());
}
} // namespace HelperFunctions
template <typename Number>
class CellwiseMaterialData
{
public:
template <int dim>
CellwiseMaterialData(
const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
const std::map<types::material_id, std::pair<double, double>>
&material_id_map)
: homogeneous(material_id_map.size() == 1)
{
Assert(material_id_map.size() > 0,
ExcMessage("No materials given to CellwiseMaterialData"));
if (homogeneous)
{
speed_of_sound_homogeneous = material_id_map.begin()->second.first;
density_homogeneous = material_id_map.begin()->second.second;
}
else
{
const auto n_cell_batches =
matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches();
speed_of_sound.resize(n_cell_batches);
density.resize(n_cell_batches);
for (unsigned int cell = 0; cell < n_cell_batches; ++cell)
{
speed_of_sound[cell] = 1.;
density[cell] = 1.;
for (unsigned int v = 0;
v < matrix_free.n_active_entries_per_cell_batch(cell);
++v)
{
const auto material_id =
matrix_free.get_cell_iterator(cell, v)->material_id();
speed_of_sound[cell][v] =
material_id_map.at(material_id).first;
density[cell][v] = material_id_map.at(material_id).second;
}
}
}
}
bool is_homogeneous() const
{
return homogeneous;
}
const AlignedVector<VectorizedArray<Number>> &get_speed_of_sound() const
{
Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()"));
return speed_of_sound;
}
const AlignedVector<VectorizedArray<Number>> &get_density() const
{
Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()"));
return density;
}
VectorizedArray<Number> get_homogeneous_speed_of_sound() const
{
Assert(homogeneous, ExcMessage("Use get_speed_of_sound()"));
return speed_of_sound_homogeneous;
}
VectorizedArray<Number> get_homogeneous_density() const
{
Assert(homogeneous, ExcMessage("Use get_density()"));
return density_homogeneous;
}
private:
const bool homogeneous;
VectorizedArray<Number> speed_of_sound_homogeneous;
VectorizedArray<Number> density_homogeneous;
};
template <int dim, typename Number>
class MaterialEvaluation
{
public:
MaterialEvaluation(
const MatrixFree<dim, Number, VectorizedArray<Number>> &matrix_free,
const CellwiseMaterialData<Number> &material_data)
: phi(matrix_free)
, phi_face(matrix_free, true)
, material_data(material_data)
{
if (material_data.is_homogeneous())
{
speed_of_sound = material_data.get_homogeneous_speed_of_sound();
density = material_data.get_homogeneous_density();
}
}
bool is_homogeneous() const
{
return material_data.is_homogeneous();
}
void reinit_cell(const unsigned int cell)
{
if (!material_data.is_homogeneous())
{
phi.reinit(cell);
speed_of_sound =
phi.read_cell_data(material_data.get_speed_of_sound());
density = phi.read_cell_data(material_data.get_density());
}
}
void reinit_face(const unsigned int face)
{
if (!material_data.is_homogeneous())
{
phi_face.reinit(face);
speed_of_sound =
phi_face.read_cell_data(material_data.get_speed_of_sound());
density = phi_face.read_cell_data(material_data.get_density());
}
}
VectorizedArray<Number> get_speed_of_sound() const
{
return speed_of_sound;
}
VectorizedArray<Number> get_density() const
{
return density;
}
private:
FEEvaluation<dim, -1, 0, 1, Number> phi;
FEFaceEvaluation<dim, -1, 0, 1, Number> phi_face;
const CellwiseMaterialData<Number> &material_data;
VectorizedArray<Number> speed_of_sound;
};
template <int dim, typename Number>
class BCEvaluationP
{
public:
BCEvaluationP(const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m)
: pressure_m(pressure_m)
{}
typename FEFaceEvaluation<dim, -1, 0, 1, Number>::value_type
get_value(const unsigned int q) const
{
return -pressure_m.get_value(q);
}
private:
const FEFaceEvaluation<dim, -1, 0, 1, Number> &pressure_m;
};
template <int dim, typename Number>
class BCEvaluationU
{
public:
BCEvaluationU(const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m)
: velocity_m(velocity_m)
{}
typename FEFaceEvaluation<dim, -1, 0, dim, Number>::value_type
get_value(const unsigned int q) const
{
return velocity_m.get_value(q);
}
private:
const FEFaceEvaluation<dim, -1, 0, dim, Number> &velocity_m;
};
template <int dim, typename Number, typename remote_value_type>
class AcousticOperator
{
static constexpr bool use_mortaring =
std::is_floating_point_v<remote_value_type>;
public:
AcousticOperator(
const MatrixFree<dim, Number> &matrix_free,
std::shared_ptr<CellwiseMaterialData<Number>> material_data,
const std::set<types::boundary_id> &remote_face_ids,
pressure_r_eval,
velocity_r_eval,
std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>> rho_r_eval,
std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>> nm_info =
nullptr)
: matrix_free(matrix_free)
, material_data(material_data)
, remote_face_ids(remote_face_ids)
, pressure_r_eval(pressure_r_eval)
, velocity_r_eval(velocity_r_eval)
, c_r_eval(c_r_eval)
, rho_r_eval(rho_r_eval)
, nm_mapping_info(nm_info)
{
if (use_mortaring)
Assert(nm_info,
ExcMessage(
"In case of Nitsche-type mortaring NonMatching::MappingInfo \
has to be provided."));
}
template <typename VectorType>
void evaluate(VectorType &dst, const VectorType &src) const
{
pressure_r_eval->gather_evaluate(src, EvaluationFlags::values);
velocity_r_eval->gather_evaluate(src, EvaluationFlags::values);
if constexpr (use_mortaring)
{
matrix_free.loop(
&AcousticOperator::local_apply_cell<VectorType>,
&AcousticOperator::local_apply_face<VectorType>,
&AcousticOperator::local_apply_boundary_face_mortaring<VectorType>,
this,
dst,
src,
true,
}
else
{
matrix_free.loop(
&AcousticOperator::local_apply_cell<VectorType>,
&AcousticOperator::local_apply_face<VectorType>,
&AcousticOperator::local_apply_boundary_face_point_to_point<
VectorType>,
this,
dst,
src,
true,
}
}
private:
template <typename VectorType>
void local_apply_cell(
const MatrixFree<dim, Number> &matrix_free,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
FEEvaluation<dim, -1, 0, 1, Number> pressure(matrix_free, 0, 0, 0);
FEEvaluation<dim, -1, 0, dim, Number> velocity(matrix_free, 0, 0, 1);
MaterialEvaluation material(matrix_free, *material_data);
for (unsigned int cell = cell_range.first; cell < cell_range.second;
++cell)
{
velocity.reinit(cell);
pressure.reinit(cell);
pressure.gather_evaluate(src, EvaluationFlags::gradients);
velocity.gather_evaluate(src, EvaluationFlags::gradients);
material.reinit_cell(cell);
const auto c = material.get_speed_of_sound();
const auto rho = material.get_density();
for (unsigned int q : pressure.quadrature_point_indices())
{
pressure.submit_value(rho * c * c * velocity.get_divergence(q),
q);
velocity.submit_value(1.0 / rho * pressure.get_gradient(q), q);
}
pressure.integrate_scatter(EvaluationFlags::values, dst);
velocity.integrate_scatter(EvaluationFlags::values, dst);
}
}
template <bool weight_neighbor,
typename InternalFaceIntegratorPressure,
typename InternalFaceIntegratorVelocity,
typename ExternalFaceIntegratorPressure,
typename ExternalFaceIntegratorVelocity>
inline DEAL_II_ALWAYS_INLINE void evaluate_face_kernel(
InternalFaceIntegratorPressure &pressure_m,
InternalFaceIntegratorVelocity &velocity_m,
ExternalFaceIntegratorPressure &pressure_p,
ExternalFaceIntegratorVelocity &velocity_p,
const typename InternalFaceIntegratorPressure::value_type c,
const typename InternalFaceIntegratorPressure::value_type rho) const
{
const auto tau = 0.5 * rho * c;
const auto gamma = 0.5 / (rho * c);
for (unsigned int q : pressure_m.quadrature_point_indices())
{
const auto n = pressure_m.normal_vector(q);
const auto pm = pressure_m.get_value(q);
const auto um = velocity_m.get_value(q);
const auto pp = pressure_p.get_value(q);
const auto up = velocity_p.get_value(q);
const auto momentum_flux =
0.5 * (pm + pp) + 0.5 * tau * (um - up) * n;
velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
if constexpr (weight_neighbor)
velocity_p.submit_value(1.0 / rho * (momentum_flux - pp) * (-n), q);
const auto mass_flux = 0.5 * (um + up) + 0.5 * gamma * (pm - pp) * n;
pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
if constexpr (weight_neighbor)
pressure_p.submit_value(rho * c * c * (mass_flux - up) * (-n), q);
}
}
template <typename InternalFaceIntegratorPressure,
typename InternalFaceIntegratorVelocity,
typename ExternalFaceIntegratorPressure,
typename ExternalFaceIntegratorVelocity,
typename MaterialIntegrator>
void evaluate_face_kernel_inhomogeneous(
InternalFaceIntegratorPressure &pressure_m,
InternalFaceIntegratorVelocity &velocity_m,
const ExternalFaceIntegratorPressure &pressure_p,
const ExternalFaceIntegratorVelocity &velocity_p,
const typename InternalFaceIntegratorPressure::value_type c,
const typename InternalFaceIntegratorPressure::value_type rho,
const MaterialIntegrator &c_r,
const MaterialIntegrator &rho_r) const
{
const auto tau_m = 0.5 * rho * c;
const auto gamma_m = 0.5 / (rho * c);
for (unsigned int q : pressure_m.quadrature_point_indices())
{
const auto c_p = c_r.get_value(q);
const auto rho_p = rho_r.get_value(q);
const auto tau_p = 0.5 * rho_p * c_p;
const auto gamma_p = 0.5 / (rho_p * c_p);
const auto tau_sum_inv = 1.0 / (tau_m + tau_p);
const auto gamma_sum_inv = 1.0 / (gamma_m + gamma_p);
const auto n = pressure_m.normal_vector(q);
const auto pm = pressure_m.get_value(q);
const auto um = velocity_m.get_value(q);
const auto pp = pressure_p.get_value(q);
const auto up = velocity_p.get_value(q);
const auto momentum_flux =
pm - tau_m * tau_sum_inv * (pm - pp) +
tau_m * tau_p * tau_sum_inv * (um - up) * n;
velocity_m.submit_value(1.0 / rho * (momentum_flux - pm) * n, q);
const auto mass_flux =
um - gamma_m * gamma_sum_inv * (um - up) +
gamma_m * gamma_p * gamma_sum_inv * (pm - pp) * n;
pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
}
}
template <typename VectorType>
void local_apply_face(
const MatrixFree<dim, Number> &matrix_free,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
matrix_free, true, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_p(
matrix_free, false, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
matrix_free, true, 0, 0, 1);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_p(
matrix_free, false, 0, 0, 1);
MaterialEvaluation material(matrix_free, *material_data);
for (unsigned int face = face_range.first; face < face_range.second;
face++)
{
velocity_m.reinit(face);
velocity_p.reinit(face);
pressure_m.reinit(face);
pressure_p.reinit(face);
pressure_m.gather_evaluate(src, EvaluationFlags::values);
pressure_p.gather_evaluate(src, EvaluationFlags::values);
velocity_m.gather_evaluate(src, EvaluationFlags::values);
velocity_p.gather_evaluate(src, EvaluationFlags::values);
material.reinit_face(face);
evaluate_face_kernel<true>(pressure_m,
velocity_m,
pressure_p,
velocity_p,
material.get_speed_of_sound(),
material.get_density());
pressure_m.integrate_scatter(EvaluationFlags::values, dst);
pressure_p.integrate_scatter(EvaluationFlags::values, dst);
velocity_m.integrate_scatter(EvaluationFlags::values, dst);
velocity_p.integrate_scatter(EvaluationFlags::values, dst);
}
}
template <typename VectorType>
void local_apply_boundary_face_point_to_point(
const MatrixFree<dim, Number> &matrix_free,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
matrix_free, true, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
matrix_free, true, 0, 0, 1);
BCEvaluationP pressure_bc(pressure_m);
BCEvaluationU velocity_bc(velocity_m);
MaterialEvaluation material(matrix_free, *material_data);
auto pressure_r = pressure_r_eval->get_data_accessor();
auto velocity_r = velocity_r_eval->get_data_accessor();
auto c_r = c_r_eval->get_data_accessor();
auto rho_r = rho_r_eval->get_data_accessor();
for (unsigned int face = face_range.first; face < face_range.second;
face++)
{
velocity_m.reinit(face);
pressure_m.reinit(face);
pressure_m.gather_evaluate(src, EvaluationFlags::values);
velocity_m.gather_evaluate(src, EvaluationFlags::values);
if (HelperFunctions::is_non_matching_face(
remote_face_ids, matrix_free.get_boundary_id(face)))
{
velocity_r.reinit(face);
pressure_r.reinit(face);
material.reinit_face(face);
if (material.is_homogeneous())
{
evaluate_face_kernel<false>(pressure_m,
velocity_m,
pressure_r,
velocity_r,
material.get_speed_of_sound(),
material.get_density());
}
else
{
c_r.reinit(face);
rho_r.reinit(face);
evaluate_face_kernel_inhomogeneous(
pressure_m,
velocity_m,
pressure_r,
velocity_r,
material.get_speed_of_sound(),
material.get_density(),
c_r,
rho_r);
}
}
else
{
material.reinit_face(face);
evaluate_face_kernel<false>(pressure_m,
velocity_m,
pressure_bc,
velocity_bc,
material.get_speed_of_sound(),
material.get_density());
}
pressure_m.integrate_scatter(EvaluationFlags::values, dst);
velocity_m.integrate_scatter(EvaluationFlags::values, dst);
}
}
template <typename VectorType>
void local_apply_boundary_face_mortaring(
const MatrixFree<dim, Number> &matrix_free,
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &face_range) const
{
FEFaceEvaluation<dim, -1, 0, 1, Number> pressure_m(
matrix_free, true, 0, 0, 0);
FEFaceEvaluation<dim, -1, 0, dim, Number> velocity_m(
matrix_free, true, 0, 0, 1);
get_thread_safe_fe_face_point_evaluation_object<1>(
thread_local_pressure_m_mortar, 0);
get_thread_safe_fe_face_point_evaluation_object<dim>(
thread_local_velocity_m_mortar, 1);
BCEvaluationP pressure_bc(pressure_m);
BCEvaluationU velocity_bc(velocity_m);
MaterialEvaluation material(matrix_free, *material_data);
auto pressure_r_mortar = pressure_r_eval->get_data_accessor();
auto velocity_r_mortar = velocity_r_eval->get_data_accessor();
auto c_r = c_r_eval->get_data_accessor();
auto rho_r = rho_r_eval->get_data_accessor();
for (unsigned int face = face_range.first; face < face_range.second;
++face)
{
if (HelperFunctions::is_non_matching_face(
remote_face_ids, matrix_free.get_boundary_id(face)))
{
material.reinit_face(face);
pressure_m.reinit(face);
velocity_m.reinit(face);
pressure_m.read_dof_values(src);
velocity_m.read_dof_values(src);
pressure_m.project_to_face(EvaluationFlags::values);
velocity_m.project_to_face(EvaluationFlags::values);
for (unsigned int v = 0;
v < matrix_free.n_active_entries_per_face_batch(face);
++v)
{
constexpr unsigned int n_lanes =
velocity_m_mortar.reinit(face * n_lanes + v);
pressure_m_mortar.reinit(face * n_lanes + v);
velocity_m_mortar.evaluate_in_face(
&velocity_m.get_scratch_data().begin()[0][v],
pressure_m_mortar.evaluate_in_face(
&pressure_m.get_scratch_data().begin()[0][v],
velocity_r_mortar.reinit(face * n_lanes + v);
pressure_r_mortar.reinit(face * n_lanes + v);
if (material.is_homogeneous())
{
evaluate_face_kernel<false>(
pressure_m_mortar,
velocity_m_mortar,
pressure_r_mortar,
velocity_r_mortar,
material.get_speed_of_sound()[v],
material.get_density()[v]);
}
else
{
c_r.reinit(face * n_lanes + v);
rho_r.reinit(face * n_lanes + v);
evaluate_face_kernel_inhomogeneous(
pressure_m_mortar,
velocity_m_mortar,
pressure_r_mortar,
velocity_r_mortar,
material.get_speed_of_sound()[v],
material.get_density()[v],
c_r,
rho_r);
}
velocity_m_mortar.integrate_in_face(
&velocity_m.get_scratch_data().begin()[0][v],
pressure_m_mortar.integrate_in_face(
&pressure_m.get_scratch_data().begin()[0][v],
}
pressure_m.collect_from_face(EvaluationFlags::values);
velocity_m.collect_from_face(EvaluationFlags::values);
pressure_m.distribute_local_to_global(dst);
velocity_m.distribute_local_to_global(dst);
}
else
{
velocity_m.reinit(face);
pressure_m.reinit(face);
pressure_m.gather_evaluate(src, EvaluationFlags::values);
velocity_m.gather_evaluate(src, EvaluationFlags::values);
material.reinit_face(face);
evaluate_face_kernel<false>(pressure_m,
velocity_m,
pressure_bc,
velocity_bc,
material.get_speed_of_sound(),
material.get_density());
pressure_m.integrate_scatter(EvaluationFlags::values, dst);
velocity_m.integrate_scatter(EvaluationFlags::values, dst);
}
}
}
const MatrixFree<dim, Number> &matrix_free;
const std::shared_ptr<CellwiseMaterialData<Number>> material_data;
const std::set<types::boundary_id> remote_face_ids;
const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
pressure_r_eval;
const std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
velocity_r_eval;
const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
c_r_eval;
const std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
rho_r_eval;
const std::shared_ptr<NonMatching::MappingInfo<dim, dim, Number>>
nm_mapping_info;
std::unique_ptr<FEFacePointEvaluation<1, dim, dim, Number>>>
thread_local_pressure_m_mortar;
std::unique_ptr<FEFacePointEvaluation<dim, dim, dim, Number>>>
thread_local_velocity_m_mortar;
template <int n_components>
get_thread_safe_fe_face_point_evaluation_object(
&fe_face_point_eval_thread_local,
unsigned int fist_selected_comp) const
{
if (fe_face_point_eval_thread_local.get() == nullptr)
{
fe_face_point_eval_thread_local = std::make_unique<
*nm_mapping_info,
matrix_free.get_dof_handler().get_fe(),
true,
fist_selected_comp);
}
return *fe_face_point_eval_thread_local.get();
}
};
template <int dim, typename Number>
class InverseMassOperator
{
public:
InverseMassOperator(const MatrixFree<dim, Number> &matrix_free)
: matrix_free(matrix_free)
{}
template <typename VectorType>
void apply(VectorType &dst, const VectorType &src) const
{
dst.zero_out_ghost_values();
matrix_free.cell_loop(&InverseMassOperator::local_apply_cell<VectorType>,
this,
dst,
src);
}
private:
template <typename VectorType>
void local_apply_cell(
VectorType &dst,
const VectorType &src,
const std::pair<unsigned int, unsigned int> &cell_range) const
{
FEEvaluation<dim, -1, 0, dim + 1, Number> phi(mf);
minv(phi);
for (unsigned int cell = cell_range.first; cell < cell_range.second;
++cell)
{
phi.reinit(cell);
phi.read_dof_values(src);
minv.apply(phi.begin_dof_values(), phi.begin_dof_values());
phi.set_dof_values(dst);
}
}
const MatrixFree<dim, Number> &matrix_free;
};
template <int dim, typename Number, typename remote_value_type>
class RungeKutta2
{
public:
RungeKutta2(
const std::shared_ptr<InverseMassOperator<dim, Number>>
inverse_mass_operator,
const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
acoustic_operator)
: inverse_mass_operator(inverse_mass_operator)
, acoustic_operator(acoustic_operator)
{}
void run(const MatrixFree<dim, Number> &matrix_free,
const double cr,
const double end_time,
const double speed_of_sound,
const Function<dim> &initial_condition,
const std::string &vtk_prefix)
{
const auto &dof_handler = matrix_free.get_dof_handler();
const auto &mapping = *matrix_free.get_mapping_info().mapping;
const auto degree = dof_handler.get_fe().degree;
VectorType solution;
matrix_free.initialize_dof_vector(solution);
VectorType solution_temp;
matrix_free.initialize_dof_vector(solution_temp);
HelperFunctions::set_initial_condition(matrix_free,
initial_condition,
solution);
double h_local_min = std::numeric_limits<double>::max();
for (const auto &cell : dof_handler.active_cell_iterators())
h_local_min =
std::min(h_local_min,
(cell->vertex(1) - cell->vertex(0)).norm_square());
h_local_min = std::sqrt(h_local_min);
const double h_min =
Utilities::MPI::min(h_local_min, dof_handler.get_communicator());
const double dt =
cr * HelperFunctions::compute_dt_cfl(h_min, degree, speed_of_sound);
double time = 0.0;
unsigned int timestep = 0;
while (time < end_time)
{
HelperFunctions::write_vtu(solution,
matrix_free.get_dof_handler(),
mapping,
degree,
"step_89-" + vtk_prefix +
std::to_string(timestep));
std::swap(solution, solution_temp);
time += dt;
timestep++;
perform_time_step(dt, solution, solution_temp);
}
}
private:
void
perform_time_step(const double dt, VectorType &dst, const VectorType &src)
{
VectorType k1 = src;
evaluate_stage(k1, src);
k1.sadd(0.5 * dt, 1.0, src);
evaluate_stage(dst, k1);
dst.sadd(dt, 1.0, src);
}
void evaluate_stage(VectorType &dst, const VectorType &src)
{
acoustic_operator->evaluate(dst, src);
dst *= -1.0;
inverse_mass_operator->apply(dst, dst);
}
const std::shared_ptr<InverseMassOperator<dim, Number>>
inverse_mass_operator;
const std::shared_ptr<AcousticOperator<dim, Number, remote_value_type>>
acoustic_operator;
};
template <int dim>
void build_non_matching_triangulation(
std::set<types::boundary_id> &non_matching_faces,
const unsigned int refinements)
{
const double length = 1.0;
const types::boundary_id non_matching_id_left = 98;
const types::boundary_id non_matching_id_right = 99;
non_matching_faces.insert(non_matching_id_left);
non_matching_faces.insert(non_matching_id_right);
Triangulation<dim> tria_left;
const unsigned int subdiv_left = 11;
{subdiv_left, 2 * subdiv_left},
{0.0, 0.0},
{0.525 * length, length});
for (const auto &cell : tria_left.active_cell_iterators())
cell->set_material_id(0);
for (const auto &face : tria_left.active_face_iterators())
if (face->at_boundary())
{
face->set_boundary_id(0);
if (face->center()[0] > 0.525 * length - 1e-6)
face->set_boundary_id(non_matching_id_left);
}
Triangulation<dim> tria_right;
const unsigned int subdiv_right = 4;
{subdiv_right, 2 * subdiv_right},
{0.525 * length, 0.0},
{length, length});
for (const auto &cell : tria_right.active_cell_iterators())
cell->set_material_id(1);
for (const auto &face : tria_right.active_face_iterators())
if (face->at_boundary())
{
face->set_boundary_id(0);
if (face->center()[0] < 0.525 * length + 1e-6)
face->set_boundary_id(non_matching_id_right);
}
tria_right,
tria,
/*tolerance*/ 0.,
/*copy_manifold_ids*/ false,
/*copy_boundary_ids*/ true);
tria.refine_global(refinements);
}
template <int dim, typename Number>
void run_with_point_to_point_interpolation(
const MatrixFree<dim, Number> &matrix_free,
const std::set<types::boundary_id> &non_matching_faces,
const std::map<types::material_id, std::pair<double, double>> &materials,
const double end_time,
const Function<dim> &initial_condition,
const std::string &vtk_prefix)
{
const auto &dof_handler = matrix_free.get_dof_handler();
const auto &tria = dof_handler.get_triangulation();
std::vector<
std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
non_matching_faces_marked_vertices;
for (const auto &nm_face : non_matching_faces)
{
auto marked_vertices = [&]() {
std::vector<bool> mask(tria.n_vertices(), true);
for (const auto &cell : tria.active_cell_iterators())
for (auto const &f : cell->face_indices())
if (cell->face(f)->at_boundary() &&
cell->face(f)->boundary_id() == nm_face)
for (const auto v : cell->vertex_indices())
mask[cell->vertex_index(v)] = false;
return mask;
};
non_matching_faces_marked_vertices.emplace_back(
std::make_pair(nm_face, marked_vertices));
}
auto remote_communicator =
matrix_free, non_matching_faces_marked_vertices);
using remote_value_type = VectorizedArray<Number>;
const auto pressure_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 0);
const auto velocity_r =
std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 1);
const auto material_data =
std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
const auto c_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator,
matrix_free.get_dof_handler().get_triangulation(),
/*first_selected_component*/ 0);
const auto rho_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator,
matrix_free.get_dof_handler().get_triangulation(),
/*first_selected_component*/ 0);
if (!material_data->is_homogeneous())
{
matrix_free.get_dof_handler().get_triangulation().n_active_cells());
matrix_free.get_dof_handler().get_triangulation().n_active_cells());
for (const auto &cell : matrix_free.get_dof_handler()
.get_triangulation()
.active_cell_iterators())
{
c[cell->active_cell_index()] =
materials.at(cell->material_id()).first;
rho[cell->active_cell_index()] =
materials.at(cell->material_id()).second;
}
c_r->gather_evaluate(c, EvaluationFlags::values);
rho_r->gather_evaluate(rho, EvaluationFlags::values);
}
const auto inverse_mass_operator =
std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
const auto acoustic_operator =
std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
matrix_free,
material_data,
non_matching_faces,
pressure_r,
velocity_r,
c_r,
rho_r);
double speed_of_sound_max = 0.0;
for (const auto &mat : materials)
speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
RungeKutta2<dim, Number, remote_value_type> time_integrator(
inverse_mass_operator, acoustic_operator);
time_integrator.run(matrix_free,
/*Cr*/ 0.2,
end_time,
speed_of_sound_max,
initial_condition,
vtk_prefix);
}
template <int dim, typename Number>
void run_with_nitsche_type_mortaring(
const MatrixFree<dim, Number> &matrix_free,
const std::set<types::boundary_id> &non_matching_faces,
const std::map<types::material_id, std::pair<double, double>> &materials,
const double end_time,
const Function<dim> &initial_condition,
const std::string &vtk_prefix)
{
#ifndef DEAL_II_WITH_CGAL
(void)matrix_free;
(void)non_matching_faces;
(void)materials;
(void)end_time;
(void)initial_condition;
(void)vtk_prefix;
std::cout, (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0));
pcout << "In this function, mortars are computed using CGAL. "
"Configure deal.II with DEAL_II_WITH_CGAL to run this function.\n";
return;
#else
const auto &dof_handler = matrix_free.get_dof_handler();
const auto &tria = dof_handler.get_triangulation();
const auto &mapping = *matrix_free.get_mapping_info().mapping;
const auto n_quadrature_pnts = matrix_free.get_quadrature().size();
std::vector<
std::pair<types::boundary_id, std::function<std::vector<bool>()>>>
non_matching_faces_marked_vertices;
for (const auto &nm_face : non_matching_faces)
{
auto marked_vertices = [&]() {
std::vector<bool> mask(tria.n_vertices(), true);
for (const auto &cell : tria.active_cell_iterators())
for (auto const &f : cell->face_indices())
if (cell->face(f)->at_boundary() &&
cell->face(f)->boundary_id() == nm_face)
for (const auto v : cell->vertex_indices())
mask[cell->vertex_index(v)] = false;
return mask;
};
non_matching_faces_marked_vertices.emplace_back(
std::make_pair(nm_face, marked_vertices));
}
additional_data;
additional_data.use_global_weights = true;
auto nm_mapping_info =
std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
additional_data);
auto remote_communicator =
matrix_free,
non_matching_faces_marked_vertices,
n_quadrature_pnts,
0,
nm_mapping_info.get());
using remote_value_type = Number;
const auto pressure_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 0);
const auto velocity_r =
std::make_shared<FERemoteEvaluation<dim, dim, remote_value_type>>(
remote_communicator, dof_handler, /*first_selected_component*/ 1);
const auto material_data =
std::make_shared<CellwiseMaterialData<Number>>(matrix_free, materials);
const auto c_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator,
matrix_free.get_dof_handler().get_triangulation(),
/*first_selected_component*/ 0);
const auto rho_r =
std::make_shared<FERemoteEvaluation<dim, 1, remote_value_type>>(
remote_communicator,
matrix_free.get_dof_handler().get_triangulation(),
/*first_selected_component*/ 0);
if (!material_data->is_homogeneous())
{
matrix_free.get_dof_handler().get_triangulation().n_active_cells());
matrix_free.get_dof_handler().get_triangulation().n_active_cells());
for (const auto &cell : matrix_free.get_dof_handler()
.get_triangulation()
.active_cell_iterators())
{
c[cell->active_cell_index()] =
materials.at(cell->material_id()).first;
rho[cell->active_cell_index()] =
materials.at(cell->material_id()).second;
}
c_r->gather_evaluate(c, EvaluationFlags::values);
rho_r->gather_evaluate(rho, EvaluationFlags::values);
}
const auto inverse_mass_operator =
std::make_shared<InverseMassOperator<dim, Number>>(matrix_free);
const auto acoustic_operator =
std::make_shared<AcousticOperator<dim, Number, remote_value_type>>(
matrix_free,
material_data,
non_matching_faces,
pressure_r,
velocity_r,
c_r,
rho_r,
nm_mapping_info);
double speed_of_sound_max = 0.0;
for (const auto &mat : materials)
speed_of_sound_max = std::max(speed_of_sound_max, mat.second.first);
RungeKutta2<dim, Number, remote_value_type> time_integrator(
inverse_mass_operator, acoustic_operator);
time_integrator.run(matrix_free,
/*Cr*/ 0.2,
end_time,
speed_of_sound_max,
initial_condition,
vtk_prefix);
#endif
}
} // namespace Step89
int main(int argc, char *argv[])
{
using namespace dealii;
constexpr int dim = 2;
using Number = double;
std::cout.precision(5);
ConditionalOStream pcout(std::cout,
0));
const unsigned int refinements = 1;
const unsigned int degree = 3;
#ifdef DEAL_II_WITH_P4EST
#else
#endif
pcout << "Create non-matching grid..." << std::endl;
std::set<types::boundary_id> non_matching_faces;
Step89::build_non_matching_triangulation(tria,
non_matching_faces,
refinements);
pcout << " - Refinement level: " << refinements << std::endl;
pcout << " - Number of cells: " << tria.n_cells() << std::endl;
pcout << "Create DoFHandler..." << std::endl;
DoFHandler<dim> dof_handler(tria);
dof_handler.distribute_dofs(FESystem<dim>(FE_DGQ<dim>(degree), dim + 1));
pcout << " - Number of DoFs: " << dof_handler.n_dofs() << std::endl;
constraints.close();
pcout << "Set up MatrixFree..." << std::endl;
matrix_free.reinit(
MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
pcout << "Run vibrating membrane test case..." << std::endl;
const double modes = 10.0;
std::map<types::material_id, std::pair<double, double>> homogeneous_material;
homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0);
const auto initial_solution_membrane =
Step89::InitialConditionVibratingMembrane<dim>(modes);
pcout << " - Point-to-point interpolation: " << std::endl;
Step89::run_with_point_to_point_interpolation(
matrix_free,
non_matching_faces,
homogeneous_material,
8.0 * initial_solution_membrane.get_period_duration(
homogeneous_material.begin()->second.first),
initial_solution_membrane,
"vm-p2p");
pcout << " - Nitsche-type mortaring: " << std::endl;
Step89::run_with_nitsche_type_mortaring(
matrix_free,
non_matching_faces,
homogeneous_material,
8.0 * initial_solution_membrane.get_period_duration(
homogeneous_material.begin()->second.first),
initial_solution_membrane,
"vm-nitsche");
pcout << "Run test case with in-homogeneous material..." << std::endl;
std::map<types::material_id, std::pair<double, double>>
inhomogeneous_material;
inhomogeneous_material[0] = std::make_pair(1.0, 1.0);
inhomogeneous_material[1] = std::make_pair(3.0, 1.0);
Step89::run_with_nitsche_type_mortaring(matrix_free,
non_matching_faces,
inhomogeneous_material,
/*runtime*/ 0.3,
Step89::GaussPulse<dim>(0.3, 0.5),
"inhomogeneous");
return 0;
}
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
void set_flags(const FlagType &flags)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition data_out.cc:1062
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
const Triangulation< dim, spacedim > & get_triangulation() const
MPI_Comm get_communicator() const
value_type get_value(const unsigned int q_point) const
void integrate_in_face(ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &integration_flags, const bool sum_into_values=false)
void evaluate_in_face(const ScalarNumber *face_dof_values, const EvaluationFlags::EvaluationFlags &evaluation_flags)
void reinit(const unsigned int cell_index, const unsigned int face_number)
const unsigned int degree
Definition fe_data.h:452
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
types::boundary_id get_boundary_id(const unsigned int face_batch_index) const
const Quadrature< dim > & get_quadrature(const unsigned int quad_index=0, const unsigned int hp_active_fe_index=0) const
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
unsigned int n_active_entries_per_face_batch(const unsigned int face_batch_index) const
void write_vtu(const std::vector< Patch< dim, spacedim > > &patches, const std::vector< std::string > &data_names, const std::vector< std::tuple< unsigned int, unsigned int, std::string, DataComponentInterpretation::DataComponentInterpretation > > &nonscalar_data_ranges, const VtkFlags &flags, std::ostream &out)
spacedim & mapping
const Triangulation< dim, spacedim > & tria
spacedim const Point< spacedim > const std::vector< bool > & marked_vertices
Definition grid_tools.h:991
if(marked_vertices.size() !=0) for(auto it
for(unsigned int j=best_vertex+1;j< vertices.size();++j) if(vertices_to_use[j])
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
long double gamma(const unsigned int n)
unsigned int material_id
Definition types.h:167
unsigned int boundary_id
Definition types.h:144
UpdateFlags mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_boundary_faces