]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Initialization of FEEvaluation like FEValues 258/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 20 Nov 2014 14:36:18 +0000 (15:36 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 22 Nov 2014 07:21:27 +0000 (08:21 +0100)
A new way to initialize an FEEvaluation object (fast tensor product evaluation routines) has been introduced.
It takes the same arguments as the initialization method of FEValues, namely a Mapping, a FiniteElement, a
one-dimensional Quadrature<1> object, and update flags. This enables using FEEvaluation on the fly just as
FEValues. Since FEEvaluation treats objects of one base element only, this class can select which components
to read.

The concept is explained in the documentation of the FEEvaluation class and an example for the Laplace matrix
is given. However, it also works for the Stokes operator as shown by the test program assemble_matrix_02.

13 files changed:
doc/news/changes.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/mapping_data_on_the_fly.h [new file with mode: 0644]
include/deal.II/matrix_free/mapping_fe_evaluation.h [deleted file]
include/deal.II/matrix_free/shape_info.templates.h
tests/matrix_free/assemble_matrix_01.cc [new file with mode: 0644]
tests/matrix_free/assemble_matrix_01.output [new file with mode: 0644]
tests/matrix_free/assemble_matrix_02.cc [new file with mode: 0644]
tests/matrix_free/assemble_matrix_02.output [new file with mode: 0644]
tests/matrix_free/assemble_matrix_03.cc [new file with mode: 0644]
tests/matrix_free/assemble_matrix_03.output [new file with mode: 0644]
tests/matrix_free/matrix_vector_15.cc
tests/matrix_free/matrix_vector_16.cc

index 98d5e1ba562444379c5234137e8d5abf8b08d742..56f80def735242b823cf5890858119231083b9f1 100644 (file)
@@ -333,6 +333,15 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: The class FEEvaluation with its fast tensor evaluation routines
+  can now be initialized from a mapping, a finite element, a quadrature, and
+  update flags on the fly similar to FEValues. This provides an alternative
+  entry point to these routines without the need to create a MatrixFree data
+  object.
+  <br>
+  (Martin Kronbichler, 2014/11/20)
+  </li>
+
   <li> New: Complete the interface for periodic boundary conditions.
   Now, it is possible to specify an additional (geometric) rotation matrix along
   with vector components that should be rotated when applying periodic boundary
index 9dfa0eabcf4ef2dee697b119081301a1bccf25a7..c63d3ec5dc2bc4ffa2b2fa68e6b89920fc283d82 100644 (file)
 #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/mapping_fe_evaluation.h>
+#include <deal.II/matrix_free/mapping_data_on_the_fly.h>
 
 
 DEAL_II_NAMESPACE_OPEN
 
+
+
+// forward declarations
 namespace parallel
 {
   namespace distributed
@@ -38,15 +41,13 @@ namespace parallel
     template <typename> class Vector;
   }
 }
-
-
-
-// forward declarations
 namespace internal
 {
   DeclException0 (ExcAccessToUninitializedField);
 }
 
+template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
+          int n_components_ = 1, typename Number = double > class FEEvaluation;
 
 
 /**
@@ -88,14 +89,42 @@ public:
    */
   //@{
   /**
-   * Initializes the operation pointer to the current cell. Unlike the
-   * FEValues::reinit function, where the information related to a particular
+   * Initializes the operation pointer to the current cell. Unlike the reinit
+   * functions taking a cell iterator as argument below and the
+   * FEValues::reinit() methods, where the information related to a particular
    * cell is generated in the reinit call, this function is very cheap since
    * all data is pre-computed in @p matrix_free, and only a few indices have
    * to be set appropriately.
    */
   void reinit (const unsigned int cell);
 
+  /**
+   * Initialize the data to the current cell using a TriaIterator object as
+   * usual in FEValues. The argument is either of type
+   * DoFHandler::active_cell_iterator or DoFHandler::level_cell_iterator. This
+   * option is only available if the FEEvaluation object was created with a
+   * finite element, quadrature formula and correct update flags and
+   * <b>without</b> a MatrixFree object. This initialization method loses the
+   * ability to use vectorization, see also the description of the
+   * FEEvaluation class. When this reinit method is used, FEEvaluation can
+   * also read from vectors (but less efficient than with data coming from
+   * MatrixFree).
+   */
+  template <class DH, bool level_dof_access>
+  void reinit (const TriaIterator<DoFCellAccessor<DH,level_dof_access> > &cell);
+
+  /**
+   * Initialize the data to the current cell using a TriaIterator object as
+   * usual in FEValues. This option is only available if the FEEvaluation
+   * object was created with a finite element, quadrature formula and correct
+   * update flags and <b>without</b> a MatrixFree object. This initialization
+   * method loses the ability to use vectorization, see also the description
+   * of the FEEvaluation class. When this reinit method is used, FEEvaluation
+   * can <b>not</b> read from vectors because no DoFHandler information is
+   * available.
+   */
+  void reinit (const typename Triangulation<dim>::cell_iterator &cell);
+
   /**
    * For the transformation information stored in MappingInfo, this function
    * returns the index which belongs to the current cell as specified in @p
@@ -148,13 +177,14 @@ public:
    * 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).
+   * If this class was constructed without a MatrixFree object and the
+   * information is acquired on the fly through a
+   * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+   * class and this function extracts the values of the underlying components
+   * on the given 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::n_array_elements same copies are worked on).
    */
   template <typename VectorType>
   void read_dof_values (const VectorType &src);
@@ -191,14 +221,14 @@ public:
    * 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.
+   * If this class was constructed without a MatrixFree object and the
+   * information is acquired on the fly through a
+   * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+   * class and this function extracts the values of the underlying components
+   * on the given 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::n_array_elements same copies are worked on).
    */
   template <typename VectorType>
   void read_dof_values_plain (const VectorType &src);
@@ -230,14 +260,14 @@ public:
    * 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.
+   * If this class was constructed without a MatrixFree object and the
+   * information is acquired on the fly through a
+   * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+   * class and this function extracts the values of the underlying components
+   * on the given 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::n_array_elements same copies are worked on).
    */
   template<typename VectorType>
   void distribute_local_to_global (VectorType &dst) const;
@@ -270,14 +300,14 @@ public:
    * 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.
+   * If this class was constructed without a MatrixFree object and the
+   * information is acquired on the fly through a
+   * DoFHandler<dim>::cell_iterator, only one single cell is used by this
+   * class and this function extracts the values of the underlying components
+   * on the given 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::n_array_elements same copies are worked on).
    */
   template<typename VectorType>
   void set_dof_values (VectorType &dst) const;
@@ -539,6 +569,14 @@ public:
    */
   VectorizedArray<Number> *begin_hessians ();
 
+  /**
+   * Returns the numbering of local degrees of freedom within the evaluation
+   * routines of FEEvaluation in terms of the standard numbering on finite
+   * elements.
+   */
+  const std::vector<unsigned int> &
+  get_internal_dof_numbering() const;
+
   //@}
 
 protected:
@@ -558,30 +596,51 @@ protected:
 
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
+   * FEValues. The arguments are similar to the ones passed to the constructor
+   * of FEValues, with the notable difference that FEEvaluation expects a
+   * one-dimensional quadrature formula, Quadrature<1>, instead of a @p dim
+   * dimensional one. The finite element can be both scalar or vector valued,
+   * but this method always only selects a scalar base element at a time (with
+   * @p n_components copies as specified by the class template argument). For
+   * vector-valued elements, the optional argument @p first_selected_component
+   * allows to specify the index of the base element to be used for
+   * evaluation. Note that the internal data structures always assume that the
+   * base element is primitive, non-primitive are not supported currently.
    *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   */
-  FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
-                    const DoFHandler<dim>                 &dof_handler,
-                    const unsigned int                     first_selected_component = 0);
-
-  /**
-   * Copy constructor
+   * As known from FEValues, a call to the reinit method with a
+   * Triangulation::cell_iterator is necessary to make the geometry and
+   * degrees of freedom of the current class known. If the iterator includes
+   * DoFHandler information (i.e., it is a DoFHandler::cell_iterator or
+   * similar), the initialization allows to also read from or write to vectors
+   * in the standard way for DoFHandler::active_cell_iterator types for
+   * one cell at a time. However, this approach is much slower than the path
+   * with MatrixFree with MPI since index translation has to be done. As only
+   * one cell at a time is used, this method does not vectorize over several
+   * elements (which is most efficient for vector operations), but only
+   * possibly within the element if the evaluate/integrate routines are
+   * combined inside user code (e.g. for computing cell matrices).
+   *
+   * The optional FEEvaluationBase object allows several FEEvaluation objects
+   * to share the geometry evaluation, i.e., the underlying mapping and
+   * quadrature points do only need to be evaluated once. This only works if
+   * the quadrature formulas are the same. Otherwise, a new evaluation object
+   * is created. Make sure to not pass an optional object around when you
+   * intend to use the FEEvaluation object in %parallel with another one
+   * because otherwise the intended sharing may create race conditions.
+   */
+  template <int n_components_other>
+  FEEvaluationBase (const Mapping<dim>       &mapping,
+                    const FiniteElement<dim> &fe,
+                    const Quadrature<1>      &quadrature,
+                    const UpdateFlags         update_flags,
+                    const unsigned int        first_selected_component,
+                    const FEEvaluationBase<dim,n_components_other,Number> *other);
+
+  /**
+   * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
+   * quadrature, and update flags, the underlying geometry evaluation based on
+   * FEValues will be deep-copied in order to allow for using in parallel with
+   * threads.
    */
   FEEvaluationBase (const FEEvaluationBase &other);
 
