From 84e4cc73aeaee4502c6075f6b2dc7d7578e24a1d Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 20 Nov 2014 15:36:18 +0100 Subject: [PATCH] Initialization of FEEvaluation like FEValues 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. --- doc/news/changes.h | 9 + include/deal.II/matrix_free/fe_evaluation.h | 876 +++++++++++------- .../matrix_free/mapping_data_on_the_fly.h | 337 +++++++ .../matrix_free/mapping_fe_evaluation.h | 313 ------- .../matrix_free/shape_info.templates.h | 24 +- tests/matrix_free/assemble_matrix_01.cc | 162 ++++ tests/matrix_free/assemble_matrix_01.output | 11 + tests/matrix_free/assemble_matrix_02.cc | 211 +++++ tests/matrix_free/assemble_matrix_02.output | 5 + tests/matrix_free/assemble_matrix_03.cc | 201 ++++ tests/matrix_free/assemble_matrix_03.output | 11 + tests/matrix_free/matrix_vector_15.cc | 10 +- tests/matrix_free/matrix_vector_16.cc | 10 +- 13 files changed, 1539 insertions(+), 641 deletions(-) create mode 100644 include/deal.II/matrix_free/mapping_data_on_the_fly.h delete mode 100644 include/deal.II/matrix_free/mapping_fe_evaluation.h create mode 100644 tests/matrix_free/assemble_matrix_01.cc create mode 100644 tests/matrix_free/assemble_matrix_01.output create mode 100644 tests/matrix_free/assemble_matrix_02.cc create mode 100644 tests/matrix_free/assemble_matrix_02.output create mode 100644 tests/matrix_free/assemble_matrix_03.cc create mode 100644 tests/matrix_free/assemble_matrix_03.output diff --git a/doc/news/changes.h b/doc/news/changes.h index 98d5e1ba56..56f80def73 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -333,6 +333,15 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/11/20) +
    2. +
  2. 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 diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 9dfa0eabcf..c63d3ec5dc 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -26,11 +26,14 @@ #include #include #include -#include +#include DEAL_II_NAMESPACE_OPEN + + +// forward declarations namespace parallel { namespace distributed @@ -38,15 +41,13 @@ namespace parallel template class Vector; } } - - - -// forward declarations namespace internal { DeclException0 (ExcAccessToUninitializedField); } +template 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 + * without 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 + void reinit (const TriaIterator > &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 without 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 not read from vectors because no DoFHandler information is + * available. + */ + void reinit (const typename Triangulation::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::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::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 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::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::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 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::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::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 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::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::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 void set_dof_values (VectorType &dst) const; @@ -539,6 +569,14 @@ public: */ VectorizedArray *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 & + 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::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 &geometry, - const DoFHandler &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 + FEEvaluationBase (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *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 > mapped_geometry; + std_cxx1x::shared_ptr > 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 *dof_handler; + std::vector 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 local_dof_indices; + + /** + * Make other FEEvaluationBase as well as FEEvaluation objects friends. + */ + template friend class FEEvaluationBase; + template 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::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 &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); + * Constructor with reduced functionality for similar usage of FEEvaluation + * as FEValues, including matrix assembly. + */ + template + FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *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::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 &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); + * Constructor with reduced functionality for similar usage of FEEvaluation + * as FEValues, including matrix assembly. + */ + template + FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *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::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 &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); + * Constructor with reduced functionality for similar usage of FEEvaluation + * as FEValues, including matrix assembly. + */ + template + FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *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::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 + 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). + *

    Usage and initialization

    + * + *

    Fast usage in combination with MatrixFree

    * - * 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. + * + *

    Usage without pre-initialized MatrixFree object

    + * + * 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 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::n_array_elements) + * { + * const unsigned int n_items = + * i+VectorizedArray::n_array_elements > dofs_per_cell ? + * (dofs_per_cell - i) : VectorizedArray::n_array_elements; + * + * // Set n_items unit vectors + * for (unsigned int j=0; j(); + * for (unsigned int v=0; vDescription of evaluation routines + * + * 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 dim (fe_degree+1)dim+1 + * operations instead of the slower + * (fe_degree+1)2*dim 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 not 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 +template class FEEvaluation : public FEEvaluationAccess { 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::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 &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); + * 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). + */ + FEEvaluation (const Mapping &mapping, + const FiniteElement &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 &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 + FEEvaluation (const FiniteElement &fe, + const FEEvaluationBase &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 &geometry, - const DoFHandler &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 &geometry, - const DoFHandler &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 &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED -: - BaseClass (geometry, dof_handler, first_selected_component) - {} }; @@ -1753,7 +1853,6 @@ FEEvaluationBase 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 template +template inline FEEvaluationBase -::FEEvaluationBase (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler_in, - const unsigned int first_selected_component) +::FEEvaluationBase (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *other) : quad_no (-1), n_fe_components (n_components_), @@ -1801,24 +1904,25 @@ FEEvaluationBase matrix_info (0), dof_info (0), mapping_info (0), - stored_shape_info (new internal::MatrixFreeFunctions::ShapeInfo(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(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 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(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 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 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(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::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::reinit (const unsigned int cell_in) +template +template +inline +void +FEEvaluationBase +::reinit (const TriaIterator > &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::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 +inline +void +FEEvaluationBase +::reinit (const typename Triangulation::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 inline unsigned int @@ -2278,13 +2448,7 @@ FEEvaluationBase // process constraints and need not care about vectorization if (matrix_info == 0) { - Assert (dof_handler != 0, ExcNotInitialized()); - typename DoFHandler::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 /*------------------------------ access to data fields ----------------------*/ +template +inline +const std::vector & +FEEvaluationBase:: +get_internal_dof_numbering() const +{ + return data->lexicographic_numbering; +} + + + template inline @@ -3620,13 +3795,18 @@ FEEvaluationAccess template +template inline FEEvaluationAccess -::FEEvaluationAccess (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) +::FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *other) : - FEEvaluationBase (geometry, dof_handler, first_selected_component) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} @@ -3661,13 +3841,18 @@ FEEvaluationAccess template +template inline FEEvaluationAccess -::FEEvaluationAccess (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) +::FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *other) : - FEEvaluationBase (geometry, dof_handler, first_selected_component) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} @@ -3895,13 +4080,18 @@ FEEvaluationAccess template +template inline FEEvaluationAccess -::FEEvaluationAccess (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) +::FEEvaluationAccess (const Mapping &mapping, + const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBase *other) : - FEEvaluationBase (geometry, dof_handler, first_selected_component) + FEEvaluationBase (mapping, fe, quadrature, update_flags, + first_selected_component, other) {} @@ -4250,13 +4440,18 @@ FEEvaluationAccess<1,1,Number> template +template 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 inline FEEvaluation -::FEEvaluation (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) +::FEEvaluation (const Mapping &mapping, + const FiniteElement &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*>(0)), + dofs_per_cell (this->data->dofs_per_cell) +{ + check_template_arguments(numbers::invalid_unsigned_int); + set_data_pointers(); +} + + + +template +inline +FEEvaluation +::FEEvaluation (const FiniteElement &fe, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component) + : + BaseClass (StaticMappingQ1::mapping, fe, quadrature, update_flags, + first_selected_component, + static_cast*>(0)), + dofs_per_cell (this->data->dofs_per_cell) +{ + check_template_arguments(numbers::invalid_unsigned_int); + set_data_pointers(); +} + + + +template +template +inline +FEEvaluation +::FEEvaluation (const FiniteElement &fe, + const FEEvaluationBase &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 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 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 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 { 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 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 index 0000000000..678b971c1b --- /dev/null +++ b/include/deal.II/matrix_free/mapping_data_on_the_fly.h @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include + + +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 + 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 &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::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::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 &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 > > & + 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 > & + 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 > > & + 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 > > & + 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::cell_iterator present_cell; + + /** + * Dummy finite element object necessary for initializing the FEValues + * object. + */ + FE_Nothing fe_dummy; + + /** + * An underlying FEValues object that performs the (scalar) evaluation. + */ + FEValues 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 > > inverse_jacobians; + + /** + * Stored Jacobian determinants and quadrature weights + */ + AlignedVector > jxw_values; + + /** + * Stored quadrature points + */ + AlignedVector > > quadrature_points; + + /** + * Stored normal vectors (for face integration) + */ + AlignedVector > > normal_vectors; + }; + + + /*----------------------- Inline functions ----------------------------------*/ + + template + inline + MappingDataOnTheFly::MappingDataOnTheFly (const Mapping &mapping, + const Quadrature<1> &quadrature, + const UpdateFlags update_flags) + : + fe_values(mapping, fe_dummy, Quadrature(quadrature), + internal::MatrixFreeFunctions::MappingInfo::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 + inline + MappingDataOnTheFly::MappingDataOnTheFly (const Quadrature<1> &quadrature, + const UpdateFlags update_flags) + : + fe_values(fe_dummy, Quadrature(quadrature), + internal::MatrixFreeFunctions::MappingInfo::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 + inline + void + MappingDataOnTheFly::reinit(typename dealii::Triangulation::cell_iterator cell) + { + if (present_cell == cell) + return; + present_cell = cell; + fe_values.reinit(present_cell); + for (unsigned int q=0; q + inline + bool + MappingDataOnTheFly::is_initialized() const + { + return present_cell != typename dealii::Triangulation::cell_iterator(); + } + + + + template + inline + typename dealii::Triangulation::cell_iterator + MappingDataOnTheFly::get_cell() const + { + return fe_values.get_cell(); + } + + + + template + inline + const FEValues & + MappingDataOnTheFly::get_fe_values() const + { + return fe_values; + } + + + + template + inline + const AlignedVector > > & + MappingDataOnTheFly::get_inverse_jacobians() const + { + return inverse_jacobians; + } + + + + template + inline + const AlignedVector > > & + MappingDataOnTheFly::get_normal_vectors() const + { + return normal_vectors; + } + + + + template + inline + const AlignedVector > > & + MappingDataOnTheFly::get_quadrature_points() const + { + return quadrature_points; + } + + + + template + inline + const AlignedVector > & + MappingDataOnTheFly::get_JxW_values() const + { + return jxw_values; + } + + + + template + inline + const Quadrature<1> & + MappingDataOnTheFly::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 index 626fc78cba..0000000000 --- a/include/deal.II/matrix_free/mapping_fe_evaluation.h +++ /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 -#include -#include -#include -#include -#include -#include -#include -#include - - -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 -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 &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 - void reinit(ITERATOR cell); - - /** - * Return a triangulation iterator to the current cell. - */ - typename Triangulation::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 &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 > > & - 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 > & - 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 > > & - 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 > > & - 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 fe_dummy; - - /** - * An underlying FEValues object that performs the (scalar) evaluation. - */ - FEValues 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 > > inverse_jacobians; - - /** - * Stored Jacobian determinants and quadrature weights - */ - AlignedVector > jxw_values; - - /** - * Stored quadrature points - */ - AlignedVector > > quadrature_points; - - /** - * Stored normal vectors (for face integration) - */ - AlignedVector > > normal_vectors; -}; - - -/*----------------------- Inline functions ----------------------------------*/ - -template -inline -MappingFEEvaluation::MappingFEEvaluation (const Mapping &mapping, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags) - : - fe_values(mapping, fe_dummy, Quadrature(quadrature), - internal::MatrixFreeFunctions::MappingInfo::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 -inline -MappingFEEvaluation::MappingFEEvaluation (const Quadrature<1> &quadrature, - const UpdateFlags update_flags) - : - fe_values(fe_dummy, Quadrature(quadrature), - internal::MatrixFreeFunctions::MappingInfo::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 -template -inline -void -MappingFEEvaluation::reinit(ITERATOR cell) -{ - fe_values.reinit(static_cast::cell_iterator>(cell)); - for (unsigned int q=0; q -inline -typename Triangulation::cell_iterator -MappingFEEvaluation::get_cell() const -{ - return fe_values.get_cell(); -} - - - -template -inline -const FEValues & -MappingFEEvaluation::get_fe_values() const -{ - return fe_values; -} - - - -template -inline -const AlignedVector > > & -MappingFEEvaluation::get_inverse_jacobians() const -{ - return inverse_jacobians; -} - - - -template -inline -const AlignedVector > > & -MappingFEEvaluation::get_normal_vectors() const -{ - return normal_vectors; -} - - - -template -inline -const AlignedVector > > & -MappingFEEvaluation::get_quadrature_points() const -{ - return quadrature_points; -} - - - -template -inline -const AlignedVector > & -MappingFEEvaluation::get_JxW_values() const -{ - return jxw_values; -} - - - -template -inline -const Quadrature<1> & -MappingFEEvaluation::get_quadrature() const -{ - return quadrature_1d; -} - - -#ifndef DOXYGEN - - -#endif // ifndef DOXYGEN - - -DEAL_II_NAMESPACE_CLOSE - -#endif diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index bcbaac2c3f..8b3aa6fa81 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -116,15 +116,27 @@ namespace internal // and invert back std::vector scalar_inv = Utilities::invert_permutation(scalar_lexicographic); - std::vector lexicographic (fe_in.dofs_per_cell); - for (unsigned int comp=0; comp lexicographic(fe_in.dofs_per_cell, + numbers::invalid_unsigned_int); + unsigned int components_before = 0; + for (unsigned int e=0; edofs_per_cell); + for (unsigned int i=0; i +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +std::ofstream logfile("output"); + + + +template +void do_test (const DoFHandler &dof) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + + MappingQ mapping(fe_degree+1); + + // compute matrix with (\nabla v, \nabla u) + (v, 10 * u) + { + QGauss quadrature_formula(fe_degree+1); + + FEValues fe_values (mapping, dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + FEEvaluation + 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 cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix test_matrix (dofs_per_cell, dofs_per_cell); + + typename DoFHandler::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_array_elements) + { + const unsigned int n_items = + i+VectorizedArray::n_array_elements > dofs_per_cell ? + (dofs_per_cell - i) : VectorizedArray::n_array_elements; + for (unsigned int j=0; j(); + for (unsigned int v=0; v +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary 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 fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + do_test (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 index 0000000000..c2eceeca3b --- /dev/null +++ b/tests/matrix_free/assemble_matrix_01.output @@ -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 index 0000000000..5965fa5343 --- /dev/null +++ b/tests/matrix_free/assemble_matrix_02.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +std::ofstream logfile("output"); + + +const unsigned int degree_p = 1; + +template +void do_test (const DoFHandler &dof) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + MappingQ mapping(degree_p+2); + + // compute matrix with (\nabla v, \nabla u) + (v, 10 * u) + { + QGauss quadrature_formula(degree_p+2); + + FEValues fe_values (mapping, dof.get_fe(), quadrature_formula, + update_values | update_gradients | + update_JxW_values); + + FEEvaluation + phi_u (mapping, dof.get_fe(), QGauss<1>(degree_p+2), + update_values | update_gradients | update_JxW_values, + 0); + FEEvaluation 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 cell_matrix (dofs_per_cell, dofs_per_cell); + FullMatrix 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 > phi_grads_u (dofs_per_cell); + std::vector div_phi_u (dofs_per_cell); + std::vector phi_pres (dofs_per_cell); + + typename DoFHandler::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_array_elements) + { + const unsigned int n_items = + i+VectorizedArray::n_array_elements > dofs_per_cell_u ? + (dofs_per_cell_u - i) : VectorizedArray::n_array_elements; + for (unsigned int j=0; j(); + for (unsigned int v=0; v 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_array_elements) + { + const unsigned int n_items = + i+VectorizedArray::n_array_elements > dofs_per_cell_p ? + (dofs_per_cell_p - i) : VectorizedArray::n_array_elements; + for (unsigned int j=0; j(); + for (unsigned int v=0; v +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary 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 fe (FE_Q(degree_p+1), dim, FE_Q(degree_p), 1); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + do_test (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 index 0000000000..80c38ff210 --- /dev/null +++ b/tests/matrix_free/assemble_matrix_02.output @@ -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 index 0000000000..8697f70bb3 --- /dev/null +++ b/tests/matrix_free/assemble_matrix_03.cc @@ -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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +std::ofstream logfile("output"); + + +namespace Assembly +{ + namespace Scratch + { + template + struct Data + { + Data (const FiniteElement &fe) + : + fe_values(fe, + QGauss(fe.degree+1), + update_values | update_gradients | + update_JxW_values), + fe_eval (1, FEEvaluation(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 fe_values; + + // Need to put FEEvaluation into an aligned vector because of member + // variable alignment requirements + AlignedVector > fe_eval; + + FullMatrix cell_matrix, test_matrix; + }; + } +} + + + +// compute matrix with (\nabla v, \nabla u) + (v, 10 * u) +template +void assemble_on_cell (const typename DoFHandler::active_cell_iterator &cell, + Assembly::Scratch::Data &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 &fe_eval = data.fe_eval[0]; + fe_eval.reinit(cell); + for (unsigned int i=0; i::n_array_elements) + { + const unsigned int n_items = + i+VectorizedArray::n_array_elements > dofs_per_cell ? + (dofs_per_cell - i) : VectorizedArray::n_array_elements; + for (unsigned int j=0; j(); + for (unsigned int v=0; v +void do_test (const DoFHandler &dof) +{ + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + unsigned int dummy = 0; + WorkStream:: + run (dof.begin_active(), dof.end(), + &assemble_on_cell, + ©_data_local_to_global, + Assembly::Scratch::Data(dof.get_fe()), + dummy, + 2*multithread_info.n_threads(), + 1); + deallog << "OK" << std::endl; +} + + + +template +void test () +{ + Triangulation tria; + GridGenerator::hyper_ball (tria); + static const HyperBallBoundary 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 fe (fe_degree); + DoFHandler dof (tria); + dof.distribute_dofs(fe); + do_test (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 index 0000000000..6bf454b955 --- /dev/null +++ b/tests/matrix_free/assemble_matrix_03.output @@ -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 diff --git a/tests/matrix_free/matrix_vector_15.cc b/tests/matrix_free/matrix_vector_15.cc index 150d5b1d8e..b1d2e2daae 100644 --- a/tests/matrix_free/matrix_vector_15.cc +++ b/tests/matrix_free/matrix_vector_15.cc @@ -21,7 +21,6 @@ #include "../tests.h" #include -#include #include #include #include @@ -56,18 +55,17 @@ public: void vmult (VECTOR &dst, const VECTOR &src) const { - MappingFEEvaluation mapping(QGauss<1>(fe_degree+1), - update_gradients | update_values | - update_JxW_values); VECTOR src_cpy = src; constraints.distribute(src_cpy); - FEEvaluation fe_eval(mapping, dof_handler); + FEEvaluation + fe_eval(dof_handler.get_fe(), QGauss<1>(fe_degree+1), + update_gradients | update_values | update_JxW_values); dst = 0; typename DoFHandler::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 -#include #include #include #include @@ -56,18 +55,17 @@ public: void vmult (VECTOR &dst, const VECTOR &src) const { - MappingFEEvaluation mapping(QGauss<1>(fe_degree+1), - update_gradients | update_values | - update_JxW_values); VECTOR src_cpy = src; constraints.distribute(src_cpy); - FEEvaluation fe_eval(mapping, dof_handler); + FEEvaluation fe_eval + (dof_handler.get_fe(), QGauss<1>(fe_degree+1), + update_values | update_gradients | update_JxW_values); dst = 0; typename DoFHandler::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