]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Convert step-12 to a simplex test 11457/head
authorPeter Munch <peterrmuench@gmail.com>
Sat, 2 Jan 2021 23:26:40 +0000 (00:26 +0100)
committerPeter Munch <peterrmuench@gmail.com>
Tue, 5 Jan 2021 22:06:01 +0000 (23:06 +0100)
tests/simplex/step-12.cc
tests/simplex/step-12.with_simplex_support=on.output
tests/simplex/step-12a.cc [new file with mode: 0644]
tests/simplex/step-12a.with_simplex_support=on.output [new file with mode: 0644]

index 29484a3aa275befba6168b0a641b17f5e1686417..94d64a085605e9751810a326d3fc70ed0baddebc 100644 (file)
 // ---------------------------------------------------------------------
 
 
-// Step-12 with tetrahedron mesh.
+// Step-02 on a simplex mesh. Following incompatible modifications had to be
+// made:
+//  - The function GridGenerator::hyper_cube() is only working for hypercube
+//    meshes. Use this function to create a temporary mesh and convert it to
+//    simplex mesh.
+//  - Local refinement is replaced by global refinement.
+
 
 #include <deal.II/base/function.h>
 #include <deal.II/base/quadrature_lib.h>
@@ -26,6 +32,7 @@
 #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_fe.h>
 #include <deal.II/fe/mapping_q1.h>
 
 #include <deal.II/grid/grid_generator.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/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>
@@ -62,8 +61,6 @@
 #include <fstream>
 #include <iostream>
 
-//#define HEX
-
 namespace Step12
 {
   using namespace dealii;
@@ -110,20 +107,21 @@ namespace Step12
     wind_field(0) = -p(1);
     wind_field(1) = p(0);
 
-    if (wind_field.norm() > 1e-6)
+    if (wind_field.norm() > 1e-10)
       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 Quadrature<dim> &    quadrature,
+                const Quadrature<dim - 1> &quadrature_face,
                 const UpdateFlags          update_flags = update_values |
                                                  update_gradients |
                                                  update_quadrature_points |
@@ -131,8 +129,11 @@ namespace Step12
                 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)
+      : fe_values(mapping, fe, quadrature, update_flags)
+      , fe_interface_values(mapping,
+                            fe,
+                            quadrature_face,
+                            interface_update_flags)
     {}
 
 
@@ -141,12 +142,10 @@ namespace Step12
                   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())
+      , fe_interface_values(scratch_data.fe_interface_values.get_mapping(),
+                            scratch_data.fe_interface_values.get_fe(),
+                            scratch_data.fe_interface_values.get_quadrature(),
+                            scratch_data.fe_interface_values.get_update_flags())
     {}
 
     FEValues<dim>          fe_values;
@@ -203,21 +202,14 @@ namespace Step12
     void
     output_results(const unsigned int cycle) const;
 
-    Triangulation<dim> triangulation;
-#ifdef HEX
-    const MappingQ1<dim> mapping;
-#else
-    Simplex::FE_P<dim>   fe_mapping;
+    Triangulation<dim>   triangulation;
     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;
+    const Simplex::FE_DGP<dim> fe;
+    DoFHandler<dim>            dof_handler;
+
+    const Simplex::QGauss<dim>     quadrature;
+    const Simplex::QGauss<dim - 1> quadrature_face;
 
     SparsityPattern      sparsity_pattern;
     SparseMatrix<double> system_matrix;
@@ -229,15 +221,11 @@ namespace Step12
 
   template <int dim>
   AdvectionProblem<dim>::AdvectionProblem()
-#ifdef HEX
-    : mapping()
-    , fe(2)
-#else
-    : fe_mapping(1)
-    , mapping(fe_mapping)
-    , fe(2)
-#endif
+    : mapping(Simplex::FE_P<dim>(1))
+    , fe(1)
     , dof_handler(triangulation)
+    , quadrature(fe.tensor_degree() + 1)
+    , quadrature_face(fe.tensor_degree() + 1)
   {}
 
 
@@ -263,10 +251,11 @@ namespace Step12
     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;
+    const auto cell_worker = [&](const Iterator &  cell,
+                                 ScratchData<dim> &scratch_data,
+                                 CopyData &        copy_data) {
+      const unsigned int n_dofs =
+        scratch_data.fe_values.get_fe().n_dofs_per_cell();
       copy_data.reinit(cell, n_dofs);
       scratch_data.fe_values.reinit(cell);
 
@@ -290,10 +279,10 @@ namespace Step12
         }
     };
 
-    auto boundary_worker = [&](const Iterator &    cell,
-                               const unsigned int &face_no,
-                               ScratchData<dim> &  scratch_data,
-                               CopyData &          copy_data) {
+    const 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);
@@ -330,14 +319,14 @@ namespace Step12
         }
     };
 
-    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) {
+    const 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();
@@ -367,9 +356,9 @@ namespace Step12
         }
     };
 
-    AffineConstraints<double> constraints;
+    const AffineConstraints<double> constraints;
 