@@ -807,17 +866,19 @@ protected:
   bool gradients_quad_submitted;
 
   /**
-   * Geometry data generated by FEValues on the fly.
+   * Geometry data that can be generated FEValues on the fly with the
+   * respective constructor.
    */
-  SmartPointer<const MappingFEEvaluation<dim,Number> > mapped_geometry;
+  std_cxx1x::shared_ptr<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> > mapped_geometry;
 
   /**
-   * A pointer to the underlying DoFHandler.
+   * For use with on-the-fly evaluation, provide a data structure to store the
+   * global dof indices on the current cell from a reinit call.
    */
-  const DoFHandler<dim> *dof_handler;
+  std::vector<types::global_dof_index> old_style_dof_indices;
 
   /**
-   * For a DoFHandler with more than one finite element, select at which
+   * For a FiniteElement with more than one finite element, select at which
    * component this data structure should start.
    */
   const unsigned int first_selected_component;
@@ -827,6 +888,12 @@ protected:
    * MatrixFree object was given at initialization.
    */
   mutable std::vector<types::global_dof_index> local_dof_indices;
+
+  /**
+   * Make other FEEvaluationBase as well as FEEvaluation objects friends.
+   */
+  template <int, int, typename> friend class FEEvaluationBase;
+  template <int, int, int, int, typename> friend class FEEvaluation;
 };
 
 
@@ -865,34 +932,16 @@ protected:
                       const unsigned int            n_q_points);
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   * With this initialization, no call to a reinit method of this
-   * class. Instead, it is enough if the geometry is initialized to a given
-   * cell iterator. Moreover, beware that a kernel using this method does not
-   * vectorize over several elements (which is most efficient for vector
-   * operations), but only possibly within the element if the
-   * evaluate/integrate routines are combined (e.g. for matrix assembly).
-   */
-  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>                 &dof_handler,
-                      const unsigned int                     first_selected_component = 0);
+   * Constructor with reduced functionality for similar usage of FEEvaluation
+   * as FEValues, including matrix assembly.
+   */
+  template <int n_components_other>
+  FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other);
 
   /**
    * Copy constructor
@@ -1022,34 +1071,16 @@ protected:
                       const unsigned int            n_q_points);
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   * With this initialization, no call to a reinit method of this
-   * class. Instead, it is enough if the geometry is initialized to a given
-   * cell iterator. Moreover, beware that a kernel using this method does not
-   * vectorize over several elements (which is most efficient for vector
-   * operations), but only possibly within the element if the
-   * evaluate/integrate routines are combined (e.g. for matrix assembly).
-   */
-  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>                 &dof_handler,
-                      const unsigned int                     first_selected_component = 0);
+   * Constructor with reduced functionality for similar usage of FEEvaluation
+   * as FEValues, including matrix assembly.
+   */
+  template <int n_components_other>
+  FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other);
 
   /**
    * Copy constructor
@@ -1188,34 +1219,16 @@ protected:
                       const unsigned int            n_q_points);
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   * With this initialization, no call to a reinit method of this
-   * class. Instead, it is enough if the geometry is initialized to a given
-   * cell iterator. Moreover, beware that a kernel using this method does not
-   * vectorize over several elements (which is most efficient for vector
-   * operations), but only possibly within the element if the
-   * evaluate/integrate routines are combined (e.g. for matrix assembly).
-   */
-  FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>                 &dof_handler,
-                      const unsigned int                     first_selected_component = 0);
+   * Constructor with reduced functionality for similar usage of FEEvaluation
+   * as FEValues, including matrix assembly.
+   */
+  template <int n_components_other>
+  FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other);
 
   /**
    * Copy constructor
@@ -1343,34 +1356,16 @@ protected:
                       const unsigned int          n_q_points);
 
   /**
-   * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
-   *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   * With this initialization, no call to a reinit method of this
-   * class. Instead, it is enough if the geometry is initialized to a given
-   * cell iterator. Moreover, beware that a kernel using this method does not
-   * vectorize over several elements (which is most efficient for vector
-   * operations), but only possibly within the element if the
-   * evaluate/integrate routines are combined (e.g. for matrix assembly).
-   */
-  FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry,
-                      const DoFHandler<1>                 &dof_handler,
-                      const unsigned int                   first_selected_component = 0);
+   * Constructor with reduced functionality for similar usage of FEEvaluation
+   * as FEValues, including matrix assembly.
+   */
+  template <int n_components_other>
+  FEEvaluationAccess (const Mapping<1>       &mapping,
+                      const FiniteElement<1> &fe,
+                      const Quadrature<1>    &quadrature,
+                      const UpdateFlags       update_flags,
+                      const unsigned int      first_selected_component,
+                      const FEEvaluationBase<1,n_components_other,Number> *other);
 
   /**
    * Copy constructor
@@ -1387,26 +1382,125 @@ protected:
  * functions that make it much faster (between 5 and 500, depending on the
  * polynomial order).
  *
- * This class can be used in two different ways. The first way is to
- * initialize it from a MatrixFree object that caches everything related to
- * the degrees of freedom and the mapping information. This way, it is
- * 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
- * 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).
+ * <h3>Usage and initialization</h3>
+ *
+ * <h4>Fast usage in combination with MatrixFree</h4>
  *
- * This class contains specialized evaluation routines for several elements,
- * including standard FE_Q or FE_DGQ elements and quadrature points symmetric
- * around 0.5 (like Gauss quadrature), FE_DGP elements based on truncated
- * tensor products as well as the faster case of Gauss-Lobatto elements with
- * Gauss-Lobatto quadrature which give diagonal mass matrices and quicker
- * evaluation internally. Note that many of the operations available through
- * this class are inherited from the base class FEEvaluationBase, in
- * particular reading from and writing to vectors. Also, the class inherits
- * from FEEvaluationAccess that implements access to values, gradients and
- * Hessians of the finite element function on quadrature points.
+ * The first and foremost way of usage is to initialize this class from a
+ * MatrixFree object that caches everything related to the degrees of freedom
+ * and the mapping information. This way, it is possible to use vectorization
+ * for applying a vector operation for several cells at once. This setting is
+ * explained in the step-37 and step-48 tutorial programs. For vector-valued
+ * problems, the deal.II test suite includes a few additional examples as
+ * well, e.g. the Stokes operator found at
+ * https://github.com/dealii/dealii/blob/master/tests/matrix_free/matrix_vector_stokes_noflux.cc
+ *
+ * For most operator evaluation tasks, this path provides the most efficient
+ * solution by combining pre-computed data for the mapping (Jacobian
+ * transformations for the geometry description) with on-the-fly evaluation of
+ * basis functions. In other words, the framework provides a trade-off between
+ * memory usage and initialization of objects that is suitable for matrix-free
+ * operator evaluation.
+ *
+ * <h4>Usage without pre-initialized MatrixFree object</h4>
+ *
+ * The second form of usage is to initialize FEEvaluation from geometry
+ * information generated by FEValues. This allows to apply the integration
+ * loops on the fly without prior initialization of MatrixFree objects. This
+ * can be useful when the memory and initialization cost of MatrixFree is not
+ * acceptable, e.g. when a different number of quadrature points should be
+ * used for one single evaluation in error computation. Also, when using the
+ * routines of this class to assemble matrices the trade-off implied by the
+ * MatrixFree class may not be desired. In such a case, the cost to initialize
+ * the necessary geometry data on the fly is comparably low and thus avoiding
+ * a global object MatrixFree can be useful. When used in this way, reinit
+ * methods reminiscent from FEValues with a cell iterator are to be
+ * used. However, note that this model results in working on a single cell at
+ * a time, with geometry data duplicated in all components of the vectorized
+ * array. Thus, vectorization is only useful when it can apply the same
+ * operation on different data, e.g. when performing matrix assembly.
+ *
+ * As an example, consider the following code to assemble the contributions to
+ * the Laplace matrix:
+ *
+ * @code
+ * FEEvaluation<dim,fe_degree> fe_eval (mapping, finite_element,
+ *                                      QGauss<1>(fe_degree+1), flags);
+ * for (cell = dof_handler.begin_active();
+ *      cell != dof_handler.end();
+ *      ++cell)
+ *   {
+ *     fe_eval.reinit(cell);
+ *     for (unsigned int i=0; i<dofs_per_cell;
+ *          i += VectorizedArray<double>::n_array_elements)
+ *       {
+ *         const unsigned int n_items =
+ *           i+VectorizedArray<double>::n_array_elements > dofs_per_cell ?
+ *           (dofs_per_cell - i) : VectorizedArray<double>::n_array_elements;
+ *
+ *         // Set n_items unit vectors
+ *         for (unsigned int j=0; j<dofs_per_cell; ++j)
+ *           fe_eval.begin_dof_values()[j]  = VectorizedArray<double>();
+ *         for (unsigned int v=0; v<n_items; ++v)
+ *           fe_eval.begin_dof_values()[i+v][v] = 1.;
+ *
+ *         // Apply operator on unit vector
+ *         fe_eval.evaluate(true, true);
+ *         for (unsigned int q=0; q<n_q_points; ++q)
+ *           {
+ *             fe_eval.submit_value(10.*fe_eval.get_value(q), q);
+ *             fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+ *           }
+ *         fe_eval.integrate(true, true);
+ *
+ *         // Insert computed entries in matrix
+ *         for (unsigned int v=0; v<n_items; ++v)
+ *           for (unsigned int j=0; j<dofs_per_cell; ++j)
+ *             cell_matrix(fe_eval.get_internal_dof_numbering()[j],
+ *                         fe_eval.get_internal_dof_numbering()[i+v])
+ *               = fe_eval.begin_dof_values()[j][v];
+ *       }
+ *     ...
+ *   }
+ * @endcode
+ *
+ * This code generates the columns of the cell matrix with the loop over @p i
+ * above. The way this is done is the following: FEEvaluation's routines focus
+ * on the evaluation of finite element operators, so the way to get a cell
+ * matrix out of an operator evaluation is to apply it to all the unit vectors
+ * on the cell. This might seem inefficient but the evaluation routines used
+ * here are so quick that they still work much faster than what is possible
+ * with FEValues.
+ *
+ * Due to vectorization, we can actually generate matrix data for several unit
+ * vectors at a time (e.g. 4). The variable @p n_items make sure that we do
+ * the last iteration where the number of cell dofs is not divisible by the
+ * vectorization length correctly. Also note that we need to get the internal
+ * dof numbering applied by fe_eval because FEEvaluation internally uses a
+ * lexicographic numbering of degrees of freedom. This is necessary to
+ * efficiently work with tensor products where all degrees of freedom along a
+ * dimension must be laid out in a regular way.
+ *
+ * <h3>Description of evaluation routines</h3>
+ *
+ * This class contains specialized evaluation routines for several elements
+ * based on tensor-product quadrature formulas and tensor-product-like shape
+ * functions, including standard FE_Q or FE_DGQ elements and quadrature points
+ * symmetric around 0.5 (like Gauss quadrature), FE_DGP elements based on
+ * truncated tensor products as well as the faster case of Gauss-Lobatto
+ * elements with Gauss-Lobatto quadrature which give diagonal mass matrices
+ * and quicker evaluation internally. The main benefit of this class is the
+ * evaluation of all shape functions in all quadrature or integration over all
+ * shape functions in <code>dim (fe_degree+1)<sup>dim+1</sup> </code>
+ * operations instead of the slower <code>
+ * (fe_degree+1)<sup>2*dim</sup></code> complexity in the evaluation routines
+ * of FEValues.
+ *
+ * Note that many of the operations available through this class are inherited
+ * from the base class FEEvaluationBase, in particular reading from and
+ * writing to vectors. Also, the class inherits from FEEvaluationAccess that
+ * implements access to values, gradients and Hessians of the finite element
+ * function on quadrature points.
  *
  * This class assumes that shape functions of the FiniteElement under
  * consideration do <em>not</em> depend on the actual shape of the cells in
@@ -1421,20 +1515,21 @@ protected:
  *                  fe_degree+1 degrees of freedom per coordinate direction
  *
  * @param n_q_points_1d Number of points in the quadrature formula in 1D,
- *                   usually chosen as fe_degree+1
+ *                   defaults to fe_degree+1
  *
  * @param n_components Number of vector components when solving a system of
  *                  PDEs. If the same operation is applied to several
  *                  components of a PDE (e.g. a vector Laplace equation), they
  *                  can be applied simultaneously with one call (and often
- *                  more efficiently)
+ *                  more efficiently). Defaults to 1.
  *
- * @param Number Number format, usually @p double or @p float
+ * @param Number Number format, usually @p double or @p float.
+ *                  Defaults to @p double
  *
  * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
  */
