Modernize step-21.

1. Ranged for-loops
2. Improve discussion relative to step-20
3. Use DGQ, not DQ, for standard discontinuous elements
4. inline function definitions

diff --git a/examples/step-21/doc/intro.dox b/examples/step-21/doc/intro.dox
index 90b63a47ca380c662e3418c1295db8f617a77d56..4c1f2c8d9611a2368b46ee34f0e153675f732ddb 100644 (file)
--- a/examples/step-21/doc/intro.dox
+++ b/examples/step-21/doc/intro.dox
@@ -238,7 +238,7 @@ cell term to get an equation as follows:
 In each time step, we then apply the mixed finite method of @ref step_20
 "step-20" to the velocity and pressure. To be well-posed, we choose
 Raviart-Thomas spaces $RT_{k}$ for $\mathbf{u}$ and discontinuous elements of
-class $DQ_{k}$ for $p$. For the saturation, we will also choose $DQ_{k}$
+class $DGQ_{k}$ for $p$. For the saturation, we will also choose $DGQ_{k}$
 spaces.
 
 Since we have discontinuous spaces, we have to think about how to evaluate
@@ -267,8 +267,8 @@ fluxes can also be found in step-12.
 <h3>Linear solvers</h3>
 
 The linear solvers used in this program are a straightforward extension of the
-ones used in step-20. Essentially, we simply have to extend
-everything from
+ones used in step-20 (but without LinearOperator). Essentially, we simply have
+to extend everything from
 two to three solution components. If we use the discrete spaces
 mentioned above and put shape functions into the bilinear forms, we
 arrive at the following linear system to be solved for time step $n+1$:
@@ -294,7 +294,7 @@ B &    0 & 0\\
 @f]
 where the individual matrices and vectors are defined as follows using
 shape functions $\mathbf v_i$ (of type Raviart Thomas $RT_k$) for
-velocities and $\phi_i$ (of type $DG_k$) for both pressures and saturations:
+velocities and $\phi_i$ (of type $DGQ_k$) for both pressures and saturations:
 @f{eqnarray*}
 M^u(S^n)_{ij} &=&
 \left((\mathbf{K}\lambda(S^n))^{-1} \mathbf{v}_i,\mathbf

diff --git a/examples/step-21/doc/results.dox b/examples/step-21/doc/results.dox
index cf2fb3232d4216e4a63c359f70759bff7d096329..66345340db832c7e77514956275f5bf40f4f4251 100644 (file)
--- a/examples/step-21/doc/results.dox
+++ b/examples/step-21/doc/results.dox
@@ -42,7 +42,7 @@ the pressure once is significantly reduced.
 The final observation concerns the number of iterations needed to solve for
 the saturation, i.e. one. This shouldn't surprise us too much: the matrix we
 have to solve with is the mass matrix. However, this is the mass matrix for
-the $DQ_0$ element of piecewise constants where no element couples with the
+the $DGQ_0$ element of piecewise constants where no element couples with the
 degrees of freedom on neighboring cells. The matrix is therefore a diagonal
 one, and it is clear that we should be able to invert this matrix in a single
 CG iteration.

diff --git a/examples/step-21/step-21.cc b/examples/step-21/step-21.cc
index 4b97da08e573ecf9f9c08955a676518f2984a65f..897de3cadf680fedac8be2dedf9fd8a373578c6e 100644 (file)
--- a/examples/step-21/step-21.cc
+++ b/examples/step-21/step-21.cc
@@ -144,21 +144,15 @@ namespace Step21
       : Function<dim>(1)
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+    virtual double value(const Point<dim> & /*p*/,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      return 0;
+    }
   };
 
 
 
-  template <int dim>
-  double
-  PressureRightHandSide<dim>::value(const Point<dim> & /*p*/,
-                                    const unsigned int /*component*/) const
-  {
-    return 0;
-  }
-
-
   // @sect4{Pressure boundary values}
 
   // The next are pressure boundary values. As mentioned in the introduction,
@@ -171,19 +165,14 @@ namespace Step21
       : Function<dim>(1)
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      return 1 - p[0];
+    }
   };
 
 
-  template <int dim>
-  double
-  PressureBoundaryValues<dim>::value(const Point<dim> &p,
-                                     const unsigned int /*component*/) const
-  {
-    return 1 - p[0];
-  }
-
 
   // @sect4{Saturation boundary values}
 
@@ -200,25 +189,18 @@ namespace Step21
       : Function<dim>(1)
     {}
 
