]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Change name of a class to better represent what it is doing, add inverse mass operato...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Feb 2014 22:09:03 +0000 (22:09 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 24 Feb 2014 22:09:03 +0000 (22:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@32546 0785d39b-7218-0410-832d-ea1e28bc413d

17 files changed:
deal.II/include/deal.II/matrix_free/fe_evaluation.h
deal.II/include/deal.II/matrix_free/mapping_fe_evaluation.h [moved from deal.II/include/deal.II/matrix_free/evaluated_geometry.h with 65% similarity]
deal.II/include/deal.II/matrix_free/matrix_free.templates.h
deal.II/include/deal.II/matrix_free/operators.h [new file with mode: 0644]
deal.II/include/deal.II/matrix_free/shape_info.h
deal.II/include/deal.II/matrix_free/shape_info.templates.h
tests/matrix_free/inverse_mass_01.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_01.output [new file with mode: 0644]
tests/matrix_free/inverse_mass_02.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_02.output [new file with mode: 0644]
tests/matrix_free/inverse_mass_03.cc [new file with mode: 0644]
tests/matrix_free/inverse_mass_03.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_14.cc
tests/matrix_free/matrix_vector_15.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_15.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_16.cc [new file with mode: 0644]
tests/matrix_free/matrix_vector_16.output [new file with mode: 0644]

index f49488d2c7a4fd05d2f2c7fa32f8dabc86ea276a..bd755ee0acfce69af4d3b3f785fe113fb86fdc1a 100644 (file)
@@ -27,7 +27,7 @@
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/matrix_free/matrix_free.h>
 #include <deal.II/matrix_free/shape_info.h>
-#include <deal.II/matrix_free/evaluated_geometry.h>
+#include <deal.II/matrix_free/mapping_fe_evaluation.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -115,6 +115,18 @@ public:
    */
   internal::MatrixFreeFunctions::CellType get_cell_type() const;
 
+  /**
+   * Returns a reference to the ShapeInfo object currently in use.
+   */
+  const internal::MatrixFreeFunctions::ShapeInfo<Number> &
+  get_shape_info() const;
+
+  /**
+   * Fills the JxW values currently used.
+   */
+  void
+  fill_JxW_values(AlignedVector<VectorizedArray<Number> > &JxW_values) const;
+
   //@}
 
   /**
@@ -136,6 +148,14 @@ public:
    * ConstraintMatrix::distribute. When accessing vector entries during the
    * solution of linear systems, the temporary solution should always have
    * homogeneous constraints and this method is the correct one.
+   *
+   * If the class was constructed through a MappingFEEvaluation object, only
+   * one single cell is used by this class and this function extracts the
+   * values of the underlying components on this cell. This call is slower
+   * than the ones done through a MatrixFree object and lead to a structure
+   * that does not effectively use vectorization in the evaluate routines
+   * based on these values (instead, VectorizedArray<Number>::n_array_elements
+   * same copies are worked on).
    */
   template <typename VectorType>
   void read_dof_values (const VectorType &src);
@@ -171,6 +191,15 @@ public:
    * function is also necessary when inhomogeneous constraints are to be used,
    * as MatrixFree can only handle homogeneous constraints. Note that if
    * vectorization is enabled, the DoF values for several cells are set.
+   *
+   * If the class was constructed through a MappingFEEvaluation object, only
+   * one single cell is used by this class and this function extracts the
+   * values of the underlying components on this cell. This call is slower
+   * than the ones done through a MatrixFree object and lead to a structure
+   * that does not effectively use vectorization in the evaluate routines
+   * based on these values (instead, VectorizedArray<Number>::n_array_elements
+   * same copies are worked on). In that case, no constraints can be
+   * processed as these are not available here.
    */
   template <typename VectorType>
   void read_dof_values_plain (const VectorType &src);
@@ -199,8 +228,17 @@ public:
    * Takes the values stored internally on dof values of the current cell and
    * sums them into the vector @p dst. The function also applies constraints
    * during the write operation. The functionality is hence similar to the
-   * function ConstraintMatrix::distribute_local_to_global.  Note that if
-   * vectorization is enabled, the DoF values for several cells are used.
+   * function ConstraintMatrix::distribute_local_to_global. If vectorization
+   * is enabled, the DoF values for several cells are used.
+   *
+   * If the class was constructed through a MappingFEEvaluation object, only
+   * one single cell is used by this class and this function extracts the
+   * values of the underlying components on this cell. This call is slower
+   * than the ones done through a MatrixFree object and lead to a structure
+   * that does not effectively use vectorization in the evaluate routines
+   * based on these values (instead, VectorizedArray<Number>::n_array_elements
+   * same copies are worked on). In that case, no constraints can be
+   * processed as these are not available here.
    */
   template<typename VectorType>
   void distribute_local_to_global (VectorType &dst) const;
@@ -211,8 +249,8 @@ public:
    * and sums them into the collection of vectors vector @p dst, starting at
    * index @p first_index. The function also applies constraints during the
    * write operation. The functionality is hence similar to the function
-   * ConstraintMatrix::distribute_local_to_global.  Note that if vectorization
-   * is enabled, the DoF values for several cells are used.
+   * ConstraintMatrix::distribute_local_to_global. If vectorization is
+   * enabled, the DoF values for several cells are used.
    */
   template<typename VectorType>
   void distribute_local_to_global (std::vector<VectorType> &dst,
@@ -232,6 +270,15 @@ public:
    * during the write operation. The functionality is hence similar to the
    * function ConstraintMatrix::distribute_local_to_global.  Note that if
    * vectorization is enabled, the DoF values for several cells are used.
+   *
+   * If the class was constructed through a MappingFEEvaluation object, only
+   * one single cell is used by this class and this function extracts the
+   * values of the underlying components on this cell. This call is slower
+   * than the ones done through a MatrixFree object and lead to a structure
+   * that does not effectively use vectorization in the evaluate routines
+   * based on these values (instead, VectorizedArray<Number>::n_array_elements
+   * same copies are worked on). In that case, no constraints can be
+   * processed as these are not available here.
    */
   template<typename VectorType>
   void set_dof_values (VectorType &dst) const;
@@ -512,7 +559,7 @@ protected:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -526,9 +573,9 @@ protected:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationBase (const EvaluatedGeometry<dim,Number> &geometry,
-                    const DoFHandler<dim>               &dof_handler,
-                    const unsigned int                   base_element = 0);
+  FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
+                    const DoFHandler<dim>                 &dof_handler,
+                    const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -638,13 +685,18 @@ protected:
   const internal::MatrixFreeFunctions::MappingInfo<dim,Number> *mapping_info;
 
   /**
-   * Stores a reference to the unit cell data, i.e., values, gradients and
+   * In case the class is initialized from MappingFEEvaluation instead of
+   * MatrixFree, this data structure holds the evaluated shape data.
+   */
+  std_cxx1x::shared_ptr<internal::MatrixFreeFunctions::ShapeInfo<Number> > stored_shape_info;
+
+  /**
+   * Stores a pointer to the unit cell shape data, i.e., values, gradients and
    * Hessians in 1D at the quadrature points that constitute the tensor
    * product. Also contained in matrix_info, but it simplifies code if we
-   * store a reference to it. If the object is initialized without MatrixFree
-   * object, the constructor creates this data structure.
+   * store a reference to it.
    */
-  std_cxx1x::shared_ptr<const internal::MatrixFreeFunctions::ShapeInfo<Number> > data;
+  const internal::MatrixFreeFunctions::ShapeInfo<Number> *data;
 
   /**
    * A pointer to the Cartesian Jacobian information of the present cell. Only
@@ -754,12 +806,24 @@ protected:
   /**
    * Geometry data generated by FEValues on the fly.
    */
-  SmartPointer<const EvaluatedGeometry<dim,Number> > evaluated_geometry;
+  SmartPointer<const MappingFEEvaluation<dim,Number> > mapped_geometry;
 
   /**
    * A pointer to the underlying DoFHandler.
    */
   const DoFHandler<dim> *dof_handler;
+
+  /**
+   * For a DoFHandler with more than one finite element, select at which
+   * component this data structure should start.
+   */
+  const unsigned int first_selected_component;
+
+  /**
+   * A temporary data structure necessary to read degrees of freedom when no
+   * MatrixFree object was given at initialization.
+   */
+   mutable std::vector<types::global_dof_index> local_dof_indices;
 };
 
 
@@ -799,7 +863,7 @@ protected:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -813,9 +877,9 @@ protected:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element = 0);
+  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+                      const DoFHandler<dim>                 &dof_handler,
+                      const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -946,7 +1010,7 @@ protected:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -960,9 +1024,9 @@ protected:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element = 0);
+  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+                      const DoFHandler<dim>                 &dof_handler,
+                      const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -1102,7 +1166,7 @@ protected:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -1116,9 +1180,9 @@ protected:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element = 0);
+  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
+                      const DoFHandler<dim>                 &dof_handler,
+                      const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -1193,7 +1257,7 @@ public:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -1207,9 +1271,9 @@ public:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationGeneral (const EvaluatedGeometry<dim,Number> &geometry,
-                       const DoFHandler<dim>               &dof_handler,
-                       const unsigned int                   base_element = 0);
+  FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
+                       const DoFHandler<dim>                 &dof_handler,
+                       const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -1309,7 +1373,7 @@ private:
  * possible to use vectorization for applying a vector operation for several
  * cells at once. The second form of usage is to initialize it from geometry
  * information generated by FEValues, which is stored in the class