-template <int dim, int fe_degree, int n_q_points_1d = fe_degree+1,
-          int n_components_ = 1, typename Number = double >
+template <int dim, int fe_degree, int n_q_points_1d, int n_components_,
+          typename Number >
 class FEEvaluation : public FEEvaluationAccess<dim,n_components_,Number>
 {
 public:
@@ -1457,34 +1552,69 @@ public:
                 const unsigned int            fe_no   = 0,
                 const unsigned int            quad_no = 0);
 
-  /**
-   * Copy constructor
-   */
-  FEEvaluation (const FEEvaluation &other);
-
   /**
    * Constructor that comes with reduced functionality and works similar as
-   * FEValues. The user has to provide a structure of type MappingFEEvaluation
-   * and a DoFHandler in order to allow for reading out the finite element
-   * data. It uses the data provided by dof_handler.get_fe(). If the element
-   * is vector-valued, the optional argument allows to specify the index of
-   * the base element (as long as the element is primitive, non-primitive are
-   * not supported currently).
+   * FEValues. The arguments are similar to the ones passed to the constructor
+   * of FEValues, with the notable difference that FEEvaluation expects a
+   * one-dimensional quadrature formula, Quadrature<1>, instead of a @p dim
+   * dimensional one. The finite element can be both scalar or vector valued,
+   * but this method always only selects a scalar base element at a time (with
+   * @p n_components copies as specified by the class template). For
+   * vector-valued elements, the optional argument @p first_selected_component
+   * allows to specify the index of the base element to be used for
+   * evaluation. Note that the internal data structures always assume that the
+   * base element is primitive, non-primitive are not supported currently.
    *
-   * With initialization from a FEValues object, no call to a reinit method of
-   * this class is necessary. Instead, it is enough if the geometry is
-   * initialized to a given cell iterator. It can also read from or write to
-   * vectors in the standard way for DoFHandler<dim>::active_cell_iterator
-   * types (which is less efficient with MPI since index translation has to be
-   * done), but of course only for one cell at a time. Hence, a kernel using
-   * this method does not vectorize over several elements (which is most
-   * efficient for vector operations), but only possibly within the element if
-   * the evaluate/integrate routines are combined (e.g. for computing cell
-   * matrices).
-   */
-  FEEvaluation (const MappingFEEvaluation<dim,Number> &geometry,
-                const DoFHandler<dim>                 &dof_handler,
-                const unsigned int                     first_selected_component = 0);
+   * As known from FEValues, a call to the reinit method with a
+   * Triangulation<dim>::cell_iterator is necessary to make the geometry and
+   * degrees of freedom of the current class known. If the iterator includes
+   * DoFHandler information (i.e., it is a DoFHandler<dim>::cell_iterator or
+   * similar), the initialization allows to also read from or write to vectors
+   * in the standard way for DoFHandler<dim>::active_cell_iterator types for
+   * one cell at a time. However, this approach is much slower than the path
+   * with MatrixFree with MPI since index translation has to be done. As only
+   * one cell at a time is used, this method does not vectorize over several
+   * elements (which is most efficient for vector operations), but only
+   * possibly within the element if the evaluate/integrate routines are
+   * combined inside user code (e.g. for computing cell matrices).
+   */
+  FEEvaluation (const Mapping<dim>       &mapping,
+                const FiniteElement<dim> &fe,
+                const Quadrature<1>      &quadrature,
+                const UpdateFlags         update_flags,
+                const unsigned int        first_selected_component = 0);
+
+  /**
+   * Constructor for the reduced functionality. Similar as the other
+   * constructor but uses MappingQ1 implicitly.
+   */
+  FEEvaluation (const FiniteElement<dim> &fe,
+                const Quadrature<1>      &quadrature,
+                const UpdateFlags         update_flags,
+                const unsigned int        first_selected_component = 0);
+
+  /**
+   * Constructor for the reduced functionality. Similar to the other
+   * constructor with FiniteElement argument but using another
+   * FEEvaluationBase object to provide info about the geometry. This allows
+   * several FEEvaluation objects to share the geometry evaluation, i.e., the
+   * underlying mapping and quadrature points do only need to be evaluated
+   * once. Make sure to not pass an optional object around when you intend to
+   * use the FEEvaluation object in %parallel to the given one because
+   * otherwise the intended sharing may create race conditions.
+   */
+  template <int n_components_other>
+  FEEvaluation (const FiniteElement<dim> &fe,
+                const FEEvaluationBase<dim,n_components_other,Number> &other,
+                const unsigned int        first_selected_component = 0);
+
+  /**
+   * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
+   * quadrature, and update flags, the underlying geometry evaluation based on
+   * FEValues will be deep-copied in order to allow for using in parallel with
+   * threads.
+   */
+  FEEvaluation (const FEEvaluation &other);
 
   /**
    * Evaluates the function values, the gradients, and the Laplacians of the
@@ -1590,16 +1720,6 @@ public:
 :
   BaseClass (matrix_free, fe_no, quad_no)
   {}
-
-  /**
-   * Constructor.
-   */
-  FEEvaluationGeneral (const MappingFEEvaluation<dim,Number> &geometry,
-                       const DoFHandler<dim>                 &dof_handler,
-                       const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
-:
-  BaseClass (geometry, dof_handler, first_selected_component)
-  {}
 };
 
 
@@ -1631,16 +1751,6 @@ public:
 :
   BaseClass (matrix_free, fe_no, quad_no)
   {}
-
-  /**
-   * Constructor.
-   */
-  FEEvaluationGL (const MappingFEEvaluation<dim,Number> &geometry,
-                  const DoFHandler<dim>                 &dof_handler,
-                  const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
-:
-  BaseClass (geometry, dof_handler, first_selected_component)
-  {}
 };
 
 
@@ -1698,16 +1808,6 @@ public:
 :
   BaseClass (matrix_free, fe_no, quad_no)
   {}
-
-  /**
-   * Constructor.
-   */
-  FEEvaluationDGP (const MappingFEEvaluation<dim,Number> &geometry,
-                   const DoFHandler<dim>                 &dof_handler,
-                   const unsigned int                     first_selected_component = 0) DEAL_II_DEPRECATED
-:
-  BaseClass (geometry, dof_handler, first_selected_component)
-  {}
 };
 
 
