--- /dev/null
+<br>
+
+<i>This program was contributed by Luca Heltai (SISSA, Trieste) and Bruno Blais (Polytechnique Montréal)
+</i>
+
+<b>Change this!!!!</b>
+
+@dealiiTutorialDOI{10.5281/zenodo.1243280,https://zenodo.org/badge/DOI/10.5281/zenodo.1243280.svg}
+
+
+<h1>Introduction</h1>
+
+<h3>Smart ass title here.</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$ ($\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$ (read <a
+href="https://en.wikipedia.org/wiki/Trace_operator">this wikipedia article</a>
+for further details on the trace operator).
+
+The co-dimension two case is a little more complicated, and in general it is not
+possible to construct a continuous 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 with respect to 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.
+
+This problem has a physical interpretation: harmonic functions, i.e., functions
+that satisfy the Laplace equation, can be thought of as the displacements of a
+membrane whose boundary values are prescribed. The current situation then
+corresponds to finding the shape of a membrane for which not only the
+displacement at the boundary, but also on $\Gamma$ is prescribed. For example,
+if $\Gamma$ is a closed curve in 2d space, then that would model a soap film
+that is held in place by a wire loop along $\partial \Omega$ as well as a second
+loop along $\Gamma$. In cases where $\Gamma$ is a whole area, you can think of
+this as a membrane that is stretched over an obstacle where $\Gamma$ is the
+contact area. (If the contact area is not known we have a different problem --
+called the "obstacle problem" -- which is modeled in step-41.)
+
+As a first example we study the zero Dirichlet boundary condition on
+$\partial\Omega$. The same equations apply if we apply zero Neumann boundary
+conditions on $\partial\Omega$ or a mix of the two.
+
+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 \qquad \forall v \in V(\Omega) \\
+(\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} \qquad \forall q \in Q(\Gamma),
+@f}
+
+where $(\cdot, \cdot)_{\Omega}$ and $(\cdot, \cdot)_{\Gamma}$ represent,
+respectively, $L^2$ scalar products in $\Omega$ and 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$ should therefore 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 discretization of the problem above. One could choose
+matching discretizations, 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.
+
+The first option is clearly more indicated for the simple problem we
+proposed above: it is sufficient to use a single Triangulation for $\Omega$ and
+then impose certain constraints depending $\Gamma$. An example of this approach
+is studied in step-40, where the solution has to stay above an obstacle and this
+is achieved imposing constraints on $\Omega$.
+
+To solve more complex problems, for example one where the domain $\Gamma$ is time
+dependent, the second option could be a more viable solution. Handling
+non aligned meshes is complex by itself: to illustrate how is done we study a
+simple problem.
+
+The technique we describe here is presented in the literature using one of many names:
+the <b>immersed finite element method</b>, the <b>fictitious boundary method</b>, the
+<b>distributed Lagrange multiplier method</b>, and others. The main principle is
+that the discretization 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 applied 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`: on
+this triangulation 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 @p false) or $\delta\psi$ (if the
+parameter `use_displacement` is set to @p 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, either 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`, using the
+FiniteElement `space_fe`, following almost verbatim the approach taken 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, the variational formulation
+of the problem above leaves us 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} &\dealcoloneq& (\nabla v_j, \nabla v_i)_\Omega \qquad i,j=1,\dots,n \\
+C_{\alpha j} &\dealcoloneq& (v_j, q_\alpha)_\Gamma \qquad j=1,\dots,n, \alpha = 1,\dots, m \\\\
+G_{\alpha} &\dealcoloneq& (g, q_\alpha)_\Gamma \qquad \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 appears in the computation of a single entry of
+$C$, is computed on $\Gamma$. As usual in finite elements we split this
+integral into contributions from all cells of the triangulation used to
+discretize $\Gamma$, we transform the integral on $K$ to an integral on the
+reference element $\hat K$, where $F_{K}$ is the mapping from $\hat K$ to $K$,
+and compute the integral on $\hat K$ using a quadrature formula:
+
+\f[
+C_{\alpha j} \dealcoloneq (v_j, 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]
+
+Computing this sum is non-trivial because we have to evaluate $(v_j \circ F_{K})
+(\hat x_i)$. In general, if $\Gamma$ and $\Omega$ are not aligned, the point
+$F_{K}(\hat x_i)$ is completely arbitrary with respect to $\Omega$, and unless
+we figure out a way to interpolate all basis functions of $V_h(\Omega)$ on an
+arbitrary point on $\Omega$, we cannot compute the integral needed for an entry
+of the matrix $C$.
+
+To evaluate $(v_j \circ F_{K}) (\hat x_i)$ the following steps needs to be
+taken (as shown in the picture below):
+
+- For a given cell $K$ in $\Gamma$ compute the real point $y_i \dealcoloneq F_{K} (\hat
+x_i)$, where $x_i$ is one of the quadrature points used for the integral on $K
+\subseteq \Gamma$.
+
+- Find the cell of $\Omega$ in which $y_i$ lies. We shall call this element $T$.
+
+- To evaluate the basis function use the inverse of the mapping $G_T$ that
+transforms the reference element $\hat T$ into the element $T$: $v_j(y_i) = \hat
+v_j \circ G^{-1}_{T} (y_i)$.
+
+<p align="center"> <img
+ src="https://www.dealii.org/images/steps/developer/step-60.C_interpolation.png"
+ alt=""> </p>
+
+The three steps above can be computed by calling, in turn,
+
+- GridTools::find_active_cell_around_point(), followed by
+
+- Mapping::transform_real_to_unit_cell(). We then
+
+- construct a custom Quadrature formula, containing the point in the reference
+ cell and then
+
+- construct an FEValues object, with the given quadrature formula, and
+ initialized with the cell obtained in the first step.
+
+This is what the deal.II function VectorTools::point_value() does when
+evaluating a finite element field (not just a single shape function) at an
+arbitrary point; but this would be inefficient in this case.
+
+A better solution is to use a convenient wrapper to perform the first three
+steps on a collection of points: GridTools::compute_point_locations(). If one is
+actually interested in computing the full coupling matrix, then it is possible
+to call the method NonMatching::create_coupling_mass_matrix(), that performs the
+above steps in an efficient way, reusing all possible data structures, and
+gathering expensive steps together. This is the function we'll be using later in
+this tutorial.
+
+We solve the final saddle point problem by an iterative solver, applied to the
+Schur complement $S$ (whose construction is described, for example, in step-20),
+and we construct $S$ using LinearOperator classes.
+
+
+<h3>The testcase</h3>
+
+The problem we solve here is identical to step-4, with the difference that we
+impose some constraints on an embedded domain $\Gamma$. The tutorial is written
+in a dimension independent way, and in the results section we show how to vary
+both `dim` and `spacedim`.
+
+The tutorial is compiled for `dim` equal to one and `spacedim` equal to two. If
+you want to run the program in embedding dimension `spacedim` equal to three,
+you will most likely want to change the reference domain for $\Gamma$ to be, for
+example, something you read from file, or a closed sphere that you later deform
+to something more interesting.
+
+In the default scenario, $\Gamma$ has co-dimension one, and this tutorial
+program implements the Fictitious Boundary Method. As it turns out, the same
+techniques are used in the Variational Immersed Finite Element Method, and
+the coupling operator $C$ defined above is the same in almost all of these
+non-matching methods.
+
+The embedded domain is assumed to be included in $\Omega$, which we take as the
+unit square $[0,1]^2$. The definition of the fictitious domain $\Gamma$ can be
+modified through the parameter file, and can be given as a mapping from the
+reference interval $[0,1]$ to a curve in $\Omega$.
+
+If the curve is closed, then the results will be similar to running the same
+problem on a grid whose boundary is $\Gamma$. The program will happily run also
+with a non-closed $\Gamma$, although in those cases the mathematical
+formulation of the problem is more difficult, since $\Gamma$ will have a
+boundary by itself that has co-dimension two with respect to the domain
+$\Omega$.
+
+<h3>References</h3>
+
+- Glowinski, R., T.-W. Pan, T.I. Hesla, and D.D. Joseph. 1999. “A Distributed
+ Lagrange Multiplier/fictitious Domain Method for Particulate Flows.”
+ International Journal of Multiphase Flow 25 (5). Pergamon: 755–94.
+
+- Boffi, D., L. Gastaldi, L. Heltai, and C.S. Peskin. 2008. “On the
+ Hyper-Elastic Formulation of the Immersed Boundary Method.” Computer Methods
+ in Applied Mechanics and Engineering 197 (25–28).
+
+- Heltai, L., and F. Costanzo. 2012. “Variational Implementation of Immersed
+ Finite Element Methods.” Computer Methods in Applied Mechanics and Engineering
+ 229–232.
+
--- /dev/null
+<h1>Results</h1>
+
+The directory in which this program is run does not contain a parameter file by
+default. On the other hand, this program wants to read its parameters from a
+file called parameters.prm -- and so, when you execute it the first time, you
+will get an exception that no such file can be found:
+
+@code
+----------------------------------------------------
+Exception on processing:
+
+--------------------------------------------------------
+An error occurred in line <74> of file <../source/base/parameter_acceptor.cc> in function
+ static void dealii::ParameterAcceptor::initialize(const std::string &, const std::string &, const ParameterHandler::OutputStyle, dealii::ParameterHandler &)
+The violated condition was:
+ false
+Additional information:
+ You specified <parameters.prm> as input parameter file, but it does not exist. We created it for you.
+--------------------------------------------------------
+
+Aborting!
+----------------------------------------------------
+@endcode
+
+However, as the error message already states, the code that triggers the
+exception will also generate a parameters.prm file that simply contains the
+default values for all parameters this program cares about. By inspection of the
+parameter file, we see the following:
+
+@code
+subsection Stokes Immersed Problem
+ set Final time = 1
+ set Homogeneous Dirichlet boundary ids = 0, 1, 2, 3
+ set Initial fluid refinement = 4
+ set Initial solid refinement = 4
+ set Particle insertion refinement = 4
+ set Nitsche penalty term = 1000
+ set Number of time steps = 501
+ set Velocity degree = 2
+ set Viscosity = 1
+ subsection Angular velocity
+ set Function constants =
+ set Function expression = t < .5 ? 5 : -5
+ set Variable names = x,y,t
+ end
+ subsection Grid generation
+ set Grid one generator = hyper_cube
+ set Grid one generator arguments = -1: 1: false
+ set Grid two generator = hyper_rectangle
+ set Grid two generator arguments = -.5, -.1: .5, .1: false
+ set Particle grid generator = hyper_ball
+ set Particle grid generator arguments = 0.3, 0.3: 0.1: false
+ end
+ subsection Right hand side
+ set Function constants =
+ set Function expression = 0; 0; 0
+ set Variable names = x,y,t
+ end
+end
+@endcode
+
+If you now run the program, you will get a file called `used_parameters.prm`,
+containing a shorter version of the above parameters (without comments and
+documentation), documenting all parameters that were used to run your program:
+@code
+subsection Stokes Immersed Problem
+ set Final time = 1
+ set Homogeneous Dirichlet boundary ids = 0, 1, 2, 3
+ set Initial fluid refinement = 4
+ set Initial solid refinement = 4
+ set Particle insertion refinement = 4
+ set Nitsche penalty term = 1000
+ set Number of time steps = 501
+ set Velocity degree = 2
+ set Viscosity = 1
+ subsection Angular velocity
+ set Function constants =
+ set Function expression = t < .5 ? 5 : -5
+ set Variable names = x,y,t
+ end
+ subsection Grid generation
+ set Grid one generator = hyper_cube
+ set Grid one generator arguments = -1: 1: false
+ set Grid two generator = hyper_rectangle
+ set Grid two generator arguments = -.5, -.1: .5, .1: false
+ set Particle grid generator = hyper_ball
+ set Particle grid generator arguments = 0.3, 0.3: 0.1: false
+ end
+ subsection Right hand side
+ set Function constants =
+ set Function expression = 0; 0; 0
+ set Variable names = x,y,t
+ end
+end
+@endcode
+
+The rationale behind creating first `parameters.prm` file (the first time the
+program is run) and then a `used_parameters.prm` (every other times you run the
+program), is because you may want to leave most parameters to their default
+values, and only modify a handful of them.
+
+For example, you could use the following (perfectly valid) parameter file with
+this tutorial program:
+@code
+SHORTER VERSION HERE
+@endcode
+
+and you would obtain exactly the same results as in test case 1 below.
+
+<h3> Test case 1: </h3>
+
+For the default problem the value of $u$ on $\Gamma$ is set to the constant $1$:
+this is like imposing a constant Dirichlet boundary condition on $\Gamma$, seen
+as boundary of the portion of $\Omega$ inside $\Gamma$. Similarly on $\partial
+\Omega$ we have zero Dirichlet boundary conditions.
+
+
+<div class="twocolumn" style="width: 80%">
+ <div class="parent">
+ <div class="img" align="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-60.1_no_grid.png"
+ alt = ""
+ width="500">
+ </div>
+ </div>
+ <div class="parent">
+ <div class="img" align="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-60.1_grid.png"
+ alt = ""
+ width="500">
+ </div>
+ </div>
+</div>
+
+The output of the program will look like the following:
+
+@code
+
+@endcode
+
+You may notice that, in terms of CPU time, assembling the coupling system is
+twice as expensive as assembling the standard Poisson system, even though the
+matrix is smaller. This is due to the non-matching nature of the discretization.
+Whether this is acceptable or not, depends on the applications.
+
+If the problem was set in a three-dimensional setting, and the immersed mesh was
+time dependent, it would be much more expensive to recreate the mesh at each
+step rather than use the technique we present here. Moreover, you may be able to
+create a very fast and optimized solver on a uniformly refined square or cubic
+grid, and embed the domain where you want to perform your computation using the
+technique presented here. This would require you to only have a surface
+representatio of your domain (a much cheaper and easier mesh to produce).
+
+To play around a little bit, we are going to complicate a little the fictitious
+domain as well as the boundary conditions we impose on it.
+
+<h3> Test case 2 and 3: </h3>
+
+If we use the following parameter file:
+@code
+subsection Distributed Lagrange<1,2>
+ set Coupling quadrature order = 3
+ set Embedded configuration finite element degree = 1
+ set Embedded space finite element degree = 1
+ set Embedding space finite element degree = 1
+ set Homogeneous Dirichlet boundary ids = 0,1,2,3
+ set Initial embedded space refinement = 8
+ set Initial embedding space refinement = 4
+ set Local refinements steps near embedded domain = 4
+ set Use displacement in embedded interface = false
+ set Verbosity level = 10
+ subsection Embedded configuration
+ set Function constants = R=.3, Cx=.5, Cy=.5, r=.1, w=12
+ set Function expression = (R+r*cos(w*pi*x))*cos(2*pi*x)+Cx; (R+r*cos(w*pi*x))*sin(2*pi*x)+Cy
+ set Variable names = x,y,t
+ end
+ subsection Embedded value
+ set Function constants =
+ set Function expression = x-.5
+ set Variable names = x,y,t
+ end
+ subsection Schur solver control
+ set Log frequency = 1
+ set Log history = false
+ set Log result = true
+ set Max steps = 100000
+ set Reduction = 1.e-12
+ set Tolerance = 1.e-12
+ end
+end
+@endcode
+
+We get a "flowery" looking domain, where we impose a linear boundary condition
+$g=x-.5$. This test shows that the method is actually quite accurate in
+recovering an exactly linear function from its boundary conditions, and even
+though the meshes are not aligned, we obtain a pretty good result.
+
+Replacing $x-.5$ with $2(x-.5)^2-2(y-.5)^2$, i.e., modifying the parameter file
+such that we have
+@code
+ ...
+ subsection Embedded value
+ set Function constants =
+ set Function expression = 2*(x-.5)^2-2*(y-.5)^2
+ set Variable names = x,y,t
+ end
+@endcode
+produces the saddle on the right.
+
+<div class="twocolumn" style="width: 80%">
+ <div class="parent">
+ <div class="img" align="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-60.3_grid.png"
+ alt = ""
+ width="500">
+ </div>
+ </div>
+ <div class="parent">
+ <div class="img" align="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-60.4_grid.png"
+ alt = ""
+ width="500">
+ </div>
+ </div>
+</div>
+
+<a name="extensions"></a>
+<h3>Possibilities for extensions</h3>
+
+<h4> Running with `spacedim` equal to three</h4>
+
+While the current tutorial program is written for `spacedim` equal to two, there
+are only minor changes you have to do in order for the program to run in
+different combinations of dimensions.
+
+If you want to run with `spacedim` equal to three and `dim` equal to two, then
+you will almost certainly want to perform the following changes:
+
+- use a different reference domain for the embedded grid, maybe reading it from
+ a file. It is not possible to construct a smooth closed surface with one
+ single parametrization of a square domain, therefore you'll most likely want
+ to use a reference domain that is topologically equivalent to a the boundary
+ of a sphere.
+
+- use a displacement instead of the deformation to map $\Gamma_0$ into $\Gamma$
+
+<h4> More general domains </h4>
+
+We have seen in other tutorials (for example in step-5 and step-54) how to read
+grids from input files. A nice generalization for this tutorial program would be
+to allow the user to select a grid to read from the parameter file itself,
+instead of hardcoding the mesh type in the tutorial program itself.
+
+<h4> Preconditioner</h4>
+
+At the moment, we have no preconditioner on the Schur complement. This is ok for
+two dimensional problems, where a few hundred iterations bring the residual down
+to the machine precision, but it's not going to work in three dimensions.
+
+It is not obvious what a good preconditioner would be here. The physical problem
+we are solving with the Schur complement, is to associate to the Dirichlet data
+$g$, the value of the Lagrange multiplier $\lambda$. $\lambda$ can be
+interpreted as the *jump* in the normal gradient that needs to be imposed on $u$
+across $\Gamma$, in order to obtain the Dirichlet data $g$.
+
+So $S$ is some sort of Neumann to Dirichlet map, and we would like to have a
+good approximation for the Dirichlet to Neumann map. A possibility would be to
+use a Boundary Element approximation of the problem on $\Gamma$, and construct a
+rough approximation of the hyper-singular operator for the Poisson problem
+associated to $\Gamma$, which is precisely a Dirichlet to Neumann map.
+
+<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 information may not
+be available, because the locally owned part of the $\Omega$ triangulation
+stored on a given processor may not be physically co-located with the locally
+owned part of the $\Gamma$ triangulation stored on the same processor.
+
+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;
+
+- use a distributed triangulation for the embedding space, and a shared
+ triangulation for the emdedded configuration.
+
+The latter strategy is clearly the easiest to implement, as most of the
+functions used in this tutorial program will work unchanged also in the parallel
+case. Of course one could use the reversal strategy (that is, have a distributed
+embedded Triangulation and a shared embedding Triangulation).
+
+However, this strategy is most likely going to be more expensive, since by
+definition the embedding grid is larger than the embedded grid, and it makes
+more sense to distribute the largest of the two grids, maintaining the smallest
+one shared among all processors.
projects / dealii.git / commitdiff