]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement FEEvaluation for FE_Q_DG0 elements
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 13 Jun 2014 12:01:42 +0000 (12:01 +0000)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 13 Jun 2014 12:01:42 +0000 (12:01 +0000)
git-svn-id: https://svn.dealii.org/trunk@33040 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/base/aligned_vector.h
deal.II/include/deal.II/fe/fe_poly.templates.h
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/matrix_free/shape_info.h
deal.II/include/deal.II/matrix_free/shape_info.templates.h
tests/matrix_free/matrix_vector_17.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_17.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_stokes_qdg0.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_stokes_qdg0.output [new file with mode: 0644]

index 84a45349c66a90b21f1548451311419283f81443..de4f3e3f36d86d7c7233ae591985823976c5b99e 100644 (file)
@@ -148,6 +148,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: There is now a class FEEvaluationQ_DG0 that does
+  optimized matrix-free evaluation for FE_Q_DG0 elements.
+  <br>
+  (Martin Kronbichler, 2014/06/13)
+  </li>
+
   <li> Bugfix: Filter libclang_rt* from the PETSc link line.
   <br>
   (Matthias Maier, 2014/06/04)
index 88576b5acba941580b6e26d5dfba85eaed12af5d..f963e8e8944ddcb4f824135d000ad9e5bb8ee437 100644 (file)
@@ -487,6 +487,7 @@ AlignedVector<T>::resize_fast (const size_type size)
 }
 
 
+
 template < class T >
 inline
 void
index bfcc19d3acb74793da4d22c01b232b9549e22901..6bdbed94009a65e45b652d39168a9512f27ca5bb 100644 (file)
@@ -18,6 +18,7 @@
 #include <deal.II/base/qprojector.h>
 #include <deal.II/base/polynomial_space.h>
 #include <deal.II/base/tensor_product_polynomials.h>
+#include <deal.II/base/tensor_product_polynomials_const.h>
 #include <deal.II/fe/fe_values.h>
 #include <deal.II/fe/fe_poly.h>
 
@@ -480,6 +481,22 @@ namespace internal
   {
     return poly.get_numbering_inverse();
   }
+
+  template <int dim>
+  inline
+  std::vector<unsigned int>
+  get_poly_space_numbering (const TensorProductPolynomialsConst<dim> &poly)
+  {
+    return poly.get_numbering();
+  }
+
+  template <int dim>
+  inline
+  std::vector<unsigned int>
+  get_poly_space_numbering_inverse (const TensorProductPolynomialsConst<dim> &poly)
+  {
+    return poly.get_numbering_inverse();
+  }
 }
 
 
index 02295fe3223ed68719541bc610878beb4202ca0d..c00b23b0b95d8b90fdb983a934031f2324be2c66 100644 (file)
@@ -1394,13 +1394,25 @@ protected:
   /**
    * Internally stored variables for the different data fields.
    */
-  VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+(dim*dim+2*dim+1)*n_q_points)];
+  VectorizedArray<Number> my_data_array[n_components*(dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)];
+
+  /**
+   * Checks the number of cell dofs with the value that we expect
+   */
+  void check_dofs_per_cell (const unsigned int dofs_per_cell) const;
 
 private:
   /**
    * Sets the pointers from the data array to values_dof, etc.
    */
   void set_data_pointers();
+
+  /**
+   * Checks that the element number of degrees of freedom given by the
+   * template arguments (via fe_degree) coincides with the number of degrees
+   * of freedom in the stored shape values.
+   */
+  void check_template_arguments(const unsigned int fe_no);
 };
 
 
@@ -1576,7 +1588,7 @@ protected:
 
 private:
   /**
-   * Fills the fields shapve_???_evenodd, called in the constructor.
+   * Fills the fields shape_???_evenodd, called in the constructor.
    */
   void compute_even_odd_factors();
 };
@@ -1828,6 +1840,98 @@ public:
 };
 
 
