]> https://gitweb.dealii.org/ - dealii.git/commitdiff
First commented version of step-60
authorLuca Heltai <luca.heltai@sissa.it>
Mon, 23 Apr 2018 15:44:31 +0000 (17:44 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Mon, 7 May 2018 12:51:59 +0000 (14:51 +0200)
Co-authored-by: Giovanni Alzetta <giovannialzetta@hotmail.it>
examples/step-60/CMakeLists.txt [new file with mode: 0644]
examples/step-60/doc/intro.dox [new file with mode: 0644]
examples/step-60/doc/kind [new file with mode: 0644]
examples/step-60/doc/results.dox [new file with mode: 0644]
examples/step-60/doc/tooltip [new file with mode: 0644]
examples/step-60/step-60.cc [new file with mode: 0644]

diff --git a/examples/step-60/CMakeLists.txt b/examples/step-60/CMakeLists.txt
new file mode 100644 (file)
index 0000000..99bfe40
--- /dev/null
@@ -0,0 +1,53 @@
+##
+#  CMake script for the step-60 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-60")
+
+# 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 2.8.12)
+
+FIND_PACKAGE(deal.II 9.0.0 QUIET
+  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_UMFPACK) # keep in one line
+  MESSAGE(FATAL_ERROR "
+Error! This tutorial requires a deal.II library that was configured with the following options:
+    DEAL_II_WITH_UMFPACK = ON
+However, the deal.II library found at ${DEAL_II_PATH} was configured with these options
+    DEAL_II_WITH_UMFPACK = ${DEAL_II_WITH_UMFPACK}
+which conflict with the requirements."
+    )
+ENDIF()
+
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox
new file mode 100644 (file)
index 0000000..7fd0508
--- /dev/null
@@ -0,0 +1,205 @@
+<br>
+
+<i>This program was contributed by Luca Heltai and Giovanni Alzetta, SISSA, Trieste.
+
+<a name="Introduction"></a><h1>Introduction</h1>
+
+
+<a name="DistributedLagrangemultipliers"></a><h3> Distributed Lagrange multipliers</h3>
+
+
+In this tutorial we consider the case of two domains, $\Omega$ in
+$R^{\text{spacedim}}$ and $\Gamma$ in $R^{\text{dim}}$, where $\Gamma$ is
+embedded in $\Omega$ and thus $\Gamma \subseteq \Omega$. We want to solve a
+partial differential equation on $\Omega$, enforcing some conditions on the
+solution of the problem *on the embedded domain* $\Gamma$.
+
+There are two interesting scenarios:
+
+- the geometrical dimension `dim` of the embedded domain $\Gamma$ is the
+same of the domain $\Omega$ (`spacedim`), that is, the spacedim-dimensional measure
+of $\Gamma$ is not zero, or
+
+- the embedded domain $\Gamma$ has an intrinsic dimension `dim` which is smaller
+than that of $\Omega$ (`spacedim`), thus its spacedim-dimensional measure is
+zero; for example it is a curve embedded in a two
+dimensional domain, or a surface embedded in a three-dimensional domain.
+
+In both cases define the restriction operator $\gamma$ as the operator
+that, given a continuous function on $\Omega$, returns its (continuous)
+restriction on $\Gamma$, i.e.,
+\f[
+\gamma : C^0(\Omega) \mapsto C^0(\Gamma), \quad \text{ s.t. } \gamma u = u|_{\Gamma} \in C^0(\Gamma),
+\quad \forall u \in C^0(\Omega).
+\f]
+
+It is well known that the operator $\gamma$ can be extended to a continuous
+operator on $H^1(\Omega)$, mapping functions in $H^1(\Omega)$ to functions in
+$H^1(\Gamma)$ when the intrinsic dimension of $\Gamma$ is the same of $\Omega$.
+
+The same is true, with a less regular range space (namely $H^{1/2}(\Gamma)$),
+when the dimension of $\Gamma$ is one less with respect to $\Omega$, and
+$\Gamma$ does not have a boundary. In this second case, the operator $\gamma$
+is also known as the *trace* operator, and it is well defined for Lipschitz
+co-dimension one curves and surfaces $\Gamma$ embedded in $\Omega$.
+
+The co-dimension two case is a little more complicated, and in general it is
+not possible to construct a continous trace operator, not even from
+$H^1(\Omega)$ to $L^2(\Gamma)$, when the dimension of $\Gamma$ is zero or one
+respectively in two and three dimensions.
+
+In this tutorial program we're not interested in further details on $\gamma$:
+we take the extension $\gamma$ for granted, assuming that the dimension of the
+embedded domain (`dim`) is always smaller by one or equal w.r.t. the dimension
+of the embedding domain $\Omega$ (`spacedim`).
+
+We are going to solve the following differential problem: given a sufficiently
+regular function $g$ on $\Gamma$, find the solution $u$ to
+
+@f{eqnarray*}
+- \Delta u + \gamma^T \lambda &=& 0 & \text{ in } \Omega\\
+\gamma u &=& g &  \text{ in } \Gamma \\
+u & = & 0 & \text{ on } \partial\Omega.
+@f}
+
+This is a constrained problem, where we are looking for a harmonic function
+$u$, that satisfies homogeneous boundary conditions on $\partial\Omega$, subject
+to the constraint $\gamma u = g$ using a Lagrange multiplier.
+
+The variational formulation can be derived by introducing two infinite
+dimensional spaces $V(\Omega)$ and $Q^*(\Gamma)$, respectively for the solution
+$u$ and for the Lagrange multiplier $\lambda$.
+
+Multiplying the first equation by $v \in V(\Omega)$ and the second by $q \in
+Q(\Gamma)$, integrating by parts when possible, and exploiting the boundary
+conditions on $\partial\Omega$, we obtain the following variational problem:
+
+Given a sufficiently regular function $g$ on $\Gamma$, find the solution $u$ to
+@f{eqnarray*}
+(\nabla u, \nabla v)_{\Omega} + (\lambda, \gamma v)_{\Gamma} &=& 0 & \forall v \in V(\Omega) \\
+(\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} &  \forall q \in Q(\Gamma),
+@f}
+
+where $(\cdot, \cdot)_{\Omega}$ (respectively $(\cdot, \cdot)_{\Gamma}$)
+represent the $L^2$ scalar product in $\Omega$ (respectively in $\Gamma$).
+
+Inspection of the variational formulation tells us that the space $V(\Omega)$
+can be taken to be $H^1_0(\Omega)$. The space $Q(\Gamma)$, in the co-dimension
+zero case, should be taken as $H^1(\Gamma)$, while in the co-dimension one case
+should be taken as $H^{1/2}(\Gamma)$.
+
+The function $g$, therefore, should be either in $H^1(\Gamma)$ (for the
+co-dimension zero case) or $H^{1/2}(\Gamma)$ (for the co-dimension one case).
+This leaves us with a Lagrange multiplier $\lambda$ in $Q^*(\Gamma)$, which is
+either $(H^1(\Gamma))^*$ or $H^{-1/2}(\Gamma)$.
+
+There are two options for the discretisation of the problem above. One could choose
+matching discretisations, where the Triangulation for $\Gamma$ is aligned with the
+Triangulation for $\Omega$, or one could choose to discretize the two domains in
+a completely independent way.
+
+While the first option is clearly more indicated for the simple problem we
+proposed above, if the domain $\Gamma$ was to be time dependent, then the
+second option could be a more viable solution.
+
+The technique we describe here is indicated in the literature with many names:
+the **immersed finite element method**, the **fictitious boundary method**, the
+**distributed Lagrange multiplier method**, and others. The main principle is
+that the discretisation of the two grids and of the two finite element spaces
+are kept completely independent. This technique is particularly efficient for
+the simulation of fluid-structure interaction problems, where the configuration
+of the embedded structure is part of the problem itself, and one solves a
+(possibly non-linear) elastic problem to determine the (time dependent)
+configuration of $\Gamma$, and a (possibly non-linear) flow problem in $\Omega
+\setminus \Gamma$, plus coupling conditions on the interface between the fluid
+and the solid.
+
+In this tutorial program we keep things a little simpler, and we assume that the
+configuration of the embedded domain is given in one of two possible ways:
+
+- as a deformation mapping $\psi: \Gamma_0 \mapsto \Gamma \subseteq \Omega$,
+defined on a continuous finite dimensional space on $\Gamma_0$ and representing,
+for any point $x \in \Gamma_0$ its coordinate $\psi(x)$ in $\Omega$;
+
+- as a displacement mapping $\delta \psi(x) = \psi(x)-x$ for $x\in \Gamma_0$,
+representing for any point $x$ the displacement vector to apply in order to
+deform $x$ to its actual configuration $\psi(x) = x +\delta\psi(x)$.
+
+We define the embedded reference domain $\Gamma_0$ `embedded_grid`, and on
+this domain, we construct a finite dimensional space (`embedded_configuration_dh`)
+to describe either the deformation or the displacement, through a FiniteElement
+system of FE_Q objects (`embedded_configuration_fe`). This finite dimensional
+space is used only to interpolate a user supplied function
+(`embedded_configuration_function`) representing either $\psi$ (if the
+parameter `use_displacement` is set to false) or $\delta\psi$ (if the parameter
+`use_displacement` is set to true).
+
+The Lagrange multiplier $\lambda$ and the user supplied function $g$ are
+defined through another finite dimensional space `embedded_dh`, and through
+another FiniteElement `embedded_fe`, using the same reference domain. In
+order to take into account the deformation of the domain, a MappingFEField or a
+MappingQEulerian object are initialized with the `embedded_configuration` vector.
+
+In the embedding space, a standard finite dimensional space `space_dh` is
+constructed on the embedding grid `space_grid` (`space_dh`), using the
+FiniteElement `space_fe`, following almost verbatim what was done in step-6.
+
+We represent the discretizations of the spaces $V$ and $Q$ with
+\f[
+V_h(\Omega) = \text{span} \{v_i\}_{i=1}^n
+\f]
+and
+\f[
+Q_h(\Gamma) = \text{span} \{q_i\}_{i=1}^m
+\f]
+respectively, where $n$ is the dimension of `space_dh`, and $m$
+the dimension of `embedded_dh`.
+
+Once all the finite dimensional spaces are defined, we are left with the following
+finite dimensional system of equations:
+
+\f[
+\begin{pmatrix}
+K & C^T \\
+C & 0
+\end{pmatrix}
+\begin{pmatrix}
+u \\
+\lambda
+\end{pmatrix}
+=
+\begin{pmatrix}
+0 \\
+G
+\end{pmatrix}
+\f]
+
+where
+
+@f{eqnarray*}
+K_{ij} &:=& (\nabla v_j, \nabla v_i)_\Omega   & i,j=1,\dots,n \\
+C_{\alpha j} &:=& (v_j, \nabla q_\alpha)_\Gamma  &j=1,\dots,n, \alpha = 1,\dots, m \\\\
+G_{\alpha} &:=& (g, q_\alpha)_\Gamma & \alpha = 1,\dots, m.
+@f}
+
+While the matrix $K$ is the standard stiffness matrix for the Poisson problem
+on $\Omega$, and the vector $G$ is a standard right-hand-side vector for a
+finite element problem with forcing term $g$ on $\Gamma$, (see, for example,
+step-3), the matrix $C$ or its transpose $C^T$ are non-standard since they
+couple information on two non-matching grids.
+
+In particular, the integral that appear in the computation of a single entry of $C$,
+is computed on $\Gamma$. As usual in finite elements, we split this integral on each
+cell of the triangulation used to discretize $\Gamma$, we tranform the integral on $K$ to
+an integral on the reference element $\hat K$, where $F_{K}$ is the corresponding
+shape function, and compute the integral there using a quadrature formula:
+
+\f[
+C_{\alpha j} := (v_j, \nabla q_\alpha)_\Gamma  = \sum_{K\in \Gamma} \int_{\hat K}
+\hat q_\alpha(\hat x) (v_j \circ F_{K}) (\hat x) J_K (\hat x) \mathrm{d} \hat x =
+\sum_{K\in \Gamma} \sum_{i=1}^{n_q}  \big(\hat q_\alpha(\hat x_i)  (v_j \circ F_{K}) (\hat x_i) J_K (\hat x_i) w_i \big)
+\f]
+
+We solve the above saddle point problem by iterating over the Schur complement
+(which is described, for example, in step-20),
+and we construct such complement using LinearOperator classes.
diff --git a/examples/step-60/doc/kind b/examples/step-60/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox
new file mode 100644 (file)
index 0000000..9c7a88b
--- /dev/null
@@ -0,0 +1,37 @@
+<h1>Results</h1>
+
+<h3> Test case 1: </h3>
+
+For the default problem the value of u on Gamma is 1 and on $\partial\Omega$
+is 0. This means we expect the following solution:
+
+
+
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+<h4> Different Parameters</h4>
+
+<h4> Parallel Code </h4>
+
+The simple code proposed here can serve as a starting point for more
+complex problems which, to be solved, need to be run on parallel
+code, possibly using distributed meshes (see step-17, step-40, and the
+documentation for parallel::shared::Triangulation and
+parallel::distributed::Triangulation ).
+
+When using non-matching grids in parallel a problem arises:
+to compute the matrix $C$ a process needs information about both
+meshes on the same portion of real space but, when working with
+distributed meshes, this might not be the case.
+Various strategies can be implemented to tackle this problem:
+
+- distribute the two meshes so that this constraint is satisfied
+- use communication for the parts of real space where the constraint is
+ not satisfied
+- make use of a shared triangulation and a distributed triangulation
+
+The latter strategy is clearly the easier to implement, as all
+the function used in this tutorial program can work lettin $\Omega$
+be distributed and $\Gamma$ be a shared triangulation.
diff --git a/examples/step-60/doc/tooltip b/examples/step-60/doc/tooltip
new file mode 100644 (file)
index 0000000..0bb91e2
--- /dev/null
@@ -0,0 +1 @@
+The fictitious domain method using distributed Lagrange multipliers
diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc
new file mode 100644 (file)
index 0000000..5483e95
--- /dev/null
@@ -0,0 +1,1056 @@
+/* ---------------------------------------------------------------------
+ *
+ * Copyright (C) 2018 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 at
+ * the top level of the deal.II distribution.
+ *
+ * ---------------------------------------------------------------------
+ *
+ * Author: Luca Heltai, Giovanni Alzetta,
+ * International School for Advanced Studies, Trieste, 2018
+ */
+
+// @sect3{Include files}
+// Most of these have been introduced elsewhere, we'll comment only on the new ones.
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/base/parameter_acceptor.h>
+
+// The parameter acceptor class is the first novelty of this tutorial program:
+// in general parameter files are used to steer the execution of a program
+// at run time. While even a simple approach saves compiling time, as the same
+// executable can be run with different parameter settings, it can become
+// difficult to handle hundreds of parameters simultaneously while maintaining
+// compatibility between different programs. This is where the class
+// ParameterAcceptor proves useful.
+//
+// This class is used to define a public interface for classes that want to use
+// a single global ParameterHandler to handle parameters. The class provides a
+// static ParameterHandler member, namely ParameterAcceptor::prm, and
+// implements the "Command design patter" (see, for example, E. Gamma, R. Helm,
+// R. Johnson, J. Vlissides, Design Patterns: Elements of Reusable
+// Object-Oriented Software, Addison-Wesley Professional, 1994.
+// https://goo.gl/FNYByc).
+//
+// ParameterAcceptor provides a global subscription mechanism. Whenever an
+// object of a class derived from ParameterAcceptor is constructed, a pointer
+// to that object-of-derived-type is registered, together with a section entry
+// in the parameter file. Such registry is traversed upon invocation of the
+// single function ParameterAcceptor::initialize(file.prm) which in turn makes
+// sure that all classes stored in the global register declare the parameters
+// they will be using, and after having declared them, it reads the content of
+// `file.prm` to parse the actual parameters.
+//
+// If you call the method ParameterHandler::add_parameter for each of the
+// parameters you want to use in your code, there is nothing else you need to
+// do. If you are using an already existing class that provides the two
+// functions `declare_parameters` and `parse_parameters`, you can still use
+// ParameterAcceptor, by encapsulating the existing class into a
+// ParameterAcceptorProxy class.
+//
+// In this example, we'll use both strategies, using ParameterAcceptorProxy for
+// deal.II classes, and deriving our own parameter classes directly from
+// ParameterAcceptor.
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/grid/grid_tools_cache.h>
+
+// The other new include file is the one that contains the GridTools::Cache
+// class. The structure of deal.II, as many modern numerical libraries, is
+// organized following a Directed Acyclic Graph (DAG). A DAG is a directed graph
+// with topological ordering: each node structurally represents an object, and
+// one or more directed edges represent how it can be used to generate new
+// objects. This has several advantages, but intrinsically creates “asymmetries”
+// as certain operations are fast unlike their inverse. For example, in deal.II
+// finding the vertices of a cell has very low computational cost, while finding
+// all the cells that share a vertex requires a non-trivial computation unless a
+// new data structure is added.
+//
+// Since inverse operations are usually not needed in a finite element code,
+// these are implemented in GridTools without the use of extra data structures
+// related to the Triangulation which would make them much faster. One such data
+// structure, for example, is a map from the vertices of a Triangulation to all
+// cells that share those vertices, which would reduce the computations needed
+// to answer to the previous example.
+//
+// Some methods, for example GridTools::find_active_cell_around_point, make
+// heavy usage of these non-standard operations. If you need to call these
+// methods more than once, it becomes convenient to store those data structures
+// somewhere. GridTools::Cache does exactly this, giving you access to
+// previously computed objects, or computing them on the fly (and then storing
+// them inside the class for later use), and making sure that whenever the
+// Triangulation is updated, also the relevant data strucutres are recomputed.
+
+#include <deal.II/fe/fe.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/fe/mapping_fe_field.h>
+
+// In this example, we will be using a reference domain to describe an embedded
+// Triangulation, deformed through a finite element vector field.
+//
+// The two include files above contain the definition of two classes that can be
+// used in these cases. MappingQEulerian allows one to describe a domain through
+// a *displacement* field, based on a FESystem[FE_Q(p)^spacedim] finite element
+// space. The second is a little more generic, and allows you to use arbitrary
+// vector FiniteElement spaces, as long as they provide a *continuous*
+// description of your domain. In this case, the description is done through the
+// actual *deformation* field, rather than a *displacement* field.
+//
+// Which one is used depends on how the user wants to specify the reference
+// domain, and/or the actual configuration. We'll provide both options, and
+// experiment a little in the results section of this tutorial program.
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/base/parsed_function.h>
+
+// The parsed function class is another new entry. It allows one to create a
+// Function object, starting from a string in a parameter file which is parsed
+// into an object that you can use anywhere deal.II accepts a Function (for
+// example, for interpolation, boundary conditions, etc.).
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#include <deal.II/non_matching/coupling.h>
+
+// This is the last new entry for this tutorial program. The namespace
+// NonMatching contains a few methods that are useful when performing
+// computations on non-matching grids, or on curves that are not aligned with
+// the underlying mesh.
+//
+// We'll discuss its use in details later on in the `setup_coupling` method.
+
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/sparse_direct.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/linear_operator_tools.h>
+
+#include <iostream>
+#include <fstream>
+
+namespace Step60
+{
+  using namespace dealii;
+
+  // In the DistributedLagrangeProblem, we need two parameters describing the
+  // dimensions of the domain $\Gamma$ (`dim`) and of the domain $\Omega$
+  // (`spacedim`).
+  //
+  // These will be used to initialize a Triangulation<dim,spacedim> (for
+  // $\Gamma$) and a Triangulation<spacedim,spacedim> (for $\Omega$).
+  //
+  // A novelty w.r.t. other tutorial programs is the heavy use of
+  // std::unique_ptr. These behave like classical pointers, with the advantage
+  // of doing automatic house-keeping: the contained object is automatically
+  // destroyed as soon as the unique_ptr goes out of scope, even if it is inside
+  // a container or there's an exception. Moreover it does not allow for
+  // duplicate pointers, which prevents ownership problems. We do this, because
+  // we want to be able to i) construct the problem, ii) read the parameters,
+  // and iii) initialize all objects according to what is specified in a
+  // parameter file.
+  //
+  // We construct the parameters of our problem in the internal class
+  // DistributedLagrangeProblemParameters, derived from ParameterAcceptor. The
+  // DistributedLagrangeProblem class takes a const reference to a
+  // DistributedLagrangeProblemParameters object, so that it is not possible to
+  // modify the parameters from within the DistributedLagrangeProblem class
+  // itself.
+  //
+  // We could have initialized the parameters first, and then pass the
+  // parameters to the DistributedLagrangeProblem assuming all entries are set to
+  // the desired values, but this has two disadvantages:
+  //
+  // - we should not make assumptions on how the user initializes a class that
+  // is not under our direct control. If the user fails to initialize the
+  // class, we should notice and throw an exception;
+  //
+  // - not all objects that need to read parameters from a parameter file may
+  // be available when we construct the DistributedLagrangeProblemParameters;
+  // this is often the case for complex programs, with multiple physics, or
+  // where we reuse existing code in some external classes. We simulate this by
+  // keeping some "complex" objects, like ParsedFunction objects, inside the
+  // DistributedLagrangeProblem instead of inside the
+  // DistributedLagrangeProblemParameters.
+  //
+  // Here we assume that upon construction, the classes that build up our
+  // problem are not usable yet. Parsing the parameter file is what ensures we have
+  // all ingredients to build up our classes, and we design them so that if parsing
+  // fails, or is not executed, the run is aborted.
+
+  template <int dim, int spacedim=dim>
+  class DistributedLagrangeProblem
+  {
+  public:
+
+    // The DistributedLagrangeProblemParameters is derived from
+    // ParameterAcceptor. This allows us to use the
+    // ParameterAcceptor::add_parameter methods in its constructor.
+    //
+    // The members of this function are all non-const, but the
+    // DistributedLagrangeProblem class takes a const reference to a
+    // DistributedLagrangeProblemParameters object: this ensures that
+    // parameters are not modified from within the DistributedLagrangeProblem
+    // class.
+
+    class DistributedLagrangeProblemParameters : public ParameterAcceptor
+    {
+    public:
+      DistributedLagrangeProblemParameters();
+
+      // The parameters now described can all be set externally
+      // using a parameter file: if no parameter file is present when running
+      // the executable, the program shall create a parameters.prm file with
+      // the default parameters here encoded, and then abort as no
+      // parameters.prm file was found.
+
+      // Initial refinement for the embedding grid, corresponding to the domain
+      // $\Omega$.
+      unsigned int initial_refinement                           = 4;
+
+      // The interaction between the embedded grid $\Omega$ and the embedding
+      // grid $\Gamma$ is handled throught the computation of $C$, which
+      // involves all cells of $\Omega$ overlapping with parts of $\Gamma$:
+      // a higher refinement of such cells might improve the results quality.
+      // For this reason we define `delta_refinement`: if it is greater
+      // than zero, then we mark each cell of the space grid that contains
+      // a vertex of the embedded grid, execute the refinement, and repeat
+      // this process `delta_refinement` times.
+      unsigned int delta_refinement                             = 3;
+
+      // Starting refinement of the embedded grid, corresponding to the domain
+      // $\Gamma$.
+      unsigned int initial_embedded_refinement                  = 7;
+
+      // A list of boundary ids where we impose homogeneous Dirichlet boundary
+      // conditions. On the remaining boundary ids (if any), we impose
+      // homogeneous Neumann boundary conditions
+      std::list<types::boundary_id> homogeneous_dirichlet_ids   {0};
+
+      // FiniteElement degree of the embedding space: $V_h(\Omega)$
+      unsigned int embedding_space_finite_element_degree        = 1;
+
+      // FiniteElement degree of the embedded space: $Q_h(\Gamma)$
+      unsigned int embedded_space_finite_element_degree         = 1;
+
+      // FiniteElement degree of the space used to describe the deformation
+      // of the embedded domain
+      unsigned int embedded_configuration_finite_element_degree = 1;
+
+      // Order of the quadrature formula used to integrate the coupling
+      unsigned int coupling_quadrature_order                    = 3;
+
+      // If set to true, then the embedded configuration function is
+      // interpreted as a displacement function
+      bool use_displacement                                     = false;
+
+      // Level of verbosity to use in the output
+      unsigned int verbosity_level                              = 10;
+
+      // A flag to keep track if we were initialized or not
+      bool initialized                                          = false;
+    };
+
+    DistributedLagrangeProblem(const DistributedLagrangeProblemParameters &parameters);
+
+    // Entry point for the DistributedLagrangeProblem
+    void run();
+
+  private:
+    // Object containing the actual parameters
+    const DistributedLagrangeProblemParameters &parameters;
+
+    // The following functions are similar to all other tutorial programs, with
+    // the exception that we now need to setup things for two different
+    // families of objects, namely the ones related to the *embedding* grids, and the
+    // ones related to the *embedded* one.
+
+    void setup_grids_and_dofs();
+
+    void setup_embedding_dofs();
+
+    void setup_embedded_dofs();
+
+    // The only unconventional function we have here is the `setup_coupling()` method,
+    // used to generate the sparsity patter for the coupling matrix $C$.
+
+    void setup_coupling();
+
+    void assemble_system();
+
+    void solve();
+
+    void output_results();
+
+
+    // First we gather all the objects related to the embedding space geometry
+
+    std::unique_ptr<Triangulation<spacedim> > space_grid;
+    std::unique_ptr<GridTools::Cache<spacedim, spacedim> > space_grid_tools_cache;
+    std::unique_ptr<FiniteElement<spacedim> > space_fe;
+    std::unique_ptr<DoFHandler<spacedim> > space_dh;
+
+    // Then the ones related to the embedded grid, with the DoFHandler associated
+    // to the Lagrange multiplier `lambda`
+
+    std::unique_ptr<Triangulation<dim, spacedim> > embedded_grid;
+    std::unique_ptr<FiniteElement<dim, spacedim> > embedded_fe;
+    std::unique_ptr<DoFHandler<dim, spacedim> > embedded_dh;
+
+    // And finally, everything that is needed to *deform* the embedded
+    // triangulation
+    std::unique_ptr<FiniteElement<dim, spacedim> > embedded_configuration_fe;
+    std::unique_ptr<DoFHandler<dim, spacedim> > embedded_configuration_dh;
+    Vector<double> embedded_configuration;
+
+    // The ParameterAcceptorProxy class, is a "transparent" wrapper, derived
+    // from both ParameterAcceptor and the type passed as its template
+    // parameter. At construction, the arguments are split into two parts: the
+    // first argument is an std::string, forwarded to the ParameterAcceptor
+    // class, and containing the name of the section that should be used for
+    // this class, while all the remaining arguments are forwarded to the
+    // constructor of the templated type, in this case, to the
+    // Functions::ParsedFunction constructor.
+    //
+    // This class allows you to use existing classes in conjunction with the
+    // ParameterAcceptor registration mechanism, provided that those classes
+    // have the members `declare_parameters()` and `parse_parameters()`.
+    //
+    // This is the case here, making it fairly easy to exploit the
+    // Functions::ParsedFunction class: instead of requiring the user to create
+    // new Function objects in its code for the RHS, boundary functions, etc.,
+    // (like it is done in most of the other tutorials), here we allow the user
+    // to use deal.II interface to muParser (http://muparser.beltoforion.de),
+    // where the specification of the function is not done at compile time, but
+    // at run time, using a string that is parsed into an actual Function object.
+    //
+    // In this case, the `embedded_configuration_function` is a vector valued
+    // Function that can be interpreted as either a *deformation* or a
+    // *displacement* according to the boolean value of
+    // parameters.use_displacement. The number of components is specified later
+    // on in the construction.
+
+    ParameterAcceptorProxy<Functions::ParsedFunction<spacedim> >
+    embedded_configuration_function;
+
+
+    // The embedded mapping. Notice that the order in which we construct these
+    // unique pointers is important. They will be destroied in the reversed
+    // order, so it is important that we respect the dependency tree. In
+    // particular, the embedded mapping will depend on both the `embedded_dh`
+    // and the `embedded_configuration`. If we declare it after the above two,
+    // we are fine, otherwise we would have do release this pointer manually in the
+    // destructor, or we'd get an error like
+    //
+    // @code
+    // --------------------------------------------------------
+    // An error occurred in line <104> of file <../source/base/subscriptor.cc> in function
+    //     void dealii::Subscriptor::check_no_subscribers() const
+    // The violated condition was:
+    //     counter == 0
+    // Additional information:
+    //     (none)
+    // @endcode
+    //
+    // at the end of the program.
+
+    std::unique_ptr<Mapping<dim,spacedim> > embedded_mapping;
+
+    // We do the same thing to specify the value of the function $g$,
+    // which is what we want our solution to be in the embedded space.
+    // In this case the Function is a scalar one.
+
+    ParameterAcceptorProxy<Functions::ParsedFunction<spacedim> >
+    embedded_value_function;
+
+    // Similarly to what we have done with the Functions::ParsedFunction class,
+    // we repeat the same for the ReductionControl class, allowing us to
+    // specify all possible stopping criterions for the Schur complement
+    // iterative solver we'll use later on.
+
+    ParameterAcceptorProxy<ReductionControl> schur_solver_control;
+
+    // Next we gather all SparsityPattern, SparseMatrix, and Vector objects
+    // we'll need
+
+    SparsityPattern           stiffness_sparsity;
+    SparsityPattern           coupling_sparsity;
+    SparsityPattern           embedded_sparsity;
+
+    SparseMatrix<double>      stiffness_matrix;
+    SparseMatrix<double>      coupling_matrix;
+    SparseMatrix<double>      embedded_stiffness_matrix;
+
+    ConstraintMatrix          constraints;
+
+    Vector<double>            solution;
+    Vector<double>            rhs;
+
+    Vector<double>            lambda;
+    Vector<double>            embedded_rhs;
+    Vector<double>            embedded_value;
+
+    // The TimerOuput class is used to provide some statistics on
+    // the performance of our program.
+    TimerOutput monitor;
+  };
+
+  // At construction time, we initialize also the ParameterAcceptor class, with
+  // the section name we want our problem to use when parsing the parameter
+  // file.
+  //
+  // Parameter files can be organized into section/subsection/etc. :
+  // this has the advantage that defined objects share parameters when
+  // sharing the same section/subsection/etc. ParameterAcceptor allows
+  // to specify the section name using unix conventions on paths.
+  // If the section name starts with a slash ("/"), then the section is
+  // interpreted as an *absolute path*, ParameterAcceptor enters a subsection
+  // for each directory in the path, using the last name it encountered as
+  // the landing subsection for the current class.
+  //
+  // For example, if you construct your class using
+  // ParameterAcceptor("/first/second/third/My Class"), the parameters will be
+  // organized as follows:
+  //
+  // @code
+  // subsection first
+  //   subsection second
+  //     subsection third
+  //       subsection My Class
+  //        ... # all the parameters
+  //       end
+  //     end
+  //   end
+  // end
+  // @endcode
+  //
+  // Internally, the *current path* stored in ParameterAcceptor, is now
+  // considered to be "/first/second/third/", i.e. when you specify an
+  // absolute path, ParameterAcceptor *changes* the current section to the
+  // current path, i.e. to the path of the section name until the *last* "/".
+  //
+  // You can now construct another class derived from ParameterAcceptor using a
+  // relative path (e.g., ParameterAcceptor("My Other Class")) instead of the
+  // absolute one (e.g. ParameterAcceptor("/first/second/third/My Other Class")),
+  // obtaining:
+  // @code
+  // subsection first
+  //   subsection second
+  //     subsection third
+  //       subsection MyClass
+  //        ... # all the parameters
+  //       end
+  //       subsection My Other Class
+  //        ... # all the parameters of MyOtherClass
+  //       end
+  //     end
+  //   end
+  // end
+  // @endcode
+  //
+  // If the section name *ends* with a slash then subsequent classes will
+  // interpret this as a full path: for example, similar to the one above, if
+  // we have two classes, one initialized with
+  // ParameterAcceptor("/first/second/third/My Class/")
+  // and the other with ParameterAcceptor("My Other Class"), then the
+  // resulting parameter file will look like:
+  //
+  // @code
+  // subsection first
+  //   subsection second
+  //     subsection third
+  //       subsection MyClass
+  //        ... # all the parameters
+  //        ... # notice MyClass subsection does not end here
+  //        subsection My Other Class
+  //         ... # all the parameters of MyOtherClass
+  //        end
+  //       end # of subsection MyClass
+  //     end
+  //   end
+  // end
+  // @endcode
+  //
+  // We are going to exploit this, by making our
+  // DistributedLagrangeProblemParameters
+  // the *parent* of all subsequently constructed classes. Since most of the other
+  // classes are members of DistributedLagrangeProblem this allows, for example,
+  // to construct two DistributedLagrangeProblem for two different dimensions, without
+  // having conflicts in the parameters for the two problems.
+
+  template<int dim, int spacedim>
+  DistributedLagrangeProblem<dim,spacedim>::DistributedLagrangeProblemParameters::
+  DistributedLagrangeProblemParameters() :
+    ParameterAcceptor("/Distributed Lagrange<" + Utilities::int_to_string(dim)
+                      + "," + Utilities::int_to_string(spacedim) +">/")
+  {
+
+    // The ParameterAcceptor::add_parameter does a few things:
+    //
+    // - enters the subsection specified at construction time to ParameterAcceptor
+    //
+    // - calls the ParameterAcceptor::prm.add_parameter
+    //
+    // - calls any signal you may have attached to
+    // ParameterAcceptor::declare_parameters_call_back
+    //
+    // - leaves the subsection
+    //
+    // In turns, ParameterAcceptor::prm.add_parameter
+    //
+    // - declares an entry in the parameter handler for the given variable;
+    //
+    // - reads the value of the variable,
+    //
+    // - transforms it to a string, used as the default value for the parameter
+    // file
+    //
+    // - attaches an *action* to ParameterAcceptor::prm that monitors when a file
+    // is parsed, or when an entry is set, and when this happens, it updates the
+    // content of the given variable to the value parsed by the string
+
+    add_parameter("Initial embedding space refinement",
+                  initial_refinement);
+
+    add_parameter("Initial embedded space refinement",
+                  initial_embedded_refinement);
+
+    add_parameter("Local refinements steps near embedded domain",
+                  delta_refinement);
+
+    add_parameter("Homogeneous Dirichlet boundary ids",
+                  homogeneous_dirichlet_ids);
+
+    add_parameter("Use displacement in embedded interface",
+                  use_displacement);
+
+    add_parameter("Embedding space finite element degree",
+                  embedding_space_finite_element_degree);
+
+    add_parameter("Embedded space finite element degree",
+                  embedded_space_finite_element_degree);
+
+    add_parameter("Embedded configuration finite element degree",
+                  embedded_configuration_finite_element_degree);
+
+    add_parameter("Coupling quadrature order",
+                  coupling_quadrature_order);
+
+    add_parameter("Verbosity level",
+                  verbosity_level);
+
+    // If we get parsed, then the parameters are good to go. Set the internal
+    // variable `initialized` to true.
+    parse_parameters_call_back.connect(
+      [&]() -> void {initialized = true;}
+    );
+  }
+
+  // The constructor is pretty standard, with the exception of the `ParameterAcceptorProxy`
+  // objects, as explained earlier on.
+
+  template<int dim, int spacedim>
+  DistributedLagrangeProblem<dim,spacedim>::DistributedLagrangeProblem(
+    const DistributedLagrangeProblemParameters &parameters) :
+    parameters(parameters),
+    embedded_configuration_function("Embedded configuration", spacedim),
+    embedded_value_function("Embedded value"),
+    schur_solver_control("Schur solver control"),
+    monitor(std::cout,
+            TimerOutput::summary,
+            TimerOutput::cpu_and_wall_times)
+  {
+    // Here is a way to set default values for a ParameterAcceptor class
+    // that was constructed using ParameterAcceptorProxy.
+    //
+    // In this case, we set the default deformation of the embedded grid to be
+    // a circle with radius `R` and center (Cx, Cy), we set the default value
+    // for the embedded_value_function to be the constant one, and specify some
+    // sensible values for the SolverControl object.
+    //
+    // It is fundamental for $\Gamma% to be embedded: from the definition of
+    // $C_{\alpha j}$ is clear that, if $\Gamma \not\subseteq \Omega$, certain
+    // rows of the matrix $C$ shall be zero. This would be a problem, as the Schur
+    // complement method requires $C$ needs to have full column rank.
+
+    embedded_configuration_function.declare_parameters_call_back.connect(
+      [] () -> void
+    {
+      ParameterAcceptor::prm.set("Function constants", "R=.3, Cx=.4, Cy=.4");
+      ParameterAcceptor::prm.set("Function expression", "R*cos(2*pi*x)+Cx; R*sin(2*pi*x)+Cy");
+    });
+
+    embedded_value_function.declare_parameters_call_back.connect(
+      [] () -> void
+    {
+      ParameterAcceptor::prm.set("Function expression", "1");
+    });
+
+    schur_solver_control.declare_parameters_call_back.connect(
+      [] () -> void
+    {
+      ParameterAcceptor::prm.set("Max steps", "1000");
+      ParameterAcceptor::prm.set("Reduction", "1.e-12");
+      ParameterAcceptor::prm.set("Tolerance", "1.e-12");
+    });
+
+  }
+
+  // Nothing special here, except that we check if parsing was done before
+  // we actually attempt to run our program.
+
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim, spacedim>::run()
+  {
+    AssertThrow(parameters.initialized,
+                ExcNotInitialized());
+    deallog.depth_console(parameters.verbosity_level);
+
+    setup_grids_and_dofs();
+    setup_coupling();
+    assemble_system();
+    solve();
+    output_results();
+  }
+
+
+  // The function DistributedLagrangeProblem::setup_grids_and_dofs
+  // is used to set up the finite element spaces. Notice how
+  // @code std_cxx14::make_unique @endcode is used to create objects
+  // wrapped inside unique pointers
+
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::setup_grids_and_dofs()
+  {
+    TimerOutput::Scope timer_section(monitor, "Setup grids and dofs");
+
+    // Initializing $\Omega$:
+    // constructing the Triangulation and wrapping it into a unique_ptr
+    space_grid = std_cxx14::make_unique<Triangulation<spacedim> >();
+    GridGenerator::hyper_cube(*space_grid);
+    // Requesting the varius values to the parameters object, which is
+    // of type DistributedLagrangeProblemParameters
+    space_grid->refine_global(parameters.initial_refinement);
+    space_grid_tools_cache =
+      std_cxx14::make_unique<GridTools::Cache<spacedim, spacedim> >(*space_grid);
+    // The refinement of $\Omega$ depends on $\Gamma$: this means we need to
+    // set up $\Gamma$ before we can finish with $\Omega$
+    embedded_grid = std_cxx14::make_unique<Triangulation<dim,spacedim> >();
+    GridGenerator::hyper_cube(*embedded_grid);
+    embedded_grid->refine_global(parameters.initial_embedded_refinement);
+
+    // At this point we need to configure the deformation of the embedded grid
+    embedded_configuration_fe =
+      std_cxx14::make_unique<FESystem<dim,spacedim> >
+      (FE_Q<dim,spacedim>(parameters.embedded_configuration_finite_element_degree), spacedim);
+
+    embedded_configuration_dh =
+      std_cxx14::make_unique<DoFHandler<dim, spacedim> >(*embedded_grid);
+
+    embedded_configuration_dh->distribute_dofs(*embedded_configuration_fe);
+    embedded_configuration.reinit(embedded_configuration_dh->n_dofs());
+
+    // Interpolating the embedded configuration function
+    VectorTools::interpolate(*embedded_configuration_dh,
+                             embedded_configuration_function,
+                             embedded_configuration);
+
+    // Once we have the embedded configuration vector stored, we can interpret
+    // it according to what the user has specified in the parameter file, as a
+    // displacement, in which case we construct a mapping that *displaces* the
+    // position of each support point of our configuration finite element space
+    // by the specified amount on the corresponding configuration vector, or as
+    // an absolution position.
+    //
+    // In the first case, the class MappingQEulerian offers its services, while
+    // in the second one, we'll use the class MappingFEField. They are in fact
+    // very similar. MappingQEulerian will only work for systems of FE_Q finite
+    // element spaces, where the displacment vector is stored in the first
+    // `spacedim` components of the FESystem, and the degree given as a
+    // parameter at construction time, must match the degree of the first
+    // `spacedim` components.
+    //
+    // The class MappingFEField is slightly more general, in that it allows you
+    // to select arbitrary FiniteElement types when constructing your
+    // approximation. Naturally some choices may (or may not) make sense,
+    // according to the type of FiniteElement you choose. MappingFEField
+    // implements the pure iso-parametric concept, and can be used, for example,
+    // to implement iso-geometric analysis codes in deal.II, by combining it
+    // with the FEBernstein finite element class. In this example, we'll use the
+    // two interchangeably, by taking into account the fact that one
+    // configuration will be a `displacment`, while the other will be an
+    // absolute `deformation` field.
+
+    if (parameters.use_displacement == true)
+      embedded_mapping =
+        std_cxx14::make_unique<MappingQEulerian<dim, Vector<double>, spacedim> >
+        (parameters.embedded_configuration_finite_element_degree,
+         *embedded_configuration_dh, embedded_configuration);
+    else
+      embedded_mapping =
+        std_cxx14::make_unique<MappingFEField<dim, spacedim, Vector<double>, DoFHandler<dim,spacedim> > >
+        (*embedded_configuration_dh,
+         embedded_configuration);
+
+    // Estimating the diameter of the largest active cell of $\Gamma$, and
+    // the smalles one of $\Omega$
+    double embedded_space_maximal_diameter =
+      GridTools::maximal_cell_diameter(*embedded_grid, *embedded_mapping);
+    double embedding_space_minimal_diameter =
+      GridTools::minimal_cell_diameter(*space_grid);
+    // Setting up the $\Gamma$'s DoFs
+    setup_embedded_dofs();
+
+    // With the mapping in place, it is now possible to query what is the
+    // location of all support points associated with the `embedded_dh`, by
+    // calling the method DoFTools::map_dofs_to_support_points.
+    //
+    // This method has two variants. One that does *not* take a Mapping, and
+    // one that takes a Mapping. If you use the second type, like we are doing
+    // in this case, the support points are computed through the specified mapping,
+    // which can manipulate them accordingly.
+    //
+    // This is precisely what the `embedded_mapping` is there for.
+
+    std::vector<Point<spacedim> > support_points(embedded_dh->n_dofs());
+    if (parameters.delta_refinement != 0)
+      DoFTools::map_dofs_to_support_points(*embedded_mapping,
+                                           *embedded_dh,
+                                           support_points);
+
+    // Once we have the support points of the embedded finite element space, we
+    // would like to identify what cells of the embedding space contain what
+    // support point, to get a chance at refining the embedding grid where it is
+    // necessary, i.e., where the embedded grid is. This can be done manually,
+    // by looping over each point, and then calling the method
+    // Mapping::tranform_real_to_unit_cell for each cell of the embedding
+    // space, until we find one that returns points in the unit reference cell,
+    // or it can be done in a more intelligent way.
+    //
+    // The GridTools::find_active_cell_around_point is a possible option, that
+    // performs the above task in a cheaper way, by first identifying the
+    // closest vertex of the embedding Triangulation to the target point, and
+    // then by calling Mapping::tranform_real_to_unit_cell only for those cells
+    // that share the found vertex.
+    //
+    // In fact, there are algorithm in the GridTools namespace that exploit a
+    // GridTools::Cache object, and possibly a KDTree object to speed up these
+    // operations as much as possible.
+    //
+    // The simplest way to exploit the maximum speed, is by calling a
+    // specialized method, GridTools::compute_point_locations, that will store a
+    // lot of useful information and data structures during the first point
+    // search, and then reuse all of this for subsequent points.
+    //
+    // GridTools::compute_point_locations returns a tuple where the first
+    // element is a vector of cells containing the input points, in this
+    // case support_points. For refinement, this is the only information we need,
+    // and this is exactly what happens now.
+    //
+    // When we need to assemble a coupling matrix, however, we'll also need the
+    // reference location of each point to evaluate the basis functions of the
+    // embedding space. The other elements of the tuple returned by
+    // GridTools::compute_point_locations allows you to reconstruct, for each
+    // point, what cell contains it, and what is the location in the reference
+    // cell of the given point. Since this information is better grouped into
+    // cells, then this is what the algorithm returns: a tuple, containing a
+    // vector of all cells that have at least one point in them, together with a
+    // list of all reference points and their corresponding index in the
+    // original vector.
+    //
+    // In the following loop, we will be ignoring all returned objects except
+    // the first, identifying all cells contain at least one support point of
+    // the embedded space.
+    //
+    // Notice that we need to do some sanity checks, in the sense that we want
+    // to have an embedding grid which is well refined around the embedded grid,
+    // but where two consecutive support points lie either in the same cell, or
+    // in neighbor embedding cells.
+    //
+    // This is only possible if we ensure that the smallest cell size of the
+    // embedding grid is nonetheless bigger than the largest cell size of the
+    // embedded grid. Since users can modify both levels of refinements, as well
+    // as the amount of local refinement they want around the embedded grid, we
+    // make sure that the resulting meshes satisfy our requirements, and if this
+    // is not the case, we bail out with an exception.
+
+    for (unsigned int i=0; i<parameters.delta_refinement; ++i)
+      {
+        const auto point_locations = GridTools::compute_point_locations(*space_grid_tools_cache,
+                                     support_points);
+        const auto &cells = std::get<0>(point_locations);
+        for (auto cell : cells)
+          cell->set_refine_flag();
+        space_grid->execute_coarsening_and_refinement();
+        embedding_space_minimal_diameter = GridTools::minimal_cell_diameter(*space_grid);
+        AssertThrow(embedded_space_maximal_diameter < embedding_space_minimal_diameter,
+                    ExcMessage("The embedding grid is too refined (or the embedded grid"
+                               "is too coarse). Adjust the parameters so that the minimal"
+                               "grid size of the embedding grid is larger "
+                               "than the maximal grid size of the embedded grid."))
+      }
+
+    deallog << "Embedding minimal diameter: " << embedding_space_minimal_diameter
+            << ", embedded maximal diameter: " << embedded_space_maximal_diameter
+            << ", ratio: "
+            << embedded_space_maximal_diameter/embedding_space_minimal_diameter << std::endl;
+
+    // $\Omega$ has been refined and we can now set up its DoFs
+    setup_embedding_dofs();
+  }
+
+  // We now set up the DoFs of $\Omega$ and $\Gamma$: since they are fundamentally
+  // independent (except for the fact that $\Omega$'s mesh is more refined "around"
+  // $\Gamma$) the procedure is standard
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::setup_embedding_dofs()
+  {
+    space_dh = std_cxx14::make_unique<DoFHandler<spacedim> >(*space_grid);
+    space_fe = std_cxx14::make_unique<FE_Q<spacedim> >
+               (parameters.embedding_space_finite_element_degree);
+    space_dh->distribute_dofs(*space_fe);
+
+    DoFTools::make_hanging_node_constraints(*space_dh, constraints);
+    for (auto id:parameters.homogeneous_dirichlet_ids)
+      {
+        VectorTools::interpolate_boundary_values(*space_dh, id,
+                                                 Functions::ZeroFunction<spacedim>(),
+                                                 constraints);
+      }
+    constraints.close();
+
+    // By definition the stiffness matrix involves only $\Omega$'s DoFs
+    DynamicSparsityPattern dsp(space_dh->n_dofs(), space_dh->n_dofs());
+    DoFTools::make_sparsity_pattern(*space_dh, dsp, constraints);
+    stiffness_sparsity.copy_from(dsp);
+    stiffness_matrix.reinit(stiffness_sparsity);
+    solution.reinit(space_dh->n_dofs());
+    rhs.reinit(space_dh->n_dofs());
+
+    deallog << "Embedding dofs: " << space_dh->n_dofs() << std::endl;
+  }
+
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::setup_embedded_dofs()
+  {
+    embedded_dh = std_cxx14::make_unique<DoFHandler<dim,spacedim> >(*embedded_grid);
+    embedded_fe = std_cxx14::make_unique<FE_Q<dim,spacedim> >
+                  (parameters.embedded_space_finite_element_degree);
+    embedded_dh->distribute_dofs(*embedded_fe);
+
+    DynamicSparsityPattern dsp(embedded_dh->n_dofs(), embedded_dh->n_dofs());
+    DoFTools::make_sparsity_pattern(*embedded_dh, dsp);
+    embedded_sparsity.copy_from(dsp);
+    embedded_stiffness_matrix.reinit(embedded_sparsity);
+    // By definition the rhs of the system we're solving involves only a zero
+    // vector and $G$, which is computed using only $\Gamma$'s DoFs
+    lambda.reinit(embedded_dh->n_dofs());
+    embedded_rhs.reinit(embedded_dh->n_dofs());
+    embedded_value.reinit(embedded_dh->n_dofs());
+
+    deallog << "Embedded dofs: " << embedded_dh->n_dofs() << std::endl;
+  }
+
+  // Creating the coupling sparsity pattern is a complex operation,
+  // but it can be easily done using the
+  // NonMatching::create_coupling_sparsity_pattern, which requires the
+  // two DoFHandlers, the quadrature points for the coupling,
+  // a DynamicSparsityPattern (which then needs to be copied into the
+  // sparsity one, as usual), the component mask for the embedding and
+  // embedded Triangulation (which we leave empty) and the mappings
+  // for both the embedding and the embedded Triangulation.
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::setup_coupling()
+  {
+    TimerOutput::Scope timer_section(monitor, "Setup coupling");
+
+    QGauss<dim> quad(parameters.coupling_quadrature_order);
+
+    DynamicSparsityPattern dsp(space_dh->n_dofs(), embedded_dh->n_dofs());
+
+    NonMatching::create_coupling_sparsity_pattern(*space_dh,
+                                                  *embedded_dh,
+                                                  quad,
+                                                  dsp, ComponentMask(), ComponentMask(),
+                                                  StaticMappingQ1<spacedim>::mapping,
+                                                  *embedded_mapping);
+    coupling_sparsity.copy_from(dsp);
+    coupling_matrix.reinit(coupling_sparsity);
+  }
+
+  // This function creates the matrices: as noted before computing
+  // the stiffness matrix and the rhs is a standard procedure
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::assemble_system()
+  {
+    {
+      TimerOutput::Scope timer_section(monitor, "Assemble system");
+
+      // Embedding stiffness matrix $K$
+      MatrixTools::create_laplace_matrix(*space_dh, QGauss<spacedim>(2*space_fe->degree+1),
+                                         stiffness_matrix, (const Function<spacedim> *) nullptr, constraints);
+
+      // Embedded stiffness matrix and rhs vector $G$
+      MatrixTools::create_laplace_matrix(*embedded_mapping,
+                                         *embedded_dh,
+                                         QGauss<dim>(2*embedded_fe->degree+1),
+                                         embedded_stiffness_matrix,
+                                         embedded_value_function,
+                                         embedded_rhs);
+    }
+    {
+      TimerOutput::Scope timer_section(monitor, "Assemble coupling system");
+
+      // To compute the coupling matrix we use the NonMatching::create_coupling_mass_matrix
+      // tool, which works similarly to NonMatching::create_coupling_sparsity_pattern,
+      // requiring only an additional parameter: a constraint matrix
+      QGauss<dim> quad(parameters.coupling_quadrature_order);
+      NonMatching::create_coupling_mass_matrix(*space_dh,
+                                               *embedded_dh,
+                                               quad,
+                                               coupling_matrix, ConstraintMatrix(),
+                                               ComponentMask(), ComponentMask(),
+                                               StaticMappingQ1<spacedim>::mapping,
+                                               *embedded_mapping);
+
+      VectorTools::interpolate(*embedded_dh, embedded_value_function, embedded_value);
+    }
+  }
+
+  // All parts have been assembled: we solve the system
+  // using the Schur complement method
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::solve()
+  {
+    TimerOutput::Scope timer_section(monitor, "Solve system");
+
+    // Start by creating the inverse stiffness matrix
+    SparseDirectUMFPACK K_inv_umfpack;
+    K_inv_umfpack.initialize(stiffness_matrix);
+
+    // Same thing, for the embedded space
+    SparseDirectUMFPACK A_inv_umfpack;
+    A_inv_umfpack.initialize(embedded_stiffness_matrix);
+    // Initializing the operators, as described in the introduction
+    auto K = linear_operator(stiffness_matrix);
+    auto A = linear_operator(embedded_stiffness_matrix);
+    auto Ct = linear_operator(coupling_matrix);
+    auto C = transpose_operator(Ct);
+
+    auto K_inv = linear_operator(K, K_inv_umfpack);
+    auto A_inv = linear_operator(A, A_inv_umfpack);
+
+    auto S = C*K_inv*Ct;
+    // Using the Schur complement method
+    SolverCG<Vector<double> > solver_cg(schur_solver_control);
+    auto S_inv = inverse_operator(S, solver_cg, A_inv);
+
+    lambda = S_inv * embedded_rhs;
+
+    solution = K_inv * Ct * lambda;
+
+    constraints.distribute(solution);
+  }
+
+  // Standard result output on two separate files, one
+  // for each mesh
+  template<int dim, int spacedim>
+  void DistributedLagrangeProblem<dim,spacedim>::output_results()
+  {
+    TimerOutput::Scope timer_section(monitor, "Output results");
+
+    DataOut<spacedim> embedding_out;
+
+    std::ofstream embedding_out_file("embedding.vtu");
+
+    embedding_out.attach_dof_handler(*space_dh);
+    embedding_out.add_data_vector(solution, "solution");
+    embedding_out.build_patches(parameters.embedding_space_finite_element_degree);
+    embedding_out.write_vtu(embedding_out_file);
+
+    // The only difference between the two output routines, is that in the
+    // second case, we want to output the data on the current configuration, and
+    // not on the reference one. This is possible by passing the actual
+    // embedded_mapping to the DataOut::build_patches function. The mapping will
+    // take care of outputting the result on the actual deformed configuration.
+
+    DataOut<dim, DoFHandler<dim,spacedim> > embedded_out;
+
+    std::ofstream embedded_out_file("embedded.vtu");
+
+    embedded_out.attach_dof_handler(*embedded_dh);
+    embedded_out.add_data_vector(lambda, "lambda");
+    embedded_out.add_data_vector(embedded_value, "g");
+    embedded_out.build_patches(*embedded_mapping,
+                               parameters.embedded_space_finite_element_degree);
+    embedded_out.write_vtu(embedded_out_file);
+  }
+}
+
+
+
+int main()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace Step60;
+
+      const unsigned int dim=1, spacedim=2;
+
+      // Differently to what happens in other tutorial programs, here we the the
+      // ParameterAcceptor style of initialization, i.e., all objects are first
+      // constructed, and then a single call to the static method
+      // ParameterAcceptor::initialize is issued to fill all parameters of the
+      // classes that are derived from ParameterAcceptor.
+
+      DistributedLagrangeProblem<dim, spacedim>::DistributedLagrangeProblemParameters parameters;
+      DistributedLagrangeProblem<dim, spacedim> problem(parameters);
+      ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm");
+      problem.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  return 0;
+}

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.