-    virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+    virtual double value(const Point<dim> &p,
+                         const unsigned int /*component*/ = 0) const override
+    {
+      if (p[0] == 0)
+        return 1;
+      else
+        return 0;
+    }
   };
 
 
 
-  template <int dim>
-  double
-  SaturationBoundaryValues<dim>::value(const Point<dim> &p,
-                                       const unsigned int /*component*/) const
-  {
-    if (p[0] == 0)
-      return 1;
-    else
-      return 0;
-  }
-
-
-
   // @sect4{Initial data}
 
   // Finally, we need initial data. In reality, we only need initial data for
@@ -241,29 +223,19 @@ namespace Step21
     {}
 
     virtual double value(const Point<dim> & p,
-                         const unsigned int component = 0) const override;
+                         const unsigned int component = 0) const override
+    {
+      return Functions::ZeroFunction<dim>(dim + 2).value(p, component);
+    }
 
     virtual void vector_value(const Point<dim> &p,
-                              Vector<double> &  value) const override;
+                              Vector<double> &  values) const override
+    {
+      Functions::ZeroFunction<dim>(dim + 2).vector_value(p, values);
+    }
   };
 
 
-  template <int dim>
-  double InitialValues<dim>::value(const Point<dim> & p,
-                                   const unsigned int component) const
-  {
-    return Functions::ZeroFunction<dim>(dim + 2).value(p, component);
-  }
-
-
-  template <int dim>
-  void InitialValues<dim>::vector_value(const Point<dim> &p,
-                                        Vector<double> &  values) const
-  {
-    Functions::ZeroFunction<dim>(dim + 2).vector_value(p, values);
-  }
-
-
 
   // @sect3{The inverse permeability tensor}
 
@@ -289,34 +261,30 @@ namespace Step21
         : TensorFunction<2, dim>()
       {}
 
-      virtual void value_list(const std::vector<Point<dim>> &points,
-                              std::vector<Tensor<2, dim>> &  values) const;
-    };
-
-
-    template <int dim>
-    void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
-                                   std::vector<Tensor<2, dim>> &  values) const
-    {
-      Assert(points.size() == values.size(),
-             ExcDimensionMismatch(points.size(), values.size()));
+      virtual void
+      value_list(const std::vector<Point<dim>> &points,
+                 std::vector<Tensor<2, dim>> &  values) const override
+      {
+        Assert(points.size() == values.size(),
+               ExcDimensionMismatch(points.size(), values.size()));
 
-      for (unsigned int p = 0; p < points.size(); ++p)
-        {
-          values[p].clear();
+        for (unsigned int p = 0; p < points.size(); ++p)
+          {
+            values[p].clear();
 
-          const double distance_to_flowline =
-            std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
+            const double distance_to_flowline =
+              std::fabs(points[p][1] - 0.5 - 0.1 * std::sin(10 * points[p][0]));
 
-          const double permeability =
-            std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
-                              (0.1 * 0.1)),
-                     0.01);
+            const double permeability =
+              std::max(std::exp(-(distance_to_flowline * distance_to_flowline) /
+                                (0.1 * 0.1)),
+                       0.01);
 
-          for (unsigned int d = 0; d < dim; ++d)
-            values[p][d][d] = 1. / permeability;
-        }
-    }
+            for (unsigned int d = 0; d < dim; ++d)
+              values[p][d][d] = 1. / permeability;
+          }
+      }
+    };
   } // namespace SingleCurvingCrack
 
 
@@ -364,60 +332,50 @@ namespace Step21
 
       virtual void
       value_list(const std::vector<Point<dim>> &points,
-                 std::vector<Tensor<2, dim>> &  values) const override;
-
-    private:
-      static std::vector<Point<dim>> centers;
+                 std::vector<Tensor<2, dim>> &  values) const override
+      {
+        Assert(points.size() == values.size(),
+               ExcDimensionMismatch(points.size(), values.size()));
 
-      static std::vector<Point<dim>> get_centers();
-    };
+        for (unsigned int p = 0; p < points.size(); ++p)
+          {
+            values[p].clear();
 
+            double permeability = 0;
+            for (unsigned int i = 0; i < centers.size(); ++i)
+              permeability += std::exp(-(points[p] - centers[i]).norm_square() /
+                                       (0.05 * 0.05));
 
+            const double normalized_permeability =
+              std::min(std::max(permeability, 0.01), 4.);
 
-    template <int dim>
-    std::vector<Point<dim>>
-      KInverse<dim>::centers = KInverse<dim>::get_centers();
+            for (unsigned int d = 0; d < dim; ++d)
+              values[p][d][d] = 1. / normalized_permeability;
+          }
+      }
 