+
+/**
+ * The class that provides all functions necessary to evaluate functions at
+ * quadrature points and cell integrations. In functionality, this class is
+ * similar to FEValues<dim>, however, it includes a lot of specialized
+ * functions that make it much faster (between 5 and 500 times as fast,
+ * depending on the polynomial order). Access to the data fields is provided
+ * through functionality in the class FEEvaluationAccess.
+ *
+ * This class is an extension of FEEvaluation to work with continuous elements
+ * supplemented with a single discontinuous degree of freedom, i.e., FE_Q_DG0.
+ *
+ * @author Martin Kronbichler, 2014
+ */
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+          int n_components_ = 1, typename Number = double >
+class FEEvaluationQ_DG0 :
+  public FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+{
+public:
+  typedef FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number> BaseClass;
+  typedef Number                            number_type;
+  typedef typename BaseClass::value_type    value_type;
+  typedef typename BaseClass::gradient_type gradient_type;
+  static const unsigned int dimension     = dim;
+  static const unsigned int n_components  = n_components_;
+  static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell + 1;
+  static const unsigned int n_q_points    = BaseClass::n_q_points;
+
+  /**
+   * Constructor. Takes all data stored in MatrixFree. If applied to problems
+   * with more than one finite element or more than one quadrature formula
+   * selected during construction of @p matrix_free, @p fe_no and @p quad_no
+   * allow to select the appropriate components.
+   */
+  FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &matrix_free,
+                     const unsigned int            fe_no   = 0,
+                     const unsigned int            quad_no = 0);
+
+  /**
+   * Constructor that comes with reduced functionality and works similar as
+   * FEValues. The user has to provide a structure of type MappingFEEvaluation
+   * and a DoFHandler in order to allow for reading out the finite element
+   * data. It uses the data provided by dof_handler.get_fe(). If the element
+   * is vector-valued, the optional argument allows to specify the index of
+   * the base element (as long as the element is primitive, non-primitive are
+   * not supported currently).
+   *
+   * With initialization from a FEValues object, no call to a reinit method of
+   * this class is necessary. Instead, it is enough if the geometry is
+   * initialized to a given cell iterator. It can also read from or write to
+   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
+   * types (which is less efficient with MPI since index translation has to be
+   * done), but of course only for one cell at a time. Hence, a kernel using
+   * this method does not vectorize over several elements (which is most
+   * efficient for vector operations), but only possibly within the element if
+   * the evaluate/integrate routines are combined (e.g. for computing cell
+   * matrices).
+   */
+  FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
+                     const DoFHandler<dim>                 &dof_handler,
+                     const unsigned int                     first_selected_component = 0);
+
+  /**
+   * Copy constructor
+   */
+  FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other);
+
+  /**
+   * Evaluates the function values, the gradients, and the Hessians of the FE
+   * function given at the DoF values in the input vector at the quadrature
+   * points of the unit cell. The function arguments specify which parts shall
+   * actually be computed. Needs to be called before the functions @p
+   * get_value(), @p get_gradient() or @p get_laplacian give useful
+   * information (unless these values have been set manually).
+   */
+  void evaluate (const bool evaluate_val,
+                 const bool evaluate_grad,
+                 const bool evaluate_lapl = false);
+
+  /**
+   * This function takes the values and/or gradients that are stored on
+   * quadrature points, tests them by all the basis functions/gradients on the
+   * cell and performs the cell integration. The two function arguments @p
+   * integrate_val and @p integrate_grad are used to enable/disable some of
+   * values or gradients.
+   */
+  void integrate (const bool integrate_val,
+                  const bool integrate_grad);
+};
+
+
 /*----------------------- Inline functions ----------------------------------*/
 
 #ifndef DOXYGEN
@@ -4355,15 +4459,53 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
                        const unsigned int quad_no)
   :
   BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points)