@@ -1753,7 +1853,6 @@ FEEvaluationBase<dim,n_components_,Number>
   cell               (numbers::invalid_unsigned_int),
   cell_type          (internal::MatrixFreeFunctions::undefined),
   cell_data_number   (0),
-  dof_handler        (0),
   first_selected_component (0)
 {
   for (unsigned int c=0; c<n_components_; ++c)
@@ -1788,11 +1887,15 @@ FEEvaluationBase<dim,n_components_,Number>
 
 
 template <int dim, int n_components_, typename Number>
+template <int n_components_other>
 inline
 FEEvaluationBase<dim,n_components_,Number>
-::FEEvaluationBase (const MappingFEEvaluation<dim,Number> &geometry,
-                    const DoFHandler<dim>                 &dof_handler_in,
-                    const unsigned int                     first_selected_component)
+::FEEvaluationBase (const Mapping<dim>       &mapping,
+                    const FiniteElement<dim> &fe,
+                    const Quadrature<1>      &quadrature,
+                    const UpdateFlags         update_flags,
+                    const unsigned int        first_selected_component,
+                    const FEEvaluationBase<dim,n_components_other,Number> *other)
   :
   quad_no            (-1),
   n_fe_components    (n_components_),
@@ -1801,24 +1904,25 @@ FEEvaluationBase<dim,n_components_,Number>
   matrix_info        (0),
   dof_info           (0),
   mapping_info       (0),
-  stored_shape_info  (new internal::MatrixFreeFunctions::ShapeInfo<Number>(geometry.get_quadrature(), dof_handler_in.get_fe(), dof_handler_in.get_fe().component_to_base_index(first_selected_component).first)),
+  // select the correct base element from the given FE component
+  stored_shape_info  (new internal::MatrixFreeFunctions::ShapeInfo<Number>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
   data               (stored_shape_info.get()),
   cartesian_data     (0),
-  jacobian           (geometry.get_inverse_jacobians().begin()),
-  J_value            (geometry.get_JxW_values().begin()),
+  jacobian           (0),
+  J_value            (0),
   quadrature_weights (0),
-  quadrature_points  (geometry.get_quadrature_points().begin()),
+  quadrature_points  (0),
   jacobian_grad      (0),
   jacobian_grad_upper(0),
   cell               (0),
   cell_type          (internal::MatrixFreeFunctions::general),
   cell_data_number   (0),
-  mapped_geometry    (&geometry),
-  dof_handler        (&dof_handler_in),
-  first_selected_component (first_selected_component)
+  // keep the number of the selected component within the current base element
+  // for reading dof values
+  first_selected_component (fe.component_to_base_index(first_selected_component).second)
 {
   const unsigned int base_element_number =
-    dof_handler_in.get_fe().component_to_base_index(first_selected_component).first;
+    fe.component_to_base_index(first_selected_component).first;
   for (unsigned int c=0; c<n_components_; ++c)
     {
       values_dofs[c] = 0;
@@ -1828,8 +1932,21 @@ 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_number) == 1 ||
-         dof_handler->get_fe().element_multiplicity(base_element_number)-first_selected_component >= n_components_,
+
+  Assert(other == 0 || other->mapped_geometry.get() != 0, ExcInternalError());
+  if (other != 0 &&
+      other->mapped_geometry->get_quadrature() == quadrature)
+    mapped_geometry = other->mapped_geometry;
+  else
+    mapped_geometry.reset(new internal::MatrixFreeFunctions::
+                          MappingDataOnTheFly<dim,Number>(mapping, quadrature,
+                                                          update_flags));
+  jacobian = mapped_geometry->get_inverse_jacobians().begin();
+  J_value = mapped_geometry->get_JxW_values().begin();
+  quadrature_points = mapped_geometry->get_quadrature_points().begin();
+
+  Assert(fe.element_multiplicity(base_element_number) == 1 ||
+         fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_,
          ExcMessage("The underlying element must at least contain as many "
                     "components as requested by this class"));
 }
@@ -1860,8 +1977,6 @@ FEEvaluationBase<dim,n_components_,Number>
   cell               (other.cell),
   cell_type          (other.cell_type),
   cell_data_number   (other.cell_data_number),
-  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)
@@ -1873,6 +1988,19 @@ FEEvaluationBase<dim,n_components_,Number>
       for (unsigned int d=0; d<(dim*dim+dim)/2; ++d)
         hessians_quad[c][d] = 0;
     }
+
+  // Create deep copy of mapped geometry for use in parallel...
+  if (other.mapped_geometry.get() != 0)
+    {
+      mapped_geometry.reset
+      (new internal::MatrixFreeFunctions::
+       MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
+                                       other.mapped_geometry->get_quadrature(),
+                                       other.mapped_geometry->get_fe_values().get_update_flags()));
+      jacobian = mapped_geometry->get_inverse_jacobians().begin();
+      J_value = mapped_geometry->get_JxW_values().begin();
+      quadrature_points = mapped_geometry->get_quadrature_points().begin();
+    }
 }
 
 
@@ -1882,7 +2010,9 @@ inline
 void
 FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
 {
-  Assert (mapped_geometry == 0, ExcMessage("FEEvaluation was initialized without a matrix-free object. Integer indexing is not possible"));
+  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());
@@ -1955,6 +2085,46 @@ FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
 
 
 
+template <int dim, int n_components_, typename Number>
+template <typename DH, bool level_dof_access>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::reinit (const TriaIterator<DoFCellAccessor<DH,level_dof_access> > &cell)
+{
+  Assert(matrix_info == 0,
+         ExcMessage("Cannot use initialization from cell iterator if "
+                    "initialized from MatrixFree object. Use variant for "
+                    "on the fly computation with arguments as for FEValues "
+                    "instead"));
+  Assert(mapped_geometry.get() != 0, ExcNotInitialized());
+  mapped_geometry->reinit(static_cast<typename Triangulation<dim>::cell_iterator>(cell));
+  local_dof_indices.resize(cell->get_fe().dofs_per_cell);
+  if (level_dof_access)
+    cell->get_mg_dof_indices(local_dof_indices);
+  else
+    cell->get_dof_indices(local_dof_indices);
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+void
+FEEvaluationBase<dim,n_components_,Number>
+::reinit (const typename Triangulation<dim>::cell_iterator &cell)
+{
+  Assert(matrix_info == 0,
+         ExcMessage("Cannot use initialization from cell iterator if "
+                    "initialized from MatrixFree object. Use variant for "
+                    "on the fly computation with arguments as for FEValues "
+                    "instead"));
+  Assert(mapped_geometry.get() != 0, ExcNotInitialized());
+  mapped_geometry->reinit(cell);
+}
+
+
+
 template <int dim, int n_components_, typename Number>
 inline
 unsigned int
@@ -2278,13 +2448,7 @@ FEEvaluationBase<dim,n_components_,Number>
   // 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);
+      Assert (!local_dof_indices.empty(), ExcNotInitialized());
 
       unsigned int index = first_selected_component * this->data->dofs_per_cell;
       for (unsigned int comp = 0; comp<n_components; ++comp)
@@ -3030,6 +3194,17 @@ FEEvaluationBase<dim,n_components_,Number>
 
 /*------------------------------ access to data fields ----------------------*/
 
+template <int dim, int n_components, typename Number>
+inline
+const std::vector<unsigned int> &
+FEEvaluationBase<dim,n_components,Number>::
+get_internal_dof_numbering() const
+{
+  return data->lexicographic_numbering;
+}
+
+
+
 
 template <int dim, int n_components, typename Number>
 inline
@@ -3620,13 +3795,18 @@ FEEvaluationAccess<dim,n_components_,Number>
 
 
 template <int dim, int n_components_, typename Number>
+template <int n_components_other>
 inline
 FEEvaluationAccess<dim,n_components_,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   first_selected_component)
+::FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other)
   :
-  FEEvaluationBase <dim,n_components_,Number> (geometry, dof_handler, first_selected_component)
+  FEEvaluationBase <dim,n_components_,Number>(mapping, fe, quadrature, update_flags,
+                                              first_selected_component, other)
 {}
 
 
@@ -3661,13 +3841,18 @@ FEEvaluationAccess<dim,1,Number>
 
 
 template <int dim, typename Number>
+template <int n_components_other>
 inline
 FEEvaluationAccess<dim,1,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   first_selected_component)
+::FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other)
   :
-  FEEvaluationBase <dim,1,Number> (geometry, dof_handler, first_selected_component)
+  FEEvaluationBase <dim,1,Number> (mapping, fe, quadrature, update_flags,
+                                   first_selected_component, other)
 {}
 
 
@@ -3895,13 +4080,18 @@ FEEvaluationAccess<dim,dim,Number>
 
 
 template <int dim, typename Number>
+template <int n_components_other>
 inline
 FEEvaluationAccess<dim,dim,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<dim,Number> &geometry,
-                      const DoFHandler<dim>               &dof_handler,
-                      const unsigned int                   first_selected_component)
+::FEEvaluationAccess (const Mapping<dim>       &mapping,
+                      const FiniteElement<dim> &fe,
+                      const Quadrature<1>      &quadrature,
+                      const UpdateFlags         update_flags,
+                      const unsigned int        first_selected_component,
+                      const FEEvaluationBase<dim,n_components_other,Number> *other)
   :
-  FEEvaluationBase <dim,dim,Number> (geometry, dof_handler, first_selected_component)
+  FEEvaluationBase <dim,dim,Number> (mapping, fe, quadrature, update_flags,
+                                     first_selected_component, other)
 {}
 
 
