@f}
This is a constrained problem, where we are looking for a harmonic function
-$u$, that satisfies homogeneous boundary conditions on $\partial\Omega$, subject
+$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
(\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$).
+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$, therefore, should be either in $H^1(\Gamma)$ (for the
+The function $g$ should, 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
+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.
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 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 discretisation of the two grids and of the two finite element spaces
+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
- 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$;
+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
+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`)
-to describe either the deformation or the displacement, through a FiniteElement
+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).
+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, a MappingFEField or a
-MappingQEulerian object are initialized with the `embedded_configuration` vector.
+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` (`space_dh`), using the
-FiniteElement `space_fe`, following almost verbatim what was done in step-6.
+FiniteElement `space_fe`, following almost verbatim the approach taken in step-6.
We represent the discretizations of the spaces $V$ and $Q$ with
\f[
@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 \\\\
+C_{\alpha j} &:=& (v_j, q_\alpha)_\Gamma &j=1,\dots,n, \alpha = 1,\dots, m \\\\
G_{\alpha} &:=& (g, q_\alpha)_\Gamma & \alpha = 1,\dots, m.
@f}
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$,
+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 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}
+C_{\alpha j} := (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]
#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
// deal.II classes, and deriving our own parameter classes directly from
// ParameterAcceptor.
+#include <deal.II/base/parameter_acceptor.h>
+
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
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
+ // grid $\Gamma$ is handled through 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.
+ // a higher refinement of such cells might improve quality of our
+ // computations.
// 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.
+ // a vertex of the embedded grid and its neighbors, 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
std::unique_ptr<DoFHandler<dim, spacedim> > embedded_configuration_dh;
Vector<double> embedded_configuration;
- // The ParameterAcceptorProxy class, is a "transparent" wrapper, derived
+ // 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
// 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
+ // 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.
// 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
+ // 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
+ // In fact, there are algorithms 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
+ // 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.
// 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
+ // GridTools::compute_point_locations allow 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
//
// 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.
+ // the embedded space. This allows for a simple adaptive refinement strategy:
+ // refining these cells and their neighbors.
//
// 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,
support_points);
const auto &cells = std::get<0>(point_locations);
for (auto cell : cells)
- cell->set_refine_flag();
+ {
+ cell->set_refine_flag();
+ for (unsigned int face_no=0; face_no<GeometryInfo<spacedim>::faces_per_cell; ++face_no)
+ if (! cell->face(face_no)->at_boundary())
+ {
+ auto neighbor = cell->neighbor(face_no);
+ neighbor->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,
K_inv_umfpack.initialize(stiffness_matrix);
// Same thing, for the embedded space
- SparseDirectUMFPACK A_inv_umfpack;
- A_inv_umfpack.initialize(embedded_stiffness_matrix);
+// 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 C = transpose_operator(Ct);
auto K_inv = linear_operator(K, K_inv_umfpack);
- auto A_inv = linear_operator(A, A_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);
+ auto S_inv = inverse_operator(S, solver_cg, PreconditionIdentity() );//A_inv);
lambda = S_inv * embedded_rhs;
const unsigned int dim=1, spacedim=2;
- // Differently to what happens in other tutorial programs, here we the the
+ // Differently to what happens in other tutorial programs, here we 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
projects / dealii.git / commitdiff