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$.
+co-dimension one curves and surfaces $\Gamma$ embedded in $\Omega$ (read on
+<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 continous trace operator, not even from
$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.
(\gamma u, q)_{\Gamma} &=& (g,q)_{\Gamma} & \forall q \in Q(\Gamma),
@f}
-where $(\cdot, \cdot){\Omega}$ and $(\cdot, \cdot){\Gamma}$ represent,
+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)$
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)$.
+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.
-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 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
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`, and on
-this domain, we construct a finite dimensional space (`embedded_configuration_dh`)
+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
vector.
In the embedding space, a standard finite dimensional space `space_dh` is
-constructed on the embedding grid `space_grid` (`space_dh`), using the
+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
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:
+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}
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 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 mapping from $\hat K$ to $K$, and compute the integral
-on $\hat K$ using a quadrature formula:
+$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 tranform 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} := (v_j, q_\alpha)_\Gamma = \sum_{K\in \Gamma} \int_{\hat K}
arbitrary point on $\Omega$, we cannot compute the integral needed for an entry
of the matrix $C$.
-The process is described by the following algorithm (and by the picture below):
+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 \coloneqq F_{K}
(\hat x_i)$, where $x_i$ is one of the quadrature points used for the integral
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"
+ 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:
+The three steps above can be computed using by calling, in turn:
- GridTools::find_active_cell_around_point, followed by
- Mapping::tranform_real_to_unit_cell
- construct a custom Quadrature formula, containing the point in the reference
-cell
+ cell
+
+- construct an FEValues object, with the given quadrature formula, and initialized
+ with the cell obtained in the first step.
-- 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 is would be inefficient in this case.
-The deal.II library offers a convenient wrapper to perform the first three steps
+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
<h1>Results</h1>
-The first time the code is run, an exception is thrown and nothing is actually
-computed. You'll get the following output:
+Currently there is no `parameters.prm` file: objects can't be initialized.
+For this reason the first time the code is run, an exception is thrown
+and nothing is actually computed. You'll get the following output:
@code
----------------------------------------------------
Exception on processing:
set Variable names = x,y,t
end
@endcode
-produce the saddle on the right.
+produces the saddle on the right.
<TABLE WIDTH="60%" ALIGN="center">
<tr>
</tr>
</table>
-To solve this problem we need to improve the way in which adaptive refinement is
-performed. Currently the vector of support points contains the vertices of
-$\Gamma$, if the finite dimensional space for the configuration is of degree
-one: these are too far away to make $\Omega$ refine properly and simply
-increasing the refinements of $\Gamma$ will not be enough, the conjugate
-gradient won't converge.
+<h4> Preconditioner </h4>
-As you modify the code you shall notice that, for instance, with complex
-geometries there shall be other convergence problems: parameters are handy, but do
-not solve all problems!
+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 need 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>
SparsityPattern stiffness_sparsity;
SparsityPattern coupling_sparsity;
- SparsityPattern embedded_sparsity;
SparseMatrix<double> stiffness_matrix;
SparseMatrix<double> coupling_matrix;
- SparseMatrix<double> embedded_stiffness_matrix;
ConstraintMatrix constraints;
(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());
{
TimerOutput::Scope timer_section(monitor, "Assemble system");
- // Embedding stiffness matrix $K$
+ // Embedding stiffness matrix $K$, and the right hand side $G$.
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);
+ VectorTools::create_right_hand_side(*embedded_mapping,
+ *embedded_dh,
+ QGauss<dim>(2*embedded_fe->degree+1),
+ 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
+ // 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,
// 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);
// Using the Schur complement method
auto S = C*K_inv*Ct;
SolverCG<Vector<double> > solver_cg(schur_solver_control);
- auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() );
+ auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity());
lambda = S_inv * embedded_rhs;
-int main()
+int main(int argc, char **argv)
{
try
{
const unsigned int dim=1, spacedim=2;
- // Differently to what happens in other tutorial programs, here we the
+ // Differently to what happens in other tutorial programs, here we use
// 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.
+ //
+ // We check if the user has specified a parameter file name to use when
+ // the program was launched. If so, try to read that parameter file,
+ // otherwise, try to read the file "parameters.prm".
+ //
+ // If the parameter file that was specified (implicitly or explicitly)
+ // does not exist, ParameterAcceptor::initialize will create one for you,
+ // and exit the program.
DistributedLagrangeProblem<dim, spacedim>::DistributedLagrangeProblemParameters parameters;
DistributedLagrangeProblem<dim, spacedim> problem(parameters);
- ParameterAcceptor::initialize("parameters.prm", "used_parameters.prm");
+ std::string parameter_file = "parameters.prm";
+ if (argc > 1)
+ parameter_file = argv[1];
+ ParameterAcceptor::initialize(parameter_file, "used_parameters.prm");
problem.run();
}
catch (std::exception &exc)
projects / dealii.git / commitdiff