@@ -4250,13 +4440,18 @@ FEEvaluationAccess<1,1,Number>
 
 
 template <typename Number>
+template <int n_components_other>
 inline
 FEEvaluationAccess<1,1,Number>
-::FEEvaluationAccess (const MappingFEEvaluation<1,Number> &geometry,
-                      const DoFHandler<1>               &dof_handler,
-                      const unsigned int                 first_selected_component)
+::FEEvaluationAccess (const Mapping<1>       &mapping,
+                      const FiniteElement<1> &fe,
+                      const Quadrature<1>    &quadrature,
+                      const UpdateFlags       update_flags,
+                      const unsigned int      first_selected_component,
+                      const FEEvaluationBase<1,n_components_other,Number> *other)
   :
-  FEEvaluationBase <1,1,Number> (geometry, dof_handler, first_selected_component)
+  FEEvaluationBase <1,1,Number> (mapping, fe, quadrature, update_flags,
+                                 first_selected_component, other)
 {}
 
 
@@ -6255,11 +6450,56 @@ 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 MappingFEEvaluation<dim,Number> &geometry,
-                const DoFHandler<dim>               &dof_handler,
-                const unsigned int                   first_selected_component)
+::FEEvaluation (const Mapping<dim>       &mapping,
+                const FiniteElement<dim> &fe,
+                const Quadrature<1>      &quadrature,
+                const UpdateFlags         update_flags,
+                const unsigned int        first_selected_component)
   :
