From c95e2cfe0da59bf05552f87a07c6f515809dda16 Mon Sep 17 00:00:00 2001 From: kanschat Date: Thu, 20 Sep 2012 12:07:16 +0000 Subject: [PATCH] add elasticity example and highlight the differences between the different approaches git-svn-id: https://svn.dealii.org/trunk@26560 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/doxygen/headers/vector_valued.h | 331 ++++++++++++++------ 1 file changed, 238 insertions(+), 93 deletions(-) diff --git a/deal.II/doc/doxygen/headers/vector_valued.h b/deal.II/doc/doxygen/headers/vector_valued.h index 367b222e5f..6045cc59fb 100644 --- a/deal.II/doc/doxygen/headers/vector_valued.h +++ b/deal.II/doc/doxygen/headers/vector_valued.h @@ -57,23 +57,151 @@ * * * @anchor VVPhilosophy - *

A philosophical view of vector-valued problems

+ *

Examples of vector-valued problems

* * The way one deals systematically with vector-valued problems is not - * fundamentally different from scalar problems in that one needs to come up - * first with a weak (variational) formulation of the problem that takes into - * account all the solution variables. To understand how to do that, let us - * consider a simple example, the mixed Laplace equations discussed in - * step-20: + * fundamentally different from scalar problems: first, we need a weak + * (variational) formulation of the problem that takes into account + * all the solution variables. After we did so, generating the system + * matrix and solving the linear system follows the same outlines that + * we are used to already. + * + *

Linear elasticity

