From: kronbichler Date: Fri, 27 Jun 2014 20:31:28 +0000 (+0000) Subject: Simplify FEEvaluation routines considerably by merging all functionality into FEEvalu... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=45cbf819fd339315507b9676b988431ce1f0c8d7;p=dealii-svn.git Simplify FEEvaluation routines considerably by merging all functionality into FEEvaluation. The internal implementation is done in internal::FEEvaluationImpl that provides the same efficient kernels as before, but frees the user from tedious selection of the appropriate kernel for the element in use. git-svn-id: https://svn.dealii.org/trunk@33094 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c22b598537..ba57ec5880 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -148,6 +148,20 @@ inconvenience this causes.

Specific improvements

    +
  1. Simplified interfaces for FEEvaluation: Previously, the user had to + select the appropriate kernel (FEEvaluation, FEEvaluationGeneral, + FEEvaluationDGP, FEEvaluationGL) for the matrix-free evaluation + routines. This made it difficult to write compact code that needs to select + between different elements. Therefore, all the functionality has been merged + into FEEvaluation and one should only use FEEvaluation, while the other + interfaces have been marked deprecated and will be removed in a future + version. The internal data structures for the various special cases have + been kept in order to provide for the most efficient routines, which is + selected at construction of FEEvaluation. +
    + (Martin Kronbichler, 2014/06/27) +
    2. +
  2. Bugfix: TBB specific include directories have to be added to the list of user include directories because of direct inclusion of header files in base/thread_local_storage.h. @@ -180,8 +194,8 @@ inconvenience this causes. (Denis Davydov, 2014/06/19)
    3. -
  3. New: There is now a class FEEvaluationQ_DG0 that does - optimized matrix-free evaluation for FE_Q_DG0 elements. +
  4. New: The class FEEvaluation now provides optimized matrix-free + evaluation routines for FE_Q_DG0 elements.
    (Martin Kronbichler, 2014/06/13)
    5. diff --git a/deal.II/examples/step-48/doc/intro.dox b/deal.II/examples/step-48/doc/intro.dox index 6a89c11ca8..a4715e6f80 100644 --- a/deal.II/examples/step-48/doc/intro.dox +++ b/deal.II/examples/step-48/doc/intro.dox @@ -127,7 +127,7 @@ variable constraint_indicator at hand that holds, for each cell, the local indices of DoFs that are constrained as well as the identifier of the type of constraint. Actually, you will not see these data structures in the example program since the class -FEEvaluationGL takes care of the constraints without user +FEEvaluation takes care of the constraints without user interaction. diff --git a/deal.II/examples/step-48/step-48.cc b/deal.II/examples/step-48/step-48.cc index 77e1f71970..1c28434238 100644 --- a/deal.II/examples/step-48/step-48.cc +++ b/deal.II/examples/step-48/step-48.cc @@ -125,7 +125,7 @@ namespace Step48 data.initialize_dof_vector (inv_mass_matrix); - FEEvaluationGL fe_eval(data); + FEEvaluation fe_eval(data); const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell=0; cell &cell_range) const { AssertDimension (src.size(), 2); - FEEvaluationGL current (data), old (data); + FEEvaluation current (data), old (data); for (unsigned int cell=cell_range.first; cell &matrix_free, const unsigned int fe_no, const unsigned int quad_no, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points); /** @@ -863,7 +862,7 @@ protected: FEEvaluationAccess (const MatrixFree &matrix_free, const unsigned int fe_no, const unsigned int quad_no, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points); /** @@ -1020,7 +1019,7 @@ protected: FEEvaluationAccess (const MatrixFree &matrix_free, const unsigned int fe_no, const unsigned int quad_no, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points); /** @@ -1227,196 +1226,6 @@ protected: -/** - * The class that provides all functions necessary to evaluate functions at - * quadrature points and cell integrations. In functionality, this class is - * similar to FEValues, however, it includes a lot of specialized - * functions that make it much faster (between 5 and 500 times as fast, - * depending on the polynomial order). Access to the data fields is provided - * through functionality in the class FEEvaluationAccess. - * - * This class is designed for general local finite element operations based on - * tensor products of 1D polynomials and quadrature points. Often, there are - * some symmetries or zeros in the unit cell data that allow for a more - * efficient operator application. FEEvaluation is specialized to standard - * FE_Q/FE_DGQ elements and quadrature points symmetric around 0.5 (like Gauss - * quadrature), and hence the most common situation. FEEvaluationGL is a - * specialization for elements where quadrature formula and support points are - * chosen so that a orthogonal relation between the values holds. This is the - * case for FE_Q based on Gauss-Lobatto support points with Gauss-Lobatto - * quadrature formula of the same order. - * - * This class has five template arguments: - * - * @param dim Dimension in which this class is to be used - * - * @param fe_degree Degree of the tensor product finite element with - * 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 - * - * @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) - * - * @param Number Number format, usually @p double or @p float - * - * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 - */ -template -class FEEvaluationGeneral : public FEEvaluationAccess -{ -public: - typedef FEEvaluationAccess BaseClass; - typedef Number number_type; - typedef typename BaseClass::value_type value_type; - typedef typename BaseClass::gradient_type gradient_type; - static const unsigned int dimension = dim; - static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::value; - - /** - * Constructor. Takes all data stored in MatrixFree. If applied to problems - * with more than one finite element or more than one quadrature formula - * selected during construction of @p matrix_free, @p fe_no and @p quad_no - * allow to select the appropriate components. - */ - FEEvaluationGeneral (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); - - /** - * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type MappingFEEvaluation - * and a DoFHandler in order to allow for reading out the finite element - * data. It uses the data provided by dof_handler.get_fe(). If the element - * is vector-valued, the optional argument allows to specify the index of - * the base element (as long as the element is primitive, non-primitive are - * not supported currently). - * - * With initialization from a FEValues object, no call to a reinit method of - * this class is necessary. Instead, it is enough if the geometry is - * initialized to a given cell iterator. It can also read from or write to - * vectors in the standard way for DoFHandler::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). - */ - FEEvaluationGeneral (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); - - /** - * Copy constructor - */ - FEEvaluationGeneral (const FEEvaluationGeneral &other); - - /** - * Evaluates the function values, the gradients, and the Laplacians of the - * FE function given at the DoF values in the input vector at the quadrature - * points. The function arguments specify which parts shall actually be - * computed. Needs to be called before the functions @p get_value(), @p - * get_gradient() or @p get_laplacian return useful information. - */ - void evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_hess = false); - - /** - * This function takes the values and/or gradients that are stored on - * quadrature points, tests them by all the basis functions/gradients on the - * cell and performs the cell integration. The two function arguments @p - * integrate_val and @p integrate_grad are used to enable/disable some of - * values or gradients. - */ - void integrate (const bool integrate_val, - const bool integrate_grad); - - /** - * Returns the q-th quadrature point stored in MappingInfo. - */ - Point > - quadrature_point (const unsigned int q_point) const; - - /** - * Internal function that applies the function values of the tensor product - * in a given coordinate direction (first template argument), from - * polynomials to values on quadrature points (second flag set to true) or - * in an integration loop from values on quadrature points to values tested - * by different test function (second flag set to false), and if the result - * is to be added to previous content in the data fields or not. - */ - template - void apply_values (const VectorizedArray in [], - VectorizedArray out []); - - /** - * Internal function that applies the gradient operation of the tensor - * product in a given coordinate direction (first template argument), from - * polynomials to values on quadrature points (second flag set to true) or - * in an integration loop from values on quadrature points to values tested - * by different test function (second flag set to false), and if the result - * is to be added to previous content in the data fields or not. - */ - template - void apply_gradients (const VectorizedArray in [], - VectorizedArray out []); - - /** - * Internal function that applies the second derivative operation (Hessian) - * of the tensor product in a given coordinate direction (first template - * argument), from polynomials to values on quadrature points (second flag - * set to true) or in an integration loop from values on quadrature points - * to values tested by different test function (second flag set to false), - * and if the result is to be added to previous content in the data fields - * or not. - */ - template - void apply_hessians (const VectorizedArray in [], - VectorizedArray out []); - -protected: - - /** - * Internally stored variables for the different data fields. - */ - VectorizedArray my_data_array[n_components*(dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)]; - - /** - * Checks the number of cell dofs with the value that we expect - */ - void check_dofs_per_cell (const unsigned int dofs_per_cell) const; - -private: - /** - * Sets the pointers from the data array to values_dof, etc. - */ - void set_data_pointers(); - - /** - * Checks that the element number of degrees of freedom given by the - * template arguments (via fe_degree) coincides with the number of degrees - * of freedom in the stored shape values. - */ - void check_template_arguments(const unsigned int fe_no); -}; - - - /** * The class that provides all functions necessary to evaluate functions at * quadrature points and cell integrations. In functionality, this class is @@ -1434,14 +1243,16 @@ private: * possibly be vectorized by combining several operations (e.g. when * performing matrix assembly). * - * This class is a specialization of FEEvaluationGeneral designed for standard - * FE_Q or FE_DGQ elements and quadrature points symmetric around 0.5 (like - * Gauss quadrature), and hence the most common situation. 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 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. * * This class assumes that shape functions of the FiniteElement under * consideration do not depend on the actual shape of the cells in @@ -1470,18 +1281,17 @@ private: */ template -class FEEvaluation : - public FEEvaluationGeneral +class FEEvaluation : public FEEvaluationAccess { public: - typedef FEEvaluationGeneral BaseClass; + typedef FEEvaluationAccess BaseClass; typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell; - static const unsigned int n_q_points = BaseClass::n_q_points; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + static const unsigned int tensor_dofs_per_cell = Utilities::fixed_int_power::value; /** * Constructor. Takes all data stored in MatrixFree. If applied to problems @@ -1545,92 +1355,104 @@ public: const bool integrate_grad); /** - * Internal function that applies the function values of the tensor product - * in a given coordinate direction (first template argument), from - * polynomials to values on quadrature points (second flag set to true) or - * in an integration loop from values on quadrature points to values tested - * by different test function (second flag set to false), and if the result - * is to be added to previous content in the data fields or not. + * Returns the q-th quadrature point stored in MappingInfo. */ - template - void apply_values (const VectorizedArray in [], - VectorizedArray out []); + Point > + quadrature_point (const unsigned int q_point) const; /** - * Internal function that applies the gradient operation of the tensor - * product in a given coordinate direction (first template argument), from - * polynomials to values on quadrature points (second flag set to true) or - * in an integration loop from values on quadrature points to values tested - * by different test function (second flag set to false), and if the result - * is to be added to previous content in the data fields or not. + * The number of scalar degrees of freedom on the cell. Usually close to + * tensor_dofs_per_cell, but depends on the actual element selected and is + * thus not static. */ - template - void apply_gradients (const VectorizedArray in [], - VectorizedArray out []); + const unsigned int dofs_per_cell; +private: /** - * Internal function that applies the second derivative operation (Hessian) - * of the tensor product in a given coordinate direction (first template - * argument), from polynomials to values on quadrature points (second flag - * set to true) or in an integration loop from values on quadrature points - * to values tested by different test function (second flag set to false), - * and if the result is to be added to previous content in the data fields - * or not. + * Internally stored variables for the different data fields. */ - template - void apply_hessians (const VectorizedArray in [], - VectorizedArray out []); + VectorizedArray my_data_array[n_components*(tensor_dofs_per_cell+1+(dim*dim+2*dim+1)*n_q_points)]; -protected: - VectorizedArray shape_val_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; - VectorizedArray shape_gra_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; - VectorizedArray shape_hes_evenodd[fe_degree+1][(n_q_points_1d+1)/2]; + /** + * Checks if the template arguments regarding degree of the element + * corresponds to the actual element used at initialization. + */ + void check_template_arguments(const unsigned int fe_no); + + /** + * Sets the pointers of the base class to my_data_array. + */ + void set_data_pointers(); + + /** + * Function pointer for the evaluate function + */ + void (*evaluate_funct) (const internal::MatrixFreeFunctions::ShapeInfo &, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl); -private: /** - * Fills the fields shape_???_evenodd, called in the constructor. + * Function pointer for the integrate function */ - void compute_even_odd_factors(); + void (*integrate_funct)(const internal::MatrixFreeFunctions::ShapeInfo &, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + const bool evaluate_val, + const bool evaluate_grad); }; /** - * The class that provides all functions necessary to evaluate functions at - * quadrature points and cell integrations. In functionality, this class is - * similar to FEValues, however, it includes a lot of specialized - * functions that make it much faster (between 5 and 500, depending on the - * polynomial order). - * - * This class is a specialization of FEEvaluation for elements where - * quadrature formula and support points are chosen so that a orthonormal - * relation between the values holds. This is the case for FE_Q based on - * Gauss-Lobatto support points with Gauss-Lobatto quadrature formula of the - * same order (QGaussLobatto). In that case, application of values is trivial - * (as the transformation is the identity matrix), and application of - * gradients is considerably simpler (since all value applications in - * directions other than the gradient direction are again identity - * operations). - * - * This class has four template arguments: - * - * @param dim Dimension in which this class is to be used - * - * @param fe_degree Degree of the tensor product finite element with - * fe_degree+1 degrees of freedom per coordinate - * direction. The quadrature formula is tied to the choice of - * the element by setting n_q_points_1d = fe_degree+1, which - * gives a diagonal mass matrix - * - * @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) - * - * @param Number Number format, usually @p double or @p float - * - * @author Katharina Kormann and Martin Kronbichler, 2010, 2011 + * Deprecated. Functionality has been merged into FEEvaluation. Use class + * FEEvaluation instead. + */ +template +class FEEvaluationGeneral : public FEEvaluation +{ +public: + typedef FEEvaluation BaseClass; + typedef Number number_type; + typedef typename BaseClass::value_type value_type; + typedef typename BaseClass::gradient_type gradient_type; + static const unsigned int dimension = dim; + static const unsigned int n_components = n_components_; + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + + /** + * Constructor. + */ + FEEvaluationGeneral (const MatrixFree &matrix_free, + const unsigned int fe_no = 0, + const unsigned int quad_no = 0) DEAL_II_DEPRECATED + : + 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) + {} +}; + + + +/** + * Deprecated. Functionality has been merged into FEEvaluation. Use class + * FEEvaluation instead. */ template class FEEvaluationGL : @@ -1643,81 +1465,28 @@ public: typedef typename BaseClass::gradient_type gradient_type; static const unsigned int dimension = dim; static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell; + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; static const unsigned int n_q_points = BaseClass::n_q_points; /** - * Constructor. Takes all data stored in MatrixFree. If applied to problems - * with more than one finite element or more than one quadrature formula - * selected during construction of @p matrix_free, @p fe_no and @p quad_no - * allow to select the appropriate components. + * Constructor. */ FEEvaluationGL (const MatrixFree &matrix_free, const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int quad_no = 0) DEAL_II_DEPRECATED + : + BaseClass (matrix_free, fe_no, quad_no) + {} /** - * 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). + * Constructor. */ FEEvaluationGL (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); - - /** - * Copy constructor - */ - FEEvaluationGL (const FEEvaluationGL &other); - - /** - * Evaluates the function values, the gradients, and the Hessians of the FE - * function given at the DoF values in the input vector at the quadrature - * points of the unit cell. The function arguments specify which parts shall - * actually be computed. Needs to be called before the functions @p - * get_value(), @p get_gradient() or @p get_laplacian give useful - * information (unless these values have been set manually). - */ - void evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl = false); - - /** - * This function takes the values and/or gradients that are stored on - * quadrature points, tests them by all the basis functions/gradients on the - * cell and performs the cell integration. The two function arguments @p - * integrate_val and @p integrate_grad are used to enable/disable some of - * values or gradients. - */ - void integrate (const bool integrate_val, - const bool integrate_grad); - - /** - * Internal function that applies the gradient operation of the tensor - * product in a given coordinate direction (first template argument), from - * polynomials to values on quadrature points (second flag set to true) or - * in an integration loop from values on quadrature points to values tested - * by different test function (second flag set to false), and if the result - * is to be added to some previous results or not. - */ - template - void apply_gradients (const VectorizedArray in [], - VectorizedArray out []); + const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED + : + BaseClass (geometry, dof_handler, first_selected_component) + {} }; @@ -1748,27 +1517,16 @@ namespace internal /** - * The class that provides all functions necessary to evaluate functions at - * quadrature points and cell integrations. In functionality, this class is - * similar to FEValues, however, it includes a lot of specialized - * functions that make it much faster (between 5 and 500 times as fast, - * depending on the polynomial order). Access to the data fields is provided - * through functionality in the class FEEvaluationAccess. - * - * This class is an extension of FEEvaluationGeneral to work with elements of - * complete polynomial degree p, FE_DGP. In this case, the polynomial basis is - * a truncated tensor product, so the evaluate and integrate routines use a - * truncation. - * - * @author Martin Kronbichler, 2014 + * Functionality has been merged into FEEvaluation. Use class FEEvaluation + * instead. */ template class FEEvaluationDGP : - public FEEvaluationGeneral + public FEEvaluation { public: - typedef FEEvaluationGeneral BaseClass; + typedef FEEvaluation BaseClass; typedef Number number_type; typedef typename BaseClass::value_type value_type; typedef typename BaseClass::gradient_type gradient_type; @@ -1778,160 +1536,28 @@ public: static const unsigned int n_q_points = BaseClass::n_q_points; /** - * Constructor. Takes all data stored in MatrixFree. If applied to problems - * with more than one finite element or more than one quadrature formula - * selected during construction of @p matrix_free, @p fe_no and @p quad_no - * allow to select the appropriate components. + * Constructor. */ FEEvaluationDGP (const MatrixFree &matrix_free, const unsigned int fe_no = 0, - const unsigned int quad_no = 0); + const unsigned int quad_no = 0) DEAL_II_DEPRECATED + : + BaseClass (matrix_free, fe_no, quad_no) + {} /** - * 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). + * Constructor. */ FEEvaluationDGP (const MappingFEEvaluation &geometry, const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); - - /** - * Copy constructor - */ - FEEvaluationDGP (const FEEvaluationDGP &other); - - /** - * Evaluates the function values, the gradients, and the Hessians of the FE - * function given at the DoF values in the input vector at the quadrature - * points of the unit cell. The function arguments specify which parts shall - * actually be computed. Needs to be called before the functions @p - * get_value(), @p get_gradient() or @p get_laplacian give useful - * information (unless these values have been set manually). - */ - void evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl = false); - - /** - * This function takes the values and/or gradients that are stored on - * quadrature points, tests them by all the basis functions/gradients on the - * cell and performs the cell integration. The two function arguments @p - * integrate_val and @p integrate_grad are used to enable/disable some of - * values or gradients. - */ - void integrate (const bool integrate_val, - const bool integrate_grad); + const unsigned int first_selected_component = 0) DEAL_II_DEPRECATED + : + BaseClass (geometry, dof_handler, first_selected_component) + {} }; -/** - * The class that provides all functions necessary to evaluate functions at - * quadrature points and cell integrations. In functionality, this class is - * similar to FEValues, however, it includes a lot of specialized - * functions that make it much faster (between 5 and 500 times as fast, - * depending on the polynomial order). Access to the data fields is provided - * through functionality in the class FEEvaluationAccess. - * - * This class is an extension of FEEvaluation to work with continuous elements - * supplemented with a single discontinuous degree of freedom, i.e., FE_Q_DG0. - * - * @author Martin Kronbichler, 2014 - */ -template -class FEEvaluationQ_DG0 : - public FEEvaluation -{ -public: - typedef FEEvaluation BaseClass; - typedef Number number_type; - typedef typename BaseClass::value_type value_type; - typedef typename BaseClass::gradient_type gradient_type; - static const unsigned int dimension = dim; - static const unsigned int n_components = n_components_; - static const unsigned int dofs_per_cell = BaseClass::dofs_per_cell + 1; - static const unsigned int n_q_points = BaseClass::n_q_points; - - /** - * Constructor. Takes all data stored in MatrixFree. If applied to problems - * with more than one finite element or more than one quadrature formula - * selected during construction of @p matrix_free, @p fe_no and @p quad_no - * allow to select the appropriate components. - */ - FEEvaluationQ_DG0 (const MatrixFree &matrix_free, - const unsigned int fe_no = 0, - const unsigned int quad_no = 0); - - /** - * Constructor that comes with reduced functionality and works similar as - * FEValues. The user has to provide a structure of type MappingFEEvaluation - * and a DoFHandler in order to allow for reading out the finite element - * data. It uses the data provided by dof_handler.get_fe(). If the element - * is vector-valued, the optional argument allows to specify the index of - * the base element (as long as the element is primitive, non-primitive are - * not supported currently). - * - * With initialization from a FEValues object, no call to a reinit method of - * this class is necessary. Instead, it is enough if the geometry is - * initialized to a given cell iterator. It can also read from or write to - * vectors in the standard way for DoFHandler::active_cell_iterator - * types (which is less efficient with MPI since index translation has to be - * done), but of course only for one cell at a time. Hence, a kernel using - * this method does not vectorize over several elements (which is most - * efficient for vector operations), but only possibly within the element if - * the evaluate/integrate routines are combined (e.g. for computing cell - * matrices). - */ - FEEvaluationQ_DG0 (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component = 0); - - /** - * Copy constructor - */ - FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other); - - /** - * Evaluates the function values, the gradients, and the Hessians of the FE - * function given at the DoF values in the input vector at the quadrature - * points of the unit cell. The function arguments specify which parts shall - * actually be computed. Needs to be called before the functions @p - * get_value(), @p get_gradient() or @p get_laplacian give useful - * information (unless these values have been set manually). - */ - void evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl = false); - - /** - * This function takes the values and/or gradients that are stored on - * quadrature points, tests them by all the basis functions/gradients on the - * cell and performs the cell integration. The two function arguments @p - * integrate_val and @p integrate_grad are used to enable/disable some of - * values or gradients. - */ - void integrate (const bool integrate_val, - const bool integrate_grad); -}; - - /*----------------------- Inline functions ----------------------------------*/ #ifndef DOXYGEN @@ -1946,13 +1572,13 @@ FEEvaluationBase ::FEEvaluationBase (const MatrixFree &data_in, const unsigned int fe_no_in, const unsigned int quad_no_in, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points) : quad_no (quad_no_in), n_fe_components (data_in.get_dof_info(fe_no_in).n_components), - active_fe_index (data_in.get_dof_info(fe_no_in).fe_index_from_dofs_per_cell - (dofs_per_cell * n_fe_components)), + active_fe_index (data_in.get_dof_info(fe_no_in).fe_index_from_degree + (fe_degree)), active_quad_index (data_in.get_mapping_info(). mapping_data_gen[quad_no_in]. quad_index_from_n_q_points(n_q_points)), @@ -2065,7 +1691,7 @@ FEEvaluationBase active_quad_index (other.active_quad_index), matrix_info (other.matrix_info), dof_info (other.dof_info), - mapping_info (other.mapping_info), + mapping_info (other.mapping_info), stored_shape_info (other.stored_shape_info), data (other.data), cartesian_data (other.cartesian_data), @@ -2601,14 +2227,15 @@ FEEvaluationBase } else { - // no constraint at all: loop bounds are known, compiler can - // unroll without checks + // no constraint at all: compiler can unroll at least the + // vectorization loop AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); - for (unsigned int j=0; j::n_array_elements) + for (unsigned int v=0; v::n_array_elements; ++v) + for (unsigned int comp=0; comp } else { - // no constraint at all: loop bounds are known, compiler can - // unroll without checks + // no constraint at all: compiler can unroll at least the + // vectorization loop AssertDimension (dof_info->end_indices(cell)-dof_indices, static_cast(n_local_dofs)); - for (unsigned int j=0; j::n_array_elements) + for (unsigned int v=0; v::n_array_elements; ++v) + operation.process_dof (dof_indices[j+v], *src[0], + local_data[j+v]); } } @@ -3826,11 +3454,11 @@ FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, const unsigned int quad_no_in, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points) : FEEvaluationBase - (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) + (data_in, fe_no, quad_no_in, fe_degree, n_q_points) {} @@ -3867,11 +3495,11 @@ FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, const unsigned int quad_no_in, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points) : FEEvaluationBase - (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) + (data_in, fe_no, quad_no_in, fe_degree, n_q_points) {} @@ -4101,11 +3729,11 @@ FEEvaluationAccess ::FEEvaluationAccess (const MatrixFree &data_in, const unsigned int fe_no, const unsigned int quad_no_in, - const unsigned int dofs_per_cell, + const unsigned int fe_degree, const unsigned int n_q_points) : FEEvaluationBase - (data_in, fe_no, quad_no_in, dofs_per_cell, n_q_points) + (data_in, fe_no, quad_no_in, fe_degree, n_q_points) {} @@ -4448,278 +4076,140 @@ FEEvaluationAccess -/*----------------------- FEEvaluationGeneral -------------------------------*/ - -template -inline -FEEvaluationGeneral -::FEEvaluationGeneral (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) - : - BaseClass (data_in, fe_no, quad_no, dofs_per_cell, n_q_points) +namespace internal { - check_template_arguments(fe_no); - set_data_pointers(); -} + /** + * In this namespace, the evaluator routines that evaluate the tensor + * products are implemented. + */ + enum EvaluatorVariant + { + evaluate_general, + evaluate_symmetric, + evaluate_evenodd + }; + /** + * Generic evaluator framework + */ + template + struct EvaluatorTensorProduct + {}; + /** + * Internal evaluator for 1d-3d shape function using the tensor product form + * of the basis functions + */ + template + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + + /** + * Empty constructor. Does nothing. Be careful when using 'values' and + * related methods because they need to be filled with the other pointer + */ + EvaluatorTensorProduct () + : + shape_values (0), + shape_gradients (0), + shape_hessians (0) + {} + + /** + * Constructor, taking the data from ShapeInfo + */ + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &shape_hessians) + : + shape_values (shape_values.begin()), + shape_gradients (shape_gradients.begin()), + shape_hessians (shape_hessians.begin()) + {} + + template + void + values (const Number in [], + Number out[]) const + { + apply(shape_values, in, out); + } -template -inline -FEEvaluationGeneral -::FEEvaluationGeneral (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) - : - BaseClass (geometry, dof_handler, first_selected_component) -{ - check_template_arguments(numbers::invalid_unsigned_int); - set_data_pointers(); -} + template + void + gradients (const Number in [], + Number out[]) const + { + apply(shape_gradients, in, out); + } + template + void + hessians (const Number in [], + Number out[]) const + { + apply(shape_hessians, in, out); + } + template + static void apply (const Number *shape_data, + const Number in [], + Number out []); -template -inline -FEEvaluationGeneral -::FEEvaluationGeneral (const FEEvaluationGeneral &other) - : - BaseClass (other) -{ - set_data_pointers(); -} + const Number * shape_values; + const Number * shape_gradients; + const Number * shape_hessians; + }; + // evaluates the given shape data in 1d-3d using the tensor product + // form. does not use a particular layout of entries in the matrices + // like the functions below and corresponds to a usual matrix-matrix + // product + template + template + inline + void + EvaluatorTensorProduct + ::apply (const Number *shape_data, + const Number in [], + Number out []) + { + AssertIndexRange (direction, dim); + const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, + nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); + const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); + const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); + const int stride = Utilities::fixed_int_power::value; -template -inline -void -FEEvaluationGeneral -::check_template_arguments(const unsigned int fe_no) -{ -#ifdef DEBUG - // print error message when the dimensions do not match. Propose a possible - // fix - if (fe_degree != this->data->fe_degree - || - n_q_points != this->data->n_q_points) - { - std::string message = - "-------------------------------------------------------\n"; - message += "Illegal arguments in constructor/wrong template arguments!\n"; - message += " Called --> FEEvaluationquad_no) + ")\n"; - - // check whether some other vector component has the correct number of - // points - unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, - proposed_quad_comp = numbers::invalid_unsigned_int; - if (fe_no != numbers::invalid_unsigned_int) - { - if (fe_degree == this->data->fe_degree) - proposed_dof_comp = fe_no; - else - for (unsigned int no=0; nomatrix_info->n_components(); ++no) - if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree - == fe_degree) - { - proposed_dof_comp = no; - break; - } - if (n_q_points == - this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index]) - proposed_quad_comp = this->quad_no; - else - for (unsigned int no=0; nomapping_info->mapping_data_gen.size(); ++no) - if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index] - == n_q_points) - { - proposed_quad_comp = no; - break; - } - } - if (proposed_dof_comp != numbers::invalid_unsigned_int && - proposed_quad_comp != numbers::invalid_unsigned_int) - { - if (proposed_dof_comp != fe_no) - message += "Wrong vector component selection:\n"; - else - message += "Wrong quadrature formula selection:\n"; - message += " Did you mean FEEvaluationquad_no) - correct_pos += " ^\n"; - else - correct_pos += " \n"; - message += " " + correct_pos; - } - // ok, did not find the numbers specified by the template arguments in - // the given list. Suggest correct template arguments - const unsigned int proposed_n_q_points_1d = static_cast(std::pow(1.001*this->data->n_q_points,1./dim)); - message += "Wrong template arguments:\n"; - message += " Did you mean FEEvaluationdata->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"; - std::string correct_pos; - if (this->data->fe_degree != fe_degree) - correct_pos = " ^"; - else - correct_pos = " "; - if (proposed_n_q_points_1d != n_q_points_1d) - correct_pos += " ^\n"; - else - correct_pos += " \n"; - message += " " + correct_pos; - - Assert (fe_degree == this->data->fe_degree && - n_q_points == this->data->n_q_points, - ExcMessage(message)); - } - if (fe_no != numbers::invalid_unsigned_int) - { - AssertDimension (n_q_points, - this->mapping_info->mapping_data_gen[this->quad_no]. - n_q_points[this->active_quad_index]); - AssertDimension (this->data->dofs_per_cell * this->n_fe_components, - this->dof_info->dofs_per_cell[this->active_fe_index]); - } -#endif -} - - - -template -inline -void -FEEvaluationGeneral -::set_data_pointers() -{ - AssertIndexRange(this->data->dofs_per_cell, dofs_per_cell+2); - const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell; - - // set the pointers to the correct position in the data array - for (unsigned int c=0; cvalues_dofs[c] = &my_data_array[c*desired_dofs_per_cell]; - this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points]; - for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+ - n_q_points) - + - (c*dim+d)*n_q_points]; - for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) - this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+ - desired_dofs_per_cell) - + - (c*(dim*dim+dim)+d)*n_q_points]; - } -} - - -template -inline -void -FEEvaluationGeneral -::check_dofs_per_cell (const unsigned int given_dofs_per_cell) const -{ -#ifdef DEBUG - if (given_dofs_per_cell != this->data->dofs_per_cell) - { - std::ostringstream str; - str << "Dofs per cell in FEEvaluation* class does not match the element (" - << given_dofs_per_cell << " != " << this->data->dofs_per_cell << "). " - << "Suggestion: "; - if (this->data->dofs_per_cell == dofs_per_cell) - str << "FEEvaluation/FEEvaluationGeneral"; - else if (this->data->dofs_per_cell == - internal::MatrixFreeFunctions::DGP_dofs_per_cell::value) - str << "FEEvaluationDGP"; - else if (this->data->dofs_per_cell == dofs_per_cell+1) - str << "FEEvaluationQ_DG0"; - else - str << "No matching suggestion found"; - - Assert(given_dofs_per_cell == this->data->dofs_per_cell, - ExcMessage(str.str().c_str())); - } -#endif -} - - - -namespace internal -{ - // evaluates the given shape data in 1d-3d using the tensor product - // form. does not use a particular layout of entries in the matrices - // like the functions below and corresponds to a usual matrix-matrix - // product - template - inline - void - apply_tensor_product (const Number *shape_data, - const Number in [], - Number out []) - { - AssertIndexRange (direction, dim); - const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, - nn = dof_to_quad ? n_q_points_1d : (fe_degree+1); - - const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; - - for (int i2=0; i2 + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + + /** + * Constructor, taking the data from ShapeInfo + */ + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &shape_hessians) + : + shape_values (shape_values.begin()), + shape_gradients (shape_gradients.begin()), + shape_hessians (shape_hessians.begin()) + {} + + template + void + values (const Number in [], + Number out[]) const; + + template + void + gradients (const Number in [], + Number out[]) const; + + template + void + hessians (const Number in [], + Number out[]) const; + + const Number * shape_values; + const Number * shape_gradients; + const Number * shape_hessians; + }; + + + + // In this case, the 1D shape values read (sorted lexicographically, rows + // run over 1D dofs, columns over quadrature points): // Q2 --> [ 0.687 0 -0.087 ] // [ 0.4 1 0.4 ] // [-0.087 0 0.687 ] @@ -4861,13 +4391,13 @@ namespace internal // In these matrices, we want to use avoid computations involving zeros and // ones and in addition use the symmetry in entries to reduce the number of // read operations. - template + template + template inline void - apply_tensor_product_values (const Number *shape_values, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::values (const Number in [], + Number out []) const { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -5041,10 +4571,6 @@ namespace internal - // evaluates the given shape data in 1d-3d using the tensor product - // form assuming the symmetries of unit cell shape gradients for - // finite elements in FEEvaluation - // For the specialized loop used for the gradient computation in // here, the 1D shape values read (sorted lexicographically, rows // run over 1D dofs, columns over quadrature points): @@ -5064,13 +4590,13 @@ namespace internal // In these matrices, we want to use avoid computations involving // zeros and ones and in addition use the symmetry in entries to // reduce the number of read operations. - template + template + template inline void - apply_tensor_product_gradients (const Number *shape_gradients, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::gradients (const Number in [], + Number out []) const { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -5208,13 +4734,13 @@ namespace internal // evaluates the given shape data in 1d-3d using the tensor product // form assuming the symmetries of unit cell shape hessians for // finite elements in FEEvaluation - template + template + template inline void - apply_tensor_product_hessians (const Number *shape_hessians, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::hessians (const Number in [], + Number out []) const { AssertIndexRange (direction, dim); const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d, @@ -5360,7 +4886,7 @@ namespace internal - // This method implements a different approach to the symmetric case for + // This class implements a different approach to the symmetric case for // values, gradients, and Hessians also treated with the above functions: It // is possible to reduce the cost per dimension from N^2 to N^2/2, where N // is the number of 1D dofs (there are only N^2/2 different entries in the @@ -5373,13 +4899,80 @@ namespace internal // experiments in the book say that the method is not efficient for N<20, it // is more efficient in the context where the loop bounds are compile-time // constants (templates). - template + template + struct EvaluatorTensorProduct + { + static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; + static const unsigned int n_q_points = Utilities::fixed_int_power::value; + + /** + * Empty constructor. Does nothing. Be careful when using 'values' and + * related methods because they need to be filled with the other pointer + */ + EvaluatorTensorProduct () + : + shape_values (0), + shape_gradients (0), + shape_hessians (0) + {} + + /** + * Constructor, taking the data from ShapeInfo (using the even-odd + * variants stored there) + */ + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &shape_hessians) + : + shape_values (shape_values.begin()), + shape_gradients (shape_gradients.begin()), + shape_hessians (shape_hessians.begin()) + {} + + template + void + values (const Number in [], + Number out[]) const + { + apply(shape_values, in, out); + } + + template + void + gradients (const Number in [], + Number out[]) const + { + apply(shape_gradients, in, out); + } + + template + void + hessians (const Number in [], + Number out[]) const + { + apply(shape_hessians, in, out); + } + + template + static void apply (const Number *shape_data, + const Number in [], + Number out []); + + const Number * shape_values; + const Number * shape_gradients; + const Number * shape_hessians; + }; + + + + template + template inline void - apply_tensor_product_evenodd (const Number shapes [][(n_q_points_1d+1)/2], - const Number in [], - Number out []) + EvaluatorTensorProduct + ::apply (const Number *shapes, + const Number in [], + Number out []) { AssertIndexRange (type, 3); AssertIndexRange (direction, dim); @@ -5392,6 +4985,8 @@ namespace internal const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); const int stride = Utilities::fixed_int_power::value; + const int offset = (n_q_points_1d+1)/2; + // this code may look very inefficient at first sight due to the many // different cases with if's at the innermost loop part, but all of the // conditionals can be evaluated at compile time because they are @@ -5421,25 +5016,25 @@ namespace internal { if (dof_to_quad == true) { - r0 = shapes[0][col] * xp[0]; - r1 = shapes[fe_degree][col] * xm[0]; + r0 = shapes[col] * xp[0]; + r1 = shapes[fe_degree*offset + col] * xm[0]; } else { - r0 = shapes[col][0] * xp[0]; - r1 = shapes[fe_degree-col][0] * xm[0]; + r0 = shapes[col*offset] * xp[0]; + r1 = shapes[(fe_degree-col)*offset] * xm[0]; } for (int ind=1; ind 0)) - r0 += shapes[col][mid] * in[stride*mid]; + r0 += shapes[col*offset+mid] * in[stride*mid]; if (add == false) { @@ -5484,14 +5079,14 @@ namespace internal Number r0; if (mid > 0) { - r0 = shapes[0][n_cols] * xp[0]; + r0 = shapes[n_cols] * xp[0]; for (int ind=1; ind [-3 -1 1 ] - // [ 4 0 -4 ] - // [-1 1 3 ] - // Q3 --> [-6 -1.618 0.618 -1 ] - // [ 8.09 0 -2.236 3.09 ] - // [-3.09 2.236 0 -8.09 ] - // [ 1 -0.618 1.618 6 ] - // Q4 --> [-10 -2.482 0.75 -0.518 1 ] - // [ 13.51 0 -2.673 1.528 -2.82 ] - // [-5.333 3.491 0 -3.491 5.333 ] - // [ 2.82 -1.528 2.673 0 -13.51 ] - // [-1 0.518 -0.75 2.482 10 ] - template + // Select evaluator type from element shape function type + template + struct EvaluatorSelector {}; + + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_general; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_symmetric; + }; + + template <> struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_evenodd; + }; + + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_general; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_symmetric; + }; + + template <> + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_evenodd; + }; + + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_evenodd; + }; + + + + // This struct performs the evaluation of function values, gradients and + // Hessians for tensor-product finite elements. The operation is used for + // both the symmetric and non-symmetric case, which use different apply + // functions 'values', 'gradients' in the individual coordinate + // directions. The apply functions for values are provided through one of + // the template classes EvaluatorTensorProduct which in turn are selected + // from the MatrixFreeFunctions::ElementType template argument. + // + // There is a specialization made for Gauss-Lobatto elements further down + // where the 'values' operation is identity, which allows us to write + // shorter code. + template + struct FEEvaluationImpl + { + static + void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl); + + static + void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + const bool evaluate_val, + const bool evaluate_grad); + }; + + + template inline void - apply_tensor_product_gradients_gl (const Number *shape_gradients, - const Number in [], - Number out []) + FEEvaluationImpl + ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) { - AssertIndexRange (direction, dim); - const int mm = fe_degree+1; - const int nn = fe_degree+1; - const int n_cols = nn / 2; - const int mid = mm / 2; + const EvaluatorVariant variant = + EvaluatorSelector4)>::variant; + typedef EvaluatorTensorProduct > Eval; + Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd : + shape_info.shape_values, + variant == evaluate_evenodd ? shape_info.shape_gra_evenodd : + shape_info.shape_gradients, + variant == evaluate_evenodd ? shape_info.shape_hes_evenodd : + shape_info.shape_hessians); + + const unsigned int temp_size = Eval::dofs_per_cell > Eval::n_q_points ? + Eval::dofs_per_cell : Eval::n_q_points; + + VectorizedArray **values_dofs = values_dofs_actual; + VectorizedArray data_array[type!=MatrixFreeFunctions::truncated_tensor ? 1 : + n_components*Eval::dofs_per_cell]; + VectorizedArray *expanded_dof_values[n_components]; + if (type == MatrixFreeFunctions::truncated_tensor) + { + for (unsigned int c=0; c 1 ? (direction > 0 ? nn : mm) : 1); - const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1); - const int stride = Utilities::fixed_int_power::value; + unsigned int count_p = 0, count_q = 0; + for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) + { + for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j) + { + for (unsigned int k=0; k(); + } + for (unsigned int j=fe_degree+1-i; j(); + } + AssertDimension(count_q, Eval::dofs_per_cell); + } - for (int i2=0; i21?1:0; + const unsigned int d2 = dim>2?2:0; + const unsigned int d3 = dim>2?3:0; + const unsigned int d4 = dim>2?4:0; + const unsigned int d5 = dim>2?5:0; + + switch (dim) { - for (int i1=0; i1 (values_dofs[c], values_quad[c]); + if (evaluate_grad == true) + eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]); + if (evaluate_lapl == true) + eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]); + } + break; + + case 2: + for (unsigned int c=0; c temp1[temp_size]; + VectorizedArray temp2[temp_size]; + + // grad x + if (evaluate_grad == true) { - Number val0, val1, in0, in1, res0, res1; - if (mid > 0) + eval.template gradients<0,true,false> (values_dofs[c], temp1); + eval.template values<1,true,false> (temp1, gradients_quad[c][0]); + } + if (evaluate_lapl == true) + { + // grad xy + if (evaluate_grad == false) + eval.template gradients<0,true,false>(values_dofs[c], temp1); + eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]); + + // grad xx + eval.template hessians<0,true,false>(values_dofs[c], temp1); + eval.template values<1,true,false> (temp1, hessians_quad[c][0]); + } + + // grad y + eval.template values<0,true,false> (values_dofs[c], temp1); + if (evaluate_grad == true) + eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]); + + // grad yy + if (evaluate_lapl == true) + eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]); + + // val: can use values applied in x + if (evaluate_val == true) + eval.template values<1,true,false> (temp1, values_quad[c]); + } + break; + + case 3: + for (unsigned int c=0; c temp1[temp_size]; + VectorizedArray temp2[temp_size]; + + if (evaluate_grad == true) + { + // grad x + eval.template gradients<0,true,false> (values_dofs[c], temp1); + eval.template values<1,true,false> (temp1, temp2); + eval.template values<2,true,false> (temp2, gradients_quad[c][0]); + } + + if (evaluate_lapl == true) + { + // grad xz + if (evaluate_grad == false) { - if (dof_to_quad == true) - { - val0 = shape_gradients[col]; - val1 = shape_gradients[nn-1-col]; - } - else - { - val0 = shape_gradients[col*mm]; - val1 = shape_gradients[(nn-col-1)*mm]; - } - in0 = in[0]; - in1 = in[stride*(mm-1)]; - if (col == 0) - { - if ((mm+dof_to_quad)%2 == 1) - { - res0 = val0 * in0; - res1 = -in0; - res0 += in1; - res1 -= val0 * in1; - } - else - { - res0 = val0 * in0; - res0 -= in1; - res1 = in0; - res1 -= val0 * in1; - } - } - else - { - res0 = val0 * in0; - res1 = val1 * in0; - res0 -= val1 * in1; - res1 -= val0 * in1; - } - for (int ind=1; ind (values_dofs[c], temp1); + eval.template values<1,true,false> (temp1, temp2); } + eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]); + + // grad xy + eval.template gradients<1,true,false> (temp1, temp2); + eval.template values<2,true,false> (temp2, hessians_quad[c][d3]); + + // grad xx + eval.template hessians<0,true,false>(values_dofs[c], temp1); + eval.template values<1,true,false> (temp1, temp2); + eval.template values<2,true,false> (temp2, hessians_quad[c][0]); } - if ( nn%2 == 1 ) + + // grad y + eval.template values<0,true,false> (values_dofs[c], temp1); + if (evaluate_grad == true) { - Number val0, res0; - if (dof_to_quad == true) - val0 = shape_gradients[n_cols]; - else - val0 = shape_gradients[n_cols*mm]; - if (mid > 0) - { - res0 = in[0] - in[stride*(mm-1)]; - res0 *= val0; - for (int ind=1; ind(temp1, temp2); + eval.template values<2,true,false> (temp2, gradients_quad[c][d1]); } - // increment: in regular case, just go to the next point in - // x-direction. for y-part in 3D and if we are at the end of one - // chunk in x-dir, need to jump over to the next layer in - // z-direction - switch (direction) + if (evaluate_lapl == true) { - case 0: - in += mm; - out += nn; - break; - case 1: - case 2: - ++in; - ++out; - break; - default: - Assert (false, ExcNotImplemented()); + // grad yz + if (evaluate_grad == false) + eval.template gradients<1,true,false>(temp1, temp2); + eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]); + + // grad yy + eval.template hessians<1,true,false> (temp1, temp2); + eval.template values<2,true,false> (temp2, hessians_quad[c][d1]); } - } - - if (direction == 1) - { - in += nn * (mm-1); - out += nn * (nn-1); - } - } - } + // grad z: can use the values applied in x direction stored in temp1 + eval.template values<1,true,false> (temp1, temp2); + if (evaluate_grad == true) + eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]); + // grad zz: can use the values applied in x and y direction stored + // in temp2 + if (evaluate_lapl == true) + eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]); - // This performs the evaluation of function values, gradients and Hessians - // for tensor-product finite elements. The operation is used for both - // FEEvaluationGeneral and FEEvaluation, which provide different functions - // apply_values, apply_gradients in the individual coordinate directions. - // use different versions for 1d, 2d, 3d to avoid nasty compiler warnings - template - inline - void - do_evaluate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][1], - VectorizedArray* hessians_quad[][1], - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl, - internal::int2type<1>) - { - const unsigned int n_components = FEEval::n_components; + // val: can use the values applied in x & y direction stored in temp2 + if (evaluate_val == true) + eval.template values<2,true,false> (temp2, values_quad[c]); + } + break; - for (unsigned int c=0; c - (values_dofs[c], fe_eval.values_quad[c]); - if (evaluate_grad == true) - fe_eval.template apply_gradients<0,true,false> - (values_dofs[c], fe_eval.gradients_quad[c][0]); - if (evaluate_lapl == true) - fe_eval.template apply_hessians<0,true,false> - (values_dofs[c], fe_eval.hessians_quad[c][0]); + default: + AssertThrow(false, ExcNotImplemented()); } + + // case additional dof for FE_Q_DG0: add values; gradients and second + // derivatives evaluate to zero + if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_val) + for (unsigned int c=0; c + + template inline void - do_evaluate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][2], - VectorizedArray* hessians_quad[][3], - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl, - internal::int2type<2>) + FEEvaluationImpl + ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs_actual[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + const bool integrate_val, + const bool integrate_grad) { - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? - FEEval::dofs_per_cell : FEEval::n_q_points; - const unsigned int n_components = FEEval::n_components; - - for (unsigned int c=0; c4)>::variant; + typedef EvaluatorTensorProduct > Eval; + Eval eval (variant == evaluate_evenodd ? shape_info.shape_val_evenodd : + shape_info.shape_values, + variant == evaluate_evenodd ? shape_info.shape_gra_evenodd : + shape_info.shape_gradients, + variant == evaluate_evenodd ? shape_info.shape_hes_evenodd : + shape_info.shape_hessians); + + const unsigned int temp_size = Eval::dofs_per_cell > Eval::n_q_points ? + Eval::dofs_per_cell : Eval::n_q_points; + VectorizedArray temp1[temp_size]; + VectorizedArray temp2[temp_size]; + + // expand dof_values to tensor product for truncated tensor products + VectorizedArray ** values_dofs = values_dofs_actual; + VectorizedArray data_array[type!=MatrixFreeFunctions::truncated_tensor ? 1 : + n_components*Eval::dofs_per_cell]; + VectorizedArray *expanded_dof_values[n_components]; + if (type == MatrixFreeFunctions::truncated_tensor) { - VectorizedArray temp1[temp_size]; - VectorizedArray temp2[temp_size]; - - // grad x - if (evaluate_grad == true) - { - fe_eval.template apply_gradients<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, gradients_quad[c][0]); - } - if (evaluate_lapl == true) - { - // grad xy - if (evaluate_grad == false) - fe_eval.template apply_gradients<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_gradients<1,true,false> - (temp1, hessians_quad[c][2]); - - // grad xx - fe_eval.template apply_hessians<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, hessians_quad[c][0]); - } - - // grad y - fe_eval.template apply_values<0,true,false> - (values_dofs[c], temp1); - if (evaluate_grad == true) - fe_eval.template apply_gradients<1,true,false> - (temp1, gradients_quad[c][1]); - - // grad yy - if (evaluate_lapl == true) - fe_eval.template apply_hessians<1,true,false> - (temp1, hessians_quad[c][1]); - - // val: can use values applied in x - if (evaluate_val == true) - fe_eval.template apply_values<1,true,false> - (temp1, values_quad[c]); + for (unsigned int c=0; c - inline - void - do_evaluate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][3], - VectorizedArray* hessians_quad[][6], - const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl, - internal::int2type<3>) - { - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? - FEEval::dofs_per_cell : FEEval::n_q_points; - const unsigned int n_components = FEEval::n_components; + // These avoid compiler errors; they are only used in sensible context but + // compilers typically cannot detect when we access something like + // gradients_quad[2] only for dim==3. + const unsigned int d1 = dim>1?1:0; + const unsigned int d2 = dim>2?2:0; - for (unsigned int c=0; c temp1[temp_size]; - VectorizedArray temp2[temp_size]; - - if (evaluate_grad == true) - { - // grad x - fe_eval.template apply_gradients<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, gradients_quad[c][0]); - } - - if (evaluate_lapl == true) + case 1: + for (unsigned int c=0; c (values_quad[c], values_dofs[c]); + if (integrate_grad == true) { - fe_eval.template apply_gradients<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); + if (integrate_val == true) + eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]); + else + eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]); } - fe_eval.template apply_gradients<2,true,false> - (temp2, hessians_quad[c][4]); - - // grad xy - fe_eval.template apply_gradients<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, hessians_quad[c][3]); - - // grad xx - fe_eval.template apply_hessians<0,true,false> - (values_dofs[c], temp1); - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, hessians_quad[c][0]); } + break; - // grad y - fe_eval.template apply_values<0,true,false> - (values_dofs[c], temp1); - if (evaluate_grad == true) + case 2: + for (unsigned int c=0; c - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, gradients_quad[c][1]); + if (integrate_val == true) + { + // val + eval.template values<0,false,false> (values_quad[c], temp1); + //grad x + if (integrate_grad == true) + eval.template gradients<0,false,true> (gradients_quad[c][0], temp1); + eval.template values<1,false,false>(temp1, values_dofs[c]); + } + if (integrate_grad == true) + { + // grad y + eval.template values<0,false,false> (gradients_quad[c][d1], temp1); + if (integrate_val == false) + { + eval.template gradients<1,false,false>(temp1, values_dofs[c]); + //grad x + eval.template gradients<0,false,false> (gradients_quad[c][0], temp1); + eval.template values<1,false,true> (temp1, values_dofs[c]); + } + else + eval.template gradients<1,false,true>(temp1, values_dofs[c]); + } } + break; - if (evaluate_lapl == true) + case 3: + for (unsigned int c=0; c - (temp1, temp2); - fe_eval.template apply_gradients<2,true,false> - (temp2, hessians_quad[c][5]); - - // grad yy - fe_eval.template apply_hessians<1,true,false> - (temp1, temp2); - fe_eval.template apply_values<2,true,false> - (temp2, hessians_quad[c][1]); + if (integrate_val == true) + { + // val + eval.template values<0,false,false> (values_quad[c], temp1); + //grad x: can sum to temporary value in temp1 + if (integrate_grad == true) + eval.template gradients<0,false,true> (gradients_quad[c][0], temp1); + eval.template values<1,false,false>(temp1, temp2); + eval.template values<0,false,false> (gradients_quad[c][d1], temp1); + if (integrate_grad == true) + eval.template gradients<1,false,true>(temp1, temp2); + eval.template values<2,false,false> (temp2, values_dofs[c]); + } + else if (integrate_grad == true) + { + eval.template gradients<0,false,false>(gradients_quad[c][0], temp1); + eval.template values<1,false,false> (temp1, temp2); + eval.template values<0,false,false> (gradients_quad[c][d1], temp1); + eval.template gradients<1,false,true>(temp1, temp2); + eval.template values<2,false,false> (temp2, values_dofs[c]); + } + if (integrate_grad == true) + { + // grad z: can sum to temporary x and y value in output + eval.template values<0,false,false> (gradients_quad[c][d2], temp1); + eval.template values<1,false,false> (temp1, temp2); + eval.template gradients<2,false,true> (temp2, values_dofs[c]); + } } + break; - // grad z: can use the values applied in x direction stored in temp1 - fe_eval.template apply_values<1,true,false> - (temp1, temp2); - if (evaluate_grad == true) - fe_eval.template apply_gradients<2,true,false> - (temp2, gradients_quad[c][2]); - - // grad zz: can use the values applied in x and y direction stored - // in temp2 - if (evaluate_lapl == true) - fe_eval.template apply_hessians<2,true,false> - (temp2, hessians_quad[c][2]); - - // val: can use the values applied in x & y direction stored in temp2 - if (evaluate_val == true) - fe_eval.template apply_values<2,true,false> - (temp2, values_quad[c]); + default: + AssertThrow(false, ExcNotImplemented()); } - } + // case FE_Q_DG0: add values, gradients and second derivatives are zero + if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0) + { + if (integrate_val) + for (unsigned int c=0; c(); + } - - template - inline - void - do_integrate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][1], - const bool integrate_val, - const bool integrate_grad, - internal::int2type<1>) - { - const unsigned int n_components = FEEval::n_components; - - for (unsigned int c=0; c - (gradients_quad[c][0], fe_eval.values_dofs[c]); - if (integrate_val == true) + unsigned int count_p = 0, count_q = 0; + for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) { - if (integrate_grad == true) - fe_eval.template apply_values<0,false,true> - (values_quad[c], values_dofs[c]); - else - fe_eval.template apply_values<0,false,false> - (values_quad[c], values_dofs[c]); + for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j) + { + for (unsigned int k=0; k + // This a specialization for Gauss-Lobatto elements where the 'values' + // operation is identity, which allows us to write shorter code. + template + struct FEEvaluationImpl + { + static + void evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl); + + static + void integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + const bool integrate_val, + const bool integrate_grad); + }; + + template inline void - do_integrate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][2], - const bool integrate_val, - const bool integrate_grad, - internal::int2type<2>) + FEEvaluationImpl + ::evaluate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + VectorizedArray *hessians_quad[][(dim*(dim+1))/2], + const bool evaluate_val, + const bool evaluate_grad, + const bool evaluate_lapl) { - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? - FEEval::dofs_per_cell : FEEval::n_q_points; - const unsigned int n_components = FEEval::n_components; - - for (unsigned int c=0; c > Eval; + Eval eval (shape_info.shape_val_evenodd, shape_info.shape_gra_evenodd, + shape_info.shape_hes_evenodd); + + // These avoid compiler errors; they are only used in sensible context but + // compilers typically cannot detect when we access something like + // gradients_quad[2] only for dim==3. + const unsigned int d1 = dim>1?1:0; + const unsigned int d2 = dim>2?2:0; + const unsigned int d3 = dim>2?3:0; + const unsigned int d4 = dim>2?4:0; + const unsigned int d5 = dim>2?5:0; + + switch (dim) { - VectorizedArray temp1[temp_size]; - VectorizedArray temp2[temp_size]; + case 1: + if (evaluate_val == true) + std::memcpy (values_quad[0], values_dofs[0], + eval.dofs_per_cell * n_components * + sizeof (values_dofs[0][0])); + for (unsigned int c=0; c(values_dofs[c], gradients_quad[c][0]); + if (evaluate_lapl == true) + eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]); + } + break; - // val - if (integrate_val == true) - fe_eval.template apply_values<0,false,false> - (values_quad[c], temp1); - if (integrate_grad == true) + case 2: + if (evaluate_val == true) { - //grad x - if (integrate_val == true) - fe_eval.template apply_gradients<0,false,true> - (gradients_quad[c][0], temp1); - else - fe_eval.template apply_gradients<0,false,false> - (gradients_quad[c][0], temp1); + std::memcpy (values_quad[0], values_dofs[0], + Eval::dofs_per_cell * n_components * + sizeof (values_dofs[0][0])); } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<1,false,false> - (temp1, values_dofs[c]); - if (integrate_grad == true) + if (evaluate_grad == true) + for (unsigned int comp=0; comp (values_dofs[comp], + gradients_quad[comp][0]); + // grad y + eval.template gradients<1,true,false> (values_dofs[comp], + gradients_quad[comp][d1]); + } + if (evaluate_lapl == true) + for (unsigned int comp=0; comp (values_dofs[comp], + hessians_quad[comp][0]); + // hess y + eval.template hessians<1,true,false> (values_dofs[comp], + hessians_quad[comp][d1]); + + VectorizedArray temp1[Eval::dofs_per_cell]; + // grad x grad y + eval.template gradients<0,true,false> (values_dofs[comp], temp1); + eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d1+d1]); + } + break; + + case 3: + if (evaluate_val == true) { - // grad y - fe_eval.template apply_values<0,false,false> - (gradients_quad[c][1], temp1); - fe_eval.template apply_gradients<1,false,true> - (temp1, values_dofs[c]); + std::memcpy (values_quad[0], values_dofs[0], + Eval::dofs_per_cell * n_components * + sizeof (values_dofs[0][0])); } + if (evaluate_grad == true) + for (unsigned int comp=0; comp (values_dofs[comp], + gradients_quad[comp][0]); + // grad y + eval.template gradients<1,true,false> (values_dofs[comp], + gradients_quad[comp][d1]); + // grad y + eval.template gradients<2,true,false> (values_dofs[comp], + gradients_quad[comp][d2]); + } + if (evaluate_lapl == true) + for (unsigned int comp=0; comp (values_dofs[comp], + hessians_quad[comp][0]); + // grad y + eval.template hessians<1,true,false> (values_dofs[comp], + hessians_quad[comp][d1]); + // grad y + eval.template hessians<2,true,false> (values_dofs[comp], + hessians_quad[comp][d2]); + + VectorizedArray temp1[Eval::dofs_per_cell]; + // grad xy + eval.template gradients<0,true,false> (values_dofs[comp], temp1); + eval.template gradients<1,true,false> (temp1, hessians_quad[comp][d3]); + // grad xz + eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d4]); + // grad yz + eval.template gradients<1,true,false> (values_dofs[comp], temp1); + eval.template gradients<2,true,false> (temp1, hessians_quad[comp][d5]); + } + break; + default: + AssertThrow(false, ExcNotImplemented()); } } - template + template inline void - do_integrate (FEEval &fe_eval, - VectorizedArray* values_dofs[], - VectorizedArray* values_quad[], - VectorizedArray* gradients_quad[][3], - const bool integrate_val, - const bool integrate_grad, - internal::int2type<3>) + FEEvaluationImpl + ::integrate (const MatrixFreeFunctions::ShapeInfo &shape_info, + VectorizedArray *values_dofs[], + VectorizedArray *values_quad[], + VectorizedArray *gradients_quad[][dim], + const bool integrate_val, + const bool integrate_grad) { - const unsigned int temp_size = FEEval::dofs_per_cell > FEEval::n_q_points ? - FEEval::dofs_per_cell : FEEval::n_q_points; - const unsigned int n_components = FEEval::n_components; - - for (unsigned int c=0; c > Eval; + Eval eval (shape_info.shape_val_evenodd, shape_info.shape_gra_evenodd, + shape_info.shape_hes_evenodd); + + // These avoid compiler errors; they are only used in sensible context but + // compilers typically cannot detect when we access something like + // gradients_quad[2] only for dim==3. + const unsigned int d1 = dim>1?1:0; + const unsigned int d2 = dim>2?2:0; + + if (integrate_val == true) + std::memcpy (values_dofs[0], values_quad[0], + Eval::dofs_per_cell * n_components * + sizeof (values_dofs[0][0])); + switch (dim) { - VectorizedArray temp1[temp_size]; - VectorizedArray temp2[temp_size]; - - if (integrate_val == true) - { - // val - fe_eval.template apply_values<0,false,false> - (values_quad[c], temp1); - } - if (integrate_grad == true) - { - // grad x: can sum to temporary value in temp1 - if (integrate_val == true) - fe_eval.template apply_gradients<0,false,true> - (gradients_quad[c][0], temp1); - else - fe_eval.template apply_gradients<0,false,false> - (gradients_quad[c][0], temp1); - } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<1,false,false> - (temp1, temp2); - if (integrate_grad == true) - { - // grad y: can sum to temporary x value in temp2 - fe_eval.template apply_values<0,false,false> - (gradients_quad[c][1], temp1); - fe_eval.template apply_gradients<1,false,true> - (temp1, temp2); - } - if (integrate_val || integrate_grad) - fe_eval.template apply_values<2,false,false> - (temp2, values_dofs[c]); - if (integrate_grad == true) + case 1: + for (unsigned int c=0; c - (gradients_quad[c][2], temp1); - fe_eval.template apply_values<1,false,false> - (temp1, temp2); - fe_eval.template apply_gradients<2,false,true> - (temp2, values_dofs[c]); + if (integrate_grad == true) + { + if (integrate_val == true) + eval.template gradients<0,false,true> (gradients_quad[c][0], + values_dofs[c]); + else + eval.template gradients<0,false,false> (gradients_quad[c][0], + values_dofs[c]); + } } - } - } - -} // end of namespace internal - + break; + case 2: + if (integrate_grad == true) + for (unsigned int comp=0; comp (gradients_quad[comp][0], + values_dofs[comp]); + else + eval.template gradients<0, false, false> (gradients_quad[comp][0], + values_dofs[comp]); -template -inline -void -FEEvaluationGeneral -::evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) -{ - this->check_dofs_per_cell(dofs_per_cell); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - internal::do_evaluate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, this->hessians_quad, - evaluate_val, evaluate_grad, evaluate_lapl, - internal::int2type()); + // grad y + eval.template gradients<1, false, true> (gradients_quad[comp][d1], + values_dofs[comp]); + } + break; -#ifdef DEBUG - if (evaluate_val == true) - this->values_quad_initialized = true; - if (evaluate_grad == true) - this->gradients_quad_initialized = true; - if (evaluate_lapl == true) - this->hessians_quad_initialized = true; -#endif -} - - - -template -inline -void -FEEvaluationGeneral -::integrate (const bool integrate_val, - const bool integrate_grad) -{ - this->check_dofs_per_cell(dofs_per_cell); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - - internal::do_integrate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, integrate_val, integrate_grad, - internal::int2type()); - -#ifdef DEBUG - this->dof_values_initialized = true; -#endif -} - - - -template -inline -Point > -FEEvaluationGeneral -::quadrature_point (const unsigned int q) const -{ - Assert (this->mapping_info->quadrature_points_initialized == true, - ExcNotInitialized()); - AssertIndexRange (q, n_q_points); - - // Cartesian mesh: not all quadrature points are stored, only the - // diagonal. Hence, need to find the tensor product index and retrieve the - // value from that - if (this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - Point > point (false); - switch (dim) - { - case 1: - return this->quadrature_points[q]; - case 2: - point[0] = this->quadrature_points[q%n_q_points_1d][0]; - point[1] = this->quadrature_points[q/n_q_points_1d][1]; - return point; - case 3: - point[0] = this->quadrature_points[q%n_q_points_1d][0]; - point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1]; - point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2]; - return point; - default: - Assert (false, ExcNotImplemented()); - return point; - } - } - // all other cases: just return the respective data as it is fully stored - else - return this->quadrature_points[q]; -} - - - -template -template -inline -void -FEEvaluationGeneral -::apply_values(const VectorizedArray in [], - VectorizedArray out []) -{ - internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data->shape_values.begin(), in, out); -} - - - -template -template -inline -void -FEEvaluationGeneral -::apply_gradients(const VectorizedArray in [], - VectorizedArray out []) -{ - internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data->shape_gradients.begin(), in, out); -} - - - -template -template -inline -void -FEEvaluationGeneral -::apply_hessians(const VectorizedArray in [], - VectorizedArray out []) -{ - internal::apply_tensor_product, direction, dof_to_quad, add> - (this->data->shape_hessians.begin(), in, out); -} - - -/*-------------------------- FEEvaluation -----------------------------------*/ - - -template -inline -FEEvaluation -::FEEvaluation (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) - : - BaseClass (data_in, fe_no, quad_no) -{ - compute_even_odd_factors(); -} - - - -template -inline -FEEvaluation -::FEEvaluation (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) - : - BaseClass (geometry, dof_handler, first_selected_component) -{ - compute_even_odd_factors(); -} - - - -template -inline -FEEvaluation -::FEEvaluation (const FEEvaluation &other) - : - BaseClass (other) -{ - compute_even_odd_factors(); -} - - - -template -inline -void -FEEvaluation -::compute_even_odd_factors() -{ - // check whether element is appropriate -#ifdef DEBUG - const double zero_tol = - types_are_equal::value==true?1e-8:1e-7; - std::string error_message = "FEEvaluation not appropriate.\n"; - error_message += " It assumes symmetry of quadrature points w.r.t. 0.5 \n"; - error_message += " and the basis functions starting from left and right.\n"; - error_message += "Try FEEvaluationGeneral<...> instead!"; - - // symmetry for values - const unsigned int n_dofs_1d = fe_degree + 1; - for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) - for (unsigned int j=0; jdata->shape_values[i*n_q_points_1d+j][0] - - this->data->shape_values[(n_dofs_1d-i)*n_q_points_1d - -j-1][0]) < zero_tol, - ExcMessage(error_message)); - - // shape values should be zero at for all basis functions except for one - // where they are one in the middle - if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1) - { - for (int i=0; i(n_dofs_1d/2); ++i) - Assert (std::fabs(this->data->shape_values[i*n_q_points_1d+ - n_q_points_1d/2][0]) < zero_tol, - ExcMessage(error_message)); - Assert (std::fabs(this->data->shape_values[(n_dofs_1d/2)*n_q_points_1d+ - n_q_points_1d/2][0]-1.)< zero_tol, - ExcMessage(error_message)); - } - - // skew-symmetry for gradient, zero of middle basis function in middle - // quadrature point - for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) - for (unsigned int j=0; jdata->shape_gradients[i*n_q_points_1d+j][0] + - this->data->shape_gradients[(n_dofs_1d-i)*n_q_points_1d- - j-1][0]) < zero_tol, - ExcMessage(error_message)); - if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1) - Assert (std::fabs(this->data->shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ - (n_q_points_1d/2)][0]) < zero_tol, - ExcMessage(error_message)); - - - // symmetry for Laplacian - for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) - for (unsigned int j=0; jdata->shape_hessians[i*n_q_points_1d+j][0] - - this->data->shape_hessians[(n_dofs_1d-i)*n_q_points_1d- - j-1][0]) < zero_tol, - ExcMessage(error_message)); -#endif - - // Compute symmetric and skew-symmetric part of shape values for even-odd - // decomposition - for (int i=0; i<(fe_degree+1)/2; ++i) - for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) - { - shape_val_evenodd[i][q] = - 0.5 * (this->data->shape_values[i*n_q_points_1d+q] + - this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); - shape_val_evenodd[fe_degree-i][q] = - 0.5 * (this->data->shape_values[i*n_q_points_1d+q] - - this->data->shape_values[i*n_q_points_1d+n_q_points_1d-1-q]); - - shape_gra_evenodd[i][q] = - 0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] + - this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); - shape_gra_evenodd[fe_degree-i][q] = - 0.5 * (this->data->shape_gradients[i*n_q_points_1d+q] - - this->data->shape_gradients[i*n_q_points_1d+n_q_points_1d-1-q]); - - shape_hes_evenodd[i][q] = - 0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] + - this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); - shape_hes_evenodd[fe_degree-i][q] = - 0.5 * (this->data->shape_hessians[i*n_q_points_1d+q] - - this->data->shape_hessians[i*n_q_points_1d+n_q_points_1d-1-q]); - } - if (fe_degree % 2 == 0) - for (unsigned int q=0; q<(n_q_points_1d+1)/2; ++q) - { - shape_val_evenodd[fe_degree/2][q] = - this->data->shape_values[(fe_degree/2)*n_q_points_1d+q]; - shape_gra_evenodd[fe_degree/2][q] = - this->data->shape_gradients[(fe_degree/2)*n_q_points_1d+q]; - shape_hes_evenodd[fe_degree/2][q] = - this->data->shape_hessians[(fe_degree/2)*n_q_points_1d+q]; - } -} - - - -template -inline -void -FEEvaluation -::evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) -{ - this->check_dofs_per_cell(dofs_per_cell); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - internal::do_evaluate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, this->hessians_quad, - evaluate_val, evaluate_grad, evaluate_lapl, - internal::int2type()); - -#ifdef DEBUG - if (evaluate_val == true) - this->values_quad_initialized = true; - if (evaluate_grad == true) - this->gradients_quad_initialized = true; - if (evaluate_lapl == true) - this->hessians_quad_initialized = true; -#endif -} - - - -template -inline -void -FEEvaluation -::integrate (bool integrate_val,bool integrate_grad) -{ - this->check_dofs_per_cell(dofs_per_cell); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - - internal::do_integrate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, integrate_val, integrate_grad, - internal::int2type()); - -#ifdef DEBUG - this->dof_values_initialized = true; -#endif -} - - - -template -template -inline -void -FEEvaluation -::apply_values (const VectorizedArray in [], - VectorizedArray out []) -{ - // for linear elements, the even-odd decomposition is slower than the plain - // evaluation (additions to create the symmetric and anti-symmetric part), - // for all other orders, we choose even-odd - if (fe_degree > 1 || n_q_points_1d > 3) - internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 0> - (&shape_val_evenodd[0], in, out); - else - internal::apply_tensor_product_values, direction, dof_to_quad, add> - (this->data->shape_values.begin(), in, out); -} - - - -template -template -inline -void -FEEvaluation -::apply_gradients (const VectorizedArray in [], - VectorizedArray out []) -{ - if (fe_degree > 1 || n_q_points_1d > 3) - internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 1> - (&shape_gra_evenodd[0], in, out); - else - internal::apply_tensor_product_gradients, direction, dof_to_quad, add> - (this->data->shape_gradients.begin(), in, out); -} - - - -// Laplacian operator application. Very similar to value application because -// the same symmetry relations hold. However, it is not possible to omit some -// values that are zero for the values -template -template -inline -void -FEEvaluation -::apply_hessians (const VectorizedArray in [], - VectorizedArray out []) -{ - if (fe_degree > 1 || n_q_points_1d > 3) - internal::apply_tensor_product_evenodd, direction, dof_to_quad, add, 2> - (&shape_hes_evenodd[0], in, out); - else - internal::apply_tensor_product_hessians, direction, dof_to_quad, add> - (this->data->shape_hessians.begin(), in, out); -} - - - -/*------------------------- FEEvaluationGL ----------------------------------*/ - -template -inline -FEEvaluationGL -::FEEvaluationGL (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) - : - BaseClass (data_in, fe_no, quad_no) -{ -#ifdef DEBUG - std::string error_mess = "FEEvaluationGL not appropriate. It assumes:\n"; - error_mess += " - identity operation for shape values\n"; - error_mess += " - zero diagonal at interior points for gradients\n"; - error_mess += " - gradient equal to unity at element boundary\n"; - error_mess += "Try FEEvaluation<...> instead!"; - - const double zero_tol = - types_are_equal::value==true?1e-9:1e-7; - - const unsigned int n_points_1d = fe_degree+1; - for (unsigned int i=0; idata->shape_values[i*n_points_1d+j][0])data->shape_values[i*n_points_1d+ - j][0]-1.)data->shape_gradients[i*n_points_1d+i][0])data->shape_gradients[n_points_1d-1][0]- - (n_points_1d%2==0 ? -1. : 1.)) < zero_tol, - ExcMessage (error_mess.c_str())); -#endif -} - - - -template -inline -FEEvaluationGL -::FEEvaluationGL (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) - : - BaseClass (geometry, dof_handler, first_selected_component) -{} - - - -template -inline -FEEvaluationGL -::FEEvaluationGL (const FEEvaluationGL &other) - : - BaseClass (other) -{} - - - -template -inline -void -FEEvaluationGL -::evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) -{ - this->check_dofs_per_cell(dofs_per_cell); - Assert (this->cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - - if (evaluate_val == true) - { - std::memcpy (&this->values_quad[0][0], &this->values_dofs[0][0], - dofs_per_cell * n_components * - sizeof (this->values_dofs[0][0])); -#ifdef DEBUG - this->values_quad_initialized = true; -#endif - } - // separate implementation here compared to the general case because the - // values are an identity operation - if (evaluate_grad == true) - { - for (unsigned int comp=0; comp (this->values_dofs[comp], - this->gradients_quad[comp][0]); - // grad y - apply_gradients<1,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][1]); - // grad y - apply_gradients<2,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][2]); - } - else if (dim == 2) - { - // grad x - apply_gradients<0,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][0]); - // grad y - apply_gradients<1,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][1]); - } - else if (dim == 1) - apply_gradients<0,true,false> (this->values_dofs[comp], - this->gradients_quad[comp][0]); - } -#ifdef DEBUG - this->gradients_quad_initialized = true; -#endif - } - if (evaluate_lapl == true) - { - for (unsigned int comp=0; comptemplate apply_hessians<0,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][0]); - // grad y - this->template apply_hessians<1,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][1]); - // grad y - this->template apply_hessians<2,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][2]); - - VectorizedArray temp1[n_q_points]; - // grad xy - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][3]); - // grad xz - apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][4]); - // grad yz - apply_gradients<1,true,false> (this->values_dofs[comp], temp1); - apply_gradients<2,true,false> (temp1, this->hessians_quad[comp][5]); - } - else if (dim == 2) - { - // grad x - this->template apply_hessians<0,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][0]); - // grad y - this->template apply_hessians<1,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][1]); - VectorizedArray temp1[n_q_points]; - // grad xy - apply_gradients<0,true,false> (this->values_dofs[comp], temp1); - apply_gradients<1,true,false> (temp1, this->hessians_quad[comp][2]); - } - else if (dim == 1) - this->template apply_hessians<0,true,false> (this->values_dofs[comp], - this->hessians_quad[comp][0]); - } -#ifdef DEBUG - this->hessians_quad_initialized = true; -#endif - } -} - - - -template -inline -void -FEEvaluationGL -::integrate (const bool integrate_val, const bool integrate_grad) -{ - this->check_dofs_per_cell(dofs_per_cell); - Assert (this->cell != numbers::invalid_unsigned_int, - ExcNotInitialized()); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_val == true) - std::memcpy (&this->values_dofs[0][0], &this->values_quad[0][0], - dofs_per_cell * n_components * - sizeof (this->values_dofs[0][0])); - if (integrate_grad == true) - { - for (unsigned int comp=0; comp (this->gradients_quad[comp][0], - this->values_dofs[comp]); - else - apply_gradients<0, false, false> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - - // grad y: can sum to temporary x value in temp2 - apply_gradients<1, false, true> (this->gradients_quad[comp][1], - this->values_dofs[comp]); - - // grad z: can sum to temporary x and y value in output - apply_gradients<2, false, true> (this->gradients_quad[comp][2], - this->values_dofs[comp]); - } - else if (dim == 2) - { - // grad x: If integrate_val == true we have to add to the previous output - if (integrate_val == true) - apply_gradients<0, false, true> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - else - apply_gradients<0, false, false> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - - // grad y: can sum to temporary x value in temp2 - apply_gradients<1, false, true> (this->gradients_quad[comp][1], - this->values_dofs[comp]); - } - else if (dim == 1) - { - if (integrate_val == true) - apply_gradients<0, false, true> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - else - apply_gradients<0, false, false> (this->gradients_quad[comp][0], - this->values_dofs[comp]); - - } - } - } + case 3: + if (integrate_grad == true) + for (unsigned int comp=0; comp (gradients_quad[comp][0], + values_dofs[comp]); + else + eval.template gradients<0, false, false> (gradients_quad[comp][0], + values_dofs[comp]); -#ifdef DEBUG - this->dof_values_initialized = true; -#endif -} + // grad y + eval.template gradients<1, false, true> (gradients_quad[comp][d1], + values_dofs[comp]); + // grad z + eval.template gradients<2, false, true> (gradients_quad[comp][d2], + values_dofs[comp]); + } + break; + + default: + AssertThrow(false, ExcNotImplemented()); + } + } +} // end of namespace internal -template -template -inline -void -FEEvaluationGL -::apply_gradients (const VectorizedArray in [], - VectorizedArray out []) -{ - internal::apply_tensor_product_gradients_gl, direction, dof_to_quad, add> - (this->data->shape_gradients.begin(), in, out); -} +/*-------------------------- FEEvaluation -----------------------------------*/ -/*------------------------- FEEvaluationDGP ---------------------------------*/ -template inline -FEEvaluationDGP -::FEEvaluationDGP (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) +FEEvaluation +::FEEvaluation (const MatrixFree &data_in, + const unsigned int fe_no, + const unsigned int quad_no) : - BaseClass (data_in, fe_no, quad_no) + BaseClass (data_in, fe_no, quad_no, fe_degree, n_q_points), + dofs_per_cell (this->data->dofs_per_cell) { - AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), - n_components * dofs_per_cell); + check_template_arguments(fe_no); + set_data_pointers(); } -template inline -FEEvaluationDGP -::FEEvaluationDGP (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) +FEEvaluation +::FEEvaluation (const MappingFEEvaluation &geometry, + const DoFHandler &dof_handler, + const unsigned int first_selected_component) : - BaseClass (geometry, dof_handler, first_selected_component) + BaseClass (geometry, dof_handler, first_selected_component), + dofs_per_cell (this->data->dofs_per_cell) { - AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), - n_components * dofs_per_cell); + check_template_arguments(numbers::invalid_unsigned_int); + set_data_pointers(); } -template inline -FEEvaluationDGP -::FEEvaluationDGP (const FEEvaluationDGP &other) +FEEvaluation +::FEEvaluation (const FEEvaluation &other) : - BaseClass (other) + BaseClass (other), + dofs_per_cell (this->data->dofs_per_cell) { + set_data_pointers(); } @@ -6844,148 +5906,249 @@ template inline void -FEEvaluationDGP -::evaluate (const bool evaluate_val, - const bool evaluate_grad, - const bool evaluate_lapl) +FEEvaluation +::set_data_pointers() { - this->check_dofs_per_cell(dofs_per_cell); - Assert (this->dof_values_initialized == true, - internal::ExcAccessToUninitializedField()); - - // expand dof_values to tensor product - VectorizedArray data_array[n_components*BaseClass::dofs_per_cell]; - VectorizedArray *expanded_dof_values[n_components]; - for (unsigned int c=0; cdata->dofs_per_cell, tensor_dofs_per_cell+2); + const unsigned int desired_dofs_per_cell = this->data->dofs_per_cell; - unsigned int count_p = 0, count_q = 0; - for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) + // set the pointers to the correct position in the data array + for (unsigned int c=0; c1?fe_degree+1-i:1); ++j) - { - for (unsigned int k=0; kvalues_dofs[c][count_p]; - for (unsigned int k=fe_degree+1-j-i; k(); - } - for (unsigned int j=fe_degree+1-i; j(); + this->values_dofs[c] = &my_data_array[c*desired_dofs_per_cell]; + this->values_quad[c] = &my_data_array[n_components*desired_dofs_per_cell+c*n_q_points]; + for (unsigned int d=0; dgradients_quad[c][d] = &my_data_array[n_components*(desired_dofs_per_cell+ + n_q_points) + + + (c*dim+d)*n_q_points]; + for (unsigned int d=0; d<(dim*dim+dim)/2; ++d) + this->hessians_quad[c][d] = &my_data_array[n_components*((dim+1)*n_q_points+ + desired_dofs_per_cell) + + + (c*(dim*dim+dim)+d)*n_q_points]; } - AssertDimension(count_q, BaseClass::dofs_per_cell); - AssertDimension(count_p, dofs_per_cell); - internal::do_evaluate (*this, expanded_dof_values, this->values_quad, - this->gradients_quad, this->hessians_quad, - evaluate_val, evaluate_grad, evaluate_lapl, - internal::int2type()); + switch (this->data->element_type) + { + case internal::MatrixFreeFunctions::tensor_symmetric: + evaluate_funct = + internal::FEEvaluationImpl::evaluate; + integrate_funct = + internal::FEEvaluationImpl::integrate; + break; -#ifdef DEBUG - if (evaluate_val == true) - this->values_quad_initialized = true; - if (evaluate_grad == true) - this->gradients_quad_initialized = true; - if (evaluate_lapl == true) - this->hessians_quad_initialized = true; -#endif -} + case internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0: + evaluate_funct = + internal::FEEvaluationImpl::evaluate; + integrate_funct = + internal::FEEvaluationImpl::integrate; + break; + case internal::MatrixFreeFunctions::tensor_general: + evaluate_funct = + internal::FEEvaluationImpl::evaluate; + integrate_funct = + internal::FEEvaluationImpl::integrate; + break; + case internal::MatrixFreeFunctions::tensor_gausslobatto: + evaluate_funct = + internal::FEEvaluationImpl::evaluate; + integrate_funct = + internal::FEEvaluationImpl::integrate; + break; -template -inline -void -FEEvaluationDGP -::integrate (bool integrate_val,bool integrate_grad) -{ - this->check_dofs_per_cell(dofs_per_cell); - if (integrate_val == true) - Assert (this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); - if (integrate_grad == true) - Assert (this->gradients_quad_submitted == true, - internal::ExcAccessToUninitializedField()); + case internal::MatrixFreeFunctions::truncated_tensor: + evaluate_funct = + internal::FEEvaluationImpl::evaluate; + integrate_funct = + internal::FEEvaluationImpl::integrate; + break; - VectorizedArray data_array[n_components*BaseClass::dofs_per_cell]; - VectorizedArray *expanded_dof_values[n_components]; - for (unsigned int c=0; cvalues_quad, - this->gradients_quad, integrate_val, integrate_grad, - internal::int2type()); - - // truncate tensor product - unsigned int count_p = 0, count_q = 0; - for (unsigned int i=0; i<(dim>2?fe_degree+1:1); ++i) - { - for (unsigned int j=0; j<(dim>1?fe_degree+1-i:1); ++j) - { - for (unsigned int k=0; kvalues_dofs[c][count_p] = expanded_dof_values[c][count_q]; - } - count_q += j+i; - } - count_q += i*(fe_degree+1); + default: + AssertThrow(false, ExcNotImplemented()); } - AssertDimension(count_q, BaseClass::dofs_per_cell); - AssertDimension(count_p, dofs_per_cell); -#ifdef DEBUG - this->dof_values_initialized = true; -#endif } -/*------------------------- FEEvaluationQ_DG0 -------------------------------*/ - -template inline -FEEvaluationQ_DG0 -::FEEvaluationQ_DG0 (const MatrixFree &data_in, - const unsigned int fe_no, - const unsigned int quad_no) - : - BaseClass (data_in, fe_no, quad_no) +void +FEEvaluation +::check_template_arguments(const unsigned int fe_no) { - AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), - n_components * dofs_per_cell); -} - +#ifdef DEBUG + // print error message when the dimensions do not match. Propose a possible + // fix + if (fe_degree != this->data->fe_degree + || + n_q_points != this->data->n_q_points) + { + std::string message = + "-------------------------------------------------------\n"; + message += "Illegal arguments in constructor/wrong template arguments!\n"; + message += " Called --> FEEvaluationquad_no) + ")\n"; + // check whether some other vector component has the correct number of + // points + unsigned int proposed_dof_comp = numbers::invalid_unsigned_int, + proposed_quad_comp = numbers::invalid_unsigned_int; + if (fe_no != numbers::invalid_unsigned_int) + { + if (fe_degree == this->data->fe_degree) + proposed_dof_comp = fe_no; + else + for (unsigned int no=0; nomatrix_info->n_components(); ++no) + if (this->matrix_info->get_shape_info(no,0,this->active_fe_index,0).fe_degree + == fe_degree) + { + proposed_dof_comp = no; + break; + } + if (n_q_points == + this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index]) + proposed_quad_comp = this->quad_no; + else + for (unsigned int no=0; nomapping_info->mapping_data_gen.size(); ++no) + if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index] + == n_q_points) + { + proposed_quad_comp = no; + break; + } + } + if (proposed_dof_comp != numbers::invalid_unsigned_int && + proposed_quad_comp != numbers::invalid_unsigned_int) + { + if (proposed_dof_comp != fe_no) + message += "Wrong vector component selection:\n"; + else + message += "Wrong quadrature formula selection:\n"; + message += " Did you mean FEEvaluationquad_no) + correct_pos += " ^\n"; + else + correct_pos += " \n"; + message += " " + correct_pos; + } + // ok, did not find the numbers specified by the template arguments in + // the given list. Suggest correct template arguments + const unsigned int proposed_n_q_points_1d = static_cast(std::pow(1.001*this->data->n_q_points,1./dim)); + message += "Wrong template arguments:\n"; + message += " Did you mean FEEvaluationdata->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"; + std::string correct_pos; + if (this->data->fe_degree != fe_degree) + correct_pos = " ^"; + else + correct_pos = " "; + if (proposed_n_q_points_1d != n_q_points_1d) + correct_pos += " ^\n"; + else + correct_pos += " \n"; + message += " " + correct_pos; -template -inline -FEEvaluationQ_DG0 -::FEEvaluationQ_DG0 (const MappingFEEvaluation &geometry, - const DoFHandler &dof_handler, - const unsigned int first_selected_component) - : - BaseClass (geometry, dof_handler, first_selected_component) -{ - AssertDimension(static_cast(this->values_quad[0]-this->values_dofs[0]), - n_components * dofs_per_cell); + Assert (fe_degree == this->data->fe_degree && + n_q_points == this->data->n_q_points, + ExcMessage(message)); + } + if (fe_no != numbers::invalid_unsigned_int) + { + AssertDimension (n_q_points, + this->mapping_info->mapping_data_gen[this->quad_no]. + n_q_points[this->active_quad_index]); + AssertDimension (this->data->dofs_per_cell * this->n_fe_components, + this->dof_info->dofs_per_cell[this->active_fe_index]); + } +#endif } -template inline -FEEvaluationQ_DG0 -::FEEvaluationQ_DG0 (const FEEvaluationQ_DG0 &other) - : - BaseClass (other) +Point > +FEEvaluation +::quadrature_point (const unsigned int q) const { + Assert (this->mapping_info->quadrature_points_initialized == true, + ExcNotInitialized()); + AssertIndexRange (q, n_q_points); + + // Cartesian mesh: not all quadrature points are stored, only the + // diagonal. Hence, need to find the tensor product index and retrieve the + // value from that + if (this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + Point > point (false); + switch (dim) + { + case 1: + return this->quadrature_points[q]; + case 2: + point[0] = this->quadrature_points[q%n_q_points_1d][0]; + point[1] = this->quadrature_points[q/n_q_points_1d][1]; + return point; + case 3: + point[0] = this->quadrature_points[q%n_q_points_1d][0]; + point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1]; + point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2]; + return point; + default: + Assert (false, ExcNotImplemented()); + return point; + } + } + // all other cases: just return the respective data as it is fully stored + else + return this->quadrature_points[q]; } @@ -6994,22 +6157,19 @@ template inline void -FEEvaluationQ_DG0 +FEEvaluation ::evaluate (const bool evaluate_val, const bool evaluate_grad, const bool evaluate_lapl) { - this->check_dofs_per_cell(dofs_per_cell); Assert (this->dof_values_initialized == true, internal::ExcAccessToUninitializedField()); - internal::do_evaluate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, this->hessians_quad, - evaluate_val, evaluate_grad, evaluate_lapl, - internal::int2type()); - if (evaluate_val) - for (unsigned int c=0; cvalues_quad[c][q] += this->values_dofs[c][dofs_per_cell-1]; + + // Select algorithm matching the element type at run time (the function + // pointer is easy to predict, so negligible in cost) + evaluate_funct (*this->data, &this->values_dofs[0], + this->values_quad, this->gradients_quad, this->hessians_quad, + evaluate_val, evaluate_grad, evaluate_lapl); #ifdef DEBUG if (evaluate_val == true) @@ -7027,10 +6187,9 @@ template inline void -FEEvaluationQ_DG0 +FEEvaluation ::integrate (bool integrate_val,bool integrate_grad) { - this->check_dofs_per_cell(dofs_per_cell); if (integrate_val == true) Assert (this->values_quad_submitted == true, internal::ExcAccessToUninitializedField()); @@ -7038,20 +6197,10 @@ FEEvaluationQ_DG0 Assert (this->gradients_quad_submitted == true, internal::ExcAccessToUninitializedField()); - if (integrate_val) - for (unsigned int c=0; cvalues_dofs[c][dofs_per_cell-1] = this->values_quad[c][0]; - for (unsigned int q=1; qvalues_dofs[c][dofs_per_cell-1] += this->values_quad[c][q]; - } - else - for (unsigned int c=0; cvalues_dofs[c][dofs_per_cell-1] = VectorizedArray(); - - internal::do_integrate (*this, this->values_dofs, this->values_quad, - this->gradients_quad, integrate_val, integrate_grad, - internal::int2type()); + // Select algorithm matching the element type at run time (the function + // pointer is easy to predict, so negligible in cost) + integrate_funct (*this->data, this->values_dofs, this->values_quad, + this->gradients_quad, integrate_val, integrate_grad); #ifdef DEBUG this->dof_values_initialized = true; diff --git a/deal.II/include/deal.II/matrix_free/operators.h b/deal.II/include/deal.II/matrix_free/operators.h index ac5670f928..d246ae4cdb 100644 --- a/deal.II/include/deal.II/matrix_free/operators.h +++ b/deal.II/include/deal.II/matrix_free/operators.h @@ -86,14 +86,9 @@ namespace MatrixFreeOperators const FEEvaluationBase &fe_eval; /** - * A structure to hold inverse shape functions, put into an aligned vector - * in order to make sure it can be allocated. + * A structure to hold inverse shape functions */ - struct Inverse1DShape - { - VectorizedArray evenodd[fe_degree+1][fe_degree/2+1]; - }; - AlignedVector inverse_shape; + AlignedVector > inverse_shape; }; @@ -112,18 +107,19 @@ namespace MatrixFreeOperators for (unsigned int j=0; j > + evaluator(inverse_shape, inverse_shape, inverse_shape); + const unsigned int shift_coefficient = inverse_coefficients.size() > dofs_per_cell ? dofs_per_cell : 0; const VectorizedArray *inv_coefficient = &inverse_coefficients[0]; @@ -183,34 +183,24 @@ namespace MatrixFreeOperators { const VectorizedArray* in = in_array+d*dofs_per_cell; VectorizedArray* out = out_array+d*dofs_per_cell; - internal::apply_tensor_product_evenodd,0,false,false,2> - (inverse_shape[0].evenodd, in, temp_data_field); - internal::apply_tensor_product_evenodd,1,false,false,2> - (inverse_shape[0].evenodd, temp_data_field, out); + // Need to select 'apply' method with hessian slot because values + // assume symmetries that do not exist in the inverse shapes + evaluator.template hessians<0,false,false> (in, temp_data_field); + evaluator.template hessians<1,false,false> (temp_data_field, out); if (dim == 3) { - internal::apply_tensor_product_evenodd,2,false,false,2> - (inverse_shape[0].evenodd, out, temp_data_field); + evaluator.template hessians<2,false,false> (out, temp_data_field); for (unsigned int q=0; q,2,true,false,2> - (inverse_shape[0].evenodd, temp_data_field, out); + evaluator.template hessians<2,true,false> (temp_data_field, out); } else if (dim == 2) for (unsigned int q=0; q,1,true,false,2> - (inverse_shape[0].evenodd, out, temp_data_field); - internal::apply_tensor_product_evenodd,0,true,false,2> - (inverse_shape[0].evenodd, temp_data_field, out); + evaluator.template hessians<1,true,false>(out, temp_data_field); + evaluator.template hessians<0,true,false>(temp_data_field, out); inv_coefficient += shift_coefficient; } diff --git a/deal.II/include/deal.II/matrix_free/shape_info.h b/deal.II/include/deal.II/matrix_free/shape_info.h index f9b51e21f6..01be04f08b 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.h @@ -34,6 +34,20 @@ namespace internal { namespace MatrixFreeFunctions { + /** + * An enum that encodes the type of element detected during + * initialization. FEEvaluation will select the most efficient algorithm + * based on the given element type. + */ + enum ElementType + { + tensor_general, + tensor_symmetric, + truncated_tensor, + tensor_symmetric_plus_dg0, + tensor_gausslobatto + }; + /** * The class that stores the shape functions, gradients and Hessians * evaluated for a tensor product finite element and tensor product @@ -76,6 +90,13 @@ namespace internal */ std::size_t memory_consumption () const; + /** + * Encodes the type of element detected at construction. FEEvaluation + * will select the most efficient algorithm based on the given element + * type. + */ + ElementType element_type; + /** * Stores the shape values of the 1D finite element evaluated on all 1D * quadrature points in vectorized format, i.e., as an array of @@ -103,6 +124,27 @@ namespace internal */ AlignedVector > shape_hessians; + /** + * Stores the shape values in a different format, namely the so-called + * even-odd scheme where the symmetries in shape_values are used for + * faster evaluation. + */ + AlignedVector > shape_val_evenodd; + + /** + * Stores the shape gradients in a different format, namely the + * so-called even-odd scheme where the symmetries in shape_gradients are + * used for faster evaluation. + */ + AlignedVector > shape_gra_evenodd; + + /** + * Stores the shape second derivatives in a different format, namely the + * so-called even-odd scheme where the symmetries in shape_hessians are + * used for faster evaluation. + */ + AlignedVector > shape_hes_evenodd; + /** * Stores the indices from cell DoFs to face DoFs. The rows go through * the 2*dim faces, and the columns the DoFs on the faces. @@ -173,6 +215,19 @@ namespace internal * Stores the number of DoFs per face in @p dim dimensions. */ unsigned int dofs_per_face; + + /** + * Checks whether we have symmetries in the shape values. In that case, + * also fill the shape_???_evenodd fields. + */ + bool check_1d_shapes_symmetric(const unsigned int n_q_points_1d); + + /** + * Checks whether symmetric 1D basis functions are such that the shape + * values form a diagonal matrix, which allows to use specialized + * algorithms that save some operations. + */ + bool check_1d_shapes_gausslobatto(); }; diff --git a/deal.II/include/deal.II/matrix_free/shape_info.templates.h b/deal.II/include/deal.II/matrix_free/shape_info.templates.h index a38ea4fe72..82e6566cf9 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.templates.h @@ -40,6 +40,7 @@ namespace internal template ShapeInfo::ShapeInfo () : + element_type (tensor_general), n_q_points (0), dofs_per_cell (0) {} @@ -85,6 +86,7 @@ namespace internal const FE_Q_DG0 *fe_q_dg0 = dynamic_cast*>(fe); + element_type = tensor_general; if (fe_poly != 0) scalar_lexicographic = fe_poly->get_poly_space_numbering_inverse(); else if (fe_poly_piece != 0) @@ -94,9 +96,13 @@ namespace internal scalar_lexicographic.resize(fe_dgp->dofs_per_cell); for (unsigned int i=0; idofs_per_cell; ++i) scalar_lexicographic[i] = i; + element_type = truncated_tensor; } else if (fe_q_dg0 != 0) - scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); + { + scalar_lexicographic = fe_q_dg0->get_poly_space_numbering_inverse(); + element_type = tensor_symmetric_plus_dg0; + } else Assert(false, ExcNotImplemented()); @@ -184,6 +190,17 @@ namespace internal this->face_gradient[1][i] = fe->shape_grad(my_i,q_point)[0]; } + if (element_type == tensor_general && + check_1d_shapes_symmetric(n_q_points_1d)) + { + if (check_1d_shapes_gausslobatto()) + element_type = tensor_gausslobatto; + else + element_type = tensor_symmetric; + } + else if (element_type == tensor_symmetric_plus_dg0) + check_1d_shapes_symmetric(n_q_points_1d); + // face information unsigned int n_faces = GeometryInfo::faces_per_cell; this->face_indices.reinit(n_faces, this->dofs_per_face); @@ -229,6 +246,137 @@ namespace internal + template + bool + ShapeInfo::check_1d_shapes_symmetric(const unsigned int n_q_points_1d) + { + const double zero_tol = + types_are_equal::value==true?1e-10:1e-7; + // symmetry for values + const unsigned int n_dofs_1d = fe_degree + 1; + for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) + for (unsigned int j=0; j zero_tol) + return false; + + // shape values should be zero at x=0.5 for all basis functions except + // for one which is one + if (n_q_points_1d%2 == 1 && n_dofs_1d%2 == 1) + { + for (unsigned int i=0; i zero_tol) + return false; + if (std::fabs(shape_values[(n_dofs_1d/2)*n_q_points_1d+ + n_q_points_1d/2][0]-1.)> zero_tol) + return false; + } + + // skew-symmetry for gradient, zero of middle basis function in middle + // quadrature point + for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) + for (unsigned int j=0; j zero_tol) + return false; + if (n_dofs_1d%2 == 1 && n_q_points_1d%2 == 1) + if (std::fabs(shape_gradients[(n_dofs_1d/2)*n_q_points_1d+ + (n_q_points_1d/2)][0]) > zero_tol) + return false; + + // symmetry for Laplacian + for (unsigned int i=0; i<(n_dofs_1d+1)/2; ++i) + for (unsigned int j=0; j zero_tol) + return false; + + const unsigned int stride = (n_q_points_1d+1)/2; + shape_val_evenodd.resize((fe_degree+1)*stride); + shape_gra_evenodd.resize((fe_degree+1)*stride); + shape_hes_evenodd.resize((fe_degree+1)*stride); + for (unsigned int i=0; i<(fe_degree+1)/2; ++i) + for (unsigned int q=0; q + bool + ShapeInfo::check_1d_shapes_gausslobatto() + { + if (dofs_per_cell != n_q_points) + return false; + + const double zero_tol = + types_are_equal::value==true?1e-10:1e-7; + // check: - identity operation for shape values + // - zero diagonal at interior points for gradients + // - gradient equal to unity at element boundary + const unsigned int n_points_1d = fe_degree+1; + for (unsigned int i=0; izero_tol) + return false; + } + else + { + if (std::fabs(shape_values[i*n_points_1d+ + j][0]-1.)>zero_tol) + return false; + } + for (unsigned int i=1; izero_tol) + return false; + if (std::fabs(shape_gradients[n_points_1d-1][0]- + (n_points_1d%2==0 ? -1. : 1.)) > zero_tol) + return false; + + return true; + } + + + template std::size_t ShapeInfo::memory_consumption () const @@ -236,7 +384,10 @@ namespace internal std::size_t memory = sizeof(*this); memory += MemoryConsumption::memory_consumption(shape_values); memory += MemoryConsumption::memory_consumption(shape_gradients); - memory += MemoryConsumption::memory_consumption(shape_hessians); + memory += MemoryConsumption::memory_consumption(shape_hessians); + memory += MemoryConsumption::memory_consumption(shape_val_evenodd); + memory += MemoryConsumption::memory_consumption(shape_gra_evenodd); + memory += MemoryConsumption::memory_consumption(shape_hes_evenodd); memory += face_indices.memory_consumption(); for (unsigned int i=0; i<2; ++i) { diff --git a/tests/matrix_free/estimate_condition_number_mass.cc b/tests/matrix_free/estimate_condition_number_mass.cc index dfd39b0910..bea4f8ac59 100644 --- a/tests/matrix_free/estimate_condition_number_mass.cc +++ b/tests/matrix_free/estimate_condition_number_mass.cc @@ -45,7 +45,7 @@ mass_operator (const MatrixFree &data, const Vector &src, const std::pair &cell_range) { - FEEvaluationGeneral fe_eval (data); + FEEvaluation fe_eval (data); const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell=cell_range.first; cell void test() { deallog << "Test " << M << " x " << N << std::endl; - double shape[M][N]; + AlignedVector shape(M*N); for (unsigned int i=0; i<(M+1)/2; ++i) for (unsigned int j=0; j evaluator(shape, shape, shape); if (type == 0) - internal::apply_tensor_product_values<1,M-1,N,double,0,false,add> - (&shape[0][0],x,y); + evaluator.template values<0,false,add> (x,y); if (type == 1) - internal::apply_tensor_product_gradients<1,M-1,N,double,0,false,add> - (&shape[0][0],x,y); + evaluator.template gradients<0,false,add> (x,y); if (type == 2) - internal::apply_tensor_product_hessians<1,M-1,N,double,0,false,add> - (&shape[0][0],x,y); + evaluator.template hessians<0,false,add> (x,y); deallog << "Errors no transpose: "; for (unsigned int i=0; i - (&shape[0][0],y,x); + evaluator.template values<0,true,add> (y,x); if (type == 1) - internal::apply_tensor_product_gradients<1,M-1,N,double,0,true,add> - (&shape[0][0],y,x); + evaluator.template gradients<0,true,add> (y,x); if (type == 2) - internal::apply_tensor_product_hessians<1,M-1,N,double,0,true,add> - (&shape[0][0],y,x); + evaluator.template hessians<0,true,add> (y,x); deallog << "Errors transpose: "; for (unsigned int i=0; i shape_sym(M*((N+1)/2)); for (unsigned int i=0; i(shape_sym,x,y); + internal::EvaluatorTensorProduct evaluator(shape_sym, shape_sym, shape_sym); + if (type == 0) + evaluator.template values<0,false,add> (x,y); + if (type == 1) + evaluator.template gradients<0,false,add> (x,y); + if (type == 2) + evaluator.template hessians<0,false,add> (x,y); + deallog << "Errors no transpose: "; for (unsigned int i=0; i(shape_sym,y,x); + if (type == 0) + evaluator.template values<0,true,add> (y,x); + if (type == 1) + evaluator.template gradients<0,true,add> (y,x); + if (type == 2) + evaluator.template hessians<0,true,add> (y,x); deallog << "Errors transpose: "; for (unsigned int i=0; i -class MatrixFreeTestGL : public MatrixFreeTest -{ -public: - - MatrixFreeTestGL(const MatrixFree &data, - const Mapping &mapping): - MatrixFreeTest(data, mapping) - {}; - - void operator() (const MatrixFree &data, - Vector &, - const Vector &src, - const std::pair &cell_range) const - { - FEEvaluationGL fe_eval (this->data); - for (unsigned int cell=cell_range.first; cell reference_values (fe_eval.n_q_points); - std::vector > reference_grads (fe_eval.n_q_points); - std::vector > reference_hess (fe_eval.n_q_points); - fe_eval.read_dof_values(src); - fe_eval.evaluate (true,true,true); - - // compare values with the ones the FEValues - // gives us. Those are seen as reference - for (unsigned int j=0; jfe_val.reinit (data.get_cell_iterator(cell,j)); - this->fe_val.get_function_values(src, reference_values); - this->fe_val.get_function_gradients(src, reference_grads); - this->fe_val.get_function_hessians(src, reference_hess); - - for (int q=0; q<(int)fe_eval.n_q_points; q++) - { - this->errors[0] += std::fabs(fe_eval.get_value(q)[j]- - reference_values[q]); - for (unsigned int d=0; derrors[1] += std::fabs(fe_eval.get_gradient(q)[d][j]- - reference_grads[q][d]); - this->errors[2] += std::fabs(fe_eval.get_laplacian(q)[j]- - trace(reference_hess[q])); - this->total[0] += std::fabs(reference_values[q]); - for (unsigned int d=0; dtotal[1] += std::fabs(reference_grads[q][d]); - this->total[2] += std::fabs(fe_eval.get_laplacian(q)[j]); - } - } - } - } -}; - - template void test () @@ -134,7 +80,7 @@ void test () data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; data.mapping_update_flags = update_gradients | update_second_derivatives; mf_data.reinit (mapping, dof, constraints, quad, data); - MatrixFreeTestGL mf (mf_data, mapping); + MatrixFreeTest mf (mf_data, mapping); mf.test_functions (solution); } diff --git a/tests/matrix_free/get_functions_q_hierarchical.cc b/tests/matrix_free/get_functions_q_hierarchical.cc index ce3c796e20..6e430c9dc7 100644 --- a/tests/matrix_free/get_functions_q_hierarchical.cc +++ b/tests/matrix_free/get_functions_q_hierarchical.cc @@ -32,61 +32,6 @@ std::ofstream logfile("output"); #include "get_functions_common.h" -template -class MatrixFreeTestGen : public MatrixFreeTest -{ -public: - MatrixFreeTestGen(const MatrixFree &data, - const Mapping &mapping): - MatrixFreeTest(data, mapping) - {}; - - void operator () (const MatrixFree &data, - Vector &, - const Vector &src, - const std::pair &cell_range) const - { - FEEvaluationGeneral fe_eval (data); - for (unsigned int cell=cell_range.first; cell reference_values (fe_eval.n_q_points); - std::vector > reference_grads (fe_eval.n_q_points); - std::vector > reference_hess (fe_eval.n_q_points); - fe_eval.read_dof_values(src); - fe_eval.evaluate (true,true,true); - - // compare values with the ones the FEValues - // gives us. Those are seen as reference - for (unsigned int j=0; jfe_val.reinit (data.get_cell_iterator(cell,j)); - this->fe_val.get_function_values(src, reference_values); - this->fe_val.get_function_gradients(src, reference_grads); - this->fe_val.get_function_hessians(src, reference_hess); - - for (int q=0; q<(int)fe_eval.n_q_points; q++) - { - this->errors[0] += std::fabs(fe_eval.get_value(q)[j]- - reference_values[q]); - for (unsigned int d=0; derrors[1] += std::fabs(fe_eval.get_gradient(q)[d][j]- - reference_grads[q][d]); - } - this->errors[2] += std::fabs(fe_eval.get_laplacian(q)[j]- - trace(reference_hess[q])); - this->total[0] += std::fabs(reference_values[q]); - for (unsigned int d=0; dtotal[1] += std::fabs(reference_grads[q][d]); - this->total[2] += std::fabs(fe_eval.get_laplacian(q)[j]); - } - } - } - } -}; - - template void test () @@ -137,7 +82,7 @@ void test () data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; data.mapping_update_flags = update_gradients | update_second_derivatives; mf_data.reinit (mapping, dof, constraints, quad, data); - MatrixFreeTestGen mf (mf_data, mapping); + MatrixFreeTest mf (mf_data, mapping); mf.test_functions (solution); } diff --git a/tests/matrix_free/matrix_vector_14.cc b/tests/matrix_free/matrix_vector_14.cc index 80126fc3ea..6d18bd6199 100644 --- a/tests/matrix_free/matrix_vector_14.cc +++ b/tests/matrix_free/matrix_vector_14.cc @@ -23,175 +23,12 @@ // approximation of the Laplacian. It only contains cell terms. #include "../tests.h" -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include std::ofstream logfile("output"); - - -template -void -helmholtz_operator_dgp (const MatrixFree &data, - VECTOR &dst, - const VECTOR &src, - const std::pair &cell_range) -{ - typedef typename VECTOR::value_type Number; - FEEvaluationDGP fe_eval (data); - const unsigned int n_q_points = fe_eval.n_q_points; - - for (unsigned int cell=cell_range.first; cell > -class MatrixFreeTest -{ -public: - typedef VectorizedArray vector_t; - - MatrixFreeTest(const MatrixFree &data_in): - data (data_in) - {}; - - void vmult (VECTOR &dst, - const VECTOR &src) const - { - dst = 0; - const std_cxx1x::function &, - VECTOR &, - const VECTOR &, - const std::pair &)> - wrap = helmholtz_operator_dgp; - data.cell_loop (wrap, dst, src); - }; - -private: - const MatrixFree &data; -}; - - - -template -void do_test (const DoFHandler &dof, - const ConstraintMatrix &constraints) -{ - - deallog << "Testing " << dof.get_fe().get_name() << std::endl; - - MatrixFree mf_data; - { - const QGauss<1> quad (fe_degree+1); - typename MatrixFree::AdditionalData data; - data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; - data.tasks_block_size = 3; - - mf_data.reinit (dof, constraints, quad, data); - } - - MatrixFreeTest mf (mf_data); - Vector in (dof.n_dofs()), out (dof.n_dofs()); - Vector in_dist (dof.n_dofs()); - Vector out_dist (in_dist); - - for (unsigned int i=0; i sparse_matrix (sparsity); - { - QGauss quadrature_formula(fe_degree+1); - - FEValues fe_values (dof.get_fe(), quadrature_formula, - update_values | update_gradients | - update_JxW_values); - - const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - cell_matrix = 0; - fe_values.reinit (cell); - - for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, - local_dof_indices, - sparse_matrix); - } - } - - sparse_matrix.vmult (out, in); - out -= out_dist; - const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); - - deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; -} +#include "matrix_vector_common.h" @@ -223,26 +60,3 @@ void test () do_test (dof, constraints); } - - -int main () -{ - deallog.attach(logfile); - deallog.depth_console(0); - - deallog << std::setprecision (3); - - { - deallog.threshold_double(5.e-11); - deallog.push("2d"); - test<2,1>(); - test<2,2>(); - deallog.pop(); - deallog.push("3d"); - test<3,1>(); - test<3,2>(); - deallog.pop(); - } -} - - diff --git a/tests/matrix_free/matrix_vector_14b.cc b/tests/matrix_free/matrix_vector_14b.cc index 04e33dd37f..65c44ec4b7 100644 --- a/tests/matrix_free/matrix_vector_14b.cc +++ b/tests/matrix_free/matrix_vector_14b.cc @@ -20,175 +20,11 @@ // degree would suggest #include "../tests.h" -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - std::ofstream logfile("output"); - - -template -void -helmholtz_operator_dgp (const MatrixFree &data, - VECTOR &dst, - const VECTOR &src, - const std::pair &cell_range) -{ - typedef typename VECTOR::value_type Number; - FEEvaluationDGP fe_eval (data); - const unsigned int n_q_points = fe_eval.n_q_points; - - for (unsigned int cell=cell_range.first; cell > -class MatrixFreeTest -{ -public: - typedef VectorizedArray vector_t; - - MatrixFreeTest(const MatrixFree &data_in): - data (data_in) - {}; - - void vmult (VECTOR &dst, - const VECTOR &src) const - { - dst = 0; - const std_cxx1x::function &, - VECTOR &, - const VECTOR &, - const std::pair &)> - wrap = helmholtz_operator_dgp; - data.cell_loop (wrap, dst, src); - }; - -private: - const MatrixFree &data; -}; - - - -template -void do_test (const DoFHandler &dof, - const ConstraintMatrix &constraints) -{ - - deallog << "Testing " << dof.get_fe().get_name() << std::endl; - - MatrixFree mf_data; - { - const QGauss<1> quad (fe_degree+2); - typename MatrixFree::AdditionalData data; - data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; - data.tasks_block_size = 3; - - mf_data.reinit (dof, constraints, quad, data); - } - - MatrixFreeTest mf (mf_data); - Vector in (dof.n_dofs()), out (dof.n_dofs()); - Vector in_dist (dof.n_dofs()); - Vector out_dist (in_dist); - - for (unsigned int i=0; i sparse_matrix (sparsity); - { - QGauss quadrature_formula(fe_degree+2); - - FEValues fe_values (dof.get_fe(), quadrature_formula, - update_values | update_gradients | - update_JxW_values); - - const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - cell_matrix = 0; - fe_values.reinit (cell); - - for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, - local_dof_indices, - sparse_matrix); - } - } - - sparse_matrix.vmult (out, in); - out -= out_dist; - const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); - - deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; -} +#include "matrix_vector_common.h" @@ -220,26 +56,3 @@ void test () do_test (dof, constraints); } - - -int main () -{ - deallog.attach(logfile); - deallog.depth_console(0); - - deallog << std::setprecision (3); - - { - deallog.threshold_double(5.e-11); - deallog.push("2d"); - test<2,1>(); - test<2,2>(); - deallog.pop(); - deallog.push("3d"); - test<3,1>(); - test<3,2>(); - deallog.pop(); - } -} - - diff --git a/tests/matrix_free/matrix_vector_17.cc b/tests/matrix_free/matrix_vector_17.cc index 4ab84fd056..3fd4d4bf4f 100644 --- a/tests/matrix_free/matrix_vector_17.cc +++ b/tests/matrix_free/matrix_vector_17.cc @@ -23,175 +23,11 @@ // in an approximation of the Laplacian. It only contains cell terms. #include "../tests.h" -#include #include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - std::ofstream logfile("output"); - - -template -void -helmholtz_operator_qdg0 (const MatrixFree &data, - VECTOR &dst, - const VECTOR &src, - const std::pair &cell_range) -{ - typedef typename VECTOR::value_type Number; - FEEvaluationQ_DG0 fe_eval (data); - const unsigned int n_q_points = fe_eval.n_q_points; - - for (unsigned int cell=cell_range.first; cell > -class MatrixFreeTest -{ -public: - typedef VectorizedArray vector_t; - - MatrixFreeTest(const MatrixFree &data_in): - data (data_in) - {}; - - void vmult (VECTOR &dst, - const VECTOR &src) const - { - dst = 0; - const std_cxx1x::function &, - VECTOR &, - const VECTOR &, - const std::pair &)> - wrap = helmholtz_operator_qdg0; - data.cell_loop (wrap, dst, src); - }; - -private: - const MatrixFree &data; -}; - - - -template -void do_test (const DoFHandler &dof, - const ConstraintMatrix &constraints) -{ - - deallog << "Testing " << dof.get_fe().get_name() << std::endl; - - MatrixFree mf_data; - { - const QGauss<1> quad (fe_degree+1); - typename MatrixFree::AdditionalData data; - data.tasks_parallel_scheme = - MatrixFree::AdditionalData::partition_color; - data.tasks_block_size = 3; - - mf_data.reinit (dof, constraints, quad, data); - } - - MatrixFreeTest mf (mf_data); - Vector in (dof.n_dofs()), out (dof.n_dofs()); - Vector in_dist (dof.n_dofs()); - Vector out_dist (in_dist); - - for (unsigned int i=0; i sparse_matrix (sparsity); - { - QGauss quadrature_formula(fe_degree+1); - - FEValues fe_values (dof.get_fe(), quadrature_formula, - update_values | update_gradients | - update_JxW_values); - - const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; - const unsigned int n_q_points = quadrature_formula.size(); - - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - std::vector local_dof_indices (dofs_per_cell); - - typename DoFHandler::active_cell_iterator - cell = dof.begin_active(), - endc = dof.end(); - for (; cell!=endc; ++cell) - { - cell_matrix = 0; - fe_values.reinit (cell); - - for (unsigned int q_point=0; q_pointget_dof_indices(local_dof_indices); - constraints.distribute_local_to_global (cell_matrix, - local_dof_indices, - sparse_matrix); - } - } - - sparse_matrix.vmult (out, in); - out -= out_dist; - const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); - - deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; -} +#include "matrix_vector_common.h" @@ -219,30 +55,13 @@ void test () DoFHandler dof (tria); dof.distribute_dofs(fe); ConstraintMatrix constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values (dof, 0, ZeroFunction(), + constraints); + constraints.close(); do_test (dof, constraints); } -int main () -{ - deallog.attach(logfile); - deallog.depth_console(0); - - deallog << std::setprecision (3); - - { - deallog.threshold_double(5.e-11); - deallog.push("2d"); - test<2,1>(); - test<2,2>(); - deallog.pop(); - deallog.push("3d"); - test<3,1>(); - test<3,2>(); - deallog.pop(); - } -} - - diff --git a/tests/matrix_free/matrix_vector_stokes_qdg0.cc b/tests/matrix_free/matrix_vector_stokes_qdg0.cc index c746da6c99..3ed9366440 100644 --- a/tests/matrix_free/matrix_vector_stokes_qdg0.cc +++ b/tests/matrix_free/matrix_vector_stokes_qdg0.cc @@ -72,8 +72,8 @@ public: const std::pair &cell_range) const { typedef VectorizedArray vector_t; - FEEvaluation velocity (data, 0); - FEEvaluationQ_DG0 pressure (data, 1); + FEEvaluation velocity (data, 0); + FEEvaluation pressure (data, 1); for (unsigned int cell=cell_range.first; cell fe_eval(data); - const unsigned int n_q_points = fe_eval.n_q_points; + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell=0; cell &cell_range) const { AssertDimension (src.size(), 2); - FEEvaluationGL current (data), old (data); + FEEvaluation current (data), old (data); for (unsigned int cell=cell_range.first; cell fe_eval(data); - const unsigned int n_q_points = fe_eval.n_q_points; + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell=0; cell &cell_range) const { AssertDimension (src.size(), 2); - FEEvaluationGL current (data), old (data); + FEEvaluation current (data), old (data); deallog << "submit / sine values: "; for (unsigned int cell=cell_range.first; cell fe_eval(data); - const unsigned int n_q_points = fe_eval.n_q_points; + FEEvaluation fe_eval(data); + const unsigned int n_q_points = fe_eval.n_q_points; for (unsigned int cell=0; cell &cell_range) const { AssertDimension (src.size(), 2); - FEEvaluationGL current (data), old (data); + FEEvaluation current (data), old (data); for (unsigned int cell=cell_range.first; cell