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 [112].
Johannes Heinz was supported by the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014-2020) under the Marie Skłodowská–Curie Grant Agreement No. [812719].
Note
If you use this program as a basis for your own work, please consider citing it in your list of references. The initial version of this work was contributed to the deal.II project by the authors listed in the following citation:
Introduction
This tutorial presents one way how to apply non-matching and/or Chimera methods within matrix-free loops in deal.II. We are following [112] 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 [112].
Acoustic conservation equations
Acoustic conservation equations are used to describe linear wave propagation. The set of equations consists of the conservation of mass and momentum
Here, \(p\) is the acoustic pressure, \(\mathbf{u}\) the acoustic particle velocity, \(c\) the speed of sound, and \(\rho\) the mean density of the fluid in which waves are propagating. As stated above, the two equations are simply a different way of writing the wave equation: If you take the time derivative of the first equation, and the divergence of the second, i.e., compute
The reason one may want to consider the form we use here (rather than the form used in step-23) is that it has the form of a hyperbolic conservation law in which only first temporal and spatial derivatives appear. Whereas both the more familiar, second order form of the wave equation and the formulation as a first-order system conserve energy, it is often easier to devise numerical schemes that have the right amount of dissipation (necessary for numerical stability) using the well-developed machinery available for first-order systems.
For the discretization of this form, we make use of discontinuous Galerkin (DG) methods. DG methods are especially attractive for the acoustic conservation equations due to their low numerical dispersion errors. More importantly for this tutorial, DG methods natively extend to non-matching Nitsche-type methods [6]. I.e., numerical fluxes are not only used on interior element faces but also as non-matching coupling conditions.
with the penalty parameters \(\tau=\frac{\rho c}{2}\) and \(\gamma=\frac{1}{2\rho c}\). In these formulas, \([a] = a^- - a^+ \) denotes the jump of an arbitrary quantity \(a\) over element faces (face between elements \(K^-\) and \(K^+\)) and \(\jump{a} = a^- \mathbf{n}^- + a^+ \mathbf{n}^+\). For homogeneous materials, the fluxes reduce to standard Lax–Friedrichs fluxes ( \(\gamma^-=\gamma^+\) and \(\tau^-=\tau^+\))
The expression \(\average{a}=\frac{a^- + a^+}{2}\) denotes the averaging operator.
Non-matching discretizations
Non-matching discretizations can be used to connect mesh regions with different element sizes without the need for a transition region. Therefore, they are highly desirable in multiphysics applications. One example is a plate that radiates sound. The plate needs a much finer discretization than the surrounding air because – due to the vastly different speeds of sound in the two media – the wavelengths of sound of the same frequently is very different in the two media, and the mesh size needs to be proportional to the wavelength. We will simulate this example later on.
A different example of the usefulness of non-matching grids is where one wants to move the mesh in parts of the domain, but not others. A typical example is the simulation of windmills: One might want to enclose the rotating wings into a co-rotating mesh (to avoid having to remesh with every time step) but of course the mesh that describes the air above the surrounding landscape and around the tower on which the windmill is located should not rotate. In a case like this, one considers sliding rotating interfaces [76] between the co-rotating part of the mesh and the stationary part of the mesh, but this also requires the ability to handle non-matching discretizations.
Besides this, non-matching methods can be extended to Chimera methods where elements overlap. Chimera methods can help reduce the pressure on mesh generation tools since different regions of a mesh (that may overlap) can be meshed independently.
Different methods exist to treat non-matching interfaces. Within this tutorial, we concentrate on two methods: Point-to-point interpolation and Nitsche-type mortaring.
Point-to-point interpolation
Point-to-point interpolation is a naive approach. Whenever you need to compute integrals over the boundary of the cell at the left, for a coupled problem you then need to evaluate the solution or shape functions on the right at quadrature points of the face on the left, i.e., of the face of element \(K^-\). You can just evaluate these be interpolating the information on the right at these points, but this is in general expensive (read, for example, the documentation of VectorTools::point_value(), which implements this kind of functionality). As can be seen from the picture this approach might be subject to aliasing in some cases.
Nitsche-type mortaring
Mortaring is the process of computing intersections and is not related to the Mortar method which enforces the coupling via Lagrange multipliers. Instead, in mortaring methods one refers to obtained intersections as "mortars". On each mortar a new integration rule is defined. The integral of the face of element \(K^-\) is computed on the intersections. The idea is that if we want to integrate something over a face \(f\subset \partial K^-\), that we break that integral into pieces:
where each of the \(f_i\) corresponds to the intersection of the original face \(f\) on the left with exactly one of the faces on the right; or, if we had multiple faces on the left, then the \(f_i\) would be the intersections of exactly one face on the left and one face on the right.
The point of this approach is first, that splitting the integral this way is exact. Secondly, and maybe more importantly, the terms we are integrating (the dots in the formula above) are now defined on one cell on each side, and consequently are smooth (whereas a finite element solution considered across multiple cells is, in general, not smooth). As a consequence, if we approximate the integrals via numerical integration (quadrature), then the result is exact as long as a sufficient number of integration points is used (at least for affine element shapes; for general curved elements, the integrand will contain rational expressions that are difficult to integrate exactly).
In this tutorial, the intersections are computed using CGAL, the Computational Geometry Algorithms Library. Therefore, deal.II has to be configured with DEAL_II_WITH_CGAL for the Nitsche-type mortaring implementation. See the deal.IIReadme file for information about installation.
FERemoteEvaluation
In practice, for integrals as those mentioned above, we need to evaluate solutions (and shape functions) from cells across the non-matching interface. This is awkward enough if the other side is on the same processor, but outright difficult if the cells on the other side of the interface are owned by a different process in a parallel computation.
On regular meshes (say, doing things as we do in step-40), this is addressed by making sure that we are only computing integrals on locally owned cells and keeping around one layer of ghost cells for which we can query information. Ghost cells are the neighbors of locally owned cells, but in cases like the picture above, where the meshes are not matching, the cells on the other side of the interface are not neighbors in the logical sense – though they happen to be geometrically located adjacently. As a consequence, we need to find a way to efficiently query information on cells that are perhaps located on a different process.
FERemoteEvaluation is a wrapper class which provides a similar interface to, e.g., the FEEvaluation and FEFaceEvaluation classes to access values over non-matching interfaces in matrix-free loops. A detailed description on how to set up the class and how to use it in actual code is given below using hands-on examples. Within this tutorial we only show the usage for non-matching discretizations. Note however, that FERemoteEvaluation can also be used in other settings such as volume coupling. Under the hood, Utilities::MPI::RemotePointEvaluation is used to query the solution or gradients at quadrature points. A detailed description how this is done can be found in step-87. The main difference in the usage of FERemoteEvaluation compared to FEEvaluation is that the interpolated values/gradients of the finite element solution at all the quadrature points are precomputed globally before the loop over the cells/faces of the mesh (i.e., near the place where the communication takes place) instead of performing the interpolation on a cell-by-cell basis. (The principal reason for this design is that MPI has a communication model where you can send messages, but you won't hear back unless the other side is actually listening. As a consequence, you can't generally write code where each process is doing its thing until it needs some information at which point it sends a message to another process to ask for something; because the other process doesn't know that there are such messages, or how many, that have been sent to it, it doesn't respond and so the first process is stuck. Instead, the programming model used with MPI is generally to collect information about everything one will need up front; then each process sends it to all the others; then each process works on these combined requests and sends the required information back to the senders; and at this point everyone has everything they need for their work and can actually do that work.)
The standard code to evaluate fluxes via FEEvaluation on interior faces between two cells reads as follows (where _m corresponds to \(K^{-}\), the current cell in minus normal direction, and _p corresponds to \(K^{+}\), the neighbor cell in plus normal direction):
In DG methods we have to evaluate fluxes over element faces. Exemplarily for an upwind-like flux \(u^*(\mathbf{x}) = u^+(\mathbf{x})\) over element face \(\partial K\) we have to compute
FEFaceEvaluation::gather_evaluate(src, EvaluationFlags::values) and FEFaceEvaluation::get_value(q) extract the value at quadrature point \(\mathbf{x}_q^{\partial K}\) from src. FEFaceEvaluation::submit_value(value, q) multiplies the value with the quadrature weight and the Jacobian determinant at \(\mathbf{x}_q^{\partial K}\). Eventually FEFaceEvaluation::integrate_scatter(EvaluationFlags::values, dst) multiplies the values queued for evaluation by FEFaceEvaluation::submit_value() by the value of the basis functions and writes the result to dst. The corresponding code reads
The code to do the same with FERemoteEvaluation is shown below. For brevity, we assume all boundary faces are somehow connected via non-conforming interfaces for FERemoteEvaluation.
// Initialize FERemoteEvaluation: Note, that FERemoteEvaluation internally manages
// the memory to store precomputed values. Therefore, FERemoteEvaluation
// should be initialized only once to avoid frequent memory
// allocation/deallocation. At this point, remote_communicator is assumed
The object remote_communicator is of type FERemoteEvaluationCommunicator and assumed to be correctly initialized prior to the above code snippet. FERemoteEvaluationCommunicator 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. As mentioned above, the update of the values/gradients happens before the actual matrix-free loop. FERemoteEvaluationCommunicator, as well as FERemoteEvaluation, behaves differently for the given template parameter value_type. If we want to access values at arbitrary points (e.g. in combination with FEPointEvaluation), then we need to choose value_type=Number. If the values are defined at quadrature points of a FEEvaluation object it is possible to get the values at the quadrature points of batches and we need to choose value_type=VectorizedArray<Number>.
Overview of the test case
In this program, we implemented both the point-to-point interpolation and Nitsche-type mortaring mentioned in the introduction.
At first we are considering the test case of a vibrating membrane, see e.g. [163]. Standing waves of length \(\lambda=2/M\) are oscillating with a time period of \(T=2 / (M \sqrt{d} c)\) where \(d\) is the dimension of the space in which our domain is located and \(M\) is the number of modes per meter, i.e. the number of half-waves per meter. The corresponding analytical solution reads as
\begin{align*}
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),
\end{align*}
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 \(\Omega=(0,1)^2\). The domain will be meshed so that the left and right parts of the domain consist of separate meshes that do not match at the interface.
As will become clear from the results, the point-to-point interpolation will result in aliasing, which can be resolved using Nitsche-type mortaring.
In a more realistic second example, we apply this implementation 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 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. A test case like this has been simulated with a different domain and different initial conditions, e.g., in [12].
The commented program
Include files
The program starts with including all the relevant header files.
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/mpi.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/operators.h>
#include <deal.II/non_matching/mapping_info.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
The following header file provides the class FERemoteEvaluation, which allows to access values and/or gradients at remote triangulations similar to FEEvaluation.
In the following, let us first define a function that provides the initial condition for the vibrating membrane test case. It implements both the initial pressure (component 0) and velocity (components 1 to 1 + dim).
There is also a function that computes the duration of one oscillation.
template <int dim>
class InitialConditionVibratingMembrane : publicFunction<dim>
Next up is a function that provides the values of a pressure Gauss pulse. As with the previous function, it implements both the initial pressure (component 0) and velocity (components 1 to 1 + dim).
The following namespace contains free helper functions that are used in the tutorial.
namespace HelperFunctions
{
Helper function to check if a boundary ID is related to a non-matching face. A std::set that contains all non-matching boundary IDs is handed over in addition to the face ID under question.
The following class stores the information if the fluid is homogeneous as well as the material properties at every cell. This class helps access the correct values without accessing a large vector of materials in the homogeneous case.
The class is provided a map from material ids to material properties (given as a pair of values for the speed of sound and the density). If the map has only one entry, the material is homogeneous – using the same values everywhere – and we can remember that fact in the homogeneous member variable and use it to optimize some code paths below. If the material is not homogeneous, we will fill a vector with the correct materials for each batch of cells; this information can then be access via FEEvaluationData::read_cell_data().
As is usual when working with the MatrixFree framework, we will not only access material parameters at a single site, but for whole "batches" of cells. As a consequence, the functions below returned VectorizedArray objects for a batch at a time.
To be able to access the material data in every cell in a thread safe way, the following class MaterialEvaluation is used. Similar to FEEvaluation, functions below will create their own instances of this class; thus, there can be no race conditions even if these functions run multiple times in parallel. For inhomogeneous materials, a reinit_cell() or reinit_face() function is used to set the correct material at the current cell batch. In the homogeneous case the _reinit() functions don't have to reset the materials.
density = material_data.get_homogeneous_density();
}
}
bool is_homogeneous() const
{
return material_data.is_homogeneous();
}
The following two functions update the data for the current cell or face, given a cell batch index. If the material is homogeneous, there is nothing to do. Otherwise, we reinit the FEEvaluation object and store the data for the current object.
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\).
The following class then 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.
At the top we define a flag that decides whether we want to use mortaring. 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. The decision is made using std::is_floating_point_v, which is a variable that is true if the template argument is a floating point type (such as float or double) and false otherwise (in particular if the template argument is a vectorized array type).
The following function then evaluates the acoustic operator. It first updates the precomputed values in corresponding the FERemoteEvaluation objects. The material parameters do not change and thus, we do not have to update precomputed values in c_r_eval and rho_r_eval.
It then either performs a matrix-free loop with Nitsche-type mortaring at non-matching faces (if use_mortaring is true) or with point-to-point interpolation at non-matching faces (in the else branch). The difference is only in the third argument to the loop function, denoting the function object that is executed at boundary faces.
In the private section of the class, we then define the functions that evaluate volume, interior face, and boundary face integrals. The concrete terms these functions evaluate are stated in the documentation at the top of this tutorial program. Each of these functions has its own object of type MaterialEvaluation that provides access to the material at each cell or face.
Get the materials on the corresponding cell. Since we introduced MaterialEvaluation we can write the code independent of whether the material is homogeneous or inhomogeneous.
material.reinit_cell(cell);
constauto c = material.get_speed_of_sound();
constauto rho = material.get_density();
for (constunsignedint q : pressure.quadrature_point_indices())
{
pressure.submit_value(rho * c * c * velocity.get_divergence(q),
This next 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).
In this and following functions, we also introduce the factors \(\tau\) and \(\gamma\) that are computed from \(\rho\) and \(c\) along interfaces and that appear in the bilinear forms. Their definitions are provided in the introduction.
pressure_m.submit_value(rho * c * c * (mass_flux - um) * n, q);
ifconstexpr (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.
If face is non-matching we have to query values via the FERemoteEvaluaton objects. This is done by passing the corresponding FERemoteEvaluaton objects to the function that evaluates the kernel. As mentioned above, each side of the non-matching interface is traversed separately and we do not have to consider the neighbor in the kernel. Note, that the values in the FERemoteEvaluaton objects are already updated at this point.
For point-to-point interpolation we simply use the corresponding FERemoteEvaluaton objects in combination with the standard FEFaceEvaluation objects.
velocity_r.reinit(face);
pressure_r.reinit(face);
material.reinit_face(face);
If we are considering a homogeneous material, do not use the inhomogeneous fluxes. While it would be possible to use the inhomogeneous fluxes they are more expensive to compute.
if (material.is_homogeneous())
{
evaluate_face_kernel<false>(pressure_m,
velocity_m,
pressure_r,
velocity_r,
material.get_speed_of_sound(),
material.get_density());
}
else
{
c_r.reinit(face);
rho_r.reinit(face);
evaluate_face_kernel_inhomogeneous(
pressure_m,
velocity_m,
pressure_r,
velocity_r,
material.get_speed_of_sound(),
material.get_density(),
c_r,
rho_r);
}
}
else
{
If face is a standard boundary face, evaluate the integral as usual in the matrix free context. To be able to use the same kernel as for inner faces we pass the boundary condition objects to the function that evaluates the kernel. As detailed above weight_neighbor=false.
For Nitsche-type mortaring we are evaluating the integrals over intersections. This is why, quadrature points are arbitrarily distributed on every face. Thus, we can not make use of face batches and FEFaceEvaluation but have to consider each face individually and make use of FEFacePointEvaluation to evaluate the integrals in the arbitrarily distributed quadrature points. Since the setup of FEFacePointEvaluation is more expensive than that of FEEvaluation we do the setup only once. For this we are using the helper function get_thread_safe_fe_face_point_evaluation_object().
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 (unsignedint v = 0;
v < matrix_free.n_active_entries_per_face_batch(face);
As above, if we are considering a homogeneous material, do not use the inhomogeneous fluxes. Since we are operating on face v we call material.get_density()[v].
FERemoteEvaluation objects are stored as shared pointers. This way, they can also be used for other operators without precomputing the values multiple times.
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.
For the time stepping methods below, we need the inverse mass operator. We apply it via a loop over all (batches of) cells as always, where the contribution of each cell is computed in a matrix-free way:
The run() function of this class contains the time loop. It starts by initializing some member variables (such as short-cuts to objects that describe the finite element, its properties, and the mapping; an by initializing vectors). It then computes the time step size via minimum element edge length. Assuming non-distorted elements, we can compute the edge length as the distance between two vertices. From this, we can then obtain the time step size via the CFL condition.
The main work of this class is done by a private member function that performs one Runge-Kutta 2 time step. Recall that this method requires two sub-steps ("stages") computing intermediate values k1 and k2, after which the intermediate values are summed with weights to obtain the new solution at the end of the time step. The RK2 method allows for the elimination of the intermediate vector k2 by utilizing the dst vector as temporary storage.
Evaluating a single Runge-Kutta stage is a straightforward step that really only requires applying the operator once, and then applying the inverse of the mass matrix.
Let us now make our way to the higher-level functions of this program.
The first of these functions creates a two dimensional square 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 elements that are three times smaller compared to the ones for the right sub-domain.
At non-matching interfaces, we need to provide different boundary IDs for the cells that make up the two parts (because, while they may be physically adjacent, they are not logically neighbors given that the faces of cells on both sides do not match, and the Triangulation class will therefore treat the interface between the two parts as a "boundary"). 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. As a consequence, we declare at the top symbolic names for these boundary indicators, and ensure that we return a set with these values to the caller for later use.
The actual mesh is then constructed by first constructing the left and right parts separately (setting material ids to zero and one, respectively), and using the appropriate boundary ids for all parts of the mesh. We then use GridGenerator::merge_triangulations() to combine them into one (non-matching) mesh. We have to pay attention that should the two sub-triangulations have vertices at the same locations, that they are not merged (connecting the two triangulations logically) since we want the interface to be an actual boundary. We achieve this by providing a tolerance of zero for the merge, see the documentation of the function GridGenerator::merge_triangulations().
Setup and running of the point-to-point interpolation scheme
We are now at the two functions that run the overall schemes (the point-to-point and the mortaring schemes). The first of these functions fills 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.
Communication objects know about the communication pattern. That is, they know about the cells and quadrature points that have to be evaluated at remote faces. This information is given via RemotePointEvaluation. Additionally, the communication objects have to be able to match the quadrature points of the remote points (that provide exterior information) to the quadrature points defined at the interior cell. In case of point-to-point interpolation a vector of pairs with face batch Ids and the number of faces in the batch is needed. FERemoteCommunicationObjectEntityBatches is a container to store this information.
The information is filled outside of the actual class since in some cases the information is available from some heuristic and it is possible to skip some expensive operations. This is for example the case for sliding rotating interfaces with equally spaced elements on both sides of the non-matching interface [76].
For the standard case of point to point-to-point interpolation without any heuristic we make use of the utility function Utilities::compute_remote_communicator_faces_point_to_point_interpolation(). Please refer to the documentation of this function to see how to manually set up the remote communicator from outside.
Among the inputs for the remote communicator we need a list of boundary ids for the non-matching faces, along with a function object for each boundary id that returns a vector of true/false flags in which exactly the vertices of cells of the triangulation are marked that have a face at the boundary id in question.
We are using point-to-point interpolation and can therefore easily access the corresponding data at face batches. This is why we use a VectorizedArray as remote_value_type:
We then set up FERemoteEvaluation objects that access the pressure and velocity at remote faces, along with an object to describe cell-wise material data.
If the domain is not homogeneous, i.e., if material parameters change from cell to cell, we initialize and fill DoF vectors that contain the material properties. Materials do not change during the simulation, therefore there is no need to ever compute the values after the first gather_evaluate() (at the end of the following block) again.
Next, we set up the inverse mass operator and the acoustic operator. Using remote_value_type=VectorizedArray<Number> makes the operator use point-to-point interpolation. These two objects are then used to create a RungeKutta2 object to perform the time integration.
We also compute the maximum speed of sound, needed for the computation of the time-step size, and then run the time integrator. For the examples considered here, 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\).
Setup and running of the Nitsche-type mortaring scheme
The alternative to the previous function is to use the mortaring scheme – implemented in the following function. This function can only be run when deal.II is configured using CGAL (but the previous function can still be used without CGAL), and so errors out if CGAL is not available.
The only parts in this function that are functionally different from the previous one follow here.
First, quadrature points are arbitrarily distributed on each non-matching face. Therefore, we have to make use of FEFacePointEvaluation. FEFacePointEvaluation needs NonMatching::MappingInfo to work at the correct quadrature points that are in sync with used FERemoteEvaluation object. Using compute_remote_communicator_faces_nitsche_type_mortaring() to reinit NonMatching::MappingInfo ensures this. In the case of mortaring, we have to use the weights provided by the quadrature rules that are used to set up NonMatching::MappingInfo. Therefore we set the flag use_global_weights.
Second, since quadrature points are arbitrarily distributed we have to consider each face in a batch separately and can not make use of VecorizedArray.
using remote_value_type = Number;
The rest of the code is then identical to what we had in the previous function (though it functions differently because of the difference in remote_value_type).
Vibrating membrane: Point-to-point interpolation vs. Nitsche-type mortaring
We compare the results of the simulations after the last time step, i.e. at \(t=8T\). The \(y\)-component of the velocity field using Nitsche-type mortaring is depicted on the left. The same field using point-to-point interpolation is depicted on the right.
Besides this, the results for the pressure and the velocity in \(y\) direction are plotted along the horizontal line that spans from (0,0.69) to (1,0.69).
While the results of the pressure are similar, \(u_y\) clearly differs. 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 becomes 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 [112].
Wave propagation through in-homogeneous fluid
The example that follows is just one example in which non-matching discretizations can be efficiently used to reduce the number of DoFs. The example is nice, since results for a similar test case are shown in multiple publications. As before, we slightly adapted the test case to be able to use the same mesh for all simulations. The pressure field at \(t=0.3\) is depicted below.
As expected, we can easily see that the wave length in the right domain is roughly three times times the wave length in the left domain. Hence, the wave can be resolved with a coarser discretization.
Using the same element size in the whole domain, we can compute a reference result. The displayed reference result is obtained by choosing the same subdivision level for both sub-domains, i.e. subdiv_right = 11. In this particular example the reference result uses \(92,928\) DoFs, while the non-matching result uses \(52,608\) DoFs. The pressure result is plotted along the horizontal line that spans from (0,0.5) to (1,0.5).
The results computed using the non-matching discretization are clearly in good agreement with the reference result.
Possibilities for extensions
All the implementations are done with overlapping triangulations in mind. In particular the intersections in the mortaring case are constructed such that they are computed correctly for overlapping triangulations. For this the intersection requests are of dimension \(dim-1\). The cells which are intersected with the intersection requests are of dimension \(dim\). For the simple case of non-conforming interfaces it would be sufficient to compute the intersections between \(dim-1\) and \(dim-1\) entities. Furthermore, the lambda could be adapted, such that cells are only marked if they are connected to a certain boundary ID (in this case, e.g., 99) instead of marking every cell that is not connected to the opposite boundary ID (in this case, e.g., 98). Marking fewer cells can reduce the setup costs significantly.
Note that the use of inhomogeneous materials in this procedure is questionable, since it is not clear which material is present in the overlapping part of the mesh.
1.11.0