From: Martin Kronbichler <kronbichler@lnm.mw.tum.de>
Date: Tue, 16 Apr 2013 15:56:15 +0000 (+0000)
Subject: Add test for mass matrix evaluation in matrix-free for some elements. Avoid deprecate... 
X-Git-Tag: v8.0.0~683
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=eb0f9bfddd0c049a6183b930c6799ffed1c20097;p=dealii.git

Add test for mass matrix evaluation in matrix-free for some elements. Avoid deprecated function.

git-svn-id: https://svn.dealii.org/trunk@29305 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/tests/matrix_free/compress_mapping.cc b/tests/matrix_free/compress_mapping.cc
index 35a2c51191..d4e3c08091 100644
--- a/tests/matrix_free/compress_mapping.cc
+++ b/tests/matrix_free/compress_mapping.cc
@@ -121,62 +121,19 @@ void test_cube ()
 
 
 
-void create_parallelogram(Triangulation<2> &tria)
-{
-  const int dim = 2;
-  std::vector<Point<dim> > points (4);
-  points[0] = Point<dim> (0, 0);
-  points[1] = Point<dim> (0, 1);
-  points[2] = Point<dim> (1 ,0.5);
-  points[3] = Point<dim> (1 ,1.5);
-
-  std::vector<CellData<dim> > cells(1);
-  cells[0].vertices[0] = 0;
-  cells[0].vertices[1] = 2;
-  cells[0].vertices[2] = 1;
-  cells[0].vertices[3] = 3;
-  cells[0].material_id = 0;
-
-  tria.create_triangulation (points, cells, SubCellData());
-}
-
-
-
-void create_parallelogram(Triangulation<3> &tria)
-{
-  const int dim = 3;
-  std::vector<Point<dim> > points (8);
-  points[0] = Point<dim> (0,0,0);
-  points[1] = Point<dim> (0,1.,0.5);
-  points[2] = Point<dim> (0,0,1);
-  points[3] = Point<dim> (0,1.,1.5);
-  points[4] = Point<dim> (1.,0,1.);
-  points[5] = Point<dim> (1.,1.,1.5);
-  points[6] = Point<dim> (1.,0,2);
-  points[7] = Point<dim> (1.,1.,2.5);
-
-  std::vector<CellData<dim> > cells(1);
-  cells[0].vertices[0] = 0;
-  cells[0].vertices[1] = 4;
-  cells[0].vertices[2] = 1;
-  cells[0].vertices[3] = 5;
-  cells[0].vertices[4] = 2;
-  cells[0].vertices[5] = 6;
-  cells[0].vertices[6] = 3;
-  cells[0].vertices[7] = 7;
-  cells[0].material_id = 0;
-
-  tria.create_triangulation (points, cells, SubCellData());
-}
-
-
-
 template <int dim>
 void test_parallelogram ()
 {
   deallog << "Parallelogram" << std::endl;
   Triangulation<dim> tria;
-  create_parallelogram(tria);
+  Point<dim> corners[dim];
+  for (unsigned int d=0; d<dim; ++d)
+    {
+      corners[d][d] = 1.;
+      if (d>0)
+        corners[d][0] = 0.5+0.5*d;
+    }
+  GridGenerator::parallelepiped(tria, corners);
   tria.refine_global(5-dim);
 
   FE_Q<dim> fe (1);
diff --git a/tests/matrix_free/estimate_condition_number_mass.cc b/tests/matrix_free/estimate_condition_number_mass.cc
new file mode 100644
index 0000000000..0a78638d37
--- /dev/null
+++ b/tests/matrix_free/estimate_condition_number_mass.cc
@@ -0,0 +1,156 @@
+//------------------  estimate_condition_number_mass.cc  ---------------------
+//    $Id$
+//    Version: $Name$
+//
+//------------------  estimate_condition_number_mass.cc  ---------------------
+
+
+// this function computes condition number estimates for the mass matrix at
+// different polynomial degrees. The mesh uses a hypercube mesh with no
+// hanging nodes and no other constraints
+
+#include "../tests.h"
+
+std::ofstream logfile("estimate_condition_number_mass/output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q_hierarchical.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/base/function_lib.h>
+
+
+template <int dim, int fe_degree, typename Number>
+void
+mass_operator (const MatrixFree<dim,Number>  &data,
+               Vector<Number>       &dst,
+               const Vector<Number> &src,
+               const std::pair<unsigned int,unsigned int> &cell_range)
+{
+  FEEvaluationGeneral<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+  const unsigned int n_q_points = fe_eval.n_q_points;
+
+  for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
+    {
+      fe_eval.reinit (cell);
+      fe_eval.read_dof_values (src);
+      fe_eval.template evaluate (true, false, false);
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          fe_eval.submit_value (fe_eval.get_value(q),q);
+        }
+      fe_eval.template integrate (true,false);
+      fe_eval.distribute_local_to_global (dst);
+    }
+}
+
+
+
+template <int dim, int fe_degree, typename Number>
+class MatrixFreeTest
+{
+ public:
+  typedef VectorizedArray<Number> vector_t;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void vmult (Vector<Number>       &dst,
+              const Vector<Number> &src) const
+  {
+    dst = 0;
+    const std_cxx1x::function<void(const MatrixFree<dim,Number>  &,
+                                   Vector<Number>       &,
+                                   const Vector<Number> &,
+                                   const std::pair<unsigned int,unsigned int>&)>
+      wrap = mass_operator<dim,fe_degree,Number>;
+    data.cell_loop (wrap, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void test (const FiniteElement<dim> &fe)
+{
+  typedef double number;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube (tria);
+  tria.refine_global(2);
+
+  DoFHandler<dim> dof(tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  constraints.close();
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+1);
+    mf_data.reinit (dof, constraints, quad);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+
+  VectorTools::create_right_hand_side(dof, QGauss<dim>(fe_degree+1),
+                                      Functions::CosineFunction<dim>(), in);
+
+  SolverControl control(10000, 1e-9*in.l2_norm());
+  typename SolverCG<>::AdditionalData data;
+  data.compute_condition_number = true;
+  SolverCG<> solver(control,data);
+  solver.solve(mf, out, in, PreconditionIdentity());
+  deallog << std::endl;
+}
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (2);
+
+  {
+    deallog.threshold_double(1.e-9);
+    deallog.push("2d");
+    test<2,1>(FE_Q<2>(1));
+    test<2,1>(FE_DGQ<2>(1));
+    test<2,2>(FE_Q<2>(2));
+    test<2,2>(FE_DGQ<2>(2));
+    test<2,4>(FE_Q<2>(4));
+    test<2,4>(FE_Q<2>(QGaussLobatto<1>(5)));
+    test<2,4>(FE_DGQ<2>(4));
+    test<2,4>(FE_Q_Hierarchical<2>(4));
+    test<2,6>(FE_Q<2>(6));
+    test<2,6>(FE_Q<2>(QGaussLobatto<1>(7)));
+    test<2,6>(FE_DGQ<2>(6));
+    test<2,6>(FE_DGQArbitraryNodes<2>(QGaussLobatto<1>(7)));
+    test<2,10>(FE_Q<2>(10));
+    test<2,10>(FE_Q<2>(QGaussLobatto<1>(11)));
+    test<2,16>(FE_Q<2>(QGaussLobatto<1>(17)));
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>(FE_Q<3>(1));
+    test<3,2>(FE_Q<3>(2));
+    test<3,5>(FE_Q<3>(5));
+    test<3,5>(FE_Q<3>(QGaussLobatto<1>(6)));
+    deallog.pop();
+  }
+}
+
+
diff --git a/tests/matrix_free/estimate_condition_number_mass/cmp/generic b/tests/matrix_free/estimate_condition_number_mass/cmp/generic
new file mode 100644
index 0000000000..9a991fe035
--- /dev/null
+++ b/tests/matrix_free/estimate_condition_number_mass/cmp/generic
@@ -0,0 +1,90 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d:cg::Starting value 0.117
+DEAL:2d:cg::Convergence step 16 value 0
+DEAL:2d:cg::Condition number estimate: 15.4
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(1)
+DEAL:2d:cg::Starting value 0.0625
+DEAL:2d:cg::Convergence step 1 value 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d:cg::Starting value 0.0686
+DEAL:2d:cg::Convergence step 38 value 0
+DEAL:2d:cg::Condition number estimate: 29.1
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(2)
+DEAL:2d:cg::Starting value 0.0625
+DEAL:2d:cg::Convergence step 3 value 0
+DEAL:2d:cg::Condition number estimate: 6.95
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d:cg::Starting value 0.0367
+DEAL:2d:cg::Convergence step 92 value 0
+DEAL:2d:cg::Condition number estimate: 148.
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(QGaussLobatto(5))
+DEAL:2d:cg::Starting value 0.0355
+DEAL:2d:cg::Convergence step 62 value 0
+DEAL:2d:cg::Condition number estimate: 74.3
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(4)
+DEAL:2d:cg::Starting value 0.0353
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d:cg::Condition number estimate: 200.
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(6)
+DEAL:2d:cg::Starting value 0.0310
+DEAL:2d:cg::Convergence step 209 value 0
+DEAL:2d:cg::Condition number estimate: 1.56e+03
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(QGaussLobatto(7))
+DEAL:2d:cg::Starting value 0.0241
+DEAL:2d:cg::Convergence step 86 value 0
+DEAL:2d:cg::Condition number estimate: 140.
+DEAL:2d::
+DEAL:2d::Testing FE_DGQ<2>(6)
+DEAL:2d:cg::Starting value 0.0305
+DEAL:2d:cg::Convergence step 15 value 0
+DEAL:2d:cg::Condition number estimate: 1.95e+03
+DEAL:2d::
+DEAL:2d::Testing FE_DGQArbitraryNodes<2>(QGaussLobatto(7))
+DEAL:2d:cg::Starting value 0.0240
+DEAL:2d:cg::Convergence step 12 value 0
+DEAL:2d:cg::Condition number estimate: 185.
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(10)
+DEAL:2d:cg::Starting value 0.173
+DEAL:2d:cg::Convergence step 1305 value 0
+DEAL:2d:cg::Condition number estimate: 2.21e+06
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(QGaussLobatto(11))
+DEAL:2d:cg::Starting value 0.0148
+DEAL:2d:cg::Convergence step 132 value 0
+DEAL:2d:cg::Condition number estimate: 328.
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(QGaussLobatto(17))
+DEAL:2d:cg::Starting value 0.00937
+DEAL:2d:cg::Convergence step 200 value 0
+DEAL:2d:cg::Condition number estimate: 759.
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d:cg::Starting value 0.102
+DEAL:3d:cg::Convergence step 12 value 0
+DEAL:3d:cg::Condition number estimate: 57.8
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d:cg::Starting value 0.0498
+DEAL:3d:cg::Convergence step 51 value 0
+DEAL:3d:cg::Condition number estimate: 161.
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(5)
+DEAL:3d:cg::Starting value 0.0121
+DEAL:3d:cg::Convergence step 475 value 0
+DEAL:3d:cg::Condition number estimate: 8.96e+03
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(QGaussLobatto(6))
+DEAL:3d:cg::Starting value 0.0137
+DEAL:3d:cg::Convergence step 167 value 0
+DEAL:3d:cg::Condition number estimate: 1.06e+03
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_05.cc b/tests/matrix_free/matrix_vector_05.cc
index 838bf0790f..481c55899d 100644
--- a/tests/matrix_free/matrix_vector_05.cc
+++ b/tests/matrix_free/matrix_vector_05.cc
@@ -20,15 +20,18 @@ std::ofstream logfile("matrix_vector_05/output");
 template <int dim, int fe_degree>
 void test ()
 {
-  if (dim == 3)
-    return;
   Triangulation<dim> tria;
-  Tensor<2,dim> points;
+  Point<dim> points[dim];
   points[0][0] = 0.25;
   points[0][1] = 0.123;
   points[1][0] = 0.09983712334;
   points[1][1] = 0.314159265358979;
-  GridGenerator::parallelogram (tria, points);
+  if (dim == 3)
+    {
+      //points[2][0] = 0.21;
+      //points[2][2] = 0.4123;
+    }
+  GridGenerator::parallelepiped (tria, points);
   typename Triangulation<dim>::active_cell_iterator
     cell = tria.begin_active (),
     endc = tria.end();
@@ -41,13 +44,12 @@ void test ()
     if (cell->center().norm()<0.2)
       cell->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  if (dim < 3 || fe_degree < 2)
+  if (dim < 3)
     tria.refine_global(2);
   tria.begin(tria.n_levels()-1)->set_refine_flag();
   tria.last()->set_refine_flag();
   tria.execute_coarsening_and_refinement();
-  cell = tria.begin_active ();
-  for (unsigned int i=0; i<5; ++i)
+  for (int i=0; i<7-2*fe_degree; ++i)
     {
       cell = tria.begin_active ();
       unsigned int counter = 0;
@@ -56,6 +58,7 @@ void test ()
           cell->set_refine_flag();
       tria.execute_coarsening_and_refinement();
     }
+  std::cout << dim << " " << tria.n_active_cells() << std::endl;
 
   FE_Q<dim> fe (fe_degree);
   DoFHandler<dim> dof (tria);
diff --git a/tests/matrix_free/matrix_vector_05/cmp/generic b/tests/matrix_free/matrix_vector_05/cmp/generic
index f2943c08e9..c8d1cd1e9f 100644
--- a/tests/matrix_free/matrix_vector_05/cmp/generic
+++ b/tests/matrix_free/matrix_vector_05/cmp/generic
@@ -5,3 +5,9 @@ DEAL:2d::
 DEAL:2d::Testing FE_Q<2>(2)
 DEAL:2d::Norm of difference: 0
 DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_common.h b/tests/matrix_free/matrix_vector_common.h
index f0518dc904..63f9199c3e 100644
--- a/tests/matrix_free/matrix_vector_common.h
+++ b/tests/matrix_free/matrix_vector_common.h
@@ -52,9 +52,6 @@ helmholtz_operator (const MatrixFree<dim,Number>  &data,
   for(unsigned int cell=cell_range.first;cell<cell_range.second;++cell)
     {
       fe_eval.reinit (cell);
-
-                                // compare values with the ones the FEValues
-                                // gives us. Those are seen as reference
       fe_eval.read_dof_values (src);
       fe_eval.template evaluate (true, true, false);
       for (unsigned int q=0; q<n_q_points; ++q)