+{
+  check_template_arguments(fe_no);
+  set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
+                       const DoFHandler<dim>                 &dof_handler,
+                       const unsigned int                     first_selected_component)
+  :
+  BaseClass (geometry, dof_handler, first_selected_component)
+{
+  check_template_arguments(numbers::invalid_unsigned_int);
+  set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationGeneral (const FEEvaluationGeneral &other)
+  :
+  BaseClass (other)
 {
   set_data_pointers();
+}
+
 
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::check_template_arguments(const unsigned int fe_no)
+{
 #ifdef DEBUG
   // print error message when the dimensions do not match. Propose a possible
   // fix
-  if ((dofs_per_cell != this->data->dofs_per_cell &&
-       internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value !=
-       this->data->dofs_per_cell)
+  if (fe_degree != this->data->fe_degree
       ||
       n_q_points != this->data->n_q_points)
     {
@@ -4376,33 +4518,36 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
       message += "," + Utilities::int_to_string(n_components);
       message += ",Number>(data, ";
       message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(quad_no) + ")\n";
+      message += Utilities::int_to_string(this->quad_no) + ")\n";
 
       // check whether some other vector component has the correct number of
       // points
       unsigned int proposed_dof_comp = numbers::invalid_unsigned_int,
                    proposed_quad_comp = numbers::invalid_unsigned_int;
-      if (dofs_per_cell == this->data->dofs_per_cell)
-        proposed_dof_comp = fe_no;
-      else
-        for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
-          if (this->matrix_info->get_dof_info(no).dofs_per_cell[this->active_fe_index]
-              == dofs_per_cell)
-            {
-              proposed_dof_comp = no;
-              break;
-            }
-      if (n_q_points ==
-          this->mapping_info->mapping_data_gen[quad_no].n_q_points[this->active_quad_index])
-        proposed_quad_comp = 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)
-            {
-              proposed_quad_comp = no;
-              break;
-            }
+      if (fe_no != numbers::invalid_unsigned_int)
+        {
+          if (fe_degree == this->data->fe_degree)
+            proposed_dof_comp = fe_no;
+          else
+            for (unsigned int no=0; no<this->matrix_info->n_components(); ++no)
+              if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree
+                  == fe_degree)
+                {
+                  proposed_dof_comp = no;
+                  break;
+                }
+          if (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)
+                {
+                  proposed_quad_comp = no;
+                  break;
+                }
+        }
       if (proposed_dof_comp  != numbers::invalid_unsigned_int &&
           proposed_quad_comp != numbers::invalid_unsigned_int)
         {
@@ -4422,7 +4567,7 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
             correct_pos = " ^ ";
           else
             correct_pos = "   ";
-          if (proposed_quad_comp != quad_no)
+          if (proposed_quad_comp != this->quad_no)
             correct_pos += " ^\n";
           else
             correct_pos += "  \n";
@@ -4430,18 +4575,17 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
         }
       // ok, did not find the numbers specified by the template arguments in
       // the given list. Suggest correct template arguments
-      const unsigned int proposed_fe_degree = static_cast<unsigned int>(std::pow(1.001*this->data->dofs_per_cell,1./dim))-1;
       const unsigned int proposed_n_q_points_1d = static_cast<unsigned int>(std::pow(1.001*this->data->n_q_points,1./dim));
       message += "Wrong template arguments:\n";
       message += "    Did you mean FEEvaluation<dim,";
-      message += Utilities::int_to_string(proposed_fe_degree) + ",";
+      message += Utilities::int_to_string(this->data->fe_degree) + ",";
       message += Utilities::int_to_string(proposed_n_q_points_1d);
       message += "," + Utilities::int_to_string(n_components);
       message += ",Number>(data, ";
       message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(quad_no) + ")?\n";
+      message += Utilities::int_to_string(this->quad_no) + ")?\n";
       std::string correct_pos;
-      if (proposed_fe_degree != fe_degree)
+      if (this->data->fe_degree != fe_degree)
         correct_pos = " ^";
       else
         correct_pos = "  ";
@@ -4451,55 +4595,23 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
         correct_pos += "  \n";
       message += "                                 " + correct_pos;
 
-      Assert (dofs_per_cell == this->data->dofs_per_cell &&
+      Assert (fe_degree == this->data->fe_degree &&
               n_q_points == this->data->n_q_points,
               ExcMessage(message));
     }
-  AssertDimension (n_q_points,
-                   this->mapping_info->mapping_data_gen[this->quad_no].
-                   n_q_points[this->active_quad_index]);
-  AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
-                   this->dof_info->dofs_per_cell[this->active_fe_index]);
+  if (fe_no != numbers::invalid_unsigned_int)
+    {
+      AssertDimension (n_q_points,
+                       this->mapping_info->mapping_data_gen[this->quad_no].
+                       n_q_points[this->active_quad_index]);
+      AssertDimension (this->data->dofs_per_cell * this->n_fe_components,
+                       this->dof_info->dofs_per_cell[this->active_fe_index]);
+    }
 #endif
 }
 
 
 
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
-                       const DoFHandler<dim>               &dof_handler,
-                       const unsigned int                   first_selected_component)
-  :
-  BaseClass (geometry, dof_handler, first_selected_component)
-{
-  set_data_pointers();
-
-  Assert ((dofs_per_cell == this->data->dofs_per_cell ||
-           internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value ==
-           this->data->dofs_per_cell
-           )
-          &&
-          n_q_points == this->data->n_q_points,
-          ExcMessage("Underlying element and template arguments do not match"));
-}
-
-
-template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
-          typename Number>
-inline
-FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationGeneral (const FEEvaluationGeneral &other)
-  :
-  BaseClass (other)
-{
-  set_data_pointers();
-}
-
-
-
 template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
@@ -4507,23 +4619,59 @@ void
 FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::set_data_pointers()
 {
+  AssertIndexRange(this->data->dofs_per_cell, dofs_per_cell+2);
+  const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell;
+
   // set the pointers to the correct position in the data array
   for (unsigned int c=0; c<n_components_; ++c)
     {
-      this->values_dofs[c] = &my_data_array[c*dofs_per_cell];
-      this->values_quad[c] = &my_data_array[n_components*dofs_per_cell+c*n_q_points];
+      this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell];
+      this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points];
       for (unsigned int d=0; d<dim; ++d)
-        this->gradients_quad[c][d] = &my_data_array[n_components*(dofs_per_cell+n_q_points)
+        this->gradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+
+                                                                  n_q_points)
                                                     +
                                                     (c*dim+d)*n_q_points];
       for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
