From 7ef46783dd00b3dcf1b29374adaa07ecd4620af0 Mon Sep 17 00:00:00 2001 From: Giovanni Alzetta Date: Wed, 2 May 2018 14:17:20 +0200 Subject: [PATCH] Corrected refinement --- examples/step-60/doc/intro.dox | 39 +++++++++++++------------- examples/step-60/doc/results.dox | 9 ++---- examples/step-60/step-60.cc | 48 +++++++++++++++++++------------- 3 files changed, 52 insertions(+), 44 deletions(-) diff --git a/examples/step-60/doc/intro.dox b/examples/step-60/doc/intro.dox index 7fd050893c..dd02773084 100644 --- a/examples/step-60/doc/intro.dox +++ b/examples/step-60/doc/intro.dox @@ -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] diff --git a/examples/step-60/doc/results.dox b/examples/step-60/doc/results.dox index 9c7a88bb4d..85f416278d 100644 --- a/examples/step-60/doc/results.dox +++ b/examples/step-60/doc/results.dox @@ -3,15 +3,12 @@

Test case 1:

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:

Possibilities for extensions

-

Different Parameters

+Add something

Parallel Code

@@ -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. diff --git a/examples/step-60/step-60.cc b/examples/step-60/step-60.cc index 5483e95e1b..2c05c8c39c 100644 --- a/examples/step-60/step-60.cc +++ b/examples/step-60/step-60.cc @@ -24,8 +24,6 @@ #include #include -#include - // 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 + #include #include #include @@ -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 > embedded_configuration_dh; Vector 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::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 > 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 -- 2.39.5