]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Tests for FEEvaluation without template on polynomial degree
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 16 Feb 2017 11:18:03 +0000 (12:18 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 17 Feb 2017 12:26:55 +0000 (13:26 +0100)
include/deal.II/matrix_free/fe_evaluation.h
tests/matrix_free/get_functions_notempl.cc [new file with mode: 0644]
tests/matrix_free/get_functions_notempl.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_06_notempl.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_06_notempl.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_common.h
tests/matrix_free/matrix_vector_stokes_notempl.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_stokes_notempl.output [new file with mode: 0644]

index 7a95e3c88130e43825d995e8a4630e1c73457b97..3ff41dcc6534667c32bcc5230383f83300bd9b1d 100644 (file)
@@ -7404,9 +7404,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 #ifdef DEBUG
   // print error message when the dimensions do not match. Propose a possible
   // fix
-  if (fe_degree != -1 && static_cast<unsigned int>(fe_degree) != this->data->fe_degree
+  if ((fe_degree != -1 && static_cast<unsigned int>(fe_degree) != this->data->fe_degree)
       ||
-      n_q_points != this->data->n_q_points)
+      (fe_degree != -1 && static_n_q_points != this->data->n_q_points))
     {
       std::string message =
         "-------------------------------------------------------\n";
@@ -7439,13 +7439,13 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
                   proposed_dof_comp = no;
                   break;
                 }
-          if (n_q_points ==
+          if (static_n_q_points ==
               this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
             proposed_quad_comp = this->quad_no;
           else
             for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
               if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
-                  == n_q_points)
+                  == static_n_q_points)
                 {
                   proposed_quad_comp = no;
                   break;
@@ -7507,7 +7507,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       message += "                                 " + correct_pos;
 
       Assert (static_cast<unsigned int>(fe_degree) == this->data->fe_degree &&
-              n_q_points == this->data->n_q_points,
+              static_n_q_points == this->data->n_q_points,
               ExcMessage(message));
     }
   if (fe_no != numbers::invalid_unsigned_int)
diff --git a/tests/matrix_free/get_functions_notempl.cc b/tests/matrix_free/get_functions_notempl.cc
new file mode 100644 (file)
index 0000000..418cd4a
--- /dev/null
@@ -0,0 +1,115 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// this function tests the correctness of the implementation of matrix free
+// operations in getting the function values, the function gradients, and the
+// function Laplacians on a hyperball mesh with different sizes in the number
+// of degrees of freedom per cell and quadrature points per cell. This is the
+// same test as get_functions_rect but without using a template parameter on
+// the degree for FEEvaluation.
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include "get_functions_common.h"
+
+
+template <int dim, int fe_degree, int n_q_points_1d, typename number>
+void sub_test (const DoFHandler<dim> &dof,
+               const ConstraintMatrix &constraints,
+               MatrixFree<dim,number> &mf_data,
+               Vector<number> &solution)
+{
+  deallog << "Test with fe_degree " << fe_degree
+          << ", n_q_points_1d: " << (n_q_points_1d) << std::endl;
+  const QGauss<1> quad (n_q_points_1d);
+  MappingQ<dim> mapping (2);
+  typename MatrixFree<dim,number>::AdditionalData data;
+  data.tasks_parallel_scheme = MatrixFree<dim,number>::AdditionalData::none;
+  data.mapping_update_flags = update_gradients | update_hessians;
+  mf_data.reinit (mapping, dof, constraints, quad, data);
+  MatrixFreeTest<dim,-1,n_q_points_1d,number> mf (mf_data, mapping);
+  mf.test_functions (solution);
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  typedef double number;
+  const SphericalManifold<dim> manifold;
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  typename Triangulation<dim>::active_cell_iterator
+  cell = tria.begin_active (),
+  endc = tria.end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+      if (cell->at_boundary(f))
+        cell->face(f)->set_all_manifold_ids(0);
+  tria.set_manifold (0, manifold);
+
+  // refine first and last cell
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global (1);
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints (dof, constraints);
+  constraints.close();
+
+
+  // in the other functions, use do_test in
+  // get_functions_common, but here we have to
+  // manually choose non-rectangular tests.
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  //std::cout << "Number of cells: " << dof.get_triangulation().n_active_cells()
+  //          << std::endl;
+  //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
+  //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
+
+  Vector<number> solution (dof.n_dofs());
+
+  // create vector with random entries
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      solution(i) = entry;
+    }
+  constraints.distribute(solution);
+
+  MatrixFree<dim,number> mf_data;
+  if (fe_degree > 1)
+    sub_test <dim,fe_degree,fe_degree-1,number> (dof, constraints, mf_data,
+                                                 solution);
+  sub_test <dim,fe_degree,fe_degree,number> (dof, constraints, mf_data,
+                                             solution);
+  sub_test <dim,fe_degree,fe_degree+2,number> (dof, constraints, mf_data,
+                                               solution);
+  if (dim == 2)
+    sub_test <dim,fe_degree,fe_degree+3,number> (dof, constraints, mf_data,
+                                                 solution);
+}
diff --git a/tests/matrix_free/get_functions_notempl.output b/tests/matrix_free/get_functions_notempl.output
new file mode 100644 (file)
index 0000000..bcea3ca
--- /dev/null
@@ -0,0 +1,147 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Test with fe_degree 1, n_q_points_1d: 1
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 1, n_q_points_1d: 3
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 1, n_q_points_1d: 4
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Test with fe_degree 2, n_q_points_1d: 1
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 2, n_q_points_1d: 2
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 2, n_q_points_1d: 4
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 2, n_q_points_1d: 5
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(3)
+DEAL:2d::Test with fe_degree 3, n_q_points_1d: 2
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 3, n_q_points_1d: 3
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 3, n_q_points_1d: 5
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 3, n_q_points_1d: 6
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Test with fe_degree 4, n_q_points_1d: 3
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 4, n_q_points_1d: 4
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 4, n_q_points_1d: 6
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:2d::Test with fe_degree 4, n_q_points_1d: 7
+DEAL:2d::Error function values: 0
+DEAL:2d::Error function gradients: 0
+DEAL:2d::Error function Laplacians: 0
+DEAL:2d::Error function diagonal of Hessian: 0
+DEAL:2d::Error function Hessians: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::Test with fe_degree 1, n_q_points_1d: 1
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Test with fe_degree 1, n_q_points_1d: 3
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::Test with fe_degree 2, n_q_points_1d: 1
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Test with fe_degree 2, n_q_points_1d: 2
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
+DEAL:3d::Test with fe_degree 2, n_q_points_1d: 4
+DEAL:3d::Error function values: 0
+DEAL:3d::Error function gradients: 0
+DEAL:3d::Error function Laplacians: 0
+DEAL:3d::Error function diagonal of Hessian: 0
+DEAL:3d::Error function Hessians: 0
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_06_notempl.cc b/tests/matrix_free/matrix_vector_06_notempl.cc
new file mode 100644 (file)
index 0000000..1772e8d
--- /dev/null
@@ -0,0 +1,71 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same as matrix_vector_06 but without setting a template parameter on the
+// polynomial degree, which requires to read it as a run time parameter in a
+// separate code path in FEEvaluation
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include "create_mesh.h"
+
+std::ofstream logfile("output");
+
+#include "matrix_vector_common.h"
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  create_mesh (tria);
+  tria.begin_active ()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator cell, endc;
+  cell = tria.begin_active ();
+  endc = tria.end();
+  for (; cell!=endc; ++cell)
+    if (cell->center().norm()<0.5)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  tria.refine_global(1);
+  cell = tria.begin_active ();
+  for (unsigned int i=0; i<10-3*dim; ++i)
+    {
+      cell = tria.begin_active ();
+      endc = tria.end();
+      unsigned int counter = 0;
+      for (; cell!=endc; ++cell, ++counter)
+        if (counter % (7-i) == 0)
+          cell->set_refine_flag();
+      tria.execute_coarsening_and_refinement();
+    }
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+  DoFTools::make_hanging_node_constraints(dof, constraints);
+  VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction<dim>(),
+                                            constraints);
+  constraints.close();
+
+  // do not enable templates in FEEvaluation
+  do_test<dim, -1, double, fe_degree+1> (dof, constraints);
+}
diff --git a/tests/matrix_free/matrix_vector_06_notempl.output b/tests/matrix_free/matrix_vector_06_notempl.output
new file mode 100644 (file)
index 0000000..c8d1cd1
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+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::
index 87c97352c3f362cd0aad17020746258af3ba1e91..1837813c7844c1876a1b49a8fbce687fe6bc5d2e 100644 (file)
@@ -61,7 +61,7 @@ void do_test (const DoFHandler<dim> &dof,
   //std::cout << "Number of degrees of freedom: " << dof.n_dofs() << std::endl;
   //std::cout << "Number of constraints: " << constraints.n_constraints() << std::endl;
 
-  MappingQGeneric<dim> mapping(fe_degree);
+  MappingQGeneric<dim> mapping(dof.get_fe().degree);
   MatrixFree<dim,number> mf_data;
   {
     const QGauss<1> quad (n_q_points_1d);
diff --git a/tests/matrix_free/matrix_vector_stokes_notempl.cc b/tests/matrix_free/matrix_vector_stokes_notempl.cc
new file mode 100644 (file)
index 0000000..0e17648
--- /dev/null
@@ -0,0 +1,347 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// same as matrix_vector_stokes_noflux but no template parameter on the
+// polynomial degree
+
+#include "../tests.h"
+
+std::ofstream logfile("output");
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/utilities.h>
+#include <deal.II/lac/block_vector.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_renumbering.h>
+#include <deal.II/lac/constraint_matrix.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_sparsity_pattern.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/fe/mapping_q.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+
+
+template <int dim, typename VectorType>
+class MatrixFreeTest
+{
+public:
+  typedef typename DoFHandler<dim>::active_cell_iterator CellIterator;
+  typedef double Number;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_apply (const MatrixFree<dim,Number> &data,
+               VectorType          &dst,
+               const VectorType    &src,
+               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    typedef VectorizedArray<Number> vector_t;
+    FEEvaluation<dim,-1,0,dim,Number> velocity (data, 0);
+    FEEvaluation<dim,-1,0,1,  Number> pressure (data, 1);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        velocity.reinit (cell);
+        velocity.read_dof_values (src.block(0));
+        velocity.evaluate (false,true,false);
+        pressure.reinit (cell);
+        pressure.read_dof_values (src.block(1));
+        pressure.evaluate (true,false,false);
+
+        for (unsigned int q=0; q<velocity.n_q_points; ++q)
+          {
+            SymmetricTensor<2,dim,vector_t> sym_grad_u =
+              velocity.get_symmetric_gradient (q);
+            vector_t pres = pressure.get_value(q);
+            vector_t div = -trace(sym_grad_u);
+            pressure.submit_value   (div, q);
+
+            // subtract p * I
+            for (unsigned int d=0; d<dim; ++d)
+              sym_grad_u[d][d] -= pres;
+
+            velocity.submit_symmetric_gradient(sym_grad_u, q);
+          }
+
+        velocity.integrate (false,true);
+        velocity.distribute_local_to_global (dst.block(0));
+        pressure.integrate (true,false);
+        pressure.distribute_local_to_global (dst.block(1));
+      }
+  }
+
+
+  void vmult (VectorType &dst,
+              const VectorType &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,VectorType>::local_apply,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim>
+void test (const unsigned int fe_degree)
+{
+  SphericalManifold<dim> manifold;
+  HyperShellBoundary<dim> boundary;
+  Triangulation<dim>   triangulation;
+  GridGenerator::hyper_shell (triangulation, Point<dim>(),
+                              0.5, 1., 96, true);
+  triangulation.set_all_manifold_ids(0);
+  triangulation.set_all_manifold_ids_on_boundary(1);
+  triangulation.set_manifold (0, manifold);
+  triangulation.set_manifold (1, boundary);
+
+  triangulation.begin_active()->set_refine_flag();
+  triangulation.last()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+  triangulation.refine_global (3-dim);
+  triangulation.last()->set_refine_flag();
+  triangulation.execute_coarsening_and_refinement();
+
+  MappingQ<dim>        mapping (3);
+  FE_Q<dim>            fe_u_scal (fe_degree+1);
+  FESystem<dim>        fe_u (fe_u_scal,dim);
+  FE_Q<dim>            fe_p (fe_degree);
+  FESystem<dim>        fe (fe_u_scal, dim, fe_p, 1);
+  DoFHandler<dim>      dof_handler_u (triangulation);
+  DoFHandler<dim>      dof_handler_p (triangulation);
+  DoFHandler<dim>      dof_handler (triangulation);
+
+  MatrixFree<dim,double> mf_data;
+
+  ConstraintMatrix     constraints, constraints_u, constraints_p;
+
+  BlockSparsityPattern      sparsity_pattern;
+  BlockSparseMatrix<double> system_matrix;
+
+  BlockVector<double> solution;
+  BlockVector<double> system_rhs;
+  BlockVector<double> mf_solution;
+
+  dof_handler.distribute_dofs (fe);
+  dof_handler_u.distribute_dofs (fe_u);
+  dof_handler_p.distribute_dofs (fe_p);
+  std::vector<unsigned int> stokes_sub_blocks (dim+1,0);
+  stokes_sub_blocks[dim] = 1;
+  DoFRenumbering::component_wise (dof_handler, stokes_sub_blocks);
+
+  std::set<types::boundary_id> no_normal_flux_boundaries;
+  no_normal_flux_boundaries.insert (0);
+  no_normal_flux_boundaries.insert (1);
+  DoFTools::make_hanging_node_constraints (dof_handler,
+                                           constraints);
+  VectorTools::compute_no_normal_flux_constraints (dof_handler, 0,
+                                                   no_normal_flux_boundaries,
+                                                   constraints, mapping);
+  constraints.close ();
+  DoFTools::make_hanging_node_constraints (dof_handler_u,
+                                           constraints_u);
+  VectorTools::compute_no_normal_flux_constraints (dof_handler_u, 0,
+                                                   no_normal_flux_boundaries,
+                                                   constraints_u, mapping);
+  constraints_u.close ();
+  DoFTools::make_hanging_node_constraints (dof_handler_p,
+                                           constraints_p);
+  constraints_p.close ();
+
+  std::vector<types::global_dof_index> dofs_per_block (2);
+  DoFTools::count_dofs_per_block (dof_handler, dofs_per_block,
+                                  stokes_sub_blocks);
+
+  //std::cout << "Number of active cells: "
+  //          << triangulation.n_active_cells()
+  //          << std::endl
+  //          << "Number of degrees of freedom: "
+  //          << dof_handler.n_dofs()
+  //          << " (" << n_u << '+' << n_p << ')'
+  //          << std::endl;
+
+  {
+    BlockDynamicSparsityPattern csp (2,2);
+
+    for (unsigned int d=0; d<2; ++d)
+      for (unsigned int e=0; e<2; ++e)
+        csp.block(d,e).reinit (dofs_per_block[d], dofs_per_block[e]);
+
+    csp.collect_sizes();
+
+    DoFTools::make_sparsity_pattern (dof_handler, csp, constraints, false);
+    sparsity_pattern.copy_from (csp);
+  }
+
+  system_matrix.reinit (sparsity_pattern);
+
+  // this is from step-22
+  {
+    QGauss<dim>   quadrature_formula(fe_degree+2);
+
+    FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                             update_values    |
+                             update_JxW_values |
+                             update_gradients);
+
+    const unsigned int   dofs_per_cell   = fe.dofs_per_cell;
+    const unsigned int   n_q_points      = quadrature_formula.size();
+
+    FullMatrix<double>   local_matrix (dofs_per_cell, dofs_per_cell);
+
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+
+    std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+    std::vector<double>                  div_phi_u   (dofs_per_cell);
+    std::vector<double>                  phi_p       (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof_handler.begin_active(),
+    endc = dof_handler.end();
+    for (; cell!=endc; ++cell)
+      {
+        fe_values.reinit (cell);
+        local_matrix = 0;
+
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            for (unsigned int k=0; k<dofs_per_cell; ++k)
+              {
+                phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                div_phi_u[k]   = fe_values[velocities].divergence (k, q);
+                phi_p[k]       = fe_values[pressure].value (k, q);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                for (unsigned int j=0; j<=i; ++j)
+                  {
+                    local_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+                                          - div_phi_u[i] * phi_p[j]
+                                          - phi_p[i] * div_phi_u[j])
+                                         * fe_values.JxW(q);
+                  }
+              }
+          }
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          for (unsigned int j=i+1; j<dofs_per_cell; ++j)
+            local_matrix(i,j) = local_matrix(j,i);
+
+        cell->get_dof_indices (local_dof_indices);
+        constraints.distribute_local_to_global (local_matrix,
+                                                local_dof_indices,
+                                                system_matrix);
+      }
+  }
+
+
+  solution.reinit (2);
+  for (unsigned int d=0; d<2; ++d)
+    solution.block(d).reinit (dofs_per_block[d]);
+  solution.collect_sizes ();
+
+  system_rhs.reinit (solution);
+  mf_solution.reinit (solution);
+
+  // fill system_rhs with random numbers
+  for (unsigned int j=0; j<system_rhs.block(0).size(); ++j)
+    if (constraints_u.is_constrained(j) == false)
+      {
+        const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+        system_rhs.block(0)(j) = val;
+      }
+  for (unsigned int j=0; j<system_rhs.block(1).size(); ++j)
+    if (constraints_p.is_constrained(j) == false)
+      {
+        const double val = -1 + 2.*(double)Testing::rand()/double(RAND_MAX);
+        system_rhs.block(1)(j) = val;
+      }
+
+  // setup matrix-free structure
+  {
+    std::vector<const DoFHandler<dim>*> dofs;
+    dofs.push_back(&dof_handler_u);
+    dofs.push_back(&dof_handler_p);
+    std::vector<const ConstraintMatrix *> constraints;
+    constraints.push_back (&constraints_u);
+    constraints.push_back (&constraints_p);
+    QGauss<1> quad(fe_degree+2);
+    // no parallelism
+    mf_data.reinit (mapping, dofs, constraints, quad,
+                    typename MatrixFree<dim>::AdditionalData
+                    (MPI_COMM_WORLD,
+                     MatrixFree<dim>::AdditionalData::none));
+  }
+
+  system_matrix.vmult (solution, system_rhs);
+
+  MatrixFreeTest<dim,BlockVector<double> > mf (mf_data);
+  mf.vmult (mf_solution, system_rhs);
+
+  // Verification
+  mf_solution -= solution;
+  const double error = mf_solution.linfty_norm();
+  const double relative = solution.linfty_norm();
+  deallog << "Verification fe degree " << fe_degree  <<  ": "
+          << error/relative << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog << std::endl << "Test with doubles" << std::endl << std::endl;
+    deallog.threshold_double(1.e-12);
+    deallog.push("2d");
+    test<2>(1);
+    test<2>(2);
+    test<2>(3);
+    deallog.pop();
+    deallog.push("3d");
+    test<3>(1);
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_stokes_notempl.output b/tests/matrix_free/matrix_vector_stokes_notempl.output
new file mode 100644 (file)
index 0000000..0434059
--- /dev/null
@@ -0,0 +1,12 @@
+
+DEAL::
+DEAL::Test with doubles
+DEAL::
+DEAL:2d::Verification fe degree 1: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 2: 0
+DEAL:2d::
+DEAL:2d::Verification fe degree 3: 0
+DEAL:2d::
+DEAL:3d::Verification fe degree 1: 0
+DEAL:3d::

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.