-        this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+dofs_per_cell)
+        this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+
+                                                                 desired_dofs_per_cell)
                                                    +
                                                    (c*(dim*dim+dim)+d)*n_q_points];
     }
 }
 
 
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+void
+FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::check_dofs_per_cell (const unsigned int given_dofs_per_cell) const
+{
+#ifdef DEBUG
+  if (given_dofs_per_cell != this->data->dofs_per_cell)
+    {
+      std::ostringstream str;
+      str << "Dofs per cell in FEEvaluation* class does not match the element ("
+          << given_dofs_per_cell << " != " << this->data->dofs_per_cell << "). "
+          << "Suggestion: ";
+      if (this->data->dofs_per_cell == dofs_per_cell)
+        str << "FEEvaluation/FEEvaluationGeneral";
+      else if (this->data->dofs_per_cell ==
+               internal::MatrixFreeFunctions::DGP_dofs_per_cell<dim,fe_degree>::value)
+        str << "FEEvaluationDGP";
+      else if (this->data->dofs_per_cell == dofs_per_cell+1)
+        str << "FEEvaluationQ_DG0";
+      else
+        str << "No matching suggestion found";
+
+      Assert(given_dofs_per_cell == this->data->dofs_per_cell,
+             ExcMessage(str.str().c_str()));
+    }
+#endif         
+}
+
+
 
 namespace internal
 {
@@ -4656,6 +4804,19 @@ namespace internal
                 out += nn;
                 break;
               case 1:
+                ++in;
+                ++out;
+                // faces 2 and 3 in 3D use local coordinate system zx, which
+                // is the other way around compared to the tensor
+                // product. Need to take that into account.
+                if (dim == 3)
+                  {
+                    if (dof_to_quad)
+                      out += fe_degree;
+                    else
+                      in += fe_degree;
+                  }
+                break;
               case 2:
                 ++in;
                 ++out;
@@ -4664,10 +4825,15 @@ namespace internal
                 Assert (false, ExcNotImplemented());
               }
           }
-        if (face_direction == 1)
+        if (face_direction == 1 && dim == 3)
           {
             in += mm*(mm-1);
             out += nn*(nn-1);
+            // adjust for local coordinate system zx
+            if (dof_to_quad)
+              out -= (fe_degree+1)*(fe_degree+1)-1;
+            else
+              in -= (fe_degree+1)*(fe_degree+1)-1;
           }
       }
   }
@@ -5954,6 +6120,7 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
             const bool evaluate_grad,
             const bool evaluate_lapl)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->dof_values_initialized == true,
           internal::ExcAccessToUninitializedField());
   internal::do_evaluate (*this, this->values_dofs, this->values_quad,
@@ -5981,6 +6148,7 @@ FEEvaluationGeneral<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::integrate (const bool integrate_val,
              const bool integrate_grad)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   if (integrate_val == true)
     Assert (this->values_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
@@ -6244,6 +6412,7 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
             const bool evaluate_grad,
             const bool evaluate_lapl)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->dof_values_initialized == true,
           internal::ExcAccessToUninitializedField());
   internal::do_evaluate (*this, this->values_dofs, this->values_quad,
@@ -6270,6 +6439,7 @@ void
 FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::integrate (bool integrate_val,bool integrate_grad)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   if (integrate_val == true)
     Assert (this->values_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
@@ -6433,6 +6603,7 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
             const bool evaluate_grad,
             const bool evaluate_lapl)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->cell != numbers::invalid_unsigned_int,
           ExcNotInitialized());
   Assert (this->dof_values_initialized == true,
@@ -6539,6 +6710,7 @@ void
 FEEvaluationGL<dim,fe_degree,n_components_,Number>
 ::integrate (const bool integrate_val, const bool integrate_grad)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->cell != numbers::invalid_unsigned_int,
           ExcNotInitialized());
   if (integrate_val == true)
@@ -6634,9 +6806,8 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
   :
   BaseClass (data_in, fe_no, quad_no)
 {
-  // reset values_dofs pointer as it has wider gaps in the base class
-  for (unsigned int c=0; c<n_components_; ++c)
-    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+                  n_components * dofs_per_cell);
 }
 
 
@@ -6651,8 +6822,8 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
   :
   BaseClass (geometry, dof_handler, first_selected_component)
 {
-  for (unsigned int c=0; c<n_components_; ++c)
-    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
+  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+                  n_components * dofs_per_cell);
 }
 
 
@@ -6665,8 +6836,6 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
   :
   BaseClass (other)
 {
-  for (unsigned int c=0; c<n_components_; ++c)
-    this->values_dofs[c] = &this->my_data_array[c*dofs_per_cell];
 }
 
 
@@ -6680,6 +6849,7 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
             const bool evaluate_grad,
             const bool evaluate_lapl)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   Assert (this->dof_values_initialized == true,
           internal::ExcAccessToUninitializedField());
 