-  BaseClass (geometry, dof_handler, first_selected_component),
+  BaseClass (mapping, fe, quadrature, update_flags,
+             first_selected_component,
+             static_cast<FEEvaluationBase<dim,1,Number>*>(0)),
+  dofs_per_cell (this->data->dofs_per_cell)
+{
+  check_template_arguments(numbers::invalid_unsigned_int);
+  set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FiniteElement<dim> &fe,
+                const Quadrature<1>      &quadrature,
+                const UpdateFlags         update_flags,
+                const unsigned int        first_selected_component)
+  :
+  BaseClass (StaticMappingQ1<dim>::mapping, fe, quadrature, update_flags,
+             first_selected_component,
+             static_cast<FEEvaluationBase<dim,1,Number>*>(0)),
+  dofs_per_cell (this->data->dofs_per_cell)
+{
+  check_template_arguments(numbers::invalid_unsigned_int);
+  set_data_pointers();
+}
+
+
+
+template <int dim, int fe_degree,  int n_q_points_1d, int n_components_,
+          typename Number>
+template <int n_components_other>
+inline
+FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
+::FEEvaluation (const FiniteElement<dim> &fe,
+                const FEEvaluationBase<dim,n_components_other,Number> &other,
+                const unsigned int        first_selected_component)
+  :
+  BaseClass (other.mapped_geometry->get_fe_values().get_mapping(),
+             fe, other.mapped_geometry->get_quadrature(),
+             other.mapped_geometry->get_fe_values().get_update_flags(),
+             first_selected_component, &other),
   dofs_per_cell (this->data->dofs_per_cell)
 {
   check_template_arguments(numbers::invalid_unsigned_int);
@@ -6395,9 +6635,13 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       message += Utilities::int_to_string(fe_degree) + ",";
       message += Utilities::int_to_string(n_q_points_1d);
       message += "," + Utilities::int_to_string(n_components);
-      message += ",Number>(data, ";
-      message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(this->quad_no) + ")\n";
+      message += ",Number>(data";
+      if (fe_no != numbers::invalid_unsigned_int)
+        {
+          message += ", " + Utilities::int_to_string(fe_no) + ", ";
+          message += Utilities::int_to_string(this->quad_no);
+        }
+      message += ")\n";
 
       // check whether some other vector component has the correct number of
       // points
@@ -6438,9 +6682,13 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
           message += Utilities::int_to_string(fe_degree) + ",";
           message += Utilities::int_to_string(n_q_points_1d);
           message += "," + Utilities::int_to_string(n_components);
-          message += ",Number>(data, ";
-          message += Utilities::int_to_string(proposed_dof_comp) + ", ";
-          message += Utilities::int_to_string(proposed_quad_comp) + ")?\n";
+          message += ",Number>(data";
+          if (fe_no != numbers::invalid_unsigned_int)
+            {
+              message += ", " + Utilities::int_to_string(proposed_dof_comp) + ", ";
+              message += Utilities::int_to_string(proposed_quad_comp);
+            }
+          message += ")?\n";
           std::string correct_pos;
           if (proposed_dof_comp != fe_no)
             correct_pos = " ^ ";
@@ -6460,9 +6708,13 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
       message += Utilities::int_to_string(this->data->fe_degree) + ",";
       message += Utilities::int_to_string(proposed_n_q_points_1d);
       message += "," + Utilities::int_to_string(n_components);
-      message += ",Number>(data, ";
-      message += Utilities::int_to_string(fe_no) + ", ";
-      message += Utilities::int_to_string(this->quad_no) + ")?\n";
+      message += ",Number>(data";
+      if (fe_no != numbers::invalid_unsigned_int)
+        {
+          message += ", " + Utilities::int_to_string(fe_no) + ", ";
+          message += Utilities::int_to_string(this->quad_no);
+        }
+      message += ")?\n";
       std::string correct_pos;
       if (this->data->fe_degree != fe_degree)
         correct_pos = " ^";
@@ -6544,6 +6796,8 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
 {
   Assert (this->dof_values_initialized == true,
           internal::ExcAccessToUninitializedField());
+  Assert(this->matrix_info != 0 ||
+         this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
   // Select algorithm matching the element type at run time (the function
   // pointer is easy to predict, so negligible in cost)
@@ -6576,6 +6830,8 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
   if (integrate_grad == true)
     Assert (this->gradients_quad_submitted == true,
             internal::ExcAccessToUninitializedField());
+  Assert(this->matrix_info != 0 ||
+         this->mapped_geometry->is_initialized(), ExcNotInitialized());
 
   // Select algorithm matching the element type at run time (the function
   // pointer is easy to predict, so negligible in cost)
diff --git a/include/deal.II/matrix_free/mapping_data_on_the_fly.h b/include/deal.II/matrix_free/mapping_data_on_the_fly.h
new file mode 100644 (file)
index 0000000..678b971
--- /dev/null
@@ -0,0 +1,337 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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_mapping_data_on_the_fly_h
+#define __deal2__matrix_free_mapping_data_on_the_fly_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/subscriptor.h>
+#include <deal.II/base/vectorization.h>
+#include <deal.II/base/aligned_vector.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
+
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    /**
+     * 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.
+     *
+     * As opposed to the Mapping classes in deal.II, this class does not
+     * actually provide a boundary description that can be used to evaluate
+     * the geometry, but it rather provides the evaluated geometry from a
+     * given deal.II mapping (as passed to the constructor of this class) in a
+     * form accessible to FEEvaluation.
+     *
+     * @author Martin Kronbichler, 2014
+     */
+    template <int dim, typename Number=double>
+    class MappingDataOnTheFly
+    {
+    public:
+      /**
+       * 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.
+       */
+      MappingDataOnTheFly (const Mapping<dim> &mapping,
+                           const Quadrature<1> &quadrature,
+                           const UpdateFlags update_flags);
+
+      /**
+       * Constructor. Instead of providing a mapping, use MappingQ1.
+       */
+      MappingDataOnTheFly (const Quadrature<1> &quadrature,
+                           const UpdateFlags update_flags);
+
+      /**
+       * Initialize with the given cell iterator.
+       */
+      void reinit(typename dealii::Triangulation<dim>::cell_iterator cell);
+
+      /**
+       * Returns whether reinit() has been called at least once, i.e., a cell
+       * has been set.
+       */
+      bool is_initialized() const;
+
+      /**
+       * Return a triangulation iterator to the current cell.
+       */
+      typename dealii::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
+       * FEEvaluation, it returns tensors of vectorized arrays, even though all
+       * components are equal.
+       */
+      const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
+      get_inverse_jacobians() const;
+
+      /**
+       * Return a vector of quadrature weights times the Jacobian determinant
+       * (JxW). For compatibility with FEEvaluation, it returns tensors of
+       * vectorized arrays, even though all components are equal.
+       */
+      const AlignedVector<VectorizedArray<Number> > &
+      get_JxW_values() const;
+
+      /**
+       * Return a vector of quadrature points in real space on the given
+       * cell. For compatibility with FEEvaluation, it returns tensors of
+       * vectorized arrays, even though all components are equal.
+       */
+      const AlignedVector<Point<dim,VectorizedArray<Number> > > &
+      get_quadrature_points() const;
+
+      /**
+       * Return a vector of quadrature points in real space on the given
+       * cell. For compatibility with FEEvaluation, it returns tensors of
+       * vectorized arrays, even though all components are equal.
+       */
+      const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
+      get_normal_vectors() const;
+
+      /**
+       * Return a reference to 1D quadrature underlying this object.
+       */
+      const Quadrature<1> &
+      get_quadrature () const;
+
+    private:
+      /**
+       * A cell iterator in case we generate the data on the fly to be able to
+       * check if we need to re-generate the information stored in this class.
+       */
+      typename dealii::Triangulation<dim>::cell_iterator present_cell;
+
+      /**
+       * Dummy finite element object necessary for initializing the FEValues
+       * object.
+       */
+      FE_Nothing<dim> fe_dummy;
+
+      /**
+       * An underlying FEValues object that performs the (scalar) evaluation.
+       */
+      FEValues<dim> fe_values;
+
+      /**
+       * Get 1D quadrature formula to be used for reinitializing shape info.
+       */
+      const Quadrature<1> quadrature_1d;
+
+      /**
+       * Inverse Jacobians, stored in vectorized array form.
+       */
+      AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > inverse_jacobians;
+
+      /**
+       * Stored Jacobian determinants and quadrature weights
+       */
+      AlignedVector<VectorizedArray<Number> > jxw_values;
+
+      /**
+       * Stored quadrature points
+       */
+      AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
+
+      /**
+       * Stored normal vectors (for face integration)
+       */
+      AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
+    };
+
+
+    /*----------------------- Inline functions ----------------------------------*/
+
+    template <int dim, typename Number>
+    inline
+    MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Mapping<dim> &mapping,
+                                                          const Quadrature<1> &quadrature,
+                                                          const UpdateFlags update_flags)
+      :
+      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()),
+      jxw_values(fe_values.get_quadrature().size()),
+      quadrature_points(fe_values.get_quadrature().size()),
+      normal_vectors(fe_values.get_quadrature().size())
+    {
+      Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
+             ExcNotImplemented());
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Quadrature<1> &quadrature,
+                                                          const UpdateFlags update_flags)
+      :
+      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()),
+      jxw_values(fe_values.get_quadrature().size()),
+      quadrature_points(fe_values.get_quadrature().size()),
+      normal_vectors(fe_values.get_quadrature().size())
+    {
+      Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
+             ExcNotImplemented());
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    void
+    MappingDataOnTheFly<dim,Number>::reinit(typename dealii::Triangulation<dim>::cell_iterator cell)
+    {
+      if (present_cell == cell)
+        return;
+      present_cell = cell;
+      fe_values.reinit(present_cell);
+      for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
+        {
+          if (fe_values.get_update_flags() & update_inverse_jacobians)
+            for (unsigned int d=0; d<dim; ++d)
+              for (unsigned int e=0; e<dim; ++e)
+                inverse_jacobians[q][d][e] = fe_values.inverse_jacobian(q)[e][d];
+          if (fe_values.get_update_flags() & update_quadrature_points)
+            for (unsigned int d=0; d<dim; ++d)
+              quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
+          if (fe_values.get_update_flags() & update_normal_vectors)
+            for (unsigned int d=0; d<dim; ++d)
+              normal_vectors[q][d] = fe_values.normal_vector(q)[d];
+          if (fe_values.get_update_flags() & update_JxW_values)
+            jxw_values[q] = fe_values.JxW(q);
+        }
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    bool
+    MappingDataOnTheFly<dim,Number>::is_initialized() const
+    {
+      return present_cell != typename dealii::Triangulation<dim>::cell_iterator();
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    typename dealii::Triangulation<dim>::cell_iterator
+    MappingDataOnTheFly<dim,Number>::get_cell() const
+    {
+      return fe_values.get_cell();
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const FEValues<dim> &
+    MappingDataOnTheFly<dim,Number>::get_fe_values() const
+    {
+      return fe_values;
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
+    MappingDataOnTheFly<dim,Number>::get_inverse_jacobians() const
+    {
+      return inverse_jacobians;
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
+    MappingDataOnTheFly<dim,Number>::get_normal_vectors() const
+    {
+      return normal_vectors;
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const AlignedVector<Point<dim,VectorizedArray<Number> > > &
+    MappingDataOnTheFly<dim,Number>::get_quadrature_points() const
+    {
+      return quadrature_points;
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const AlignedVector<VectorizedArray<Number> > &
+    MappingDataOnTheFly<dim,Number>::get_JxW_values() const
+    {
+      return jxw_values;
+    }
+
+
+
+    template <int dim, typename Number>
+    inline
+    const Quadrature<1> &
+    MappingDataOnTheFly<dim,Number>::get_quadrature() const
+    {
+      return quadrature_1d;
+    }
+
+
+  } // end of namespace MatrixFreeFunctions
+} // end of namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/matrix_free/mapping_fe_evaluation.h b/include/deal.II/matrix_free/mapping_fe_evaluation.h
deleted file mode 100644 (file)
index 626fc78..0000000
+++ /dev/null
@@ -1,313 +0,0 @@
-// ---------------------------------------------------------------------
-//
-// 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_mapping_fe_evaluation_h
-#define __deal2__matrix_free_mapping_fe_evaluation_h
-
-
-#include <deal.II/base/config.h>
-#include <deal.II/base/exceptions.h>
-#include <deal.II/base/subscriptor.h>
-#include <deal.II/base/vectorization.h>
-#include <deal.II/base/aligned_vector.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
-
-
-/**
- * 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.
- *
- * As opposed to the standard Mapping classes in deal.II, this class does not
- * actually provide a boundary description that can be used to evaluate the
- * geometry, but it rather provides the evaluated geometry from a given
- * deal.II mapping (as passed to the constructor of this class) in a form
- * accessible to FEEvaluation.
- *
- * @author Martin Kronbichler, 2014
- */
-template <int dim, typename Number=double>
-class MappingFEEvaluation : public Subscriptor
-{
-public:
-  /**
-   * 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.
-   */
-  MappingFEEvaluation (const Mapping<dim> &mapping,
-                       const Quadrature<1> &quadrature,
-                       const UpdateFlags update_flags);
-
-  /**
-   * Constructor. Instead of providing a mapping, use MappingQ1.
-   */
-  MappingFEEvaluation (const Quadrature<1> &quadrature,
-                       const UpdateFlags update_flags);
-
-  /**
-   * 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);
-
-  /**
-   * 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
-   * FEEvaluation, it returns tensors of vectorized arrays, even though all
-   * components are equal.
-   */
-  const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
-  get_inverse_jacobians() const;
-
-  /**
-   * Return a vector of quadrature weights times the Jacobian determinant
-   * (JxW). For compatibility with FEEvaluation, it returns tensors of
-   * vectorized arrays, even though all components are equal.
-   */
-  const AlignedVector<VectorizedArray<Number> > &
-  get_JxW_values() const;
-
-  /**
-   * Return a vector of quadrature points in real space on the given
-   * cell. For compatibility with FEEvaluation, it returns tensors of
-   * vectorized arrays, even though all components are equal.
-   */
-  const AlignedVector<Point<dim,VectorizedArray<Number> > > &
-  get_quadrature_points() const;
-
-  /**
-   * Return a vector of quadrature points in real space on the given
-   * cell. For compatibility with FEEvaluation, it returns tensors of
-   * vectorized arrays, even though all components are equal.
-   */
-  const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
-  get_normal_vectors() const;
-
-  /**
-   * Return a reference to 1D quadrature underlying this object.
-   */
-  const Quadrature<1> &
-  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.
-   */
-  FEValues<dim> fe_values;
-
-  /**
-   * Get 1D quadrature formula to be used for reinitializing shape info.
-   */
-  const Quadrature<1> quadrature_1d;
-
-  /**
-   * Inverse Jacobians, stored in vectorized array form.
-   */
-  AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > inverse_jacobians;
-
-  /**
-   * Stored Jacobian determinants and quadrature weights
-   */
-  AlignedVector<VectorizedArray<Number> > jxw_values;
-
-  /**
-   * Stored quadrature points
-   */
-  AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
-
-  /**
-   * Stored normal vectors (for face integration)
-   */
-  AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
-};
-
-
-/*----------------------- Inline functions ----------------------------------*/
-
-template <int dim, typename Number>
-inline
-MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Mapping<dim> &mapping,
-                                                      const Quadrature<1> &quadrature,
-                                                      const UpdateFlags update_flags)
-  :
-  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()),
-  jxw_values(fe_values.get_quadrature().size()),
-  quadrature_points(fe_values.get_quadrature().size()),
-  normal_vectors(fe_values.get_quadrature().size())
-{
-  Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
-         ExcNotImplemented());
-}
-
-
-
-template <int dim, typename Number>
-inline
-MappingFEEvaluation<dim,Number>::MappingFEEvaluation (const Quadrature<1> &quadrature,
-                                                      const UpdateFlags update_flags)
-  :
-  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()),
-  jxw_values(fe_values.get_quadrature().size()),
-  quadrature_points(fe_values.get_quadrature().size()),
-  normal_vectors(fe_values.get_quadrature().size())
-{
-  Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
-         ExcNotImplemented());
-}
-
-
-
-template <int dim, typename Number>
-template <typename ITERATOR>
-inline
-void
-MappingFEEvaluation<dim,Number>::reinit(ITERATOR 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)
-        for (unsigned int d=0; d<dim; ++d)
-          for (unsigned int e=0; e<dim; ++e)
-            inverse_jacobians[q][d][e] = fe_values.inverse_jacobian(q)[e][d];
-      if (fe_values.get_update_flags() & update_quadrature_points)
-        for (unsigned int d=0; d<dim; ++d)
-          quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
-      if (fe_values.get_update_flags() & update_normal_vectors)
-        for (unsigned int d=0; d<dim; ++d)
-          normal_vectors[q][d] = fe_values.normal_vector(q)[d];
-      if (fe_values.get_update_flags() & update_JxW_values)
-        jxw_values[q] = fe_values.JxW(q);
-    }
-}
-
-
-
-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> > > &
-MappingFEEvaluation<dim,Number>::get_inverse_jacobians() const
-{
-  return inverse_jacobians;
-}
-
-
-
-template <int dim, typename Number>
-inline
-const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
-MappingFEEvaluation<dim,Number>::get_normal_vectors() const
-{
-  return normal_vectors;
-}
-
-
-
-template <int dim, typename Number>
-inline
-const AlignedVector<Point<dim,VectorizedArray<Number> > > &
-MappingFEEvaluation<dim,Number>::get_quadrature_points() const
-{
-  return quadrature_points;
-}
-
-
-
-template <int dim, typename Number>
-inline
-const AlignedVector<VectorizedArray<Number> > &
-MappingFEEvaluation<dim,Number>::get_JxW_values() const
-{
-  return jxw_values;
-}
-
-
-
-template <int dim, typename Number>
-inline
-const Quadrature<1> &
-MappingFEEvaluation<dim,Number>::get_quadrature() const
-{
-  return quadrature_1d;
-}
-
-
-#ifndef DOXYGEN
-
-
-#endif  // ifndef DOXYGEN
-
-
-DEAL_II_NAMESPACE_CLOSE
-
-#endif
index bcbaac2c3fc9e5ee4c7c45e3bedb5e272a99e5a1..8b3aa6fa8171ca0b58ab62109a067b023bd083ed 100644 (file)
@@ -116,15 +116,27 @@ namespace internal
             // 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)
