]> https://gitweb.dealii.org/ - dealii.git/commitdiff
add tests/feinterface/stokes_simplex 17211/head
authorTimo Heister <timo.heister@gmail.com>
Fri, 5 Jul 2024 10:13:20 +0000 (06:13 -0400)
committerTimo Heister <timo.heister@gmail.com>
Fri, 5 Jul 2024 10:22:45 +0000 (06:22 -0400)
This is adding a meshworker simplex test for DG Stokes, extracted from
#13868

tests/feinterface/stokes_simplex.cc [new file with mode: 0644]
tests/feinterface/stokes_simplex.with_umfpack=on.output [new file with mode: 0644]

diff --git a/tests/feinterface/stokes_simplex.cc b/tests/feinterface/stokes_simplex.cc
new file mode 100644 (file)
index 0000000..ae3189c
--- /dev/null
@@ -0,0 +1,923 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2019 - 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+// test FEInterfaceValues for a DG Stokes problem on a
+// simplex mesh.
+
+#include <deal.II/base/flow_function.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/dofs/block_info.h>
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_bdm.h>
+#include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_raviart_thomas.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/sparse_direct.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+#include <deal.II/meshworker/scratch_data.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <sstream>
+
+#include "../tests.h"
+
+#define USE_SIMPLEX
+
+namespace StokesTests
+{
+  using namespace dealii;
+
+  struct CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+    std::array<double, 2>                values;
+    std::array<unsigned int, 2>          cell_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;
+    double                               value;
+    unsigned int                         cell_index;
+
+
+
+    template <class Iterator>
+    void
+    reinit(const Iterator &cell, const 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 <class MatrixType, class VectorType>
+  inline void
+  copy(const CopyData                  &c,
+       const AffineConstraints<double> &constraints,
+       MatrixType                      &system_matrix,
+       VectorType                      &system_rhs)
+  {
+    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)
+      {
+        const unsigned int dofs_per_cell = cdf.joint_dof_indices.size();
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          for (unsigned int k = 0; k < dofs_per_cell; ++k)
+            system_matrix.add(cdf.joint_dof_indices[i],
+                              cdf.joint_dof_indices[k],
+                              cdf.cell_matrix(i, k));
+      }
+  }
+
+  template <int dim>
+  class Solution : public Function<dim>
+  {
+  public:
+    Solution()
+      : Function<dim>(dim + 1)
+    {}
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const;
+    virtual Tensor<1, dim>
+    gradient(const Point<dim> &p, const unsigned int component = 0) const;
+  };
+
+  template <>
+  double
+  Solution<2>::value(const Point<2> &p, const unsigned int component) const
+  {
+    using numbers::PI;
+    const double x = p(0);
+    const double y = p(1);
+    // zero on BD's
+    if (component == 0)
+      return PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y);
+    if (component == 1)
+      return -PI * sin(PI * y) * sin(PI * y) * sin(2.0 * PI * x);
+    if (component == 2)
+      return cos(PI * x) * sin(PI * y);
+
+    return 0;
+  }
+
+  template <>
+  double
+  Solution<3>::value(const Point<3> &p, const unsigned int component) const
+  {
+    Assert(component <= 3 + 1, ExcIndexRange(component, 0, 3 + 1));
+
+    using numbers::PI;
+    const double x = p(0);
+    const double y = p(1);
+    const double z = p(2);
+
+    if (component == 0)
+      return 2. * PI * sin(PI * x) * sin(PI * x) * sin(2.0 * PI * y) *
+             sin(2.0 * PI * z);
+    if (component == 1)
+      return -PI * sin(PI * y) * sin(PI * y) * sin(2.0 * PI * x) *
+             sin(2.0 * PI * z);
+    if (component == 2)
+      return -PI * sin(PI * z) * sin(PI * z) * sin(2.0 * PI * x) *
+             sin(2.0 * PI * y);
+    if (component == 3)
+      return sin(PI * x) * cos(PI * y) * sin(PI * z);
+
+    return 0;
+  }
+
+  // Note that for the gradient we need to return a Tensor<1,dim>
+  template <>
+  Tensor<1, 2>
+  Solution<2>::gradient(const Point<2> &p, const unsigned int component) const
+  {
+    Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1));
+
+    using numbers::PI;
+    const double x = p(0);
+    const double y = p(1);
+
+    Tensor<1, 2> return_value;
+    if (component == 0)
+      {
+        return_value[0] = PI * PI * sin(2.0 * PI * y) * sin(2.0 * PI * x);
+        return_value[1] =
+          2.0 * PI * PI * sin(PI * x) * sin(PI * x) * cos(2.0 * PI * y);
+      }
+    else if (component == 1)
+      {
+        return_value[0] =
+          -2.0 * PI * PI * sin(PI * y) * sin(PI * y) * cos(2.0 * PI * x);
+        return_value[1] = -PI * PI * sin(2.0 * PI * y) * sin(2.0 * PI * x);
+      }
+    else if (component == 2)
+      {
+        return_value[0] = PI * cos(PI * x) * cos(PI * y);
+        return_value[1] = -PI * sin(PI * x) * sin(PI * y);
+      }
+
+    return return_value;
+  }
+
+  template <>
+  Tensor<1, 3>
+  Solution<3>::gradient(const Point<3> &p, const unsigned int component) const
+  {
+    Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1));
+
+    using numbers::PI;
+    const double x = p(0);
+    const double y = p(1);
+    const double z = p(2);
+
+    Tensor<1, 3> return_value;
+    if (component == 0)
+      {
+        return_value[0] =
+          2 * PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z);
+        return_value[1] = 4 * PI * PI * sin(PI * x) * sin(PI * x) *
+                          cos(2 * PI * y) * sin(2 * PI * z);
+        return_value[2] = 4 * PI * PI * sin(PI * x) * sin(PI * x) *
+                          cos(2 * PI * z) * sin(2 * PI * y);
+      }
+    else if (component == 1)
+      {
+        return_value[0] = -2 * PI * PI * sin(PI * y) * sin(PI * y) *
+                          cos(2 * PI * x) * sin(2 * PI * z);
+        return_value[1] =
+          -PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z);
+        return_value[2] = -2 * PI * PI * sin(PI * y) * sin(PI * y) *
+                          cos(2 * PI * z) * sin(2 * PI * x);
+      }
+    else if (component == 2)
+      {
+        return_value[0] = -2 * PI * PI * sin(PI * z) * sin(PI * z) *
+                          cos(2 * PI * x) * sin(2 * PI * y);
+        return_value[1] = -2 * PI * PI * sin(PI * z) * sin(PI * z) *
+                          cos(2 * PI * y) * sin(2 * PI * x);
+        return_value[2] =
+          -PI * PI * sin(2 * PI * x) * sin(2 * PI * y) * sin(2 * PI * z);
+      }
+    else if (component == 3)
+      {
+        return_value[0] = PI * cos(PI * x) * cos(PI * y) * sin(PI * z);
+        return_value[1] = -PI * sin(PI * x) * sin(PI * y) * sin(PI * z);
+        return_value[2] = PI * sin(PI * x) * cos(PI * y) * cos(PI * z);
+      }
+
+    return return_value;
+  }
+
+
+
+  // Implementation of $f$. See the introduction for more information.
+  template <int dim>
+  class RightHandSide : public Function<dim>
+  {
+  public:
+    RightHandSide()
+      : Function<dim>(dim + 1)
+    {}
+
+    virtual double
+    value(const Point<dim> &p, const unsigned int component = 0) const;
+  };
+
+  template <>
+  double
+  RightHandSide<2>::value(const Point<2> &p, const unsigned int component) const
+  {
+    Assert(component <= 2, ExcIndexRange(component, 0, 2 + 1));
+
+    using numbers::PI;
+    double x  = p(0);
+    double y  = p(1);
+    double nu = 1.0;
+
+    // RHS for 0 BD's
+    if (component == 0)
+      return -nu * 2.0 * PI * PI * PI *
+               (-2.0 * sin(PI * x) * sin(PI * x) + cos(2. * PI * x)) *
+               sin(2.0 * PI * y) -
+             PI * sin(PI * x) * sin(PI * y);
+    if (component == 1)
+      return nu * 2.0 * PI * PI * PI * (2.0 * cos(2.0 * PI * y) - 1) *
+               sin(2.0 * PI * x) +
+             PI * cos(PI * x) * cos(PI * y);
+    if (component == 2)
+      return 0.0;
+
+    return 0.0;
+  }
+
+  template <>
+  double
+  RightHandSide<3>::value(const Point<3> &p, const unsigned int component) const
+  {
+    Assert(component <= 3, ExcIndexRange(component, 0, 3 + 1));
+
+    using numbers::PI;
+    double x = p(0);
+    double y = p(1);
+    double z = p(2);
+
+    if (component == 0)
+      return 4. * PI * PI * PI *
+               (4. * sin(PI * x) * sin(PI * x) - cos(2. * PI * x)) *
+               sin(2. * PI * y) * sin(2. * PI * z) +
+             PI * cos(PI * x) * cos(PI * y) * sin(PI * z);
+    if (component == 1)
+      return -2. * PI * PI * PI *
+               (4. * sin(PI * y) * sin(PI * y) - cos(2. * PI * y)) *
+               sin(2. * PI * x) * sin(2. * PI * z) +
+             PI * (-1) * sin(PI * y) * sin(PI * x) * sin(PI * z);
+    if (component == 2)
+      return -2. * PI * PI * PI *
+               (4. * sin(PI * z) * sin(PI * z) - cos(2. * PI * z)) *
+               sin(2. * PI * x) * sin(2. * PI * y) +
+             PI * cos(PI * z) * sin(PI * x) * cos(PI * y);
+    if (component == 3)
+      return 0.0;
+
+    return 0.0;
+  }
+
+
+
+  // @sect3{The StokesProblem class}
+  //
+  // This is the main class of the problem.
+  template <int dim>
+  class StokesProblem
+  {
+  public:
+    StokesProblem(FiniteElement<dim> &fe, const unsigned int pressure_degree);
+    void
+    run();
+
+  private:
+    void
+    setup_dofs();
+    void
+    assemble_system_mesh_loop();
+    void
+    solve();
+    void
+    compute_errors(unsigned int k);
+
+    const unsigned int pressure_degree;
+
+#ifdef USE_SIMPLEX
+    const QGaussSimplex<dim>     quadrature;
+    const QGaussSimplex<dim - 1> face_quadrature;
+    const QGaussSimplex<dim>     quadrature_2;
+    const QGaussSimplex<dim - 1> face_quadrature_2;
+    const MappingFE<dim>         mapping;
+#else
+    const QGauss<dim>     quadrature;
+    const QGauss<dim - 1> face_quadrature;
+    const QGauss<dim>     quadrature_2;
+    const QGauss<dim - 1> face_quadrature_2;
+    const MappingQ1<dim>  mapping;
+#endif
+
+    Triangulation<dim>  triangulation;
+    FiniteElement<dim> &fe;
+    DoFHandler<dim>     dof_handler;
+
+    AffineConstraints<> constraints;
+
+    BlockSparsityPattern      sparsity_pattern;
+    BlockSparseMatrix<double> system_matrix;
+
+    BlockVector<double> solution;
+    BlockVector<double> system_rhs;
+
+    double last_l2_error;
+    double last_H1_error;
+    double last_Hdiv_error1;
+    double last_Hdiv_error2;
+  };
+
+
+
+  template <int dim>
+  StokesProblem<dim>::StokesProblem(FiniteElement<dim> &fe,
+                                    const unsigned int  pressure_degree)
+    : pressure_degree(pressure_degree)
+    , quadrature(pressure_degree + 1)
+    , face_quadrature(pressure_degree + 1)
+    , quadrature_2(pressure_degree + 2)
+    , face_quadrature_2(pressure_degree + 2)
+#ifdef USE_SIMPLEX
+    , mapping(FE_SimplexP<dim>(1))
+#else
+    , mapping()
+#endif
+    , fe(fe)
+    , dof_handler(triangulation)
+  {}
+
+  template <int dim>
+  void
+  StokesProblem<dim>::setup_dofs()
+  {
+    system_matrix.clear();
+
+    // The main DoFHandler only needs active DoFs, so we are not calling
+    // distribute_mg_dofs() here
+    dof_handler.distribute_dofs(fe);
+
+    // This block structure separates the dim velocity components from
+    // the pressure component (used for reordering). Note that we have
+    // 2 instead of dim+1 blocks like in step-22, because our FESystem
+    // is nested and the dim velocity components appear as one block.
+    std::vector<unsigned int> block_component(2);
+    block_component[0] = 0;
+    block_component[1] = 1;
+
+    // Velocities start at component 0:
+    const FEValuesExtractors::Vector velocities(0);
+    DoFRenumbering::block_wise(dof_handler);
+
+    std::vector<types::global_dof_index> dofs_per_block =
+      DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
+    const unsigned int n_u = dofs_per_block[0], n_p = dofs_per_block[1];
+
+    deallog << "\tNumber of active cells: " << triangulation.n_active_cells()
+            << std::endl
+            << "\tNumber of degrees of freedom: " << dof_handler.n_dofs()
+            << " (" << n_u << '+' << n_p << ')' << std::endl;
+
+    {
+      constraints.reinit();
+      DoFTools::make_hanging_node_constraints(dof_handler, constraints);
+      VectorTools::interpolate_boundary_values(dof_handler,
+                                               0,
+                                               Solution<dim>(),
+                                               constraints,
+                                               fe.component_mask(velocities));
+      constraints.close();
+    }
+
+    {
+      BlockDynamicSparsityPattern csp(dofs_per_block, dofs_per_block);
+      DoFTools::make_flux_sparsity_pattern(dof_handler, csp, constraints);
+      sparsity_pattern.copy_from(csp);
+    }
+    system_matrix.reinit(sparsity_pattern);
+
+    solution.reinit(dofs_per_block);
+    system_rhs.reinit(dofs_per_block);
+  }
+
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::assemble_system_mesh_loop()
+  {
+    system_matrix = 0;
+    system_rhs    = 0;
+
+    typedef decltype(dof_handler.begin_active()) Iterator;
+    const RightHandSide<dim>                     rhs_function;
+    const Solution<dim>                          boundary_function;
+    const FEValuesExtractors::Vector             velocities(0);
+    const FEValuesExtractors::Scalar             pressure(dim);
+
+    const double nu = 1.0;
+
+    auto penalty_parameter = [](const double degree,
+                                const double extent1,
+                                const double extent2) -> double {
+      return 4.0 * degree * (degree + 1.0) * 0.5 *
+             (1.0 / extent1 + 1.0 / extent2);
+    };
+
+    auto cell_worker = [&](const Iterator               &cell,
+                           MeshWorker::ScratchData<dim> &scratch_data,
+                           CopyData                     &copy_data) {
+      const FEValues<dim> &fe_v = scratch_data.reinit(cell);
+
+      const unsigned int dofs_per_cell = fe_v.dofs_per_cell;
+      const unsigned int n_q_points    = fe_v.get_quadrature().size();
+
+      copy_data.reinit(cell, dofs_per_cell);
+
+      const std::vector<double> &JxW = fe_v.get_JxW_values();
+
+      const double                nu = 1.0;
+      std::vector<Vector<double>> rhs_values(n_q_points,
+                                             Vector<double>(dim + 1));
+      rhs_function.vector_value_list(fe_v.get_quadrature_points(), rhs_values);
+      Tensor<1, dim> force_f;
+
+      for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+        {
+          for (unsigned int d = 0; d < dim; ++d)
+            force_f[d] = rhs_values[point](d);
+          for (unsigned int i = 0; i < fe_v.dofs_per_cell; ++i)
+            {
+              for (unsigned int j = 0; j < fe_v.dofs_per_cell; ++j)
+                copy_data.cell_matrix(i, j) +=
+                  (
+                    // nu \nabla v : \nabla u
+                    nu * scalar_product(fe_v[velocities].gradient(i, point),
+                                        fe_v[velocities].gradient(j, point))
+                    // -q, div u
+                    - fe_v[pressure].value(i, point) *
+                        fe_v[velocities].divergence(j, point)
+                    // -p, div v
+                    - fe_v[pressure].value(j, point) *
+                        fe_v[velocities].divergence(i, point)) *
+                  JxW[point];
+
+              copy_data.cell_rhs(i) +=
+                // f,v
+                (force_f * fe_v[velocities].value(i, point)) * JxW[point];
+            }
+        }
+    };
+
+    auto boundary_worker = [&](const Iterator               &cell,
+                               const unsigned int           &face_no,
+                               MeshWorker::ScratchData<dim> &scratch_data,
+                               CopyData                     &copy_data) {
+      const FEFaceValuesBase<dim> &fe_fv = scratch_data.reinit(cell, face_no);
+
+      const auto &q_points = fe_fv.get_quadrature_points();
+
+      const std::vector<double>         &JxW     = fe_fv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
+
+      std::vector<Vector<double>> g_values(q_points.size(),
+                                           Vector<double>(dim + 1));
+      boundary_function.vector_value_list(q_points, g_values);
+      Tensor<1, dim> g;
+
+      const double degree =
+        std::max(1.0, static_cast<double>(fe_fv.get_fe().degree));
+      const double extent1 = cell->measure() / cell->face(face_no)->measure();
+      const double penalty = penalty_parameter(degree, extent1, extent1);
+
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          for (unsigned int d = 0; d < dim; ++d)
+            g[d] = g_values[point](d);
+
+          for (unsigned int i = 0; i < fe_fv.dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < fe_fv.dofs_per_cell; ++j)
+              copy_data.cell_matrix(i, j) +=
+                (
+                  // - nu (\nabla u n) . v
+                  -nu *
+                    ((fe_fv[velocities].gradient(j, point) * normals[point]) *
+                     fe_fv[velocities].value(i, point))
+
+                  // - nu u . (\nabla v n)  // NIPG: use +
+                  -
+                  nu * (fe_fv[velocities].value(j, point) *
+                        (fe_fv[velocities].gradient(i, point) * normals[point]))
+
+                  // + nu * penalty u . v
+                  + nu * penalty *
+                      (fe_fv[velocities].value(j, point) *
+                       fe_fv[velocities].value(i, point))
+
+                  // p (v.n)
+                  + fe_fv[pressure].value(j, point) *
+                      scalar_product(fe_fv[velocities].value(i, point),
+                                     normals[point])
+
+                  // q (u.n)
+                  + fe_fv[pressure].value(i, point) *
+                      scalar_product(fe_fv[velocities].value(j, point),
+                                     normals[point])
+
+                    ) *
+                JxW[point];
+
+          for (unsigned int i = 0; i < fe_fv.dofs_per_cell; ++i)
+            copy_data.cell_rhs(i) +=
+              (
+                // -nu g . (\nabla v n) // NIPG: use +
+                -nu * scalar_product(g,
+                                     (fe_fv[velocities].gradient(i, point) *
+                                      normals[point]))
+
+                // +nu penalty g . v
+                + nu * penalty *
+                    scalar_product(g, fe_fv[velocities].value(i, point))
+
+                // q (g.n) (weak normal component of boundary condition)
+                + fe_fv[pressure].value(i, point) *
+                    scalar_product(g, normals[point])) *
+              JxW[point];
+        }
+    };
+
+    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,
+                           MeshWorker::ScratchData<dim> &scratch_data,
+                           CopyData                     &copy_data) {
+      const FEInterfaceValues<dim> &fe_fv =
+        scratch_data.reinit(cell, f, sf, ncell, nf, nsf);
+      FEInterfaceViews::Scalar<dim> interface_scalar(fe_fv, pressure.component);
+      FEInterfaceViews::Vector<dim> interface_vector(
+        fe_fv, velocities.first_vector_component);
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace      &copy_data_face = copy_data.face_data.back();
+      const unsigned int dofs_per_cell  = fe_fv.n_current_interface_dofs();
+
+      copy_data_face.joint_dof_indices = fe_fv.get_interface_dof_indices();
+      copy_data_face.cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+
+      const std::vector<double>         &JxW     = fe_fv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_fv.get_normal_vectors();
+      const auto &q_points = fe_fv.get_quadrature_points();
+
+      double       nu = 1.0;
+      const double degree =
+        std::max(1.0, static_cast<double>(fe_fv.get_fe().degree));
+      const double extent1 = cell->measure() / cell->face(f)->measure();
+      const double extent2 = ncell->measure() / ncell->face(nf)->measure();
+      const double penalty = penalty_parameter(degree, extent1, extent2);
+
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              copy_data_face.cell_matrix(i, j) +=
+                (
+                  // - nu {\nabla u}n . [v] (consistency)
+                  -nu *
+                    (fe_fv[velocities].average_of_gradients(j, point) *
+                     normals[point]) *
+                    fe_fv[velocities].jump_in_values(i, point)
+
+                  // - nu [u] . {\nabla v}n  (symmetry) // NIPG: use +
+                  - nu * fe_fv[velocities].jump_in_values(j, point) *
+                      (fe_fv[velocities].average_of_gradients(i, point) *
+                       normals[point])
+
+                  // nu sigma [u].[v] (penalty)
+                  + nu * penalty *
+                      scalar_product(fe_fv[velocities].jump_in_values(j, point),
+                                     fe_fv[velocities].jump_in_values(i, point))
+
+                  // {p} ([v].n)
+                  + fe_fv[pressure].average_of_values(j, point) *
+                      scalar_product(fe_fv[velocities].jump_in_values(i, point),
+                                     normals[point])
+
+                  // {q} ([u].n)
+                  + fe_fv[pressure].average_of_values(i, point) *
+                      scalar_product(fe_fv[velocities].jump_in_values(j, point),
+                                     normals[point])) *
+                JxW[point];
+        }
+    };
+
+    auto copier = [&](const CopyData &c) {
+      copy(c, constraints, system_matrix, system_rhs);
+    };
+
+    const UpdateFlags cell_flags = update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values;
+    const UpdateFlags face_flags = update_values | update_gradients |
+                                   update_quadrature_points |
+                                   update_normal_vectors | update_JxW_values;
+
+
+    MeshWorker::ScratchData<dim> scratch_data(
+      mapping, fe, quadrature, cell_flags, face_quadrature, face_flags);
+    CopyData cd;
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          cd,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_boundary_faces |
+                            MeshWorker::assemble_own_interior_faces_once,
+                          boundary_worker,
+                          face_worker);
+  }
+
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::solve()
+  {
+    solution = 0;
+    SparseDirectUMFPACK A_direct;
+    A_direct.initialize(system_matrix);
+    A_direct.vmult(solution, system_rhs);
+
+    constraints.distribute(solution);
+
+    const double mean_pressure =
+      VectorTools::compute_mean_value(dof_handler, quadrature_2, solution, dim);
+    VectorTools::add_constant(solution, dof_handler, dim, -mean_pressure);
+  }
+
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::compute_errors(unsigned int k)
+  {
+    const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
+                                                     dim + 1);
+    const ComponentSelectFunction<dim> pressure_mask(dim, 1.0, dim + 1);
+
+    Vector<float> difference_per_cell(triangulation.n_active_cells());
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      quadrature_2,
+                                      VectorTools::L2_norm,
+                                      &velocity_mask);
+
+    const double Velocity_L2_error = difference_per_cell.l2_norm();
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      quadrature_2,
+                                      VectorTools::H1_norm,
+                                      &velocity_mask);
+
+    const double Velocity_H1_error = difference_per_cell.l2_norm();
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Functions::ZeroFunction<dim>(dim + 1),
+                                      difference_per_cell,
+                                      quadrature_2,
+                                      VectorTools::Hdiv_seminorm,
+                                      &velocity_mask);
+
+    const double Velocity_Hdiv_error1 = difference_per_cell.l2_norm();
+
+    static double last_Pressure_L2_error = 0;
+
+    VectorTools::integrate_difference(mapping,
+                                      dof_handler,
+                                      solution,
+                                      Solution<dim>(),
+                                      difference_per_cell,
+                                      quadrature_2,
+                                      VectorTools::L2_norm,
+                                      &pressure_mask);
+    const double Pressure_L2_error =
+      VectorTools::compute_global_error(triangulation,
+                                        difference_per_cell,
+                                        VectorTools::L2_norm);
+
+    deallog << " At " << k + 1 << "th mesh" << std::endl
+            << " L2 error:     " << std::setw(12) << Velocity_L2_error
+            << std::setw(0) << " L2_Conv_rate: " << std::setw(6)
+            << (k == 0 ? 0 : last_l2_error / Velocity_L2_error) << std::endl
+            << " H1 error:     " << std::setw(12) << Velocity_H1_error
+            << std::setw(0) << " H1_Conv_rate: " << std::setw(6)
+            << (k == 0 ? 0 : last_H1_error / Velocity_H1_error) << std::endl
+            << " Hdiv error1:  " << std::setw(12) << Velocity_Hdiv_error1
+            << std::setw(0) << " Hdiv_Conv_rate1: " << std::setw(6)
+            << (k == 0 ? 0 : last_Hdiv_error1 / Velocity_Hdiv_error1)
+            << std::endl
+            << " L2 pressure:  " << std::setw(12) << Pressure_L2_error
+            << std::setw(0) << " rate: " << std::setw(6)
+            << (k == 0 ? 0 : last_Pressure_L2_error / Pressure_L2_error)
+            << std::endl
+            << std::setw(0) << std::endl;
+    last_l2_error          = Velocity_L2_error;
+    last_H1_error          = Velocity_H1_error;
+    last_Hdiv_error1       = Velocity_Hdiv_error1;
+    last_Pressure_L2_error = Pressure_L2_error;
+  }
+
+
+
+  template <int dim>
+  void
+  StokesProblem<dim>::run()
+  {
+#ifdef USE_SIMPLEX
+    Triangulation<dim> temp;
+    GridGenerator::hyper_cube(temp);
+    GridGenerator::convert_hypercube_to_simplex_mesh(temp, triangulation);
+#else
+    GridGenerator::hyper_cube(triangulation);
+#endif
+
+    triangulation.refine_global(1);
+
+    deallog << "  Now running with " << fe.get_name() << std::endl;
+
+    for (unsigned int refinement_cycle = 0; refinement_cycle < 4;
+         ++refinement_cycle)
+      {
+        if (refinement_cycle > 0)
+          triangulation.refine_global();
+
+        setup_dofs();
+
+        assemble_system_mesh_loop();
+
+        solve();
+
+        compute_errors(refinement_cycle);
+      }
+  }
+} // namespace StokesTests
+
+
+template <int dim>
+void
+test_stokes(FiniteElement<dim> &fe, const unsigned int pressure_degree)
+{
+  deallog << fe.get_name() << ": degree=" << fe.degree
+          << " tensor_degree=" << fe.tensor_degree() << std::endl;
+  StokesTests::StokesProblem<dim> flow_problem(fe, pressure_degree);
+  flow_problem.run();
+}
+
+int
+main()
+{
+  try
+    {
+      using namespace dealii;
+      using namespace StokesTests;
+      const int dim = 2;
+
+      initlog();
+
+      std::unique_ptr<FiniteElement<dim>> fe;
+      const int                           degree = 2;
+
+      Assert(degree >= 1, ExcMessage("invalid degree!"));
+      {
+#ifdef USE_SIMPLEX
+        fe = std::make_unique<FESystem<dim>>(
+          FESystem<dim>(FE_SimplexDGP<dim>(degree), dim),
+          1,
+          FE_SimplexDGP<dim>(degree - 1),
+          1);
+#else
+        fe = std::make_unique<FESystem<dim>>(FESystem<dim>(FE_DGQ<dim>(degree),
+                                                           dim),
+                                             1,
+                                             FE_DGQ<dim>(degree - 1),
+                                             1);
+#endif
+        test_stokes(*fe, degree);
+      }
+    }
+  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;
+}
diff --git a/tests/feinterface/stokes_simplex.with_umfpack=on.output b/tests/feinterface/stokes_simplex.with_umfpack=on.output
new file mode 100644 (file)
index 0000000..54396e3
--- /dev/null
@@ -0,0 +1,35 @@
+
+DEAL::FESystem<2>[FESystem<2>[FE_SimplexDGP<2>(2)^2]-FE_SimplexDGP<2>(1)]: degree=2 tensor_degree=2
+DEAL::  Now running with FESystem<2>[FESystem<2>[FE_SimplexDGP<2>(2)^2]-FE_SimplexDGP<2>(1)]
+DEAL:: Number of active cells: 32
+DEAL:: Number of degrees of freedom: 480 (384+96)
+DEAL:: At 1th mesh
+DEAL:: L2 error:     0.117546     L2_Conv_rate: 0.00000
+DEAL:: H1 error:     2.20951      H1_Conv_rate: 0.00000
+DEAL:: Hdiv error1:  0.430647     Hdiv_Conv_rate1: 0.00000
+DEAL:: L2 pressure:  4.15183      rate: 0.00000
+DEAL::
+DEAL:: Number of active cells: 128
+DEAL:: Number of degrees of freedom: 1920 (1536+384)
+DEAL:: At 2th mesh
+DEAL:: L2 error:     0.0138896    L2_Conv_rate: 8.46286
+DEAL:: H1 error:     0.577807     H1_Conv_rate: 3.82396
+DEAL:: Hdiv error1:  0.108836     Hdiv_Conv_rate1: 3.95685
+DEAL:: L2 pressure:  1.41599      rate: 2.93210
+DEAL::
+DEAL:: Number of active cells: 512
+DEAL:: Number of degrees of freedom: 7680 (6144+1536)
+DEAL:: At 3th mesh
+DEAL:: L2 error:     0.00130976   L2_Conv_rate: 10.6047
+DEAL:: H1 error:     0.137071     H1_Conv_rate: 4.21539
+DEAL:: Hdiv error1:  0.0288825    Hdiv_Conv_rate1: 3.76823
+DEAL:: L2 pressure:  0.380977     rate: 3.71674
+DEAL::
+DEAL:: Number of active cells: 2048
+DEAL:: Number of degrees of freedom: 30720 (24576+6144)
+DEAL:: At 4th mesh
+DEAL:: L2 error:     0.000134320  L2_Conv_rate: 9.75105
+DEAL:: H1 error:     0.0329769    H1_Conv_rate: 4.15658
+DEAL:: Hdiv error1:  0.00739000   Hdiv_Conv_rate1: 3.90832
+DEAL:: L2 pressure:  0.0977908    rate: 3.89584
+DEAL::

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.