@@ -6733,6 +6903,7 @@ void
 FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
 ::integrate (bool integrate_val,bool integrate_grad)
 {
+  this->check_dofs_per_cell(dofs_per_cell);
   if (integrate_val == true)
     Assert (this->values_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
@@ -6773,6 +6944,122 @@ FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
 
 
 
+/*------------------------- FEEvaluationQ_DG0 -------------------------------*/
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const MatrixFree<dim,Number> &data_in,
+                     const unsigned int fe_no,
+                     const unsigned int quad_no)
+  :
+  BaseClass (data_in, fe_no, quad_no)
+{
+  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+                  n_components * dofs_per_cell);
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const MappingFEEvaluation<dim,Number> &geometry,
+                     const DoFHandler<dim>                  &dof_handler,
+                     const unsigned int                      first_selected_component)
+  :
+  BaseClass (geometry, dof_handler, first_selected_component)
+{
+  AssertDimension(static_cast<unsigned int>(this->values_quad[0]-this->values_dofs[0]),
+                  n_components * dofs_per_cell);
+}
+
+
+
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other)
+  :
+  BaseClass (other)
+{
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+void
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::evaluate (const bool evaluate_val,
+            const bool evaluate_grad,
+            const bool evaluate_lapl)
+{
+  this->check_dofs_per_cell(dofs_per_cell);
+  Assert (this->dof_values_initialized == true,
+          internal::ExcAccessToUninitializedField());
+  internal::do_evaluate (*this, this->values_dofs, this->values_quad,
+                         this->gradients_quad, this->hessians_quad,
+                         evaluate_val, evaluate_grad, evaluate_lapl,
+                         internal::int2type<dim>());
+  if (evaluate_val)
+    for (unsigned int c=0; c<n_components; ++c)
+      for (unsigned int q=0; q<n_q_points; ++q)
+        this->values_quad[c][q] += this->values_dofs[c][dofs_per_cell-1];
+
+#ifdef DEBUG
+  if (evaluate_val == true)
+    this->values_quad_initialized = true;
+  if (evaluate_grad == true)
+    this->gradients_quad_initialized = true;
+  if (evaluate_lapl == true)
+    this->hessians_quad_initialized  = true;
+#endif
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+void
+FEEvaluationQ_DG0<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::integrate (bool integrate_val,bool integrate_grad)
+{
+  this->check_dofs_per_cell(dofs_per_cell);
+  if (integrate_val == true)
+    Assert (this->values_quad_submitted == true,
+            internal::ExcAccessToUninitializedField());
+  if (integrate_grad == true)
+    Assert (this->gradients_quad_submitted == true,
+            internal::ExcAccessToUninitializedField());
+
+  if (integrate_val)
+    for (unsigned int c=0; c<n_components; ++c)
+      {
+        this->values_dofs[c][dofs_per_cell-1] = this->values_quad[c][0];
+        for (unsigned int q=1; q<n_q_points; ++q)
+          this->values_dofs[c][dofs_per_cell-1] += this->values_quad[c][q];
+      }
+  else
+    for (unsigned int c=0; c<n_components; ++c)
+      this->values_dofs[c][dofs_per_cell-1] = VectorizedArray<Number>();
+
+  internal::do_integrate (*this, this->values_dofs, this->values_quad,
+                          this->gradients_quad, integrate_val, integrate_grad,
+                          internal::int2type<dim>());
+
+#ifdef DEBUG
+  this->dof_values_initialized = true;
+#endif
+}
+
+
+
 #endif  // ifndef DOXYGEN
 
 
index 3164b4679ce354c25291a109945f2eda59aa610b..f9b51e21f6c6f1787e3d337d5e2ed6e972fc8605 100644 (file)
@@ -72,55 +72,40 @@ namespace internal
                    const unsigned int base_element = 0);
 
       /**
-       * Returns the memory consumption of this
-       * class in bytes.
+       * Returns the memory consumption of this class in bytes.
        */
       std::size_t memory_consumption () const;
 
       /**
-       * Stores the shape values of the 1D finite
-       * element evaluated on all 1D quadrature
-       * points in vectorized format, i.e., as an
-       * array of
-       * VectorizedArray<dim>::n_array_elements
-       * equal elements. The length of this array is
-       * <tt>n_dofs_1d * n_q_points_1d</tt> and
-       * quadrature points are the index running
-       * fastest.
+       * Stores the shape values of the 1D finite element evaluated on all 1D
+       * quadrature points in vectorized format, i.e., as an array of
+       * VectorizedArray<dim>::n_array_elements equal elements. The length of
+       * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+       * points are the index running fastest.
        */
       AlignedVector<VectorizedArray<Number> > shape_values;
 
       /**
-       * Stores the shape gradients of the 1D finite
-       * element evaluated on all 1D quadrature
-       * points in vectorized format, i.e., as an
-       * array of
-       * VectorizedArray<dim>::n_array_elements
-       * equal elements. The length of this array is
-       * <tt>n_dofs_1d * n_q_points_1d</tt> and
-       * quadrature points are the index running
-       * fastest.
+       * Stores the shape gradients of the 1D finite element evaluated on all
+       * 1D quadrature points in vectorized format, i.e., as an array of
+       * VectorizedArray<dim>::n_array_elements equal elements. The length of
+       * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+       * points are the index running fastest.
        */
       AlignedVector<VectorizedArray<Number> > shape_gradients;
 
       /**
-       * Stores the shape Hessians of the 1D finite
-       * element evaluated on all 1D quadrature
-       * points in vectorized format, i.e., as an
-       * array of
-       * VectorizedArray<dim>::n_array_elements
-       * equal elements. The length of this array is
-       * <tt>n_dofs_1d * n_q_points_1d</tt> and
-       * quadrature points are the index running
-       * fastest.
+       * Stores the shape Hessians of the 1D finite element evaluated on all
+       * 1D quadrature points in vectorized format, i.e., as an array of
+       * VectorizedArray<dim>::n_array_elements equal elements. The length of
+       * this array is <tt>n_dofs_1d * n_q_points_1d</tt> and quadrature
+       * points are the index running fastest.
        */
       AlignedVector<VectorizedArray<Number> > shape_hessians;
 
       /**
-       * Stores the indices from cell DoFs to face
-       * DoFs. The rows go through the
-       * <tt>2*dim</tt> faces, and the columns the
-       * DoFs on the faces.
+       * Stores the indices from cell DoFs to face DoFs. The rows go through
+       * the <tt>2*dim</tt> faces, and the columns the DoFs on the faces.
        */
       dealii::Table<2,unsigned int>  face_indices;
 
@@ -143,15 +128,14 @@ namespace internal
       std::vector<Number>    subface_value[2];
 
       /**
-       * Non-vectorized version of shape
-       * values. Needed when evaluating face info.
+       * Non-vectorized version of shape values. Needed when evaluating face
+       * info.
        */
       std::vector<Number>    shape_values_number;
 
       /**
-       * Non-vectorized version of shape
-       * gradients. Needed when evaluating face
-       * info.
+       * Non-vectorized version of shape gradients. Needed when evaluating
+       * face info.
        */
       std::vector<Number>    shape_gradient_number;
 
@@ -165,26 +149,28 @@ namespace internal
       std::vector<unsigned int> lexicographic_numbering;
 
       /**
-       * Stores the number of quadrature points in
-       * @p dim dimensions for a cell.
+       * Stores the degree of the element.
+       */
+      unsigned int fe_degree;
+
+      /**
+       * Stores the number of quadrature points in @p dim dimensions for a
+       * cell.
        */
       unsigned int n_q_points;
 
       /**
-       * Stores the number of DoFs per cell in @p
-       * dim dimensions.
+       * Stores the number of DoFs per cell in @p dim dimensions.
        */
       unsigned int dofs_per_cell;
 
       /**
-       * Stores the number of quadrature points per
-       * face in @p dim dimensions.
+       * Stores the number of quadrature points per face in @p dim dimensions.
        */
       unsigned int n_q_points_face;
 
       /**
-       * Stores the number of DoFs per face in @p
-       * dim dimensions.
+       * Stores the number of DoFs per face in @p dim dimensions.
        */
       unsigned int dofs_per_face;
     };
