]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Corrected refinement
authorGiovanni Alzetta <giovannialzetta@hotmail.it>
Wed, 2 May 2018 12:17:20 +0000 (14:17 +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 7fd050893cbcf6db6b6eb2b7e74bc20a8ce0aca5..dd02773084d32af52621e2b4e769ac6bf9624b01 100644 (file)
@@ -63,7 +63,7 @@ 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
+$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
@@ -80,21 +80,21 @@ 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}$ (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.
 
@@ -102,10 +102,10 @@ 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 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
@@ -119,30 +119,31 @@ 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$;
+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[
@@ -178,7 +179,7 @@ where
 
 @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}
 
@@ -188,14 +189,14 @@ 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 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]
index 9c7a88bb4d56397434bee67b75e5dfbf78148e01..85f416278d32d9c4e8c98765536f320021666e5a 100644 (file)
@@ -3,15 +3,12 @@
 <h3> Test case 1: </h3>
 
 For the default problem the value of u on Gamma is 1 and on $\partial\Omega$
-is 0. This means we expect the following solution:
-
-
-
+is 0. In fact this is the solution:
 
 <a name="extensions"></a>
 <h3>Possibilities for extensions</h3>
 
-<h4> Different Parameters</h4>
+Add something
 
 <h4> Parallel Code </h4>
 
@@ -33,5 +30,5 @@ Various strategies can be implemented to tackle this problem:
 - make use of a shared triangulation and a distributed triangulation
 
 The latter strategy is clearly the easier to implement, as all
-the function used in this tutorial program can work lettin $\Omega$
+the function used in this tutorial program can work letting $\Omega$
 be distributed and $\Gamma$ be a shared triangulation.
index 5483e95e1b56ee4806d6b1d6ce1ee78b1256dc5c..2c05c8c39cc77288c924f4593acf368ee7f9fe31 100644 (file)
@@ -24,8 +24,6 @@
 #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
@@ -62,6 +60,8 @@
 // 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>
@@ -229,13 +229,14 @@ namespace Step60
       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
@@ -323,7 +324,7 @@ namespace Step60
     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
@@ -588,7 +589,7 @@ namespace Step60
     // 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.
@@ -747,17 +748,17 @@ namespace Step60
     // 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.
@@ -770,7 +771,7 @@ namespace Step60
     // 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
@@ -780,7 +781,8 @@ namespace Step60
     //
     // 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,
@@ -800,7 +802,15 @@ namespace Step60
                                      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,
@@ -949,8 +959,8 @@ namespace Step60
     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);
@@ -958,12 +968,12 @@ namespace Step60
     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;
 
@@ -1018,7 +1028,7 @@ int main()
 
       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

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.