]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added possibility to read parameter from command line.
authorLuca Heltai <luca.heltai@sissa.it>
Wed, 2 May 2018 23:09:29 +0000 (01:09 +0200)
committerMatthias Maier <tamiko@43-1.org>
Tue, 8 May 2018 04:44:55 +0000 (23:44 -0500)
examples/step-60/doc/intro.dox
examples/step-60/doc/results.dox
examples/step-60/step-60.cc

index 32443edac50fa121d6ace5fcb21f7ec260ff3eb8..3fb137029e3a4b2f344298b9cf3a9a5e1a7a6e65 100644 (file)
@@ -41,7 +41,9 @@ 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$.
+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
@@ -66,6 +68,19 @@ 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.
@@ -84,7 +99,7 @@ Given a sufficiently regular function $g$ on $\Gamma$, find the solution $u$ to
 (\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)$
@@ -95,16 +110,23 @@ 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)$.
+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
@@ -129,8 +151,8 @@ for any point $x \in \Gamma_0$, its coordinate $\psi(x)$ in $\Omega$;
 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
@@ -146,7 +168,7 @@ 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
+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
@@ -160,8 +182,9 @@ Q_h(\Gamma) = \text{span} \{q_i\}_{i=1}^m
 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}
@@ -194,11 +217,11 @@ 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 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}
@@ -213,7 +236,8 @@ 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$.
 
-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
@@ -226,21 +250,26 @@ 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"
+  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
index ddaafdff0c8193cd1d4cbf470711de406f8b8068..d582d0ac3285b8ee7e6bdf35b9d38874f149f7ea 100644 (file)
@@ -1,7 +1,8 @@
 <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: 
@@ -321,7 +322,7 @@ such that we have
     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>
@@ -368,16 +369,23 @@ of this hot air balloon shape. Why? What's the difference with the following?
   </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>
 
index c3c24f016c9095f6999371d7bf4e29f7f6599e6b..404e1444fc4e51af100c765d5aa1b3d8b27b6443 100644 (file)
@@ -398,11 +398,9 @@ namespace Step60
 
     SparsityPattern           stiffness_sparsity;
     SparsityPattern           coupling_sparsity;
-    SparsityPattern           embedded_sparsity;
 
     SparseMatrix<double>      stiffness_matrix;
     SparseMatrix<double>      coupling_matrix;
-    SparseMatrix<double>      embedded_stiffness_matrix;
 
     ConstraintMatrix          constraints;
 
@@ -876,10 +874,6 @@ namespace Step60
                   (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());
@@ -924,24 +918,23 @@ namespace Step60
     {
       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,
@@ -968,7 +961,6 @@ namespace Step60
 
     // 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);
 
@@ -977,7 +969,7 @@ namespace Step60
     // 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;
 
@@ -1023,7 +1015,7 @@ namespace Step60
 
 
 
-int main()
+int main(int argc, char **argv)
 {
   try
     {
@@ -1032,15 +1024,26 @@ int main()
 
       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)

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.