+            std::vector<unsigned int> lexicographic(fe_in.dofs_per_cell,
+                                                    numbers::invalid_unsigned_int);
+            unsigned int components_before = 0;
+            for (unsigned int e=0; e<base_element_number; ++e)
+              components_before += fe_in.element_multiplicity(e);
+            for (unsigned int comp=0;
+                 comp<fe_in.element_multiplicity(base_element_number); ++comp)
               for (unsigned int i=0; i<scalar_inv.size(); ++i)
-                lexicographic[fe_in.component_to_system_index(comp,i)]
+                lexicographic[fe_in.component_to_system_index(comp+components_before,i)]
                   = scalar_inv.size () * comp + scalar_inv[i];
 
-            // invert numbering again
-            lexicographic_numbering =
-              Utilities::invert_permutation(lexicographic);
+            // invert numbering again. Need to do it manually because we might
+            // have undefined blocks
+            lexicographic_numbering.resize(fe_in.element_multiplicity(base_element_number)*fe->dofs_per_cell);
+            for (unsigned int i=0; i<lexicographic.size(); ++i)
+              if (lexicographic[i] != numbers::invalid_unsigned_int)
+                {
+                  AssertIndexRange(lexicographic[i],
+                                   lexicographic_numbering.size());
+                  lexicographic_numbering[lexicographic[i]] = i;
+                }
           }
 
         // to evaluate 1D polynomials, evaluate along the line where y=z=0,