+ * + * Let us take for example the elasticity problem from step-8 and even + * simplify it by choosing $\lambda = 0$ and $\mu = 1$ to highlight the important concepts. Therefore, + * let consider the following weak formulation: find $\mathbf u \in + * \mathbf V = H^1_0(\Omega; \mathbb R^3)$ such that for all $\mathbf + * v\in V$ holds + * @f[ + * a(u,v) \equiv 2\int_{\Omega} \mathbf D\mathbf u : \mathbf D\mathbf + * v\,dx = \int_\Omega \mathbf f\cdot \mathbf v \,dx. + * @f] + * Here, D denotes the symmetric gradient defined by + * $\mathbf Du = \tfrac12 (\nabla \mathbf u + (\nabla \mathbf u)^T)$ + * and the colon indicates double contraction of two tensors of rank 2 + * (the Frobenius inner product). This bilinear form looks indeed very + * much like the bilinear form of the Poisson problem in step-3. The + * only differences are + *
    + *
  1. We replaced the gradient operator by the symmetric gradient; + * this is actually not a significant difference, and everything said + * here is true if your replace $\mathbf D$ by $\nabla$. Indeed, let + * us do this to simplify the discussion: + * @f[ + * a(u,v) \equiv \int_{\Omega} \nabla\mathbf u : \nabla\mathbf + * v\,dx = \int_\Omega \mathbf f\cdot \mathbf v \,dx. + * @f] + * Note though, that this system is not very exciting, since we could + * solve for the three components of u separately. + * + *
  2. The trial and test functions are now from the space + * $H^1_0(\Omega; \mathbb R^3)$, which can be viewed as three copies + * of the scalar space $H^1_0(\Omega)$. And this is exactly, how we + * are going to implement this space below, using FESystem. + * + * But for now, let us look at the system a little more + * closely. First, let us exploit that + * u=(u1,u2,u3)T + * and v accordingly. Then, we can write the simplified equation in + * coordinates as + * @f[ + * a(u,v) = \int_\Omega \bigl(\nabla u_1\cdot \nabla v_1 + +\nabla u_2\cdot \nabla v_2+\nabla u_3\cdot \nabla v_3\bigr)\,dx + = \int_\Omega \bigl(f_1v_1 + f_2 v_2 + f_3 v_3\bigr)\,dx. + * @f] + * We see, that this is just three copies of the bilinear form of the + * Laplacian, one applied to each component (this is where the + * formulation with the $\mathbf D$ is more exciting, and we want to + * derive a framework that applies to that one as well). We can make + * this weak form a system of differential equations again by choosing + * special test functions: first, choose + * v=(v1,0,0)T, then + * v=(0,v2,0)T, and finally + * v=(0,0,v3)T. writing the outcomes below + * each other, we obtain the system + * @f[ + * \arraycolsep1pt + * \begin{matrix} + * (\nabla u_1,\nabla v_1) &&& = (f_1, v_1) + * \\ + * & (\nabla u_2,\nabla v_2) && = (f_2, v_2) + * \\ + * && (\nabla u_3,\nabla v_3) & = (f_3, v_3) + * @f] + * where we used the standard inner product notation $(\mathbf + * f,\mathbf g) = + * \int_\Omega \mathbf f \cdot \mathbf g \,dx$. It is important for our understanding, that + * we keep in mind that the latter form as a system of PDE is + * completely equivalent to the original definition of the bilinear + * form a(u,v), which does not immediately + * exhibit this system structure. Let us close by writing the full + * system of the elastic equation with symmetric gradient D: + * @f[ + * \arraycolsep1pt + * \begin{matrix} + * (\nabla u_1,\nabla v_1) + (\partial_1 u_1,\partial_1 v_1) + * & (\partial_1 u_2,\partial_2 v_1) + * & (\partial_1 u_3,\partial_3 v_1) + * & = (f_1, v_1) + * \\ + * (\partial_2 u_1,\partial_1 v_2) + * & (\nabla u_2,\nabla v_2) + (\partial_2 u_2,\partial_2 v_2) + * & (\partial_2 u_3,\partial_3 v_2) + * & = (f_2, v_2) + * \\ + * (\partial_3 u_1,\partial_1 v_3) + * & (\partial_3 u_2,\partial_2 v_3) + * & (\nabla u_3,\nabla v_3) + (\partial_3 u_3,\partial_3 v_3) + * & = (f_3, v_3) + * \end{matrix}. + * @f] + * Very formally, if we believe in operator valued matrices, we can + * rewrite this in the form vTAu = vTf or + * @f[ + * \begin{pmatrix} + * v_1 \\ v_2 \\ v_3 + * \end{pmatrix}^T + * \begin{pmatrix} + * (\nabla \cdot,\nabla \cdot) + (\partial_1 \cdot,\partial_1 \cdot) + * & (\partial_1 \cdot,\partial_2 \cdot) + * & (\partial_1 \cdot,\partial_3 \cdot) + * \\ + * (\partial_2 \cdot,\partial_1 \cdot) + * & (\nabla \cdot,\nabla \cdot) + (\partial_2 \cdot,\partial_2 \cdot) + * & (\partial_2 \cdot,\partial_3 \cdot) + * \\ + * (\partial_3 \cdot,\partial_1 \cdot) + * & (\partial_3 \cdot,\partial_2 \cdot) + * & (\nabla \cdot,\nabla \cdot) + (\partial_3 \cdot,\partial_3 \cdot) + * \end{pmatrix} + * \begin{pmatrix} + * u_1 \\ u_2 \\ u_3 + * \end{pmatrix} + * = + * \begin{pmatrix} + * v_1 \\ v_2 \\ v_3 + * \end{pmatrix}^T + * \begin{pmatrix} f_1 \\ f_2 \\ f_3\end{pmatrix} + * @f] + * + *

    Mixed elliptic problems

    + * Now, let us + * consider a more complex example, the mixed Laplace equations discussed in + * step-20 in three dimensions: @f{eqnarray*} \textbf{u} + \nabla p &=& 0, \\ -\textrm{div}\; \textbf{u} &=& f, @f} * - * Here, we have dim+1 solution components: the scalar pressure - * $p$ and the vector-valued velocity $\textbf u$ with dim vector - * components. + * Here, we have four solution components: the scalar pressure + * $p \in L^2(\Omega)$ and the vector-valued velocity $\mathbf u \in + * \mathbf V + * = H^{\text{div}}_0(\Omega)$ with three vector + * components. Note as important difference to the previous example, + * that the vector space V is not just simply a copy of three + * identical spaces/ * * A systematic way to get a weak or variational form for this and other * vector problems is to first consider it as a problem where the operators @@ -99,7 +227,7 @@ \begin{array}{c} \mathbf u \\ p \end{array} \right) @f} - * indeed has dim+1 components. We note that we could change the + * indeed has four components. We note that we could change the * ordering of the solution components $\textbf u$ and $p$ inside $U$ if we * also change columns of the matrix operator. * @@ -155,13 +283,8 @@ = (q,f), @f} - * where we have defined the scalar product - * $(\mathbf v, \mathbf u) = \int_\Omega \mathbf v(x) - * \cdot \mathbf u(x) \; dx$, and similarly if both parts of the scalar - * product are scalar-valued functions (e.g. the pressure) rather - * than vector-valued onces (like the velocity). - * - * We get the final form by integrating by part the second term: +* +* We get the final form by integrating by part the second term: @f{eqnarray*} (\mathbf v, \mathbf u) - @@ -181,83 +304,122 @@ * @anchor VVFEs *

    Describing finite element spaces

    * - * Once we have settled on this description, we need to find a way to describe - * the vector-valued finite element space from which we draw solution and test + * Once we have settled on a bilinear form and a functional setting, we need to find a way to describe + * the vector-valued finite element spaces from which we draw solution and test * functions. This is where the FESystem class comes in: it composes - * vector-valued finite element spaces from simpler ones. For example, if we - * were to attempt to use $Q_1$ elements for all dim components - * of $\mathbf u$ and the one pressure component $p$, we could use the - * following object: + * vector-valued finite element spaces from simpler ones. + * In the example of the elasticity problem, we need dim + * copies of the same element, for instance * @code - * FESystem finite_element (FE_Q(1), dim, - * FE_Q(1), 1); + * FESystem elasticity_element (FE_Q(1), dim); * @endcode - * - * This means that the final finite element will consist of dim - * components made up of FE_Q elements of degree 1, and another one also of - * degree 1. Of course, a simpler (and more efficient) way to achieve the same - * is to use the following form instead: + * This will generate a vector valued space of dimension + * dim, where each component is a + * continuous bilinear element of type FE_Q. It will have dim times + * as many basis functions as the corresponding FE_Q, and each of + * these basis functions is a basis function of FE_Q, lifted into one + * of the components of the vector. + * + * For the mixed Laplacian, the situation is more complex. First, we + * have to settle on a pair of discrete spaces $\mathbf V_h \times Q_h + * \subset H^{\text{div}}_0(\Omega) \times L^2_0(\Omega)$. One option + * would be the stable Raviart-Thomas pair * @code - * FESystem finite_element (FE_Q(1), dim+1); + * FESystem rt_element (FE_RaviartThomas(1), 1, + * FE_DGQ(1), 1); * @endcode - * Another way (also not efficient, but making it clear which components - * belong to which element) would be + * The first element in this system is already a vector valued + * element of dimension dim, while the second is a + * regular scalar element. + * + * Alternatively to using the stable Raviart-Thomas pair, we could + * consider a stabilized formulation for the mixed Laplacian, for + * instance the LDG method. There, we have the option of using the + * same spaces for velocity components and pressure, namely * @code - * FESystem finite_element (FESystem(FE_Q(1), dim), 1, - * FE_Q(1), 1); + * FESystem ldg_convoluted_element_1 (FE_DGQ(1), dim+1); * @endcode - * meaning that we first create a vector element with dim - * components, each consisting of FE_Q elements of order 1. We then couple - * one copy of this already vector-valued element to a single scalar - * copy of the FE_Q element of order 1 which will describe the pressure - * component. - * - * As it turns out these (equivalent) choices do not lead to a stable scheme - * for the mixed Laplace equation. In step-20, we therefore use - * a Raviart-Thomas element for the velocities. What exactly this means may be - * of less concern for now except that the FE_RaviartThomas class describes - * elements that already have dim components. For the pressure, - * we use piecewise bi-/trilinear elements that may be discontinuous between - * cells; this is done using the FE_DGQ class. The combined element will then - * be described by + * This system just has dim+1 equal copies of the same + * discontinuous element, which not really reflects the structure of + * the system. Therefore, we prefer * @code - * FESystem finite_element (FE_RaviartThomas(1), 1, - * FE_DGQ(1), 1); + * FESystem ldg_equal_element (FESystem(FE_DGQ(1), dim), 1, + * FE_DGQ(1), 1); * @endcode - * i.e. we combine a single copy of the Raviart-Thomas element with a single - * copy of the element used for the pressure $p$. - * - * In this manner, we can combine the whole vector-valued element from its - * individual components. However, it is worth pointing out that whichever - * way you construct your finite element object, this object doesn't really - * know what it represents. For example, for + * Here, we have a system of two elements, one vector-valued and one + * scalar, very much like with the rt_element. Indeed, in + * many codes, the two can be interchanged. This element also allows + * us easily to switch to an LDG method with lower order approximation + * in the velocity, namely * @code - * FESystem finite_element (FE_Q(1), dim, - * FE_Q(1), 1); + * FESystem ldg_unequal_element (FESystem(FE_DGQ(1), dim), 1, + * FE_DGQ(2), 1); * @endcode - * the constructed object really only knows that it has dim+1 - * vector components. It has no notion, however, which of these components - * represent scalar fields (e.g. temperature, pressure, concentration) and/or - * if any of its components are parts of vector fields (velocities, - * displacements) or tensors (e.g. stresses). As a consequence, the FEValues - * object we use below to evaluate finite element shape functions at - * quadrature points only knows that it has a finite element with a number - * of vector components, but doesn't know how to group them. We will show - * how to give these components a logical connection using the - * FEValuesExtractors classes. - * + * It must be pointed out, + * that this element is different from + * @code + * FESystem ldg_convoluted_element_2 (FE_DGQ(1), dim, + * FE_DGQ(2), 1); + * @endcode + * While the constructor call is very similar to + * rt_element, the result actually resembles more + * ldg_convoluted_element_1 in that this element produces + * dim+1 independent components. + * A more detailed comparison of the resulting FESystem objects is below. + * + *

    Internal structure of FESystem

    + * + * FESystem has a few internal variables which reflect the internal + * structure set up by the constructor. These can then also be used by + * application programs to give structure to matrix assembling and + * linear algebra. We give the names and values of these variables for + * the examples above in the following table: + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
    System ElementFiniteElementData::n_blocks()FiniteElementData::n_components()FiniteElement::n_base_elements()
    elasticity_elementdimdim1
    rt_element2dim+12
    ldg_equal_element2dim+12
    ldg_convoluted_element_1dim+1dim+11
    ldg_convoluted_element_2dim+1dim+12
    + * + * From this table, it is clear that the FESystem reflects a lot of + * the structure of the system of differential equations in the cases + * of the rt_element and the + * ldg_equal_element, in that we have a vector valued and + * a scalar variable. On the other hand, the convoluted elements do + * not have this structure and we have to reconstruct it somehow when + * assembling systems, as described below. + * + * At this point, it is important to note that nesting of two FESystem + * object can give the whole FESystem a richer structure than just + * concatenating them. This structure can be exploited by application + * programs, but is not automatically so. * * @anchor VVAssembling *

    Assembling linear systems

    - * * The next step is to assemble the linear system. How to do this for the * simple case of a scalar problem has been shown in many tutorial programs, * starting with step-3. Here we will show how to do it for - * vector problems. + * vector problems. Corresponding to the different characterizations + * of weak formulations above and the different system elements + * created, we have a few options which are outlined below. * * The whole concept is probably best explained by showing an example * illustrating how the local contribution of a cell to the weak form of above - * mixed Laplace equations could be assembled. This is essentially how + * mixed Laplace equations could be assembled. + * + *

    A single FEValues and FEValuesExtractors

    + * This is essentially how * step-20 does it: * @code const FEValuesExtractors::Vector velocities (0); @@ -450,33 +612,16 @@ phi_i_div * phi_j_div + mu_values[q_point] * - scalar_product(phi_i_grad, phi_j_grad) + double_contract(phi_i_grad, phi_j_grad) + mu_values[q_point] * - scalar_product(phi_i_grad, transpose(phi_j_grad)) + double_contract(phi_i_grad, transpose(phi_j_grad)) ) * fe_values.JxW(q_point); } } * @endcode * - * The scalar product between two tensors used in this bilinear form is - * implemented as follows: - * - * @code -template -double -scalar_product (const Tensor<2,dim> &u, - const Tensor<2,dim> &v) -{ - double tmp = 0; - for (unsigned int i=0; i