From: Johannes Heinz <43043310+jh66637@users.noreply.github.com> Date: Fri, 2 Feb 2024 15:32:21 +0000 (+0100) Subject: Add tutorial on Nitsche-type mortaring X-Git-Tag: v9.6.0-rc1~398^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=73194f7530d8c42c16f46ded94531a5bb1985c52;p=dealii.git Add tutorial on Nitsche-type mortaring --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index b82b2a5d7b..85c081c077 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1856,6 +1856,71 @@ } +%------------------------------------------------------------------------------- +% 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 index 0000000000..f8f70b4659 --- /dev/null +++ b/doc/news/changes/major/20231127HeinzBergbauerFederMunch @@ -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. +
+(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 index 0000000000..9795f5d7a3 --- /dev/null +++ b/examples/step-89/CMakeLists.txt @@ -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 index 0000000000..809da9b9ad --- /dev/null +++ b/examples/step-89/doc/builds-on @@ -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 index 0000000000..666eef93e2 --- /dev/null +++ b/examples/step-89/doc/intro.dox @@ -0,0 +1,221 @@ +
+ + +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]. + + +@dealiiTutorialDOI{10.5281/zenodo.10033975,https://zenodo.org/badge/DOI/10.5281/zenodo.10033975.svg} + +

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 @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. + +

%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 +@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. + +

%Non-matching discretizations

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

%Point-to-point interpolation

+ +@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. + +

%Nitsche-type mortaring

+ +@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. + +

%FERemoteEvaluation

+ +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 before +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({}, 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 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({}, {}, 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 +before 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 batches and @c value_type=VectorizedArray. + +

Overview

+ +In the following, point-to-point interpolation and Nitsche-type mortaring is implemented. + +At first we are considering the test case of a vibrating membrane, see e.g. @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 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-89/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-89/doc/results.dox b/examples/step-89/doc/results.dox new file mode 100644 index 0000000000..141fba9a97 --- /dev/null +++ b/examples/step-89/doc/results.dox @@ -0,0 +1,80 @@ +

Results

+ +

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

+ +We compare the results of the simulations after the last time step, i.e. at $t=8T$. +The $y$-component of the velocity field using Nitsche-type mortaring is depicted on the left. +The same field using point-to-point interpolation is depicted on the right. + + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_velocity_Y.png "" width=60% + + @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_point_to_point_velocity_Y.png "" width=60% +
+ +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). + + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_pressure.svg "" width=100% + + @image html https://www.dealii.org/images/steps/developer/step_89_membrane_test_case_mortaring_vs_point_to_point_velocity_Y.svg "" width=100% +
+ +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. + +

Wave propagation through in-homogeneous fluid

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

Possibilities for extensions

+ +All the implementations are done with overlapping triangulations in mind. In particular the +intersections in the mortaring case are constructed such that they are computed correctly +for overlapping triangulations. For this the intersection requests are of dimension $dim-1$. +The cells which are intersected with the intersection requests are of dimension $dim$. For the +simple case of non-conforming interfaces it would be sufficient to compute the intersections +between $dim-1$ and $dim-1$ entities. Furthermore, the lambda could be adapted, such that cells are +only marked if they are connected to a certain boundary ID (in this case, e.g. 99) instead of +marking every cell that is not connected to the opposite boundary ID (in this case, e.g. 98). +Marking less cells can reduce the setup costs significantly. + +Note that for in-homogeneous material in this procedure is questionable, since it is not clear which +material is present in the overlapping part of the mesh. diff --git a/examples/step-89/doc/tooltip b/examples/step-89/doc/tooltip new file mode 100644 index 0000000000..10f51f3381 --- /dev/null +++ b/examples/step-89/doc/tooltip @@ -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 index 0000000000..b39ffb885c --- /dev/null +++ b/examples/step-89/step-89.cc @@ -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 +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include +// The following header file provides the class FERemoteEvaluation, which allows +// to access values and/or gradients at remote triangulations similar to +// FEEvaluation. +#include + +// 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 + class InitialConditionVibratingMembrane : public Function + { + 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 &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 + InitialConditionVibratingMembrane::InitialConditionVibratingMembrane( + const double modes) + : Function(dim + 1, 0.0) + , M(modes) + { + static_assert(dim == 2, "Only implemented for dim==2"); + } + + template + double + InitialConditionVibratingMembrane::value(const Point &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 + double InitialConditionVibratingMembrane::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 + class GaussPulse : public Function + { + 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 &p, const unsigned int comp) const final; + + private: + const double shift_x; + const double shift_y; + }; + + template + GaussPulse::GaussPulse(const double shift_x, const double shift_y) + : Function(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 + double GaussPulse::value(const Point &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 &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 + void set_initial_condition(MatrixFree matrix_free, + const Function &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 + void write_vtu(const VectorType &solution, + const DoFHandler &dof_handler, + const Mapping &mapping, + const unsigned int degree, + const std::string &name_prefix) + { + DataOut data_out; + DataOutBase::VtkFlags flags; + flags.write_higher_order_cells = true; + data_out.set_flags(flags); + + std::vector + interpretation( + dim + 1, DataComponentInterpretation::component_is_part_of_vector); + std::vector 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::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 + class CellwiseMaterialData + { + public: + template + CellwiseMaterialData( + const MatrixFree> &matrix_free, + const std::map> + &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> &get_speed_of_sound() const + { + Assert(!homogeneous, ExcMessage("Use get_homogeneous_speed_of_sound()")); + return speed_of_sound; + } + + const AlignedVector> &get_density() const + { + Assert(!homogeneous, ExcMessage("Use get_homogeneous_density()")); + return density; + } + + VectorizedArray get_homogeneous_speed_of_sound() const + { + Assert(homogeneous, ExcMessage("Use get_speed_of_sound()")); + return speed_of_sound_homogeneous; + } + + VectorizedArray get_homogeneous_density() const + { + Assert(homogeneous, ExcMessage("Use get_density()")); + return density_homogeneous; + } + + private: + const bool homogeneous; + + // Materials in the in-homogeneous case. + AlignedVector> speed_of_sound; + AlignedVector> density; + + // Materials in the homogeneous case. + VectorizedArray speed_of_sound_homogeneous; + VectorizedArray 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 + class MaterialEvaluation + { + public: + MaterialEvaluation( + const MatrixFree> &matrix_free, + const CellwiseMaterialData &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 get_speed_of_sound() const + { + return speed_of_sound; + } + + // Return the density at the current cell batch. + VectorizedArray get_density() const + { + return density; + } + + private: + // Members needed for the in-homogeneous case. + FEEvaluation phi; + FEFaceEvaluation phi_face; + + // Material defined at every cell. + const CellwiseMaterialData &material_data; + + // Materials at current cell. + VectorizedArray speed_of_sound; + VectorizedArray 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 + class BCEvaluationP + { + public: + BCEvaluationP(const FEFaceEvaluation &pressure_m) + : pressure_m(pressure_m) + {} + + typename FEFaceEvaluation::value_type + get_value(const unsigned int q) const + { + return -pressure_m.get_value(q); + } + + private: + const FEFaceEvaluation &pressure_m; + }; + + // We don't have to apply boundary conditions for the velocity, i.e. + // $\mathbf{u}_h^+=\mathbf{u}_h^-$. + template + class BCEvaluationU + { + public: + BCEvaluationU(const FEFaceEvaluation &velocity_m) + : velocity_m(velocity_m) + {} + + typename FEFaceEvaluation::value_type + get_value(const unsigned int q) const + { + return velocity_m.get_value(q); + } + + private: + const FEFaceEvaluation &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 + 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; + + public: + // In case of Nitsche-type mortaring, `NonMatching::MappingInfo` has to + // be provided in the constructor. + AcousticOperator( + const MatrixFree &matrix_free, + std::shared_ptr> material_data, + const std::set &remote_face_ids, + std::shared_ptr> + pressure_r_eval, + std::shared_ptr> + velocity_r_eval, + std::shared_ptr> c_r_eval, + std::shared_ptr> rho_r_eval, + std::shared_ptr> 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 + 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, + &AcousticOperator::local_apply_face, + &AcousticOperator::local_apply_boundary_face_mortaring, + this, + dst, + src, + true, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + } + else + { + // Perform matrix free loop with point-to-point interpolation at + // non-matching faces. + matrix_free.loop( + &AcousticOperator::local_apply_cell, + &AcousticOperator::local_apply_face, + &AcousticOperator::local_apply_boundary_face_point_to_point< + VectorType>, + this, + dst, + src, + true, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + } + } + + private: + // This function evaluates the volume integrals. + template + void local_apply_cell( + const MatrixFree &matrix_free, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + FEEvaluation pressure(matrix_free, 0, 0, 0); + FEEvaluation 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 + 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 + 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 + void local_apply_face( + const MatrixFree &matrix_free, + VectorType &dst, + const VectorType &src, + const std::pair &face_range) const + { + FEFaceEvaluation pressure_m( + matrix_free, true, 0, 0, 0); + FEFaceEvaluation pressure_p( + matrix_free, false, 0, 0, 0); + FEFaceEvaluation velocity_m( + matrix_free, true, 0, 0, 1); + FEFaceEvaluation 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(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 + void local_apply_boundary_face_point_to_point( + const MatrixFree &matrix_free, + VectorType &dst, + const VectorType &src, + const std::pair &face_range) const + { + // Standard face evaluators. + FEFaceEvaluation pressure_m( + matrix_free, true, 0, 0, 0); + FEFaceEvaluation 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(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(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 + void local_apply_boundary_face_mortaring( + const MatrixFree &matrix_free, + VectorType &dst, + const VectorType &src, + const std::pair &face_range) const + { + // Standard face evaluators for BCs. + FEFaceEvaluation pressure_m( + matrix_free, true, 0, 0, 0); + FEFaceEvaluation 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 &velocity_m_mortar = + get_thread_safe_fe_face_point_evaluation_object( + 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::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( + 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(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 &matrix_free; + + // CellwiseMaterialData is stored as shared pointer with the same + // argumentation. + const std::shared_ptr> material_data; + + const std::set 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> + pressure_r_eval; + const std::shared_ptr> + velocity_r_eval; + + const std::shared_ptr> + c_r_eval; + const std::shared_ptr> + rho_r_eval; + + const std::shared_ptr> + 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>> + thread_local_pressure_m_mortar; + + mutable Threads::ThreadLocalStorage< + std::unique_ptr>> + 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 + FEFacePointEvaluation & + get_thread_safe_fe_face_point_evaluation_object( + Threads::ThreadLocalStorage< + std::unique_ptr>> + &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>( + *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 + class InverseMassOperator + { + public: + InverseMassOperator(const MatrixFree &matrix_free) + : matrix_free(matrix_free) + {} + + // Function to apply the inverse mass operator. + template + void apply(VectorType &dst, const VectorType &src) const + { + dst.zero_out_ghost_values(); + matrix_free.cell_loop(&InverseMassOperator::local_apply_cell, + this, + dst, + src); + } + + private: + // Apply the inverse mass operator onto every cell batch. + template + void local_apply_cell( + const MatrixFree &mf, + VectorType &dst, + const VectorType &src, + const std::pair &cell_range) const + { + FEEvaluation phi(mf); + MatrixFreeOperators::CellwiseInverseMassMatrix + 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 &matrix_free; + }; + + //@sect3{Runge-Kutta time-stepping} + // + // This class implements a Runge-Kutta scheme of order 2. + template + class RungeKutta2 + { + using VectorType = LinearAlgebra::distributed::Vector; + + public: + RungeKutta2( + const std::shared_ptr> + inverse_mass_operator, + const std::shared_ptr> + acoustic_operator) + : inverse_mass_operator(inverse_mass_operator) + , acoustic_operator(acoustic_operator) + {} + + // Set up and run time loop. + void run(const MatrixFree &matrix_free, + const double cr, + const double end_time, + const double speed_of_sound, + const Function &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::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> + inverse_mass_operator; + const std::shared_ptr> + 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 + void build_non_matching_triangulation( + Triangulation &tria, + std::set &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 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 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 + void run_with_point_to_point_interpolation( + const MatrixFree &matrix_free, + const std::set &non_matching_faces, + const std::map> &materials, + const double end_time, + const Function &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()>>> + 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 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; + + // Set up FERemoteEvaluation object that accesses the pressure + // at remote faces. + const auto pressure_r = + std::make_shared>( + 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>( + remote_communicator, dof_handler, /*first_selected_component*/ 1); + + // Set up cell-wise material data. + const auto material_data = + std::make_shared>(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>( + remote_communicator, + matrix_free.get_dof_handler().get_triangulation(), + /*first_selected_component*/ 0); + const auto rho_r = + std::make_shared>( + 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 c( + matrix_free.get_dof_handler().get_triangulation().n_active_cells()); + Vector 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>(matrix_free); + + // Set up the acoustic operator. Using + // `remote_value_type=VectorizedArray` makes the operator use + // point-to-point interpolation. + const auto acoustic_operator = + std::make_shared>( + 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 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 + void run_with_nitsche_type_mortaring( + const MatrixFree &matrix_free, + const std::set &non_matching_faces, + const std::map> &materials, + const double end_time, + const Function &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()>>> + 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 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::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>( + 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>( + remote_communicator, dof_handler, /*first_selected_component*/ 0); + + const auto velocity_r = + std::make_shared>( + remote_communicator, dof_handler, /*first_selected_component*/ 1); + + const auto material_data = + std::make_shared>(matrix_free, materials); + + const auto c_r = + std::make_shared>( + remote_communicator, + matrix_free.get_dof_handler().get_triangulation(), + /*first_selected_component*/ 0); + const auto rho_r = + std::make_shared>( + remote_communicator, + matrix_free.get_dof_handler().get_triangulation(), + /*first_selected_component*/ 0); + + if (!material_data->is_homogeneous()) + { + Vector c( + matrix_free.get_dof_handler().get_triangulation().n_active_cells()); + Vector 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>(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>( + 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 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 tria(MPI_COMM_WORLD); +#else + parallel::shared::Triangulation tria(MPI_COMM_WORLD); +#endif + + pcout << "Create non-matching grid..." << std::endl; + + std::set 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 dof_handler(tria); + dof_handler.distribute_dofs(FESystem(FE_DGQ(degree), dim + 1)); + pcout << " - Number of DoFs: " << dof_handler.n_dofs() << std::endl; + + AffineConstraints constraints; + constraints.close(); + + pcout << "Set up MatrixFree..." << std::endl; + typename MatrixFree::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 matrix_free; + matrix_free.reinit( + MappingQ1(), dof_handler, constraints, QGauss(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> homogeneous_material; + homogeneous_material[numbers::invalid_material_id] = std::make_pair(1.0, 1.0); + const auto initial_solution_membrane = + Step89::InitialConditionVibratingMembrane(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> + 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(0.3, 0.5), + "inhomogeneous"); + + + return 0; +}