@@ -200,6 +186,7 @@ namespace internal
                                   const FiniteElement<dim> &fe_in,
                                   const unsigned int base_element_number)
       :
+      fe_degree (0),
       n_q_points (0),
       dofs_per_cell (0)
     {
index eb8e2bda753e9f9ade4acea0970635a58eacb904..826649644fe8b780d69a3900e2ef6334041c0e0f 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/polynomials_piecewise.h>
 #include <deal.II/fe/fe_poly.h>
 #include <deal.II/fe/fe_dgp.h>
+#include <deal.II/fe/fe_q_dg0.h>
 
 #include <deal.II/matrix_free/shape_info.h>
 
@@ -59,8 +60,9 @@ namespace internal
       Assert (fe->n_components() == 1,
               ExcMessage("FEEvaluation only works for scalar finite elements."));
 
+      fe_degree = fe->degree;
 
-      const unsigned int n_dofs_1d = fe->degree+1,
+      const unsigned int n_dofs_1d = fe_degree+1,
                          n_q_points_1d = quad.size();
 
       // renumber (this is necessary for FE_Q, for example, since there the
@@ -82,6 +84,8 @@ namespace internal
 
         const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(fe);
 
+        const FE_Q_DG0<dim> *fe_q_dg0 = dynamic_cast<const FE_Q_DG0<dim>*>(fe);
+
         if (fe_poly != 0)
           scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
         else if (fe_poly_piece != 0)
@@ -89,9 +93,11 @@ namespace internal
         else if (fe_dgp != 0)
           {
             scalar_lexicographic.resize(fe_dgp->dofs_per_cell);
-          for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
-            scalar_lexicographic[i] = i;
+            for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
+              scalar_lexicographic[i] = i;
           }
+        else if (fe_q_dg0 != 0)
+          scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse();
         else
           Assert(false, ExcNotImplemented());
 
diff --git a/tests/matrix_free/matrix_vector_17.cc b/tests/matrix_free/matrix_vector_17.cc
new file mode 100644 (file)
index 0000000..4ab84fd
--- /dev/null
@@ -0,0 +1,248 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 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 test is similar to matrix_vector_06, but implements the operations for
+// FE_QDG0 instead of FE_Q (where there is an additional discontinuous degree
+// of freedom and different routines need to be used). The data is still not
+// very useful because the matrix does not include face terms actually present
+// in an approximation of the Laplacian. It only contains cell terms.
+
+#include "../tests.h"
+#include <deal.II/base/function.h>
+#include <deal.II/fe/fe_q_dg0.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.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/lac/constraint_matrix.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/compressed_simple_sparsity_pattern.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include <deal.II/lac/vector.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename VECTOR>
+void
+helmholtz_operator_qdg0 (const MatrixFree<dim,typename VECTOR::value_type>  &data,
+                         VECTOR       &dst,
+                         const VECTOR &src,
+                         const std::pair<unsigned int,unsigned int> &cell_range)
+{
+  typedef typename VECTOR::value_type Number;
+  FEEvaluationQ_DG0<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.evaluate (true, true, false);
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          fe_eval.submit_value (Number(10)*fe_eval.get_value(q),q);
+          fe_eval.submit_gradient (fe_eval.get_gradient(q),q);
+        }
+      fe_eval.integrate (true,true);
+      fe_eval.distribute_local_to_global (dst);
+    }
+}
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+  typedef VectorizedArray<Number> vector_t;
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    const std_cxx1x::function<void(const MatrixFree<dim,typename VECTOR::value_type> &,
+                                   VECTOR &,
+                                   const VECTOR &,
+                                   const std::pair<unsigned int,unsigned int> &)>
+      wrap = helmholtz_operator_qdg0<dim,fe_degree,VECTOR>;
+    data.cell_loop (wrap, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+              const ConstraintMatrix &constraints)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MatrixFree<dim,number> mf_data;
+  {
+    const QGauss<1> quad (fe_degree+1);
+    typename MatrixFree<dim,number>::AdditionalData data;
+    data.tasks_parallel_scheme =
+      MatrixFree<dim,number>::AdditionalData::partition_color;
+    data.tasks_block_size = 3;
+
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), out (dof.n_dofs());
+  Vector<number> in_dist (dof.n_dofs());
+  Vector<number> out_dist (in_dist);
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      if (constraints.is_constrained(i))
+        continue;
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+      in_dist(i) = entry;
+    }
+
+  mf.vmult (out_dist, in_dist);
+
+
+  // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  SparsityPattern sparsity;
+  {
+    CompressedSimpleSparsityPattern csp(dof.n_dofs(), dof.n_dofs());
+    DoFTools::make_sparsity_pattern (dof, csp, constraints, true);
+    sparsity.copy_from(csp);
+  }
+  SparseMatrix<double> sparse_matrix (sparsity);
+  {
+    QGauss<dim>  quadrature_formula(fe_degree+1);
+
+    FEValues<dim> fe_values (dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    const unsigned int   dofs_per_cell = dof.get_fe().dofs_per_cell;
+    const unsigned int   n_q_points    = quadrature_formula.size();
+
+    FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+    std::vector<types::global_dof_index> local_dof_indices (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        fe_values.reinit (cell);
+
+        for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            {
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+                                      fe_values.shape_grad(j,q_point)
+                                      +
+                                      10. *
+                                      fe_values.shape_value(i,q_point) *
+                                      fe_values.shape_value(j,q_point)) *
+                                     fe_values.JxW(q_point));
+            }
+
+        cell->get_dof_indices(local_dof_indices);
+        constraints.distribute_local_to_global (cell_matrix,
+                                                local_dof_indices,
+                                                sparse_matrix);
+      }
+  }
+
+  sparse_matrix.vmult (out, in);
+  out -= out_dist;
+  const double diff_norm = out.linfty_norm() / out_dist.linfty_norm();
+
+  deallog << "Norm of difference: " << diff_norm << std::endl << std::endl;
+}
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  static const HyperBallBoundary<dim> boundary;
+  tria.set_boundary (0, boundary);
+  if (dim < 3 || fe_degree < 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+  typename Triangulation<dim>::active_cell_iterator
+  cell = tria.begin_active (),
+  endc = tria.end();
+  for (; cell!=endc; ++cell)
+    if (cell->center().norm()<1e-8)
+      cell->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_Q_DG0<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  ConstraintMatrix constraints;
+
+  do_test<dim, fe_degree, double> (dof, constraints);
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  deallog << std::setprecision (3);
+
+  {
+    deallog.threshold_double(5.e-11);
+    deallog.push("2d");
+    test<2,1>();
+    test<2,2>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
+
+
diff --git a/tests/matrix_free/matrix_vector_17.output b/tests/matrix_free/matrix_vector_17.output
new file mode 100644 (file)
index 0000000..2bbf553
--- /dev/null
@@ -0,0 +1,13 @@
+
+DEAL:2d::Testing FE_Q_DG0<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FE_Q_DG0<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FE_Q_DG0<3>(1)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_Q_DG0<3>(2)
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0.cc b/tests/matrix_free/matrix_vector_stokes_qdg0.cc
new file mode 100644 (file)
index 0000000..c746da6
--- /dev/null
@@ -0,0 +1,333 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2013 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
+// matrix-vector products by comparing with the result of deal.II sparse
+// matrix. No hanging nodes and no other
+// constraints for a vector-valued problem (stokes equations).
+
+#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/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_q_dg0.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "create_mesh.h"
+
+#include <fstream>
+#include <iostream>
+#include <complex>
+
+
+
+template <int dim, int degree_p, 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,degree_p+1,degree_p+2,dim,Number> velocity (data, 0);
+    FEEvaluationQ_DG0<dim,degree_p  ,degree_p+2,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, 0);
+        velocity.evaluate (false,true,false);
+        pressure.reinit (cell);
+        pressure.read_dof_values (src, dim);
+        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 = -velocity.get_divergence(q);
+            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, 0);
+        pressure.integrate (true,false);
+        pressure.distribute_local_to_global (dst, dim);
+      }
+  }
+
+
+  void vmult (VectorType &dst,
+              const VectorType &src) const
+  {
+    AssertDimension (dst.size(), dim+1);
+    for (unsigned int d=0; d<dim+1; ++d)
+      dst[d] = 0;
+    data.cell_loop (&MatrixFreeTest<dim,degree_p,VectorType>::local_apply,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree>
+void test ()
+{
+  Triangulation<dim>   triangulation;
+  create_mesh (triangulation);
+  if (fe_degree == 1)
+    triangulation.refine_global (4-dim);
+  else
+    triangulation.refine_global (3-dim);
+
+  FE_Q<dim>            fe_u (fe_degree+1);
+  FE_Q_DG0<dim>        fe_p (fe_degree);
+  FESystem<dim>        fe (fe_u, 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;
+
+  BlockSparsityPattern      sparsity_pattern;
+  BlockSparseMatrix<double> system_matrix;
+
+  BlockVector<double> solution;
+  BlockVector<double> system_rhs;
+  std::vector<Vector<double> > vec1, vec2;
+
+  dof_handler.distribute_dofs (fe);
+  dof_handler_u.distribute_dofs (fe_u);
+  dof_handler_p.distribute_dofs (fe_p);
+  DoFRenumbering::component_wise (dof_handler);
+
+  constraints.close ();
+
+  std::vector<types::global_dof_index> dofs_per_block (dim+1);
+  DoFTools::count_dofs_per_component (dof_handler, dofs_per_block);
+
+  //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;
+
+  {
+    BlockCompressedSimpleSparsityPattern csp (dim+1,dim+1);
+
+    for (unsigned int d=0; d<dim+1; ++d)
+      for (unsigned int e=0; e<dim+1; ++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);
+
+  solution.reinit (dim+1);
+  for (unsigned int i=0; i<dim+1; ++i)
+    solution.block(i).reinit (dofs_per_block[i]);
+  solution.collect_sizes ();
+
+  system_rhs.reinit (solution);
+
+  vec1.resize (dim+1);
+  vec2.resize (dim+1);
+  vec1[0].reinit (dofs_per_block[0]);
+  vec2[0].reinit (vec1[0]);
+  for (unsigned int i=1; i<dim; ++i)
+    {
+      vec1[i].reinit (vec1[0]);
+      vec2[i].reinit (vec1[0]);
+    }
+  vec1[dim].reinit (dofs_per_block[dim]);
+  vec2[dim].reinit (vec1[dim]);
+
+  // this is from step-22
+  {
+    QGauss<dim>   quadrature_formula(fe_degree+2);
+
+    FEValues<dim> fe_values (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);
+      }
+  }
+
+  // first system_rhs with random numbers
+  for (unsigned int i=0; i<dim+1; ++i)
+    for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+      {
+        const double val = -1. + 2.*(double)Testing::rand()/double(RAND_MAX);
+        system_rhs.block(i)(j) = val;
+        vec1[i](j) = val;
+      }
+
+  // setup matrix-free structure
+  {
+    std::vector<const DoFHandler<dim>*> dofs;
+    dofs.push_back(&dof_handler_u);
+    dofs.push_back(&dof_handler_p);
+    ConstraintMatrix dummy_constraints;
+    dummy_constraints.close();
+    std::vector<const ConstraintMatrix *> constraints;
+    constraints.push_back (&dummy_constraints);
+    constraints.push_back (&dummy_constraints);
+    QGauss<1> quad(fe_degree+2);
+    mf_data.reinit (dofs, constraints, quad,
+                    typename MatrixFree<dim>::AdditionalData
+                    (MPI_COMM_WORLD,
+                     MatrixFree<dim>::AdditionalData::none));
+  }
+
+  system_matrix.vmult (solution, system_rhs);
+
+  typedef std::vector<Vector<double> > VectorType;
+  MatrixFreeTest<dim,fe_degree,VectorType> mf (mf_data);
+  mf.vmult (vec2, vec1);
+
+  // Verification
+  double error = 0.;
+  for (unsigned int i=0; i<dim+1; ++i)
+    for (unsigned int j=0; j<system_rhs.block(i).size(); ++j)
+      error += std::fabs (solution.block(i)(j)-vec2[i](j));
+  double relative = solution.block(0).l1_norm();
+  deallog << "  Verification fe degree " << fe_degree  <<  ": "
+          << error/relative << std::endl << std::endl;
+}
+
+
+
+int main ()
+{
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  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>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
+
diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0.output b/tests/matrix_free/matrix_vector_stokes_qdg0.output
new file mode 100644 (file)
index 0000000..db6c7de
--- /dev/null
@@ -0,0 +1,14 @@
+
+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::
+DEAL:3d::  Verification fe degree 2: 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.