+    private:
+      static std::vector<Point<dim>> centers;
 
-    template <int dim>
-    std::vector<Point<dim>> KInverse<dim>::get_centers()
-    {
-      const unsigned int N =
-        (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
+      static std::vector<Point<dim>> get_centers()
+      {
+        const unsigned int N =
+          (dim == 2 ? 40 : (dim == 3 ? 100 : throw ExcNotImplemented()));
 
-      std::vector<Point<dim>> centers_list(N);
-      for (unsigned int i = 0; i < N; ++i)
-        for (unsigned int d = 0; d < dim; ++d)
-          centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
+        std::vector<Point<dim>> centers_list(N);
+        for (unsigned int i = 0; i < N; ++i)
+          for (unsigned int d = 0; d < dim; ++d)
+            centers_list[i][d] = static_cast<double>(rand()) / RAND_MAX;
 
-      return centers_list;
-    }
+        return centers_list;
+      }
+    };
 
 
 
     template <int dim>
-    void KInverse<dim>::value_list(const std::vector<Point<dim>> &points,
-                                   std::vector<Tensor<2, dim>> &  values) const
-    {
-      Assert(points.size() == values.size(),
-             ExcDimensionMismatch(points.size(), values.size()));
-
-      for (unsigned int p = 0; p < points.size(); ++p)
-        {
-          values[p].clear();
-
-          double permeability = 0;
-          for (unsigned int i = 0; i < centers.size(); ++i)
-            permeability +=
-              std::exp(-(points[p] - centers[i]).norm_square() / (0.05 * 0.05));
-
-          const double normalized_permeability =
-            std::min(std::max(permeability, 0.01), 4.);
-
-          for (unsigned int d = 0; d < dim; ++d)
-            values[p][d][d] = 1. / normalized_permeability;
-        }
-    }
+    std::vector<Point<dim>>
+      KInverse<dim>::centers = KInverse<dim>::get_centers();
   } // namespace RandomMedium
 
 
@@ -458,44 +416,44 @@ namespace Step21
   class InverseMatrix : public Subscriptor
   {
   public:
-    InverseMatrix(const MatrixType &m);
+    InverseMatrix(const MatrixType &m)
+      : matrix(&m)
+    {}
 
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+    void vmult(Vector<double> &dst, const Vector<double> &src) const
+    {
+      SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
+                                   1e-8 * src.l2_norm());
+      SolverCG<>    cg(solver_control);
+
+      dst = 0;
+
+      cg.solve(*matrix, dst, src, PreconditionIdentity());
+    }
 
   private:
     const SmartPointer<const MatrixType> matrix;
   };
 
 
-  template <class MatrixType>
-  InverseMatrix<MatrixType>::InverseMatrix(const MatrixType &m)
-    : matrix(&m)
-  {}
-
-
-
-  template <class MatrixType>
-  void InverseMatrix<MatrixType>::vmult(Vector<double> &      dst,
-                                        const Vector<double> &src) const
-  {
-    SolverControl solver_control(std::max<unsigned int>(src.size(), 200),
-                                 1e-8 * src.l2_norm());
-    SolverCG<>    cg(solver_control);
-
-    dst = 0;
-
-    cg.solve(*matrix, dst, src, PreconditionIdentity());
-  }
-
-
 
   class SchurComplement : public Subscriptor
   {
   public:
     SchurComplement(const BlockSparseMatrix<double> &          A,
-                    const InverseMatrix<SparseMatrix<double>> &Minv);
+                    const InverseMatrix<SparseMatrix<double>> &Minv)
+      : system_matrix(&A)
+      , m_inverse(&Minv)
+      , tmp1(A.block(0, 0).m())
+      , tmp2(A.block(0, 0).m())
+    {}
 
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+    void vmult(Vector<double> &dst, const Vector<double> &src) const
+    {
+      system_matrix->block(0, 1).vmult(tmp1, src);
+      m_inverse->vmult(tmp2, tmp1);
+      system_matrix->block(1, 0).vmult(dst, tmp2);
+    }
 
   private:
     const SmartPointer<const BlockSparseMatrix<double>>           system_matrix;
@@ -506,32 +464,21 @@ namespace Step21
 
 
 
