]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Enable mixed meshes for MF and DG 11416/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 19 Dec 2020 14:47:41 +0000 (15:47 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Wed, 10 Mar 2021 09:12:23 +0000 (10:12 +0100)
source/fe/mapping_fe.cc
tests/simplex/matrix_free_face_integral_01.cc [new file with mode: 0644]
tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/step-12b.cc [new file with mode: 0644]
tests/simplex/step-12b.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/step-12c.cc [new file with mode: 0644]
tests/simplex/step-12c.with_simplex_support=on.output [new file with mode: 0644]

index 3a4cc76024e54aa8639a0d96475440e6c5e2bdc5..e0e68633339713398e9465bea3c9724ed63089fb 100644 (file)
@@ -1697,7 +1697,13 @@ namespace internal
         const typename Mapping<dim, spacedim>::InternalDataBase &mapping_data,
         const ArrayView<Tensor<rank, spacedim>> &                output)
       {
-        AssertDimension(input.size(), output.size());
+        // In the case of wedges and pyramids, faces might have different
+        // numbers of quadrature points on each face with the result
+        // that input and output have different sizes, since input has
+        // the correct size but the size of output is the maximum of
+        // all possible sizes.
+        AssertIndexRange(input.size(), output.size() + 1);
+
         Assert(
           (dynamic_cast<
              const typename dealii::MappingFE<dim, spacedim>::InternalData *>(
@@ -1717,7 +1723,7 @@ namespace internal
                   typename FEValuesBase<dim>::ExcAccessToUninitializedField(
                     "update_contravariant_transformation"));
 
-                for (unsigned int i = 0; i < output.size(); ++i)
+                for (unsigned int i = 0; i < input.size(); ++i)
                   output[i] =
                     apply_transformation(data.contravariant[i], input[i]);
 
@@ -1738,7 +1744,7 @@ namespace internal
                 if (rank != 1)
                   return;
 
-                for (unsigned int i = 0; i < output.size(); ++i)
+                for (unsigned int i = 0; i < input.size(); ++i)
                   {
                     output[i] =
                       apply_transformation(data.contravariant[i], input[i]);
@@ -1756,7 +1762,7 @@ namespace internal
                   typename FEValuesBase<dim>::ExcAccessToUninitializedField(
                     "update_covariant_transformation"));
 
-                for (unsigned int i = 0; i < output.size(); ++i)
+                for (unsigned int i = 0; i < input.size(); ++i)
                   output[i] = apply_transformation(data.covariant[i], input[i]);
 
                 return;
diff --git a/tests/simplex/matrix_free_face_integral_01.cc b/tests/simplex/matrix_free_face_integral_01.cc
new file mode 100644 (file)
index 0000000..418ce55
--- /dev/null
@@ -0,0 +1,238 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Compute the right-hand side vector due to a Neumann boundary condition on
+// different meshes (pure simplex, wedge, pyramid mesh).
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/solver_control.h>
+#include <deal.II/lac/sparse_matrix.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "./simplex_grids.h"
+
+using namespace dealii;
+
+
+template <int dim>
+void
+test(const unsigned int v, const unsigned int degree)
+{
+  Triangulation<dim> tria;
+
+  std::shared_ptr<FiniteElement<dim>> fe;
+  std::shared_ptr<Quadrature<dim>>    quad;
+  std::shared_ptr<FiniteElement<dim>> fe_mapping;
+  hp::QCollection<dim - 1>            face_quad;
+
+  if (v == 0)
+    {
+      GridGenerator::subdivided_hyper_cube_with_simplices(tria,
+                                                          dim == 2 ? 16 : 8);
+      fe         = std::make_shared<FE_SimplexP<dim>>(degree);
+      quad       = std::make_shared<QGaussSimplex<dim>>(degree + 1);
+      fe_mapping = std::make_shared<FE_SimplexP<dim>>(1);
+      face_quad  = hp::QCollection<dim - 1>{QGaussSimplex<dim - 1>(degree + 1)};
+    }
+  else if (v == 1)
+    {
+      GridGenerator::subdivided_hyper_cube_with_wedges(tria, dim == 2 ? 16 : 8);
+      fe         = std::make_shared<FE_WedgeP<dim>>(degree);
+      quad       = std::make_shared<QGaussWedge<dim>>(degree + 1);
+      fe_mapping = std::make_shared<FE_WedgeP<dim>>(1);
+      face_quad  = hp::QCollection<dim - 1>{QGaussSimplex<dim - 1>(degree + 1),
+                                           QGaussSimplex<dim - 1>(degree + 1),
+                                           QGauss<dim - 1>(degree + 1),
+                                           QGauss<dim - 1>(degree + 1),
+                                           QGauss<dim - 1>(degree + 1)};
+    }
+  else if (v == 2)
+    {
+      GridGenerator::subdivided_hyper_cube_with_pyramids(tria,
+                                                         dim == 2 ? 16 : 8);
+      fe         = std::make_shared<FE_PyramidP<dim>>(degree);
+      quad       = std::make_shared<QGaussPyramid<dim>>(degree + 1);
+      fe_mapping = std::make_shared<FE_PyramidP<dim>>(1);
+      face_quad  = hp::QCollection<dim - 1>{QGauss<dim - 1>(degree + 1),
+                                           QGaussSimplex<dim - 1>(degree + 1),
+                                           QGaussSimplex<dim - 1>(degree + 1),
+                                           QGaussSimplex<dim - 1>(degree + 1),
+                                           QGaussSimplex<dim - 1>(degree + 1)};
+    }
+  else
+    Assert(false, ExcNotImplemented());
+
+  MappingFE<dim> mapping(*fe_mapping);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(*fe);
+
+  AffineConstraints<double> constraints;
+
+  using VectorType = LinearAlgebra::distributed::Vector<double>;
+  using number     = double;
+
+  VectorType vec_mf, vec_mb;
+
+  vec_mf.reinit(dof_handler.n_dofs());
+  vec_mb.reinit(dof_handler.n_dofs());
+
+  // compute vector with FEFaceEvaluation
+  {
+    typename MatrixFree<dim, double>::AdditionalData additional_data;
+    additional_data.mapping_update_flags = update_gradients | update_values;
+    additional_data.mapping_update_flags_boundary_faces =
+      update_gradients | update_values;
+
+    MatrixFree<dim, double> matrix_free;
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, *quad, additional_data);
+
+    matrix_free.template loop<VectorType, VectorType>(
+      [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
+        (void)data;
+        (void)dst;
+        (void)src;
+        (void)cell_range;
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        (void)data;
+        (void)dst;
+        (void)src;
+        (void)face_range;
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, number> fe_eval(data, face_range, true);
+        for (unsigned int face = face_range.first; face < face_range.second;
+             face++)
+          {
+            fe_eval.reinit(face);
+            for (unsigned int q = 0; q < fe_eval.n_q_points; ++q)
+              fe_eval.submit_value(1.0, q);
+
+            fe_eval.integrate_scatter(true, false, dst);
+          }
+      },
+      vec_mf,
+      vec_mf);
+  }
+
+  // compute vector with FEFaceValues
+  {
+    const UpdateFlags flag = update_JxW_values | update_values |
+                             update_gradients | update_quadrature_points;
+    FEValues<dim> fe_values(mapping, *fe, *quad, flag);
+
+    FEFaceValues<dim> fe_face_values(mapping, *fe, face_quad, flag);
+
+    const unsigned int dofs_per_cell = fe->dofs_per_cell;
+
+    std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+    Vector<double>                       cell_rhs(dofs_per_cell);
+
+    for (const auto &cell : dof_handler.cell_iterators())
+      {
+        if (!cell->is_locally_owned())
+          continue;
+
+        fe_values.reinit(cell);
+        cell_rhs = 0;
+
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary())
+            {
+              fe_face_values.reinit(cell, face);
+              for (const auto q : fe_face_values.quadrature_point_indices())
+                for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                  cell_rhs(i) +=
+                    (1.0 *                              // 1.0
+                     fe_face_values.shape_value(i, q) * // phi_i(x_q)
+                     fe_face_values.JxW(q));            // dx
+            }
+
+        cell->get_dof_indices(dof_indices);
+
+        constraints.distribute_local_to_global(cell_rhs, dof_indices, vec_mb);
+      }
+  }
+
+#if false
+  vec_mf.print(deallog.get_file_stream());
+  vec_mb.print(deallog.get_file_stream());
+#endif
+
+  for (unsigned int i = 0; i < vec_mf.size(); ++i)
+    Assert(std::abs(vec_mf[i] - vec_mb[i]) < 1e-8, ExcNotImplemented());
+
+  deallog << " dim=" << dim << " degree=" << degree << ": OK" << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  deallog.depth_file(2);
+
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+
+  for (unsigned int i = 0; i <= 2; ++i)
+    {
+      if (i == 0)
+        deallog.push("SIMPLEX");
+      else if (i == 1)
+        deallog.push("WEDGE  ");
+      else if (i == 2)
+        deallog.push("PYRAMID");
+      else
+        Assert(false, ExcNotImplemented());
+
+      if (i == 0) // 2D makes only sense for simplex
+        {
+          test<2>(i, /*degree=*/1);
+          test<2>(i, /*degree=*/2);
+        }
+
+      test<3>(i, /*degree=*/1);
+
+      if (i != 2) // for pyramids: no quadratic elements implemented
+        test<3>(i, /*degree=*/2);
+
+      deallog.pop();
+    }
+}
diff --git a/tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output b/tests/simplex/matrix_free_face_integral_01.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..e55e763
--- /dev/null
@@ -0,0 +1,8 @@
+
+DEAL:SIMPLEX:: dim=2 degree=1: OK
+DEAL:SIMPLEX:: dim=2 degree=2: OK
+DEAL:SIMPLEX:: dim=3 degree=1: OK
+DEAL:SIMPLEX:: dim=3 degree=2: OK
+DEAL:WEDGE  :: dim=3 degree=1: OK
+DEAL:WEDGE  :: dim=3 degree=2: OK
+DEAL:PYRAMID:: dim=3 degree=1: OK
diff --git a/tests/simplex/step-12b.cc b/tests/simplex/step-12b.cc
new file mode 100644 (file)
index 0000000..3673d19
--- /dev/null
@@ -0,0 +1,1151 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// a copy of tests/hp/step-12 testing different reference-cell types
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_values.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/derivative_approximation.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+#include "./simplex_grids.h"
+
+template <int dim>
+struct ScratchData
+{
+  hp::MappingCollection<dim>            mapping;
+  hp::FECollection<dim>                 fe;
+  hp::QCollection<dim>                  quadrature;
+  std::vector<hp::QCollection<dim - 1>> face_quadrature;
+
+  std::function<void(Triangulation<dim> &)> mesh_generator;
+};
+
+
+template <int dim>
+class RHS : public Function<dim>
+{
+public:
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int             component = 0) const;
+};
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int             component = 0) const;
+};
+
+
+template <int dim>
+class Beta
+{
+public:
+  Beta()
+  {}
+  void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<Point<dim>> &      values) const;
+};
+
+
+template <int dim>
+void
+RHS<dim>::value_list(const std::vector<Point<dim>> &points,
+                     std::vector<double> &          values,
+                     const unsigned int) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i)
+    values[i] = 0;
+}
+
+
+template <int dim>
+void
+Beta<dim>::value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Point<dim>> &      values) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i)
+    {
+      const Point<dim> &p    = points[i];
+      Point<dim> &      beta = values[i];
+
+      beta(0) = -p(1);
+      beta(1) = p(0);
+
+      if (beta.norm() > 1e-10)
+        beta /= std::sqrt(beta.square());
+    }
+}
+
+
+template <int dim>
+void
+BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+                                std::vector<double> &          values,
+                                const unsigned int) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i)
+    {
+      if (points[i](0) < 0.5)
+        values[i] = 1.;
+      else
+        values[i] = 0.;
+    }
+}
+
+
+template <int dim>
+class DGTransportEquation
+{
+public:
+  DGTransportEquation();
+
+  void
+  assemble_cell_term(const hp::FEValues<dim> &fe_v,
+                     FullMatrix<double> &     ui_vi_matrix,
+                     Vector<double> &         cell_vector) const;
+
+  void
+  assemble_boundary_term(const hp::FEFaceValues<dim> &fe_v,
+                         FullMatrix<double> &         ui_vi_matrix,
+                         Vector<double> &             cell_vector) const;
+
+  template <class X, class Y>
+  void
+  assemble_face_term1(const X &           fe_v,
+                      const Y &           fe_v_neighbor,
+                      FullMatrix<double> &ui_vi_matrix,
+                      FullMatrix<double> &ue_vi_matrix) const;
+
+  template <class X, class Y>
+  void
+  assemble_face_term2(const X &           fe_v,
+                      const Y &           fe_v_neighbor,
+                      FullMatrix<double> &ui_vi_matrix,
+                      FullMatrix<double> &ue_vi_matrix,
+                      FullMatrix<double> &ui_ve_matrix,
+                      FullMatrix<double> &ue_ve_matrix) const;
+
+private:
+  const Beta<dim>           beta_function;
+  const RHS<dim>            rhs_function;
+  const BoundaryValues<dim> boundary_function;
+};
+
+
+template <int dim>
+DGTransportEquation<dim>::DGTransportEquation()
+  : beta_function()
+  , rhs_function()
+  , boundary_function()
+{}
+
+
+template <int dim>
+void
+DGTransportEquation<dim>::assemble_cell_term(const hp::FEValues<dim> &fe_v,
+                                             FullMatrix<double> &ui_vi_matrix,
+                                             Vector<double> &cell_vector) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> rhs(fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+  rhs_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                          rhs);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+         ++i)
+      {
+        for (unsigned int j = 0; j < fe_v.get_present_fe_values().dofs_per_cell;
+             ++j)
+          ui_vi_matrix(i, j) -=
+            beta[point] * fe_v.get_present_fe_values().shape_grad(i, point) *
+            fe_v.get_present_fe_values().shape_value(j, point) * JxW[point];
+
+        cell_vector(i) += rhs[point] *
+                          fe_v.get_present_fe_values().shape_value(i, point) *
+                          JxW[point];
+      }
+}
+
+
+template <int dim>
+void
+DGTransportEquation<dim>::assemble_boundary_term(
+  const hp::FEFaceValues<dim> &fe_v,
+  FullMatrix<double> &         ui_vi_matrix,
+  Vector<double> &             cell_vector) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> g(fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+  boundary_function.value_list(
+    fe_v.get_present_fe_values().get_quadrature_points(), g);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int j = 0;
+               j < fe_v.get_present_fe_values().dofs_per_cell;
+               ++j)
+            ui_vi_matrix(i, j) +=
+              beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+      else
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          cell_vector(i) -= beta_n * g[point] *
+                            fe_v.get_present_fe_values().shape_value(i, point) *
+                            JxW[point];
+    }
+}
+
+
+template <int dim>
+template <class X, class Y>
+void
+DGTransportEquation<dim>::assemble_face_term1(
+  const X &           fe_v,
+  const Y &           fe_v_neighbor,
+  FullMatrix<double> &ui_vi_matrix,
+  FullMatrix<double> &ue_vi_matrix) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int j = 0;
+               j < fe_v.get_present_fe_values().dofs_per_cell;
+               ++j)
+            ui_vi_matrix(i, j) +=
+              beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+      else
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            ue_vi_matrix(i, k) +=
+              beta_n *
+              fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+    }
+}
+
+
+template <int dim>
+template <class X, class Y>
+void
+DGTransportEquation<dim>::assemble_face_term2(
+  const X &           fe_v,
+  const Y &           fe_v_neighbor,
+  FullMatrix<double> &ui_vi_matrix,
+  FullMatrix<double> &ue_vi_matrix,
+  FullMatrix<double> &ui_ve_matrix,
+  FullMatrix<double> &ue_ve_matrix) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        {
+          for (unsigned int i = 0;
+               i < fe_v.get_present_fe_values().dofs_per_cell;
+               ++i)
+            for (unsigned int j = 0;
+                 j < fe_v.get_present_fe_values().dofs_per_cell;
+                 ++j)
+              ui_vi_matrix(i, j) +=
+                beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+                fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            for (unsigned int j = 0;
+                 j < fe_v.get_present_fe_values().dofs_per_cell;
+                 ++j)
+              ui_ve_matrix(k, j) -=
+                beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+                fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+                JxW[point];
+        }
+      else
+        {
+          for (unsigned int i = 0;
+               i < fe_v.get_present_fe_values().dofs_per_cell;
+               ++i)
+            for (unsigned int l = 0;
+                 l < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+                 ++l)
+              ue_vi_matrix(i, l) +=
+                beta_n *
+                fe_v_neighbor.get_present_fe_values().shape_value(l, point) *
+                fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            for (unsigned int l = 0;
+                 l < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+                 ++l)
+              ue_ve_matrix(k, l) -=
+                beta_n *
+                fe_v_neighbor.get_present_fe_values().shape_value(l, point) *
+                fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+                JxW[point];
+        }
+    }
+}
+
+
+template <int dim>
+class DGMethod
+{
+public:
+  DGMethod(const ScratchData<dim> &scratch_data);
+  ~DGMethod();
+
+  void
+  run();
+
+private:
+  void
+  setup_system();
+  void
+  assemble_system1();
+  void
+  assemble_system2();
+  void
+  solve(Vector<double> &solution);
+  void
+  refine_grid();
+  void
+  output_results(const unsigned int cycle) const;
+
+  Triangulation<dim>               triangulation;
+  const hp::MappingCollection<dim> mapping;
+
+  hp::FECollection<dim> fe;
+  DoFHandler<dim>       dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  const hp::QCollection<dim>                  quadrature;
+  const std::vector<hp::QCollection<dim - 1>> face_quadrature;
+  const hp::QCollection<dim - 1>              face_quadrature_dummy;
+
+  Vector<double> solution1;
+  Vector<double> solution2;
+  Vector<double> right_hand_side;
+
+  const std::function<void(Triangulation<dim> &)> mesh_generator;
+
+  const DGTransportEquation<dim> dg;
+};
+
+
+template <int dim>
+DGMethod<dim>::DGMethod(const ScratchData<dim> &scratch_data)
+  : mapping(scratch_data.mapping)
+  , fe(scratch_data.fe)
+  , dof_handler(triangulation)
+  , quadrature(scratch_data.quadrature)
+  , face_quadrature(scratch_data.face_quadrature)
+  , mesh_generator(scratch_data.mesh_generator)
+  , dg()
+{}
+
+
+template <int dim>
+DGMethod<dim>::~DGMethod()
+{
+  dof_handler.clear();
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  sparsity_pattern.reinit(dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          (GeometryInfo<dim>::faces_per_cell *
+                             GeometryInfo<dim>::max_children_per_face +
+                           1) *
+                            fe[0].dofs_per_cell);
+
+  DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit(sparsity_pattern);
+
+  solution1.reinit(dof_handler.n_dofs());
+  solution2.reinit(dof_handler.n_dofs());
+  right_hand_side.reinit(dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::assemble_system1()
+{
+  const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell;
+  std::vector<types::global_dof_index> dofs(dofs_per_cell);
+  std::vector<types::global_dof_index> dofs_neighbor(dofs_per_cell);
+
+  const UpdateFlags update_flags = update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values;
+
+  const UpdateFlags face_update_flags =
+    update_values | update_quadrature_points | update_JxW_values |
+    update_normal_vectors;
+
+  const UpdateFlags neighbor_face_update_flags = update_values;
+
+  hp::FEValues<dim> fe_v(mapping, fe, quadrature, update_flags);
+
+  hp::FEFaceValues<dim>    fe_v_face(mapping,
+                                  fe,
+                                  face_quadrature,
+                                  face_update_flags);
+  hp::FESubfaceValues<dim> fe_v_subface(mapping,
+                                        fe,
+                                        face_quadrature_dummy,
+                                        face_update_flags);
+  hp::FEFaceValues<dim>    fe_v_face_neighbor(mapping,
+                                           fe,
+                                           face_quadrature,
+                                           neighbor_face_update_flags);
+  hp::FESubfaceValues<dim> fe_v_subface_neighbor(mapping,
+                                                 fe,
+                                                 face_quadrature_dummy,
+                                                 neighbor_face_update_flags);
+
+  FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
+
+  Vector<double> cell_vector(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+
+  for (; cell != endc; ++cell)
+    {
+      ui_vi_matrix = 0;
+      cell_vector  = 0;
+
+      fe_v.reinit(cell);
+
+      dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
+
+      cell->get_dof_indices(dofs);
+
+      for (const unsigned int face_no : cell->face_indices())
+        {
+          typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+
+          ue_vi_matrix = 0;
+
+          if (face->at_boundary())
+            {
+              fe_v_face.reinit(cell, face_no);
+
+              dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
+            }
+          else
+            {
+              typename DoFHandler<dim>::cell_iterator neighbor =
+                cell->neighbor(face_no);
+              ;
+
+              if (face->has_children())
+                {
+                  const unsigned int neighbor2 =
+                    cell->neighbor_of_neighbor(face_no);
+
+
+                  for (unsigned int subface_no = 0;
+                       subface_no < face->n_children();
+                       ++subface_no)
+                    {
+                      typename DoFHandler<dim>::active_cell_iterator
+                        neighbor_child =
+                          cell->neighbor_child_on_subface(face_no, subface_no);
+
+                      Assert(neighbor_child->face(neighbor2) ==
+                               face->child(subface_no),
+                             ExcInternalError());
+                      Assert(!neighbor_child->has_children(),
+                             ExcInternalError());
+
+                      ue_vi_matrix = 0;
+
+                      fe_v_subface.reinit(cell, face_no, subface_no);
+                      fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
+
+                      dg.assemble_face_term1(fe_v_subface,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix);
+
+                      neighbor_child->get_dof_indices(dofs_neighbor);
+
+                      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                        for (unsigned int k = 0; k < dofs_per_cell; ++k)
+                          system_matrix.add(dofs[i],
+                                            dofs_neighbor[k],
+                                            ue_vi_matrix(i, k));
+                    }
+                }
+              else
+                {
+                  if (neighbor->level() == cell->level())
+                    {
+                      const unsigned int neighbor2 =
+                        cell->neighbor_of_neighbor(face_no);
+
+                      fe_v_face.reinit(cell, face_no);
+                      fe_v_face_neighbor.reinit(neighbor, neighbor2);
+
+                      dg.assemble_face_term1(fe_v_face,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix);
+                    }
+                  else
+                    {
+                      Assert(neighbor->level() < cell->level(),
+                             ExcInternalError());
+
+                      const std::pair<unsigned int, unsigned int>
+                        faceno_subfaceno =
+                          cell->neighbor_of_coarser_neighbor(face_no);
+                      const unsigned int neighbor_face_no =
+                                           faceno_subfaceno.first,
+                                         neighbor_subface_no =
+                                           faceno_subfaceno.second;
+
+                      Assert(neighbor->neighbor_child_on_subface(
+                               neighbor_face_no, neighbor_subface_no) == cell,
+                             ExcInternalError());
+
+                      fe_v_face.reinit(cell, face_no);
+                      fe_v_subface_neighbor.reinit(neighbor,
+                                                   neighbor_face_no,
+                                                   neighbor_subface_no);
+
+                      dg.assemble_face_term1(fe_v_face,
+                                             fe_v_subface_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix);
+                    }
+
+                  neighbor->get_dof_indices(dofs_neighbor);
+
+                  for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                    for (unsigned int k = 0; k < dofs_per_cell; ++k)
+                      system_matrix.add(dofs[i],
+                                        dofs_neighbor[k],
+                                        ue_vi_matrix(i, k));
+                }
+            }
+        }
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+          system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j));
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        right_hand_side(dofs[i]) += cell_vector(i);
+    }
+}
+
+
+
+template <int dim>
+void
+DGMethod<dim>::assemble_system2()
+{
+  const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell;
+  std::vector<types::global_dof_index> dofs(dofs_per_cell);
+  std::vector<types::global_dof_index> dofs_neighbor(dofs_per_cell);
+
+  const UpdateFlags update_flags = update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values;
+
+  const UpdateFlags face_update_flags =
+    update_values | update_quadrature_points | update_JxW_values |
+    update_normal_vectors;
+
+  const UpdateFlags neighbor_face_update_flags = update_values;
+
+  hp::FEValues<dim>        fe_v(mapping, fe, quadrature, update_flags);
+  hp::FEFaceValues<dim>    fe_v_face(mapping,
+                                  fe,
+                                  face_quadrature,
+                                  face_update_flags);
+  hp::FESubfaceValues<dim> fe_v_subface(mapping,
+                                        fe,
+                                        face_quadrature_dummy,
+                                        face_update_flags);
+  hp::FEFaceValues<dim>    fe_v_face_neighbor(mapping,
+                                           fe,
+                                           face_quadrature,
+                                           neighbor_face_update_flags);
+
+
+  FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
+
+  FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);
+
+  Vector<double> cell_vector(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  for (; cell != endc; ++cell)
+    {
+      ui_vi_matrix = 0;
+      cell_vector  = 0;
+
+      fe_v.reinit(cell);
+
+      dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
+
+      cell->get_dof_indices(dofs);
+
+      for (const unsigned int face_no : cell->face_indices())
+        {
+          typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+
+          if (face->at_boundary())
+            {
+              fe_v_face.reinit(cell, face_no);
+
+              dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
+            }
+          else
+            {
+              Assert(cell->neighbor(face_no).state() == IteratorState::valid,
+                     ExcInternalError());
+              typename DoFHandler<dim>::cell_iterator neighbor =
+                cell->neighbor(face_no);
+              if (face->has_children())
+                {
+                  const unsigned int neighbor2 =
+                    cell->neighbor_of_neighbor(face_no);
+
+                  for (unsigned int subface_no = 0;
+                       subface_no < face->n_children();
+                       ++subface_no)
+                    {
+                      typename DoFHandler<dim>::cell_iterator neighbor_child =
+                        cell->neighbor_child_on_subface(face_no, subface_no);
+                      Assert(neighbor_child->face(neighbor2) ==
+                               face->child(subface_no),
+                             ExcInternalError());
+                      Assert(!neighbor_child->has_children(),
+                             ExcInternalError());
+
+                      ue_vi_matrix = 0;
+                      ui_ve_matrix = 0;
+                      ue_ve_matrix = 0;
+
+                      fe_v_subface.reinit(cell, face_no, subface_no);
+                      fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
+
+                      dg.assemble_face_term2(fe_v_subface,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix,
+                                             ui_ve_matrix,
+                                             ue_ve_matrix);
+
+                      neighbor_child->get_dof_indices(dofs_neighbor);
+
+                      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            system_matrix.add(dofs[i],
+                                              dofs_neighbor[j],
+                                              ue_vi_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs[j],
+                                              ui_ve_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs_neighbor[j],
+                                              ue_ve_matrix(i, j));
+                          }
+                    }
+                }
+              else
+                {
+                  if (neighbor->level() == cell->level() &&
+                      neighbor->index() > cell->index())
+                    {
+                      const unsigned int neighbor2 =
+                        cell->neighbor_of_neighbor(face_no);
+
+                      ue_vi_matrix = 0;
+                      ui_ve_matrix = 0;
+                      ue_ve_matrix = 0;
+
+                      fe_v_face.reinit(cell, face_no);
+                      fe_v_face_neighbor.reinit(neighbor, neighbor2);
+
+                      dg.assemble_face_term2(fe_v_face,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix,
+                                             ui_ve_matrix,
+                                             ue_ve_matrix);
+
+                      neighbor->get_dof_indices(dofs_neighbor);
+
+                      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            system_matrix.add(dofs[i],
+                                              dofs_neighbor[j],
+                                              ue_vi_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs[j],
+                                              ui_ve_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs_neighbor[j],
+                                              ue_ve_matrix(i, j));
+                          }
+                    }
+                }
+            }
+        }
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+          system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j));
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        right_hand_side(dofs[i]) += cell_vector(i);
+    }
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::solve(Vector<double> &solution)
+{
+  ReductionControl   solver_control(10000, 1e-4, 1e-3, false, false);
+  SolverRichardson<> solver(solver_control);
+
+  PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+  preconditioner.initialize(system_matrix, fe[0].dofs_per_cell);
+
+  deallog << "  " << right_hand_side.l2_norm() << std::endl;
+  solver.solve(system_matrix, solution, right_hand_side, preconditioner);
+  deallog << "  " << solution.l2_norm() << std::endl;
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::refine_grid()
+{
+  Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+  DerivativeApproximation::approximate_gradient(mapping[0],
+                                                dof_handler,
+                                                solution2,
+                                                gradient_indicator);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  for (unsigned int cell_no = 0; cell != endc; ++cell, ++cell_no)
+    gradient_indicator(cell_no) *=
+      std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+  GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                  gradient_indicator,
+                                                  0.3,
+                                                  0.1);
+
+  triangulation.execute_coarsening_and_refinement();
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::output_results(const unsigned int cycle) const
+{
+  return;
+
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert(cycle < 10, ExcInternalError());
+
+  filename += ".eps";
+  deallog << "Writing grid to <" << filename << ">..." << std::endl;
+
+  GridOut grid_out;
+  // grid_out.write_eps(triangulation, deallog.get_file_stream());
+
+  filename = "sol-";
+  filename += ('0' + cycle);
+  Assert(cycle < 10, ExcInternalError());
+
+  filename += ".gnuplot";
+  deallog << "Writing solution to <" << filename << ">..." << std::endl
+          << std::endl;
+
+  DataOut<dim, DoFHandler<dim>> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution2, "u");
+
+  data_out.build_patches(mapping, 2);
+
+#if false
+  data_out.write_vtk(deallog.get_file_stream());
+#else
+  static unsigned int counter = 0;
+
+  std::ofstream output("test." + std::to_string(dim) + "." +
+                       std::to_string(counter++) + ".vtk");
+  data_out.write_vtk(output);
+#endif
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::run()
+{
+  for (unsigned int cycle = 0; cycle < 6; ++cycle)
+    {
+      deallog << "Cycle " << cycle << ':' << std::endl;
+
+      mesh_generator(triangulation);
+
+
+      deallog << "   Number of active cells:       "
+              << triangulation.n_active_cells() << std::endl;
+
+      setup_system();
+
+      deallog << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+              << std::endl;
+
+      Timer assemble_timer;
+      assemble_system1();
+      solve(solution1);
+
+      system_matrix   = 0;
+      right_hand_side = 0;
+      assemble_timer.reset();
+
+      assemble_timer.start();
+      assemble_system2();
+      solve(solution2);
+
+      solution1 -= solution2;
+      const double difference = solution1.linfty_norm();
+      if (difference > 1e-13)
+        deallog << "solution1 and solution2 differ!!" << std::endl;
+      else
+        deallog << "solution1 and solution2 coincide." << std::endl;
+
+      output_results(cycle);
+
+      break;
+    }
+}
+
+int
+main()
+{
+  try
+    {
+      initlog();
+      deallog.get_file_stream().precision(2);
+
+      // triangle
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<2> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<2>(MappingFE<2>(FE_SimplexP<2>(1)));
+          scratch_data.fe         = hp::FECollection<2>(FE_SimplexDGP<2>(i));
+          scratch_data.quadrature = hp::QCollection<2>(QGaussSimplex<2>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<1>>{
+            hp::QCollection<1>(QGaussSimplex<1>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<2> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                                16);
+          };
+
+          DGMethod<2> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // quadrilateral
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<2> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<2>(MappingQGeneric<2>(1));
+          scratch_data.fe              = hp::FECollection<2>(FE_DGQ<2>(i));
+          scratch_data.quadrature      = hp::QCollection<2>(QGauss<2>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<1>>{
+            hp::QCollection<1>(QGauss<1>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<2> &triangulation) -> void {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(3);
+          };
+
+          DGMethod<2> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // tetrahedron
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_SimplexP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_SimplexDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussSimplex<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGaussSimplex<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                                8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // pyramid
+      for (unsigned int i = 1; i <= 1; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_PyramidP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_PyramidDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussPyramid<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGauss<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_pyramids(triangulation,
+                                                               8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // wedge
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_WedgeP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_WedgeDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussWedge<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGauss<2>(i + 1),
+                               QGauss<2>(i + 1),
+                               QGauss<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_wedges(triangulation, 8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // hexahedron
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingQGeneric<3>(1));
+          scratch_data.fe              = hp::FECollection<3>(FE_DGQ<3>(i));
+          scratch_data.quadrature      = hp::QCollection<3>(QGauss<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGauss<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube(triangulation, 8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.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;
+}
diff --git a/tests/simplex/step-12b.with_simplex_support=on.output b/tests/simplex/step-12b.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..89ed3cd
--- /dev/null
@@ -0,0 +1,89 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 1536
+DEAL::  0.125000
+DEAL::  17.2926
+DEAL::  0.125000
+DEAL::  17.2926
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 3072
+DEAL::  0.125000
+DEAL::  24.3239
+DEAL::  0.125000
+DEAL::  24.3239
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 256
+DEAL::  0.176777
+DEAL::  7.04972
+DEAL::  0.176777
+DEAL::  7.04972
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 576
+DEAL::  0.176777
+DEAL::  10.4560
+DEAL::  0.176777
+DEAL::  10.4560
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       2560
+DEAL::   Number of degrees of freedom: 10240
+DEAL::  0.0360844
+DEAL::  45.1474
+DEAL::  0.0360844
+DEAL::  45.1474
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       2560
+DEAL::   Number of degrees of freedom: 25600
+DEAL::  0.0360844
+DEAL::  69.8629
+DEAL::  0.0360844
+DEAL::  69.8629
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       3072
+DEAL::   Number of degrees of freedom: 15360
+DEAL::  0.0441942
+DEAL::  54.8297
+DEAL::  0.0441942
+DEAL::  54.8297
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 6144
+DEAL::  0.0441942
+DEAL::  34.6857
+DEAL::  0.0441942
+DEAL::  34.6857
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 18432
+DEAL::  0.0441942
+DEAL::  59.5766
+DEAL::  0.0441942
+DEAL::  59.5766
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 4096
+DEAL::  0.0441942
+DEAL::  28.1989
+DEAL::  0.0441942
+DEAL::  28.1989
+DEAL::solution1 and solution2 coincide.
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 13824
+DEAL::  0.0441942
+DEAL::  51.2235
+DEAL::  0.0441942
+DEAL::  51.2235
+DEAL::solution1 and solution2 coincide.
diff --git a/tests/simplex/step-12c.cc b/tests/simplex/step-12c.cc
new file mode 100644 (file)
index 0000000..086c091
--- /dev/null
@@ -0,0 +1,1151 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// a copy of tests/hp/step-12 testing different reference-cell types
+
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/timer.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_pyramid_p.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_wedge_p.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/hp/fe_values.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/precondition_block.h>
+#include <deal.II/lac/solver_richardson.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/derivative_approximation.h>
+
+#include <iostream>
+
+#include "../tests.h"
+
+#include "./simplex_grids.h"
+
+template <int dim>
+struct ScratchData
+{
+  hp::MappingCollection<dim>            mapping;
+  hp::FECollection<dim>                 fe;
+  hp::QCollection<dim>                  quadrature;
+  std::vector<hp::QCollection<dim - 1>> face_quadrature;
+
+  std::function<void(Triangulation<dim> &)> mesh_generator;
+};
+
+
+template <int dim>
+class RHS : public Function<dim>
+{
+public:
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int             component = 0) const;
+};
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+public:
+  virtual void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<double> &          values,
+             const unsigned int             component = 0) const;
+};
+
+
+template <int dim>
+class Beta
+{
+public:
+  Beta()
+  {}
+  void
+  value_list(const std::vector<Point<dim>> &points,
+             std::vector<Point<dim>> &      values) const;
+};
+
+
+template <int dim>
+void
+RHS<dim>::value_list(const std::vector<Point<dim>> &points,
+                     std::vector<double> &          values,
+                     const unsigned int) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i)
+    values[i] = 0;
+}
+
+
+template <int dim>
+void
+Beta<dim>::value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Point<dim>> &      values) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(points.size(), values.size()); ++i)
+    {
+      const Point<dim> &p    = points[i];
+      Point<dim> &      beta = values[i];
+
+      beta(0) = -p(1);
+      beta(1) = p(0);
+
+      if (beta.norm() > 1e-10)
+        beta /= std::sqrt(beta.square());
+    }
+}
+
+
+template <int dim>
+void
+BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+                                std::vector<double> &          values,
+                                const unsigned int) const
+{
+  // Assert(values.size() == points.size(),
+  //       ExcDimensionMismatch(values.size(), points.size()));
+
+  for (unsigned int i = 0; i < std::min(values.size(), points.size()); ++i)
+    {
+      if (points[i](0) < 0.5)
+        values[i] = 1.;
+      else
+        values[i] = 0.;
+    }
+}
+
+
+template <int dim>
+class DGTransportEquation
+{
+public:
+  DGTransportEquation();
+
+  void
+  assemble_cell_term(const hp::FEValues<dim> &fe_v,
+                     FullMatrix<double> &     ui_vi_matrix,
+                     Vector<double> &         cell_vector) const;
+
+  void
+  assemble_boundary_term(const hp::FEFaceValues<dim> &fe_v,
+                         FullMatrix<double> &         ui_vi_matrix,
+                         Vector<double> &             cell_vector) const;
+
+  template <class X, class Y>
+  void
+  assemble_face_term1(const X &           fe_v,
+                      const Y &           fe_v_neighbor,
+                      FullMatrix<double> &ui_vi_matrix,
+                      FullMatrix<double> &ue_vi_matrix) const;
+
+  template <class X, class Y>
+  void
+  assemble_face_term2(const X &           fe_v,
+                      const Y &           fe_v_neighbor,
+                      FullMatrix<double> &ui_vi_matrix,
+                      FullMatrix<double> &ue_vi_matrix,
+                      FullMatrix<double> &ui_ve_matrix,
+                      FullMatrix<double> &ue_ve_matrix) const;
+
+private:
+  const Beta<dim>           beta_function;
+  const RHS<dim>            rhs_function;
+  const BoundaryValues<dim> boundary_function;
+};
+
+
+template <int dim>
+DGTransportEquation<dim>::DGTransportEquation()
+  : beta_function()
+  , rhs_function()
+  , boundary_function()
+{}
+
+
+template <int dim>
+void
+DGTransportEquation<dim>::assemble_cell_term(const hp::FEValues<dim> &fe_v,
+                                             FullMatrix<double> &ui_vi_matrix,
+                                             Vector<double> &cell_vector) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> rhs(fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+  rhs_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                          rhs);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+         ++i)
+      {
+        for (unsigned int j = 0; j < fe_v.get_present_fe_values().dofs_per_cell;
+             ++j)
+          ui_vi_matrix(i, j) -=
+            beta[point] * fe_v.get_present_fe_values().shape_grad(i, point) *
+            fe_v.get_present_fe_values().shape_value(j, point) * JxW[point];
+
+        cell_vector(i) += rhs[point] *
+                          fe_v.get_present_fe_values().shape_value(i, point) *
+                          JxW[point];
+      }
+}
+
+
+template <int dim>
+void
+DGTransportEquation<dim>::assemble_boundary_term(
+  const hp::FEFaceValues<dim> &fe_v,
+  FullMatrix<double> &         ui_vi_matrix,
+  Vector<double> &             cell_vector) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  std::vector<double> g(fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+  boundary_function.value_list(
+    fe_v.get_present_fe_values().get_quadrature_points(), g);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int j = 0;
+               j < fe_v.get_present_fe_values().dofs_per_cell;
+               ++j)
+            ui_vi_matrix(i, j) +=
+              beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+      else
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          cell_vector(i) -= beta_n * g[point] *
+                            fe_v.get_present_fe_values().shape_value(i, point) *
+                            JxW[point];
+    }
+}
+
+
+template <int dim>
+template <class X, class Y>
+void
+DGTransportEquation<dim>::assemble_face_term1(
+  const X &           fe_v,
+  const Y &           fe_v_neighbor,
+  FullMatrix<double> &ui_vi_matrix,
+  FullMatrix<double> &ue_vi_matrix) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int j = 0;
+               j < fe_v.get_present_fe_values().dofs_per_cell;
+               ++j)
+            ui_vi_matrix(i, j) +=
+              beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+      else
+        for (unsigned int i = 0; i < fe_v.get_present_fe_values().dofs_per_cell;
+             ++i)
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            ue_vi_matrix(i, k) +=
+              beta_n *
+              fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+              fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+    }
+}
+
+
+template <int dim>
+template <class X, class Y>
+void
+DGTransportEquation<dim>::assemble_face_term2(
+  const X &           fe_v,
+  const Y &           fe_v_neighbor,
+  FullMatrix<double> &ui_vi_matrix,
+  FullMatrix<double> &ue_vi_matrix,
+  FullMatrix<double> &ui_ve_matrix,
+  FullMatrix<double> &ue_ve_matrix) const
+{
+  const std::vector<double> &JxW =
+    fe_v.get_present_fe_values().get_JxW_values();
+  const std::vector<Tensor<1, dim>> &normals =
+    fe_v.get_present_fe_values().get_normal_vectors();
+
+  std::vector<Point<dim>> beta(
+    fe_v.get_present_fe_values().n_quadrature_points);
+
+  beta_function.value_list(fe_v.get_present_fe_values().get_quadrature_points(),
+                           beta);
+
+  for (unsigned int point = 0;
+       point < fe_v.get_present_fe_values().n_quadrature_points;
+       ++point)
+    {
+      const double beta_n = beta[point] * normals[point];
+      if (beta_n > 0)
+        {
+          for (unsigned int i = 0;
+               i < fe_v.get_present_fe_values().dofs_per_cell;
+               ++i)
+            for (unsigned int j = 0;
+                 j < fe_v.get_present_fe_values().dofs_per_cell;
+                 ++j)
+              ui_vi_matrix(i, j) +=
+                beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+                fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            for (unsigned int j = 0;
+                 j < fe_v.get_present_fe_values().dofs_per_cell;
+                 ++j)
+              ui_ve_matrix(k, j) -=
+                beta_n * fe_v.get_present_fe_values().shape_value(j, point) *
+                fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+                JxW[point];
+        }
+      else
+        {
+          for (unsigned int i = 0;
+               i < fe_v.get_present_fe_values().dofs_per_cell;
+               ++i)
+            for (unsigned int l = 0;
+                 l < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+                 ++l)
+              ue_vi_matrix(i, l) +=
+                beta_n *
+                fe_v_neighbor.get_present_fe_values().shape_value(l, point) *
+                fe_v.get_present_fe_values().shape_value(i, point) * JxW[point];
+
+          for (unsigned int k = 0;
+               k < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+               ++k)
+            for (unsigned int l = 0;
+                 l < fe_v_neighbor.get_present_fe_values().dofs_per_cell;
+                 ++l)
+              ue_ve_matrix(k, l) -=
+                beta_n *
+                fe_v_neighbor.get_present_fe_values().shape_value(l, point) *
+                fe_v_neighbor.get_present_fe_values().shape_value(k, point) *
+                JxW[point];
+        }
+    }
+}
+
+
+template <int dim>
+class DGMethod
+{
+public:
+  DGMethod(const ScratchData<dim> &scratch_data);
+  ~DGMethod();
+
+  void
+  run();
+
+private:
+  void
+  setup_system();
+  void
+  assemble_system2();
+  void
+  solve(Vector<double> &solution);
+  void
+  refine_grid();
+  void
+  output_results(const unsigned int cycle) const;
+
+  Triangulation<dim>               triangulation;
+  const hp::MappingCollection<dim> mapping;
+
+  hp::FECollection<dim> fe;
+  DoFHandler<dim>       dof_handler;
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  const hp::QCollection<dim>                  quadrature;
+  const std::vector<hp::QCollection<dim - 1>> face_quadrature;
+  const hp::QCollection<dim - 1>              face_quadrature_dummy;
+
+  Vector<double> solution1;
+  Vector<double> solution2;
+  Vector<double> right_hand_side;
+
+  const std::function<void(Triangulation<dim> &)> mesh_generator;
+
+  const DGTransportEquation<dim> dg;
+};
+
+
+template <int dim>
+DGMethod<dim>::DGMethod(const ScratchData<dim> &scratch_data)
+  : mapping(scratch_data.mapping)
+  , fe(scratch_data.fe)
+  , dof_handler(triangulation)
+  , quadrature(scratch_data.quadrature)
+  , face_quadrature(scratch_data.face_quadrature)
+  , mesh_generator(scratch_data.mesh_generator)
+  , dg()
+{}
+
+
+template <int dim>
+DGMethod<dim>::~DGMethod()
+{
+  dof_handler.clear();
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(fe);
+
+  sparsity_pattern.reinit(dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          (GeometryInfo<dim>::faces_per_cell *
+                             GeometryInfo<dim>::max_children_per_face +
+                           1) *
+                            fe[0].dofs_per_cell);
+
+  DoFTools::make_flux_sparsity_pattern(dof_handler, sparsity_pattern);
+
+  sparsity_pattern.compress();
+
+  system_matrix.reinit(sparsity_pattern);
+
+  solution1.reinit(dof_handler.n_dofs());
+  solution2.reinit(dof_handler.n_dofs());
+  right_hand_side.reinit(dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void
+DGMethod<dim>::assemble_system2()
+{
+  const unsigned int dofs_per_cell = dof_handler.get_fe(0).dofs_per_cell;
+  std::vector<types::global_dof_index> dofs(dofs_per_cell);
+  std::vector<types::global_dof_index> dofs_neighbor(dofs_per_cell);
+
+  const UpdateFlags update_flags = update_values | update_gradients |
+                                   update_quadrature_points | update_JxW_values;
+
+  const UpdateFlags face_update_flags =
+    update_values | update_quadrature_points | update_JxW_values |
+    update_normal_vectors;
+
+  const UpdateFlags neighbor_face_update_flags = update_values;
+
+  hp::FEValues<dim>        fe_v(mapping, fe, quadrature, update_flags);
+  hp::FEFaceValues<dim>    fe_v_face(mapping,
+                                  fe,
+                                  face_quadrature,
+                                  face_update_flags);
+  hp::FESubfaceValues<dim> fe_v_subface(mapping,
+                                        fe,
+                                        face_quadrature_dummy,
+                                        face_update_flags);
+  hp::FEFaceValues<dim>    fe_v_face_neighbor(mapping,
+                                           fe,
+                                           face_quadrature,
+                                           neighbor_face_update_flags);
+
+
+  FullMatrix<double> ui_vi_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ue_vi_matrix(dofs_per_cell, dofs_per_cell);
+
+  FullMatrix<double> ui_ve_matrix(dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> ue_ve_matrix(dofs_per_cell, dofs_per_cell);
+
+  Vector<double> cell_vector(dofs_per_cell);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  for (; cell != endc; ++cell)
+    {
+      ui_vi_matrix = 0;
+      cell_vector  = 0;
+
+      fe_v.reinit(cell);
+
+      dg.assemble_cell_term(fe_v, ui_vi_matrix, cell_vector);
+
+      cell->get_dof_indices(dofs);
+
+      for (const unsigned int face_no : cell->face_indices())
+        {
+          typename DoFHandler<dim>::face_iterator face = cell->face(face_no);
+
+          if (face->at_boundary())
+            {
+              fe_v_face.reinit(cell, face_no);
+
+              dg.assemble_boundary_term(fe_v_face, ui_vi_matrix, cell_vector);
+            }
+          else
+            {
+              Assert(cell->neighbor(face_no).state() == IteratorState::valid,
+                     ExcInternalError());
+              typename DoFHandler<dim>::cell_iterator neighbor =
+                cell->neighbor(face_no);
+              if (face->has_children())
+                {
+                  const unsigned int neighbor2 =
+                    cell->neighbor_of_neighbor(face_no);
+
+                  for (unsigned int subface_no = 0;
+                       subface_no < face->n_children();
+                       ++subface_no)
+                    {
+                      typename DoFHandler<dim>::cell_iterator neighbor_child =
+                        cell->neighbor_child_on_subface(face_no, subface_no);
+                      Assert(neighbor_child->face(neighbor2) ==
+                               face->child(subface_no),
+                             ExcInternalError());
+                      Assert(!neighbor_child->has_children(),
+                             ExcInternalError());
+
+                      ue_vi_matrix = 0;
+                      ui_ve_matrix = 0;
+                      ue_ve_matrix = 0;
+
+                      fe_v_subface.reinit(cell, face_no, subface_no);
+                      fe_v_face_neighbor.reinit(neighbor_child, neighbor2);
+
+                      dg.assemble_face_term2(fe_v_subface,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix,
+                                             ui_ve_matrix,
+                                             ue_ve_matrix);
+
+                      neighbor_child->get_dof_indices(dofs_neighbor);
+
+                      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            system_matrix.add(dofs[i],
+                                              dofs_neighbor[j],
+                                              ue_vi_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs[j],
+                                              ui_ve_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs_neighbor[j],
+                                              ue_ve_matrix(i, j));
+                          }
+                    }
+                }
+              else
+                {
+                  if (neighbor->level() == cell->level() &&
+                      neighbor->index() > cell->index())
+                    {
+                      const unsigned int neighbor2 =
+                        cell->neighbor_of_neighbor(face_no);
+
+                      ue_vi_matrix = 0;
+                      ui_ve_matrix = 0;
+                      ue_ve_matrix = 0;
+
+                      fe_v_face.reinit(cell, face_no);
+                      fe_v_face_neighbor.reinit(neighbor, neighbor2);
+
+                      dg.assemble_face_term2(fe_v_face,
+                                             fe_v_face_neighbor,
+                                             ui_vi_matrix,
+                                             ue_vi_matrix,
+                                             ui_ve_matrix,
+                                             ue_ve_matrix);
+
+                      neighbor->get_dof_indices(dofs_neighbor);
+
+                      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+                        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+                          {
+                            system_matrix.add(dofs[i],
+                                              dofs_neighbor[j],
+                                              ue_vi_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs[j],
+                                              ui_ve_matrix(i, j));
+                            system_matrix.add(dofs_neighbor[i],
+                                              dofs_neighbor[j],
+                                              ue_ve_matrix(i, j));
+                          }
+                    }
+                }
+            }
+        }
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int j = 0; j < dofs_per_cell; ++j)
+          system_matrix.add(dofs[i], dofs[j], ui_vi_matrix(i, j));
+
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        right_hand_side(dofs[i]) += cell_vector(i);
+    }
+}
+
+
+
+template <int dim>
+class Operator
+{
+public:
+  using VectorType = Vector<double>;
+  using number     = double;
+
+  Operator(const MatrixFree<dim, double> &matrix_free)
+    : matrix_free(matrix_free)
+  {}
+
+  void
+  initialize_dof_vector(VectorType &vec)
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  void
+  rhs(VectorType &vec) const
+  {
+    const int dummy = 0;
+
+    matrix_free.template loop<VectorType, int>(
+      [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
+        FEEvaluation<dim, -1, 0, 1, double> phi(data, cell_range);
+        for (unsigned int cell = cell_range.first; cell < cell_range.second;
+             ++cell)
+          {
+            phi.reinit(cell);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_value(0.0, q);
+
+            phi.integrate_scatter(true, false, dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        (void)data;
+        (void)dst;
+        (void)src;
+        (void)face_range;
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, double> phi(data, face_range);
+        for (unsigned int face = face_range.first; face < face_range.second;
+             ++face)
+          {
+            phi.reinit(face);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              {
+                const auto beta_n =
+                  beta(phi.quadrature_point(q)) * phi.get_normal_vector(q);
+                const auto beta_n_m = (-std::abs(beta_n) + beta_n) * 0.5;
+                phi.submit_value(-beta_n_m *
+                                   boundary_values(phi.quadrature_point(q)),
+                                 q);
+              }
+
+            phi.integrate_scatter(true, false, dst);
+          }
+      },
+      vec,
+      dummy,
+      true);
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    matrix_free.template loop<VectorType, VectorType>(
+      [&](const auto &data, auto &dst, const auto &src, const auto cell_range) {
+        FEEvaluation<dim, -1, 0, 1, double> phi(data, cell_range);
+        for (unsigned int cell = cell_range.first; cell < cell_range.second;
+             ++cell)
+          {
+            phi.reinit(cell);
+            phi.gather_evaluate(src, true, false);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              phi.submit_gradient(-beta(phi.quadrature_point(q)) *
+                                    phi.get_value(q),
+                                  q);
+            phi.integrate_scatter(false, true, dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, double> phi_m(data, face_range, true);
+        FEFaceEvaluation<dim, -1, 0, 1, double> phi_p(data, face_range, false);
+        for (unsigned int cell = face_range.first; cell < face_range.second;
+             ++cell)
+          {
+            phi_m.reinit(cell);
+            phi_m.gather_evaluate(src, true, false);
+            phi_p.reinit(cell);
+            phi_p.gather_evaluate(src, true, false);
+            for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
+              {
+                const auto beta_n =
+                  beta(phi_m.quadrature_point(q)) * phi_m.get_normal_vector(q);
+
+                const auto u_m = phi_m.get_value(q);
+                const auto u_p = phi_p.get_value(q);
+                const auto temp =
+                  0.5 * ((u_m + u_p) * beta_n + (u_m - u_p) * std::abs(beta_n));
+
+                phi_m.submit_value(+temp, q);
+                phi_p.submit_value(-temp, q);
+              }
+
+            phi_m.integrate_scatter(true, false, dst);
+            phi_p.integrate_scatter(true, false, dst);
+          }
+      },
+      [&](const auto &data, auto &dst, const auto &src, const auto face_range) {
+        FEFaceEvaluation<dim, -1, 0, 1, double> phi(data, face_range);
+        for (unsigned int cell = face_range.first; cell < face_range.second;
+             ++cell)
+          {
+            phi.reinit(cell);
+            phi.gather_evaluate(src, true, false);
+            for (unsigned int q = 0; q < phi.n_q_points; ++q)
+              {
+                const auto beta_n =
+                  beta(phi.quadrature_point(q)) * phi.get_normal_vector(q);
+                const auto beta_n_p = (std::abs(beta_n) + beta_n) * 0.5;
+                phi.submit_value(beta_n_p * phi.get_value(q), q);
+              }
+
+            phi.integrate_scatter(true, false, dst);
+          }
+      },
+      dst,
+      src,
+      true);
+  }
+
+  Tensor<1, dim, VectorizedArray<double>>
+  beta(const Point<dim, VectorizedArray<double>> &points) const
+  {
+    Tensor<1, dim, VectorizedArray<double>> betas;
+
+    for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+      {
+        Point<dim, double>     p;
+        Tensor<1, dim, double> beta;
+
+        for (int d = 0; d < dim; ++d)
+          p[d] = points[d][v];
+
+        beta[0] = -p[1];
+        beta[1] = +p[0];
+
+        if (beta.norm() > 1e-10)
+          beta /= beta.norm();
+
+        for (int d = 0; d < dim; ++d)
+          betas[d][v] = beta[d];
+      }
+
+    return betas;
+  }
+
+  VectorizedArray<double>
+  boundary_values(const Point<dim, VectorizedArray<double>> &points) const
+  {
+    VectorizedArray<double> betas;
+
+    for (unsigned int v = 0; v < VectorizedArray<double>::size(); ++v)
+      betas[v] = (points[0][v] < 0.5) ? 1.0 : 0.0;
+
+    return betas;
+  }
+
+private:
+  const MatrixFree<dim, double> &matrix_free;
+};
+
+
+template <int dim>
+void
+DGMethod<dim>::solve(Vector<double> &solution)
+{
+  ReductionControl   solver_control(10000, 1e-4, 1e-3, false, false);
+  SolverRichardson<> solver(solver_control);
+
+  PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+  preconditioner.initialize(system_matrix, fe[0].dofs_per_cell);
+
+  typename MatrixFree<dim, double>::AdditionalData additional_data;
+  additional_data.mapping_update_flags =
+    update_gradients | update_values | update_quadrature_points;
+  additional_data.mapping_update_flags_inner_faces =
+    update_gradients | update_values | update_quadrature_points;
+  additional_data.mapping_update_flags_boundary_faces =
+    update_gradients | update_values | update_quadrature_points;
+  additional_data.tasks_parallel_scheme =
+    MatrixFree<dim, double>::AdditionalData::none;
+
+  AffineConstraints<double> constraints;
+  constraints.close();
+
+  MatrixFree<dim, double> matrix_free;
+  matrix_free.reinit(
+    mapping, dof_handler, constraints, quadrature, additional_data);
+
+  Operator<dim> op(matrix_free);
+
+  op.initialize_dof_vector(solution);
+  op.initialize_dof_vector(right_hand_side);
+
+  op.rhs(right_hand_side);
+
+  deallog << "  " << right_hand_side.l2_norm() << std::endl;
+  solver.solve(op, solution, right_hand_side, preconditioner);
+  deallog << "  " << solution.l2_norm() << std::endl;
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::refine_grid()
+{
+  Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+  DerivativeApproximation::approximate_gradient(mapping[0],
+                                                dof_handler,
+                                                solution2,
+                                                gradient_indicator);
+
+  typename DoFHandler<dim>::active_cell_iterator cell =
+                                                   dof_handler.begin_active(),
+                                                 endc = dof_handler.end();
+  for (unsigned int cell_no = 0; cell != endc; ++cell, ++cell_no)
+    gradient_indicator(cell_no) *=
+      std::pow(cell->diameter(), 1 + 1.0 * dim / 2);
+
+  GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                  gradient_indicator,
+                                                  0.3,
+                                                  0.1);
+
+  triangulation.execute_coarsening_and_refinement();
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::output_results(const unsigned int cycle) const
+{
+  return;
+
+  std::string filename = "grid-";
+  filename += ('0' + cycle);
+  Assert(cycle < 10, ExcInternalError());
+
+  filename += ".eps";
+  deallog << "Writing grid to <" << filename << ">..." << std::endl;
+
+  GridOut grid_out;
+  // grid_out.write_eps(triangulation, deallog.get_file_stream());
+
+  filename = "sol-";
+  filename += ('0' + cycle);
+  Assert(cycle < 10, ExcInternalError());
+
+  filename += ".gnuplot";
+  deallog << "Writing solution to <" << filename << ">..." << std::endl
+          << std::endl;
+
+  DataOut<dim, DoFHandler<dim>> data_out;
+  data_out.attach_dof_handler(dof_handler);
+  data_out.add_data_vector(solution1, "u");
+
+  data_out.build_patches(mapping, 2);
+
+#if false
+  data_out.write_vtk(deallog.get_file_stream());
+#else
+  static unsigned int counter = 0;
+
+  std::ofstream output("test." + std::to_string(dim) + "." +
+                       std::to_string(counter++) + ".vtk");
+  data_out.write_vtk(output);
+#endif
+}
+
+
+template <int dim>
+void
+DGMethod<dim>::run()
+{
+  deallog << "Cycle " << 0 << ':' << std::endl;
+  mesh_generator(triangulation);
+  deallog << "   Number of active cells:       "
+          << triangulation.n_active_cells() << std::endl;
+
+  setup_system();
+  assemble_system2();
+
+  deallog << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+          << std::endl;
+  solve(solution1);
+
+  output_results(0);
+}
+
+int
+main()
+{
+  try
+    {
+      initlog();
+      deallog.get_file_stream().precision(2);
+
+      // triangle
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<2> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<2>(MappingFE<2>(FE_SimplexP<2>(1)));
+          scratch_data.fe         = hp::FECollection<2>(FE_SimplexDGP<2>(i));
+          scratch_data.quadrature = hp::QCollection<2>(QGaussSimplex<2>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<1>>{
+            hp::QCollection<1>(QGaussSimplex<1>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<2> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                                16);
+          };
+
+          DGMethod<2> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // quadrilateral
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<2> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<2>(MappingQGeneric<2>(1));
+          scratch_data.fe              = hp::FECollection<2>(FE_DGQ<2>(i));
+          scratch_data.quadrature      = hp::QCollection<2>(QGauss<2>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<1>>{
+            hp::QCollection<1>(QGauss<1>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<2> &triangulation) -> void {
+            GridGenerator::hyper_cube(triangulation);
+            triangulation.refine_global(3);
+          };
+
+          DGMethod<2> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // tetrahedron
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_SimplexP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_SimplexDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussSimplex<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGaussSimplex<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                                8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // pyramid
+      for (unsigned int i = 1; i <= 1; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_PyramidP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_PyramidDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussPyramid<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGauss<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_pyramids(triangulation,
+                                                               8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // wedge
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingFE<3>(FE_WedgeP<3>(1)));
+          scratch_data.fe         = hp::FECollection<3>(FE_WedgeDGP<3>(i));
+          scratch_data.quadrature = hp::QCollection<3>(QGaussWedge<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGaussSimplex<2>(i + 1),
+                               QGaussSimplex<2>(i + 1),
+                               QGauss<2>(i + 1),
+                               QGauss<2>(i + 1),
+                               QGauss<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube_with_wedges(triangulation, 8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.run();
+        }
+
+      // hexahedron
+      for (unsigned int i = 1; i <= 2; ++i)
+        {
+          ScratchData<3> scratch_data;
+
+          scratch_data.mapping =
+            hp::MappingCollection<3>(MappingQGeneric<3>(1));
+          scratch_data.fe              = hp::FECollection<3>(FE_DGQ<3>(i));
+          scratch_data.quadrature      = hp::QCollection<3>(QGauss<3>(i + 1));
+          scratch_data.face_quadrature = std::vector<hp::QCollection<2>>{
+            hp::QCollection<2>(QGauss<2>(i + 1))};
+          scratch_data.mesh_generator =
+            [](Triangulation<3> &triangulation) -> void {
+            GridGenerator::subdivided_hyper_cube(triangulation, 8);
+          };
+
+          DGMethod<3> dgmethod(scratch_data);
+          dgmethod.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;
+}
diff --git a/tests/simplex/step-12c.with_simplex_support=on.output b/tests/simplex/step-12c.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..b8e96aa
--- /dev/null
@@ -0,0 +1,56 @@
+
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 1536
+DEAL::  0.125000
+DEAL::  17.2926
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 3072
+DEAL::  0.125000
+DEAL::  24.3239
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 256
+DEAL::  0.176777
+DEAL::  7.04972
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       64
+DEAL::   Number of degrees of freedom: 576
+DEAL::  0.176777
+DEAL::  10.4560
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       2560
+DEAL::   Number of degrees of freedom: 10240
+DEAL::  0.0360844
+DEAL::  45.1474
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       2560
+DEAL::   Number of degrees of freedom: 25600
+DEAL::  0.0360844
+DEAL::  69.8629
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       3072
+DEAL::   Number of degrees of freedom: 15360
+DEAL::  0.0441942
+DEAL::  54.8297
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 6144
+DEAL::  0.0441942
+DEAL::  34.6857
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       1024
+DEAL::   Number of degrees of freedom: 18432
+DEAL::  0.0441942
+DEAL::  59.5766
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 4096
+DEAL::  0.0441942
+DEAL::  28.1989
+DEAL::Cycle 0:
+DEAL::   Number of active cells:       512
+DEAL::   Number of degrees of freedom: 13824
+DEAL::  0.0441942
+DEAL::  51.2235

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.