]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Initial version of step-71 for the biharmonic equation.
authorWolfgang Bangerth <bangerth@colostate.edu>
Sat, 30 Nov 2019 04:47:22 +0000 (21:47 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 14 Jan 2020 00:52:45 +0000 (17:52 -0700)
examples/step-71/CMakeLists.txt [new file with mode: 0644]
examples/step-71/doc/builds-on [new file with mode: 0644]
examples/step-71/doc/intro.dox [new file with mode: 0644]
examples/step-71/doc/kind [new file with mode: 0644]
examples/step-71/doc/results.dox [new file with mode: 0644]
examples/step-71/doc/tooltip [new file with mode: 0644]
examples/step-71/step-71.cc [new file with mode: 0644]

diff --git a/examples/step-71/CMakeLists.txt b/examples/step-71/CMakeLists.txt
new file mode 100644 (file)
index 0000000..4773496
--- /dev/null
@@ -0,0 +1,39 @@
+##
+#  CMake script for the step-71 tutorial program:
+##
+
+# Set the name of the project and target:
+SET(TARGET "step-71")
+
+# Declare all source files the target consists of. Here, this is only
+# the one step-X.cc file, but as you expand your project you may wish
+# to add other source files as well. If your project becomes much larger,
+# you may want to either replace the following statement by something like
+#    FILE(GLOB_RECURSE TARGET_SRC  "source/*.cc")
+#    FILE(GLOB_RECURSE TARGET_INC  "include/*.h")
+#    SET(TARGET_SRC ${TARGET_SRC}  ${TARGET_INC}) 
+# or switch altogether to the large project CMakeLists.txt file discussed
+# in the "CMake in user projects" page accessible from the "User info"
+# page of the documentation.
+SET(TARGET_SRC
+  ${TARGET}.cc
+  )
+
+# Usually, you will not need to modify anything beyond this point...
+
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
+
+FIND_PACKAGE(deal.II 9.2.0 QUIET
+  HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
+  )
+IF(NOT ${deal.II_FOUND})
+  MESSAGE(FATAL_ERROR "\n"
+    "*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
+    "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
+    "or set an environment variable \"DEAL_II_DIR\" that contains this path."
+    )
+ENDIF()
+
+DEAL_II_INITIALIZE_CACHED_VARIABLES()
+PROJECT(${TARGET})
+DEAL_II_INVOKE_AUTOPILOT()
diff --git a/examples/step-71/doc/builds-on b/examples/step-71/doc/builds-on
new file mode 100644 (file)
index 0000000..01f2099
--- /dev/null
@@ -0,0 +1 @@
+step-4 step-12
diff --git a/examples/step-71/doc/intro.dox b/examples/step-71/doc/intro.dox
new file mode 100644 (file)
index 0000000..83ec640
--- /dev/null
@@ -0,0 +1,12 @@
+<br>
+
+<i>
+This program was contributed by Natasha Sharma, Guido Kanschat, Timo
+Heister, Wolfgang Bangerth, and Zhuoran Wang.
+</i>
+
+<a name="Intro"></a>
+<h1>Introduction</h1>
+
+
+<h2>The testcase</h2>
diff --git a/examples/step-71/doc/kind b/examples/step-71/doc/kind
new file mode 100644 (file)
index 0000000..c1d9154
--- /dev/null
@@ -0,0 +1 @@
+techniques
diff --git a/examples/step-71/doc/results.dox b/examples/step-71/doc/results.dox
new file mode 100644 (file)
index 0000000..b5eaba9
--- /dev/null
@@ -0,0 +1,2 @@
+<h1>Results</h1>
+
diff --git a/examples/step-71/doc/tooltip b/examples/step-71/doc/tooltip
new file mode 100644 (file)
index 0000000..5166747
--- /dev/null
@@ -0,0 +1 @@
+Solving the fourth-order biharmonic equation
diff --git a/examples/step-71/step-71.cc b/examples/step-71/step-71.cc
new file mode 100644 (file)
index 0000000..3478fcb
--- /dev/null
@@ -0,0 +1,712 @@
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/function.h>
+
+#include <deal.II/lac/sparse_matrix.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/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/fe/fe_interface_values.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_accessor.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+
+#include <fstream>
+#include <iostream>
+#include <cmath>
+
+
+namespace StepBiharmonic
+{
+  using namespace dealii;
+
+
+  namespace ExactSolution
+  {
+    using numbers::PI;
+
+    /**
+     * An exact solution of the form
+     * $ u(x,y) = \sin(\pi x) \sin(\pi y) $.
+     *
+     * Note that this solution has zero boundary values for the *value*
+     * of the solution, but not for its Laplacian. Consequently, the
+     * boundary contribution to the penalty terms is not zero.
+     */
+    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
+      {
+        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
+      {
+        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
+      {
+        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);
+          }
+      }
+    };
+
+
+    /**
+     * The corresponding right hand side.
+     */
+    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
+
+      {
+        return 4 * std::pow(PI, 4.0) * std::sin(PI * p[0]) *
+               std::sin(PI * p[1]);
+      }
+    };
+  } // namespace ExactSolution
+
+
+
+  /*************************************************************/
+  // @sect3{The main class}
+  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;
+    const MappingQ<dim>       mapping;
+    const 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)
+  {}
+
+
+
+  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 c_sparsity(dof_handler.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof_handler,
+                                         c_sparsity,
+                                         constraints,
+                                         true);
+    sparsity_pattern.copy_from(c_sparsity);
+    system_matrix.reinit(sparsity_pattern);
+
+    solution.reinit(dof_handler.n_dofs());
+    system_rhs.reinit(dof_handler.n_dofs());
+  }
+
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData(const Mapping<dim> &      mapping,
+                const FiniteElement<dim> &fe,
+                const unsigned int        quadrature_degree,
+                const UpdateFlags         update_flags = update_values |
+                                                 update_gradients |
+                                                 update_quadrature_points |
+                                                 update_JxW_values,
+                const UpdateFlags interface_update_flags =
+                  update_values | update_gradients | update_quadrature_points |
+                  update_JxW_values | update_normal_vectors)
+      : 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 CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+  };
+
+
+
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<CopyDataFace>            face_data;
+
+    template <class Iterator>
+    void reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+    }
+  };
+
+
+
+  template <int dim>
+  void BiharmonicProblem<dim>::assemble_system()
+  {
+    typedef decltype(dof_handler.begin_active()) Iterator;
+    const ExactSolution::RightHandSide<dim>      right_hand_side;
+
+    auto cell_worker = [&](const Iterator &  cell,
+                           ScratchData<dim> &scratch_data,
+                           CopyData &        copy_data) {
+      const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+      copy_data.reinit(cell, n_dofs);
+      scratch_data.fe_values.reinit(cell);
+
+      const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+      const FEValues<dim> &      fe_v = scratch_data.fe_values;
+      const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+      // scalar_product(fe.shape_hessian_component(j,k,d),
+      // fe.shape_hessian_component(i,k,d));
+      const double nu = 1.0;
+
+      for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+        {
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            {
+              for (unsigned int j = 0; j < n_dofs; ++j)
+                {
+                  // \int_Z \nu \nabla^2 u \cdot \nabla^2 v \, dx.
+                  copy_data.cell_matrix(i, j) +=
+                    nu *
+                    scalar_product(fe_v.shape_hessian(i, point),
+                                   fe_v.shape_hessian(j, point)) *
+                    JxW[point]; // dx
+                }
+
+
+              copy_data.cell_rhs(i) += fe_v.shape_value(i, point) *
+                                       right_hand_side.value(q_points[point]) *
+                                       JxW[point]; // dx
+            }
+        }
+    };
+
+
+    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_i = scratch_data.fe_interface_values;
+      fe_i.reinit(cell, f, sf, ncell, nf, nsf);
+      const auto &q_points = fe_i.get_quadrature_points();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      const unsigned int n_dofs        = fe_i.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices = fe_i.get_interface_dof_indices();
+
+      copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+      const std::vector<double> &        JxW     = fe_i.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_i.get_normal_vectors();
+
+      // eta = 1/2 + 2C_2
+      // gamma = eta/|e|
+
+      double gamma = 1.0; // TODO:
+
+      {
+        int                degree = fe.tensor_degree();
+        const unsigned int normal1 =
+          GeometryInfo<dim>::unit_normal_direction[f];
+        const unsigned int normal2 =
+          GeometryInfo<dim>::unit_normal_direction[nf];
+        const unsigned int deg1sq =
+          degree * (degree + 1); //(deg1 == 0) ? 1 : deg1 * (deg1+1);
+        const unsigned int deg2sq =
+          degree * (degree + 1); //(deg2 == 0) ? 1 : deg2 * (deg2+1);
+
+        double penalty1 = deg1sq / cell->extent_in_direction(normal1);
+        double penalty2 = deg2sq / ncell->extent_in_direction(normal2);
+        if (cell->has_children() ^ ncell->has_children())
+          {
+            penalty1 *= 8;
+          }
+        gamma = 0.5 * (penalty1 + penalty2);
+      }
+
+
+      for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+        {
+          // \int_F -{grad^2 u n n } [grad v n]
+          //   - {grad^2 v n n } [grad u n]
+          //   +  gamma [grad u n ][grad v n]
+          const auto &n = normals[qpoint];
+
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                Assert((fe_i.average_hessian(i, qpoint) * n * n) ==
+                         contract3(n, fe_i.average_hessian(i, qpoint), n),
+                       ExcInternalError());
+
+                Assert((fe_i.jump_gradient(j, qpoint) * n) ==
+                         (n * fe_i.jump_gradient(j, qpoint)),
+                       ExcInternalError());
+
+
+                copy_data_face.cell_matrix(i, j) +=
+                  (-(fe_i.average_hessian(i, qpoint) * n *
+                     n)                                    // - {grad^2 v n n }
+                     * (fe_i.jump_gradient(j, qpoint) * n) // [grad u n]
+                   - (fe_i.average_hessian(j, qpoint) * n *
+                      n) // - {grad^2 u n n }
+                       * (fe_i.jump_gradient(i, qpoint) * n) // [grad v n]
+                   // gamma [grad u n ][grad v n]:
+                   + gamma * (fe_i.jump_gradient(i, qpoint) * n) *
+                       (fe_i.jump_gradient(j, qpoint) * n)) *
+                  JxW[qpoint]; // dx
+              }
+        }
+    };
+
+
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData<dim> &  scratch_data,
+                               CopyData &          copy_data) {
+      // return;
+      FEInterfaceValues<dim> &fe_i = scratch_data.fe_interface_values;
+      fe_i.reinit(cell, face_no);
+      const auto &q_points = fe_i.get_quadrature_points();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      const unsigned int n_dofs        = fe_i.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices = fe_i.get_interface_dof_indices();
+
+      copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+      const std::vector<double> &        JxW     = fe_i.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_i.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);
+
+
+      // eta = 1/2 + 2C_2
+      // gamma = eta/|e|
+
+      double gamma = 1.0;
+
+      {
+        int                degree = fe.tensor_degree();
+        const unsigned int normal1 =
+          GeometryInfo<dim>::unit_normal_direction[face_no];
+        const unsigned int deg1sq =
+          degree * (degree + 1); //(deg1 == 0) ? 1 : deg1 * (deg1+1);
+
+        gamma = deg1sq / cell->extent_in_direction(normal1);
+        //      gamma = 0.5*(penalty1 + penalty2);
+      }
+
+      for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+        {
+          const auto &n = normals[qpoint];
+
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            {
+              for (unsigned int j = 0; j < n_dofs; ++j)
+                copy_data_face.cell_matrix(i, j) +=
+                  (-(fe_i.average_hessian(i, qpoint) * n *
+                     n)                                    // - {grad^2 v n n }
+                     * (fe_i.jump_gradient(j, qpoint) * n) // [grad u n]
+                   //
+                   - (fe_i.average_hessian(j, qpoint) * n *
+                      n) // - {grad^2 u n n }
+                       * (fe_i.jump_gradient(i, qpoint) * n) //  [grad v n]
+                                                             //
+                   + 2.0 * gamma *
+                       (fe_i.jump_gradient(i, qpoint) * n) // 2 gamma [grad v n]
+                       * (fe_i.jump_gradient(j, qpoint) * n)) // [grad u n]
+                  * JxW[qpoint];                              // dx
+
+
+              copy_data.cell_rhs(i) +=
+                (-(fe_i.average_hessian(i, qpoint) * n *
+                   n) *                                    //  - {grad^2 v n n }
+                   (exact_gradients[qpoint] * n)           // (grad u_exact n)
+                 + 2.0 * gamma                             //
+                     * (fe_i.jump_gradient(i, qpoint) * n) // [grad v n]
+                     * (exact_gradients[qpoint] * n)       // (grad u_exact n)
+                 ) *
+                JxW[qpoint]; // dx
+            }
+        }
+    };
+
+    auto copier = [&](const CopyData &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             system_rhs);
+
+      for (auto &cdf : c.face_data)
+        {
+          constraints.distribute_local_to_global(cdf.cell_matrix,
+                                                 cdf.joint_dof_indices,
+                                                 system_matrix);
+        }
+    };
+
+    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;
+    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);
+  }
+
+
+  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);
+  }
+
+
+
+  template <int dim>
+  void BiharmonicProblem<dim>::compute_errors()
+  {
+    const unsigned int n_gauss_points =
+      dof_handler.get_fe().tensor_degree() + 1;
+
+    {
+      Vector<float> norm_per_cell(triangulation.n_active_cells());
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        ExactSolution::Solution<dim>(),
+                                        norm_per_cell,
+                                        QGauss<dim>(n_gauss_points + 1),
+                                        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>(n_gauss_points + 1),
+                                        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;
+    }
+
+    // Now also compute the H2 seminorm error, integrating over the interiors
+    // of the cells but not taking into account the interface jump terms.
+    // This is *not* equivalent to the energy error for the problem.
+    {
+      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 diff = 0;
+          for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+            {
+              diff +=
+                ((exact_hessians[q_point] - hessians[q_point]).norm_square() *
+                 fe_values.JxW(q_point));
+            }
+          error_per_cell[cell->active_cell_index()] = std::sqrt(diff);
+        }
+      const double error_norm = error_per_cell.l2_norm();
+      std::cout << "   Error in the broken H2 seminorm: " << error_norm
+                << std::endl;
+    }
+  }
+
+
+  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, "u");
+    Vector<double>                     exact  = solution;
+    unsigned int                       degree = fe.tensor_degree();
+    const ExactSolution::Solution<dim> exact_solution;
+    VectorTools::project(mapping,
+                         dof_handler,
+                         constraints,
+                         QGauss<dim>(degree + 1),
+                         exact_solution,
+                         exact);
+    data_out.add_data_vector(exact, "exact");
+
+    data_out.build_patches();
+
+    std::ofstream output_vtk(
+      ("output_" + Utilities::int_to_string(iteration, 6) + ".vtk").c_str());
+    data_out.write_vtk(output_vtk);
+  }
+
+
+
+  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 StepBiharmonic
+
+
+
+int main(int argc, char *argv[])
+{
+  try
+    {
+      using namespace dealii;
+      using namespace StepBiharmonic;
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
+
+      unsigned int degree = 2; // minimum degree 2
+
+      // If provided on the command line, override the polynomial degree
+      // by the one given there.
+      if (argc > 1)
+        degree = Utilities::string_to_int(argv[1]);
+
+      BiharmonicProblem<2> my_bi(degree);
+      my_bi.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;
+}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.