-  SchurComplement::SchurComplement(
-    const BlockSparseMatrix<double> &          A,
-    const InverseMatrix<SparseMatrix<double>> &Minv)
-    : system_matrix(&A)
-    , m_inverse(&Minv)
-    , tmp1(A.block(0, 0).m())
-    , tmp2(A.block(0, 0).m())
-  {}
-
-
-  void SchurComplement::vmult(Vector<double> &      dst,
-                              const Vector<double> &src) const
-  {
-    system_matrix->block(0, 1).vmult(tmp1, src);
-    m_inverse->vmult(tmp2, tmp1);
-    system_matrix->block(1, 0).vmult(dst, tmp2);
-  }
-
-
-
   class ApproximateSchurComplement : public Subscriptor
   {
   public:
-    ApproximateSchurComplement(const BlockSparseMatrix<double> &A);
+    ApproximateSchurComplement(const BlockSparseMatrix<double> &A)
+      : system_matrix(&A)
+      , tmp1(A.block(0, 0).m())
+      , tmp2(A.block(0, 0).m())
+    {}
 
-    void vmult(Vector<double> &dst, const Vector<double> &src) const;
+    void vmult(Vector<double> &dst, const Vector<double> &src) const
+    {
+      system_matrix->block(0, 1).vmult(tmp1, src);
+      system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
+      system_matrix->block(1, 0).vmult(dst, tmp2);
+    }
 
   private:
     const SmartPointer<const BlockSparseMatrix<double>> system_matrix;
@@ -540,23 +487,6 @@ namespace Step21
   };
 
 
-  ApproximateSchurComplement::ApproximateSchurComplement(
-    const BlockSparseMatrix<double> &A)
-    : system_matrix(&A)
-    , tmp1(A.block(0, 0).m())
-    , tmp2(A.block(0, 0).m())
-  {}
-
-
-  void ApproximateSchurComplement::vmult(Vector<double> &      dst,
-                                         const Vector<double> &src) const
-  {
-    system_matrix->block(0, 1).vmult(tmp1, src);
-    system_matrix->block(0, 0).precondition_Jacobi(tmp2, tmp1);
-    system_matrix->block(1, 0).vmult(dst, tmp2);
-  }
-
-
 
   // @sect3{<code>TwoPhaseFlowProblem</code> class implementation}
 
@@ -720,10 +650,7 @@ namespace Step21
     const FEValuesExtractors::Scalar pressure(dim);
     const FEValuesExtractors::Scalar saturation(dim + 1);
 
-    typename DoFHandler<dim>::active_cell_iterator cell =
-                                                     dof_handler.begin_active(),
-                                                   endc = dof_handler.end();
-    for (; cell != endc; ++cell)
+    for (const auto &cell : dof_handler.active_cell_iterators())
       {
         fe_values.reinit(cell);
         local_matrix = 0;
@@ -878,10 +805,7 @@ namespace Step21
 
     const FEValuesExtractors::Scalar saturation(dim + 1);
 
-    typename DoFHandler<dim>::active_cell_iterator cell =
-                                                     dof_handler.begin_active(),
-                                                   endc = dof_handler.end();
-    for (; cell != endc; ++cell)
+    for (const auto &cell : dof_handler.active_cell_iterators())
       {
         local_rhs = 0;
         fe_values.reinit(cell);
@@ -936,8 +860,7 @@ namespace Step21
                 fe_face_values.get_quadrature_points(), neighbor_saturation);
             else
               {
-                const typename DoFHandler<dim>::active_cell_iterator neighbor =
-                  cell->neighbor(face_no);
+                const auto         neighbor = cell->neighbor(face_no);
                 const unsigned int neighbor_face =
                   cell->neighbor_of_neighbor(face_no);
 
@@ -1093,18 +1016,11 @@ namespace Step21
     switch (dim)
       {
         case 2:
-          solution_names.emplace_back("u");
-          solution_names.emplace_back("v");
-          solution_names.emplace_back("p");
-          solution_names.emplace_back("S");
+          solution_names = {"u", "v", "p", "S"};
           break;
 
         case 3:
-          solution_names.emplace_back("u");
-          solution_names.emplace_back("v");
-          solution_names.emplace_back("w");
-          solution_names.emplace_back("p");
-          solution_names.emplace_back("S");
+          solution_names = {"u", "v", "w", "p", "S"};
           break;
 
         default:
@@ -1169,10 +1085,7 @@ namespace Step21
                                                 Vector<double>(dim + 2));
     double                      max_velocity = 0;
 
-    typename DoFHandler<dim>::active_cell_iterator cell =
-                                                     dof_handler.begin_active(),
-                                                   endc = dof_handler.end();
-    for (; cell != endc; ++cell)
+    for (const auto &cell : dof_handler.active_cell_iterators())
       {
         fe_values.reinit(cell);
         fe_values.get_function_values(solution, solution_values);