- * EvaluatedGeometry. Here, the operations can only work on a single cell, but
+ * MappingFEEvaluation. Here, the operations can only work on a single cell, but
  * possibly be vectorized by combining several operations (e.g. when
  * performing matrix assembly).
  *
@@ -1379,7 +1443,7 @@ public:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -1393,9 +1457,9 @@ public:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluation (const EvaluatedGeometry<dim,Number> &geometry,
-                const DoFHandler<dim>               &dof_handler,
-                const unsigned int                   base_element = 0);
+  FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
+                const DoFHandler<dim>                 &dof_handler,
+                const unsigned int                     first_selected_component = 0);
 
   /**
    * Evaluates the function values, the gradients, and the Laplacians of the
@@ -1533,7 +1597,7 @@ public:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -1547,9 +1611,9 @@ public:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationGL (const EvaluatedGeometry<dim,Number> &geometry,
-                  const DoFHandler<dim>               &dof_handler,
-                  const unsigned int                   base_element = 0);
+  FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
+                  const DoFHandler<dim>                 &dof_handler,
+                  const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -1660,7 +1724,7 @@ public:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type EvaluatedGeometry
+   * 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
@@ -1674,9 +1738,9 @@ public:
    * operations), but only possibly within the element if the
    * evaluate/integrate routines are combined (e.g. for matrix assembly).
    */
-  FEEvaluationDGP (const EvaluatedGeometry<dim,Number> &geometry,
-                   const DoFHandler<dim>               &dof_handler,
-                   const unsigned int                   base_element = 0);
+  FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
+                   const DoFHandler<dim>                 &dof_handler,
+                   const unsigned int                     first_selected_component = 0);
 
   /**
    * Copy constructor
@@ -1712,33 +1776,6 @@ public:
 #ifndef DOXYGEN
 
 
-namespace internal
-{
-  namespace MatrixFreeFunctions
-  {
-    // a small class that gives control over the delete behavior of
-    // std::shared_ptr: we need to disable it when we initialize a pointer
-    // from another structure.
-    template <typename CLASS>
-    struct DummyDeleter
-    {
-      DummyDeleter (const bool do_delete = false)
-        :
-        do_delete(do_delete)
-      {}
-
-      void operator () (CLASS *pointer)
-      {
-        if (do_delete)
-          delete pointer;
-      }
-
-      const bool do_delete;
-    };
-  }
-}
-
-
 
 /*----------------------- FEEvaluationBase ----------------------------------*/
 
@@ -1763,9 +1800,7 @@ FEEvaluationBase<dim,n_components_,Number>
   mapping_info       (&data_in.get_mapping_info()),
   data               (&data_in.get_shape_info
                       (fe_no_in, quad_no_in, active_fe_index,
-                       active_quad_index),
-                      internal::MatrixFreeFunctions::DummyDeleter
-                      <const internal::MatrixFreeFunctions::ShapeInfo<Number> >(false)),
+                       active_quad_index)),
   cartesian_data     (0),
   jacobian           (0),
   J_value            (0),
@@ -1776,7 +1811,9 @@ FEEvaluationBase<dim,n_components_,Number>
   jacobian_grad_upper(0),
   cell               (numbers::invalid_unsigned_int),
   cell_type          (internal::MatrixFreeFunctions::undefined),
-  cell_data_number   (0)
+  cell_data_number   (0),
+  dof_handler        (0),
+  first_selected_component (0)
 {
   for (unsigned int c=0; c<n_components_; ++c)
     {
@@ -1812,9 +1849,9 @@ FEEvaluationBase<dim,n_components_,Number>
 template <int dim, int n_components_, typename Number>
 inline
 FEEvaluationBase<dim,n_components_,Number>
-::FEEvaluationBase (const EvaluatedGeometry<dim,Number> &geometry,
-                    const DoFHandler<dim>               &dof_handler_in,
-                    const unsigned int                   base_element)
+::FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
+                    const DoFHandler<dim>                 &dof_handler_in,
+                    const unsigned int                     first_selected_component)
   :
   quad_no            (-1),
   n_fe_components    (n_components_),
@@ -1823,7 +1860,8 @@ FEEvaluationBase<dim,n_components_,Number>
   matrix_info        (0),
   dof_info           (0),
   mapping_info       (0),
-  data               (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), base_element)),
+  stored_shape_info  (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), 0)),
+  data               (stored_shape_info.get()),
   cartesian_data     (0),
   jacobian           (geometry.get_inverse_jacobians().begin()),
   J_value            (geometry.get_JxW_values().begin()),
@@ -1834,8 +1872,9 @@ FEEvaluationBase<dim,n_components_,Number>
   cell               (0),
   cell_type          (internal::MatrixFreeFunctions::general),
   cell_data_number   (0),
-  evaluated_geometry (&geometry),
-  dof_handler        (&dof_handler_in)
+  mapped_geometry    (&geometry),
+  dof_handler        (&dof_handler_in),
+  first_selected_component (first_selected_component)
 {
   for (unsigned int c=0; c<n_components_; ++c)
     {
@@ -1846,8 +1885,8 @@ FEEvaluationBase<dim,n_components_,Number>
       for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
         hessians_quad[c][d] = 0;
     }
-  Assert(dof_handler->get_fe().element_multiplicity(base_element) == 1 ||
-         dof_handler->get_fe().element_multiplicity(base_element) >= n_components_,
+  Assert(dof_handler->get_fe().element_multiplicity(0) == 1 ||
+         dof_handler->get_fe().element_multiplicity(0)-first_selected_component >= n_components_,
          ExcMessage("The underlying element must at least contain as many "
                     "components as requested by this class"));
 }
@@ -1877,8 +1916,9 @@ FEEvaluationBase<dim,n_components_,Number>
   cell               (other.cell),
   cell_type          (other.cell_type),
   cell_data_number   (other.cell_data_number),
-  evaluated_geometry (other.evaluated_geometry),
-  dof_handler        (other.dof_handler)
+  mapped_geometry    (other.mapped_geometry),
+  dof_handler        (other.dof_handler),
+  first_selected_component (other.first_selected_component)
 {
   for (unsigned int c=0; c<n_components_; ++c)
     {
@@ -1893,15 +1933,13 @@ FEEvaluationBase<dim,n_components_,Number>
 
 
 
-
 template <int dim, int n_components_, typename Number>
 inline
 void
-FEEvaluationBase<dim,n_components_,Number>
-::reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
 {
-  Assert (evaluated_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible"));
-  if (evaluated_geometry != 0)
+  Assert (mapped_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible"));
+  if (mapped_geometry != 0)
     return;
   Assert (dof_info != 0, ExcNotInitialized());
   Assert (mapping_info != 0, ExcNotInitialized());
@@ -1988,8 +2026,7 @@ FEEvaluationBase<dim,n_components_,Number>
 template <int dim, int n_components_, typename Number>
 inline
 internal::MatrixFreeFunctions::CellType
-FEEvaluationBase<dim,n_components_,Number>
-::get_cell_type () const
+FEEvaluationBase<dim,n_components_,Number>::get_cell_type () const
 {
   Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   return cell_type;
@@ -1997,6 +2034,40 @@ FEEvaluationBase<dim,n_components_,Number>
 
 
 
+template <int dim, int n_components_, typename Number>
+inline
+const internal::MatrixFreeFunctions::ShapeInfo<Number> &
+FEEvaluationBase<dim,n_components_,Number>::get_shape_info() const
+{
+  Assert(data != 0, ExcInternalError());
+  return *data;
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::fill_JxW_values(AlignedVector<VectorizedArray<Number> > &JxW_values) const
+{
+  AssertDimension(JxW_values.size(), data->n_q_points);
+  Assert (this->J_value != 0, ExcNotImplemented());
+  if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
+      this->cell_type == internal::MatrixFreeFunctions::affine)
+    {
+      Assert (this->mapping_info != 0, ExcNotImplemented());
+      VectorizedArray<Number> J = this->J_value[0];
+      for (unsigned int q=0; q<this->data->n_q_points; ++q)
+        JxW_values[q] = J * this->quadrature_weights[q];
+    }
+  else
+    for (unsigned int q=0; q<data->n_q_points; ++q)
+      JxW_values[q] = this->J_value[q];
+}
+
+
+
 namespace internal
 {
   // write access to generic vectors that have operator ().
@@ -2086,6 +2157,14 @@ namespace internal
       res = vector_access (const_cast<const VectorType &>(vec), index);
     }
 
+    template <typename VectorType>
+    void process_dof_global (const types::global_dof_index index,
+                             VectorType         &vec,
+                             Number             &res) const
+    {
+      res = const_cast<const VectorType &>(vec)(index);
+    }
+
     void pre_constraints (const Number &,
                           Number       &res) const
     {
@@ -2125,6 +2204,14 @@ namespace internal
       vector_access (vec, index) += res;
     }
 
+    template <typename VectorType>
+    void process_dof_global (const types::global_dof_index index,
+                             VectorType         &vec,
+                             Number             &res) const
+    {
+      vec(index) += res;
+    }
+
     void pre_constraints (const Number &input,
                           Number       &res) const
     {
@@ -2163,6 +2250,14 @@ namespace internal
       vector_access (vec, index) = res;
     }
 
+    template <typename VectorType>
+    void process_dof_global (const types::global_dof_index index,
+                             VectorType         &vec,
+                             Number             &res) const
+    {
+      vec(index) = res;
+    }
+
     void pre_constraints (const Number &,
                           Number &) const
     {
@@ -2235,7 +2330,32 @@ FEEvaluationBase<dim,n_components_,Number>
   // into the local data field or write local data into the vector. Certain
   // operations are no-ops for the given use case.
 
-  Assert (matrix_info != 0, ExcNotInitialized());
+  // Case 1: No MatrixFree object given, simple case because we do not need to
+  // process constraints and need not care about vectorization
+  if (matrix_info == 0)
+    {
+      Assert (dof_handler != 0, ExcNotInitialized());
+      typename DoFHandler<dim>::cell_iterator cell (&dof_handler->get_tria(),
+                                                    mapped_geometry->get_cell()->level(),
+                                                    mapped_geometry->get_cell()->index(),
+                                                    dof_handler);
+      local_dof_indices.resize(dof_handler->get_fe().dofs_per_cell);
+      cell->get_dof_indices(local_dof_indices);
+
+      unsigned int index = first_selected_component * this->data->dofs_per_cell;
+      for (unsigned int comp = 0; comp<n_components; ++comp)
+        {
+          for (unsigned int i=0; i<this->data->dofs_per_cell; ++i, ++index)
+            {
+              operation.process_dof_global(local_dof_indices[this->data->lexicographic_numbering[index]],
+                                           *src[0], values_dofs[comp][i][0]);
+              for (unsigned int v=1; v<VectorizedArray<Number>::n_array_elements; ++v)
+                operation.process_empty(values_dofs[comp][i][v]);
+            }
+        }
+      return;
+    }
+
   Assert (dof_info != 0, ExcNotInitialized());
   Assert (matrix_info->indices_initialized() == true,
           ExcNotInitialized());
@@ -2843,9 +2963,16 @@ void
 FEEvaluationBase<dim,n_components_,Number>
 ::read_dof_values_plain (const VectorType *src[])
 {
+  // Case without MatrixFree initialization object
+  if (matrix_info == 0)
+    {
+      internal::VectorReader<Number> reader;
+      read_write_operation (reader, src);
+      return;
+    }
+
   // this is different from the other three operations because we do not use
   // constraints here, so this is a separate function.
-  Assert (matrix_info != 0, ExcNotInitialized());
   Assert (dof_info != 0, ExcNotInitialized());
   Assert (matrix_info->indices_initialized() == true,
           ExcNotInitialized());
@@ -3549,11 +3676,11 @@ FEEvaluationAccess<dim,n_components_,Number>
 template <int dim, int n_components_, typename Number>
 inline
 FEEvaluationAccess<dim,n_components_,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element)
+                      const unsigned int                   first_selected_component)
   :
-  FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, base_element)
+  FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, first_selected_component)
 {}
 
 
@@ -3590,11 +3717,11 @@ FEEvaluationAccess<dim,1,Number>
 template <int dim, typename Number>
 inline
 FEEvaluationAccess<dim,1,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element)
+                      const unsigned int                   first_selected_component)
   :
-  FEEvaluationBase <dim,1,Number> (geometry, dof_handler, base_element)
+  FEEvaluationBase <dim,1,Number> (geometry, dof_handler, first_selected_component)
 {}
 
 
@@ -3824,11 +3951,11 @@ FEEvaluationAccess<dim,dim,Number>
 template <int dim, typename Number>
 inline
 FEEvaluationAccess<dim,dim,Number>
-::FEEvaluationAccess (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
                       const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   base_element)
+                      const unsigned int                   first_selected_component)
   :
-  FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, base_element)
+  FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, first_selected_component)
 {}
 
 
