]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add tutorial on Nitsche-type mortaring 16299/head
authorJohannes Heinz <43043310+jh66637@users.noreply.github.com>
Fri, 2 Feb 2024 15:32:21 +0000 (16:32 +0100)
committerJohannes Heinz <43043310+jh66637@users.noreply.github.com>
Fri, 22 Mar 2024 18:50:31 +0000 (19:50 +0100)
doc/doxygen/references.bib
doc/news/changes/major/20231127HeinzBergbauerFederMunch [new file with mode: 0644]
examples/step-89/CMakeLists.txt [new file with mode: 0644]
examples/step-89/doc/builds-on [new file with mode: 0644]
examples/step-89/doc/intro.dox [new file with mode: 0644]
examples/step-89/doc/kind [new file with mode: 0644]
examples/step-89/doc/results.dox [new file with mode: 0644]
examples/step-89/doc/tooltip [new file with mode: 0644]
examples/step-89/step-89.cc [new file with mode: 0644]

index b82b2a5d7ba49f396ef4564ce4a5428a155685c2..85c081c07751d4c2db88adb76e609cd0be106aec 100644 (file)
 }
 
 
+%-------------------------------------------------------------------------------
+% Step 89
+%-------------------------------------------------------------------------------
+
+@article{hochbruck2014efficient,
+  author    = {Marlis Hochbruck and Tomislav Pa{\v{z}}ur and Andreas Schulz and Ekkachai Thawinan and Christian Wieners},
+  journal   = {{ZAMM} - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik},
+  title     = {Efficient time integration for discontinuous {G}alerkin approximations of linear wave equations},
+  year      = {2014},
+  month     = {sep},
+  number    = {3},
+  pages     = {237--259},
+  volume    = {95},
+  doi       = {10.1002/zamm.201300306},
+  publisher = {Wiley},
+}
+
+@article{arnold2002unified,
+  author  = {Arnold, D. and Brezzi, F. and Cockburn, B. and Marini, L.},
+  journal = {SIAM Journal on Numerical Analysis},
+  title   = {Unified analysis of discontinuous {Galerkin} methods for elliptic problems},
+  year    = {2002},
+  number  = {5},
+  pages   = {1749-1779},
+  volume  = {39},
+  doi     = {10.1137/S0036142901384162},
+  eprint  = {https://doi.org/10.1137/S0036142901384162},
+}
+
+@article{nguyen2011high,
+  author   = {N.C. Nguyen and J. Peraire and B. Cockburn},
+  journal  = {Journal of Computational Physics},
+  title    = {High-order implicit hybridizable discontinuous {G}alerkin methods for acoustics and elastodynamics},
+  year     = {2011},
+  issn     = {0021-9991},
+  number   = {10},
+  pages    = {3695 - 3718},
+  volume   = {230},
+  doi      = {https://doi.org/10.1016/j.jcp.2011.01.035},
+}
+
+@article{bangerth2010adaptive,
+  author    = {W. Bangerth and M. Geiger and R. Rannacher},
+  journal   = {Computational Methods in Applied Mathematics},
+  title     = {Adaptive {G}alerkin finite element methods for the wave equation},
+  year      = {2010},
+  number    = {1},
+  pages     = {3--48},
+  volume    = {10},
+  doi       = {10.2478/cmam-2010-0001},
+  publisher = {Walter de Gruyter {GmbH}},
+}
+
+@article{duerrwaechter2021an,
+  author    = {Jakob Dürrwächter and Marius Kurz and Patrick Kopper and Daniel Kempf and Claus-Dieter Munz and Andrea Beck},
+  journal   = {Computers \& Fluids},
+  title     = {An efficient sliding mesh interface method for high-order discontinuous {G}alerkin schemes},
+  year      = {2021},
+  month     = {mar},
+  pages     = {104825},
+  volume    = {217},
+  doi       = {10.1016/j.compfluid.2020.104825},
+  publisher = {Elsevier {BV}},
+}
+
 
 %-------------------------------------------------------------------------------
 % References used elsewhere
diff --git a/doc/news/changes/major/20231127HeinzBergbauerFederMunch b/doc/news/changes/major/20231127HeinzBergbauerFederMunch
new file mode 100644 (file)
index 0000000..f8f70b4
--- /dev/null
@@ -0,0 +1,6 @@
+New: The new tutorial step-89 presents the use of FERemoteEvaluation during matrix-free operator
+evaluation for non-matching and Chimera methods. The acoustic conservation equations are solved
+using Nitsche-type mortaring and point-to-point interpolation to demonstrate that
+a simple point-to-point interpolation approach is sometimes not sufficient.
+<br>
+(Johannes Heinz, Marco Feder, Peter Munch, 2023/11/27)
diff --git a/examples/step-89/CMakeLists.txt b/examples/step-89/CMakeLists.txt
new file mode 100644 (file)
index 0000000..9795f5d
--- /dev/null
@@ -0,0 +1,55 @@
+##
+#  CMake script
+##
+
+# Set the name of the project and target:
+set(TARGET "step-89")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#  file(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#  file(GLOB_RECURSE TARGET_INC  "include/*.h")
+#  set(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC})
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+set(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+cmake_minimum_required(VERSION 3.13.4)
+
+find_package(deal.II 9.6.0
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+if(NOT ${deal.II_FOUND})
+  message(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+endif()
+
+#
+# Are all dependencies fulfilled?
+#
+if(NOT DEAL_II_WITH_MPI) # keep in one line
+  message(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following option:
+    DEAL_II_WITH_MPI = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options:
+    DEAL_II_WITH_MPI = ${DEAL_II_WITH_MPI}
+This conflicts with the requirements."
+    )
+endif()
+
+
+deal_ii_initialize_cached_variables()
+set(CLEAN_UP_FILES *.log *.gmv *.gnuplot *.gpl *.eps *.pov *.ucd *.d2 *.vtu *.pvtu)
+project(${TARGET})
+deal_ii_invoke_autopilot()
+-
diff --git a/examples/step-89/doc/builds-on b/examples/step-89/doc/builds-on
new file mode 100644 (file)
index 0000000..809da9b
--- /dev/null
@@ -0,0 +1 @@
+step-67 step-87
diff --git a/examples/step-89/doc/intro.dox b/examples/step-89/doc/intro.dox
new file mode 100644 (file)
index 0000000..666eef9
--- /dev/null
@@ -0,0 +1,221 @@
+<br>
+
+<i>
+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 @cite heinz2022high.
+
+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].
+</i>
+
+@dealiiTutorialDOI{10.5281/zenodo.10033975,https://zenodo.org/badge/DOI/10.5281/zenodo.10033975.svg}
+
+<h1>Introduction</h1>
+
+This tutorial presents one way how to apply non-matching and/or Chimera methods
+within matrix-free loops in deal.II.
+We are following @cite heinz2022high 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 @cite heinz2022high.
+
+<h3>%Acoustic conservation equations</h3>
+
+Acoustic conservation equations are used to describe linear wave propagation.
+The set of equations consists of the conservation of mass and momentum
+@f[
+  \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}.
+@f]
+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.
+For the discretization 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 @cite arnold2002unified. I.e., numerical fluxes are not only used on inner
+element faces but also as non-matching coupling conditions.
+
+The discretized equations read
+@f[
+  \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},
+@f]
+where $\mathbf{w}_h$ and $q_h$ are test functions. The numerical fluxes are
+defined as follows @cite hochbruck2014efficient
+@f[
+  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},
+@f]
+with the penalty parameters $\tau=\frac{\rho c}{2}$ and $\gamma=\frac{1}{2\rho c}$.
+$[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^+$)
+@f[
+  p_h^*=\average{p_h}+\tau\jump{\mathbf{u}_h},\\
+  \mathbf{u}_h^*=\average{\mathbf{u}_h}+\gamma\jump{p_h}.
+@f]
+$\average{a}=\frac{a^- + a^+}{2}$ denots the averaging operator.
+
+<h3>%Non-matching discretizations</h3>
+
+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).
+
+Considering sliding rotating interfaces @cite duerrwaechter2021an also requires the ability to
+handle non-matching discretizations: A cylindrical mesh is embedded in a surrounding mesh with a
+cylindrical hole. Every time step, the cylinder rotates while the outer mesh remains at the
+same position. In this situation it is hardly possible to construct a conforming mesh in every
+time step.
+
+Besides this, non-matching methods can be extended to Chimera methods (elements overlap).
+Chimera methods can help to reduce the pressure on mesh generation tools since different regions
+of a mesh (that may overlap) can be considered 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.
+
+<h4>%Point-to-point interpolation</h4>
+
+@image html https://www.dealii.org/images/steps/developer/step_89_intro_point_to_point.svg "" width=25%
+
+Point-to-point interpolation is a naive approach. The points in which values/gradients are
+queried in the coupling terms are defined by the quadrature points on the element face of
+element $K^-$. As it can be seen from the picture this approach might be subject to aliasing
+in some cases. In the picture, information from element $K_1^+$ is completely neglected.
+
+<h4>%Nitsche-type mortaring</h4>
+
+@image html https://www.dealii.org/images/steps/developer/step_89_intro_mortaring.svg "" width=25%
+
+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. This way, the numerical integration is exact as long
+as a sufficient number of integration points is used.
+In this tutorial, the intersections are computed using `CGAL`. Therefore, `deal.II` has
+to be configured with `DEAL_II_WITH_CGAL` for the Nitsche-type mortaring implementation.
+
+<h3>%FERemoteEvaluation</h3>
+
+FERemoteEvaluation is a wrapper class which provides a similar interface to, e.g., FEEvaluation 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 on 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 <i>before</i>
+the loop that passes through 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:
+@code
+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>({}, face_function, {}, dst, src);
+@endcode
+
+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.
+
+@code
+// 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);
+@endcode
+The object @c 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
+<i>before</i> the actual matrix-free loop. FERemoteCommunicator, as well as FERemoteEvaluation behaves differently for
+the given template parameter @c value_type. If we want to access values at arbitrary points (e.g. in combination with
+@c FEPointEvaluation) @c value_type=Number. If the values are defined at quadrature points of a @c FEEvaluation object
+it is possible to get the values at the quadrature points of <i>batches</i> and @c value_type=VectorizedArray<Number>.
+
+<h3>Overview</h3>
+
+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. @cite nguyen2011high.
+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
+
+@f[
+  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),
+@f]
+
+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 @cite bangerth2010adaptive.
diff --git a/examples/step-89/doc/kind b/examples/step-89/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-89/doc/results.dox b/examples/step-89/doc/results.dox
new file mode 100644 (file)
index 0000000..141fba9
--- /dev/null
@@ -0,0 +1,80 @@
+<h1>Results</h1>
+
+<h3>Vibrating membrane: Point-to-point interpolation vs. Nitsche-type mortaring</h3>
+
+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.
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+        @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_velocity_Y.png "" width=60%
+    </td>
+    <td>
+        @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_point_to_point_velocity_Y.png "" width=60%
+    </td>
+  </tr>
+</table>
+
+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).
+
+<table align="center" class="doxtable">
+  <tr>
+    <td>
+        @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_pressure.svg "" width=100%
+    </td>
+    <td>
+        @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_velocity_Y.svg "" width=100%
+    </td>
+  </tr>
+</table>
+
+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 @cite heinz2022high.
+
+<h3>Wave propagation through in-homogeneous fluid</h3>
+
+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.
+
+@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_pressure.png "" width=30%
+
+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. @c 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).
+
+@image html https://www.dealii.org/images/steps/developer/step_89_inhomogenous_test_case_conforming_vs_nonmatching.svg "" width=60%
+
+The results, obtained with the non-matching discretization is in good agreement with
+the reference result.
+
+<h3>Possibilities for extensions</h3>
+
+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 <i>not</i> 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.
diff --git a/examples/step-89/doc/tooltip b/examples/step-89/doc/tooltip
new file mode 100644 (file)
index 0000000..10f51f3
--- /dev/null
@@ -0,0 +1 @@
+Matrix-free operator evaluation for non-matching and Chimera methods with application to acoustic conservation equations.
diff --git a/examples/step-89/step-89.cc b/examples/step-89/step-89.cc
new file mode 100644 (file)
index 0000000..b39ffb8
--- /dev/null
@@ -0,0 +1,1796 @@
+/* ---------------------------------------------------------------------
+ *
+ * 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
+ */
+
+// @sect3{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;
+
+  // @sect3{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;
+
+    // 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);
+  }
+
+  // @sect3{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;
+  }
+
+  // @sect3{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 @c 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);
+    }
+
+    // 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;
+      DataOutBase::VtkFlags flags;
+      flags.write_higher_order_cells = true;
+      data_out.set_flags(flags);
+
+      std::vector<DataComponentInterpretation::DataComponentInterpretation>
+        interpretation(
+          dim + 1, DataComponentInterpretation::component_is_part_of_vector);
+      std::vector<std::string> names(dim + 1, "U");
+
+      interpretation[0] = DataComponentInterpretation::component_is_scalar;
+      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
+
+  //@sect3{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)
+        {
+          // 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
+          // @c 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.
+    AlignedVector<VectorizedArray<Number>> speed_of_sound;
+    AlignedVector<VectorizedArray<Number>> density;
+
+    // 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
+  // @c MaterialEvaluation is used. Similar to @c FEEvaluation, every thread
+  // creates its own instance and thus, there are no race conditions. For
+  // in-homogeneous materials, a @c reinit_cell() or @c reinit_face() function
+  // is used to set the correct material at the current cell batch. In the
+  // homogeneous case the @c _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;
+    VectorizedArray<Number> density;
+  };
+
+
+  //@sect3{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;
+  };
+
+  //@sect3{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,
+      std::shared_ptr<FERemoteEvaluation<dim, 1, remote_value_type>>
+        pressure_r_eval,
+      std::shared_ptr<FERemoteEvaluation<dim, dim, remote_value_type>>
+        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 c_r_eval and @c 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,
+            MatrixFree<dim, Number>::DataAccessOnFaces::values,
+            MatrixFree<dim, Number>::DataAccessOnFaces::values);
+        }
+      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,
+            MatrixFree<dim, Number>::DataAccessOnFaces::values,
+            MatrixFree<dim, Number>::DataAccessOnFaces::values);
+        }
+    }
+
+  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 @c 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
+    {
+      // 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);
+        }
+    }
+
+
+    //@sect4{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 @c 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 @c 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);
+        }
+    }
+
+    //@sect4{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 @c 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 @c get_thread_safe_fe_face_point_evaluation_object().
+      FEFacePointEvaluation<1, dim, dim, Number> &pressure_m_mortar =
+        get_thread_safe_fe_face_point_evaluation_object<1>(
+          thread_local_pressure_m_mortar, 0);
+      FEFacePointEvaluation<dim, dim, dim, Number> &velocity_m_mortar =
+        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 =
+                    VectorizedArray<Number>::size();
+                  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],
+                    EvaluationFlags::values);
+
+                  pressure_m_mortar.evaluate_in_face(
+                    &pressure_m.get_scratch_data().begin()[0][v],
+                    EvaluationFlags::values);
+
+                  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 @c
+                      // v we call @c 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],
+                    EvaluationFlags::values);
+
+                  pressure_m_mortar.integrate_in_face(
+                    &pressure_m.get_scratch_data().begin()[0][v],
+                    EvaluationFlags::values);
+                }
+
+              // Collect the contributions from the face DoFs to
+              // the internal cell DoFs to be able to use the
+              // member function @c 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 @c 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.
+    mutable Threads::ThreadLocalStorage<
+      std::unique_ptr<FEFacePointEvaluation<1, dim, dim, Number>>>
+      thread_local_pressure_m_mortar;
+
+    mutable Threads::ThreadLocalStorage<
+      std::unique_ptr<FEFacePointEvaluation<dim, dim, dim, Number>>>
+      thread_local_velocity_m_mortar;
+
+    // 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>
+    FEFacePointEvaluation<n_components, dim, dim, Number> &
+    get_thread_safe_fe_face_point_evaluation_object(
+      Threads::ThreadLocalStorage<
+        std::unique_ptr<FEFacePointEvaluation<n_components, dim, dim, Number>>>
+                  &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<
+            FEFacePointEvaluation<n_components, dim, dim, Number>>(
+            *nm_mapping_info,
+            matrix_free.get_dof_handler().get_fe(),
+            true,
+            fist_selected_comp);
+        }
+      return *fe_face_point_eval_thread_local.get();
+    }
+  };
+
+  //@sect3{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:
+    // Apply the inverse mass operator onto every cell batch.
+    template <typename VectorType>
+    void local_apply_cell(
+      const MatrixFree<dim, Number>               &mf,
+      VectorType                                  &dst,
+      const VectorType                            &src,
+      const std::pair<unsigned int, unsigned int> &cell_range) const
+    {
+      FEEvaluation<dim, -1, 0, dim + 1, Number> phi(mf);
+      MatrixFreeOperators::CellwiseInverseMassMatrix<dim, -1, dim + 1, Number>
+        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;
+  };
+
+  //@sect3{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
+  {
+    using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  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);
+
+      // 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());
+
+      // 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;
+  };
+
+
+  // @sect3{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(
+    Triangulation<dim>           &tria,
+    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;
+    GridGenerator::subdivided_hyper_rectangle(tria_left,
+                                              {subdiv_left, 2 * subdiv_left},
+                                              {0.0, 0.0},
+                                              {0.525 * length, length});
+
+    // 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;
+    GridGenerator::subdivided_hyper_rectangle(tria_right,
+                                              {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
+    // @c merge_triangulations().
+    GridGenerator::merge_triangulations(tria_left,
+                                        tria_right,
+                                        tria,
+                                        /*tolerance*/ 0.,
+                                        /*copy_manifold_ids*/ false,
+                                        /*copy_boundary_ids*/ true);
+    tria.refine_global(refinements);
+  }
+
+  // @sect3{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();
+
+    // 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. @c 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 @cite duerrwaechter2021an.
+    //
+    // For the standard case of point to point-to-point interpolation without
+    // any heuristic we make use of the utility function
+    // @c 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
+          // @c 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 =
+      Utilities::compute_remote_communicator_faces_point_to_point_interpolation(
+        matrix_free, non_matching_faces_marked_vertices);
+
+    // We are using point-to-point interpolation and can therefore
+    // easily access the corresponding data at face batches. This
+    // is why we use a @c VectorizedArray as @c 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.
+        Vector<Number> c(
+          matrix_free.get_dof_handler().get_triangulation().n_active_cells());
+        Vector<Number> rho(
+          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;
+          }
+
+        // Materials do not change during the simulation, therefore
+        // there is no need to precompute the values after
+        // the first @c 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);
+
+    // 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);
+  }
+
+  // @sect3{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;
+
+    ConditionalOStream pcout(
+      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();
+
+    // In case of Nitsche-type mortaring a vector of pairs with cell iterator
+    // and face number is needed as communication object.
+    // @c 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
+    // @c 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
+          // @c 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 @c
+    // use_global_weights.
+    typename NonMatching::MappingInfo<dim, dim, Number>::AdditionalData
+      additional_data;
+    additional_data.use_global_weights = true;
+
+    // Set up NonMatching::MappingInfo with needed update flags and
+    // @c additional_data.
+    auto nm_mapping_info =
+      std::make_shared<NonMatching::MappingInfo<dim, dim, Number>>(
+        mapping,
+        update_values | update_JxW_values | update_normal_vectors |
+          update_quadrature_points,
+        additional_data);
+
+    auto remote_communicator =
+      Utilities::compute_remote_communicator_faces_nitsche_type_mortaring(
+        matrix_free,
+        non_matching_faces_marked_vertices,
+        n_quadrature_pnts,
+        0,
+        nm_mapping_info.get());
+
+    // 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 @c 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())
+      {
+        Vector<Number> c(
+          matrix_free.get_dof_handler().get_triangulation().n_active_cells());
+        Vector<Number> rho(
+          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
+
+
+// @sect3{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;
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv);
+  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
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+#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;
+
+  // 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;
+
+  AffineConstraints<Number> constraints;
+  constraints.close();
+
+  pcout << "Set up MatrixFree..." << std::endl;
+  typename MatrixFree<dim, Number>::AdditionalData data;
+  data.mapping_update_flags             = update_gradients | update_values;
+  data.mapping_update_flags_inner_faces = update_values;
+  data.mapping_update_flags_boundary_faces =
+    update_quadrature_points | update_values;
+
+  MatrixFree<dim, Number> matrix_free;
+  matrix_free.reinit(
+    MappingQ1<dim>(), dof_handler, constraints, QGauss<dim>(degree + 1), data);
+
+
+  //@sect4{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;
+  // 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");
+
+  //@sect4{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;
+}

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.