-    auto copier = [&](const CopyData &c) {
+    const auto copier = [&](const CopyData &c) {
       constraints.distribute_local_to_global(c.cell_matrix,
                                              c.cell_rhs,
                                              c.local_dof_indices,
@@ -384,19 +373,7 @@ namespace Step12
         }
     };
 
-    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::QGauss<dim> quad(degree + 1);
-
-    Simplex::QGauss<dim - 1> face_quad(degree + 1);
-#endif
-
-    ScratchData<dim> scratch_data(mapping, fe, quad, face_quad);
+    ScratchData<dim> scratch_data(mapping, fe, quadrature, quadrature_face);
     CopyData         copy_data;
 
     MeshWorker::mesh_loop(dof_handler.begin_active(),
@@ -421,12 +398,12 @@ namespace Step12
 
     PreconditionBlockSSOR<SparseMatrix<double>> preconditioner;
 
-    preconditioner.initialize(system_matrix, fe.dofs_per_cell);
+    preconditioner.initialize(system_matrix, fe.n_dofs_per_cell());
 
     solver.solve(system_matrix, solution, right_hand_side, preconditioner);
 
-    deallog << "  Solver converged in " << solver_control.last_step()
-            << " iterations." << std::endl;
+    std::cout << "  Solver converged in " << solver_control.last_step()
+              << " iterations." << std::endl;
   }
 
 
@@ -434,24 +411,7 @@ namespace Step12
   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();
+    triangulation.refine_global();
   }
 
 
@@ -459,27 +419,17 @@ namespace Step12
   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;
+    const std::string filename = "solution-" + std::to_string(cycle) + ".vtk";
+    std::cout << "  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.build_patches(mapping);
 
     data_out.write_vtk(output);
-#endif
 
     {
       Vector<float> values(triangulation.n_active_cells());
@@ -488,13 +438,13 @@ namespace Step12
                                         solution,
                                         Functions::ZeroFunction<dim>(),
                                         values,
-                                        QGauss<dim>(fe.degree + 1),
+                                        quadrature,
                                         VectorTools::Linfty_norm);
       const double l_infty =
         VectorTools::compute_global_error(triangulation,
                                           values,
                                           VectorTools::Linfty_norm);
-      deallog << "  L-infinity norm: " << l_infty << std::endl;
+      std::cout << "  L-infinity norm: " << l_infty << std::endl;
     }
   }
 
@@ -503,36 +453,27 @@ namespace Step12
   void
   AdvectionProblem<dim>::run()
   {
-    //#ifdef HEX
-    //    for (unsigned int cycle = 0; cycle < 6; ++cycle)
-    //#else
-    for (unsigned int cycle = 0; cycle < 1; ++cycle)
-      //#endif
+    for (unsigned int cycle = 0; cycle < 4; ++cycle)
       {
-        deallog << "Cycle " << cycle << std::endl;
+        std::cout << "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
+            Triangulation<dim> temp;
+            GridGenerator::hyper_cube(temp);
+            GridGenerator::convert_hypercube_to_simplex_mesh(temp,
+                                                             triangulation);
           }
         else
           refine_grid();
 
-        deallog << "  Number of active cells:       "
-                << triangulation.n_active_cells() << std::endl;
+        std::cout << "  Number of active cells:       "
+                  << triangulation.n_active_cells() << std::endl;
 
         setup_system();
 
-        deallog << "  Number of degrees of freedom: " << dof_handler.n_dofs()
-                << std::endl;
+        std::cout << "  Number of degrees of freedom: " << dof_handler.n_dofs()
+                  << std::endl;
 
         assemble_system();
         solve();
@@ -546,16 +487,36 @@ namespace Step12
 int
 main()
 {
-  initlog();
-
-  {
-    Step12::AdvectionProblem<2> dgmethod;
-    dgmethod.run();
-  }
-  {
-    Step12::AdvectionProblem<3> dgmethod;
-    dgmethod.run();
-  }
+  try
+    {
+      Step12::AdvectionProblem<2> dgmethod;
+      dgmethod.run();
+    }
+  catch (std::exception &exc)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Exception on processing: " << std::endl
+                << exc.what() << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
+  catch (...)
+    {
+      std::cerr << std::endl
+                << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      std::cerr << "Unknown exception!" << std::endl
+                << "Aborting!" << std::endl
+                << "----------------------------------------------------"
+                << std::endl;
+      return 1;
+    }
 
   return 0;
 }
index 7fbe5eb0f73c1d1509cc6a9351f22300b593e771..f6149cf94475aab37b29bbdd0037e6cbf7b5a980 100644 (file)
@@ -1,15 +1,24 @@
-
-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
+Cycle 0
+  Number of active cells:       8
+  Number of degrees of freedom: 24
+  Solver converged in 7 iterations.
+  Writing solution to <solution-0.vtk>
+  L-infinity norm: 1.00799
+Cycle 1
+  Number of active cells:       32
+  Number of degrees of freedom: 96
+  Solver converged in 8 iterations.
+  Writing solution to <solution-1.vtk>
+  L-infinity norm: 1.14371
+Cycle 2
+  Number of active cells:       128
+  Number of degrees of freedom: 384
+  Solver converged in 11 iterations.
+  Writing solution to <solution-2.vtk>
+  L-infinity norm: 1.09928
+Cycle 3
+  Number of active cells:       512
+  Number of degrees of freedom: 1536
+  Solver converged in 21 iterations.
+  Writing solution to <solution-3.vtk>
+  L-infinity norm: 1.11521
diff --git a/tests/simplex/step-12a.cc b/tests/simplex/step-12a.cc
new file mode 100644 (file)
index 0000000..29484a3
--- /dev/null
@@ -0,0 +1,561 @@
+// ---------------------------------------------------------------------
+//
+// 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::QGauss<dim> quad(degree + 1);
+
+    Simplex::QGauss<dim - 1> face_quad(degree + 1);
+#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-12a.with_simplex_support=on.output b/tests/simplex/step-12a.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

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.