From: Peter Munch Date: Thu, 10 Aug 2023 14:29:42 +0000 (+0200) Subject: Add step-87 X-Git-Tag: relicensing~432^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16072%2Fhead;p=dealii.git Add step-87 --- diff --git a/doc/doxygen/references.bib b/doc/doxygen/references.bib index fdaa7f4a4e..f443f79b47 100644 --- a/doc/doxygen/references.bib +++ b/doc/doxygen/references.bib @@ -1697,6 +1697,152 @@ } +%------------------------------------------------------------------------------- +% Step 87 +%------------------------------------------------------------------------------- + +@article{henri2022geometrical, + title={Geometrical level set reinitialization using closest point method and kink detection for thin filaments, topology changes and two-phase flows}, + author={Henri, F{\'e}lix and Coquerelle, Mathieu and Lubin, Pierre}, + journal={Journal of Computational Physics}, + volume={448}, + pages={110704}, + year={2022}, + publisher={Elsevier} +} + +@article{coquerelle2016fourth, + title={A fourth-order accurate curvature computation in a level set framework for two-phase flows subjected to surface tension forces}, + author={Coquerelle, Mathieu and Glockner, St{\'e}phane}, + journal={Journal of Computational Physics}, + volume={305}, + pages={838--876}, + year={2016}, + publisher={Elsevier} +} + +@article{brackbill1992continuum, + title={A continuum method for modeling surface tension}, + author={Brackbill, Jeremiah U and Kothe, Douglas B and Zemach, Charles}, + journal={Journal of computational physics}, + volume={100}, + number={2}, + pages={335--354}, + year={1992}, + publisher={Elsevier} +} + +@article{olsson2005conservative, + title={A conservative level set method for two phase flow}, + author={Olsson, Elin and Kreiss, Gunilla}, + journal={Journal of computational physics}, + volume={210}, + number={1}, + pages={225--246}, + year={2005}, + publisher={Elsevier} +} + +@article{kronbichler2018fast, + title={A fast massively parallel two-phase flow solver for microfluidic chip simulation}, + author={Kronbichler, Martin and Diagne, Ababacar and Holmgren, Hanna}, + journal={The International Journal of High Performance Computing Applications}, + volume={32}, + number={2}, + pages={266--287}, + year={2018}, + publisher={SAGE Publications Sage UK: London, England} +} + +@article{peskin1977numerical, + title={Numerical analysis of blood flow in the heart}, + author={Peskin, Charles S}, + journal={Journal of computational physics}, + volume={25}, + number={3}, + pages={220--252}, + year={1977}, + publisher={Elsevier} +} + +@article{unverdi1992front, + title={A front-tracking method for viscous, incompressible, multi-fluid flows}, + author={Unverdi, Salih Ozen and Tryggvason, Gr{\'e}tar}, + journal={Journal of computational physics}, + volume={100}, + number={1}, + pages={25--37}, + year={1992}, + publisher={Elsevier} +} + +@article{bungartz2016precice, + title={{preCICE} -- a fully parallel library for multi-physics surface coupling}, + author={Bungartz, Hans-Joachim and Lindner, Florian and Gatzhammer, Bernhard and Mehl, Miriam and Scheufele, Klaudius and Shukaev, Alexander and Uekermann, Benjamin}, + journal={Computers \& Fluids}, + volume={141}, + pages={250--258}, + year={2016}, + publisher={Elsevier} +} + +@article{chourdakis2021precice, + title={{preCICE} v2: A sustainable and user-friendly coupling library}, + author={Chourdakis, Gerasimos and Davis, Kyle and Rodenberg, Benjamin and Schulte, Miriam and Simonis, Fr{\'e}d{\'e}ric and Uekermann, Benjamin and Abrams, Georg and Bungartz, Hans-Joachim and Yau, Lucia Cheung and Desai, Ishaan and others}, + journal={arXiv preprint arXiv:2109.14470}, + year={2021} +} + +@article{lebrun2020arborx, + title={ArborX: A performance portable geometric search library}, + author={Lebrun-Grandi{\'e}, Damien and Prokopenko, Andrey and Turcksin, Bruno and Slattery, Stuart R}, + journal={ACM Transactions on Mathematical Software (TOMS)}, + volume={47}, + number={1}, + pages={1--15}, + year={2020}, + publisher={ACM New York, NY, USA} +} + +@article{heinz2022high, + title={High-Order Non-Conforming Discontinuous {G}alerkin Methods for the Acoustic Conservation Equations}, + author={Heinz, Johannes and Munch, Peter and Kaltenbacher, Manfred}, + journal={International Journal for Numerical Methods in Engineering}, + year={2022}, + volume={124}, + number={9}, + pages={2034-2049}, + publisher={Wiley Online Library} +} + + + +%%% sparse communication pattern + +@article{hoefler2010scalable, + title={Scalable communication protocols for dynamic sparse data exchange}, + author={Hoefler, Torsten and Siebert, Christian and Lumsdaine, Andrew}, + journal={ACM Sigplan Notices}, + volume={45}, + number={5}, + pages={159--168}, + year={2010}, + publisher={ACM New York, NY, USA} +} + +@article{burman2015cutfem, + title={CutFEM: discretizing geometry and partial differential equations}, + author={Burman, Erik and Claus, Susanne and Hansbo, Peter and Larson, Mats G and Massing, Andr{\'e}}, + journal={International Journal for Numerical Methods in Engineering}, + volume={104}, + number={7}, + pages={472--501}, + year={2015}, + publisher={Wiley Online Library} +} + + + %------------------------------------------------------------------------------- % References used elsewhere %------------------------------------------------------------------------------- diff --git a/doc/doxygen/tutorial/tutorial.h.in b/doc/doxygen/tutorial/tutorial.h.in index 7818674238..57b6e60f2a 100644 --- a/doc/doxygen/tutorial/tutorial.h.in +++ b/doc/doxygen/tutorial/tutorial.h.in @@ -684,6 +684,13 @@ *
Keywords: FEInterfaceValues, NonMatching::FEImmersedSurfaceValues * * + * + * step-87 + * Evaluation of finite element solutions at arbitrary points within a distributed + * mesh with application to two-phase flow. + *
Keywords: Utilities::MPI::RemotePointEvaluation, VectorTools::point_values() + * + * * * * @@ -1027,7 +1034,8 @@ * * step-60, * step-70, - * step-85 + * step-85, + * step-87 * * * diff --git a/doc/news/changes/major/20230905SchreterMunch b/doc/news/changes/major/20230905SchreterMunch new file mode 100644 index 0000000000..3216cf34ec --- /dev/null +++ b/doc/news/changes/major/20230905SchreterMunch @@ -0,0 +1,5 @@ +New: The new tutorial step-87 presents the advanced point evaluation +functionalities of deal.II, specifically useful for evaluating +finite element solutions at arbitrary points on distributed meshes. +
+(Magdalena Schreter-Fleischhacker, Peter Munch, 2023/09/05) diff --git a/examples/step-87/CMakeLists.txt b/examples/step-87/CMakeLists.txt new file mode 100644 index 0000000000..735c52b1d0 --- /dev/null +++ b/examples/step-87/CMakeLists.txt @@ -0,0 +1,54 @@ +## +# CMake script +## + +# Set the name of the project and target: +set(TARGET "step-87") + +# 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-87/doc/builds-on b/examples/step-87/doc/builds-on new file mode 100644 index 0000000000..000aea605e --- /dev/null +++ b/examples/step-87/doc/builds-on @@ -0,0 +1 @@ +step-40 diff --git a/examples/step-87/doc/intro.dox b/examples/step-87/doc/intro.dox new file mode 100644 index 0000000000..4d6b3250a5 --- /dev/null +++ b/examples/step-87/doc/intro.dox @@ -0,0 +1,312 @@ +
+ + +This program was contributed by Magdalena Schreter-Fleischhacker +and Peter Munch. Many ideas presented here are the result of common code +development with Maximilian Bergbauer, Marco Feder, +Niklas Fehn, Johannes Heinz, Luca Heltai, Martin Kronbichler, Nils Much, +and Judith Pauen. + +This tutorial is loosely based on chapter 3.4 of the submitted +PhD thesis "Matrix-free finite element computations at extreme scale and for +challenging applications" by Peter Munch. Magdalena Schreter-Fleischhacker +is funded by the Austrian Science Fund (FWF) Schrödinger Fellowship (J4577). + + +

Introduction

+ +This tutorial presents the advanced point-evaluation functionalities of +deal.II, specifically useful for evaluating finite element solutions at +arbitrary points. The underlying finite element mesh can be distributed +among processes, which makes the operations more involved due to communication. +In the examples discussed in this tutorial, we focus on point evaluation for +MPI-parallel computations, like parallel::distributed::Triangulation. +Nevertheless, the application to non-distributed meshes is also possible. + +

%Point evaluation

+ +In the context of the finite element method (FEM), it is a common task to +query the solution $u$ at an arbitrary point $\boldsymbol{x}_q$ in the +domain of interest $\Omega$ +@f[ +u(\boldsymbol{x}_q) = \sum_{i} N_i(\boldsymbol{x}_q) u_i \quad\text{with}\quad +i\in[0,n_{\text{dofs}}), +@f] +by evaluating the shape functions $N_i$ at this point together with the +corresponding solution coefficients $u_i$. +After identifying the cell $K$ where the point $\boldsymbol{x}_q$ resides, the +transformation between $\boldsymbol{x}_q$ and the corresponding coordinates in +the reference cell $\hat{\boldsymbol{x}}_q$ is obtained by the mapping +$\boldsymbol{x}_q=\boldsymbol{F}_K(\hat{\boldsymbol{x}}_q)$. In this setting, +the evaluation of the solution at an arbitrary point boils down to a cell-local +evaluation +@f[ +u(\boldsymbol{x}_q) = \sum_{i} \hat{N}^K_i(\hat{\boldsymbol{x}}_q) u_i^K +\quad\text{with}\quad i\in[0,n_{\text{dofs_per_cell}}), +@f] +with $\hat{N}^K_i$ being the shape functions defined on the reference cell and +$u_i^{K}$ the solution coefficients +restricted to the cell $K$. + +Alternatively to point evaluation, evaluating weak-form (integration) operations +of the type +@f[ +u_i += +\left(N_i(\boldsymbol{x}), u(\boldsymbol{x})\right)_\Omega += +\int_{\Omega} N_i(\boldsymbol{x}) u(\boldsymbol{x}) dx += +\sum_q N_i\left(\boldsymbol{x}_q\right) u\left(\boldsymbol{x}_q\right) |J\left( +\boldsymbol{x}_q\right)| w\left(\boldsymbol{x}_q\right) \quad\text{with}\quad +i\in[0,n_{\text{dofs}}) +@f] +is possible, with $\boldsymbol{x}_q$ being quadrature points at arbitrary +positions. +After the values at the quadrature points have been multiplied by the +integration weights, this operation can be interpreted as the transpose of the +evaluation. Not surprisingly, such an operation can be also implemented as a +cell loop. + +

Setup and communication

+ +To determine the cell $K$ and the reference position $\hat{\boldsymbol x}_q$ +within the cell for a given point $\boldsymbol{x}_q$ on distributed meshes, +deal.II performs a two-level-search approach. First, all processes whose portion +of the global mesh might contain the point are determined ("coarse search"). +For this purpose, e.g., a distributed tree based on bounding boxes around +locally owned domains using "ArborX" @cite lebrun2020arborx is applied. After +the potentially owning processes have been determined and the points have been +sent to them as a request, one can start to find the cells that surround the +points among locally owned cells ("fine search"). In order to accelerate this +search, an R-tree from "boost::geometry" built around the vertices of the mesh +is used. + +Once the cell $K$ that surrounds point $\boldsymbol{x}_q$ has been found, +the reference position $\hat{\boldsymbol{x}}_q$ is obtained by performing the +minimization: +@f[ +\min_{\hat{\boldsymbol{x}}_q}(| \boldsymbol{F}_K(\hat{\boldsymbol{x}}_q) +- {\boldsymbol{x}_q}|) +\quad +\text{with} +\quad +\hat{\boldsymbol{x}}_q\in[0,1]^{dim}. +@f] +With the determined pieces of information, the desired evaluation can be +performed by the process that owns the cell. The result can now be communicated +to the requesting process. + +In summary, the coarse search determines, for each point, a list of processes +that might own it. The subsequent fine search by each process determines whether +the processes actually own these points by the sequence of request +("Does the process own the point?") and answer ("Yes."/"No."). +Processes might post any number of point requests and communicate with +any process. We propose to collect the point requests to a process to +use the dynamic, sparse, scalable consensus-based communication +algorithms @cite hoefler2010scalable, and to consider the obtained information +to set up point-to-point communication patterns. + +

Implementation: Utilities::MPI::RemotePointEvaluation

+ +The algorithm described above is implemented in +Utilities::MPI::RemotePointEvaluation (short: ``rpe'') and related +classes/functions. In this section, basic functionalities are briefly +summarized. Their advanced capabilities will be shown subsequently based on +concrete application cases. + +The following code snippet shows the setup steps for the communication pattern: +@code +std::vector> points; // ... (filling of points not shown) + +RemotePointEvaluation rpe; +rpe.reinit(points, triangulation, mapping); +@endcode + +All what is needed is a list of evaluation points and the mesh with a mapping. + +The following code snippet shows the evaluation steps: +@code +const std::function &, const CellData &)> + evaluation_function; + +std::vector output; +rpe.evaluate_and_process(output, evaluation_function); + +@endcode + +The user provides a function that processes the locally owned points. +These values are communicated by Utilities::MPI::RemotePointEvaluation. + +The relevant class during the local evaluation is +Utilities::MPI::RemotePointEvaluation::CellData. It allows to loop over +cells that surround the points. On these cells, +a cell iterator and the positions in the reference cell of the +requested points can be queried. Furthermore, this class provides +controlled access to the output vector of the +Utilities::MPI::RemotePointEvaluation::evaluate_and_process() function. +@code +for (const auto cell_index : cell_data.cell_indices()) + { + const auto cell = cell_data.get_active_cell_iterator(cell_index); + const auto unit_points = cell_data.get_unit_points(cell_index); + const auto local_output = cell_data.get_data_view(cell_index, output); + } +@endcode + +The functions +@code +const auto evaluated_values = + VectorTools::point_values(rpe, dof_handler, vector); + +const auto evaluated_gradients = + VectorTools::point_gradients(rpe, dof_handler, vector); +@endcode + +evaluate the values and gradients of a solution defined by DoFHandler and a +vector at the requested points. Internally, a lambda function is passed to +Utilities::MPI::RemotePointEvaluation. +Additionally it handles the special case if points belong to multiple cells +by taking, e.g., the average, the minimum, or +the maximum via an optional argument of type EvaluationFlags::EvaluationFlags. +This occurs when a point lies on a cell boundary or within a small tolerance +around it and might be relevant for discontinuous solution quantities, such +as values of discontinuous Galerkin methods or gradients in continuous finite +element methods. + + +

Motivation: two-phase flow

+ +The minimal code examples (short "mini examples") presented in this tutorial +are motivated by the application of two-phase-flow simulations formulated in +a one-fluid setting using a Eulerian framework. In diffuse interface methods, +the two phases may be implicitly described by a level-set function, here +chosen as a signed distance function $\phi(\boldsymbol{ x})$ in +$\Omega$ and illustrated for a popular benchmark case of a rising +bubble in the following figure. + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_87_rising_bubble.png +
+ +The discrete interface $\Gamma$ is +represented implicitly through a certain isosurface of the level-set function +e.g. for the signed-distance function $\Gamma=\{\boldsymbol x \in \Omega~|~\phi +(\boldsymbol{x})=0\}$. +We would like to note that deal.II provides a set of analytical signed distance +functions for simple geometries in the Functions::SignedDistance namespace. +Those can be combined via Boolean operations to describe more complex +geometries @cite burman2015cutfem. The temporal evolution of the level-set field +is obtained by the transport equation +@f[ + \frac{\partial \, \phi}{\partial \, t} + \boldsymbol{u}\vert_\Gamma \cdot + \nabla \phi = 0 +@f] +with the transport velocity at the interface $\boldsymbol{u}\vert_\Gamma$, +which might be approximated by the local fluid velocity $\boldsymbol{u} +\vert_\Gamma\approx\boldsymbol{u}(\boldsymbol{x})$. To reobtain the +signed-distance property of the level-set field throughout the numerical +solution procedure, PDE-based or, alternatively, also geometric reinitialization +methods are used. For the latter, an algorithm for computing the distance from +the support points to the discrete interface, e.g., via closest-point point +projection @cite henri2022geometrical, is needed. This will be part of one of +the mini examples, where we describe how to obtain the closest point +$\boldsymbol{x}^*$ to the interface $\Gamma$ for an arbitrary point +$\boldsymbol{x}$. For the simplest case, the former can be computed from the +following equation +@f[ + \boldsymbol{x}^{*} = \boldsymbol{x} - \boldsymbol{n}(\boldsymbol{x}) + \phi(\boldsymbol{x}), +@f] +assuming that the interface normal vector $\boldsymbol{n}(\boldsymbol{x})$ and +$\phi(\boldsymbol{x})$ represent exact quantities. Typically, this projection +is only performed for points located within a narrow band region around the +interface, indicated in the right panel of the figure above. + +Alternatively to the implicit representation of the interface, in sharp +interface methods, e.g., via front tracking, the interface $\Gamma$ is +explicitly represented by a surface mesh. The latter is immersed into a +background mesh, from which the local velocity at the support points of the +surface mesh is extracted and leads to a movement of the support points of the +immersed mesh as +@f[ + \boldsymbol{x}_q^{(i + 1)} = \boldsymbol{x}_q^{(i)} + \Delta t \cdot + \boldsymbol{u}(\boldsymbol{x}_q^{(i)}) \quad \text{ for } \boldsymbol{x}_q + \in \Gamma +@f] +which considers an explicit Euler time integration scheme from time step $i$ to +$i+1$ with time step-size $\Delta t$. + +For a two-phase-flow model considering the incompressible Navier-Stokes +equations, the two phases are usually coupled by a singular surface-tension +force $\boldsymbol{F}_S$, which results, together with the difference in fluid +properties, in discontinuities across the interface: +@f[ + \boldsymbol{F}_S(\boldsymbol{x})= \sigma \kappa(\boldsymbol{x}) + \boldsymbol{n}(\boldsymbol{x}) \delta_{\Gamma}(\boldsymbol{x}). +@f] +Here $\sigma$ represents the surface-tension coefficient, +$\boldsymbol{n}(\boldsymbol{x})$ the interface normal vector +and $\kappa(\boldsymbol{x})$ the interface mean curvature field. +The singularity at the interface is imposed by the Dirac delta function +@f[ +\delta_{\Gamma}(\boldsymbol{x}) = \begin{cases} +1 & \text{on } \Gamma \\ +0 & \text{else}\end{cases} +@f] +with support on the interface $\Gamma$. +In a finite element context, the weak form of the surface-tension force +is needed. The latter can be applied as a sharp surface-tension force model +@f[ + (\boldsymbol v, \boldsymbol F_S)_{\Omega} = \left( \boldsymbol{v}, \sigma + \kappa \boldsymbol{n} \right)_\Gamma, +@f] +exploiting the property of the Dirac delta function for any smooth +function $v$, i.e., +$\int_\Omega\delta_{\Gamma}\,v\,\text{d}x=\int_\Gamma v\,\text{d}y$. For +front-tracking methods, the curvature and the normal vector are directly +computed from the surface mesh. + +Alternatively, in regularized surface-tension-force models +@cite brackbill1992continuum @cite olsson2005conservative +@cite kronbichler2018fast, the Dirac delta function is approximated by a smooth +ansatz +@f[ +(\boldsymbol v, \boldsymbol F_S)_{\Omega} \approx \left(\boldsymbol v, \sigma +\kappa \boldsymbol{n} \|\nabla H\| \right)_\Omega +@f] +considering the absolute value of the gradient of a regularized indicator +function $\|\nabla H\|$, which is related to the level-set field. In such +models, the interface normal vector +@f[ + \boldsymbol{n}(\boldsymbol{x}) = \nabla \phi(\boldsymbol{x}), +@f] +and the interface curvature field +@f[ + \kappa(\boldsymbol{x}) = \nabla \cdot \boldsymbol{n}(\boldsymbol{x})= + \Delta \phi(\boldsymbol{x}) \,. +@f] +are derived from the level-set function. + + +

Overview

+ +In the following, we present three simple use cases of +Utilities::MPI::RemotePointEvaluation. +We start with discussing a serial code in mini example 0. +In the subsequent mini examples, advanced problems are solved on distributed +meshes: +
    +
  • mini example 1: we evaluate values and user quantities along a line;
    • +
  • mini example 2: we perform a closest-point projection within a narrow band, +based on a level-set function, use the information to update the distance and +to perform an extrapolation from the interface;
    • +
  • mini example 3: we compute the surface-tension term sharply +with the interface given by an codim-1 mesh, which is advected by +the velocity from the background mesh (front tracking; +solution transfer between a background mesh and an immersed surface mesh).
    • +
diff --git a/examples/step-87/doc/kind b/examples/step-87/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-87/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-87/doc/results.dox b/examples/step-87/doc/results.dox new file mode 100644 index 0000000000..7f6679681c --- /dev/null +++ b/examples/step-87/doc/results.dox @@ -0,0 +1,272 @@ +

Results

+ +

Mini example 0

+ +We present a part of the terminal output. It shows, for each point, the +determined cell and reference position. Also, one can see that +the values evaluated with FEValues, FEPointEvaluation, and +VectorTools::point_values() are identical, as expected. + +@verbatim +Running: example 0 + - Found point with real coordinates: 0 0.5 + - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6) + - with coordinates on the unit cell: (0 0.5) + - Values at point: + - 0.25002 (w. FEValues) + - 0.25002 (w. FEPointEvaluation) + - 0.25002 (w. VectorTools::point_value()) + + - Found point with real coordinates: 0.05 0.5 + - in cell with vertices: (0 0.4) (0.2 0.4) (0 0.6) (0.2 0.6) + - with coordinates on the unit cell: (0.25 0.5) + - Values at point: + - 0.20003 (w. FEValues) + - 0.20003 (w. FEPointEvaluation) + - 0.20003 (w. VectorTools::point_value()) + +... + + - Found point with real coordinates: 1 0.5 + - in cell with vertices: (0.8 0.4) (1 0.4) (0.8 0.6) (1 0.6) + - with coordinates on the unit cell: (1 0.5) + - Values at point: + - 0.25002 (w. FEValues) + - 0.25002 (w. FEPointEvaluation) + - 0.25002 (w. VectorTools::point_value()) + + - writing csv file +@endverbatim + +The CSV output is as follows and contains, in the +first column, the distances with respect to the first point, +the second and the third column represent the coordinates +of the points and the fourth column the evaluated solution +values at those points. + +@verbatim +0.000 0.000 0.500 0.250 +0.050 0.050 0.500 0.200 +0.100 0.100 0.500 0.150 +0.150 0.150 0.500 0.100 +0.200 0.200 0.500 0.050 +0.250 0.250 0.500 0.000 +0.300 0.300 0.500 -0.050 +0.350 0.350 0.500 -0.100 +0.400 0.400 0.500 -0.149 +0.450 0.450 0.500 -0.200 +0.500 0.500 0.500 -0.222 +0.550 0.550 0.500 -0.200 +0.600 0.600 0.500 -0.149 +0.650 0.650 0.500 -0.100 +0.700 0.700 0.500 -0.050 +0.750 0.750 0.500 0.000 +0.800 0.800 0.500 0.050 +0.850 0.850 0.500 0.100 +0.900 0.900 0.500 0.150 +0.950 0.950 0.500 0.200 +1.000 1.000 0.500 0.250 +@endverbatim + +

Mini example 1

+ +We show the terminal output. + +@verbatim +Running: example 1 + - writing csv file +@endverbatim + +The CSV output is as follows and identical to the results +of the serial case presented in mini example 0. +The fifth column shows the +user quantity evaluated additionally in this mini example. + +@verbatim +0.000 0.000 0.500 0.250 0.000 +0.050 0.050 0.500 0.200 0.050 +0.100 0.100 0.500 0.150 0.100 +0.150 0.150 0.500 0.100 0.150 +0.200 0.200 0.500 0.050 0.200 +0.250 0.250 0.500 0.000 0.250 +0.300 0.300 0.500 -0.050 0.300 +0.350 0.350 0.500 -0.100 0.350 +0.400 0.400 0.500 -0.149 0.400 +0.450 0.450 0.500 -0.200 0.450 +0.500 0.500 0.500 -0.222 0.500 +0.550 0.550 0.500 -0.200 0.550 +0.600 0.600 0.500 -0.149 0.600 +0.650 0.650 0.500 -0.100 0.650 +0.700 0.700 0.500 -0.050 0.700 +0.750 0.750 0.500 0.000 0.750 +0.800 0.800 0.500 0.050 0.800 +0.850 0.850 0.500 0.100 0.850 +0.900 0.900 0.500 0.150 0.900 +0.950 0.950 0.500 0.200 0.950 +1.000 1.000 0.500 0.250 1.000 +@endverbatim + + +

Mini example 2

+ +We show the terminal output. +@verbatim +Running: example 2 + - create system + - determine narrow band + - determine closest point iteratively + - iteration 0: 7076 -> 7076 + - iteration 1: 7076 -> 104 + - iteration 2: 104 -> 0 + - determine distance in narrow band + - perform extrapolation in narrow band + - output results +@endverbatim + +The following three plots, representing the performed iterations of the +closest-point projection, show the current position of the closest +points exceeding the required tolerance of the discrete interface +of the circle and still need to +be corrected. +It can be seen that the numbers of points to be processed decrease +from iteration to iteration. + + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_0.png + + @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_1.png + + @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_p_2.png +
+ +The output visualized in Paraview looks like the following: On the +left, the original distance function is shown as the light gray surface. +In addition, the contour values refer to the distance values determined +from calculation of the distance to the closest points at the interface +in the narrow band. It can be seen that the two functions coincide. +Similarly, on the right, the original solution and the extrapolated +solution from the interface is shown. + + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_0.png + + @image html https://www.dealii.org/images/steps/developer/step_87_ex_2_res_1.png +
+ +

Mini example 3

+ +We show a shortened version of the terminal output. + +@verbatim +Running: example 3 + - creating background mesh + - creating immersed mesh +time: 0 + - compute to be tested values (immersed mesh) + - test values (background mesh) + - write data (background mesh) + - write mesh (immersed mesh) + +time: 0.01 + - move support points (immersed mesh) + - compute to be tested values (immersed mesh) + - test values (background mesh) + +time: 0.02 + - move support points (immersed mesh) + - compute to be tested values (immersed mesh) + - test values (background mesh) + +... + +time: 2 + - move support points (immersed mesh) + - compute to be tested values (immersed mesh) + - test values (background mesh) + - write data (background mesh) + - write mesh (immersed mesh) +@endverbatim + +The output visualized in Paraview looks like the following: The deformation of +the immersed mesh by the reversible vortex flow can be seen. Due to +discretization errors, the shape is not exactly circular at the end, illustrated +in the right figure. The sharp nature of the surface-tension force vector, shown +as vector plots, can be seen by its restriction to cells that are intersected by +the immersed mesh. + + + + + + + +
+ @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0000.png + + @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0010.png + + @image html https://www.dealii.org/images/steps/developer/step_87_ex_3_force.0020.png +
+ +

Possibilities for extension

+ +This program highlights some of the main capabilities +of the distributed evaluation routines in deal.II. However, there are many +related topics worth mentioning: +- Performing a distributed search is an expensive step. That is why we suggest +to provide hints to Utilities::MPI::RemotePointEvaluation and to reuse +Utilities::MPI::RemotePointEvaluation +instances in the case that the communication pattern has not changed. +Furthermore, there are instances where no search is needed and the points are +already sorted into the right cells. This is the case if the points are +generated on the cell level (see step-85; CutFEM) or the points are +automatically sorted into the correct (neighboring) cell (see step-68; PIC with +Particles::ParticleHandler). Having said that, the +Particles::ParticleHandler::insert_global_particles() function uses +the described infrastructure to perform the initial sorting of particles into +cells. +- We concentrated on parallelization aspects in this tutorial. However, we would +like to point out the need for fast evaluation on cell level. +The task for this in deal.II is FEPointEvaluation. It exploits the structure of +@f[ +\hat{u}(\hat{\boldsymbol{x}}) = \sum_i \hat{N}_i(\hat{\boldsymbol{x}}) \hat{u}_i +@f] +to derive fast implementations, e.g., for tensor-product elements +@f[ +\hat{u}(\hat{x}_0, \hat{x}_1, \hat{x}_2) = +\sum_k \hat{N}^{\text{1D}}_k(\hat{x}_2) +\sum_j \hat{N}^{\text{1D}}_j(\hat{x}_1) +\sum_i \hat{N}^{\text{1D}}_i(\hat{x}_0) +\hat{u}_{ijk}. +@f] +Since only 1D shape functions are queried and are re-used in re-occurring terms, +this formulation is faster than without exploitation of the structure. +- Utilities::MPI::RemotePointEvaluation is used in multiple places in deal.II. +The class DataOutResample allows to output results on a different mesh than +the computational mesh. This is useful if one wants to output the results +on a coarser mesh or one does not want to output 3D results but instead 2D +slices. In addition, MGTwoLevelTransferNonNested allows to prolongate solutions +and restrict residuals between two independent meshes. By passing a sequence +of such two-level transfer operators to MGTransferMF and, finally, to Multigrid, +non-nested multigrid can be computed. +- Utilities::MPI::RemotePointEvaluation can be used to couple non-matching +grids via surfaces (example: fluid-structure interaction, independently created +grids). The evaluation points can come from any side (pointwise interpolation) +or are defined on intersected meshes (Nitsche-type mortaring +@cite heinz2022high). Concerning the creation of such intersected meshes and the +corresponding evaluation points, see +GridTools::internal::distributed\_compute_intersection_locations(). +- Alternatively to the coupling via Utilities::MPI::RemotePointEvaluation, +preCICE @cite bungartz2016precice @cite chourdakis2021precice can be used. The +code-gallery program "Laplace equation coupled to an external simulation +program" shows how to use this library with deal.II. diff --git a/examples/step-87/doc/tooltip b/examples/step-87/doc/tooltip new file mode 100644 index 0000000000..aba8ccdf88 --- /dev/null +++ b/examples/step-87/doc/tooltip @@ -0,0 +1 @@ +Evaluation of finite element solutions at arbitrary points within a distributed mesh with application to two-phase flow. diff --git a/examples/step-87/step-87.cc b/examples/step-87/step-87.cc new file mode 100644 index 0000000000..cf6fad791a --- /dev/null +++ b/examples/step-87/step-87.cc @@ -0,0 +1,1252 @@ +/* --------------------------------------------------------------------- + * + * 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: Magdalena Schreter-Fleischhacker, Technical University of + * Munich, 2023 + * Peter Munch, University of Augsburg, 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 + +#include +#include + +// The files most relevant for this tutorial are the ones that +// contain Utilities::MPI::RemotePointEvaluation and the distributed evaluation +// functions in the VectorTools namespace, which use +// Utilities::MPI::RemotePointEvaluation. +#include +#include + +// The following header file provides the class FEPointEvaluation, which allows +// us to evaluate values of a local solution vector at arbitrary unit points of +// a cell. +#include + +// We pack everything that is specific for this program into a namespace +// of its own. + +namespace Step87 +{ + using namespace dealii; + +// @sect3{Utility functions (declaration)} +// +// In the following, we declare utility functions that are used +// in the mini examples below. You find the definitions at the end +// of the tutorial. +// +// The minimum requirement of this tutorial is MPI. If deal.II is built +// with p4est, we use parallel::distributed::Triangulation as +// distributed mesh. The class parallel::shared::Triangulation is +// used if deal.II is built without p4est or if the dimension of the +// triangulation is 1D, e.g., in the case of codim-1 meshes. +#ifdef DEAL_II_WITH_P4EST + template + using DistributedTriangulation = typename std::conditional< + dim == 1, + parallel::shared::Triangulation, + parallel::distributed::Triangulation>::type; +#else + template + using DistributedTriangulation = + parallel::shared::Triangulation; +#endif + + // A list of points along a line is created by definition of a + // start point @p p0, an end point @p p1, and the number of subdivisions + // @p n_subdivisions. + template + std::vector> + create_points_along_line(const Point &p0, + const Point &p1, + const unsigned int n_subdivisions); + + // A given list of @p points and the corresponding values @p values_0 + // and @p values_1 (optional) are printed column-wise to a file @p + // file_name. In addition, the first column represents the distance + // of the points from the first point. + template + void print_along_line(const std::string &file_name, + const std::vector> &points, + const std::vector &values_0, + const std::vector &values_1 = {}); + + // Create a unique list of the real coordinates of support points into + // @p support_points from the provided Mapping @p mapping and the DoFHandler + // @p dof_handler. + template + void collect_support_points( + const Mapping &mapping, + const DoFHandler &dof_handler, + LinearAlgebra::distributed::Vector &support_points); + + // From the provided Mapping @p mapping and the DoFHandler @p dof_handler collect + // the global DoF indices and corresponding support points within a narrow + // band around the zero-level-set isosurface. Thereto, + // the value of the finite element function @p signed_distance corresponding to + // the DoFHandler @p dof_handler_support_points is evaluated at each support point. + // A support point is only collected if the absolute value is below the value + // for + // the @p narrow_band_threshold. + template + std::tuple>, std::vector> + collect_support_points_with_narrow_band( + const Mapping &mapping, + const DoFHandler &dof_handler_signed_distance, + const LinearAlgebra::distributed::Vector &signed_distance, + const DoFHandler &dof_handler_support_points, + const double narrow_band_threshold); + + // Convert a distributed vector of support points (@p support_points_unrolled) + // with a sequential order of the coordinates per point into a list of points. + template + std::vector> convert( + const LinearAlgebra::distributed::Vector &support_points_unrolled); + + // @sect3{Mini example 0: Evaluation at given points for a serial mesh} + // + // In this introductory example, we demonstrate basic functionalities + // available in deal.II to evaluate solution quantities at arbitrary + // points on a serial mesh. The same functionalities are used directly + // or indirecly in the distributed case to evaluate solution on locally + // owned cells. This, however, needs to be augmented by communication, + // as presented in following examples. + // + // We first create the typical objects needed + // for a finite element discretization (defined by mapping, triangulation, + // and finite element) and a vector containing finite element solution + // coefficients. + void example_0() + { + std::cout << "Running: example 0" << std::endl; + + constexpr unsigned int dim = 2; + constexpr unsigned int fe_degree = 3; + + MappingQ1 mapping; + Triangulation tria; + GridGenerator::subdivided_hyper_cube(tria, 7); + + FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + Vector vector(dof_handler.n_dofs()); + VectorTools::interpolate(mapping, + dof_handler, + Functions::SignedDistance::Sphere( + (dim == 1) ? Point(0.5) : + (dim == 2) ? Point(0.5, 0.5) : + Point(0.5, 0.5, 0.5), + 0.25), + vector); + + // We create a list of points inside the domain at which we + // would like to evaluate the finite element interpolant. + const auto points_line = + create_points_along_line((dim == 1) ? Point(0.0) : + (dim == 2) ? Point(0.0, 0.5) : + Point(0.0, 0.5, 0.5), + (dim == 1) ? Point(1.0) : + (dim == 2) ? Point(1.0, 0.5) : + Point(1.0, 0.5, 0.5), + 20); + + // Now, we loop over all evaluation points. In the first step, we determine + // via GridTools::find_active_cell_around_point() the cell $K$ that + // surrounds the point and translate the given real coordinate + // $\boldsymbol{x}$ to the corresponding coordinate on the unit cell + // $\hat{\boldsymbol{x}}_K$ according to the provided mapping. + // The resulting information is printed to the screen. + std::vector values_line; + values_line.reserve(points_line.size()); + + for (const auto &p_real : points_line) + { + const auto [cell, p_unit] = + GridTools::find_active_cell_around_point(mapping, + dof_handler, + p_real); + + { + AssertThrow(cell != dof_handler.end(), ExcInternalError()); + std::cout << " - Found point with real coordinates: " << p_real + << std::endl; + std::cout << " - in cell with vertices:"; + for (const auto &v : cell->vertex_indices()) + std::cout << " (" << cell->vertex(v) << ")"; + std::cout << std::endl; + std::cout << " - with coordinates on the unit cell: (" << p_unit + << ")" << std::endl; + } + + // Having determined $K$ and $\hat{\boldsymbol{x}}_K$, we can + // perform the evaluation of the finite element solution at this + // point. In the following, we show three approaches for this + // purpose. In the first approach, we follow a traditional technique + // by using FEValues based on a cell-specific quadrature rule + // consisting of the unit point. + std::cout << " - Values at point:" << std::endl; + + { + FEValues fe_values(mapping, + fe, + Quadrature(p_unit), + update_values); + fe_values.reinit(cell); + + std::vector quad_values(1); + fe_values.get_function_values(vector, quad_values); + const double value_0 = quad_values[0]; + std::cout << " - " << value_0 << " (w. FEValues)" << std::endl; + values_line.push_back(value_0); + } + + // The second approach considers FEPointEvaluation, which directly + // takes a list of unit points for the subsequent evaluation. + // The class FEPointEvaluation is a class optimized for the evaluation + // on cell level at arbitrary points and should be favored for such + // tasks. + { + std::vector cell_vector(fe.n_dofs_per_cell()); + cell->get_dof_values(vector, cell_vector.begin(), cell_vector.end()); + + FEPointEvaluation<1, dim> fe_point(mapping, fe, update_values); + fe_point.reinit(cell, ArrayView>(p_unit)); + fe_point.evaluate(cell_vector, EvaluationFlags::values); + const auto value_1 = fe_point.get_value(0); + std::cout << " - " << value_1 << " (w. FEPointEvaluation)" + << std::endl; + } + + // Finally, in the third approach, the function + // VectorTools::point_value() is considered. It performs both + // the search of the surrounding cell and the evaluation at the + // requested point. However, its application is limited + // to a serial run of the code. + { + const auto value_2 = + VectorTools::point_value(dof_handler, vector, p_real); + std::cout << " - " << value_2 << " (w. VectorTools::point_value())" + << std::endl; + std::cout << std::endl; + } + } + + // We output the requested points together with the corresponding + // evaluated solution to a CSV file. + std::cout << " - writing csv file" << std::endl; + print_along_line("example_0_profile.csv", points_line, values_line); + } + + // Obviously, the code above cannot work for distributed meshes, since + // the search (which might require communication) is called within a for-loop + // with loop bounds possibly different on each process. In the following + // code examples, we present the usage of arbitrary point evaluation + // in a parallel computation. + // + // @sect3{Mini example 1: Evaluation at given points on a distributed mesh} + // + // Just like in the introductory example, we evaluate the solution + // along a line, however, on a distributed mesh. We again start with + // setting up the objects needed for a finite element discretization. + void example_1() + { + constexpr unsigned int dim = 2; + constexpr unsigned int fe_degree = 3; + + ConditionalOStream pcout(std::cout, + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == + 0); + + pcout << "Running: example 1" << std::endl; + + MappingQ1 mapping; + DistributedTriangulation tria(MPI_COMM_WORLD); + GridGenerator::subdivided_hyper_cube(tria, 7); + + FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + // We determine a finite element solution representing implicitly + // the geometry of a sphere with a radius of $r=0.25$ and the center at + // $(0.5,0.5)$ via a signed distance function. + LinearAlgebra::distributed::Vector signed_distance; + signed_distance.reinit(dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + MPI_COMM_WORLD); + + VectorTools::interpolate(mapping, + dof_handler, + Functions::SignedDistance::Sphere( + (dim == 1) ? Point(0.5) : + (dim == 2) ? Point(0.5, 0.5) : + Point(0.5, 0.5, 0.5), + 0.25), + signed_distance); + + // Next, we fill a vector from an arbitrary function that should represent + // a possible finite element solution, which we would like to evaluate. + LinearAlgebra::distributed::Vector solution; + solution.reinit(dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + MPI_COMM_WORLD); + + VectorTools::interpolate(mapping, + dof_handler, + Functions::SignedDistance::Plane( + Point(), Point::unit_vector(0)), + solution); + + // We create a list of arbitrary (evaluation) points along a horizontal + // line, which intersects the center of the sphere. We do this only + // on the root rank, since we intend to output the results to a CSV file + // by the root rank. + std::vector> profile; + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + profile = create_points_along_line((dim == 1) ? Point(0.0) : + (dim == 2) ? Point(0.0, 0.5) : + Point(0.0, 0.5, 0.5), + (dim == 1) ? Point(1.0) : + (dim == 2) ? Point(1.0, 0.5) : + Point(1.0, 0.5, 0.5), + 20); + + // Now, we can evaluate the results, e.g., for the signed distance + // function at all evaluation points in one go. First, we create a + // modifiable Utilities::MPI::RemotePointEvaluation object. We use + // VectorTools::point_values() by specifying the number of components of the + // solution vector (1 for the present example) as a template parameter. + // Within this function, the provided object for + // Utilities::MPI::RemotePointEvaluation is automatically reinitialized with + // the given points (profile). The ghost values of the solution vector + // need to be updated from the user. + Utilities::MPI::RemotePointEvaluation rpe; + + signed_distance.update_ghost_values(); + const std::vector profile_signed_distance = + VectorTools::point_values<1>( + mapping, dof_handler, signed_distance, profile, rpe); + + // In addition to VectorTools::point_values(), function gradients can be + // evaluated via VectorTools::point_gradient(). However, for the computation + // of user-derived quantities, one might need to fall back to the direct + // usage of Utilities::MPI::RemotePointEvaluation::evaluate_and_process() or + // Utilities::MPI::RemotePointEvaluation::process_and_evaluate(). For the + // sake of demonstration, we use the former to evaluate the values + // of the solution vector at the requested points. First, we define a + // lambda function for the operation on the surrounding cells. Using the + // CellData object, we can create a FEPointEvaluation object to evaluate the + // solution values at the cell-specific unit coordinates of the requested + // points. Then, we assign the values to the result vector. + const auto evaluate_function = [&](const ArrayView &values, + const auto &cell_data) { + FEPointEvaluation<1, dim> fe_point(mapping, fe, update_values); + + std::vector local_values; + std::vector local_dof_indices; + + for (const auto cell : cell_data.cell_indices()) + { + const auto cell_dofs = + cell_data.get_active_cell_iterator(cell)->as_dof_handler_iterator( + dof_handler); + + const auto unit_points = cell_data.get_unit_points(cell); + const auto local_value = cell_data.get_data_view(cell, values); + + local_values.resize(cell_dofs->get_fe().n_dofs_per_cell()); + cell_dofs->get_dof_values(solution, + local_values.begin(), + local_values.end()); + + fe_point.reinit(cell_dofs, unit_points); + fe_point.evaluate(local_values, EvaluationFlags::values); + + for (unsigned int q = 0; q < unit_points.size(); ++q) + local_value[q] = fe_point.get_value(q); + } + }; + + // The lambda function is passed to + // Utilities::MPI::RemotePointEvaluation::evaluate_and_process(), where + // the values are processed accordingly and stored within the created + // output vector. Again, the ghost values of the vector to be read + // need to be updated by the user. + solution.update_ghost_values(); + + const std::vector output = + rpe.evaluate_and_process(evaluate_function); + + // Finally, we output all results: the mesh as a VTU file and the + // results along the line as a CSV file. You can import the CSV file in + // ParaView and compare the output with the native line plot of + // ParaView based on the VTU file. + DataOut data_out; + data_out.add_data_vector(dof_handler, signed_distance, "signed_distance"); + data_out.build_patches(mapping); + data_out.write_vtu_in_parallel("example_1.vtu", MPI_COMM_WORLD); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::cout << " - writing csv file" << std::endl; + print_along_line("example_1_profile.csv", + profile, + profile_signed_distance, + output); + } + } + + // @sect3{Mini example 2: Closest-point evaluation of a distributed mesh} + // + // In this mini example, we perform a closest-point projection for each + // support point of a mesh within a narrow band by iteratively solving for + // @f[ + //\boldsymbol{x}^{(i+1)} = \boldsymbol{x}^{(i)} - + //\boldsymbol{n}(\boldsymbol{x}^{(i)})\phi(\boldsymbol{x}^{(i)}). + // @f] + // Once the closest point is determined, we can compute the distance and + // extrapolate the values from the interface. Note that the demonstrated + // algorithm does not guarantee that the closest points are collinear + // (see discussion in @cite coquerelle2016fourth). For the latter, one + // might also need to perform a tangential correction, which we omit + // here to keep the discussion concise. + // + // We start with creating the objects for the finite element representation + // of the background mesh. + void example_2() + { + constexpr unsigned int dim = 2; + constexpr unsigned int fe_degree = 3; + + ConditionalOStream pcout(std::cout, + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == + 0); + + pcout << "Running: example 2" << std::endl; + pcout << " - create system" << std::endl; + + FE_Q fe(fe_degree); + MappingQ1 mapping; + DistributedTriangulation tria(MPI_COMM_WORLD); + GridGenerator::subdivided_hyper_cube(tria, 50); + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + // We compute finite element solution vector, + // based on an arbitrary function. In addition, a finite element + // function computed from a signed distance function represents + // the geometry of a sphere implicitly. + LinearAlgebra::distributed::Vector solution; + solution.reinit(dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + MPI_COMM_WORLD); + + VectorTools::interpolate(mapping, + dof_handler, + Functions::SignedDistance::Plane( + Point(), Point::unit_vector(0)), + solution); + + LinearAlgebra::distributed::Vector signed_distance; + signed_distance.reinit(dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + MPI_COMM_WORLD); + + VectorTools::interpolate(mapping, + dof_handler, + Functions::SignedDistance::Sphere( + (dim == 1) ? Point(0.5) : + (dim == 2) ? Point(0.5, 0.5) : + Point(0.5, 0.5, 0.5), + 0.25), + signed_distance); + signed_distance.update_ghost_values(); + + // In the next step, we collect the points in the narrow band around + // the zero-level-set isosurface for which we would like to perform + // a closest point projection. To this end, we loop over all support + // points and collect the coordinates and the DoF indices of those + // with a maximum distance of 0.1 from the zero-level-set isosurface. + pcout << " - determine narrow band" << std::endl; + + const auto [support_points, support_points_idx] = + collect_support_points_with_narrow_band(mapping, + dof_handler, + signed_distance, + dof_handler, + 0.1 /*narrow_band_threshold*/); + + // For the iterative solution procedure of the closest-point projection, + // the maximum number of iterations and the tolerance for the maximum + // absolute acceptable change in the distance in one iteration are set. + pcout << " - determine closest point iteratively" << std::endl; + constexpr int max_iter = 30; + constexpr double tol_distance = 1e-6; + + // Now, we are ready to perform the algorithm by setting an initial guess + // for the projection points simply corresponding to the collected support + // points. We collect the global indices of the support points and the + // total number of points that need to be processed and do not + // fulfill the required tolerance. Those will be gradually reduced + // upon the iterative process. + std::vector> closest_points = support_points; // initial guess + + std::vector unmatched_points_idx(closest_points.size()); + std::iota(unmatched_points_idx.begin(), unmatched_points_idx.end(), 0); + + int n_unmatched_points = + Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD); + + // Now, we create a Utilities::MPI::RemotePointEvaluation cache object and + // start the loop for the fix-point iteration. We update the list of points + // that still need to be processed and subsequently pass this information + // to the Utilities::MPI::RemotePointEvaluation object. For the sake of + // illustration, we export the coordinates of the points to be updated for + // each iteration to a CSV file. Next, we can evaluate the signed distance + // function and the gradient at those points to update the current solution + // for the closest points. We perform the update only if the signed + // distance of the closest point is not already within the tolerance + // and register those points that still need to be processed. + Utilities::MPI::RemotePointEvaluation rpe; + + for (int it = 0; it < max_iter && n_unmatched_points > 0; ++it) + { + pcout << " - iteration " << it << ": " << n_unmatched_points; + + std::vector> unmatched_points(unmatched_points_idx.size()); + for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i) + unmatched_points[i] = closest_points[unmatched_points_idx[i]]; + + const auto all_unmatched_points = + Utilities::MPI::reduce>>( + unmatched_points, MPI_COMM_WORLD, [](const auto &a, const auto &b) { + auto result = a; + result.insert(result.end(), b.begin(), b.end()); + return result; + }); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::ofstream file("example_2_" + std::to_string(it) + ".csv"); + for (const auto &p : all_unmatched_points) + file << p << std::endl; + file.close(); + } + + rpe.reinit(unmatched_points, tria, mapping); + + AssertThrow(rpe.all_points_found(), + ExcMessage("Processed point is outside domain.")); + + const auto eval_values = + VectorTools::point_values<1>(rpe, dof_handler, signed_distance); + + const auto eval_gradient = + VectorTools::point_gradients<1>(rpe, dof_handler, signed_distance); + + std::vector unmatched_points_idx_next; + + for (unsigned int i = 0; i < unmatched_points_idx.size(); ++i) + if (std::abs(eval_values[i]) > tol_distance) + { + closest_points[unmatched_points_idx[i]] -= + eval_values[i] * eval_gradient[i]; + + unmatched_points_idx_next.emplace_back(unmatched_points_idx[i]); + } + + unmatched_points_idx.swap(unmatched_points_idx_next); + + n_unmatched_points = + Utilities::MPI::sum(unmatched_points_idx.size(), MPI_COMM_WORLD); + + pcout << " -> " << n_unmatched_points << std::endl; + } + + // We print a warning message if we exceed the maximum number of allowed + // iterations and if there are still projection points with a distance + // value exceeding the tolerance. + if (n_unmatched_points > 0) + pcout << "WARNING: The tolerance of " << n_unmatched_points + << " points is not yet attained." << std::endl; + + // As a result, we obtain a list of support points and corresponding + // closest points at the zero-isosurface level set. This information + // can be used to update the signed distance function, i.e., the + // reinitialization the values of the level-set function to maintain + // the signed distance property @cite henri2022geometrical. + pcout << " - determine distance in narrow band" << std::endl; + LinearAlgebra::distributed::Vector solution_distance; + solution_distance.reinit(solution); + + for (unsigned int i = 0; i < closest_points.size(); ++i) + solution_distance[support_points_idx[i]] = + support_points[i].distance(closest_points[i]); + + // In addition, we use the information of the closest point to + // extrapolate values from the interface, i.e., the zero-level + // set isosurface, to the support points within the narrow band. + // This might be helpful to improve accuracy, e.g., for + // diffuse interface fluxes where certain quantities are only + // accurately determined at the interface (e.g. curvature + // for surface tension @cite coquerelle2016fourth). + pcout << " - perform extrapolation in narrow band" << std::endl; + rpe.reinit(closest_points, tria, mapping); + solution.update_ghost_values(); + const auto vals = VectorTools::point_values<1>(rpe, dof_handler, solution); + + LinearAlgebra::distributed::Vector solution_extrapolated; + solution_extrapolated.reinit(solution); + + for (unsigned int i = 0; i < closest_points.size(); ++i) + solution_extrapolated[support_points_idx[i]] = vals[i]; + + // Finally, we output the results to a VTU file. + pcout << " - output results" << std::endl; + DataOut data_out; + data_out.add_data_vector(dof_handler, signed_distance, "signed_distance"); + data_out.add_data_vector(dof_handler, solution, "solution"); + data_out.add_data_vector(dof_handler, + solution_distance, + "solution_distance"); + data_out.add_data_vector(dof_handler, + solution_extrapolated, + "solution_extrapolated"); + data_out.build_patches(mapping); + data_out.write_vtu_in_parallel("example_2.vtu", MPI_COMM_WORLD); + + pcout << std::endl; + } + + // @sect3{Mini example 3: Sharp interface method on the example of surface tension for front tracking} + // + // The final mini example presents a basic implementation of + // front tracking @cite peskin1977numerical, @cite unverdi1992front + // of a surface mesh $\mathbb{T}_\Gamma$ immersed + // in a Eulerian background fluid mesh $\mathbb{T}_\Omega$. + // + // We assume that the immersed surface is transported according to a + // prescribed velocity field from the background mesh. Subsequently, + // we perform a sharp computation of the surface-tension force: + // @f[ + // (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S + // (\boldsymbol{x}))_{\Omega} + // = + // \left( \boldsymbol{v}_i (\boldsymbol{x}), \sigma (\boldsymbol{x}) \kappa + // (\boldsymbol{x}) \boldsymbol{n} (\boldsymbol{x}) \right)_\Gamma \approx + // \sum_{q\in\mathbb{T}_\Gamma} \boldsymbol{v}_i^T (\boldsymbol{x}_q) + // \sigma (\boldsymbol{x}_q) \kappa (\boldsymbol{x}_q) \boldsymbol{n} + // (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)| w(\boldsymbol{x}_q) \quad \forall + // i\in\mathbb{T}_\Omega + // . + // @f] + // We decompose this operation into two steps. In the first step, we evaluate + // the force contributions $\sigma (\boldsymbol{x}_q) \kappa + // (\boldsymbol{x}_q) \boldsymbol{n} + // (\boldsymbol{x}_q)$ at the quadrature points defined on the immersed mesh + // and multiply them with the mapped quadrature weight $|J(\boldsymbol{x}_q)| + // w_q$: + // @f[ + // \boldsymbol{F}_S (\boldsymbol{x}_q) \gets \sigma (\boldsymbol{x}_q) \kappa + // (\boldsymbol{x}_q) \boldsymbol{n} (\boldsymbol{x}_q) |J(\boldsymbol{x}_q)| + // w_q \quad \forall q\in\mathbb{T}_\Gamma. + // @f] + // In the second step, we compute the discretized weak form by multiplying + // with test functions on the background mesh: + // @f[ + // (\boldsymbol v_i (\boldsymbol{x}), \boldsymbol F_S + // (\boldsymbol{x}))_{\Omega} \gets \sum_q \boldsymbol{v}_i^T + // (\boldsymbol{x}_q) \boldsymbol{F}_S + // (\boldsymbol{x}_q) + // \quad \forall i\in\mathbb{T}_\Omega + // . + // @f] + // Obviously, we need to communicate between the two steps. The second step + // can be handled completely by Utilities::MPI::RemotePointEvaluation, which + // provides the function + // Utilities::MPI::RemotePointEvaluation::process_and_evaluate() for this + // purpose. + // + // We start with setting the parameters consisting of the polynomial degree of + // the shape functions, the dimension of the background mesh, the time-step + // size to be considered for transporting the surface mesh and the number of + // time steps. + + void example_3() + { + constexpr unsigned int degree = 3; + constexpr unsigned int dim = 2; + const double dt = 0.01; + const unsigned int n_time_steps = 200; + + // This program is intended to be executed in 2D or 3D. + static_assert(dim == 2 || dim == 3, "Only implemented for 2D or 3D."); + + ConditionalOStream pcout(std::cout, + Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == + 0); + + pcout << "Running: example 3" << std::endl; + + // Next, we create the standard objects necessary for the finite element + // representation of the background mesh + pcout << " - creating background mesh" << std::endl; + DistributedTriangulation tria_background(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria_background); + tria_background.refine_global(5); + + MappingQ1 mapping_background; + FESystem fe_background(FE_Q(degree), dim); + DoFHandler dof_handler_background(tria_background); + dof_handler_background.distribute_dofs(fe_background); + + // and, similarly, for the immersed surface mesh. + // We use a sphere with radius $r=0.75$ which is + // placed in the center of the top half of the cubic background domain. + pcout << " - creating immersed mesh" << std::endl; + const Point center((dim == 2) ? Point(0.5, 0.75) : + Point(0.5, 0.75, 0.5)); + const double radius = 0.15; + + DistributedTriangulation tria_immersed(MPI_COMM_WORLD); + GridGenerator::hyper_sphere(tria_immersed, center, radius); + tria_immersed.refine_global(4); + + // Two different mappings are considered for the immersed + // surface mesh: one valid for the initial configuration and one + // that is updated in every time step according to the nodal + // displacements. Two types of finite elements are used to + // represent scalar and vector-valued DoF values. + MappingQ mapping_immersed_base(3); + MappingQCache mapping_immersed(3); + mapping_immersed.initialize(mapping_immersed_base, tria_immersed); + QGauss quadrature_immersed(degree + 1); + + FE_Q fe_scalar_immersed(degree); + FESystem fe_immersed(fe_scalar_immersed, dim); + DoFHandler dof_handler_immersed(tria_immersed); + dof_handler_immersed.distribute_dofs(fe_immersed); + + // We renumber the DoFs related to the vector-valued problem to + // simplify access to the individual components. + DoFRenumbering::support_point_wise(dof_handler_immersed); + + // We fill a DoF vector on the background mesh with an analytical + // velocity field considering the Rayleigh-Kothe vortex flow and + // initialize a DoF vector for the weak form of the surface-tension force. + LinearAlgebra::distributed::Vector velocity; + velocity.reinit(dof_handler_background.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs( + dof_handler_background), + MPI_COMM_WORLD); + Functions::RayleighKotheVortex vortex(2); + + LinearAlgebra::distributed::Vector force_vector( + dof_handler_background.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler_background), + MPI_COMM_WORLD); + + // Next, we collect the real positions $\boldsymbol{x}_q$ of the quadrature + // points of the surface mesh in a vector. + LinearAlgebra::distributed::Vector immersed_support_points; + collect_support_points(mapping_immersed, + dof_handler_immersed, + immersed_support_points); + + // We initialize a Utilities::MPI::RemotePointEvaluation object and start + // the time loop. For any other step than the initial one, we first move the + // support points of the surface mesh according to the fluid velocity of the + // background mesh. Thereto, we first update the time of the velocity + // function. Then, we update the internal data structures of the + // Utilities::MPI::RemotePointEvaluation object with the collected support + // points of the immersed mesh. We throw an exception if one of the points + // cannot be found within the domain of the background mesh. Next, we + // evaluate the velocity at the surface-mesh support points and compute the + // resulting update of the coordinates. Finally, we update the mapping of + // the immersed surface mesh to the current position. + Utilities::MPI::RemotePointEvaluation rpe; + double time = 0.0; + for (unsigned int it = 0; it <= n_time_steps; ++it, time += dt) + { + pcout << "time: " << time << std::endl; + if (it > 0) + { + pcout << " - move support points (immersed mesh)" << std::endl; + vortex.set_time(time); + VectorTools::interpolate(mapping_background, + dof_handler_background, + vortex, + velocity); + rpe.reinit(convert(immersed_support_points), + tria_background, + mapping_background); + + AssertThrow(rpe.all_points_found(), + ExcMessage( + "Immersed domain leaves background domain!")); + + velocity.update_ghost_values(); + const auto immersed_velocity = + VectorTools::point_values(rpe, + dof_handler_background, + velocity); + + for (unsigned int i = 0, c = 0; + i < immersed_support_points.locally_owned_size() / dim; + ++i) + for (unsigned int d = 0; d < dim; ++d, ++c) + immersed_support_points.local_element(c) += + immersed_velocity[i][d] * dt; + + mapping_immersed.initialize(mapping_immersed_base, + dof_handler_immersed, + immersed_support_points, + false); + } + + // Next, we loop over all locally owned cells of the immersed mesh and + // its quadrature points to compute the value for the local surface + // tension force contribution $\boldsymbol{F}_S(\boldsymbol{x}_q)$. We + // store the real coordinates of the quadrature points and the + // corresponding force contributions in two individual vectors. For + // computation of the latter, the normal vector + // $\boldsymbol{n}(\boldsymbol{x}_q)$ can be directly extracted from the + // surface mesh via FEValues and, for the curvature, we use the + // following approximation: + // @f[ + // \kappa(\boldsymbol{x}_q) + // = + // \nabla \cdot \boldsymbol{n}(\boldsymbol{x}_q) + // = + // \text{tr}\left({\nabla \boldsymbol{n}(\boldsymbol{x}_q)}\right) + // \approx + // \text{tr}\left({\nabla \sum_i \boldsymbol{N}_i (\boldsymbol{x}_q) + // \boldsymbol n_i}\right) + // = + // \sum_i\text{tr}\left({\nabla \boldsymbol{N}_i (\boldsymbol{x}_q) + // \boldsymbol n_i}\right) + // \;\text{with}\; i\in[0,n_{\text{dofs_per_cell}}), + // @f] + // which we can apply since the immersed mesh is consistently + // orientated. The surface tension coefficient is set to 1 for the + // sake of demonstration. + pcout << " - compute to be tested values (immersed mesh)" << std::endl; + using value_type = Tensor<1, dim, double>; + + std::vector> integration_points; + std::vector integration_values; + + FEValues fe_values(mapping_immersed, + fe_immersed, + quadrature_immersed, + update_JxW_values | update_gradients | + update_normal_vectors | + update_quadrature_points); + + FEValues fe_values_co( + mapping_immersed, + fe_scalar_immersed, + fe_scalar_immersed.get_unit_support_points(), + update_JxW_values | update_normal_vectors); + + std::vector component_to_system_index( + fe_immersed.n_dofs_per_cell()); + + for (unsigned int i = 0, c = 0; + i < fe_scalar_immersed.n_dofs_per_cell(); + ++i) + for (unsigned int d = 0; d < dim; ++d, ++c) + component_to_system_index[c] = + fe_immersed.component_to_system_index(d, i); + + for (const auto &cell : tria_immersed.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + fe_values.reinit(cell); + fe_values_co.reinit(cell); + + for (const auto &q : fe_values.quadrature_point_indices()) + { + integration_points.emplace_back(fe_values.quadrature_point(q)); + + const auto sigma = 1.0; // surface tension coefficient + + const auto normal = fe_values.normal_vector(q); + double curvature = 0; + for (unsigned int i = 0, c = 0; + i < fe_scalar_immersed.n_dofs_per_cell(); + ++i) + for (unsigned int d = 0; d < dim; ++d, ++c) + curvature += fe_values.shape_grad_component( + component_to_system_index[c], q, d)[d] * + fe_values_co.normal_vector(i)[d]; + + const auto FxJxW = + sigma * curvature * normal * fe_values.JxW(q); + + integration_values.emplace_back(FxJxW); + } + } + + // Before we evaluate the weak form of the surface-tension force, the + // communication pattern of Utilities::MPI::RemotePointEvaluation is + // set up from the quadrature points of the immersed mesh, determining + // the surrounding cells on the background mesh. + pcout << " - test values (background mesh)" << std::endl; + + rpe.reinit(integration_points, tria_background, mapping_background); + + // In preparation for utilizing + // Utilities::MPI::RemotePointEvaluation::process_and_evaluate that + // performs the + // multiplication with the test function, we set up a callback function + // that contains the operation on the intersected cells of the + // background mesh. Within this function, we initialize a + // FEPointEvaluation object that allows us to integrate values at + // arbitrary points within a cell. We loop over the cells that surround + // quadrature points of the immersed mesh -- provided by the callback + // function. From the provided CellData object, we retrieve the unit + // points, i.e., the quadrature points of the immersed mesh that lie + // within the current cell and a pointer to the stored values on the + // current cell (local surface-tension force) for convenience. We + // reinitialize the data structure of FEPointEvaluation on every cell + // according to the unit points. Next, we loop over the quadrature + // points attributed to the cell and submit the local surface-tension + // force to the FEPointEvaluation object. Via + // FEPointEvaluation::test_and_sum(), the submitted values are + // multiplied by the values of the test function and a summation over + // all given points is performed. Subsequently, the contributions are + // assembled into the global vector containing the weak form of the + // surface-tension force. + const auto integration_function = [&](const auto &values, + const auto &cell_data) { + FEPointEvaluation phi_force(mapping_background, + fe_background, + update_values); + + std::vector local_values; + std::vector local_dof_indices; + + for (const auto cell : cell_data.cell_indices()) + { + const auto cell_dofs = + cell_data.get_active_cell_iterator(cell) + ->as_dof_handler_iterator(dof_handler_background); + + const auto unit_points = cell_data.get_unit_points(cell); + const auto FxJxW = cell_data.get_data_view(cell, values); + + phi_force.reinit(cell_dofs, unit_points); + + for (const auto q : phi_force.quadrature_point_indices()) + phi_force.submit_value(FxJxW[q], q); + + local_values.resize(cell_dofs->get_fe().n_dofs_per_cell()); + phi_force.test_and_sum(local_values, EvaluationFlags::values); + + local_dof_indices.resize(cell_dofs->get_fe().n_dofs_per_cell()); + cell_dofs->get_dof_indices(local_dof_indices); + AffineConstraints().distribute_local_to_global( + local_values, local_dof_indices, force_vector); + } + }; + + // The callback function is passed together with the vector holding the + // surface-tension force contribution at each quadrature point of the + // immersed mesh to + // Utilities::MPI::RemotePointEvaluation::process_and_evaluate. The only + // missing step is to compress the distributed force vector. + rpe.process_and_evaluate(integration_values, + integration_function); + force_vector.compress(VectorOperation::add); + + // After every tenth step or at the beginning/end of the time loop, we + // output the force vector and the velocity of the background mesh to + // a VTU file. In addition, we also export the geometry of the + // (deformed) immersed surface mesh to a separate VTU file. + if (it % 10 == 0 || it == n_time_steps) + { + std::vector< + DataComponentInterpretation::DataComponentInterpretation> + vector_component_interpretation( + dim, DataComponentInterpretation::component_is_part_of_vector); + pcout << " - write data (background mesh)" << std::endl; + DataOut data_out_background; + DataOutBase::VtkFlags flags_backround; + flags_backround.write_higher_order_cells = true; + data_out_background.set_flags(flags_backround); + data_out_background.add_data_vector( + dof_handler_background, + force_vector, + std::vector(dim, "force"), + vector_component_interpretation); + data_out_background.add_data_vector( + dof_handler_background, + velocity, + std::vector(dim, "velocity"), + vector_component_interpretation); + data_out_background.build_patches(mapping_background, 3); + data_out_background.write_vtu_in_parallel("example_3_background_" + + std::to_string(it) + + ".vtu", + MPI_COMM_WORLD); + + pcout << " - write mesh (immersed mesh)" << std::endl; + DataOut data_out_immersed; + data_out_immersed.attach_triangulation(tria_immersed); + data_out_immersed.build_patches(mapping_immersed, 3); + data_out_immersed.write_vtu_in_parallel("example_3_immersed_" + + std::to_string(it) + + ".vtu", + MPI_COMM_WORLD); + } + pcout << std::endl; + } + } + + // @sect3{Utility functions (definition)} + template + std::vector> + create_points_along_line(const Point &p0, + const Point &p1, + const unsigned int n_subdivisions) + { + Assert(n_subdivisions >= 1, ExcInternalError()); + + std::vector> points; + points.reserve(n_subdivisions + 1); + + points.emplace_back(p0); + for (unsigned int i = 1; i < n_subdivisions; ++i) + points.emplace_back(p0 + (p1 - p0) * static_cast(i) / + static_cast(n_subdivisions)); + points.emplace_back(p1); + + return points; + } + + template + void print_along_line(const std::string &file_name, + const std::vector> &points, + const std::vector &values_0, + const std::vector &values_1) + { + AssertThrow(points.size() == values_0.size() && + (values_1.size() == points.size() || values_1.empty()), + ExcMessage("The provided vectors must have the same length.")); + + std::ofstream file(file_name); + + for (unsigned int i = 0; i < points.size(); ++i) + { + file << std::fixed << std::right << std::setw(5) << std::setprecision(3) + << points[0].distance(points[i]); + + for (unsigned int d = 0; d < spacedim; ++d) + file << std::fixed << std::right << std::setw(10) + << std::setprecision(3) << points[i][d]; + + file << std::fixed << std::right << std::setw(10) + << std::setprecision(3) << values_0[i]; + + if (!values_1.empty()) + file << std::fixed << std::right << std::setw(10) + << std::setprecision(3) << values_1[i]; + file << std::endl; + } + } + + template + void collect_support_points( + const Mapping &mapping, + const DoFHandler &dof_handler, + LinearAlgebra::distributed::Vector &support_points) + { + support_points.reinit(dof_handler.locally_owned_dofs(), + DoFTools::extract_locally_active_dofs(dof_handler), + dof_handler.get_communicator()); + + const auto &fe = dof_handler.get_fe(); + + FEValues fe_values( + mapping, + fe, + fe.base_element(0).get_unit_support_points(), + update_quadrature_points); + + std::vector local_dof_indices( + fe.n_dofs_per_cell()); + + std::vector component_to_system_index( + fe_values.n_quadrature_points * spacedim); + + for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q) + for (unsigned int d = 0; d < spacedim; ++d, ++c) + component_to_system_index[c] = fe.component_to_system_index(d, q); + + for (const auto &cell : dof_handler.active_cell_iterators() | + IteratorFilters::LocallyOwnedCell()) + { + fe_values.reinit(cell); + cell->get_dof_indices(local_dof_indices); + + for (unsigned int q = 0, c = 0; q < fe_values.n_quadrature_points; ++q) + for (unsigned int d = 0; d < spacedim; ++d, ++c) + support_points[local_dof_indices[component_to_system_index[c]]] = + fe_values.quadrature_point(q)[d]; + } + } + + template + std::tuple>, std::vector> + collect_support_points_with_narrow_band( + const Mapping &mapping, + const DoFHandler &dof_handler_signed_distance, + const LinearAlgebra::distributed::Vector &signed_distance, + const DoFHandler &dof_handler_support_points, + const double narrow_band_threshold) + { + AssertThrow(narrow_band_threshold >= 0, + ExcMessage("The narrow band threshold" + " must be larger than or equal to 0.")); + const auto &tria = dof_handler_signed_distance.get_triangulation(); + const Quadrature quad(dof_handler_support_points.get_fe() + .base_element(0) + .get_unit_support_points()); + + FEValues distance_values(mapping, + dof_handler_signed_distance.get_fe(), + quad, + update_values); + + FEValues req_values(mapping, + dof_handler_support_points.get_fe(), + quad, + update_quadrature_points); + + std::vector temp_distance(quad.size()); + std::vector local_dof_indices( + dof_handler_support_points.get_fe().n_dofs_per_cell()); + + std::vector> support_points; + std::vector support_points_idx; + + const bool has_ghost_elements = signed_distance.has_ghost_elements(); + + const auto &locally_owned_dofs_req = + dof_handler_support_points.locally_owned_dofs(); + std::vector flags(locally_owned_dofs_req.n_elements(), false); + + if (has_ghost_elements == false) + signed_distance.update_ghost_values(); + + for (const auto &cell : + tria.active_cell_iterators() | IteratorFilters::LocallyOwnedCell()) + { + const auto cell_distance = + cell->as_dof_handler_iterator(dof_handler_signed_distance); + distance_values.reinit(cell_distance); + distance_values.get_function_values(signed_distance, temp_distance); + + const auto cell_req = + cell->as_dof_handler_iterator(dof_handler_support_points); + req_values.reinit(cell_req); + cell_req->get_dof_indices(local_dof_indices); + + for (const auto q : req_values.quadrature_point_indices()) + if (std::abs(temp_distance[q]) < narrow_band_threshold) + { + const auto idx = local_dof_indices[q]; + + if (locally_owned_dofs_req.is_element(idx) == false || + flags[locally_owned_dofs_req.index_within_set(idx)]) + continue; + + flags[locally_owned_dofs_req.index_within_set(idx)] = true; + + support_points_idx.emplace_back(idx); + support_points.emplace_back(req_values.quadrature_point(q)); + } + } + + if (has_ghost_elements == false) + signed_distance.zero_out_ghost_values(); + + return {support_points, support_points_idx}; + } + + template + std::vector> convert( + const LinearAlgebra::distributed::Vector &support_points_unrolled) + { + const unsigned int n_points = + support_points_unrolled.locally_owned_size() / spacedim; + + std::vector> points(n_points); + + for (unsigned int i = 0, c = 0; i < n_points; ++i) + for (unsigned int d = 0; d < spacedim; ++d, ++c) + points[i][d] = support_points_unrolled.local_element(c); + + return points; + } + +} // namespace Step87 + +// @sect3{Driver} +// +// Finally, the driver of the program executes the four mini examples. +int main(int argc, char **argv) +{ + using namespace dealii; + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + std::cout.precision(5); + + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + Step87::example_0(); // only run on root process + + Step87::example_1(); + Step87::example_2(); + Step87::example_3(); +}