diff --git a/tests/matrix_free/assemble_matrix_01.cc b/tests/matrix_free/assemble_matrix_01.cc
new file mode 100644 (file)
index 0000000..fb862d1
--- /dev/null
@@ -0,0 +1,162 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test FEEvaluation for assembling the Laplace matrix. It is enough to just
+// consider the resulting element matrices element by element
+
+#include "../tests.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/lac/constraint_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+
+std::ofstream logfile("output");
+
+
+
+template <int dim, int fe_degree>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+
+  MappingQ<dim> mapping(fe_degree+1);
+
+  // compute matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  {
+    QGauss<dim>  quadrature_formula(fe_degree+1);
+
+    FEValues<dim> fe_values (mapping, dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    FEEvaluation<dim,fe_degree,fe_degree+1>
+      fe_eval (mapping, dof.get_fe(), QGauss<1>(fe_degree+1),
+               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);
+    FullMatrix<double>   test_matrix (dofs_per_cell, dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        test_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));
+            }
+
+        fe_eval.reinit(cell);
+        for (unsigned int i=0; i<dofs_per_cell;
+             i += VectorizedArray<double>::n_array_elements)
+          {
+            const unsigned int n_items =
+              i+VectorizedArray<double>::n_array_elements > dofs_per_cell ?
+              (dofs_per_cell - i) : VectorizedArray<double>::n_array_elements;
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              fe_eval.begin_dof_values()[j]  = VectorizedArray<double>();
+            for (unsigned int v=0; v<n_items; ++v)
+              fe_eval.begin_dof_values()[i+v][v] = 1.;
+
+            fe_eval.evaluate(true, true);
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                fe_eval.submit_value(10.*fe_eval.get_value(q), q);
+                fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+              }
+            fe_eval.integrate(true, true);
+
+            for (unsigned int v=0; v<n_items; ++v)
+              for (unsigned int j=0; j<dofs_per_cell; ++j)
+                test_matrix(fe_eval.get_internal_dof_numbering()[j],
+                            fe_eval.get_internal_dof_numbering()[i+v])
+                  = fe_eval.begin_dof_values()[j][v];
+          }
+        test_matrix.add(-1., cell_matrix);
+        deallog << test_matrix.frobenius_norm() << " ";
+      }
+    deallog << 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();
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  do_test<dim, fe_degree> (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/assemble_matrix_01.output b/tests/matrix_free/assemble_matrix_01.output
new file mode 100644 (file)
index 0000000..c2eceec
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
diff --git a/tests/matrix_free/assemble_matrix_02.cc b/tests/matrix_free/assemble_matrix_02.cc
new file mode 100644 (file)
index 0000000..5965fa5
--- /dev/null
@@ -0,0 +1,211 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test FEEvaluation for assembling the Stokes matrix with Q2-Q1 finite elements
+
+#include "../tests.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/lac/constraint_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+
+std::ofstream logfile("output");
+
+
+const unsigned int degree_p = 1;
+
+template <int dim>
+void do_test (const DoFHandler<dim> &dof)
+{
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  MappingQ<dim> mapping(degree_p+2);
+
+  // compute matrix with (\nabla v, \nabla u) + (v, 10 * u)
+  {
+    QGauss<dim>  quadrature_formula(degree_p+2);
+
+    FEValues<dim> fe_values (mapping, dof.get_fe(), quadrature_formula,
+                             update_values    |  update_gradients |
+                             update_JxW_values);
+
+    FEEvaluation<dim,degree_p+1,degree_p+2,dim>
+      phi_u (mapping, dof.get_fe(), QGauss<1>(degree_p+2),
+             update_values | update_gradients | update_JxW_values,
+             0);
+    FEEvaluation<dim,degree_p,degree_p+2> phi_p (dof.get_fe(), phi_u, dim);
+
+    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);
+    FullMatrix<double>   test_matrix (dofs_per_cell, dofs_per_cell);
+
+    // The matrix part is from step-22
+    const FEValuesExtractors::Vector velocities (0);
+    const FEValuesExtractors::Scalar pressure (dim);
+    std::vector<SymmetricTensor<2,dim> > phi_grads_u (dofs_per_cell);
+    std::vector<double>                  div_phi_u   (dofs_per_cell);
+    std::vector<double>                  phi_pres    (dofs_per_cell);
+
+    typename DoFHandler<dim>::active_cell_iterator
+    cell = dof.begin_active(),
+    endc = dof.end();
+    for (; cell!=endc; ++cell)
+      {
+        cell_matrix = 0;
+        test_matrix = 0;
+        fe_values.reinit (cell);
+
+        for (unsigned int q=0; q<n_q_points; ++q)
+          {
+            for (unsigned int k=0; k<dofs_per_cell; ++k)
+              {
+                phi_grads_u[k] = fe_values[velocities].symmetric_gradient (k, q);
+                div_phi_u[k]   = fe_values[velocities].divergence (k, q);
+                phi_pres[k]    = fe_values[pressure].value (k, q);
+              }
+
+            for (unsigned int i=0; i<dofs_per_cell; ++i)
+              {
+                for (unsigned int j=0; j<dofs_per_cell; ++j)
+                  {
+                    cell_matrix(i,j) += (phi_grads_u[i] * phi_grads_u[j]
+                                         - div_phi_u[i] * phi_pres[j]
+                                         - phi_pres[i] * div_phi_u[j])
+                                        * fe_values.JxW(q);
+                  }
+              }
+          }
+
+        phi_u.reinit(cell);
+        phi_p.reinit(cell);
+        const unsigned int dofs_per_cell_u = phi_u.dofs_per_cell*dim;
+        const unsigned int dofs_per_cell_p = phi_p.dofs_per_cell;
+        for (unsigned int i=0; i<dofs_per_cell_u;
+             i += VectorizedArray<double>::n_array_elements)
+          {
+            const unsigned int n_items =
+              i+VectorizedArray<double>::n_array_elements > dofs_per_cell_u ?
+              (dofs_per_cell_u - i) : VectorizedArray<double>::n_array_elements;
+            for (unsigned int j=0; j<dofs_per_cell_u; ++j)
+              phi_u.begin_dof_values()[j]  = VectorizedArray<double>();
+            for (unsigned int v=0; v<n_items; ++v)
+              phi_u.begin_dof_values()[i+v][v] = 1.;
+
+            phi_u.evaluate(false, true);
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                VectorizedArray<double> div_u = phi_u.get_divergence(q);
+                phi_u.submit_symmetric_gradient(phi_u.get_symmetric_gradient(q), q);
+                phi_p.submit_value(-div_u, q);
+              }
+            phi_u.integrate(false, true);
+            phi_p.integrate(true, false);
+
+            for (unsigned int v=0; v<n_items; ++v)
+              {
+                for (unsigned int j=0; j<dofs_per_cell_u; ++j)
+                  test_matrix(phi_u.get_internal_dof_numbering()[j],
+                              phi_u.get_internal_dof_numbering()[i+v])
+                    = phi_u.begin_dof_values()[j][v];
+                for (unsigned int j=0; j<dofs_per_cell_p; ++j)
+                  test_matrix(phi_p.get_internal_dof_numbering()[j],
+                              phi_u.get_internal_dof_numbering()[i+v])
+                    = phi_p.begin_dof_values()[j][v];
+              }
+          }
+
+        for (unsigned int i=0; i<dofs_per_cell_p;
+             i += VectorizedArray<double>::n_array_elements)
+          {
+            const unsigned int n_items =
+              i+VectorizedArray<double>::n_array_elements > dofs_per_cell_p ?
+              (dofs_per_cell_p - i) : VectorizedArray<double>::n_array_elements;
+            for (unsigned int j=0; j<dofs_per_cell_p; ++j)
+              phi_p.begin_dof_values()[j]  = VectorizedArray<double>();
+            for (unsigned int v=0; v<n_items; ++v)
+              phi_p.begin_dof_values()[i+v][v] = 1.;
+
+            phi_p.evaluate(true, false);
+            for (unsigned int q=0; q<n_q_points; ++q)
+              {
+                phi_u.submit_divergence(-phi_p.get_value(q), q);
+              }
+            phi_u.integrate(false, true);
+
+            for (unsigned int v=0; v<n_items; ++v)
+              for (unsigned int j=0; j<dofs_per_cell_u; ++j)
+                test_matrix(phi_u.get_internal_dof_numbering()[j],
+                            phi_p.get_internal_dof_numbering()[i+v])
+                  = phi_u.begin_dof_values()[j][v];
+          }
+        test_matrix.add(-1., cell_matrix);
+        deallog << test_matrix.frobenius_norm() << " ";
+      }
+    deallog << std::endl;
+  }
+}
+
+
+
+template <int dim>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_ball (tria);
+  static const HyperBallBoundary<dim> boundary;
+  tria.set_boundary (0, boundary);
+  if (dim == 2)
+    tria.refine_global(1);
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FESystem<dim> fe (FE_Q<dim>(degree_p+1), dim, FE_Q<dim>(degree_p), 1);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  do_test<dim> (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>();
+    deallog.pop();
+    deallog.push("3d");
+    test<3>();
+    deallog.pop();
+  }
+}
diff --git a/tests/matrix_free/assemble_matrix_02.output b/tests/matrix_free/assemble_matrix_02.output
new file mode 100644 (file)
index 0000000..80c38ff
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:2d::Testing FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(1)]
+DEAL:2d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL:3d::Testing FESystem<3>[FE_Q<3>(2)^3-FE_Q<3>(1)]
+DEAL:3d::0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
diff --git a/tests/matrix_free/assemble_matrix_03.cc b/tests/matrix_free/assemble_matrix_03.cc
new file mode 100644 (file)
index 0000000..8697f70
--- /dev/null
@@ -0,0 +1,201 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// test FEEvaluation for assembling the Laplace matrix. Similar to
+// assemble_matrix_01 but doing the whole thing in parallel with WorkStream
+
+
+#include "../tests.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/lac/constraint_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/base/work_stream.h>
+
+
+std::ofstream logfile("output");
+
+
+namespace Assembly
+{
+  namespace Scratch
+  {
+    template <int dim, int fe_degree>
+    struct Data
+    {
+      Data (const FiniteElement<dim> &fe)
+        :
+        fe_values(fe,
+                  QGauss<dim>(fe.degree+1),
+                  update_values    |  update_gradients |
+                  update_JxW_values),
+        fe_eval (1, FEEvaluation<dim,fe_degree>(fe, QGauss<1>(fe.degree+1), fe_values.get_update_flags())),
+        cell_matrix (fe.dofs_per_cell, fe.dofs_per_cell),
+        test_matrix (cell_matrix)
+      {}
+
+      Data (const Data &data)
+        :
+        fe_values(data.fe_values.get_mapping(),
+                  data.fe_values.get_fe(),
+                  data.fe_values.get_quadrature(),
+                  data.fe_values.get_update_flags()),
+        fe_eval (data.fe_eval),
+        cell_matrix (data.cell_matrix),
+        test_matrix (data.test_matrix)
+      {}
+
+      FEValues<dim> fe_values;
+
+      // Need to put FEEvaluation into an aligned vector because of member
+      // variable alignment requirements
+      AlignedVector<FEEvaluation<dim,fe_degree> > fe_eval;
+
+      FullMatrix<double> cell_matrix, test_matrix;
+    };
+  }
+}
+
+
+
+// compute matrix with (\nabla v, \nabla u) + (v, 10 * u)
+template <int dim, int fe_degree>
+void assemble_on_cell (const typename DoFHandler<dim>::active_cell_iterator &cell,
+                       Assembly::Scratch::Data<dim,fe_degree> &data,
+                       unsigned int &)
+{
+  const unsigned int   dofs_per_cell = cell->get_fe().dofs_per_cell;
+  const unsigned int   n_q_points    = data.fe_values.get_quadrature().size();
+  data.cell_matrix = 0;
+  data.test_matrix = 0;
+  data.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)
+          data.cell_matrix(i,j) += ((data.fe_values.shape_grad(i,q_point) *
+                                     data.fe_values.shape_grad(j,q_point)
+                                     +
+                                     10. *
+                                     data.fe_values.shape_value(i,q_point) *
+                                     data.fe_values.shape_value(j,q_point)) *
+                                    data.fe_values.JxW(q_point));
+      }
+
+  FEEvaluation<dim,fe_degree> &fe_eval = data.fe_eval[0];
+  fe_eval.reinit(cell);
+  for (unsigned int i=0; i<dofs_per_cell;
+       i += VectorizedArray<double>::n_array_elements)
+    {
+      const unsigned int n_items =
+        i+VectorizedArray<double>::n_array_elements > dofs_per_cell ?
+        (dofs_per_cell - i) : VectorizedArray<double>::n_array_elements;
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+        fe_eval.begin_dof_values()[j]  = VectorizedArray<double>();
+      for (unsigned int v=0; v<n_items; ++v)
+        fe_eval.begin_dof_values()[i+v][v] = 1.;
+
+      fe_eval.evaluate(true, true);
+      for (unsigned int q=0; q<n_q_points; ++q)
+        {
+          fe_eval.submit_value(10.*fe_eval.get_value(q), q);
+          fe_eval.submit_gradient(fe_eval.get_gradient(q), q);
+        }
+      fe_eval.integrate(true, true);
+
+      for (unsigned int v=0; v<n_items; ++v)
+        for (unsigned int j=0; j<dofs_per_cell; ++j)
+          data.test_matrix(fe_eval.get_internal_dof_numbering()[j],
+                           fe_eval.get_internal_dof_numbering()[i+v])
+            = fe_eval.begin_dof_values()[j][v];
+    }
+  data.test_matrix.add(-1., data.cell_matrix);
+  AssertThrow (data.test_matrix.frobenius_norm() < 1e-10, ExcInternalError());
+}
+
+
+void copy_data_local_to_global (const unsigned int &)
+{}
+
+
+template <int dim, int fe_degree>
+void do_test (const DoFHandler<dim> &dof)
+{
+
+  deallog << "Testing " << dof.get_fe().get_name() << std::endl;
+
+  unsigned int dummy = 0;
+  WorkStream::
+    run (dof.begin_active(), dof.end(),
+         &assemble_on_cell<dim,fe_degree>,
+         &copy_data_local_to_global,
+         Assembly::Scratch::Data<dim,fe_degree>(dof.get_fe()),
+         dummy,
+         2*multithread_info.n_threads(),
+         1);
+  deallog << "OK" << 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);
+  tria.refine_global(5-dim);
+  tria.begin(tria.n_levels()-1)->set_refine_flag();
+  tria.last()->set_refine_flag();
+  tria.execute_coarsening_and_refinement();
+
+  FE_Q<dim> fe (fe_degree);
+  DoFHandler<dim> dof (tria);
+  dof.distribute_dofs(fe);
+  do_test<dim, fe_degree> (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/assemble_matrix_03.output b/tests/matrix_free/assemble_matrix_03.output
new file mode 100644 (file)
index 0000000..6bf454b
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL:2d::Testing FE_Q<2>(1)
+DEAL:2d::OK
+DEAL:2d::Testing FE_Q<2>(2)
+DEAL:2d::OK
+DEAL:2d::Testing FE_Q<2>(4)
+DEAL:2d::OK
+DEAL:3d::Testing FE_Q<3>(1)
+DEAL:3d::OK
+DEAL:3d::Testing FE_Q<3>(2)
+DEAL:3d::OK
index 150d5b1d8e7f9eef81c804e5bcbaa2d9ef853dcc..b1d2e2daae9a56cba3f7bd0dde556daefca9753a 100644 (file)
@@ -21,7 +21,6 @@
 
 #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>
@@ -56,18 +55,17 @@ public:
   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);
+    FEEvaluation<dim,fe_degree,fe_degree+1,1,Number>
+      fe_eval(dof_handler.get_fe(), QGauss<1>(fe_degree+1),
+              update_gradients | update_values | update_JxW_values);
     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.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)
index 5c982cda4c0394a8f46336eed6fe365e66b5fdf7..f964918e9980b28d48d8c7d11d09e34d8b3ff4e6 100644 (file)
@@ -20,7 +20,6 @@
 
 #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>
@@ -56,18 +55,17 @@ public:
   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);
+    FEEvaluation<dim,fe_degree,fe_degree+1,dim,Number> fe_eval
+      (dof_handler.get_fe(), QGauss<1>(fe_degree+1),
+       update_values | update_gradients | update_JxW_values);
     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.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)

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.