 |
Reference documentation for deal.II version 9.3.3
|
|
Reference documentation for deal.II version 9.3.3
|
This tutorial depends on step-12.
| Table of contents | |
|---|---|
This program was contributed by Natasha Sharma, Guido Kanschat, Timo Heister, Wolfgang Bangerth, and Zhuoran Wang.
The first author would like to acknowledge the support of NSF Grant No. DMS-1520862. Timo Heister and Wolfgang Bangerth acknowledge support through NSF awards DMS-1821210, EAR-1550901, and OAC-1835673.
This program deals with the biharmonic equation,
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega. \end{align*}
This equation appears in the modeling of thin structures such as roofs of stadiums. These objects are of course in reality three-dimensional with a large aspect ratio of lateral extent to perpendicular thickness, but one can often very accurately model these structures as two dimensional by making assumptions about how internal forces vary in the perpendicular direction. These assumptions lead to the equation above.
The model typically comes in two different kinds, depending on what kinds of boundary conditions are imposed. The first case,
\begin{align*} u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \end{align*}
corresponds to the edges of the thin structure attached to the top of a wall of height \(g(\mathbf x)\) in such a way that the bending forces that act on the structure are \(h(\mathbf x)\); in most physical situations, one will have \(h=0\), corresponding to the structure simply sitting atop the wall.
In the second possible case of boundary values, one would have
\begin{align*} u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \frac{\partial u(\mathbf x)}{\partial \mathbf n} &= j(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
This corresponds to a "clamped" structure for which a nonzero \(j(\mathbf x)\) implies a certain angle against the horizontal.
As with Dirichlet and Neumann boundary conditions for the Laplace equation, it is of course possible to have one kind of boundary conditions on one part of the boundary, and the other on the remainder.
The fundamental issue with the equation is that it takes four derivatives of the solution. In the case of the Laplace equation we treated in step-3, step-4, and several other tutorial programs, one multiplies by a test function, integrates, integrates by parts, and ends up with only one derivative on both the test function and trial function – something one can do with functions that are continuous globally, but may have kinks at the interfaces between cells: The derivative may not be defined at the interfaces, but that is on a lower-dimensional manifold (and so doesn't show up in the integrated value).
But for the biharmonic equation, if one followed the same procedure using integrals over the entire domain (i.e., the union of all cells), one would end up with two derivatives on the test functions and trial functions each. If one were to use the usual piecewise polynomial functions with their kinks on cell interfaces, the first derivative would yield a discontinuous gradient, and the second derivative with delta functions on the interfaces – but because both the second derivatives of the test functions and of the trial functions yield a delta function, we would try to integrate the product of two delta functions. For example, in 1d, where \(\varphi_i\) are the usual piecewise linear "hat functions", we would get integrals of the sort
\begin{align*} \int_0^L (\Delta \varphi_i) (\Delta \varphi_j) = \int_0^L \frac 1h \left[\delta(x-x_{i-1}) - 2\delta(x-x_i) + \delta(x-x_{i+1})\right] \frac 1h \left[\delta(x-x_{j-1}) - 2\delta(x-x_j) + \delta(x-x_{j+1})\right] \end{align*}
where \(x_i\) is the node location at which the shape function \(\varphi_i\) is defined, and \(h\) is the mesh size (assumed uniform). The problem is that delta functions in integrals are defined using the relationship
\begin{align*} \int_0^L \delta(x-\hat x) f(x) \; dx = f(\hat x). \end{align*}
But that only works if (i) \(f(\cdot)\) is actually well defined at \(\hat x\), and (ii) if it is finite. On the other hand, an integral of the form
\begin{align*} \int_0^L \delta(x-x_i) \delta (x-x_i) \end{align*}
does not make sense. Similar reasoning can be applied for 2d and 3d situations.
In other words: This approach of trying to integrate over the entire domain and then integrating by parts can't work.
Historically, numerical analysts have tried to address this by inventing finite elements that are "C<sup>1</sup> continuous", i.e., that use shape functions that are not just continuous but also have continuous first derivatives. This is the realm of elements such as the Argyris element, the Clough-Tocher element and others, all developed in the late 1960s. From a twenty-first century perspective, they can only be described as bizarre in their construction. They are also exceedingly cumbersome to implement if one wants to use general meshes. As a consequence, they have largely fallen out of favor and deal.II currently does not contain implementations of these shape functions.
So how does one approach solving such problems then? That depends a bit on the boundary conditions. If one has the first set of boundary conditions, i.e., if the equation is
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \Delta u(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \end{align*}
then the following trick works (at least if the domain is convex, see below): In the same way as we obtained the mixed Laplace equation of step-20 from the regular Laplace equation by introducing a second variable, we can here introduce a variable \(v=\Delta u\) and can then replace the equations above by the following, "mixed" system:
\begin{align*} -\Delta u(\mathbf x) +v(\mathbf x) &= 0 \qquad \qquad &&\forall \mathbf x \in \Omega, \\ -\Delta v(\mathbf x) &= -f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ v(\mathbf x) &= h(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
In other words, we end up with what is in essence a system of two coupled Laplace equations for \(u,v\), each with Dirichlet-type boundary conditions. We know how to solve such problems, and it should not be very difficult to construct good solvers and preconditioners for this system either using the techniques of step-20 or step-22. So this case is pretty simple to deal with.
The more complicated situation is if we have the "clamped" boundary conditions, i.e., if the equation looks like this:
\begin{align*} \Delta^2 u(\mathbf x) &= f(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \Omega, \\ u(\mathbf x) &= g(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega, \\ \frac{\partial u(\mathbf x)}{\partial \mathbf n} &= j(\mathbf x) \qquad \qquad &&\forall \mathbf x \in \partial\Omega. \end{align*}
The same trick with the mixed system does not work here, because we would end up with both Dirichlet and Neumann boundary conditions for \(u\), but none for \(v\).
The solution to this conundrum arrived with the Discontinuous Galerkin method wave in the 1990s and early 2000s: In much the same way as one can use discontinuous shape functions for the Laplace equation by penalizing the size of the discontinuity to obtain a scheme for an equation that has one derivative on each shape function, we can use a scheme that uses continuous (but not \(C^1\) continuous) shape functions and penalize the jump in the derivative to obtain a scheme for an equation that has two derivatives on each shape function. In analogy to the Interior Penalty (IP) method for the Laplace equation, this scheme for the biharmonic equation is typically called the \(C^0\) IP (or C0IP) method, since it uses \(C^0\) (continuous but not continuously differentiable) shape functions with an interior penalty formulation.
We base this program on the \(C^0\) IP method presented by Susanne Brenner and Li-Yeng Sung in the paper "C \_form#5084 Interior Penalty Method for Linear Fourth Order Boundary Value Problems on polygonal domains'' @cite Brenner2005 , where the method is derived for the biharmonic equation with "clamped" boundary conditions. As mentioned, this method relies on the use of \_form#1029 Lagrange finite elements where the \_form#442 continuity requirement is relaxed and has been replaced with interior penalty techniques. To derive this method, we consider a \_form#1029 shape function \_form#2525 which vanishes on \_form#2935. We introduce notation \_form#5085 as the set of all faces of \_form#5086, \_form#5087 as the set of boundary faces, and \_form#5088 as the set of interior faces for use further down below. Since the higher order derivatives of \_form#2525 have two values on each interface \_form#5089 (shared by the two cells \_form#5090), we cope with this discontinuity by defining the following single-valued functions on \_form#2082: \_form#5091@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl for \_form#5092 (i.e., for the gradient and the matrix of second derivatives), and where \_form#2764 denotes a unit vector normal to \_form#2082 pointing from \_form#5093 to \_form#5094. In the literature, these functions are referred to as the "jump" and "average" operations, respectively. To obtain the \_form#1029 IP approximation \_form#1172, we left multiply the biharmonic equation by \_form#2525, and then integrate over \_form#264. As explained above, we can't do the integration by parts on all of \_form#264 with these shape functions, but we can do it on each cell individually since the shape functions are just polynomials on each cell. Consequently, we start by using the following integration-by-parts formula on each mesh cell \_form#5095: \_form#5096@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl At this point, we have two options: We can integrate the domain term's \_form#5097 one more time to obtain \_form#5098@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl For a variety of reasons, this turns out to be a variation that is not useful for our purposes. Instead, what we do is recognize that \_form#5099, and we can re-sort these operations as \_form#5100 where we typically write \_form#5101 to indicate that this is the "Hessian" matrix of second derivatives. With this re-ordering, we can now integrate the divergence, rather than the gradient operator, and we get the following instead: \_form#5102@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl Here, the colon indicates a double-contraction over the indices of the matrices to its left and right, i.e., the scalar product between two tensors. The outer product of two vectors \_form#5103 yields the matrix \_form#5104. Then, we sum over all cells \_form#5105, and take into account that this means that every interior face appears twice in the sum. If we therefore split everything into a sum of integrals over cell interiors and a separate sum over cell interfaces, we can use the jump and average operators defined above. There are two steps left: First, because our shape functions are continuous, the gradients of the shape functions may be discontinuous, but the continuity guarantees that really only the normal component of the gradient is discontinuous across faces whereas the tangential component(s) are continuous. Second, the discrete formulation that results is not stable as the mesh size goes to zero, and to obtain a stable formulation that converges to the correct solution, we need to add the following terms: \_form#5106@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl Then, after making cancellations that arise, we arrive at the following C0IP formulation of the biharmonic equation: find \_form#1172 such that \_form#5107@_fakenl on \_form#2651 and \_form#5108@_fakenl@_fakenl where \_form#5109@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl and \_form#5110@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl Here, \_form#1463 is the penalty parameter which both weakly enforces the boundary condition \_form#5111@_fakenl@_fakenl on the boundary interfaces \_form#5112, and also ensures that in the limit \_form#4069, \_form#1172 converges to a \_form#442 continuous function. \_form#1463 is chosen to be large enough to guarantee the stability of the method. We will discuss our choice in the program below. @anchor ConvergenceRates <a></a><h4>Convergence Rates </h4> On polygonal domains, the weak solution \_form#314 to the biharmonic equation lives in \_form#5113 where \_form#5114 is determined by the interior angles at the corners of \_form#264. For instance, whenever \_form#264 is convex, \_form#2647; \_form#143 may be less than one if the domain has re-entrant corners but \_form#143 is close to \_form#1261 if one of all interior angles is close to \_form#1271. Now suppose that the \_form#1029 IP solution \_form#1172 is approximated by \_form#1029 shape functions with polynomial degree \_form#5115. Then the discretization outlined above yields the convergence rates as discussed below. <b>Convergence in the \_form#1029 IP-norm</b> Ideally, we would like to measure convergence in the "energy norm" \_form#5116. However, this does not work because, again, the discrete solution \_form#1172 does not have two (weak) derivatives. Instead, one can define a discrete ( \_form#1029 IP) seminorm that is "equivalent" to the energy norm, as follows: \_form#5117@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl In this seminorm, the theory in the paper mentioned above yields that we can expect \_form#5118@_fakenl@_fakenl much as one would expect given the convergence rates we know are true for the usual discretizations of the Laplace equation. Of course, this is true only if the exact solution is sufficiently smooth. Indeed, if \_form#5119 with \_form#5120, \_form#5121 where \_form#5122, then the convergence rate of the \_form#1029 IP method is \_form#5123. In other words, the optimal convergence rate can only be expected if the solution is so smooth that \_form#5124; this can only happen if (i) the domain is convex with a sufficiently smooth boundary, and (ii) \_form#5125. In practice, of course, the solution is what it is (independent of the polynomial degree we choose), and the last condition can then equivalently be read as saying that there is definitely no point in choosing \_form#158 large if \_form#311 is not also large. In other words, the only reasonably choices for \_form#158 are \_form#5126@_fakenl because larger polynomial degrees do not result in higher convergence orders. For the purposes of this program, we're a bit too lazy to actually implement this equivalent seminorm – though it's not very difficult and would make for a good exercise. Instead, we'll simply check in the program what the "broken" \_form#2546 seminorm \_form#5127@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl yields. The convergence rate in this norm can, from a theoretical perspective, of course not be <i>worse</i> than the one for \_form#5128 because it contains only a subset of the necessary terms, but it could at least conceivably be better. It could also be the case that we get the optimal convergence rate even though there is a bug in the program, and that that bug would only show up in sub-optimal rates for the additional terms present in \_form#5128. But, one might hope that if we get the optimal rate in the broken norm and the norms discussed below, then the program is indeed correct. The results section will demonstrate that we obtain optimal rates in all norms shown. <b>Convergence in the \_form#993-norm</b> The optimal convergence rate in the \_form#993-norm is \_form#5129 provided \_form#5130. More details can be found in Theorem 4.6 of @cite Engel2002 . The default in the program below is to choose \_form#5131. In that case, the theorem does not apply, and indeed one only gets \_form#5132 instead of \_form#5133 as we will show in the results section. <b>Convergence in the \_form#69-seminorm</b> Given that we expect \_form#5134 in the best of cases for a norm equivalent to the \_form#2546 seminorm, and \_form#5129 for the \_form#993 norm, one may ask about what happens in the \_form#69 seminorm that is intermediate to the two others. A reasonable guess is that one should expect \_form#5135. There is probably a paper somewhere that proves this, but we also verify that this conjecture is experimentally true below. @anchor OtherBoundaryConditions <a></a><h3>Other Boundary Conditions</h3> We remark that the derivation of the \_form#1029 IP method for the biharmonic equation with other boundary conditions – for instance, for the first set of boundary conditions namely \_form#5136@_fakenl and \_form#5137 on \_form#2935 – can be obtained with suitable modifications to \_form#5138 and \_form#5139 described in the book chapter @cite Brenner2011 . @anchor Thetestcase <a></a><h3>The testcase</h3> The last step that remains to describe is what this program solves for. As always, a trigonometric function is both a good and a bad choice because it does not lie in any polynomial space in which we may seek the solution while at the same time being smoother than real solutions typically are (here, it is in \_form#3439 while real solutions are typically only in \_form#5140 or so on convex polygonal domains, or somewhere between \_form#2546 and \_form#5140 if the domain is not convex). But, since we don't have the means to describe solutions of realistic problems in terms of relatively simple formulas, we just go with the following, on the unit square for the domain \_form#264: \_form#5141@_fakenl@_fakenl As a consequence, we then need choose as boundary conditions the following: \_form#5142@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl The right hand side is easily computed as \_form#5143@_fakenl@_fakenl The program has classes <tt>ExactSolution::Solution</tt> and <tt>ExactSolution::RightHandSide</tt> that encode this information. @anchor CommProg <a></a> <h1> The commented program</h1> @anchor Includefiles <a></a> <h3>Include files</h3> The first few include files have already been used in the previous example, so we will not explain their meaning here again. The principal structure of the program is very similar to that of, for example, @ref step_4 "step-4" and so we include many of the same header files. @code #include <deal.II/base/quadrature_lib.h> #include <deal.II/base/function.h> #include <deal.II/lac/vector.h> #include <deal.II/lac/full_matrix.h> #include <deal.II/lac/dynamic_sparsity_pattern.h> #include <deal.II/lac/solver_cg.h> #include <deal.II/lac/sparse_matrix.h> #include <deal.II/lac/sparse_direct.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_values.h> #include <deal.II/fe/mapping_q.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_tools.h> #include <deal.II/numerics/vector_tools.h> #include <deal.II/numerics/data_out.h> @endcode The two most interesting header files will be these two: @code #include <deal.II/fe/fe_interface_values.h> #include <deal.II/meshworker/mesh_loop.h> @endcode The first of these is responsible for providing the class FEInterfaceValues that can be used to evaluate quantities such as the jump or average of shape functions (or their gradients) across interfaces between cells. This class will be quite useful in evaluating the penalty terms that appear in the C0IP formulation. @code #include <fstream> #include <iostream> #include <cmath> namespace Step47 { using namespace dealii; @endcode In the following namespace, let us define the exact solution against which we will compare the numerically computed one. It has the form \_form#5144 (only the 2d case is implemented), and the namespace also contains a class that corresponds to the right hand side that produces this solution. @code namespace ExactSolution { using numbers::PI; template <int dim> class Solution : public Function<dim> { public: static_assert(dim == 2, "Only dim==2 is implemented."); virtual double value(const Point<dim> &p, const unsigned int /*component*/ = 0) const override { return std::sin(PI * p[0]) * std::sin(PI * p[1]); } virtual Tensor<1, dim> gradient(const Point<dim> &p, const unsigned int /*component*/ = 0) const override { Tensor<1, dim> r; r[0] = PI * std::cos(PI * p[0]) * std::sin(PI * p[1]); r[1] = PI * std::cos(PI * p[1]) * std::sin(PI * p[0]); return r; } virtual void hessian_list(const std::vector<Point<dim>> & points, std::vector<SymmetricTensor<2, dim>> &hessians, const unsigned int /*component*/ = 0) const override { for (unsigned i = 0; i < points.size(); ++i) { const double x = points[i][0]; const double y = points[i][1]; hessians[i][0][0] = -PI * PI * std::sin(PI * x) * std::sin(PI * y); hessians[i][0][1] = PI * PI * std::cos(PI * x) * std::cos(PI * y); hessians[i][1][1] = -PI * PI * std::sin(PI * x) * std::sin(PI * y); } } }; template <int dim> class RightHandSide : public Function<dim> { public: static_assert(dim == 2, "Only dim==2 is implemented"); virtual double value(const Point<dim> &p, const unsigned int /*component*/ = 0) const override { return 4 * std::pow(PI, 4.0) * std::sin(PI * p[0]) * std::sin(PI * p[1]); } }; } // namespace ExactSolution @endcode @anchor Themainclass <a></a> <h3>The main class</h3> The following is the principal class of this tutorial program. It has the structure of many of the other tutorial programs and there should really be nothing particularly surprising about its contents or the constructor that follows it. @code template <int dim> class BiharmonicProblem { public: BiharmonicProblem(const unsigned int fe_degree); void run(); private: void make_grid(); void setup_system(); void assemble_system(); void solve(); void compute_errors(); void output_results(const unsigned int iteration) const; Triangulation<dim> triangulation; MappingQ<dim> mapping; FE_Q<dim> fe; DoFHandler<dim> dof_handler; AffineConstraints<double> constraints; SparsityPattern sparsity_pattern; SparseMatrix<double> system_matrix; Vector<double> solution; Vector<double> system_rhs; }; template <int dim> BiharmonicProblem<dim>::BiharmonicProblem(const unsigned int fe_degree) : mapping(1) , fe(fe_degree) , dof_handler(triangulation) {} @endcode Next up are the functions that create the initial mesh (a once refined unit square) and set up the constraints, vectors, and matrices on each mesh. Again, both of these are essentially unchanged from many previous tutorial programs. @code template <int dim> void BiharmonicProblem<dim>::make_grid() { GridGenerator::hyper_cube(triangulation, 0., 1.); triangulation.refine_global(1); std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl << "Total number of cells: " << triangulation.n_cells() << std::endl; } template <int dim> void BiharmonicProblem<dim>::setup_system() { dof_handler.distribute_dofs(fe); std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; constraints.clear(); DoFTools::make_hanging_node_constraints(dof_handler, constraints); VectorTools::interpolate_boundary_values(dof_handler, 0, ExactSolution::Solution<dim>(), constraints); constraints.close(); DynamicSparsityPattern dsp(dof_handler.n_dofs()); DoFTools::make_flux_sparsity_pattern(dof_handler, dsp, constraints, true); sparsity_pattern.copy_from(dsp); system_matrix.reinit(sparsity_pattern); solution.reinit(dof_handler.n_dofs()); system_rhs.reinit(dof_handler.n_dofs()); } @endcode @anchor Assemblingthelinearsystem <a></a> <h4>Assembling the linear system</h4> The following pieces of code are more interesting. They all relate to the assembly of the linear system. While assembling the cell-interior terms is not of great difficulty – that works in essence like the assembly of the corresponding terms of the Laplace equation, and you have seen how this works in @ref step_4 "step-4" or @ref step_6 "step-6", for example – the difficulty is with the penalty terms in the formulation. These require the evaluation of gradients of shape functions at interfaces of cells. At the least, one would therefore need to use two FEFaceValues objects, but if one of the two sides is adaptively refined, then one actually needs an FEFaceValues and one FESubfaceValues objects; one also needs to keep track which shape functions live where, and finally we need to ensure that every face is visited only once. All of this is a substantial overhead to the logic we really want to implement (namely the penalty terms in the bilinear form). As a consequence, we will make use of the FEInterfaceValues class – a helper class in deal.II that allows us to abstract away the two FEFaceValues or FESubfaceValues objects and directly access what we really care about: jumps, averages, etc. But this doesn't yet solve our problem of having to keep track of which faces we have already visited when we loop over all cells and all of their faces. To make this process simpler, we use the MeshWorker::mesh_loop() function that provides a simple interface for this task: Based on the ideas outlined in the WorkStream namespace documentation, MeshWorker::mesh_loop() requires three functions that do work on cells, interior faces, and boundary faces. These functions work on scratch objects for intermediate results, and then copy the result of their computations into copy data objects from where a copier function copies them into the global matrix and right hand side objects. The following structures then provide the scratch and copy objects that are necessary for this approach. You may look up the WorkStream namespace as well as the @ref threads "Parallel computing with multiple processors" module for more information on how they typically work. @code template <int dim> struct ScratchData { ScratchData(const Mapping<dim> & mapping, const FiniteElement<dim> &fe, const unsigned int quadrature_degree, const UpdateFlags update_flags, const UpdateFlags interface_update_flags) : fe_values(mapping, fe, QGauss<dim>(quadrature_degree), update_flags) , fe_interface_values(mapping, fe, QGauss<dim - 1>(quadrature_degree), interface_update_flags) {} ScratchData(const ScratchData<dim> &scratch_data) : fe_values(scratch_data.fe_values.get_mapping(), scratch_data.fe_values.get_fe(), scratch_data.fe_values.get_quadrature(), scratch_data.fe_values.get_update_flags()) , fe_interface_values(scratch_data.fe_values.get_mapping(), scratch_data.fe_values.get_fe(), scratch_data.fe_interface_values.get_quadrature(), scratch_data.fe_interface_values.get_update_flags()) {} FEValues<dim> fe_values; FEInterfaceValues<dim> fe_interface_values; }; struct CopyData { CopyData(const unsigned int dofs_per_cell) : cell_matrix(dofs_per_cell, dofs_per_cell) , cell_rhs(dofs_per_cell) , local_dof_indices(dofs_per_cell) {} CopyData(const CopyData &) = default; CopyData(CopyData &&) = default; ~CopyData() = default; CopyData &operator=(const CopyData &) = default; CopyData &operator=(CopyData &&) = default; struct FaceData { FullMatrix<double> cell_matrix; std::vector<types::global_dof_index> joint_dof_indices; }; FullMatrix<double> cell_matrix; Vector<double> cell_rhs; std::vector<types::global_dof_index> local_dof_indices; std::vector<FaceData> face_data; }; @endcode The more interesting part is where we actually assemble the linear system. Fundamentally, this function has five parts: - The definition of the <tt>cell_worker</tt> lambda function, a small function that is defined within the <tt>assemble_system()</tt> function and that will be responsible for computing the local integrals on an individual cell. It will work on a copy of the <tt>ScratchData</tt> class and put its results into the corresponding <tt>CopyData</tt> object. - The definition of the <tt>face_worker</tt> lambda function that does the integration of all terms that live on the interfaces between cells. - The definition of the <tt>boundary_worker</tt> function that does the same but for cell faces located on the boundary of the domain. - The definition of the <tt>copier</tt> function that is responsible for copying all of the data the previous three functions have put into copy objects for a single cell, into the global matrix and right hand side. The fifth part is the one where we bring all of this together. Let us go through each of these pieces necessary for the assembly in turns. @code template <int dim> void BiharmonicProblem<dim>::assemble_system() { using Iterator = typename DoFHandler<dim>::active_cell_iterator; @endcode The first piece is the <tt>cell_worker</tt> that does the assembly on the cell interiors. It is a (lambda) function that takes a cell (input), a scratch object, and a copy object (output) as arguments. It looks like the assembly functions of many other of the tutorial programs, or at least the body of the loop over all cells. The terms we integrate here are the cell contribution \_form#5145@_fakenl@_fakenl to the global matrix, and \_form#5146@_fakenl@_fakenl to the right hand side vector. We use the same technique as used in the assembly of @ref step_22 "step-22" to accelerate the function: Instead of calling <tt>fe_values.shape_hessian(i, qpoint)</tt> in the innermost loop, we create a variable <tt>hessian_i</tt> that evaluates this value once in the loop over <tt>i</tt> and re-use the so-evaluated value in the loop over <tt>j</tt>. For symmetry, we do the same with a variable <tt>hessian_j</tt>, although it is indeed only used once and we could have left the call to <tt>fe_values.shape_hessian(j,qpoint)</tt> in the instruction that computes the scalar product between the two terms. @code auto cell_worker = [&](const Iterator & cell, ScratchData<dim> &scratch_data, CopyData & copy_data) { copy_data.cell_matrix = 0; copy_data.cell_rhs = 0; FEValues<dim> &fe_values = scratch_data.fe_values; fe_values.reinit(cell); cell->get_dof_indices(copy_data.local_dof_indices); const ExactSolution::RightHandSide<dim> right_hand_side; const unsigned int dofs_per_cell = scratch_data.fe_values.get_fe().n_dofs_per_cell(); for (unsigned int qpoint = 0; qpoint < fe_values.n_quadrature_points; ++qpoint) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { const Tensor<2, dim> &hessian_i = fe_values.shape_hessian(i, qpoint); for (unsigned int j = 0; j < dofs_per_cell; ++j) { const Tensor<2, dim> &hessian_j = fe_values.shape_hessian(j, qpoint); copy_data.cell_matrix(i, j) += scalar_product(hessian_i, // nabla^2 phi_i(x) hessian_j) * // nabla^2 phi_j(x) fe_values.JxW(qpoint); // dx } copy_data.cell_rhs(i) += fe_values.shape_value(i, qpoint) * // phi_i(x) right_hand_side.value( fe_values.quadrature_point(qpoint)) * // f(x) fe_values.JxW(qpoint); // dx } } }; @endcode The next building block is the one that assembles penalty terms on each of the interior faces of the mesh. As described in the documentation of MeshWorker::mesh_loop(), this function receives arguments that denote a cell and its neighboring cell, as well as (for each of the two cells) the face (and potentially sub-face) we have to integrate over. Again, we also get a scratch object, and a copy object for putting the results in. The function has three parts itself. At the top, we initialize the FEInterfaceValues object and create a new <tt>CopyData::FaceData</tt> object to store our input in. This gets pushed to the end of the <tt>copy_data.face_data</tt> variable. We need to do this because the number of faces (or subfaces) over which we integrate for a given cell differs from cell to cell, and the sizes of these matrices also differ, depending on what degrees of freedom are adjacent to the face or subface. As discussed in the documentation of MeshWorker::mesh_loop(), the copy object is reset every time a new cell is visited, so that what we push to the end of <tt>copy_data.face_data()</tt> is really all that the later <tt>copier</tt> function gets to see when it copies the contributions of each cell to the global matrix and right hand side objects. @code auto face_worker = [&](const Iterator & cell, const unsigned int &f, const unsigned int &sf, const Iterator & ncell, const unsigned int &nf, const unsigned int &nsf, ScratchData<dim> & scratch_data, CopyData & copy_data) { FEInterfaceValues<dim> &fe_interface_values = scratch_data.fe_interface_values; fe_interface_values.reinit(cell, f, sf, ncell, nf, nsf); copy_data.face_data.emplace_back(); CopyData::FaceData ©_data_face = copy_data.face_data.back(); copy_data_face.joint_dof_indices = fe_interface_values.get_interface_dof_indices(); const unsigned int n_interface_dofs = fe_interface_values.n_current_interface_dofs(); copy_data_face.cell_matrix.reinit(n_interface_dofs, n_interface_dofs); @endcode The second part deals with determining what the penalty parameter should be. By looking at the units of the various terms in the bilinear form, it is clear that the penalty has to have the form \_form#5147 (i.e., one over length scale), but it is not a priori obvious how one should choose the dimension-less number \_form#1463. From the discontinuous Galerkin theory for the Laplace equation, one might conjecture that the right choice is \_form#5148 is the right choice, where \_form#158 is the polynomial degree of the finite element used. We will discuss this choice in a bit more detail in the results section of this program. In the formula above, \_form#1417 is the size of cell \_form#87. But this is not quite so straightforward either: If one uses highly stretched cells, then a more involved theory says that \_form#1170 should be replaced by the diameter of cell \_form#87 normal to the direction of the edge in question. It turns out that there is a function in deal.II for that. Secondly, \_form#1417 may be different when viewed from the two different sides of a face. To stay on the safe side, we take the maximum of the two values. We will note that it is possible that this computation has to be further adjusted if one were to use hanging nodes resulting from adaptive mesh refinement. @code const unsigned int p = fe.degree; const double gamma_over_h = std::max((1.0 * p * (p + 1) / cell->extent_in_direction( GeometryInfo<dim>::unit_normal_direction[f])), (1.0 * p * (p + 1) / ncell->extent_in_direction( GeometryInfo<dim>::unit_normal_direction[nf]))); @endcode Finally, and as usual, we loop over the quadrature points and indices <tt>i</tt> and <tt>j</tt> to add up the contributions of this face or sub-face. These are then stored in the <tt>copy_data.face_data</tt> object created above. As for the cell worker, we pull the evaluation of averages and jumps out of the loops if possible, introducing local variables that store these results. The assembly then only needs to use these local variables in the innermost loop. Regarding the concrete formula this code implements, recall that the interface terms of the bilinear form were as follows: \_form#5149@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl @code for (unsigned int qpoint = 0; qpoint < fe_interface_values.n_quadrature_points; ++qpoint) { const auto &n = fe_interface_values.normal(qpoint); for (unsigned int i = 0; i < n_interface_dofs; ++i) { const double av_hessian_i_dot_n_dot_n = (fe_interface_values.average_hessian(i, qpoint) * n * n); const double jump_grad_i_dot_n = (fe_interface_values.jump_gradient(i, qpoint) * n); for (unsigned int j = 0; j < n_interface_dofs; ++j) { const double av_hessian_j_dot_n_dot_n = (fe_interface_values.average_hessian(j, qpoint) * n * n); const double jump_grad_j_dot_n = (fe_interface_values.jump_gradient(j, qpoint) * n); copy_data_face.cell_matrix(i, j) += (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n } * jump_grad_j_dot_n // [grad u n] - av_hessian_j_dot_n_dot_n // - {grad^2 u n n } * jump_grad_i_dot_n // [grad v n] + // + gamma_over_h * // gamma/h jump_grad_i_dot_n * // [grad v n] jump_grad_j_dot_n) * // [grad u n] fe_interface_values.JxW(qpoint); // dx } } } }; @endcode The third piece is to do the same kind of assembly for faces that are at the boundary. The idea is the same as above, of course, with only the difference that there are now penalty terms that also go into the right hand side. As before, the first part of the function simply sets up some helper objects: @code auto boundary_worker = [&](const Iterator & cell, const unsigned int &face_no, ScratchData<dim> & scratch_data, CopyData & copy_data) { FEInterfaceValues<dim> &fe_interface_values = scratch_data.fe_interface_values; fe_interface_values.reinit(cell, face_no); const auto &q_points = fe_interface_values.get_quadrature_points(); copy_data.face_data.emplace_back(); CopyData::FaceData ©_data_face = copy_data.face_data.back(); const unsigned int n_dofs = fe_interface_values.n_current_interface_dofs(); copy_data_face.joint_dof_indices = fe_interface_values.get_interface_dof_indices(); copy_data_face.cell_matrix.reinit(n_dofs, n_dofs); const std::vector<double> &JxW = fe_interface_values.get_JxW_values(); const std::vector<Tensor<1, dim>> &normals = fe_interface_values.get_normal_vectors(); const ExactSolution::Solution<dim> exact_solution; std::vector<Tensor<1, dim>> exact_gradients(q_points.size()); exact_solution.gradient_list(q_points, exact_gradients); @endcode Positively, because we now only deal with one cell adjacent to the face (as we are on the boundary), the computation of the penalty factor \_form#1463 is substantially simpler: @code const unsigned int p = fe.degree; const double gamma_over_h = (1.0 * p * (p + 1) / cell->extent_in_direction( GeometryInfo<dim>::unit_normal_direction[face_no])); @endcode The third piece is the assembly of terms. This is now slightly more involved since these contains both terms for the matrix and for the right hand side. The former is exactly the same as for the interior faces stated above if one just defines the jump and average appropriately (which is what the FEInterfaceValues class does). The latter requires us to evaluate the boundary conditions \_form#5065, which in the current case (where we know the exact solution) we compute from \_form#5150@_fakenl. The term to be added to the right hand side vector is then \_form#5151@_fakenl. @code for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) { const auto &n = normals[qpoint]; for (unsigned int i = 0; i < n_dofs; ++i) { const double av_hessian_i_dot_n_dot_n = (fe_interface_values.average_hessian(i, qpoint) * n * n); const double jump_grad_i_dot_n = (fe_interface_values.jump_gradient(i, qpoint) * n); for (unsigned int j = 0; j < n_dofs; ++j) { const double av_hessian_j_dot_n_dot_n = (fe_interface_values.average_hessian(j, qpoint) * n * n); const double jump_grad_j_dot_n = (fe_interface_values.jump_gradient(j, qpoint) * n); copy_data_face.cell_matrix(i, j) += (-av_hessian_i_dot_n_dot_n // - {grad^2 v n n} * jump_grad_j_dot_n // [grad u n] - av_hessian_j_dot_n_dot_n // - {grad^2 u n n} * jump_grad_i_dot_n // [grad v n] + gamma_over_h // gamma/h * jump_grad_i_dot_n // [grad v n] * jump_grad_j_dot_n // [grad u n] ) * JxW[qpoint]; // dx } copy_data.cell_rhs(i) += (-av_hessian_i_dot_n_dot_n * // - {grad^2 v n n } (exact_gradients[qpoint] * n) // (grad u_exact . n) + // + gamma_over_h // gamma/h * jump_grad_i_dot_n // [grad v n] * (exact_gradients[qpoint] * n) // (grad u_exact . n) ) * JxW[qpoint]; // dx } } }; @endcode Part 4 is a small function that copies the data produced by the cell, interior, and boundary face assemblers above into the global matrix and right hand side vector. There really is not very much to do here: We distribute the cell matrix and right hand side contributions as we have done in almost all of the other tutorial programs using the constraints objects. We then also have to do the same for the face matrix contributions that have gained content for the faces (interior and boundary) and that the <tt>face_worker</tt> and <tt>boundary_worker</tt> have added to the <tt>copy_data.face_data</tt> array. @code auto copier = [&](const CopyData ©_data) { constraints.distribute_local_to_global(copy_data.cell_matrix, copy_data.cell_rhs, copy_data.local_dof_indices, system_matrix, system_rhs); for (auto &cdf : copy_data.face_data) { constraints.distribute_local_to_global(cdf.cell_matrix, cdf.joint_dof_indices, system_matrix); } }; @endcode Having set all of this up, what remains is to just create a scratch and copy data object and call the MeshWorker::mesh_loop() function that then goes over all cells and faces, calls the respective workers on them, and then the copier function that puts things into the global matrix and right hand side. As an additional benefit, MeshWorker::mesh_loop() does all of this in parallel, using as many processor cores as your machine happens to have. @code const unsigned int n_gauss_points = dof_handler.get_fe().degree + 1; ScratchData<dim> scratch_data(mapping, fe, n_gauss_points, update_values | update_gradients | update_hessians | update_quadrature_points | update_JxW_values, update_values | update_gradients | update_hessians | update_quadrature_points | update_JxW_values | update_normal_vectors); CopyData copy_data(dof_handler.get_fe().n_dofs_per_cell()); MeshWorker::mesh_loop(dof_handler.begin_active(), dof_handler.end(), cell_worker, copier, scratch_data, copy_data, MeshWorker::assemble_own_cells | MeshWorker::assemble_boundary_faces | MeshWorker::assemble_own_interior_faces_once, boundary_worker, face_worker); } @endcode @anchor Solvingthelinearsystemandpostprocessing <a></a> <h4>Solving the linear system and postprocessing</h4> The show is essentially over at this point: The remaining functions are not overly interesting or novel. The first one simply uses a direct solver to solve the linear system (see also @ref step_29 "step-29"): @code template <int dim> void BiharmonicProblem<dim>::solve() { std::cout << " Solving system..." << std::endl; SparseDirectUMFPACK A_direct; A_direct.initialize(system_matrix); A_direct.vmult(solution, system_rhs); constraints.distribute(solution); } @endcode The next function evaluates the error between the computed solution and the exact solution (which is known here because we have chosen the right hand side and boundary values in a way so that we know the corresponding solution). In the first two code blocks below, we compute the error in the \_form#993 norm and the \_form#69 semi-norm. @code template <int dim> void BiharmonicProblem<dim>::compute_errors() { { Vector<float> norm_per_cell(triangulation.n_active_cells()); VectorTools::integrate_difference(mapping, dof_handler, solution, ExactSolution::Solution<dim>(), norm_per_cell, QGauss<dim>(fe.degree + 2), VectorTools::L2_norm); const double error_norm = VectorTools::compute_global_error(triangulation, norm_per_cell, VectorTools::L2_norm); std::cout << " Error in the L2 norm : " << error_norm << std::endl; } { Vector<float> norm_per_cell(triangulation.n_active_cells()); VectorTools::integrate_difference(mapping, dof_handler, solution, ExactSolution::Solution<dim>(), norm_per_cell, QGauss<dim>(fe.degree + 2), VectorTools::H1_seminorm); const double error_norm = VectorTools::compute_global_error(triangulation, norm_per_cell, VectorTools::H1_seminorm); std::cout << " Error in the H1 seminorm : " << error_norm << std::endl; } @endcode Now also compute an approximation to the \_form#2546 seminorm error. The actual \_form#2546 seminorm would require us to integrate second derivatives of the solution \_form#1172, but given the Lagrange shape functions we use, \_form#1172 of course has kinks at the interfaces between cells, and consequently second derivatives are singular at interfaces. As a consequence, we really only integrate over the interior of cells and ignore the interface contributions. This is <em>not</em> an equivalent norm to the energy norm for the problem, but still gives us an idea of how fast the error converges. We note that one could address this issue by defining a norm that is equivalent to the energy norm. This would involve adding up not only the integrals over cell interiors as we do below, but also adding penalty terms for the jump of the derivative of \_form#1172 across interfaces, with an appropriate scaling of the two kinds of terms. We will leave this for later work. @code { const QGauss<dim> quadrature_formula(fe.degree + 2); ExactSolution::Solution<dim> exact_solution; Vector<double> error_per_cell(triangulation.n_active_cells()); FEValues<dim> fe_values(mapping, fe, quadrature_formula, update_values | update_hessians | update_quadrature_points | update_JxW_values); FEValuesExtractors::Scalar scalar(0); const unsigned int n_q_points = quadrature_formula.size(); std::vector<SymmetricTensor<2, dim>> exact_hessians(n_q_points); std::vector<Tensor<2, dim>> hessians(n_q_points); for (auto &cell : dof_handler.active_cell_iterators()) { fe_values.reinit(cell); fe_values[scalar].get_function_hessians(solution, hessians); exact_solution.hessian_list(fe_values.get_quadrature_points(), exact_hessians); double local_error = 0; for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { local_error += ((exact_hessians[q_point] - hessians[q_point]).norm_square() * fe_values.JxW(q_point)); } error_per_cell[cell->active_cell_index()] = std::sqrt(local_error); } const double error_norm = error_per_cell.l2_norm(); std::cout << " Error in the broken H2 seminorm: " << error_norm << std::endl; } } @endcode Equally uninteresting is the function that generates graphical output. It looks exactly like the one in @ref step_6 "step-6", for example. @code template <int dim> void BiharmonicProblem<dim>::output_results(const unsigned int iteration) const { std::cout << " Writing graphical output..." << std::endl; DataOut<dim> data_out; data_out.attach_dof_handler(dof_handler); data_out.add_data_vector(solution, "solution"); data_out.build_patches(); const std::string filename = ("output_" + Utilities::int_to_string(iteration, 6) + ".vtu"); std::ofstream output_vtu(filename); data_out.write_vtu(output_vtu); } @endcode The same is true for the <tt>run()</tt> function: Just like in previous programs. @code template <int dim> void BiharmonicProblem<dim>::run() { make_grid(); const unsigned int n_cycles = 4; for (unsigned int cycle = 0; cycle < n_cycles; ++cycle) { std::cout << "Cycle " << cycle << " of " << n_cycles << std::endl; triangulation.refine_global(1); setup_system(); assemble_system(); solve(); output_results(cycle); compute_errors(); std::cout << std::endl; } } } // namespace Step47 @endcode @anchor Themainfunction <a></a> <h3>The main() function</h3> Finally for the <tt>main()</tt> function. There is, again, not very much to see here: It looks like the ones in previous tutorial programs. There is a variable that allows selecting the polynomial degree of the element we want to use for solving the equation. Because the C0IP formulation we use requires the element degree to be at least two, we check with an assertion that whatever one sets for the polynomial degree actually makes sense. @code int main() { try { using namespace dealii; using namespace Step47; const unsigned int fe_degree = 2; Assert(fe_degree >= 2, ExcMessage("The C0IP formulation for the biharmonic problem " "only works if one uses elements of polynomial " "degree at least 2.")); BiharmonicProblem<2> biharmonic_problem(fe_degree); biharmonic_problem.run(); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; } @endcode @anchor Results <a></a><h1>Results</h1> We run the program with right hand side and boundary values as discussed in the introduction. These will produce the solution \_form#5152 on the domain \_form#5153. We test this setup using \_form#157, \_form#892, and \_form#1092 elements, which one can change via the <tt>fe_degree</tt> variable in the <tt>main()</tt> function. With mesh refinement, the \_form#993 convergence rates, \_form#69-seminorm rate, and \_form#2546-seminorm convergence of \_form#314 should then be around 2, 2, 1 for \_form#157 (with the \_form#993 norm sub-optimal as discussed in the introduction); 4, 3, 2 for \_form#892; and 5, 4, 3 for \_form#1092, respectively. From the literature, it is not immediately clear what the penalty parameter \_form#1463 should be. For example, @cite Brenner2009 state that it needs to be larger than one, and choose \_form#5154. The FEniCS/Dolphin tutorial chooses it as \_form#5155, see https://fenicsproject.org/docs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html . @cite Wells2007 uses a value for \_form#1463 larger than the number of edges belonging to an element for Kirchhoff plates (see their Section 4.2). This suggests that maybe \_form#5156, \_form#1266, are too small; on the other hand, a value \_form#5157 would be reasonable, where \_form#158 is the degree of polynomials. The last of these choices is the one one would expect to work by comparing to the discontinuous Galerkin formulations for the Laplace equation (see, for example, the discussions in @ref step_39 "step-39" and @ref step_74 "step-74"), and it will turn out to also work here. But we should check what value of \_form#1463 is right, and we will do so below; changing \_form#1463 is easy in the two <tt>face_worker</tt> and <tt>boundary_worker</tt> functions defined in <tt>assemble_system()</tt>. @anchor TestresultsoniQsub2subiwithigammapp1i <a></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = p(p+1)</i> </h3> We run the code with differently refined meshes and get the following convergence rates. <table align="center" class="doxtable"> <tr> <th>Number of refinements </th><th> \_form#5158 </th><th> Conv. rates </th><th> \_form#5159 </th><th> Conv. rates </th><th> \_form#5160 </th><th> Conv. rates </th> </tr> <tr> <td> 2 </td><td> 8.780e-03 </td><td> </td><td> 7.095e-02 </td><td> </td><td> 1.645 </td><td> </td> </tr> <tr> <td> 3 </td><td> 3.515e-03 </td><td> 1.32 </td><td> 2.174e-02 </td><td> 1.70 </td><td> 8.121e-01 </td><td> 1.018 </td> </tr> <tr> <td> 4 </td><td> 1.103e-03 </td><td> 1.67 </td><td> 6.106e-03 </td><td> 1.83 </td><td> 4.015e-01 </td><td> 1.016 </td> </tr> <tr> <td> 5 </td><td> 3.084e-04 </td><td> 1.83 </td><td> 1.622e-03 </td><td> 1.91 </td><td> 1.993e-01 </td><td> 1.010 </td> </tr> </table> We can see that the \_form#993 convergence rates are around 2, \_form#69-seminorm convergence rates are around 2, and \_form#2546-seminorm convergence rates are around 1. The latter two match the theoretically expected rates; for the former, we have no theorem but are not surprised that it is sub-optimal given the remark in the introduction. @anchor TestresultsoniQsub3subiwithigammapp1i <a></a><h3>Test results on <i>Q<sub>3</sub></i> with <i>γ = p(p+1)</i> </h3> <table align="center" class="doxtable"> <tr> <th>Number of refinements </th><th> \_form#5158 </th><th> Conv. rates </th><th> \_form#5159 </th><th> Conv. rates </th><th> \_form#5160 </th><th> Conv. rates </th> </tr> <tr> <td> 2 </td><td> 2.045e-04 </td><td> </td><td> 4.402e-03 </td><td> </td><td> 1.641e-01 </td><td> </td> </tr> <tr> <td> 3 </td><td> 1.312e-05 </td><td> 3.96 </td><td> 5.537e-04 </td><td> 2.99 </td><td> 4.096e-02 </td><td> 2.00 </td> </tr> <tr> <td> 4 </td><td> 8.239e-07 </td><td> 3.99 </td><td> 6.904e-05 </td><td> 3.00 </td><td> 1.023e-02 </td><td> 2.00 </td> </tr> <tr> <td> 5 </td><td> 5.158e-08 </td><td> 3.99 </td><td> 8.621e-06 </td><td> 3.00 </td><td> 2.558e-03 </td><td> 2.00 </td> </tr> </table> We can see that the \_form#993 convergence rates are around 4, \_form#69-seminorm convergence rates are around 3, and \_form#2546-seminorm convergence rates are around 2. This, of course, matches our theoretical expectations. @anchor TestresultsoniQsub4subiwithigammapp1i <a></a><h3>Test results on <i>Q<sub>4</sub></i> with <i>γ = p(p+1)</i> </h3> <table align="center" class="doxtable"> <tr> <th>Number of refinements </th><th> \_form#5158 </th><th> Conv. rates </th><th> \_form#5159 </th><th> Conv. rates </th><th> \_form#5160 </th><th> Conv. rates </th> </tr> <tr> <td> 2 </td><td> 6.510e-06 </td><td> </td><td> 2.215e-04 </td><td> </td><td> 1.275e-02 </td><td> </td> </tr> <tr> <td> 3 </td><td> 2.679e-07 </td><td> 4.60 </td><td> 1.569e-05 </td><td> 3.81 </td><td> 1.496e-03 </td><td> 3.09 </td> </tr> <tr> <td> 4 </td><td> 9.404e-09 </td><td> 4.83 </td><td> 1.040e-06 </td><td> 3.91 </td><td> 1.774e-04 </td><td> 3.07 </td> </tr> <tr> <td> 5 </td><td> 7.943e-10 </td><td> 3.56 </td><td> 6.693e-08 </td><td> 3.95 </td><td> 2.150e-05 </td><td> 3.04 </td> </tr> </table> We can see that the \_form#993 norm convergence rates are around 5, \_form#69-seminorm convergence rates are around 4, and \_form#2546-seminorm convergence rates are around 3. On the finest mesh, the \_form#993 norm convergence rate is much smaller than our theoretical expectations because the linear solver becomes the limiting factor due to round-off. Of course the \_form#993 error is also very small already in that case. @anchor TestresultsoniQsub2subiwithigamma1i <a></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = 1</i> </h3> For comparison with the results above, let us now also consider the case where we simply choose \_form#5161: <table align="center" class="doxtable"> <tr> <th>Number of refinements </th><th> \_form#5158 </th><th> Conv. rates </th><th> \_form#5159 </th><th> Conv. rates </th><th> \_form#5160 </th><th> Conv. rates </th> </tr> <tr> <td> 2 </td><td> 7.350e-02 </td><td> </td><td> 7.323e-01 </td><td> </td><td> 10.343 </td><td> </td> </tr> <tr> <td> 3 </td><td> 6.798e-03 </td><td> 3.43 </td><td> 1.716e-01 </td><td> 2.09 </td><td>4.836 </td><td> 1.09 </td> </tr> <tr> <td> 4 </td><td> 9.669e-04 </td><td> 2.81 </td><td> 6.436e-02 </td><td> 1.41 </td><td> 3.590 </td><td> 0.430 </td> </tr> <tr> <td> 5 </td><td> 1.755e-04 </td><td> 2.46 </td><td> 2.831e-02 </td><td> 1.18 </td><td>3.144 </td><td> 0.19 </td> </tr> </table> Although \_form#993 norm convergence rates of \_form#314 more or less follows the theoretical expectations, the \_form#69-seminorm and \_form#2546-seminorm do not seem to converge as expected. Comparing results from \_form#5156 and \_form#5162, it is clear that \_form#5162 is a better penalty. Given that \_form#5161 is already too small for \_form#157 elements, it may not be surprising that if one repeated the experiment with a \_form#892 element, the results are even more disappointing: One again only obtains convergence rates of 2, 1, zero – i.e., no better than for the \_form#157 element (although the errors are smaller in magnitude). Maybe surprisingly, however, one obtains more or less the expected convergence orders when using \_form#1092 elements. Regardless, this uncertainty suggests that \_form#5161 is at best a risky choice, and at worst an unreliable one and that we should choose \_form#1463 larger. @anchor TestresultsoniQsub2subiwithigamma2i <a></a><h3>Test results on <i>Q<sub>2</sub></i> with <i>γ = 2</i> </h3> Since \_form#5161 is clearly too small, one might conjecture that \_form#5163 might actually work better. Here is what one obtains in that case: <table align="center" class="doxtable"> <tr> <th>Number of refinements </th><th> \_form#5158 </th><th> Conv. rates </th><th> \_form#5159 </th><th> Conv. rates </th><th> \_form#5160 </th><th> Conv. rates </th> </tr> <tr> <td> 2 </td><td> 4.133e-02 </td><td> </td><td> 2.517e-01 </td><td> </td><td> 3.056 </td><td> </td> </tr> <tr> <td> 3 </td><td> 6.500e-03 </td><td>2.66 </td><td> 5.916e-02 </td><td> 2.08 </td><td>1.444 </td><td> 1.08 </td> </tr> <tr> <td> 4 </td><td> 6.780e-04 </td><td> 3.26 </td><td> 1.203e-02 </td><td> 2.296 </td><td> 6.151e-01 </td><td> 1.231 </td> </tr> <tr> <td> 5 </td><td> 1.622e-04 </td><td> 2.06 </td><td> 2.448e-03 </td><td> 2.297 </td><td> 2.618e-01 </td><td> 1.232 </td> </tr> </table> In this case, the convergence rates more or less follow the theoretical expectations, but, compared to the results from \_form#5164@_fakenl, are more variable. Again, we could repeat this kind of experiment for \_form#892 and \_form#1092 elements. In both cases, we will find that we obtain roughly the expected convergence rates. Of more interest may then be to compare the absolute size of the errors. While in the table above, for the \_form#157 case, the errors on the finest grid are comparable between the \_form#5148 and \_form#5163 case, for \_form#892 the errors are substantially larger for \_form#5163 than for \_form#5148. The same is true for the \_form#1092 case. @anchor Conclusionsforthechoiceofthepenaltyparameter <a></a><h3> Conclusions for the choice of the penalty parameter </h3> The conclusions for which of the "reasonable" choices one should use for the penalty parameter is that \_form#5148 yields the expected results. It is, consequently, what the code uses as currently written. @anchor Possibilitiesforextensions <a></a><h3> Possibilities for extensions </h3> There are a number of obvious extensions to this program that would make sense: - The program uses a square domain and a uniform mesh. Real problems don't come this way, and one should verify convergence also on domains with other shapes and, in particular, curved boundaries. One may also be interested in resolving areas of less regularity by using adaptive mesh refinement. - From a more theoretical perspective, the convergence results above only used the "broken" \_form#2546 seminorm \_form#5165 instead of the "equivalent" norm \_form#5128. This is good enough to convince ourselves that the program isn't fundamentally broken. However, it might be interesting to measure the error in the actual norm for which we have theoretical results. Implementing this addition should not be overly difficult using, for example, the FEInterfaceValues class combined with MeshWorker::mesh_loop() in the same spirit as we used for the assembly of the linear system. @anchor Derivationforthesimplysupportedplates <a></a> <h4> Derivation for the simply supported plates </h4> Similar to the "clamped" boundary condition addressed in the implementation, we will derive the \_form#1029 IP finite element scheme for simply supported plates: \_form#5166@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl We multiply the biharmonic equation by the test function \_form#2525 and integrate over \_form#5167 and get: \_form#5168@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl Summing up over all cells \_form#5105,since normal directions of \_form#5169 are pointing at opposite directions on each interior edge shared by two cells and \_form#5170 on \_form#2651, \_form#5171@_fakenl@_fakenl and by the definition of jump over cell interfaces, \_form#5172@_fakenl@_fakenl We separate interior faces and boundary faces of the domain, \_form#5173@_fakenl@_fakenl@_fakenl where \_form#5174 is the set of interior faces. This leads us to \_form#5175@_fakenl@_fakenl@_fakenl In order to symmetrize and stabilize the discrete problem, we add symmetrization and stabilization term. We finally get the \_form#1029 IP finite element scheme for the biharmonic equation: find \_form#1172 such that \_form#5176 on \_form#2651 and \_form#5177@_fakenl@_fakenl where \_form#5178@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl and \_form#5179@_fakenl@_fakenl@_fakenl@_fakenl@_fakenl The implementation of this boundary case is similar to the "clamped" version except that <tt>boundary_worker</tt> is no longer needed for system assembling and the right hand side is changed according to the formulation. @anchor PlainProg <a></a> <h1> The plain program</h1> @include "step-47.cc"
1.9.4