From 6d8dc7b85fb8b3776a59831aa6b542ed8d3874c0 Mon Sep 17 00:00:00 2001
From: Jiaqi Zhang <jiaqi2@clemson.edu>
Date: Tue, 16 Mar 2021 17:43:27 +0000
Subject: [PATCH] add some descriptions

---
 tests/simplex/step-67.cc | 69 ++++++++++------------------------------
 1 file changed, 17 insertions(+), 52 deletions(-)

diff --git a/tests/simplex/step-67.cc b/tests/simplex/step-67.cc
index 1a72141570..a4c08cc59f 100644
--- a/tests/simplex/step-67.cc
+++ b/tests/simplex/step-67.cc
@@ -13,11 +13,24 @@
  *
  * ---------------------------------------------------------------------
  *
- * Author: Martin Kronbichler, 2020
+ * step-67 on a simplex mesh (only testcase = 0 is tested). The major
+ * modifications are:
+ * - Since MatrixFreeOperators::CellwiseInverseMassMatrix cannot be used, the
+ *   following two places are changed:
+ *   - EulerOperator::project() is replaced by VectorTools::project().
+ *   - The linear system is solved globally by first using
+ * MatrixFreeOperators::MassOperator to generate the matrix and then using CG to
+ * solve it, instead of using EulerOperator::local_apply_inverse_mass_matrix()
+ *
+ * - The function GridGenerator::hyper_rectangle() is only working for hypercube
+ * meshes. Use it to create a temporary mesh and then convert it to a simplex
+ * mesh.
+ *
+ *
  */
 
-// The include files are similar to the previous matrix-free tutorial programs
-// step-37, step-48, and step-59
+
+
 //#define HEX
 
 #include <deal.II/base/conditional_ostream.h>
@@ -801,13 +814,6 @@ namespace Euler_DG
       const LinearAlgebra::distributed::Vector<Number> &src,
       const std::pair<unsigned int, unsigned int> &     cell_range) const;
 
-    void
-    local_apply_vmult(
-      const MatrixFree<dim, Number> &                   data,
-      LinearAlgebra::distributed::Vector<Number> &      dst,
-      const LinearAlgebra::distributed::Vector<Number> &src,
-      const std::pair<unsigned int, unsigned int> &     cell_range) const;
-
     void
     local_apply_face(
       const MatrixFree<dim, Number> &                   data,
@@ -1123,33 +1129,6 @@ namespace Euler_DG
       }
   }
 
-  template <int dim, int degree, int n_points_1d>
-  void
-  EulerOperator<dim, degree, n_points_1d>::local_apply_vmult(
-    const MatrixFree<dim, Number> &,
-    LinearAlgebra::distributed::Vector<Number> &      dst,
-    const LinearAlgebra::distributed::Vector<Number> &src,
-    const std::pair<unsigned int, unsigned int> &     cell_range) const
-  {
-#ifdef HEX
-    FEEvaluation<dim, degree, n_points_1d, dim + 2, Number> phi(data);
-#else
-    FEEvaluation<dim, -1, n_points_1d, dim + 2, Number>     phi(data);
-#endif
-
-    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
-      {
-        phi.reinit(cell);
-        phi.gather_evaluate(src, EvaluationFlags::values);
-
-        for (unsigned int q = 0; q < phi.n_q_points; ++q)
-          {
-            phi.submit_value(phi.get_value(q), q);
-          }
-        phi.integrate_scatter(EvaluationFlags::values, dst);
-      }
-  }
-
   // The next function concerns the computation of integrals on interior
   // faces, where we need evaluators from both cells adjacent to the face. We
   // associate the variable `phi_m` with the solution component $\mathbf{w}^-$
@@ -1486,19 +1465,6 @@ namespace Euler_DG
     }
   }
 
-  template <int dim, int degree, int n_points_1d>
-  void
-  EulerOperator<dim, degree, n_points_1d>::vmult(
-    LinearAlgebra::distributed::Vector<Number> &      dst,
-    const LinearAlgebra::distributed::Vector<Number> &src) const
-  {
-    {
-      TimerOutput::Scope t(timer, "apply - vmult");
-
-      data.cell_loop(&EulerOperator::local_apply_vmult, this, dst, src, true);
-    }
-  }
-
   // Let us move to the function that does an entire stage of a Runge--Kutta
   // update. It calls EulerOperator::apply() followed by some updates
   // to the vectors, namely `next_ri = solution + factor_ai * k_i` and
@@ -1614,12 +1580,11 @@ namespace Euler_DG
       std::shared_ptr<MatrixFree<dim, Number>> matrix_free(
         new MatrixFree<dim, Number>(data));
       mass_matrix.initialize(matrix_free);
-      mass_matrix.compute_diagonal();
 
       ReductionControl control(6 * vec_ki.size(), 0., 1e-12, false, false);
       SolverCG<LinearAlgebra::distributed::Vector<Number>> cg(control);
 
-      cg.solve(*this, next_ri, vec_ki, PreconditionIdentity());
+      cg.solve(mass_matrix, next_ri, vec_ki, PreconditionIdentity());
 
       const Number ai = factor_ai;
       const Number bi = factor_solution;
-- 
2.39.5