]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add mini simplex application tests 10758/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 28 Jul 2020 09:59:26 +0000 (11:59 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 3 Aug 2020 12:35:18 +0000 (14:35 +0200)
13 files changed:
source/base/qprojector.cc
tests/simplex/data_out_write_vtk_01.cc
tests/simplex/poisson_01.cc [new file with mode: 0644]
tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/poisson_02.cc [new file with mode: 0644]
tests/simplex/poisson_02.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/step-12.cc [new file with mode: 0644]
tests/simplex/step-12.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/step-18.cc [new file with mode: 0644]
tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output [new file with mode: 0644]
tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output [new file with mode: 0644]

index d9c87099c0dd082ab7ec0ecc7722e0c200d883bc..e3e7b4e1337122a8efe715d5f007cf8de72dd661 100644 (file)
@@ -906,8 +906,11 @@ QProjector<2>::project_to_all_subfaces(
   const ReferenceCell::Type reference_cell_type,
   const SubQuadrature &     quadrature)
 {
+  if (reference_cell_type == ReferenceCell::Type::Tri ||
+      reference_cell_type == ReferenceCell::Type::Tet)
+    return Quadrature<2>(); // nothing to do
+
   Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
-  (void)reference_cell_type;
 
   const unsigned int dim = 2;
 
@@ -964,8 +967,11 @@ QProjector<3>::project_to_all_subfaces(
   const ReferenceCell::Type reference_cell_type,
   const SubQuadrature &     quadrature)
 {
+  if (reference_cell_type == ReferenceCell::Type::Tri ||
+      reference_cell_type == ReferenceCell::Type::Tet)
+    return Quadrature<3>(); // nothing to do
+
   Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
-  (void)reference_cell_type;
 
   const unsigned int dim         = 3;
   SubQuadrature      q_reflected = reflect(quadrature);
index 027654e519a76d08eac59ff0d036d8f0e775bebe..f830e658a20cd8e318c95d1081b057adbf327eed 100644 (file)
@@ -56,7 +56,7 @@ void
 test(const FiniteElement<dim, spacedim> &fe, const unsigned int n_components)
 {
   Triangulation<dim, spacedim> tria;
-  Simplex::GridGenerator::subdivided_hyper_cube(tria, dim == 2 ? 4 : 2);
+  GridGenerator::subdivided_hyper_cube_with_simplices(tria, dim == 2 ? 4 : 2);
 
   DoFHandler<dim> dof_handler(tria);
 
diff --git a/tests/simplex/poisson_01.cc b/tests/simplex/poisson_01.cc
new file mode 100644 (file)
index 0000000..8441033
--- /dev/null
@@ -0,0 +1,455 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem on a tet mesh and on a quad mesh with the same number
+// of subdivisions.
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/shared_tria.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.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/sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+struct Parameters
+{
+  unsigned int degree = 2;
+
+
+  // GridGenerator
+  bool                      use_grid_generator = true;
+  std::vector<unsigned int> repetitions;
+  Point<dim>                p1;
+  Point<dim>                p2;
+
+  // GridIn
+  std::string file_name_in = "";
+
+  // GridOut
+  std::string file_name_out = "";
+};
+
+template <int dim, int spacedim>
+MPI_Comm
+get_communicator(const Triangulation<dim, spacedim> &tria)
+{
+  if (auto tria_ =
+        dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(&tria))
+    return tria_->get_communicator();
+
+  return MPI_COMM_SELF;
+}
+
+template <int dim, int spacedim = dim>
+void
+test(const Triangulation<dim, spacedim> &tria,
+     const FiniteElement<dim, spacedim> &fe,
+     const Quadrature<dim> &             quad,
+     const Quadrature<dim - 1> &         face_quad,
+     const Mapping<dim, spacedim> &      mapping,
+     const double                        r_boundary)
+{
+  std::string label =
+    (dynamic_cast<const parallel::shared::Triangulation<dim, spacedim> *>(
+       &tria) ?
+       "parallel::shared::Triangulation" :
+       (dynamic_cast<
+          const parallel::fullydistributed::Triangulation<dim, spacedim> *>(
+          &tria) ?
+          "parallel::fullydistributed::Triangulation" :
+          (dynamic_cast<
+             const parallel::distributed::Triangulation<dim, spacedim> *>(
+             &tria) ?
+             "parallel::distributed::Triangulation" :
+             "Triangulation")));
+
+  deallog << "   on " << label << std::endl;
+
+
+  for (const auto &cell : tria.active_cell_iterators())
+    for (const auto &face : cell->face_iterators())
+      if (face->at_boundary() &&
+          (std::abs(face->center()[0] - r_boundary) < 1e-6))
+        face->set_boundary_id(1);
+      else if (face->at_boundary() && face->center()[1] == 0.0)
+        face->set_boundary_id(2);
+      else if (face->at_boundary() && face->center()[1] == 1.0)
+        face->set_boundary_id(2);
+      else if (dim == 3 && face->at_boundary() && face->center()[2] == 0.0)
+        face->set_boundary_id(2);
+      else if (dim == 3 && face->at_boundary() && face->center()[2] == 1.0)
+        face->set_boundary_id(2);
+      else if (face->at_boundary())
+        face->set_boundary_id(0);
+
+
+  DoFHandler<dim, spacedim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<double> constraint_matrix;
+  DoFTools::make_zero_boundary_constraints(dof_handler, 0, constraint_matrix);
+  constraint_matrix.close();
+
+  // constraint_matrix.print(std::cout);
+
+  const MPI_Comm comm = get_communicator(dof_handler.get_triangulation());
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+  using VectorType = LinearAlgebra::distributed::Vector<double>;
+  TrilinosWrappers::SparseMatrix system_matrix;
+  VectorType                     solution;
+  VectorType                     system_rhs;
+#else
+  using VectorType = Vector<double>;
+  SparsityPattern        sparsity_pattern;
+  SparseMatrix<double>   system_matrix;
+  VectorType             solution(dof_handler.n_dofs());
+  VectorType             system_rhs(dof_handler.n_dofs());
+#endif
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+  TrilinosWrappers::SparsityPattern dsp(dof_handler.locally_owned_dofs(), comm);
+#else
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+#endif
+  DoFTools::make_sparsity_pattern(dof_handler, dsp, constraint_matrix);
+#ifdef DEAL_II_WITH_TRILINOS
+  dsp.compress();
+  system_matrix.reinit(dsp);
+
+
+  solution.reinit(dof_handler.locally_owned_dofs(),
+                  locally_relevant_dofs,
+                  comm);
+  system_rhs.reinit(dof_handler.locally_owned_dofs(),
+                    locally_relevant_dofs,
+                    comm);
+#else
+  sparsity_pattern.copy_from(dsp);
+  system_matrix.reinit(sparsity_pattern);
+#endif
+
+  const UpdateFlags flag = update_JxW_values | update_values |
+                           update_gradients | update_quadrature_points;
+  FEValues<dim, spacedim> fe_values(mapping, fe, quad, flag);
+
+  std::shared_ptr<FEFaceValues<dim, spacedim>> fe_face_values;
+
+  fe_face_values.reset(
+    new FEFaceValues<dim, spacedim>(mapping, fe, face_quad, flag));
+
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points    = quad.size();
+
+  std::vector<types::global_dof_index> dof_indices(dofs_per_cell);
+  FullMatrix<double> cell_matrix(dofs_per_cell, 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_matrix = 0;
+      cell_rhs    = 0;
+
+      for (unsigned int q_index = 0; q_index < n_q_points; ++q_index)
+        for (unsigned int i = 0; i < dofs_per_cell; ++i)
+          {
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              cell_matrix(i, j) +=
+                (fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
+                 fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
+                 fe_values.JxW(q_index));           // dx
+            cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
+                            1.0 *                               // 1.0
+                            fe_values.JxW(q_index));            // dx
+          }
+
+      if (fe_face_values)
+        for (const auto &face : cell->face_iterators())
+          if (face->at_boundary() && (face->boundary_id() == 1))
+            {
+              fe_face_values->reinit(cell, face);
+              for (unsigned int q = 0; q < face_quad.size(); ++q)
+                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);
+
+      constraint_matrix.distribute_local_to_global(
+        cell_matrix, cell_rhs, dof_indices, system_matrix, system_rhs);
+    }
+
+  system_matrix.compress(VectorOperation::add);
+  system_rhs.compress(VectorOperation::add);
+
+  SolverControl        solver_control(1000, 1e-12);
+  SolverCG<VectorType> solver(solver_control);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+  deallog << "   with " << solver_control.last_step()
+          << " CG iterations needed to obtain convergence" << std::endl;
+
+  // system_rhs.print(std::cout);
+  // solution.print(std::cout);
+
+  bool hex_mesh = true;
+
+  for (const auto &cell : tria.active_cell_iterators())
+    hex_mesh &= (cell->n_vertices() == GeometryInfo<dim>::vertices_per_cell);
+
+  deallog << std::endl;
+}
+
+template <int dim, int spacedim = dim>
+void
+test_tet(const MPI_Comm &comm, const Parameters<dim> &params)
+{
+  const unsigned int tria_type = 2;
+
+  // 1) Create triangulation...
+  Triangulation<dim, spacedim> *tria;
+
+  // a) serial triangulation
+  Triangulation<dim, spacedim> tr_1;
+
+  // b) shared triangulation (with artificial cells)
+  parallel::shared::Triangulation<dim> tr_2(
+    MPI_COMM_WORLD,
+    ::Triangulation<dim>::none,
+    true,
+    parallel::shared::Triangulation<dim>::partition_custom_signal);
+
+  tr_2.signals.create.connect([&]() {
+    GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+                                       tr_2);
+  });
+
+  // c) distributed triangulation
+  parallel::fullydistributed::Triangulation<dim> tr_3(comm);
+
+
+  // ... choose the right triangulation
+  if (tria_type == 0 || tria_type == 2)
+    tria = &tr_1;
+  else if (tria_type == 1)
+    tria = &tr_2;
+
+  // ... create triangulation
+  if (params.use_grid_generator)
+    {
+      // ...via Simplex::GridGenerator
+      GridGenerator::subdivided_hyper_rectangle_with_simplices(
+        *tria, params.repetitions, params.p1, params.p2, false);
+    }
+  else
+    {
+      // ...via GridIn
+      GridIn<dim, spacedim> grid_in;
+      grid_in.attach_triangulation(*tria);
+      std::ifstream input_file(params.file_name_in);
+      grid_in.read_ucd(input_file);
+      // std::ifstream input_file("test_tet_geometry.unv");
+      // grid_in.read_unv(input_file);
+    }
+
+  // ... partition serial triangulation and create distributed triangulation
+  if (tria_type == 0 || tria_type == 2)
+    {
+      GridTools::partition_triangulation(Utilities::MPI::n_mpi_processes(comm),
+                                         tr_1);
+
+      auto construction_data = TriangulationDescription::Utilities::
+        create_description_from_triangulation(tr_1, comm);
+
+      tr_3.create_triangulation(construction_data);
+
+      tria = &tr_3;
+    }
+
+  // 2) Output generated triangulation via GridOut
+  GridOut       grid_out;
+  std::ofstream out(params.file_name_out + "." +
+                    std::to_string(Utilities::MPI::this_mpi_process(comm)) +
+                    ".vtk");
+  grid_out.write_vtk(*tria, out);
+
+  // 3) Select components
+  Simplex::FE_P<dim> fe(params.degree);
+
+  Simplex::PGauss<dim> quad(dim == 2 ? (params.degree == 1 ? 3 : 7) :
+                                       (params.degree == 1 ? 4 : 10));
+
+  Simplex::PGauss<dim - 1> face_quad(dim == 2 ? (params.degree == 1 ? 2 : 3) :
+                                                (params.degree == 1 ? 3 : 7));
+
+  Simplex::FE_P<dim> fe_mapping(1);
+  MappingFE<dim>     mapping(fe_mapping);
+
+  // 4) Perform test (independent of mesh type)
+  test(*tria, fe, quad, face_quad, mapping, params.p2[0]);
+}
+
+template <int dim, int spacedim = dim>
+void
+test_hex(const MPI_Comm &comm, const Parameters<dim> &params)
+{
+  // 1) Create triangulation...
+  parallel::distributed::Triangulation<dim, spacedim> tria(comm);
+
+  if (params.use_grid_generator)
+    {
+      // ...via GridGenerator
+      GridGenerator::subdivided_hyper_rectangle(
+        tria, params.repetitions, params.p1, params.p2, false);
+    }
+  else
+    {
+      // ...via GridIn
+      GridIn<dim, spacedim> grid_in;
+      grid_in.attach_triangulation(tria);
+      std::ifstream input_file(params.file_name_in);
+      grid_in.read_ucd(input_file);
+    }
+
+  // 2) Output generated triangulation via GridOut
+  GridOut       grid_out;
+  std::ofstream out(params.file_name_out + "." +
+                    std::to_string(Utilities::MPI::this_mpi_process(comm)) +
+                    ".vtk");
+  grid_out.write_vtk(tria, out);
+
+  // 3) Select components
+  FE_Q<dim> fe(params.degree);
+
+  QGauss<dim> quad(params.degree + 1);
+
+  QGauss<dim - 1> quad_face(params.degree + 1);
+
+  MappingQ<dim, spacedim> mapping(1);
+
+  // 4) Perform test (independent of mesh type)
+  test(tria, fe, quad, quad_face, mapping, params.p2[0]);
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+  initlog();
+
+  const MPI_Comm comm = MPI_COMM_WORLD;
+
+  // 2D
+  {
+    Parameters<2> params;
+    params.use_grid_generator = true;
+    params.repetitions        = std::vector<unsigned int>{10, 10};
+
+    // test TRI
+    {
+      deallog << "Solve problem on TRI mesh:" << std::endl;
+
+      params.file_name_out = "mesh-tri";
+      params.p1            = Point<2>(0, 0);
+      params.p2            = Point<2>(1, 1);
+      test_tet(comm, params);
+    }
+
+    // test QUAD
+    {
+      deallog << "Solve problem on QUAD mesh:" << std::endl;
+
+      params.file_name_out = "mesh-quad";
+      params.p1            = Point<2>(1.1, 0); // shift to the right for
+      params.p2            = Point<2>(2.1, 1); // visualization purposes
+      test_hex(comm, params);
+    }
+  }
+
+  // 3D
+  {
+    Parameters<3> params;
+    params.use_grid_generator = true;
+    params.repetitions        = std::vector<unsigned int>{10, 10, 10};
+
+    // test TET
+    {
+      deallog << "Solve problem on TET mesh:" << std::endl;
+
+      params.file_name_out = "mesh-tet";
+      params.p1            = Point<3>(0, 0, 0);
+      params.p2            = Point<3>(1, 1, 1);
+      test_tet(comm, params);
+    }
+
+    // test HEX
+    {
+      deallog << "Solve problem on HEX mesh:" << std::endl;
+
+      params.file_name_out = "mesh-hex";
+      params.p1            = Point<3>(1.1, 0, 0);
+      params.p2            = Point<3>(2.1, 1, 1);
+      test_hex(comm, params);
+    }
+  }
+}
diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=false.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..02c1cc8
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL:0::Solve problem on TRI mesh:
+DEAL:0::   on parallel::fullydistributed::Triangulation
+DEAL:0:cg::Starting value 0.245798
+DEAL:0:cg::Convergence step 114 value 6.62438e-13
+DEAL:0::   with 114 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on QUAD mesh:
+DEAL:0::   on parallel::distributed::Triangulation
+DEAL:0:cg::Starting value 0.244628
+DEAL:0:cg::Convergence step 98 value 8.55703e-13
+DEAL:0::   with 98 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on TET mesh:
+DEAL:0::   on parallel::fullydistributed::Triangulation
+DEAL:0:cg::Starting value 0.0607616
+DEAL:0:cg::Convergence step 156 value 9.07624e-13
+DEAL:0::   with 156 CG iterations needed to obtain convergence
+DEAL:0::
+DEAL:0::Solve problem on HEX mesh:
+DEAL:0::   on parallel::distributed::Triangulation
+DEAL:0:cg::Starting value 0.0573704
+DEAL:0:cg::Convergence step 134 value 8.23917e-13
+DEAL:0::   with 134 CG iterations needed to obtain convergence
+DEAL:0::
diff --git a/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=1.with_trilinos=true.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..a45012d
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL::Solve problem on TRI mesh:
+DEAL::   on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.245798
+DEAL:cg::Convergence step 114 value 6.95996e-13
+DEAL::   with 114 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on QUAD mesh:
+DEAL::   on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.244628
+DEAL:cg::Convergence step 98 value 8.55703e-13
+DEAL::   with 98 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on TET mesh:
+DEAL::   on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0607616
+DEAL:cg::Convergence step 156 value 9.07635e-13
+DEAL::   with 156 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on HEX mesh:
+DEAL::   on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.0573704
+DEAL:cg::Convergence step 134 value 8.23917e-13
+DEAL::   with 134 CG iterations needed to obtain convergence
+DEAL::
diff --git a/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output b/tests/simplex/poisson_01.mpirun=4.with_trilinos=true.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..0e41dea
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL::Solve problem on TRI mesh:
+DEAL::   on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.245798
+DEAL:cg::Convergence step 114 value 6.45793e-13
+DEAL::   with 114 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on QUAD mesh:
+DEAL::   on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.244628
+DEAL:cg::Convergence step 98 value 8.55703e-13
+DEAL::   with 98 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on TET mesh:
+DEAL::   on parallel::fullydistributed::Triangulation
+DEAL:cg::Starting value 0.0607616
+DEAL:cg::Convergence step 156 value 9.07640e-13
+DEAL::   with 156 CG iterations needed to obtain convergence
+DEAL::
+DEAL::Solve problem on HEX mesh:
+DEAL::   on parallel::distributed::Triangulation
+DEAL:cg::Starting value 0.0573704
+DEAL:cg::Convergence step 134 value 8.23917e-13
+DEAL::   with 134 CG iterations needed to obtain convergence
+DEAL::
diff --git a/tests/simplex/poisson_02.cc b/tests/simplex/poisson_02.cc
new file mode 100644 (file)
index 0000000..03110d8
--- /dev/null
@@ -0,0 +1,670 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Solve Poisson problem on a tet mesh with DG.
+
+#include <deal.II/base/convergence_table.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.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_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_fe.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+//#define HEX
+
+using namespace dealii;
+
+template <int dim>
+struct ScratchData
+{
+  ScratchData(const Mapping<dim> &       mapping,
+              const FiniteElement<dim> & fe,
+              const Quadrature<dim> &    quad,
+              const Quadrature<dim - 1> &quad_face,
+              const UpdateFlags          update_flags = update_values |
+                                               update_gradients |
+                                               update_quadrature_points |
+                                               update_JxW_values,
+              const UpdateFlags interface_update_flags =
+                update_values | update_gradients | update_quadrature_points |
+                update_JxW_values | update_normal_vectors)
+    : fe_values(mapping, fe, quad, update_flags)
+    , fe_interface_values(mapping, fe, quad_face, interface_update_flags)
+  {}
+
+
+  ScratchData(const ScratchData<dim> &scratch_data)
+    : fe_values(scratch_data.fe_values.get_mapping(),
+                scratch_data.fe_values.get_fe(),
+                scratch_data.fe_values.get_quadrature(),
+                scratch_data.fe_values.get_update_flags())
+    , fe_interface_values(scratch_data.fe_values.get_mapping(),
+                          scratch_data.fe_values.get_fe(),
+                          scratch_data.fe_interface_values.get_quadrature(),
+                          scratch_data.fe_interface_values.get_update_flags())
+  {}
+
+  FEValues<dim>          fe_values;
+  FEInterfaceValues<dim> fe_interface_values;
+};
+
+
+
+struct CopyDataFace
+{
+  FullMatrix<double>                   cell_matrix;
+  std::vector<types::global_dof_index> joint_dof_indices;
+};
+
+
+
+struct CopyData
+{
+  FullMatrix<double>                   cell_matrix;
+  Vector<double>                       cell_rhs;
+  std::vector<types::global_dof_index> local_dof_indices;
+  std::vector<CopyDataFace>            face_data;
+
+  template <class Iterator>
+  void
+  reinit(const Iterator &cell, unsigned int dofs_per_cell)
+  {
+    cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+    cell_rhs.reinit(dofs_per_cell);
+
+    local_dof_indices.resize(dofs_per_cell);
+    cell->get_dof_indices(local_dof_indices);
+  }
+};
+
+template <int dim>
+class RightHandSideFunction : public Function<dim>
+{
+public:
+  RightHandSideFunction()
+  {}
+
+  virtual double
+  value(const Point<dim> &p, const unsigned int /*component*/ = 0) const
+  {
+    if (dim == 2)
+      return -2. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1));
+    else /* if(dim == 3)*/
+      return -3. * M_PI * M_PI * std::sin(M_PI * p(0)) * std::sin(M_PI * p(1)) *
+             std::sin(M_PI * p(2));
+  }
+};
+
+template <int dim>
+class DGHeat
+{
+public:
+  DGHeat(const bool           hex,
+         FiniteElement<dim> * fe,
+         Mapping<dim> *       mapping,
+         Quadrature<dim> *    quad,
+         Quadrature<dim - 1> *face_quad,
+         unsigned int         initial_refinement,
+         unsigned int         number_refinement)
+    : hex(hex)
+    , fe(fe)
+    , mapping(mapping)
+    , quad(quad)
+    , face_quad(face_quad)
+    , dof_handler(triangulation)
+    , initial_refinement_level(initial_refinement)
+    , number_refinement(number_refinement)
+  {}
+
+  static std::unique_ptr<DGHeat<dim>>
+  HEX(unsigned int degree,
+      unsigned int initial_refinement,
+      unsigned int number_refinement)
+  {
+    return std::make_unique<DGHeat<dim>>(true,
+                                         new FE_DGQ<dim>(degree),
+                                         new MappingQ<dim>(1),
+                                         new QGauss<dim>(degree + 1),
+                                         new QGauss<dim - 1>(degree + 1),
+                                         initial_refinement,
+                                         number_refinement);
+  }
+
+  static std::unique_ptr<DGHeat<dim>>
+  TET(unsigned int degree,
+      unsigned int initial_refinement,
+      unsigned int number_refinement)
+  {
+    return std::make_unique<DGHeat<dim>>(
+      false,
+      new Simplex::FE_DGP<dim>(degree),
+      new MappingFE<dim>(Simplex::FE_P<dim>(1)),
+      new Simplex::PGauss<dim>(dim == 2 ? (degree == 1 ? 3 : 7) :
+                                          (degree == 1 ? 4 : 10)),
+      new Simplex::PGauss<dim - 1>(dim == 2 ? (degree == 1 ? 2 : 3) :
+                                              (degree == 1 ? 3 : 7)),
+      initial_refinement,
+      number_refinement);
+  }
+
+
+
+  void
+  run();
+
+
+
+private:
+  void
+  make_grid(int refinements = -1);
+  void
+  setup_system();
+  void
+  assemble_system();
+  void
+  solve();
+  void
+  output_results(unsigned int it) const;
+  void
+  calculateL2Error();
+
+  Triangulation<dim> triangulation;
+
+  bool hex;
+
+  std::unique_ptr<FiniteElement<dim>>  fe;
+  const std::unique_ptr<Mapping<dim>>  mapping;
+  std::unique_ptr<Quadrature<dim>>     quad;
+  std::unique_ptr<Quadrature<dim - 1>> face_quad;
+
+  DoFHandler<dim> dof_handler;
+
+
+  RightHandSideFunction<dim> right_hand_side;
+
+
+  SparsityPattern      sparsity_pattern;
+  SparseMatrix<double> system_matrix;
+
+  Vector<double> solution;
+  Vector<double> system_rhs;
+  Point<dim>     center;
+
+  ConvergenceTable error_table;
+
+  unsigned int initial_refinement_level;
+  unsigned int number_refinement;
+};
+
+
+template <int dim>
+void
+DGHeat<dim>::make_grid(int refinements)
+{
+  triangulation.clear();
+
+  const unsigned int ref =
+    refinements == -1 ? initial_refinement_level : refinements;
+
+  if (hex)
+    GridGenerator::subdivided_hyper_cube(triangulation,
+                                         Utilities::pow(2, ref),
+                                         -1.0,
+                                         +1.0);
+  else
+    GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                        Utilities::pow(2, ref),
+                                                        -1.0,
+                                                        +1.0);
+
+  // deallog << "   Number of active cells: " <<
+  // triangulation.n_active_cells()
+  //          << std::endl
+  //          << "   Total number of cells: " << triangulation.n_cells()
+  //          << std::endl;
+}
+
+
+template <int dim>
+void
+DGHeat<dim>::setup_system()
+{
+  dof_handler.distribute_dofs(*fe);
+
+  // deallog << "   Number of degrees of freedom: " << dof_handler.n_dofs()
+  //          << std::endl;
+
+  DynamicSparsityPattern dsp(dof_handler.n_dofs());
+  DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+  sparsity_pattern.copy_from(dsp);
+
+  system_matrix.reinit(sparsity_pattern);
+
+  solution.reinit(dof_handler.n_dofs());
+  system_rhs.reinit(dof_handler.n_dofs());
+}
+
+
+template <int dim>
+void
+DGHeat<dim>::assemble_system()
+{
+  using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+
+  auto cell_worker = [&](const Iterator &  cell,
+                         ScratchData<dim> &scratch_data,
+                         CopyData &        copy_data) {
+    const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+    copy_data.reinit(cell, n_dofs);
+    scratch_data.fe_values.reinit(cell);
+
+    const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+    const FEValues<dim> &      fe_v = scratch_data.fe_values;
+    const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+    std::vector<double> f(q_points.size());
+    right_hand_side.value_list(q_points, f);
+
+    for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+      {
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data.cell_matrix(i, j) +=
+                  fe_v.shape_grad(i, point)   // \nabla \phi_i
+                  * fe_v.shape_grad(j, point) // \nabla \phi_j
+                  * JxW[point];               // dx
+              }
+
+            // Right Hand Side
+            copy_data.cell_rhs(i) +=
+              (fe_v.shape_value(i, point) * f[point] * JxW[point]);
+          }
+      }
+  };
+
+  auto boundary_worker = [&](const Iterator &    cell,
+                             const unsigned int &face_no,
+                             ScratchData<dim> &  scratch_data,
+                             CopyData &          copy_data) {
+    scratch_data.fe_interface_values.reinit(cell, face_no);
+
+    const FEFaceValuesBase<dim> &fe_face =
+      scratch_data.fe_interface_values.get_fe_face_values(0);
+
+    const auto &       q_points     = fe_face.get_quadrature_points();
+    const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+    const std::vector<double> &JxW  = fe_face.get_JxW_values();
+
+    const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+    double h;
+    if (dim == 2)
+      {
+        if (hex)
+          h = std::sqrt(4. * cell->measure() / M_PI);
+        else
+          h = std::sqrt(4. * (4.0 / triangulation.n_cells()) / M_PI);
+      }
+    else if (dim == 3)
+      {
+        if (hex)
+          h = pow(6 * cell->measure() / M_PI, 1. / 3.);
+        else
+          h = pow(6 * (8.0 / triangulation.n_cells()) / M_PI, 1. / 3.);
+      }
+
+
+
+    const double beta = 10.;
+
+    for (unsigned int point = 0; point < q_points.size(); ++point)
+      for (unsigned int i = 0; i < n_facet_dofs; ++i)
+        for (unsigned int j = 0; j < n_facet_dofs; ++j)
+          {
+            copy_data.cell_matrix(i, j) +=
+              -normals[point] * fe_face.shape_grad(i, point) // n*\nabla \phi_i
+              * fe_face.shape_value(j, point)                // \phi_j
+              * JxW[point];                                  // dx
+
+            copy_data.cell_matrix(i, j) +=
+              -fe_face.shape_value(i, point)                  // \phi_i
+              * fe_face.shape_grad(j, point) * normals[point] // n*\nabla \phi_j
+              * JxW[point];                                   // dx
+
+            copy_data.cell_matrix(i, j) +=
+              beta * 1. / h * fe_face.shape_value(i, point) // \phi_i
+              * fe_face.shape_value(j, point) * JxW[point]; // dx
+          }
+  };
+
+  auto face_worker = [&](const Iterator &    cell,
+                         const unsigned int &f,
+                         const unsigned int &sf,
+                         const Iterator &    ncell,
+                         const unsigned int &nf,
+                         const unsigned int &nsf,
+                         ScratchData<dim> &  scratch_data,
+                         CopyData &          copy_data) {
+    FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+
+    fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+
+    const auto &q_points = fe_iv.get_quadrature_points();
+
+    copy_data.face_data.emplace_back();
+    CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+    const unsigned int n_dofs        = fe_iv.n_current_interface_dofs();
+    copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+
+    copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+    const std::vector<double> &        JxW     = fe_iv.get_JxW_values();
+    const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+
+    double h;
+    if (dim == 2)
+      {
+        if (hex)
+          h = std::sqrt(4. * cell->measure() / M_PI);
+        else
+          h = std::sqrt(4. * (4.0 / triangulation.n_cells()) / M_PI);
+      }
+    else if (dim == 3)
+      {
+        if (hex)
+          h = pow(6 * cell->measure() / M_PI, 1. / 3.);
+        else
+          h = pow(6 * (8.0 / triangulation.n_cells()) / M_PI, 1. / 3.);
+      }
+
+    const double beta = 10.;
+
+    for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+      {
+        for (unsigned int i = 0; i < n_dofs; ++i)
+          {
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data_face.cell_matrix(i, j) +=
+                  -normals[qpoint] * fe_iv.average_gradient(i, qpoint) *
+                  fe_iv.jump(j, qpoint) * JxW[qpoint];
+
+                copy_data_face.cell_matrix(i, j) +=
+                  -fe_iv.jump(i, qpoint) // \phi_i
+                  * fe_iv.average_gradient(j, qpoint) *
+                  normals[qpoint] // n*\nabla \phi_j
+                  * JxW[qpoint];  // dx
+
+                copy_data_face.cell_matrix(i, j) +=
+                  beta * 1. / h * fe_iv.jump(i, qpoint) *
+                  fe_iv.jump(j, qpoint) * JxW[qpoint];
+              }
+          }
+      }
+  };
+
+  AffineConstraints<double> constraints;
+
+  auto copier = [&](const CopyData &c) {
+    constraints.distribute_local_to_global(c.cell_matrix,
+                                           c.cell_rhs,
+                                           c.local_dof_indices,
+                                           system_matrix,
+                                           system_rhs);
+
+    for (auto &cdf : c.face_data)
+      {
+        constraints.distribute_local_to_global(cdf.cell_matrix,
+                                               cdf.joint_dof_indices,
+                                               system_matrix);
+      }
+  };
+
+
+  ScratchData<dim> scratch_data(*mapping, *fe, *quad, *face_quad);
+  CopyData         copy_data;
+
+  MeshWorker::mesh_loop(dof_handler.begin_active(),
+                        dof_handler.end(),
+                        cell_worker,
+                        copier,
+                        scratch_data,
+                        copy_data,
+                        MeshWorker::assemble_own_cells |
+                          MeshWorker::assemble_boundary_faces |
+                          MeshWorker::assemble_own_interior_faces_once,
+                        boundary_worker,
+                        face_worker);
+}
+
+template <int dim>
+void
+DGHeat<dim>::solve()
+{
+  SolverControl solver_control(10000, 1e-8);
+  SolverCG<>    solver(solver_control);
+  solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity());
+
+  // We have made one addition, though: since we suppress output from the
+  // linear solvers, we have to print the number of iterations by hand.
+  // deallog << "   " << solver_control.last_step()
+  //          << " CG iterations needed to obtain convergence." << std::endl;
+
+  // error_table.add_value("iterations", solver_control.last_step());
+}
+
+template <int dim>
+void
+DGHeat<dim>::output_results(unsigned int it) const
+{
+  return;
+
+  std::string type = hex ? "hex" : "tet";
+
+  std::string dimension(dim == 2 ? "solution-2d-" + type + "-case-" :
+                                   "solution-3d-" + type + "-case-");
+
+  std::string fname = dimension + Utilities::int_to_string(it) + ".vtk";
+
+  deallog << "  Writing solution to <" << fname << ">" << std::endl;
+
+  std::ofstream output(fname.c_str());
+
+  if (false)
+    {
+      DataOut<dim> data_out;
+
+      data_out.attach_dof_handler(dof_handler);
+      data_out.add_data_vector(solution, "solution");
+
+      data_out.build_patches(*mapping);
+      data_out.write_vtk(output);
+    }
+}
+
+// Find the l2 norm of the error between the finite element sol'n and the exact
+// sol'n
+template <int dim>
+void
+DGHeat<dim>::calculateL2Error()
+{
+  FEValues<dim> fe_values(*mapping,
+                          *fe,
+                          *quad,
+                          update_values | update_gradients |
+                            update_quadrature_points | update_JxW_values);
+
+  const unsigned int dofs_per_cell =
+    fe->dofs_per_cell; // This gives you dofs per cell
+  std::vector<types::global_dof_index> local_dof_indices(
+    dofs_per_cell); //  Local connectivity
+
+  const unsigned int n_q_points = quad->size();
+
+  double l2error = 0.;
+
+  // loop over elements
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    {
+      fe_values.reinit(cell);
+
+      cell->get_dof_indices(local_dof_indices);
+
+      for (unsigned int q = 0; q < n_q_points; q++)
+        {
+          const double u_exact =
+            dim == 2 ? -std::sin(M_PI * fe_values.quadrature_point(q)[0]) *
+                         std::sin(M_PI * fe_values.quadrature_point(q)[1]) :
+                       -std::sin(M_PI * fe_values.quadrature_point(q)[0]) *
+                         std::sin(M_PI * fe_values.quadrature_point(q)[1]) *
+                         std::sin(M_PI * fe_values.quadrature_point(q)[2]);
+
+          double u_sim = 0;
+
+          // Find the values of x and u_h (the finite element solution) at the
+          // quadrature points
+          for (unsigned int i = 0; i < dofs_per_cell; i++)
+            {
+              u_sim +=
+                fe_values.shape_value(i, q) * solution[local_dof_indices[i]];
+            }
+          l2error += (u_sim - u_exact) * (u_sim - u_exact) * fe_values.JxW(q);
+          //       deallog << " x = " << x << " y = " << y <<  " r = " << r <<
+          //       "   u_exact = " << u_exact << "   u_sim=" << u_sim <<
+          //       std::endl;
+        }
+    }
+
+
+  // deallog << "L2Error is : " << std::sqrt(l2error) << std::endl;
+  error_table.add_value("error", std::sqrt(l2error));
+  error_table.add_value("cells", triangulation.n_global_active_cells());
+  error_table.add_value("dofs", dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void
+DGHeat<dim>::run()
+{
+  for (unsigned int it = 0; it < number_refinement; ++it)
+    {
+      make_grid(initial_refinement_level + it);
+      setup_system();
+      assemble_system();
+      solve();
+      output_results(it);
+      calculateL2Error();
+    }
+
+  // error_table.omit_column_from_convergence_rate_evaluation("iterations");
+  error_table.omit_column_from_convergence_rate_evaluation("cells");
+  error_table.evaluate_all_convergence_rates(
+    ConvergenceTable::reduction_rate_log2);
+
+  error_table.set_scientific("error", true);
+
+  error_table.write_text(deallog.get_file_stream());
+  deallog << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  deallog.depth_file(1);
+
+
+  {
+    auto problem = DGHeat<2>::TET(1 /*=degree*/, 2, 3);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<2>::TET(2 /*=degree*/, 2, 3);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<3>::TET(1 /*=degree*/, 2, 2);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<3>::TET(2 /*=degree*/, 2, 2);
+    problem->run();
+  }
+
+  {
+    auto problem = DGHeat<2>::HEX(1 /*=degree*/, 2, 3);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<2>::HEX(2 /*=degree*/, 2, 3);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<3>::HEX(1 /*=degree*/, 2, 2);
+    problem->run();
+  }
+  {
+    auto problem = DGHeat<3>::HEX(2 /*=degree*/, 2, 2);
+    problem->run();
+  }
+
+  return 0;
+}
diff --git a/tests/simplex/poisson_02.with_simplex_support=on.output b/tests/simplex/poisson_02.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..825396d
--- /dev/null
@@ -0,0 +1,37 @@
+
+     error      cells    dofs    
+3.0924e-01 -    32    96   -     
+1.1363e-01 1.44 128   384  -2.00 
+3.2531e-02 1.80 512   1536 -2.00 
+DEAL::
+     error      cells    dofs    
+3.2745e-02 -    32    192  -     
+3.6158e-03 3.18 128   768  -2.00 
+4.1566e-04 3.12 512   3072 -2.00 
+DEAL::
+     error      cells    dofs     
+1.6346e-01 -    320   1280  -     
+1.0816e-01 0.60 2560  10240 -3.00 
+DEAL::
+     error      cells    dofs     
+8.4685e-02 -    320   3200  -     
+9.4765e-03 3.16 2560  25600 -3.00 
+DEAL::
+     error      cells    dofs    
+1.6226e-01 -    16    64   -     
+4.8634e-02 1.74 64    256  -2.00 
+1.2684e-02 1.94 256   1024 -2.00 
+DEAL::
+     error      cells    dofs    
+1.6600e-02 -    16    144  -     
+1.6552e-03 3.33 64    576  -2.00 
+1.7799e-04 3.22 256   2304 -2.00 
+DEAL::
+     error      cells    dofs    
+1.6067e-01 -    64    512  -     
+4.8569e-02 1.73 512   4096 -3.00 
+DEAL::
+     error      cells    dofs     
+2.2671e-02 -    64    1728  -     
+2.2931e-03 3.31 512   13824 -3.00 
+DEAL::
diff --git a/tests/simplex/step-12.cc b/tests/simplex/step-12.cc
new file mode 100644 (file)
index 0000000..e1f1233
--- /dev/null
@@ -0,0 +1,563 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-12 with tetrahedron mesh.
+
+#include <deal.II/base/function.h>
+#include <deal.II/base/quadrature_lib.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_interface_values.h>
+#include <deal.II/fe/fe_values.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/lac/dynamic_sparsity_pattern.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 <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+// Finally, the new include file for using the mesh_loop from the MeshWorker
+// framework
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/meshworker/mesh_loop.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+#include <fstream>
+#include <iostream>
+
+//#define HEX
+
+namespace Step12
+{
+  using namespace dealii;
+
+  template <int dim>
+  class BoundaryValues : public Function<dim>
+  {
+  public:
+    BoundaryValues() = default;
+    virtual void
+    value_list(const std::vector<Point<dim>> &points,
+               std::vector<double> &          values,
+               const unsigned int             component = 0) const override;
+  };
+
+  template <int dim>
+  void
+  BoundaryValues<dim>::value_list(const std::vector<Point<dim>> &points,
+                                  std::vector<double> &          values,
+                                  const unsigned int component) const
+  {
+    (void)component;
+    AssertIndexRange(component, 1);
+    Assert(values.size() == points.size(),
+           ExcDimensionMismatch(values.size(), points.size()));
+
+    for (unsigned int i = 0; i < values.size(); ++i)
+      {
+        if (points[i](0) < 0.5)
+          values[i] = 1.;
+        else
+          values[i] = 0.;
+      }
+  }
+
+
+  template <int dim>
+  Tensor<1, dim>
+  beta(const Point<dim> &p)
+  {
+    Assert(dim >= 2, ExcNotImplemented());
+
+    Point<dim> wind_field;
+    wind_field(0) = -p(1);
+    wind_field(1) = p(0);
+
+    if (wind_field.norm() > 1e-6)
+      wind_field /= wind_field.norm();
+
+    return wind_field;
+  }
+
+
+  template <int dim>
+  struct ScratchData
+  {
+    ScratchData(const Mapping<dim> &       mapping,
+                const FiniteElement<dim> & fe,
+                const Quadrature<dim> &    quad,
+                const Quadrature<dim - 1> &quad_face,
+                const UpdateFlags          update_flags = update_values |
+                                                 update_gradients |
+                                                 update_quadrature_points |
+                                                 update_JxW_values,
+                const UpdateFlags interface_update_flags =
+                  update_values | update_gradients | update_quadrature_points |
+                  update_JxW_values | update_normal_vectors)
+      : fe_values(mapping, fe, quad, update_flags)
+      , fe_interface_values(mapping, fe, quad_face, interface_update_flags)
+    {}
+
+
+    ScratchData(const ScratchData<dim> &scratch_data)
+      : fe_values(scratch_data.fe_values.get_mapping(),
+                  scratch_data.fe_values.get_fe(),
+                  scratch_data.fe_values.get_quadrature(),
+                  scratch_data.fe_values.get_update_flags())
+      , fe_interface_values(
+          scratch_data.fe_values
+            .get_mapping(), // TODO: implement for fe_interface_values
+          scratch_data.fe_values.get_fe(),
+          scratch_data.fe_interface_values.get_quadrature(),
+          scratch_data.fe_interface_values.get_update_flags())
+    {}
+
+    FEValues<dim>          fe_values;
+    FEInterfaceValues<dim> fe_interface_values;
+  };
+
+
+
+  struct CopyDataFace
+  {
+    FullMatrix<double>                   cell_matrix;
+    std::vector<types::global_dof_index> joint_dof_indices;
+  };
+
+
+
+  struct CopyData
+  {
+    FullMatrix<double>                   cell_matrix;
+    Vector<double>                       cell_rhs;
+    std::vector<types::global_dof_index> local_dof_indices;
+    std::vector<CopyDataFace>            face_data;
+
+    template <class Iterator>
+    void
+    reinit(const Iterator &cell, unsigned int dofs_per_cell)
+    {
+      cell_matrix.reinit(dofs_per_cell, dofs_per_cell);
+      cell_rhs.reinit(dofs_per_cell);
+
+      local_dof_indices.resize(dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+    }
+  };
+
+
+  template <int dim>
+  class AdvectionProblem
+  {
+  public:
+    AdvectionProblem();
+    void
+    run();
+
+  private:
+    void
+    setup_system();
+    void
+    assemble_system();
+    void
+    solve();
+    void
+    refine_grid();
+    void
+    output_results(const unsigned int cycle) const;
+
+    Triangulation<dim> triangulation;
+#ifdef HEX
+    const MappingQ1<dim> mapping;
+#else
+    Simplex::FE_P<dim>   fe_mapping;
+    const MappingFE<dim> mapping;
+#endif
+
+    // Furthermore we want to use DG elements.
+#ifdef HEX
+    FE_DGQ<dim> fe;
+#else
+    Simplex::FE_DGP<dim> fe;
+#endif
+    DoFHandler<dim> dof_handler;
+
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+
+    Vector<double> solution;
+    Vector<double> right_hand_side;
+  };
+
+
+  template <int dim>
+  AdvectionProblem<dim>::AdvectionProblem()
+#ifdef HEX
+    : mapping()
+    , fe(2)
+#else
+    : fe_mapping(1)
+    , mapping(fe_mapping)
+    , fe(2)
+#endif
+    , dof_handler(triangulation)
+  {}
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+
+    DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    DoFTools::make_flux_sparsity_pattern(dof_handler, dsp);
+    sparsity_pattern.copy_from(dsp);
+
+    system_matrix.reinit(sparsity_pattern);
+    solution.reinit(dof_handler.n_dofs());
+    right_hand_side.reinit(dof_handler.n_dofs());
+  }
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::assemble_system()
+  {
+    using Iterator = typename DoFHandler<dim>::active_cell_iterator;
+    const BoundaryValues<dim> boundary_function;
+
+    auto cell_worker = [&](const Iterator &  cell,
+                           ScratchData<dim> &scratch_data,
+                           CopyData &        copy_data) {
+      const unsigned int n_dofs = scratch_data.fe_values.get_fe().dofs_per_cell;
+      copy_data.reinit(cell, n_dofs);
+      scratch_data.fe_values.reinit(cell);
+
+      const auto &q_points = scratch_data.fe_values.get_quadrature_points();
+
+      const FEValues<dim> &      fe_v = scratch_data.fe_values;
+      const std::vector<double> &JxW  = fe_v.get_JxW_values();
+
+      for (unsigned int point = 0; point < fe_v.n_quadrature_points; ++point)
+        {
+          auto beta_q = beta(q_points[point]);
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              {
+                copy_data.cell_matrix(i, j) +=
+                  -beta_q                      // -\beta
+                  * fe_v.shape_grad(i, point)  // \nabla \phi_i
+                  * fe_v.shape_value(j, point) // \phi_j
+                  * JxW[point];                // dx
+              }
+        }
+    };
+
+    auto boundary_worker = [&](const Iterator &    cell,
+                               const unsigned int &face_no,
+                               ScratchData<dim> &  scratch_data,
+                               CopyData &          copy_data) {
+      scratch_data.fe_interface_values.reinit(cell, face_no);
+      const FEFaceValuesBase<dim> &fe_face =
+        scratch_data.fe_interface_values.get_fe_face_values(0);
+
+      const auto &q_points = fe_face.get_quadrature_points();
+
+      const unsigned int n_facet_dofs = fe_face.get_fe().n_dofs_per_cell();
+      const std::vector<double> &        JxW     = fe_face.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_face.get_normal_vectors();
+
+      std::vector<double> g(q_points.size());
+      boundary_function.value_list(q_points, g);
+
+      for (unsigned int point = 0; point < q_points.size(); ++point)
+        {
+          const double beta_dot_n = beta(q_points[point]) * normals[point];
+
+          if (beta_dot_n > 0)
+            {
+              for (unsigned int i = 0; i < n_facet_dofs; ++i)
+                for (unsigned int j = 0; j < n_facet_dofs; ++j)
+                  copy_data.cell_matrix(i, j) +=
+                    fe_face.shape_value(i, point)   // \phi_i
+                    * fe_face.shape_value(j, point) // \phi_j
+                    * beta_dot_n                    // \beta . n
+                    * JxW[point];                   // dx
+            }
+          else
+            for (unsigned int i = 0; i < n_facet_dofs; ++i)
+              copy_data.cell_rhs(i) += -fe_face.shape_value(i, point) // \phi_i
+                                       * g[point]                     // g
+                                       * beta_dot_n  // \beta . n
+                                       * JxW[point]; // dx
+        }
+    };
+
+    auto face_worker = [&](const Iterator &    cell,
+                           const unsigned int &f,
+                           const unsigned int &sf,
+                           const Iterator &    ncell,
+                           const unsigned int &nf,
+                           const unsigned int &nsf,
+                           ScratchData<dim> &  scratch_data,
+                           CopyData &          copy_data) {
+      FEInterfaceValues<dim> &fe_iv = scratch_data.fe_interface_values;
+      fe_iv.reinit(cell, f, sf, ncell, nf, nsf);
+      const auto &q_points = fe_iv.get_quadrature_points();
+
+      copy_data.face_data.emplace_back();
+      CopyDataFace &copy_data_face = copy_data.face_data.back();
+
+      const unsigned int n_dofs        = fe_iv.n_current_interface_dofs();
+      copy_data_face.joint_dof_indices = fe_iv.get_interface_dof_indices();
+
+      copy_data_face.cell_matrix.reinit(n_dofs, n_dofs);
+
+      const std::vector<double> &        JxW     = fe_iv.get_JxW_values();
+      const std::vector<Tensor<1, dim>> &normals = fe_iv.get_normal_vectors();
+
+      for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint)
+        {
+          const double beta_dot_n = beta(q_points[qpoint]) * normals[qpoint];
+          for (unsigned int i = 0; i < n_dofs; ++i)
+            for (unsigned int j = 0; j < n_dofs; ++j)
+              copy_data_face.cell_matrix(i, j) +=
+                fe_iv.jump(i, qpoint) // [\phi_i]
+                *
+                fe_iv.shape_value((beta_dot_n > 0), j, qpoint) // phi_j^{upwind}
+                * beta_dot_n                                   // (\beta . n)
+                * JxW[qpoint];                                 // dx
+        }
+    };
+
+    AffineConstraints<double> constraints;
+
+    auto copier = [&](const CopyData &c) {
+      constraints.distribute_local_to_global(c.cell_matrix,
+                                             c.cell_rhs,
+                                             c.local_dof_indices,
+                                             system_matrix,
+                                             right_hand_side);
+
+      for (auto &cdf : c.face_data)
+        {
+          constraints.distribute_local_to_global(cdf.cell_matrix,
+                                                 cdf.joint_dof_indices,
+                                                 system_matrix);
+        }
+    };
+
+    const unsigned int degree = dof_handler.get_fe().degree;
+
+#ifdef HEX
+    QGauss<dim> quad(degree + 1);
+
+    QGauss<dim - 1> face_quad(degree + 1);
+#else
+    Simplex::PGauss<dim> quad(dim == 2 ? (degree == 1 ? 3 : 7) :
+                                         (degree == 1 ? 4 : 10));
+
+    Simplex::PGauss<dim - 1> face_quad(dim == 2 ? (degree == 1 ? 2 : 3) :
+                                                  (degree == 1 ? 3 : 7));
+#endif
+
+    ScratchData<dim> scratch_data(mapping, fe, quad, face_quad);
+    CopyData         copy_data;
+
+    MeshWorker::mesh_loop(dof_handler.begin_active(),
+                          dof_handler.end(),
+                          cell_worker,
+                          copier,
+                          scratch_data,
+                          copy_data,
+                          MeshWorker::assemble_own_cells |
+                            MeshWorker::assemble_boundary_faces |
+                            MeshWorker::assemble_own_interior_faces_once,
+                          boundary_worker,
+                          face_worker);
+  }
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::solve()
+  {
+    SolverControl                    solver_control(1000, 1e-12);
+    SolverRichardson<Vector<double>> solver(solver_control);
+
+    PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
+
+    preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+
+    solver.solve(system_matrix, solution, right_hand_side, preconditioner);
+
+    deallog << "  Solver converged in " << solver_control.last_step()
+            << " iterations." << std::endl;
+  }
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::refine_grid()
+  {
+    Vector<float> gradient_indicator(triangulation.n_active_cells());
+
+    DerivativeApproximation::approximate_gradient(mapping,
+                                                  dof_handler,
+                                                  solution,
+                                                  gradient_indicator);
+
+    unsigned int cell_no = 0;
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      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
+  AdvectionProblem<dim>::output_results(const unsigned int cycle) const
+  {
+#if false
+#  ifdef HEX
+    const std::string filename =
+      dim == 2 ? ("step12-quad-" + std::to_string(cycle) + ".vtk") :
+                 ("step12-hex-" + std::to_string(cycle) + ".vtk");
+#  else
+    const std::string filename =
+      dim == 2 ? ("step12-tri-" + std::to_string(cycle) + ".vtk") :
+                 ("step12-tet-" + std::to_string(cycle) + ".vtk");
+#  endif
+    deallog << "  Writing solution to <" << filename << ">" << std::endl;
+    std::ofstream output(filename);
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+    data_out.add_data_vector(solution, "u", DataOut<dim>::type_dof_data);
+
+    data_out.build_patches(mapping, 2);
+
+    data_out.write_vtk(output);
+#endif
+
+    {
+      Vector<float> values(triangulation.n_active_cells());
+      VectorTools::integrate_difference(mapping,
+                                        dof_handler,
+                                        solution,
+                                        Functions::ZeroFunction<dim>(),
+                                        values,
+                                        QGauss<dim>(fe.degree + 1),
+                                        VectorTools::Linfty_norm);
+      const double l_infty =
+        VectorTools::compute_global_error(triangulation,
+                                          values,
+                                          VectorTools::Linfty_norm);
+      deallog << "  L-infinity norm: " << l_infty << std::endl;
+    }
+  }
+
+
+  template <int dim>
+  void
+  AdvectionProblem<dim>::run()
+  {
+    //#ifdef HEX
+    //    for (unsigned int cycle = 0; cycle < 6; ++cycle)
+    //#else
+    for (unsigned int cycle = 0; cycle < 1; ++cycle)
+      //#endif
+      {
+        deallog << "Cycle " << cycle << std::endl;
+
+        if (cycle == 0)
+          {
+#ifdef HEX
+            // GridGenerator::hyper_cube(triangulation);
+            // triangulation.refine_global(3);
+            GridGenerator::subdivided_hyper_cube(triangulation, 16);
+#else
+            GridGenerator::subdivided_hyper_cube_with_simplices(triangulation,
+                                                                dim == 2 ? 32 :
+                                                                           8);
+#endif
+          }
+        else
+          refine_grid();
+
+        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;
+
+        assemble_system();
+        solve();
+
+        output_results(cycle);
+      }
+  }
+} // namespace Step12
+
+
+int
+main()
+{
+  initlog();
+
+  {
+    Step12::AdvectionProblem<2> dgmethod;
+    dgmethod.run();
+  }
+  {
+    Step12::AdvectionProblem<3> dgmethod;
+    dgmethod.run();
+  }
+
+  return 0;
+}
diff --git a/tests/simplex/step-12.with_simplex_support=on.output b/tests/simplex/step-12.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..7fbe5eb
--- /dev/null
@@ -0,0 +1,15 @@
+
+DEAL::Cycle 0
+DEAL::  Number of active cells:       2048
+DEAL::  Number of degrees of freedom: 12288
+DEAL:Richardson::Starting value 0.0883883
+DEAL:Richardson::Convergence step 16 value 9.71637e-17
+DEAL::  Solver converged in 16 iterations.
+DEAL::  L-infinity norm: 2.33280
+DEAL::Cycle 0
+DEAL::  Number of active cells:       2560
+DEAL::  Number of degrees of freedom: 25600
+DEAL:Richardson::Starting value 0.0360844
+DEAL:Richardson::Convergence step 13 value 2.52800e-14
+DEAL::  Solver converged in 13 iterations.
+DEAL::  L-infinity norm: 13.9260
diff --git a/tests/simplex/step-18.cc b/tests/simplex/step-18.cc
new file mode 100644 (file)
index 0000000..a146186
--- /dev/null
@@ -0,0 +1,1053 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Step-18 with tetrahedron mesh.
+
+#include <deal.II/base/conditional_ostream.h>
+#include <deal.II/base/function.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/petsc_precondition.h>
+#include <deal.II/lac/petsc_solver.h>
+#include <deal.II/lac/petsc_sparse_matrix.h>
+#include <deal.II/lac/petsc_vector.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparsity_tools.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/error_estimator.h>
+#include <deal.II/numerics/matrix_tools.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/physics/transformations.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+// simplex
+#include <deal.II/fe/mapping_fe.h>
+
+#include <deal.II/simplex/fe_lib.h>
+#include <deal.II/simplex/grid_generator.h>
+#include <deal.II/simplex/quadrature_lib.h>
+
+//#define HEX
+
+const unsigned int degree = 1;
+
+namespace Step18
+{
+  using namespace dealii;
+
+  template <int dim>
+  struct PointHistory
+  {
+    SymmetricTensor<2, dim> old_stress;
+  };
+
+  template <int dim>
+  SymmetricTensor<4, dim>
+  get_stress_strain_tensor(const double lambda, const double mu)
+  {
+    SymmetricTensor<4, dim> tmp;
+    for (unsigned int i = 0; i < dim; ++i)
+      for (unsigned int j = 0; j < dim; ++j)
+        for (unsigned int k = 0; k < dim; ++k)
+          for (unsigned int l = 0; l < dim; ++l)
+            tmp[i][j][k][l] = (((i == k) && (j == l) ? mu : 0.0) +
+                               ((i == l) && (j == k) ? mu : 0.0) +
+                               ((i == j) && (k == l) ? lambda : 0.0));
+    return tmp;
+  }
+
+  template <int dim>
+  inline SymmetricTensor<2, dim>
+  get_strain(const FEValues<dim> &fe_values,
+             const unsigned int   shape_func,
+             const unsigned int   q_point)
+  {
+    SymmetricTensor<2, dim> tmp;
+
+    for (unsigned int i = 0; i < dim; ++i)
+      tmp[i][i] = fe_values.shape_grad_component(shape_func, q_point, i)[i];
+
+    for (unsigned int i = 0; i < dim; ++i)
+      for (unsigned int j = i + 1; j < dim; ++j)
+        tmp[i][j] =
+          (fe_values.shape_grad_component(shape_func, q_point, i)[j] +
+           fe_values.shape_grad_component(shape_func, q_point, j)[i]) /
+          2;
+
+    return tmp;
+  }
+
+
+  template <int dim>
+  inline SymmetricTensor<2, dim>
+  get_strain(const std::vector<Tensor<1, dim>> &grad)
+  {
+    Assert(grad.size() == dim, ExcInternalError());
+
+    SymmetricTensor<2, dim> strain;
+    for (unsigned int i = 0; i < dim; ++i)
+      strain[i][i] = grad[i][i];
+
+    for (unsigned int i = 0; i < dim; ++i)
+      for (unsigned int j = i + 1; j < dim; ++j)
+        strain[i][j] = (grad[i][j] + grad[j][i]) / 2;
+
+    return strain;
+  }
+
+
+  Tensor<2, 2>
+  get_rotation_matrix(const std::vector<Tensor<1, 2>> &grad_u)
+  {
+    const double curl = (grad_u[1][0] - grad_u[0][1]);
+
+    const double angle = std::atan(curl);
+
+    return Physics::Transformations::Rotations::rotation_matrix_2d(-angle);
+  }
+
+
+  Tensor<2, 3>
+  get_rotation_matrix(const std::vector<Tensor<1, 3>> &grad_u)
+  {
+    const Point<3> curl(grad_u[2][1] - grad_u[1][2],
+                        grad_u[0][2] - grad_u[2][0],
+                        grad_u[1][0] - grad_u[0][1]);
+
+    const double tan_angle = std::sqrt(curl * curl);
+    const double angle     = std::atan(tan_angle);
+
+    if (std::abs(angle) < 1e-9)
+      {
+        static const double rotation[3][3] = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
+        static const Tensor<2, 3> rot(rotation);
+        return rot;
+      }
+
+    const Point<3> axis = curl / tan_angle;
+    return Physics::Transformations::Rotations::rotation_matrix_3d(axis,
+                                                                   -angle);
+  }
+
+
+
+  template <int dim>
+  class TopLevel
+  {
+  public:
+    TopLevel();
+    ~TopLevel();
+    void
+    run();
+
+  private:
+    void
+    create_coarse_grid();
+
+    void
+    setup_system();
+
+    void
+    assemble_system();
+
+    void
+    solve_timestep();
+
+    unsigned int
+    solve_linear_problem();
+
+    void
+    output_results() const;
+
+    void
+    do_initial_timestep(const bool do_output = false);
+
+    void
+    do_timestep(const bool do_output = false);
+
+    void
+    refine_initial_grid();
+
+    void
+    move_mesh();
+
+    void
+    setup_quadrature_point_history();
+
+    void
+    update_quadrature_point_history();
+
+    Triangulation<dim> triangulation;
+
+    FESystem<dim> fe;
+
+    DoFHandler<dim> dof_handler;
+
+    AffineConstraints<double> hanging_node_constraints;
+
+    const Quadrature<dim> quadrature_formula;
+
+#ifdef HEX
+    MappingQGeneric<dim, dim> mapping;
+#else
+    MappingFE<dim, dim>  mapping;
+#endif
+
+    std::vector<PointHistory<dim>> quadrature_point_history;
+
+#ifdef DEAL_II_WITH_PETSC
+    PETScWrappers::MPI::SparseMatrix system_matrix;
+    PETScWrappers::MPI::Vector       system_rhs;
+#else
+    SparsityPattern      sparsity_pattern;
+    SparseMatrix<double> system_matrix;
+    Vector<double>       system_rhs;
+#endif
+
+    Vector<double> incremental_displacement;
+
+    double       present_time;
+    double       present_timestep;
+    double       end_time;
+    unsigned int timestep_no;
+
+    MPI_Comm mpi_communicator;
+
+    const unsigned int n_mpi_processes;
+
+    const unsigned int this_mpi_process;
+
+    ConditionalOStream pcout;
+
+    IndexSet locally_owned_dofs;
+    IndexSet locally_relevant_dofs;
+
+    static const SymmetricTensor<4, dim> stress_strain_tensor;
+  };
+
+
+  template <int dim>
+  class BodyForce : public Function<dim>
+  {
+  public:
+    BodyForce();
+
+    virtual void
+    vector_value(const Point<dim> &p, Vector<double> &values) const override;
+
+    virtual void
+    vector_value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Vector<double>> &  value_list) const override;
+  };
+
+
+  template <int dim>
+  BodyForce<dim>::BodyForce()
+    : Function<dim>(dim)
+  {}
+
+
+  template <int dim>
+  inline void
+  BodyForce<dim>::vector_value(const Point<dim> & /*p*/,
+                               Vector<double> &values) const
+  {
+    Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
+
+    const double g   = 9.81;
+    const double rho = 7700;
+
+    values          = 0;
+    values(dim - 1) = -rho * g;
+  }
+
+
+
+  template <int dim>
+  void
+  BodyForce<dim>::vector_value_list(
+    const std::vector<Point<dim>> &points,
+    std::vector<Vector<double>> &  value_list) const
+  {
+    const unsigned int n_points = points.size();
+
+    Assert(value_list.size() == n_points,
+           ExcDimensionMismatch(value_list.size(), n_points));
+
+    for (unsigned int p = 0; p < n_points; ++p)
+      BodyForce<dim>::vector_value(points[p], value_list[p]);
+  }
+
+
+
+  template <int dim>
+  class IncrementalBoundaryValues : public Function<dim>
+  {
+  public:
+    IncrementalBoundaryValues(const double present_time,
+                              const double present_timestep);
+
+    virtual void
+    vector_value(const Point<dim> &p, Vector<double> &values) const override;
+
+    virtual void
+    vector_value_list(const std::vector<Point<dim>> &points,
+                      std::vector<Vector<double>> &  value_list) const override;
+
+  private:
+    const double velocity;
+    const double present_time;
+    const double present_timestep;
+  };
+
+
+  template <int dim>
+  IncrementalBoundaryValues<dim>::IncrementalBoundaryValues(
+    const double present_time,
+    const double present_timestep)
+    : Function<dim>(dim)
+    , velocity(.08)
+    , present_time(present_time)
+    , present_timestep(present_timestep)
+  {}
+
+
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::vector_value(const Point<dim> & /*p*/,
+                                               Vector<double> &values) const
+  {
+    Assert(values.size() == dim, ExcDimensionMismatch(values.size(), dim));
+
+    values    = 0;
+    values(2) = -present_timestep * velocity;
+  }
+
+
+
+  template <int dim>
+  void
+  IncrementalBoundaryValues<dim>::vector_value_list(
+    const std::vector<Point<dim>> &points,
+    std::vector<Vector<double>> &  value_list) const
+  {
+    const unsigned int n_points = points.size();
+
+    Assert(value_list.size() == n_points,
+           ExcDimensionMismatch(value_list.size(), n_points));
+
+    for (unsigned int p = 0; p < n_points; ++p)
+      IncrementalBoundaryValues<dim>::vector_value(points[p], value_list[p]);
+  }
+
+
+  template <int dim>
+  const SymmetricTensor<4, dim> TopLevel<dim>::stress_strain_tensor =
+    get_stress_strain_tensor<dim>(/*lambda = */ 9.695e10,
+                                  /*mu     = */ 7.617e10);
+
+
+#ifdef HEX
+  template <int dim>
+  TopLevel<dim>::TopLevel()
+    : triangulation()
+    , fe(FE_Q<dim>(degree), dim)
+    , dof_handler(triangulation)
+    , quadrature_formula(QGauss<dim>(fe.degree + 1))
+    , mapping(1)
+    , present_time(0.0)
+    , present_timestep(1.0)
+    , end_time(10.0)
+    , timestep_no(0)
+    , mpi_communicator(MPI_COMM_WORLD)
+    , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
+    , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
+    , pcout(std::cout, this_mpi_process == 0)
+  {}
+#else
+  template <int dim>
+  TopLevel<dim>::TopLevel()
+    : triangulation()
+    , fe(Simplex::FE_P<dim>(degree), dim)
+    , dof_handler(triangulation)
+    , quadrature_formula(Simplex::PGauss<dim>(fe.degree == 1 ? 4 : 10))
+    , mapping(Simplex::FE_P<dim>(1))
+    , present_time(0.0)
+    , present_timestep(1.0)
+    , end_time(10.0)
+    , timestep_no(0)
+    , mpi_communicator(MPI_COMM_WORLD)
+    , n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator))
+    , this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator))
+    , pcout(std::cout, this_mpi_process == 0)
+  {}
+#endif
+
+
+
+  template <int dim>
+  TopLevel<dim>::~TopLevel()
+  {
+    dof_handler.clear();
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::run()
+  {
+    do_initial_timestep(false);
+
+    while (present_time < end_time)
+      do_timestep(std::abs(end_time - present_time - present_timestep) < 10e-5);
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::create_coarse_grid()
+  {
+    const unsigned int n = 5;
+
+#ifdef HEX
+    GridGenerator::subdivided_hyper_rectangle(triangulation,
+                                              {1 * n, 1 * n, 3 * n},
+                                              {-0.5, -0.5, 0},
+                                              {+0.5, +0.5, +3});
+#else
+    GridGenerator::subdivided_hyper_rectangle_with_simplices(
+      triangulation, {1 * n, 1 * n, 3 * n}, {-0.5, -0.5, 0}, {+0.5, +0.5, +3});
+#endif
+
+    for (const auto &cell : triangulation.active_cell_iterators())
+      for (const auto &face : cell->face_iterators())
+        if (face->at_boundary())
+          {
+            const Point<dim> face_center = face->center();
+
+            if (face_center[2] == 0)
+              face->set_boundary_id(0);
+            else if (face_center[2] == 3)
+              face->set_boundary_id(1);
+            else
+              face->set_boundary_id(2);
+          }
+
+    setup_quadrature_point_history();
+  }
+
+  template <int dim>
+  void
+  TopLevel<dim>::setup_system()
+  {
+    dof_handler.distribute_dofs(fe);
+    locally_owned_dofs = dof_handler.locally_owned_dofs();
+    DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+
+    // The next step is to set up constraints due to hanging nodes. This has
+    // been handled many times before:
+    hanging_node_constraints.clear();
+    DoFTools::make_hanging_node_constraints(dof_handler,
+                                            hanging_node_constraints);
+    hanging_node_constraints.close();
+
+#ifdef DEAL_II_WITH_PETSC
+    DynamicSparsityPattern sparsity_pattern(locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    sparsity_pattern,
+                                    hanging_node_constraints,
+                                    /*keep constrained dofs*/ false);
+    SparsityTools::distribute_sparsity_pattern(sparsity_pattern,
+                                               locally_owned_dofs,
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+
+    system_matrix.reinit(locally_owned_dofs,
+                         locally_owned_dofs,
+                         sparsity_pattern,
+                         mpi_communicator);
+
+    system_rhs.reinit(locally_owned_dofs, mpi_communicator);
+    incremental_displacement.reinit(dof_handler.n_dofs());
+
+#else
+    // DynamicSparsityPattern dsp(dof_handler.n_dofs());
+    // DoFTools::make_sparsity_pattern(dof_handler, dsp,
+    // hanging_node_constraints, false);
+    //  sparsity_pattern.copy_from(dsp);
+    //  system_matrix.reinit(sparsity_pattern);
+    //  sparsity_pattern.copy_from(dsp);
+    //  system_matrix.reinit(sparsity_pattern);
+
+    DynamicSparsityPattern dsp(locally_relevant_dofs);
+    DoFTools::make_sparsity_pattern(dof_handler,
+                                    dsp,
+                                    hanging_node_constraints,
+                                    /*keep constrained dofs*/ false);
+    SparsityTools::distribute_sparsity_pattern(dsp,
+                                               locally_owned_dofs,
+                                               mpi_communicator,
+                                               locally_relevant_dofs);
+    sparsity_pattern.copy_from(dsp);
+    system_matrix.reinit(sparsity_pattern);
+
+    system_rhs.reinit(dof_handler.n_dofs());
+    incremental_displacement.reinit(dof_handler.n_dofs());
+#endif
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::assemble_system()
+  {
+    system_rhs    = 0;
+    system_matrix = 0;
+
+    FEValues<dim> fe_values(mapping,
+                            fe,
+                            quadrature_formula,
+                            update_values | update_gradients |
+                              update_quadrature_points | update_JxW_values);
+
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const unsigned int n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
+    Vector<double>     cell_rhs(dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
+
+    BodyForce<dim>              body_force;
+    std::vector<Vector<double>> body_force_values(n_q_points,
+                                                  Vector<double>(dim));
+
+    for (const auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell_matrix = 0;
+          cell_rhs    = 0;
+
+          fe_values.reinit(cell);
+
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            for (unsigned int j = 0; j < dofs_per_cell; ++j)
+              for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+                {
+                  const SymmetricTensor<2, dim>
+                    eps_phi_i = get_strain(fe_values, i, q_point),
+                    eps_phi_j = get_strain(fe_values, j, q_point);
+
+                  cell_matrix(i, j) += (eps_phi_i *            //
+                                        stress_strain_tensor * //
+                                        eps_phi_j              //
+                                        ) *                    //
+                                       fe_values.JxW(q_point); //
+                }
+
+
+          const PointHistory<dim> *local_quadrature_points_data =
+            reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+
+          body_force.vector_value_list(fe_values.get_quadrature_points(),
+                                       body_force_values);
+
+          for (unsigned int i = 0; i < dofs_per_cell; ++i)
+            {
+              const unsigned int component_i =
+                fe.system_to_component_index(i).first;
+
+              for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
+                {
+                  const SymmetricTensor<2, dim> &old_stress =
+                    local_quadrature_points_data[q_point].old_stress;
+
+                  cell_rhs(i) +=
+                    (body_force_values[q_point](component_i) *
+                       fe_values.shape_value(i, q_point) -
+                     old_stress * get_strain(fe_values, i, q_point)) *
+                    fe_values.JxW(q_point);
+                }
+            }
+
+          cell->get_dof_indices(local_dof_indices);
+
+          hanging_node_constraints.distribute_local_to_global(cell_matrix,
+                                                              cell_rhs,
+                                                              local_dof_indices,
+                                                              system_matrix,
+                                                              system_rhs);
+        }
+
+        // Now compress the vector and the system matrix:
+#ifdef DEAL_II_WITH_PETSC
+    system_matrix.compress(VectorOperation::add);
+    system_rhs.compress(VectorOperation::add);
+#endif
+
+
+    FEValuesExtractors::Scalar                z_component(dim - 1);
+    std::map<types::global_dof_index, double> boundary_values;
+    VectorTools::interpolate_boundary_values(mapping,
+                                             dof_handler,
+                                             0,
+                                             Functions::ZeroFunction<dim>(dim),
+                                             boundary_values);
+    VectorTools::interpolate_boundary_values(
+      mapping,
+      dof_handler,
+      1,
+      IncrementalBoundaryValues<dim>(present_time, present_timestep),
+      boundary_values,
+      fe.component_mask(z_component));
+
+#ifdef DEAL_II_WITH_PETSC
+    PETScWrappers::MPI::Vector tmp(locally_owned_dofs, mpi_communicator);
+#else
+    Vector<double> tmp(dof_handler.n_dofs());
+#endif
+    MatrixTools::apply_boundary_values(
+      boundary_values, system_matrix, tmp, system_rhs, false);
+    incremental_displacement = tmp;
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::solve_timestep()
+  {
+    deallog << "    Assembling system..." << std::flush;
+    assemble_system();
+    deallog << " norm of rhs is " << system_rhs.l2_norm() << std::endl;
+
+    const unsigned int n_iterations = solve_linear_problem();
+
+    deallog << "    Solver converged in " << n_iterations << " iterations."
+            << std::endl;
+
+    deallog << "    Updating quadrature point data..." << std::flush;
+    update_quadrature_point_history();
+    deallog << std::endl;
+  }
+
+
+  template <int dim>
+  unsigned int
+  TopLevel<dim>::solve_linear_problem()
+  {
+#ifdef DEAL_II_WITH_PETSC
+    PETScWrappers::MPI::Vector distributed_incremental_displacement(
+      locally_owned_dofs, mpi_communicator);
+    distributed_incremental_displacement = incremental_displacement;
+#else
+    Vector<double> distributed_incremental_displacement(dof_handler.n_dofs());
+    distributed_incremental_displacement = incremental_displacement;
+#endif
+
+    SolverControl solver_control(dof_handler.n_dofs(),
+                                 1e-16 * system_rhs.l2_norm());
+
+#ifdef DEAL_II_WITH_PETSC
+    PETScWrappers::SolverCG cg(solver_control, mpi_communicator);
+
+    PETScWrappers::PreconditionBlockJacobi preconditioner(system_matrix);
+
+    cg.solve(system_matrix,
+             distributed_incremental_displacement,
+             system_rhs,
+             preconditioner);
+#else
+    SolverCG<Vector<double>> solver(solver_control);
+    solver.solve(system_matrix,
+                 distributed_incremental_displacement,
+                 system_rhs,
+                 PreconditionIdentity());
+#endif
+
+    deallog << "norm: " << distributed_incremental_displacement.linfty_norm()
+            << " " << distributed_incremental_displacement.l1_norm() << " "
+            << distributed_incremental_displacement.l2_norm() << std::endl;
+
+    incremental_displacement = distributed_incremental_displacement;
+
+    hanging_node_constraints.distribute(incremental_displacement);
+
+    return solver_control.last_step();
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::output_results() const
+  {
+    return;
+
+    DataOut<dim> data_out;
+    data_out.attach_dof_handler(dof_handler);
+
+    std::vector<std::string> solution_names;
+    std::vector<DataComponentInterpretation::DataComponentInterpretation>
+      solution_interpretation;
+
+    solution_names.assign(dim, "delta");
+    solution_interpretation.assign(
+      dim,
+      DataComponentInterpretation::DataComponentInterpretation::
+        component_is_part_of_vector);
+
+    data_out.add_data_vector(incremental_displacement,
+                             solution_names,
+                             DataOut_DoFData<DoFHandler<dim, dim>, dim, dim>::
+                               DataVectorType::type_automatic,
+                             solution_interpretation);
+
+
+    Vector<double> norm_of_stress(triangulation.n_active_cells());
+    {
+      // Loop over all the cells...
+      for (auto &cell : triangulation.active_cell_iterators())
+        if (cell->is_locally_owned())
+          {
+            // On these cells, add up the stresses over all quadrature
+            // points...
+            SymmetricTensor<2, dim> accumulated_stress;
+            for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
+              accumulated_stress +=
+                reinterpret_cast<PointHistory<dim> *>(cell->user_pointer())[q]
+                  .old_stress;
+
+            // ...then write the norm of the average to their destination:
+            norm_of_stress(cell->active_cell_index()) =
+              (accumulated_stress / quadrature_formula.size()).norm();
+          }
+        else
+          norm_of_stress(cell->active_cell_index()) = -1e+20;
+    }
+
+    data_out.add_data_vector(norm_of_stress, "norm_of_stress");
+
+    std::vector<types::subdomain_id> partition_int(
+      triangulation.n_active_cells());
+    GridTools::get_subdomain_association(triangulation, partition_int);
+    const Vector<double> partitioning(partition_int.begin(),
+                                      partition_int.end());
+    data_out.add_data_vector(partitioning, "partitioning");
+
+    data_out.build_patches(mapping, 2);
+
+#if false
+    std::ofstream output("step18." + std::to_string(this_mpi_process) + "." +
+                         std::to_string(timestep_no) + ".vtk");
+    data_out.write_vtk(output);
+#endif
+  }
+
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::do_initial_timestep(const bool do_output)
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    deallog << "Timestep " << timestep_no << " at time " << present_time
+            << std::endl;
+
+    for (unsigned int cycle = 0; cycle < 1; ++cycle)
+      {
+        deallog << "  Cycle " << cycle << ':' << std::endl;
+
+        if (cycle == 0)
+          create_coarse_grid();
+        // else
+        //  refine_initial_grid();
+
+        deallog << "    Number of active cells:       "
+                << triangulation.n_active_cells() << " (by partition:";
+        for (unsigned int p = 0; p < n_mpi_processes; ++p)
+          deallog << (p == 0 ? ' ' : '+')
+                  << (GridTools::count_cells_with_subdomain_association(
+                       triangulation, p));
+        deallog << ")" << std::endl;
+
+        setup_system();
+
+        deallog << "    Number of degrees of freedom: " << dof_handler.n_dofs()
+                << " (by partition:";
+        for (unsigned int p = 0; p < n_mpi_processes; ++p)
+          deallog << (p == 0 ? ' ' : '+')
+                  << (DoFTools::count_dofs_with_subdomain_association(
+                       dof_handler, p));
+        deallog << ")" << std::endl;
+
+        solve_timestep();
+      }
+
+    move_mesh();
+
+    if (do_output)
+      output_results();
+
+    deallog << std::endl;
+  }
+
+  template <int dim>
+  void
+  TopLevel<dim>::do_timestep(const bool do_output)
+  {
+    present_time += present_timestep;
+    ++timestep_no;
+    deallog << "Timestep " << timestep_no << " at time " << present_time
+            << std::endl;
+    if (present_time > end_time)
+      {
+        present_timestep -= (present_time - end_time);
+        present_time = end_time;
+      }
+
+
+    solve_timestep();
+
+    move_mesh();
+
+    if (do_output)
+      output_results();
+
+    deallog << std::endl;
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::refine_initial_grid()
+  {
+    // First, let each process compute error indicators for the cells it owns:
+    Vector<float> error_per_cell(triangulation.n_active_cells());
+    KellyErrorEstimator<dim>::estimate(
+      dof_handler,
+      QGauss<dim - 1>(fe.degree + 1),
+      std::map<types::boundary_id, const Function<dim> *>(),
+      incremental_displacement,
+      error_per_cell,
+      ComponentMask(),
+      nullptr,
+      MultithreadInfo::n_threads(),
+      this_mpi_process);
+
+    const unsigned int n_local_cells = triangulation.n_active_cells();
+
+#ifdef DEAL_II_WITH_PETSC
+    PETScWrappers::MPI::Vector distributed_error_per_cell(
+      mpi_communicator, triangulation.n_active_cells(), n_local_cells);
+#else
+    Vector<double> distributed_error_per_cell(n_local_cells);
+#endif
+
+    for (unsigned int i = 0; i < error_per_cell.size(); ++i)
+      if (error_per_cell(i) != 0)
+        distributed_error_per_cell(i) = error_per_cell(i);
+    distributed_error_per_cell.compress(VectorOperation::insert);
+
+    error_per_cell = distributed_error_per_cell;
+    GridRefinement::refine_and_coarsen_fixed_number(triangulation,
+                                                    error_per_cell,
+                                                    0.35,
+                                                    0.03);
+    triangulation.execute_coarsening_and_refinement();
+
+    setup_quadrature_point_history();
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::move_mesh()
+  {
+    deallog << "    Moving mesh..." << std::endl;
+
+    std::vector<bool> vertex_touched(triangulation.n_vertices(), false);
+    for (auto &cell : dof_handler.active_cell_iterators())
+      for (unsigned int v = 0; v < cell->n_vertices(); ++v)
+        if (vertex_touched[cell->vertex_index(v)] == false)
+          {
+            vertex_touched[cell->vertex_index(v)] = true;
+
+            Point<dim> vertex_displacement;
+            for (unsigned int d = 0; d < dim; ++d)
+              vertex_displacement[d] =
+                incremental_displacement(cell->vertex_dof_index(v, d));
+
+            cell->vertex(v) += vertex_displacement;
+          }
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::setup_quadrature_point_history()
+  {
+    triangulation.clear_user_data();
+
+    {
+      std::vector<PointHistory<dim>> tmp;
+      quadrature_point_history.swap(tmp);
+    }
+    quadrature_point_history.resize(triangulation.n_active_cells() *
+                                    quadrature_formula.size());
+
+    unsigned int history_index = 0;
+    for (auto &cell : triangulation.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          cell->set_user_pointer(&quadrature_point_history[history_index]);
+          history_index += quadrature_formula.size();
+        }
+
+    Assert(history_index == quadrature_point_history.size(),
+           ExcInternalError());
+  }
+
+
+  template <int dim>
+  void
+  TopLevel<dim>::update_quadrature_point_history()
+  {
+    FEValues<dim> fe_values(mapping,
+                            fe,
+                            quadrature_formula,
+                            update_values | update_gradients);
+
+    std::vector<std::vector<Tensor<1, dim>>> displacement_increment_grads(
+      quadrature_formula.size(), std::vector<Tensor<1, dim>>(dim));
+
+    for (auto &cell : dof_handler.active_cell_iterators())
+      if (cell->is_locally_owned())
+        {
+          PointHistory<dim> *local_quadrature_points_history =
+            reinterpret_cast<PointHistory<dim> *>(cell->user_pointer());
+          Assert(local_quadrature_points_history >=
+                   &quadrature_point_history.front(),
+                 ExcInternalError());
+          Assert(local_quadrature_points_history <=
+                   &quadrature_point_history.back(),
+                 ExcInternalError());
+
+          fe_values.reinit(cell);
+          fe_values.get_function_gradients(incremental_displacement,
+                                           displacement_increment_grads);
+
+          for (unsigned int q = 0; q < quadrature_formula.size(); ++q)
+            {
+              const SymmetricTensor<2, dim> new_stress =
+                (local_quadrature_points_history[q].old_stress +
+                 (stress_strain_tensor *
+                  get_strain(displacement_increment_grads[q])));
+
+              const Tensor<2, dim> rotation =
+                get_rotation_matrix(displacement_increment_grads[q]);
+
+              const SymmetricTensor<2, dim> rotated_new_stress =
+                symmetrize(transpose(rotation) *
+                           static_cast<Tensor<2, dim>>(new_stress) * rotation);
+
+              local_quadrature_points_history[q].old_stress =
+                rotated_new_stress;
+            }
+        }
+  }
+} // namespace Step18
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+
+  try
+    {
+      using namespace dealii;
+      using namespace Step18;
+
+      Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+      TopLevel<3> elastic_problem;
+      elastic_problem.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+
+  return 0;
+}
diff --git a/tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output b/tests/simplex/step-18.mpirun=1.with_petsc=false.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..0672e69
--- /dev/null
@@ -0,0 +1,94 @@
+
+DEAL::Timestep 1 at time 1.00000
+DEAL::  Cycle 0:
+DEAL::    Number of active cells:       1875 (by partition: 1875)
+DEAL::    Number of degrees of freedom: 1728 (by partition: 1728)
+DEAL::    Assembling system... norm of rhs is 2.45750e+10
+DEAL:cg::Starting value 1.95453e+10
+DEAL:cg::Convergence step 126 value 2.13799e-06
+DEAL::norm: 0.0800000 25.2034 1.12390
+DEAL::    Solver converged in 126 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 2 at time 2.00000
+DEAL::    Assembling system... norm of rhs is 2.49844e+10
+DEAL:cg::Starting value 2.03490e+10
+DEAL:cg::Convergence step 127 value 1.74037e-06
+DEAL::norm: 0.0800000 25.3187 1.12498
+DEAL::    Solver converged in 127 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 3 at time 3.00000
+DEAL::    Assembling system... norm of rhs is 2.54414e+10
+DEAL:cg::Starting value 2.12110e+10
+DEAL:cg::Convergence step 127 value 2.31040e-06
+DEAL::norm: 0.0800000 25.4405 1.12611
+DEAL::    Solver converged in 127 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 4 at time 4.00000
+DEAL::    Assembling system... norm of rhs is 2.59506e+10
+DEAL:cg::Starting value 2.21374e+10
+DEAL:cg::Convergence step 127 value 2.50495e-06
+DEAL::norm: 0.0800000 25.5693 1.12729
+DEAL::    Solver converged in 127 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 5 at time 5.00000
+DEAL::    Assembling system... norm of rhs is 2.65175e+10
+DEAL:cg::Starting value 2.31346e+10
+DEAL:cg::Convergence step 128 value 2.04181e-06
+DEAL::norm: 0.0800000 25.7059 1.12851
+DEAL::    Solver converged in 128 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 6 at time 6.00000
+DEAL::    Assembling system... norm of rhs is 2.71482e+10
+DEAL:cg::Starting value 2.42104e+10
+DEAL:cg::Convergence step 129 value 2.29400e-06
+DEAL::norm: 0.0800000 25.8512 1.12979
+DEAL::    Solver converged in 129 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 7 at time 7.00000
+DEAL::    Assembling system... norm of rhs is 2.78496e+10
+DEAL:cg::Starting value 2.53733e+10
+DEAL:cg::Convergence step 131 value 2.19345e-06
+DEAL::norm: 0.0800000 26.0062 1.13114
+DEAL::    Solver converged in 131 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 8 at time 8.00000
+DEAL::    Assembling system... norm of rhs is 2.86298e+10
+DEAL:cg::Starting value 2.66332e+10
+DEAL:cg::Convergence step 132 value 2.63526e-06
+DEAL::norm: 0.0800000 26.1724 1.13256
+DEAL::    Solver converged in 132 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 9 at time 9.00000
+DEAL::    Assembling system... norm of rhs is 2.94979e+10
+DEAL:cg::Starting value 2.80016e+10
+DEAL:cg::Convergence step 135 value 2.78060e-06
+DEAL::norm: 0.0800000 26.3512 1.13406
+DEAL::    Solver converged in 135 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 10 at time 10.0000
+DEAL::    Assembling system... norm of rhs is 3.04646e+10
+DEAL:cg::Starting value 2.94916e+10
+DEAL:cg::Convergence step 137 value 2.61758e-06
+DEAL::norm: 0.0800000 26.5447 1.13566
+DEAL::    Solver converged in 137 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
diff --git a/tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output b/tests/simplex/step-18.mpirun=1.with_petsc=true.with_simplex_support=on.output
new file mode 100644 (file)
index 0000000..8fb41aa
--- /dev/null
@@ -0,0 +1,94 @@
+
+DEAL::Timestep 1 at time 1.00000
+DEAL::  Cycle 0:
+DEAL::    Number of active cells:       1875 (by partition: 1875)
+DEAL::    Number of degrees of freedom: 1728 (by partition: 1728)
+DEAL::    Assembling system... norm of rhs is 6.42608e+09
+DEAL::Starting value 0.406456
+DEAL::Convergence step 49 value 2.89749e-07
+DEAL::norm: 0.0800000 25.2034 1.12390
+DEAL::    Solver converged in 49 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 2 at time 2.00000
+DEAL::    Assembling system... norm of rhs is 6.38314e+09
+DEAL::Starting value 0.413459
+DEAL::Convergence step 48 value 2.33590e-07
+DEAL::norm: 0.0800000 25.3187 1.12498
+DEAL::    Solver converged in 48 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 3 at time 3.00000
+DEAL::    Assembling system... norm of rhs is 6.34811e+09
+DEAL::Starting value 0.420747
+DEAL::Convergence step 47 value 2.56013e-07
+DEAL::norm: 0.0800000 25.4405 1.12611
+DEAL::    Solver converged in 47 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 4 at time 4.00000
+DEAL::    Assembling system... norm of rhs is 6.32248e+09
+DEAL::Starting value 0.428346
+DEAL::Convergence step 45 value 6.05128e-07
+DEAL::norm: 0.0800000 25.5693 1.12729
+DEAL::    Solver converged in 45 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 5 at time 5.00000
+DEAL::    Assembling system... norm of rhs is 6.30808e+09
+DEAL::Starting value 0.436284
+DEAL::Convergence step 45 value 1.99224e-07
+DEAL::norm: 0.0800000 25.7059 1.12851
+DEAL::    Solver converged in 45 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 6 at time 6.00000
+DEAL::    Assembling system... norm of rhs is 6.30720e+09
+DEAL::Starting value 0.444594
+DEAL::Convergence step 43 value 5.57039e-07
+DEAL::norm: 0.0800000 25.8512 1.12979
+DEAL::    Solver converged in 43 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 7 at time 7.00000
+DEAL::    Assembling system... norm of rhs is 6.32271e+09
+DEAL::Starting value 0.453310
+DEAL::Convergence step 42 value 3.35520e-07
+DEAL::norm: 0.0800000 26.0062 1.13114
+DEAL::    Solver converged in 42 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 8 at time 8.00000
+DEAL::    Assembling system... norm of rhs is 6.35827e+09
+DEAL::Starting value 0.462469
+DEAL::Convergence step 41 value 2.38274e-07
+DEAL::norm: 0.0800000 26.1724 1.13256
+DEAL::    Solver converged in 41 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 9 at time 9.00000
+DEAL::    Assembling system... norm of rhs is 6.41862e+09
+DEAL::Starting value 0.472114
+DEAL::Convergence step 40 value 2.20795e-07
+DEAL::norm: 0.0800000 26.3512 1.13406
+DEAL::    Solver converged in 40 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::
+DEAL::Timestep 10 at time 10.0000
+DEAL::    Assembling system... norm of rhs is 6.50998e+09
+DEAL::Starting value 0.482291
+DEAL::Convergence step 38 value 4.09099e-07
+DEAL::norm: 0.0800000 26.5447 1.13566
+DEAL::    Solver converged in 38 iterations.
+DEAL::    Updating quadrature point data...
+DEAL::    Moving mesh...
+DEAL::

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.