From 3ef1d75b6147fc16da95d37fbd2c345448712955 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Wed, 15 Apr 2020 01:09:21 +0200 Subject: [PATCH] Copy doc from step-60. --- examples/step-70/doc/builds-on | 1 + examples/step-70/doc/intro.dox | 334 +++++++++++++++++++++++++++++++ examples/step-70/doc/kind | 1 + examples/step-70/doc/results.dox | 305 ++++++++++++++++++++++++++++ examples/step-70/doc/tooltip | 1 + 5 files changed, 642 insertions(+) create mode 100644 examples/step-70/doc/builds-on create mode 100644 examples/step-70/doc/intro.dox create mode 100644 examples/step-70/doc/kind create mode 100644 examples/step-70/doc/results.dox create mode 100644 examples/step-70/doc/tooltip diff --git a/examples/step-70/doc/builds-on b/examples/step-70/doc/builds-on new file mode 100644 index 0000000000..16313b3a57 --- /dev/null +++ b/examples/step-70/doc/builds-on @@ -0,0 +1 @@ +step-6, step-22, step-32, step-55, step-60 diff --git a/examples/step-70/doc/intro.dox b/examples/step-70/doc/intro.dox new file mode 100644 index 0000000000..4cd9edc183 --- /dev/null +++ b/examples/step-70/doc/intro.dox @@ -0,0 +1,334 @@ +
+ +This program was contributed by Luca Heltai (SISSA, Trieste) and Bruno Blais (Polytechnique Montréal) + + +Change this!!!! + +@dealiiTutorialDOI{10.5281/zenodo.1243280,https://zenodo.org/badge/DOI/10.5281/zenodo.1243280.svg} + + +

Introduction

+ +

Smart ass title here.

+ + +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 this wikipedia article +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 immersed finite element method, the fictitious boundary method, the +distributed Lagrange multiplier method, 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)$. + +

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

The testcase

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

References

+ +- 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. + diff --git a/examples/step-70/doc/kind b/examples/step-70/doc/kind new file mode 100644 index 0000000000..c1d9154931 --- /dev/null +++ b/examples/step-70/doc/kind @@ -0,0 +1 @@ +techniques diff --git a/examples/step-70/doc/results.dox b/examples/step-70/doc/results.dox new file mode 100644 index 0000000000..7e62cede86 --- /dev/null +++ b/examples/step-70/doc/results.dox @@ -0,0 +1,305 @@ +

Results

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

Test case 1:

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

Test case 2 and 3:

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

Possibilities for extensions

+ +

Running with `spacedim` equal to three

+ +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$ + +

More general domains

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

Preconditioner

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

Parallel Code

+ +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. diff --git a/examples/step-70/doc/tooltip b/examples/step-70/doc/tooltip new file mode 100644 index 0000000000..5d847a3c3d --- /dev/null +++ b/examples/step-70/doc/tooltip @@ -0,0 +1 @@ +A fluid structure interaction problem, using a penalty term -- 2.39.5