@@ -4284,11 +4411,11 @@ 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 EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
                        const DoFHandler<dim>               &dof_handler,
-                       const unsigned int                   base_element)
+                       const unsigned int                   first_selected_component)
   :
-  BaseClass (geometry, dof_handler, base_element)
+  BaseClass (geometry, dof_handler, first_selected_component)
 {
   set_data_pointers();
 }
@@ -5917,11 +6044,11 @@ template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
           typename Number>
 inline
 FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluation (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
                 const DoFHandler<dim>               &dof_handler,
-                const unsigned int                   base_element)
+                const unsigned int                   first_selected_component)
   :
-  BaseClass (geometry, dof_handler, base_element)
+  BaseClass (geometry, dof_handler, first_selected_component)
 {
   compute_even_odd_factors();
 }
@@ -6213,11 +6340,11 @@ FEEvaluationGL<dim,fe_degree,n_components_,Number>
 template <int dim, int fe_degree, int n_components_, typename Number>
 inline
 FEEvaluationGL<dim,fe_degree,n_components_,Number>
-::FEEvaluationGL (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
                   const DoFHandler<dim>               &dof_handler,
-                  const unsigned int                   base_element)
+                  const unsigned int                   first_selected_component)
   :
-  BaseClass (geometry, dof_handler, base_element)
+  BaseClass (geometry, dof_handler, first_selected_component)
 {}
 
 
@@ -6448,11 +6575,11 @@ template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
           typename Number>
 inline
 FEEvaluationDGP<dim,fe_degree,n_q_points_1d,n_components_,Number>
-::FEEvaluationDGP (const EvaluatedGeometry<dim,Number> &geometry,
+::FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
                 const DoFHandler<dim>                  &dof_handler,
-                const unsigned int                      base_element)
+                const unsigned int                      first_selected_component)
   :
-  BaseClass (geometry, dof_handler, base_element)
+  BaseClass (geometry, dof_handler, first_selected_component)
 {}
 
 
similarity index 65%
rename from deal.II/include/deal.II/matrix_free/evaluated_geometry.h
rename to deal.II/include/deal.II/matrix_free/mapping_fe_evaluation.h
index 28d486ca53e10a8d67b3910f1d05002b8ee317e9..2c42f1e98fd5d4760473fadcba1356abf974617c 100644 (file)
@@ -15,8 +15,8 @@
 // ---------------------------------------------------------------------
 
 
-#ifndef __deal2__matrix_free_evaluated_geometry_h
-#define __deal2__matrix_free_evaluated_geometry_h
+#ifndef __deal2__matrix_free_mapping_fe_evaluation_h
+#define __deal2__matrix_free_mapping_fe_evaluation_h
 
 
 #include <deal.II/base/config.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/mapping_info.h>
 #include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/fe_nothing.h>
 
 
 DEAL_II_NAMESPACE_OPEN
 
 
 /**
- * The class makes FEEvaluation with the mapping information generated by
- * FEValues.
+ * This class provides evaluated mapping information using standard deal.II
+ * information in a form that FEEvaluation and friends can use for vectorized
+ * access. Since no vectorization over cells is available with the
+ * DoFHandler/Triangulation cell iterators, the interface to FEEvaluation's
+ * vectorization model is to use @p VectorizedArray::n_array_element copies of
+ * the same element. This interface is thus primarily useful for evaluating
+ * several operators on the same cell, e.g., when assembling cell matrices.
+ *
+ * @author Martin Kronbichler
  */
 template <int dim, typename Number=double>
-class EvaluatedGeometry : public Subscriptor
+class MappingFEEvaluation : public Subscriptor
 {
 public:
   /**
-   * Constructor, similar to FEValues.
+   * Constructor, similar to FEValues. Since this class only evaluates the
+   * geometry, no finite element has to be specified and the simplest element,
+   * FE_Nothing, is used internally for the underlying FEValues object.
    */
-  EvaluatedGeometry (const Mapping<dim> &mapping,
-                     const FiniteElement<dim> &fe,
+  MappingFEEvaluation (const Mapping<dim> &mapping,
                      const Quadrature<1> &quadrature,
                      const UpdateFlags update_flags);
 
   /**
    * Constructor. Instead of providing a mapping, use MappingQ1.
    */
-  EvaluatedGeometry (const FiniteElement<dim> &fe,
-                     const Quadrature<1> &quadrature,
+  MappingFEEvaluation (const Quadrature<1> &quadrature,
                      const UpdateFlags update_flags);
 
   /**
-   * Initialize with the given cell iterator.
+   * Initialize with the given cell iterator. This works for all type of
+   * iterators that can be converted into a Triangulation::cell_iterator
+   * (DoFHandler::cell_iterator or Triangulation::cell_iterator).
    */
   template <typename ITERATOR>
-  void reinit(ITERATOR &cell);
+  void reinit(ITERATOR cell);
+
+  /**
+   * Return a triangulation iterator to the current cell.
+   */
+  typename Triangulation<dim>::cell_iterator get_cell () const;
+
+  /**
+   * Return a reference to the underlying FEValues object that evaluates
+   * certain quantities (only mapping-related ones like Jacobians or mapped
+   * quadrature points are accessible, as no finite element data is actually
+   * used).
+   */
+  const FEValues<dim>& get_fe_values () const;
 
   /**
    * Return a vector of inverse transpose Jacobians. For compatibility with
@@ -100,6 +123,12 @@ public:
   get_quadrature () const;
 
 private:
+  /**
+   * Dummy finite element object necessary for initializing the FEValues
+   * object.
+   */
+  FE_Nothing<dim> fe_dummy;
+
   /**
    * An underlying FEValues object that performs the (scalar) evaluation.
    */
@@ -136,12 +165,11 @@ private:
 
 template <int dim, typename Number>
 inline
-EvaluatedGeometry<dim,Number>::EvaluatedGeometry (const Mapping<dim> &mapping,
-                                                  const FiniteElement<dim> &fe,
-                                                  const Quadrature<1> &quadrature,
-                                                  const UpdateFlags update_flags)
+MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Mapping<dim> &mapping,
+                                                      const Quadrature<1> &quadrature,
+                                                      const UpdateFlags update_flags)
   :
-  fe_values(mapping, fe, Quadrature<dim>(quadrature),
+  fe_values(mapping, fe_dummy, Quadrature<dim>(quadrature),
             internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
   quadrature_1d(quadrature),
   inverse_jacobians(fe_values.get_quadrature().size()),
@@ -157,11 +185,10 @@ EvaluatedGeometry<dim,Number>::EvaluatedGeometry (const Mapping<dim> &mapping,
 
 template <int dim, typename Number>
 inline
-EvaluatedGeometry<dim,Number>::EvaluatedGeometry (const FiniteElement<dim> &fe,
-                                                  const Quadrature<1> &quadrature,
-                                                  const UpdateFlags update_flags)
+MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Quadrature<1> &quadrature,
+                                                      const UpdateFlags update_flags)
   :
-  fe_values(fe, Quadrature<dim>(quadrature),
+  fe_values(fe_dummy, Quadrature<dim>(quadrature),
             internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
   quadrature_1d(quadrature),
   inverse_jacobians(fe_values.get_quadrature().size()),
@@ -179,9 +206,9 @@ template <int dim, typename Number>
 template <typename ITERATOR>
 inline
 void
-EvaluatedGeometry<dim,Number>::reinit(ITERATOR &cell)
+MappingFEEvaluation<dim,Number>::reinit(ITERATOR cell)
 {
-  fe_values.reinit(cell);
+  fe_values.reinit(static_cast<typename Triangulation<dim>::cell_iterator>(cell));
   for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
     {
       if (fe_values.get_update_flags() & update_inverse_jacobians)
@@ -201,10 +228,30 @@ EvaluatedGeometry<dim,Number>::reinit(ITERATOR &cell)
 
 
 
+template <int dim, typename Number>
+inline
+typename Triangulation<dim>::cell_iterator
+MappingFEEvaluation<dim,Number>::get_cell() const
+{
+  return fe_values.get_cell();
+}
+
+
+
+template <int dim, typename Number>
+inline
+const FEValues<dim>&
+MappingFEEvaluation<dim,Number>::get_fe_values() const
+{
+  return fe_values;
+}
+
+
+
 template <int dim, typename Number>
 inline
 const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_inverse_jacobians() const
+MappingFEEvaluation<dim,Number>::get_inverse_jacobians() const
 {
   return inverse_jacobians;
 }
@@ -214,7 +261,7 @@ EvaluatedGeometry<dim,Number>::get_inverse_jacobians() const
 template <int dim, typename Number>
 inline
 const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_normal_vectors() const
+MappingFEEvaluation<dim,Number>::get_normal_vectors() const
 {
   return normal_vectors;
 }
@@ -224,7 +271,7 @@ EvaluatedGeometry<dim,Number>::get_normal_vectors() const
 template <int dim, typename Number>
 inline
 const AlignedVector<Point<dim,VectorizedArray<Number> > >&
-EvaluatedGeometry<dim,Number>::get_quadrature_points() const
+MappingFEEvaluation<dim,Number>::get_quadrature_points() const
 {
   return quadrature_points;
 }
@@ -234,7 +281,7 @@ EvaluatedGeometry<dim,Number>::get_quadrature_points() const
 template <int dim, typename Number>
 inline
 const AlignedVector<VectorizedArray<Number> >&
-EvaluatedGeometry<dim,Number>::get_JxW_values() const
+MappingFEEvaluation<dim,Number>::get_JxW_values() const
 {
   return jxw_values;
 }
@@ -244,7 +291,7 @@ EvaluatedGeometry<dim,Number>::get_JxW_values() const
 template <int dim, typename Number>
 inline
 const Quadrature<1>&
-EvaluatedGeometry<dim,Number>::get_quadrature() const
+MappingFEEvaluation<dim,Number>::get_quadrature() const
 {
   return quadrature_1d;
 }
index b50d23bc01887c999ea679c39fd8fa97fdc8647e..ad08fa64aaed6e5fd6a713a26fad90bd1830d09b 100644 (file)
@@ -111,6 +111,21 @@ internal_reinit(const Mapping<dim>                          &mapping,
                 const std::vector<hp::QCollection<1> >      &quad,
                 const typename MatrixFree<dim,Number>::AdditionalData additional_data)
 {
+
+  // Reads out the FE information and stores the shape function values,
+  // gradients and Hessians for quadrature points.
+  {
+    const unsigned int n_fe   = dof_handler.size();
+    const unsigned int n_quad = quad.size();
+    shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1));
+    for (unsigned int no=0; no<n_fe; no++)
+      for (unsigned int nq =0; nq<n_quad; nq++)
+        {
+          AssertDimension (quad[nq].size(), 1);
+          shape_info(no,nq,0,0).reinit(quad[nq][0], dof_handler[no]->get_fe());
+        }
+  }
+
   if (additional_data.initialize_indices == true)
     {
       clear();
@@ -191,21 +206,6 @@ internal_reinit(const Mapping<dim>                          &mapping,
         }
     }
 
-  // Reads out the FE information and stores the shape function values,
-  // gradients and Hessians for quadrature points.
-  const unsigned int n_fe   = dof_handler.size();
-  const unsigned int n_quad = quad.size();
-  shape_info.reinit (TableIndices<4>(n_fe, n_quad, 1, 1));
-  for (unsigned int no=0; no<n_fe; no++)
-    {
-      const FiniteElement<dim> &fe = dof_handler[no]->get_fe();
-      for (unsigned int nq =0; nq<n_quad; nq++)
-        {
-          AssertDimension (quad[nq].size(), 1);
-          shape_info(no,nq,0,0).reinit(quad[nq][0], fe.base_element(0));
-        }
-    }
-
   // Evaluates transformations from unit to real cell, Jacobian determinants,
   // quadrature points in real space, based on the ordering of the cells
   // determined in @p extract_local_to_global_indices. The algorithm assumes
@@ -233,6 +233,29 @@ internal_reinit(const Mapping<dim>                            &mapping,
                 const std::vector<hp::QCollection<1> >        &quad,
                 const typename MatrixFree<dim,Number>::AdditionalData additional_data)
 {
+  // Reads out the FE information and stores the shape function values,
+  // gradients and Hessians for quadrature points.
+  {
+    const unsigned int n_components = dof_handler.size();
+    const unsigned int n_quad       = quad.size();
+    unsigned int n_fe_in_collection = 0;
+    for (unsigned int i=0; i<n_components; ++i)
+      n_fe_in_collection = std::max (n_fe_in_collection,
+                                     dof_handler[i]->get_fe().size());
+    unsigned int n_quad_in_collection = 0;
+    for (unsigned int q=0; q<n_quad; ++q)
+      n_quad_in_collection = std::max (n_quad_in_collection, quad[q].size());
+    shape_info.reinit (TableIndices<4>(n_components, n_quad,
+                                       n_fe_in_collection,
+                                       n_quad_in_collection));
+    for (unsigned int no=0; no<n_components; no++)
+      for (unsigned int fe_no=0; fe_no<dof_handler[no]->get_fe().size(); ++fe_no)
+        for (unsigned int nq =0; nq<n_quad; nq++)
+          for (unsigned int q_no=0; q_no<quad[nq].size(); ++q_no)
+            shape_info(no,nq,fe_no,q_no).reinit (quad[nq][q_no],
+                                                 dof_handler[no]->get_fe()[fe_no]);
+  }
+
   if (additional_data.initialize_indices == true)
     {
       clear();
@@ -314,30 +337,6 @@ internal_reinit(const Mapping<dim>                            &mapping,
         }
     }
 
-  // Reads out the FE information and stores the shape function values,
-  // gradients and Hessians for quadrature points.
-  const unsigned int n_components = dof_handler.size();
-  const unsigned int n_quad       = quad.size();
-  unsigned int n_fe_in_collection = 0;
-  for (unsigned int i=0; i<n_components; ++i)
-    n_fe_in_collection = std::max (n_fe_in_collection,
-                                   dof_handler[i]->get_fe().size());
-  unsigned int n_quad_in_collection = 0;
-  for (unsigned int q=0; q<n_quad; ++q)
-    n_quad_in_collection = std::max (n_quad_in_collection, quad[q].size());
-  shape_info.reinit (TableIndices<4>(n_components, n_quad,
-                                     n_fe_in_collection,
-                                     n_quad_in_collection));
-  for (unsigned int no=0; no<n_components; no++)
-    for (unsigned int fe_no=0; fe_no<dof_handler[no]->get_fe().size(); ++fe_no)
-      {
-        const FiniteElement<dim> &fe = dof_handler[no]->get_fe()[fe_no];
-        for (unsigned int nq =0; nq<n_quad; nq++)
-          for (unsigned int q_no=0; q_no<quad[nq].size(); ++q_no)
-            shape_info(no,nq,fe_no,q_no).reinit (quad[nq][q_no],
-                                                 fe.base_element(0));
-      }
-
   // Evaluates transformations from unit to real cell, Jacobian determinants,
   // quadrature points in real space, based on the ordering of the cells
   // determined in @p extract_local_to_global_indices.
@@ -461,7 +460,7 @@ initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_hand
 template <int dim, typename Number>
 void MatrixFree<dim,Number>::initialize_indices
 (const std::vector<const ConstraintMatrix *> &constraint,
- const std::vector<IndexSet>                &locally_owned_set)
+ const std::vector<IndexSet>                 &locally_owned_set)
 {
   const unsigned int n_fe = dof_handlers.n_dof_handlers;
   const unsigned int n_active_cells = cell_level_index.size();
@@ -507,12 +506,11 @@ void MatrixFree<dim,Number>::initialize_indices
                                                  fes.back()->dofs_per_cell);
         }
 
-      lexicographic_inv[no].resize (fes.size());
       for (unsigned int fe_index = 0; fe_index<fes.size(); ++fe_index)
         {
           const FiniteElement<dim> &fe = *fes[fe_index];
           Assert (fe.n_base_elements() == 1,
-                  ExcMessage ("MatrixFree only works for DoFHandler with one base element"));
+                  ExcMessage ("MatrixFree currently only works for DoFHandler with one base element"));
           const unsigned int n_fe_components = fe.element_multiplicity (0);
 
           // cache number of finite elements and dofs_per_cell
@@ -521,40 +519,7 @@ void MatrixFree<dim,Number>::initialize_indices
           dof_info[no].dimension    = dim;
           dof_info[no].n_components = n_fe_components;
 
-          // get permutation that gives lexicographic renumbering of the cell
-          // dofs renumber (this is necessary for FE_Q, for example, since
-          // there the vertex DoFs come first, which is incompatible with the
-          // lexicographic ordering necessary to apply tensor products
-          // efficiently)
-
-          if (n_fe_components == 1)
-            {
-              lexicographic_inv[no][fe_index] =
-                internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe);
-              AssertDimension (lexicographic_inv[no][fe_index].size(),
-                               dof_info[no].dofs_per_cell[fe_index]);
-            }
-          else
-            {
-              // ok, we have more than one component, get the inverse
-              // permutation, invert it, sort the components one after one,
-              // and invert back
-              Assert (n_fe_components > 1, ExcInternalError());
-              std::vector<unsigned int> scalar_lex =
-                Utilities::invert_permutation(internal::MatrixFreeFunctions::get_lexicographic_numbering_inverse(fe.base_element(0)));
-              AssertDimension (scalar_lex.size() * n_fe_components,
-                               dof_info[no].dofs_per_cell[fe_index]);
-              std::vector<unsigned int> lexicographic (dof_info[no].dofs_per_cell[fe_index]);
-              for (unsigned int comp=0; comp<n_fe_components; ++comp)
-                for (unsigned int i=0; i<scalar_lex.size(); ++i)
-                  lexicographic[fe.component_to_system_index(comp,i)]
-                    = scalar_lex.size () * comp + scalar_lex[i];
-
-              // invert numbering
-              lexicographic_inv[no][fe_index] =
-                Utilities::invert_permutation(lexicographic);
-            }
-          AssertDimension (lexicographic_inv[no][fe_index].size(),
+          AssertDimension (shape_info(no,0,fe_index,0).lexicographic_numbering.size(),
                            dof_info[no].dofs_per_cell[fe_index]);
         }
 
@@ -604,7 +569,7 @@ void MatrixFree<dim,Number>::initialize_indices
               local_dof_indices.resize (dof_info[no].dofs_per_cell[0]);
               cell_it->get_dof_indices(local_dof_indices);
               dof_info[no].read_dof_indices (local_dof_indices,
-                                             lexicographic_inv[no][0],
+                                             shape_info(no,0,0,0).lexicographic_numbering,
                                              *constraint[no], counter,
                                              constraint_values,
                                              cell_at_boundary);
@@ -623,7 +588,7 @@ void MatrixFree<dim,Number>::initialize_indices
               local_dof_indices.resize (dof_info[no].dofs_per_cell[0]);
               cell_it->get_mg_dof_indices(local_dof_indices);
               dof_info[no].read_dof_indices (local_dof_indices,
-                                             lexicographic_inv[no][0],
+                                             shape_info(no,0,0,0).lexicographic_numbering,
                                              *constraint[no], counter,
                                              constraint_values,
                                              cell_at_boundary);
@@ -643,7 +608,7 @@ void MatrixFree<dim,Number>::initialize_indices
               local_dof_indices.resize (cell_it->get_fe().dofs_per_cell);
               cell_it->get_dof_indices(local_dof_indices);
               dof_info[no].read_dof_indices (local_dof_indices,
-                                             lexicographic_inv[no][cell_it->active_fe_index()],
+                                             shape_info(no,0,cell_it->active_fe_index(),0).lexicographic_numbering,
                                              *constraint[no], counter,
                                              constraint_values,
                                              cell_at_boundary);
diff --git a/deal.II/include/deal.II/matrix_free/operators.h b/deal.II/include/deal.II/matrix_free/operators.h
new file mode 100644 (file)
index 0000000..7e6821e
--- /dev/null
@@ -0,0 +1,224 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2011 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef __deal2__matrix_free_operators_h
+#define __deal2__matrix_free_operators_h
+
+
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace MatrixFreeOperators
+{
+  /**
+   * This class implements the operation of the action of the inverse of a mass
+   * matrix on an element for the special case of an evaluation object with as
+   * many quadrature points as there are cell degrees of freedom. It uses
+   * algorithms from FEEvaluation and produces the exact mass matrix for DGQ
+   * elements. This algorithm uses tensor products of inverse 1D shape matrices
+   * over quadrature points, so the inverse operation is exactly as expensive as
+   * applying the forward operator on each cell. Of course, for continuous
+   * finite elements this operation does not produce the inverse of a mass
+   * operation as the coupling between the elements cannot be considered by this
+   * operation.
+   *
+   * The equation may contain variable coefficients, so the user is required to
+   * provide an array for the inverse of the local coefficient (this class
+   * provide a helper method 'fill_inverse_JxW_values' to get the inverse of a
+   * constant-coefficient operator).
+   *
+   * @author Martin Kronbichler, 2014
+   */
+  template <int dim, int fe_degree, int n_components = 1, typename Number = double>
+  class CellwiseInverseMassMatrix
+  {
+  public:
+    /**
+     * Constructor. Initializes the shape information from the ShapeInfo field
+     * in the class FEEval.
+     */
+    CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval);
+
+    /**
+     * Applies the inverse mass matrix operation on an input array. It is
+     * assumed that the passed input and output arrays are of correct size,
+     * namely FEEval::dofs_per_cell * n_components long. The inverse of the
+     * local coefficient (also containing the inverse JxW values) must be passed
+     * as first argument. Passing more than one component in the coefficient is
+     * allowed.
+     */
+    void apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficient,
+               const unsigned int             n_actual_components,
+               const VectorizedArray<Number> *in_array,
+               VectorizedArray<Number>       *out_array) const;
+
+    /**
+     * Fills the given array with the inverse of the JxW values, i.e., a mass
+     * matrix with coefficient 1. Non-unit coefficients must be multiplied (in
+     * inverse form) to this array.
+     */
+    void fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const;
+
+  private:
+    /**
+     * A reference to the FEEvaluation object for getting the JxW_values.
+     */
+    const FEEvaluationBase<dim,n_components,Number> &fe_eval;
+
+    /**
+     * A structure to hold inverse shape functions, put into an aligned vector
+     * in order to make sure it can be allocated.
+     */
+    struct Inverse1DShape
+    {
+      VectorizedArray<Number> evenodd[fe_degree+1][fe_degree/2+1];
+    };
+    AlignedVector<Inverse1DShape> inverse_shape;
+  };
+
+
+
+  // ------------------------------------ inline functions ---------------------
+
+  template <int dim, int fe_degree, int n_components, typename Number>
+  inline
+  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+  ::CellwiseInverseMassMatrix (const FEEvaluationBase<dim,n_components,Number> &fe_eval)
+    :
+    fe_eval (fe_eval)
+  {
+    FullMatrix<double> shapes_1d(fe_degree+1, fe_degree+1);
+    for (unsigned int i=0, c=0; i<shapes_1d.m(); ++i)
+      for (unsigned int j=0; j<shapes_1d.n(); ++j, ++c)
+        shapes_1d(i,j) = fe_eval.get_shape_info().shape_values_number[c];
+    shapes_1d.gauss_jordan();
+    inverse_shape.resize(1);
+    for (int i=0; i<(fe_degree+1)/2; ++i)
+      for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
+        {
+          inverse_shape[0].evenodd[i][q] =
+            0.5 * (shapes_1d(i,q) + shapes_1d(i,fe_degree-q));
+          inverse_shape[0].evenodd[fe_degree-i][q] =
+            0.5 * (shapes_1d(i,q) - shapes_1d(i,fe_degree-q));
+        }
+    if (fe_degree % 2 == 0)
+      for (unsigned int q=0; q<(fe_degree+2)/2; ++q)
+        inverse_shape[0].evenodd[fe_degree/2][q] = shapes_1d(fe_degree/2,q);
+  }
+
+
+
+  template <int dim, int fe_degree, int n_components, typename Number>
+  inline
+  void
+  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+  ::fill_inverse_JxW_values(AlignedVector<VectorizedArray<Number> > &inverse_jxw) const
+  {
+    const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+    Assert(inverse_jxw.size() > 0 &&
+           inverse_jxw.size() % dofs_per_cell == 0,
+           ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));
+
+    // temporarily reduce size of inverse_jxw to dofs_per_cell to get JxW values
+    // from fe_eval (will not reallocate any memory)
+    const unsigned int previous_size = inverse_jxw.size();
+    inverse_jxw.resize(dofs_per_cell);
+    fe_eval.fill_JxW_values(inverse_jxw);
+
+    // invert
+    inverse_jxw.resize_fast(previous_size);
+    for (unsigned int q=0; q<dofs_per_cell; ++q)
+      inverse_jxw[q] = 1./inverse_jxw[q];
+    // copy values to rest of vector
+    for (unsigned int q=dofs_per_cell; q<previous_size; )
+      for (unsigned int i=0; i<dofs_per_cell; ++i, ++q)
+        inverse_jxw[q] = inverse_jxw[i];
+  }
+
+
+
+  template <int dim, int fe_degree, int n_components, typename Number>
+  inline
+  void
+  CellwiseInverseMassMatrix<dim,fe_degree,n_components,Number>
+  ::apply(const AlignedVector<VectorizedArray<Number> > &inverse_coefficients,
+          const unsigned int             n_actual_components,
+          const VectorizedArray<Number> *in_array,
+          VectorizedArray<Number>       *out_array) const
+  {
+    const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
+    Assert(inverse_coefficients.size() > 0 &&
+           inverse_coefficients.size() % dofs_per_cell == 0,
+           ExcMessage("Expected diagonal to be a multiple of scalar dof per cells"));
+    if (inverse_coefficients.size() != dofs_per_cell)
+      AssertDimension(n_actual_components * dofs_per_cell, inverse_coefficients.size());
+
+    Assert(dim == 2 || dim == 3, ExcNotImplemented());
+
+    const unsigned int shift_coefficient =
+      inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0;
+    const VectorizedArray<Number> *inv_coefficient = &inverse_coefficients[0];
+    VectorizedArray<Number> temp_data_field[dofs_per_cell];
+    for (unsigned int d=0; d<n_components; ++d)
+      {
+        const VectorizedArray<Number>* in = in_array+d*dofs_per_cell;
+        VectorizedArray<Number>* out = out_array+d*dofs_per_cell;
+        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                               VectorizedArray<Number>,0,false,false,2>
+          (inverse_shape[0].evenodd, in, temp_data_field);
+        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                               VectorizedArray<Number>,1,false,false,2>
+          (inverse_shape[0].evenodd, temp_data_field, out);
+
+        if (dim == 3)
+          {
+            internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                                   VectorizedArray<Number>,2,false,false,2>
+              (inverse_shape[0].evenodd, out, temp_data_field);
+            for (unsigned int q=0; q<dofs_per_cell; ++q)
+              temp_data_field[q] *= inv_coefficient[q];
+            internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                                   VectorizedArray<Number>,2,true,false,2>
+              (inverse_shape[0].evenodd, temp_data_field, out);
+          }
+        else if (dim == 2)
+          for (unsigned int q=0; q<dofs_per_cell; ++q)
+            out[q] *= inv_coefficient[q];
+
+        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                               VectorizedArray<Number>,1,true,false,2>
+          (inverse_shape[0].evenodd, out, temp_data_field);
+        internal::apply_tensor_product_evenodd<dim,fe_degree,fe_degree+1,
+                                               VectorizedArray<Number>,0,true,false,2>
+          (inverse_shape[0].evenodd, temp_data_field, out);
+
+        inv_coefficient += shift_coefficient;
+      }
+  }
+
+} // end of namespace MatrixFreeOperators
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 7b6057a43b6617f6d68d01c0222fdf7336736403..3164b4679ce354c25291a109945f2eda59aa610b 100644 (file)
@@ -155,6 +155,15 @@ namespace internal
        */
       std::vector<Number>    shape_gradient_number;
 
+      /**
+       * Renumbering from deal.II's numbering of cell degrees of freedom to
+       * lexicographic numbering used inside the FEEvaluation schemes of the
+       * underlying element in the DoFHandler. For vector-valued elements, the
+       * renumbering starts with a lexicographic numbering of the first
+       * component, then everything of the second component, and so on.
+       */
+      std::vector<unsigned int> lexicographic_numbering;
+
       /**
        * Stores the number of quadrature points in
        * @p dim dimensions for a cell.
@@ -181,16 +190,6 @@ namespace internal
     };
 
 
-    /**
-     * Extracts the mapping between the actual numbering inside a scalar
-     * element and lexicographic numbering as required by the evaluation
-     * routines.
-     */
-    template <int dim>
-    std::vector<unsigned int>
-    get_lexicographic_numbering_inverse (const FiniteElement<dim> &fe);
-
-
 
     // ------------------------------------------ inline functions
 
index d11c8654e07b83da00fb42e4217bb8e7d2523055..eb8e2bda753e9f9ade4acea0970635a58eacb904 100644 (file)
@@ -34,43 +34,6 @@ namespace internal
   namespace MatrixFreeFunctions
   {
 
-    // helper function
-    template <int dim>
-    std::vector<unsigned int>
-    get_lexicographic_numbering_inverse(const FiniteElement<dim> &fe)
-    {
-      Assert(fe.n_components() == 1,
-             ExcMessage("Expected a scalar element"));
-
-      const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
-        dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(&fe);
-
-      const FE_Poly<TensorProductPolynomials<dim,Polynomials::
-        PiecewisePolynomial<double> >,dim,dim> *fe_poly_piece =
-        dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
-        Polynomials::PiecewisePolynomial<double> >,dim,dim>*> (&fe);
-
-      const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(&fe);
-
-      std::vector<unsigned int> lexicographic;
-      if (fe_poly != 0)
-        lexicographic = fe_poly->get_poly_space_numbering_inverse();
-      else if (fe_poly_piece != 0)
-        lexicographic = fe_poly_piece->get_poly_space_numbering_inverse();
-      else if (fe_dgp != 0)
-        {
-          lexicographic.resize(fe_dgp->dofs_per_cell);
-          for (unsigned int i=0; i<fe_dgp->dofs_per_cell; ++i)
-            lexicographic[i] = i;
-        }
-      else
-        Assert(false, ExcNotImplemented());
-
-      return lexicographic;
-    }
-
-
-
     // ----------------- actual ShapeInfo functions --------------------
 
     template <typename Number>
@@ -99,19 +62,67 @@ namespace internal
 
       const unsigned int n_dofs_1d = fe->degree+1,
                          n_q_points_1d = quad.size();
-      std::vector<unsigned int> lexicographic (fe->dofs_per_cell);
 
       // renumber (this is necessary for FE_Q, for example, since there the
       // vertex DoFs come first, which is incompatible with the lexicographic
       // ordering necessary to apply tensor products efficiently)
+      std::vector<unsigned int> scalar_lexicographic;
       {
-        lexicographic = get_lexicographic_numbering_inverse(*fe);
+        // find numbering to lexicographic
+        Assert(fe->n_components() == 1,
+               ExcMessage("Expected a scalar element"));
+
+        const FE_Poly<TensorProductPolynomials<dim>,dim,dim> *fe_poly =
+          dynamic_cast<const FE_Poly<TensorProductPolynomials<dim>,dim,dim>*>(fe);
+
+        const FE_Poly<TensorProductPolynomials<dim,Polynomials::
+          PiecewisePolynomial<double> >,dim,dim> *fe_poly_piece =
+          dynamic_cast<const FE_Poly<TensorProductPolynomials<dim,
+          Polynomials::PiecewisePolynomial<double> >,dim,dim>*> (fe);
+
+        const FE_DGP<dim> *fe_dgp = dynamic_cast<const FE_DGP<dim>*>(fe);
+
+        if (fe_poly != 0)
+          scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse();
+        else if (fe_poly_piece != 0)
+          scalar_lexicographic = fe_poly_piece->get_poly_space_numbering_inverse();
+        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;
+          }
+        else
+          Assert(false, ExcNotImplemented());
+
+        // Finally store the renumbering into the member variable of this
+        // class
+        if (fe_in.n_components() == 1)
+          lexicographic_numbering = scalar_lexicographic;
+        else
+          {
+            // have more than one component, get the inverse
+            // permutation, invert it, sort the components one after one,
+            // and invert back
+            std::vector<unsigned int> scalar_inv =
+              Utilities::invert_permutation(scalar_lexicographic);
+            std::vector<unsigned int> lexicographic (fe_in.dofs_per_cell);
+            for (unsigned int comp=0; comp<fe_in.n_components(); ++comp)
+              for (unsigned int i=0; i<scalar_inv.size(); ++i)
+                lexicographic[fe_in.component_to_system_index(comp,i)]
+                  = scalar_inv.size () * comp + scalar_inv[i];
+
+            // invert numbering again
+            lexicographic_numbering =
+              Utilities::invert_permutation(lexicographic);
+          }
 
         // to evaluate 1D polynomials, evaluate along the line where y=z=0,
         // assuming that shape_value(0,Point<dim>()) == 1. otherwise, need
         // other entry point (e.g. generating a 1D element by reading the
         // name, as done before r29356)
-        Assert(std::fabs(fe->shape_value(lexicographic[0], Point<dim>())-1) < 1e-13,
+        Assert(std::fabs(fe->shape_value(scalar_lexicographic[0],
+                                         Point<dim>())-1) < 1e-13,
                ExcInternalError());
       }
 
@@ -138,7 +149,7 @@ namespace internal
         {
           // need to reorder from hierarchical to lexicographic to get the
           // DoFs correct
-          const unsigned int my_i = lexicographic[i];
+          const unsigned int my_i = scalar_lexicographic[i];
           for (unsigned int q=0; q<n_q_points_1d; ++q)
             {
               // fill both vectors with
diff --git a/tests/matrix_free/inverse_mass_01.cc b/tests/matrix_free/inverse_mass_01.cc
new file mode 100644 (file)
index 0000000..77faf47
--- /dev/null
@@ -0,0 +1,216 @@
+// ---------------------------------------------------------------------
+// $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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on DG elements by comparing its action on a
+// random vector to a CG solver
+
+#include "../tests.h"
+#include <deal.II/base/function.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/fe/mapping_q.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_mass_operator (const MatrixFree<dim,Number>  &data,
+                       VECTOR       &dst,
+                       const VECTOR &src,
+                       const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<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, false);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          fe_eval.submit_value (fe_eval.get_value(q),q);
+        fe_eval.integrate (true, false);
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void
+  local_inverse_mass_operator (const MatrixFree<dim,Number>  &data,
+                               VECTOR       &dst,
+                               const VECTOR &src,
+                               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval (data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,1,Number> mass_inv(fe_eval);
+    const unsigned int n_q_points = fe_eval.n_q_points;
+    AlignedVector<VectorizedArray<Number> > inverse_coefficients(n_q_points);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        fe_eval.reinit (cell);
+        mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+        fe_eval.read_dof_values (src);
+        mass_inv.apply(inverse_coefficients, 1, fe_eval.begin_dof_values(),
+                       fe_eval.begin_dof_values());
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+                    this, dst, src);
+  };
+
+  void apply_inverse (VECTOR       &dst,
+                      const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  MappingQ<dim> mapping(4);
+
+  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;
+    ConstraintMatrix constraints;
+
+    mf_data.reinit (mapping, dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+    }
+
+  mf.apply_inverse (inverse, in);
+
+  SolverControl control(1000, 1e-12);
+  SolverCG<Vector<number> > solver(control);
+  solver.solve (mf, reference, in, PreconditionIdentity());
+
+  inverse -= reference;
+  const double diff_norm = inverse.linfty_norm() / reference.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_DGQ<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  do_test<dim, fe_degree, double> (dof);
+
+  if (dim == 2)
+    {
+      deallog.push("float");
+      deallog.threshold_double(1.e-6);
+      do_test<dim, fe_degree, float> (dof);
+      deallog.threshold_double(5.e-11);
+      deallog.pop();
+    }
+}
+
+
+
+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>();
+    test<2,4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/inverse_mass_01.output b/tests/matrix_free/inverse_mass_01.output
new file mode 100644 (file)
index 0000000..98c2aaf
--- /dev/null
@@ -0,0 +1,41 @@
+
+DEAL:2d::Testing FE_DGQ<2>(1)
+DEAL:2d:cg::Starting value 6.32
+DEAL:2d:cg::Convergence step 51 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(1)
+DEAL:2d:float:cg::Starting value 6.15
+DEAL:2d:float:cg::Convergence step 66 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_DGQ<2>(2)
+DEAL:2d:cg::Starting value 8.62
+DEAL:2d:cg::Convergence step 116 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(2)
+DEAL:2d:float:cg::Starting value 9.05
+DEAL:2d:float:cg::Convergence step 158 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_DGQ<2>(4)
+DEAL:2d:cg::Starting value 14.6
+DEAL:2d:cg::Convergence step 354 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_DGQ<2>(4)
+DEAL:2d:float:cg::Starting value 14.7
+DEAL:2d:float:cg::Convergence step 437 value 0
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:3d::Testing FE_DGQ<3>(1)
+DEAL:3d:cg::Starting value 13.7
+DEAL:3d:cg::Convergence step 131 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FE_DGQ<3>(2)
+DEAL:3d:cg::Starting value 13.5
+DEAL:3d:cg::Convergence step 255 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/inverse_mass_02.cc b/tests/matrix_free/inverse_mass_02.cc
new file mode 100644 (file)
index 0000000..e3efb75
--- /dev/null
@@ -0,0 +1,205 @@
+// ---------------------------------------------------------------------
+// $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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on vector DG elements, otherwise the same
+// as inverse_mass_01.cc
+
+#include "../tests.h"
+#include <deal.II/base/function.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_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_mass_operator (const MatrixFree<dim,Number>  &data,
+                       VECTOR       &dst,
+                       const VECTOR &src,
+                       const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,dim,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, false);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          fe_eval.submit_value (fe_eval.get_value(q),q);
+        fe_eval.integrate (true, false);
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void
+  local_inverse_mass_operator (const MatrixFree<dim,Number>  &data,
+                               VECTOR       &dst,
+                               const VECTOR &src,
+                               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval (data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,dim,Number> mass_inv(fe_eval);
+    const unsigned int n_q_points = fe_eval.n_q_points;
+    AlignedVector<VectorizedArray<Number> > inverse_coefficients(n_q_points);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        fe_eval.reinit (cell);
+        mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+        fe_eval.read_dof_values (src);
+        mass_inv.apply(inverse_coefficients, dim, fe_eval.begin_dof_values(),
+                       fe_eval.begin_dof_values());
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+                    this, dst, src);
+  };
+
+  void apply_inverse (VECTOR       &dst,
+                      const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+  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;
+    ConstraintMatrix constraints;
+
+    mf_data.reinit (dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+    }
+
+  mf.apply_inverse (inverse, in);
+
+  SolverControl control(1000, 1e-12);
+  SolverCG<Vector<number> > solver(control);
+  solver.solve (mf, reference, in, PreconditionIdentity());
+
+  inverse -= reference;
+  const double diff_norm = inverse.linfty_norm() / reference.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_cube (tria, -1, 1);
+  tria.refine_global(1);
+  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();
+
+  FESystem<dim> fe (FE_DGQ<dim>(fe_degree), dim);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  do_test<dim, fe_degree, double> (dof);
+}
+
+
+
+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>();
+    test<2,4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/inverse_mass_02.output b/tests/matrix_free/inverse_mass_02.output
new file mode 100644 (file)
index 0000000..67b9f69
--- /dev/null
@@ -0,0 +1,26 @@
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^2]
+DEAL:2d:cg::Starting value 7.90
+DEAL:2d:cg::Convergence step 7 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^2]
+DEAL:2d:cg::Starting value 11.6
+DEAL:2d:cg::Convergence step 15 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^2]
+DEAL:2d:cg::Starting value 19.1
+DEAL:2d:cg::Convergence step 46 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d:cg::Starting value 25.0
+DEAL:3d:cg::Convergence step 10 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d:cg::Starting value 24.3
+DEAL:3d:cg::Convergence step 32 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
diff --git a/tests/matrix_free/inverse_mass_03.cc b/tests/matrix_free/inverse_mass_03.cc
new file mode 100644 (file)
index 0000000..02abac7
--- /dev/null
@@ -0,0 +1,220 @@
+// ---------------------------------------------------------------------
+// $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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Tests CellwiseInverseMassMatrix on vector DG elements, similar test as
+// inverse_mass_02 but using different coefficients on different components
+
+#include "../tests.h"
+#include <deal.II/base/function.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/fe/mapping_q.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/precondition.h>
+
+#include <deal.II/matrix_free/operators.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree, typename Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+
+  MatrixFreeTest(const MatrixFree<dim,Number> &data_in):
+    data (data_in)
+  {};
+
+  void
+  local_mass_operator (const MatrixFree<dim,Number>  &data,
+                       VECTOR       &dst,
+                       const VECTOR &src,
+                       const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,3,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, false);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            Tensor<1,3,VectorizedArray<double> > val = fe_eval.get_value(q);
+            val[0] *= make_vectorized_array(0.8314159);
+            val[1] *= make_vectorized_array(2.3);
+            val[2] *= make_vectorized_array(1.98);
+            fe_eval.submit_value (val, q);
+          }
+        fe_eval.integrate (true, false);
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void
+  local_inverse_mass_operator (const MatrixFree<dim,Number>  &data,
+                               VECTOR       &dst,
+                               const VECTOR &src,
+                               const std::pair<unsigned int,unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim,fe_degree,fe_degree+1,3,Number> fe_eval (data);
+    MatrixFreeOperators::CellwiseInverseMassMatrix<dim,fe_degree,3,Number> mass_inv(fe_eval);
+    const unsigned int n_q_points = fe_eval.n_q_points;
+    AlignedVector<VectorizedArray<Number> > inverse_coefficients(3*n_q_points);
+
+    for (unsigned int cell=cell_range.first; cell<cell_range.second; ++cell)
+      {
+        fe_eval.reinit (cell);
+        mass_inv.fill_inverse_JxW_values(inverse_coefficients);
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            inverse_coefficients[q] *= make_vectorized_array(1./0.8314159);
+            inverse_coefficients[n_q_points+q] *= make_vectorized_array(1./2.3);
+            inverse_coefficients[2*n_q_points+q] *= make_vectorized_array(1./1.98);
+          }
+        fe_eval.read_dof_values (src);
+        mass_inv.apply(inverse_coefficients, 3, fe_eval.begin_dof_values(),
+                       fe_eval.begin_dof_values());
+        fe_eval.distribute_local_to_global (dst);
+      }
+  }
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_mass_operator,
+                    this, dst, src);
+  };
+
+  void apply_inverse (VECTOR       &dst,
+                      const VECTOR &src) const
+  {
+    dst = 0;
+    data.cell_loop (&MatrixFreeTest<dim,fe_degree,Number,VECTOR>::local_inverse_mass_operator,
+                    this, dst, src);
+  };
+
+private:
+  const MatrixFree<dim,Number> &data;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MappingQ<dim> mapping(4);
+  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;
+    ConstraintMatrix constraints;
+
+    mf_data.reinit (mapping, dof, constraints, quad, data);
+  }
+
+  MatrixFreeTest<dim,fe_degree,number> mf (mf_data);
+  Vector<number> in (dof.n_dofs()), inverse (dof.n_dofs()), reference(dof.n_dofs());
+
+  for (unsigned int i=0; i<dof.n_dofs(); ++i)
+    {
+      const double entry = Testing::rand()/(double)RAND_MAX;
+      in(i) = entry;
+    }
+
+  mf.apply_inverse (inverse, in);
+
+  SolverControl control(10000, 1e-12);
+  SolverCG<Vector<number> > solver(control);
+  solver.solve (mf, reference, in, PreconditionIdentity());
+
+  inverse -= reference;
+  const double diff_norm = inverse.linfty_norm() / reference.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();
+
+  FESystem<dim> fe (FE_DGQ<dim>(fe_degree),3);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+
+  do_test<dim, fe_degree, double> (dof);
+}
+
+
+
+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>();
+    test<2,4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/inverse_mass_03.output b/tests/matrix_free/inverse_mass_03.output
new file mode 100644 (file)
index 0000000..3579912
--- /dev/null
@@ -0,0 +1,26 @@
+
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(1)^3]
+DEAL:2d:cg::Starting value 10.6
+DEAL:2d:cg::Convergence step 108 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(2)^3]
+DEAL:2d:cg::Starting value 15.1
+DEAL:2d:cg::Convergence step 252 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_DGQ<2>(4)^3]
+DEAL:2d:cg::Starting value 25.4
+DEAL:2d:cg::Convergence step 699 value 0
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(1)^3]
+DEAL:3d:cg::Starting value 23.7
+DEAL:3d:cg::Convergence step 412 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_DGQ<3>(2)^3]
+DEAL:3d:cg::Starting value 23.7
+DEAL:3d:cg::Convergence step 1019 value 0
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
index 7f0bd5ac0a193af2846769880effb21ed424a7c1..80126fc3ea8315a3e038629304acc3c6bb285390 100644 (file)
@@ -104,16 +104,10 @@ private:
 
 template <int dim, int fe_degree, typename number>
 void do_test (const DoFHandler<dim> &dof,
-              const ConstraintMatrix &constraints,
-              const unsigned int     parallel_option = 0)
+              const ConstraintMatrix &constraints)
 {
 
   deallog << "Testing " << dof.get_fe().get_name() << std::endl;
-  if (parallel_option > 0)
-    deallog << "Parallel option: " << parallel_option << std::endl;
-  //std::cout << "Number of cells: " << dof.get_tria().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;
 
   MatrixFree<dim,number> mf_data;
   {
diff --git a/tests/matrix_free/matrix_vector_15.cc b/tests/matrix_free/matrix_vector_15.cc
new file mode 100644 (file)
index 0000000..0c8f9aa
--- /dev/null
@@ -0,0 +1,242 @@
+// ---------------------------------------------------------------------
+// $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 on
+// the fly from an FEValues implementation instead of data cached in
+// MatrixFree.
+
+#include "../tests.h"
+#include <deal.II/base/function.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 Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const DoFHandler<dim> &dof_handler,
+                 const ConstraintMatrix &constraints)
+    :
+    dof_handler (dof_handler),
+    constraints (constraints)
+  {}
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    MappingFEEvaluation<dim,Number> mapping(QGauss<1>(fe_degree+1),
+                                            update_gradients | update_values |
+                                            update_JxW_values);
+    VECTOR src_cpy = src;
+    constraints.distribute(src_cpy);
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,Number> fe_eval(mapping, dof_handler);
+    dst = 0;
+    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for ( ; cell != endc; ++cell)
+      {
+        mapping.reinit(cell);
+        fe_eval.read_dof_values (src_cpy);
+        fe_eval.evaluate (true, true, false);
+        for (unsigned int q=0; q<fe_eval.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);
+      }
+    constraints.condense(dst);
+  };
+
+private:
+  const DoFHandler<dim> &dof_handler;
+  const ConstraintMatrix &constraints;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim> &dof,
+              const ConstraintMatrix &constraints,
+              const unsigned int     parallel_option = 0)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  if (parallel_option > 0)
+    deallog << "Parallel option: " << parallel_option << std::endl;
+  //std::cout << "Number of cells: " << dof.get_tria().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;
+
+  MatrixFreeTest<dim,fe_degree,number> mf (dof, constraints);
+  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<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_test<dim, fe_degree, double> (dof, constraints);
+
+  if (dim == 2)
+    {
+      deallog.push("float");
+      deallog.threshold_double(1.e-6);
+      do_test<dim, fe_degree, float> (dof, constraints);
+      deallog.threshold_double(5.e-11);
+      deallog.pop();
+    }
+}
+
+
+
+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>();
+    test<2,4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_15.output b/tests/matrix_free/matrix_vector_15.output
new file mode 100644 (file)
index 0000000..0cd9474
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(1)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(2)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d:float::Testing FE_Q<2>(4)
+DEAL:2d:float::Norm of difference: 0
+DEAL:2d:float::
+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_16.cc b/tests/matrix_free/matrix_vector_16.cc
new file mode 100644 (file)
index 0000000..5b2e37f
--- /dev/null
@@ -0,0 +1,237 @@
+// ---------------------------------------------------------------------
+// $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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// evaluating the geometry on the fly as in matrix_vector_15, but using
+// vector-valued finite elements.
+
+#include "../tests.h"
+#include <deal.II/base/function.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_system.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 Number, typename VECTOR=Vector<Number> >
+class MatrixFreeTest
+{
+public:
+  MatrixFreeTest(const DoFHandler<dim> &dof_handler,
+                 const ConstraintMatrix &constraints)
+    :
+    dof_handler (dof_handler),
+    constraints (constraints)
+  {}
+
+  void vmult (VECTOR       &dst,
+              const VECTOR &src) const
+  {
+    MappingFEEvaluation<dim,Number> mapping(QGauss<1>(fe_degree+1),
+                                            update_gradients | update_values |
+                                            update_JxW_values);
+    VECTOR src_cpy = src;
+    constraints.distribute(src_cpy);
+    FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval(mapping, dof_handler);
+    dst = 0;
+    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+      endc = dof_handler.end();
+    for ( ; cell != endc; ++cell)
+      {
+        mapping.reinit(cell);
+        fe_eval.read_dof_values (src_cpy);
+        fe_eval.evaluate (true, true, false);
+        for (unsigned int q=0; q<fe_eval.n_q_points; ++q)
+          {
+            fe_eval.submit_value (make_vectorized_array<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);
+      }
+    constraints.condense(dst);
+  };
+
+private:
+  const DoFHandler<dim> &dof_handler;
+  const ConstraintMatrix &constraints;
+};
+
+
+
+template <int dim, int fe_degree, typename number>
+void do_test (const DoFHandler<dim>  &dof,
+              const ConstraintMatrix &constraints,
+              const unsigned int      parallel_option = 0)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+  if (parallel_option > 0)
+    deallog << "Parallel option: " << parallel_option << std::endl;
+  //std::cout << "Number of cells: " << dof.get_tria().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;
+
+  MatrixFreeTest<dim,fe_degree,number> mf (dof, constraints);
+  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 vector-valued form of (\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)
+                if (dof.get_fe().system_to_component_index(i).first ==
+                    dof.get_fe().system_to_component_index(j).first)
+                  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();
+
+  FESystem<dim> fe(FE_Q<dim>(fe_degree), dim);
+  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>(dim),
+                                            constraints);
+  constraints.close();
+
+  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>();
+    test<2,4>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3,1>();
+    test<3,2>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/matrix_vector_16.output b/tests/matrix_free/matrix_vector_16.output
new file mode 100644 (file)
index 0000000..b1dbf5d
--- /dev/null
@@ -0,0 +1,16 @@
+
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(1)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(2)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(4)^2]
+DEAL:2d::Norm of difference: 0
+DEAL:2d::
+DEAL:3d::Testing FESystem<3>[FE_Q<3>(1)^3]
+DEAL:3d::Norm of difference: 0
+DEAL:3d::
+DEAL:3d::Testing FESystem<3>[FE_Q<3>(2)^3]
+DEAL:3d::Norm of difference: 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.