From e1067b66573b2fb673888d929ad49e617faaebc8 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 22 Mar 2024 17:22:18 +0100 Subject: [PATCH] FEEvaluation: Reduced specialized Access class --- include/deal.II/matrix_free/fe_evaluation.h | 2781 ++++++------------- 1 file changed, 884 insertions(+), 1897 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 2f59de1e17..afff5148bd 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -95,11 +95,24 @@ class FEEvaluationBase { public: using number_type = Number; - using value_type = Tensor<1, n_components_, VectorizedArrayType>; - using gradient_type = - Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>; - using hessian_type = - Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>; + using value_type = + std::conditional_t>; + using gradient_type = std::conditional_t< + n_components_ == 1, + Tensor<1, dim, VectorizedArrayType>, + std::conditional_t< + n_components_ == dim, + Tensor<2, dim, VectorizedArrayType>, + Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>>>; + using hessian_type = std::conditional_t< + n_components_ == 1, + Tensor<2, dim, VectorizedArrayType>, + std::conditional_t< + n_components_ == dim, + Tensor<3, dim, VectorizedArrayType>, + Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>>>>; static constexpr unsigned int dimension = dim; static constexpr unsigned int n_components = n_components_; static constexpr unsigned int n_lanes = VectorizedArrayType::size(); @@ -292,10 +305,6 @@ public: * the one set there. If @p integrate was called last, it instead * corresponds to the value of the integrated function with the test * function of the given index. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ value_type get_dof_value(const unsigned int dof) const; @@ -305,10 +314,6 @@ public: * component @p dof. Writes to the same field as is accessed through @p * get_dof_value. Therefore, the original data that was read from a vector * is overwritten as soon as a value is submitted. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ void submit_dof_value(const value_type val_in, const unsigned int dof); @@ -320,10 +325,6 @@ public: * a call to FEEvaluationBase::submit_value(). If the object is * vector-valued, a vector-valued return argument is given. Note that when * vectorization is enabled, values from several cells are grouped together. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ value_type get_value(const unsigned int q_point) const; @@ -335,10 +336,6 @@ public: * with EvaluationFlags::values set is called, this specifies the value * which is tested by all basis function on the current cell and integrated * over. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ void submit_value(const value_type val_in, const unsigned int q_point); @@ -348,10 +345,6 @@ public: * number @p q_point after a call to FEEvaluation::evaluate() with * EvaluationFlags::gradients, or the value that has been stored there with * a call to FEEvaluationBase::submit_gradient(). - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ gradient_type get_gradient(const unsigned int q_point) const; @@ -365,10 +358,6 @@ public: * * This call is equivalent to calling get_gradient() * normal_vector() * but will use a more efficient internal representation of data. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ value_type get_normal_derivative(const unsigned int q_point) const; @@ -380,10 +369,6 @@ public: * function FEEvaluation::integrate(EvaluationFlags::gradients) is called, * this specifies what is tested by all basis function gradients on the * current cell and integrated over. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ void submit_gradient(const gradient_type grad_in, const unsigned int q_point); @@ -401,52 +386,25 @@ public: * submit_gradient(). As a consequence, only one of these two can be * used. Usually, the contribution of a potential call to this function must * be added into the contribution for submit_gradient(). - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ void submit_normal_derivative(const value_type grad_in, const unsigned int q_point); - /** - * Write a contribution that is tested by the Hessian to the field - * containing the values at quadrature points with component @p q_point. - * Access to the same field as through get_hessian(). If applied before the - * function FEEvaluation::integrate(EvaluationFlags::hessians) is called, - * this specifies what is tested by the Hessians of all basis functions on the - * current cell and integrated over. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). - */ - void - submit_hessian(const hessian_type hessian_in, const unsigned int q_point); - /** * Return the Hessian of a finite element function at quadrature point * number @p q_point after a call to * FEEvaluation::evaluate(EvaluationFlags::hessians). If only the diagonal * or even the trace of the Hessian, the Laplacian, is needed, use the other * functions below. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ - Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>> + hessian_type get_hessian(const unsigned int q_point) const; /** * Return the diagonal of the Hessian of a finite element function at * quadrature point number @p q_point after a call to * FEEvaluation::evaluate(EvaluationFlags::hessians). - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ gradient_type get_hessian_diagonal(const unsigned int q_point) const; @@ -457,25 +415,29 @@ public: * FEEvaluation::evaluate(EvaluationFlags::hessians). Compared to the case * when computing the full Hessian, some operations can be saved when only * the Laplacian is requested. - * - * @note The derived class FEEvaluationAccess overloads this operation - * with specializations for the scalar case (n_components == 1) and for the - * vector-valued case (n_components == dim). */ value_type get_laplacian(const unsigned int q_point) const; -#ifdef DOXYGEN - // doxygen does not anyhow mention functions coming from partial template - // specialization of the base class, in this case FEEvaluationAccess. - // For now, hack in those functions manually only to fix documentation: + /** + * Write a contribution that is tested by the Hessian to the field + * containing the values at quadrature points with component @p q_point. + * Access to the same field as through get_hessian(). If applied before the + * function FEEvaluation::integrate(EvaluationFlags::hessians) is called, + * this specifies what is tested by the Hessians of all basis functions on the + * current cell and integrated over. + */ + void + submit_hessian(const hessian_type hessian_in, const unsigned int q_point); /** * Return the divergence of a vector-valued finite element at quadrature - * point number @p q_point after a call to @p evaluate(EvaluationFlags::gradients). + * point number @p q_point after a call to + * @p evaluate(EvaluationFlags::gradients). * * @note Only available for the vector-valued case (n_components == dim). */ + template > VectorizedArrayType get_divergence(const unsigned int q_point) const; @@ -487,6 +449,7 @@ public: * * @note Only available for the vector-valued case (n_components == dim). */ + template > SymmetricTensor<2, dim, VectorizedArrayType> get_symmetric_gradient(const unsigned int q_point) const; @@ -496,16 +459,19 @@ public: * * @note Only available for the vector-valued case (n_components == dim). */ + template > Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> get_curl(const unsigned int q_point) const; /** * Write a contribution that is tested by the divergence to the field * containing the values on quadrature points with component @p q_point. - * Access to the same field as through @p get_gradient. If applied before - * the function @p integrate(EvaluationFlags::gradients) is called, this specifies what is - * tested by all basis function gradients on the current cell and integrated - * over. + * Access to the same field as through get_gradient() and + * submit_gradient(). If applied before the function @p + * integrate(EvaluationFlags::gradients) is called, this quantity specifies + * what is tested by all basis function gradients on the current cell and + * integrated over. * * @note Only available for the vector-valued case (n_components == dim). * @@ -514,6 +480,7 @@ public: * used. Usually, the contribution of a potential call to this function must * be added into the diagonal of the contribution for submit_gradient(). */ + template > void submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point); @@ -521,10 +488,10 @@ public: /** * Write a contribution that is tested by the symmetric gradient to the field * containing the values on quadrature points with component @p q_point. - * Access to the same field as through @p get_symmetric_gradient. If applied before - * the function @p integrate(EvaluationFlags::gradients) is called, this specifies the - * symmetric gradient which is tested by all basis function symmetric - * gradients on the current cell and integrated over. + * Access to the same field as through @p get_symmetric_gradient. If applied + * before the function @p integrate(EvaluationFlags::gradients) is called, + * this specifies the symmetric gradient which is tested by all basis + * function symmetric gradients on the current cell and integrated over. * * @note Only available for the vector-valued case (n_components == dim). * @@ -532,8 +499,9 @@ public: * submit_gradient(). As a consequence, only one of these two can be * used. Usually, the contribution of a potential call to this function must * be added to the respective entries of the rank-2 tensor for - * submit_gradient(). + * submit_gradient(), in order not to overwrite information. */ + template > void submit_symmetric_gradient( const SymmetricTensor<2, dim, VectorizedArrayType> grad_in, @@ -551,12 +519,12 @@ public: * be added to the respective entries of the rank-2 tensor for * submit_gradient(). */ + template > void submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in, const unsigned int q_point); -#endif - /** * Take values collected at quadrature points via the submit_value() function, * multiply by the Jacobian determinant @@ -744,434 +712,24 @@ protected: -/** - * This class provides access to the data fields of the FEEvaluation classes. - * Generic access is achieved through the base class, and specializations for - * scalar and vector-valued elements are defined separately. - * - * @ingroup matrixfree - */ +// backward compatibility template , - typename Enable = void> -class FEEvaluationAccess : public FEEvaluationBase -{ - static_assert( - std::is_same_v, - "Type of Number and of VectorizedArrayType do not match."); - -public: - using number_type = Number; - using value_type = Tensor<1, n_components_, VectorizedArrayType>; - using gradient_type = - Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>; - static constexpr unsigned int dimension = dim; - static constexpr unsigned int n_components = n_components_; - using BaseClass = - FEEvaluationBase; - -protected: - /** - * Constructor. Made protected to prevent initialization in user code. 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 first_selected_component and @p - * quad_no allow to select the appropriate components. - */ - FEEvaluationAccess( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face = true, - const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int, - const unsigned int face_type = numbers::invalid_unsigned_int); - - /** - * Constructor with reduced functionality for similar usage of FEEvaluation - * as FEValues, including matrix assembly. - */ - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other); - - /** - * Copy constructor - */ - FEEvaluationAccess(const FEEvaluationAccess &other); - - /** - * Copy assignment operator - */ - FEEvaluationAccess & - operator=(const FEEvaluationAccess &other); -}; - - + typename VectorizedArrayType = VectorizedArray> +using FEEvaluationAccess DEAL_II_DEPRECATED = + FEEvaluationBase; /** - * This class provides access to the data fields of the FEEvaluation classes. - * Partial specialization for scalar fields that defines access with simple - * data fields, i.e., scalars for the values and Tensor<1,dim> for the - * gradients. + * 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 + * degree). For evaluation of face terms in DG, see the class + * FEFaceEvaluation. * - * @ingroup matrixfree - */ -template -class FEEvaluationAccess - : public FEEvaluationBase -{ - static_assert( - std::is_same_v, - "Type of Number and of VectorizedArrayType do not match."); - -public: - using number_type = Number; - using value_type = VectorizedArrayType; - using gradient_type = Tensor<1, dim, VectorizedArrayType>; - using hessian_type = Tensor<2, dim, VectorizedArrayType>; - static constexpr unsigned int dimension = dim; - using BaseClass = - FEEvaluationBase; - - /** - * @copydoc FEEvaluationBase::get_dof_value() - */ - value_type - get_dof_value(const unsigned int dof) const; - - /** - * @copydoc FEEvaluationBase::submit_dof_value() - */ - void - submit_dof_value(const value_type val_in, const unsigned int dof); - - /** - * @copydoc FEEvaluationBase::get_value() - */ - value_type - get_value(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::submit_value() - */ - void - submit_value(const value_type val_in, const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::submit_value() - */ - void - submit_value(const Tensor<1, 1, VectorizedArrayType> val_in, - const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::get_gradient() - */ - gradient_type - get_gradient(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::get_normal_derivative() - */ - value_type - get_normal_derivative(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::submit_gradient() - */ - void - submit_gradient(const gradient_type grad_in, const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::submit_normal_derivative() - */ - void - submit_normal_derivative(const value_type grad_in, - const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::get_hessian() - */ - hessian_type - get_hessian(unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::get_hessian_diagonal() - */ - gradient_type - get_hessian_diagonal(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::submit_hessian() - */ - void - submit_hessian(const hessian_type hessian_in, const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::get_laplacian() - */ - value_type - get_laplacian(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::integrate_value() - */ - value_type - integrate_value() const; - -protected: - /** - * Constructor. Made protected to avoid initialization in user code. 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 first_selected_component and @p - * quad_no allow to select the appropriate components. - */ - FEEvaluationAccess( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face = true, - const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int, - const unsigned int face_type = numbers::invalid_unsigned_int); - - /** - * Constructor with reduced functionality for similar usage of FEEvaluation - * as FEValues, including matrix assembly. - */ - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other); - - /** - * Copy constructor - */ - FEEvaluationAccess(const FEEvaluationAccess &other); - - /** - * Copy assignment operator - */ - FEEvaluationAccess & - operator=(const FEEvaluationAccess &other); -}; - - - -/** - * This class provides access to the data fields of the FEEvaluation classes. - * Partial specialization for fields with as many components as the underlying - * space dimension, i.e., values are of type Tensor<1,dim> and gradients of - * type Tensor<2,dim>. Provides some additional functions for access, like the - * symmetric gradient and divergence. - * - * @ingroup matrixfree - */ -template -class FEEvaluationAccess> - : public FEEvaluationBase -{ - static_assert( - std::is_same_v, - "Type of Number and of VectorizedArrayType do not match."); - -public: - using number_type = Number; - using value_type = Tensor<1, dim, VectorizedArrayType>; - using gradient_type = Tensor<2, dim, VectorizedArrayType>; - static constexpr unsigned int dimension = dim; - static constexpr unsigned int n_components = dim; - using BaseClass = - FEEvaluationBase; - - /** - * @copydoc FEEvaluationBase::get_value() - */ - value_type - get_value(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::get_gradient() - */ - gradient_type - get_gradient(const unsigned int q_point) const; - - /** - * Return the divergence of a vector-valued finite element at quadrature - * point number @p q_point after a call to @p evaluate(EvaluationFlags::gradients). - */ - VectorizedArrayType - get_divergence(const unsigned int q_point) const; - - /** - * Return the symmetric gradient of a vector-valued finite element at - * quadrature point number @p q_point after a call to @p - * evaluate(EvaluationFlags::gradients). It corresponds to 0.5 - * (grad+gradT). - */ - SymmetricTensor<2, dim, VectorizedArrayType> - get_symmetric_gradient(const unsigned int q_point) const; - - /** - * Return the curl of the vector field, $\nabla \times v$ after a call to @p - * evaluate(EvaluationFlags::gradients). - */ - Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> - get_curl(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::get_hessian() - */ - Tensor<3, dim, VectorizedArrayType> - get_hessian(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::get_hessian_diagonal() - */ - gradient_type - get_hessian_diagonal(const unsigned int q_point) const; - - /** - * @copydoc FEEvaluationBase::submit_value() - */ - void - submit_value(const Tensor<1, dim, VectorizedArrayType> val_in, - const unsigned int q_point); - - /** - * @copydoc FEEvaluationBase::submit_gradient() - */ - void - submit_gradient(const gradient_type grad_in, const unsigned int q_point); - - /** - * Write a contribution that is tested by the gradient to the field - * containing the values on quadrature points with component @p q_point. - * This function is an alternative to the other submit_gradient function - * when using a system of fixed number of equations which happens to - * coincide with the dimension for some dimensions, but not all. To allow - * for dimension-independent programming, this function can be used instead. - */ - void - submit_gradient( - const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in, - const unsigned int q_point); - - /** - * Write a contribution that is tested by the divergence to the field - * containing the values on quadrature points with component @p q_point. - * Access to the same field as through @p get_gradient. If applied before - * the function @p integrate(EvaluationFlags::gradients) is called, this specifies what is - * tested by all basis function gradients on the current cell and integrated - * over. - */ - void - submit_divergence(const VectorizedArrayType div_in, - const unsigned int q_point); - - /** - * Write a contribution that is tested by the symmetric gradient to the field - * containing the values on quadrature points with component @p q_point. - * Access to the same field as through @p get_symmetric_gradient. If applied before - * the function @p integrate(EvaluationFlags::gradients) is called, this specifies the - * symmetric gradient which is tested by all basis function symmetric - * gradients on the current cell and integrated over. - */ - void - submit_symmetric_gradient( - const SymmetricTensor<2, dim, VectorizedArrayType> grad_in, - const unsigned int q_point); - - /** - * Write the components of a curl containing the values on quadrature point - * @p q_point. Access to the same data field as through @p get_gradient. - */ - void - submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl_in, - const unsigned int q_point); - -protected: - /** - * Constructor. Made protected to avoid initialization in user code. 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 first_selected_component and @p - * quad_no allow to select the appropriate components. - */ - FEEvaluationAccess( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int dofs_per_cell, - const unsigned int n_q_points, - const bool is_interior_face = true, - const unsigned int active_fe_index = numbers::invalid_unsigned_int, - const unsigned int active_quad_index = numbers::invalid_unsigned_int, - const unsigned int face_type = numbers::invalid_unsigned_int); - - /** - * Constructor with reduced functionality for similar usage of FEEvaluation - * as FEValues, including matrix assembly. - */ - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other); - - /** - * Copy constructor - */ - FEEvaluationAccess(const FEEvaluationAccess &other); - - /** - * Copy assignment operator - */ - FEEvaluationAccess & - operator=(const FEEvaluationAccess &other); -}; - - - -/** - * 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 - * degree). For evaluation of face terms in DG, see the class - * FEFaceEvaluation. - * - *

Usage and initialization

+ *

Usage and initialization

* *

Fast usage in combination with MatrixFree

* @@ -1466,9 +1024,9 @@ protected: * * 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. + * writing to vectors. Furthermore, functionality to access to values, + * gradients and Hessians of the finite element function at quadrature points + * is inherited. * * This class assumes that the shape functions of the FiniteElement under * consideration do not depend on the geometry of the cells in real @@ -1582,14 +1140,13 @@ protected: * @endcode * * In a similar vein, the submit_value() and submit_gradient() calls take - * tensors of values. Note that there exist specializations for @p - * n_components=1 and @p n_components=dim, which are provided through the base - * class FEEvaluationAccess. In the scalar case, these provide the scalar - * return types described above. In the vector-valued case, the gradient is - * converted from Tensor@<1,dim,Tensor@<1,dim,VectorizedArray@ - * @> @> to Tensor@<2,dim,VectorizedArray@ - * @>. Furthermore, additional operations such as the diveregence or - * curl are available. + * tensors of values. Note that there exist specializations of these types for + * @p n_components=1 and @p n_components=dim. In the scalar case, these + * provide the scalar return types described above. In the vector-valued case, + * the gradient is converted from + * Tensor@<1,dim,Tensor@<1,dim,VectorizedArray@ @> @> to + * Tensor@<2,dim,VectorizedArray@ @>. Furthermore, + * additional operations such as the divergence or curl are available. * * In case different shape functions are combined, for example mixed finite * element formulations in Stokes flow, two FEEvaluation objects are created, @@ -1679,7 +1236,8 @@ protected: * * This observation also translates to the case when different differential * operators are implemented in a program, for example the action of a mass - * matrix for one phase of the algorithm and the action of a @ref GlossStiffnessMatrix "stiffness matrix" + * matrix for one phase of the algorithm and the action of a + * @ref GlossStiffnessMatrix "stiffness matrix" * in another one. Only a single MatrixFree object is necessary, maintaining * full efficiency by using different local functions with the respective * implementation in separate FEEvaluation objects. In other words, a user @@ -1732,11 +1290,11 @@ template -class FEEvaluation : public FEEvaluationAccess +class FEEvaluation : public FEEvaluationBase { static_assert( std::is_same_v, @@ -1747,7 +1305,7 @@ public: * An alias to the base class. */ using BaseClass = - FEEvaluationAccess; + FEEvaluationBase; /** * An underlying number type specified as template argument. @@ -2160,7 +1718,7 @@ private: * quadrature points and face integrations. The design of the class is similar * to FEEvaluation and most of the interfaces are shared with that class, in * particular most access functions that come from the common base classes - * FEEvaluationAccess and FEEvaluationBase. Furthermore, the relation of this + * FEEvaluationData and FEEvaluationBase. Furthermore, the relation of this * class to FEEvaluation is similar to the relation between FEValues and * FEFaceValues. * @@ -2196,11 +1754,11 @@ template > -class FEFaceEvaluation : public FEEvaluationAccess +class FEFaceEvaluation : public FEEvaluationBase { static_assert( std::is_same_v, @@ -2211,7 +1769,7 @@ public: * An alias to the base class. */ using BaseClass = - FEEvaluationAccess; + FEEvaluationBase; /** * A underlying number type specified as template argument. @@ -4390,16 +3948,26 @@ template -inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType> -FEEvaluationBase:: - get_dof_value(const unsigned int dof) const +inline DEAL_II_ALWAYS_INLINE + typename FEEvaluationBase::value_type + FEEvaluationBase:: + get_dof_value(const unsigned int dof) const { AssertIndexRange(dof, this->data->dofs_per_component_on_cell); - const std::size_t dofs = this->data->dofs_per_component_on_cell; - Tensor<1, n_components_, VectorizedArrayType> return_value; - for (unsigned int comp = 0; comp < n_components; ++comp) - return_value[comp] = this->values_dofs[comp * dofs + dof]; - return return_value; + if constexpr (n_components == 1) + return this->values_dofs[dof]; + else + { + const std::size_t dofs = this->data->dofs_per_component_on_cell; + Tensor<1, n_components_, VectorizedArrayType> return_value; + for (unsigned int comp = 0; comp < n_components; ++comp) + return_value[comp] = this->values_dofs[comp * dofs + dof]; + return return_value; + } } @@ -4409,9 +3977,14 @@ template -inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType> -FEEvaluationBase:: - get_value(const unsigned int q_point) const +inline DEAL_II_ALWAYS_INLINE + typename FEEvaluationBase::value_type + FEEvaluationBase:: + get_value(const unsigned int q_point) const { # ifdef DEBUG Assert(this->values_quad_initialized == true, @@ -4419,11 +3992,81 @@ FEEvaluationBase:: # endif AssertIndexRange(q_point, this->n_quadrature_points); - const std::size_t nqp = this->n_quadrature_points; - Tensor<1, n_components_, VectorizedArrayType> return_value; - for (unsigned int comp = 0; comp < n_components; ++comp) - return_value[comp] = this->values_quad[comp * nqp + q_point]; - return return_value; + if constexpr (n_components == 1) + return this->values_quad[q_point]; + else + { + if (n_components == dim && + this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { + // Piola transform is required +# ifdef DEBUG + Assert(this->values_quad_initialized == true, + internal::ExcAccessToUninitializedField()); +# endif + + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_values")); + const std::size_t nqp = this->n_quadrature_points; + Tensor<1, n_components, VectorizedArrayType> value_out; + + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + // Cartesian cell + const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; + const VectorizedArrayType inv_det = + (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] : + this->jacobian[0][0][0] * this->jacobian[0][1][1] * + this->jacobian[0][2][2]; + + // J * u * det(J^-1) + for (unsigned int comp = 0; comp < n_components; ++comp) + value_out[comp] = this->values_quad[comp * nqp + q_point] * + jac[comp][comp] * inv_det; + } + else + { + // Affine or general cell + const Tensor<2, dim, VectorizedArrayType> inv_t_jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->jacobian[q_point] : + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + transpose(invert(inv_t_jac)) : + this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into + // account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac); + // J * u * det(J^-1) + for (unsigned int comp = 0; comp < n_components; ++comp) + { + value_out[comp] = this->values_quad[q_point] * jac[comp][0]; + for (unsigned int e = 1; e < dim; ++e) + value_out[comp] += + this->values_quad[e * nqp + q_point] * jac[comp][e]; + value_out[comp] *= inv_det; + } + } + return value_out; + } + else + { + const std::size_t nqp = this->n_quadrature_points; + Tensor<1, n_components_, VectorizedArrayType> return_value; + for (unsigned int comp = 0; comp < n_components; ++comp) + return_value[comp] = this->values_quad[comp * nqp + q_point]; + return return_value; + } + } } @@ -4434,7 +4077,11 @@ template inline DEAL_II_ALWAYS_INLINE - Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> + typename FEEvaluationBase::gradient_type FEEvaluationBase:: get_gradient(const unsigned int q_point) const { @@ -4448,6 +4095,192 @@ inline DEAL_II_ALWAYS_INLINE internal::ExcMatrixFreeAccessToUninitializedMappingField( "update_gradients")); const std::size_t nqp = this->n_quadrature_points; + + if constexpr (n_components == dim && dim > 1) + { + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { + // Piola transform is required +# ifdef DEBUG + Assert(this->gradients_quad_initialized == true, + internal::ExcAccessToUninitializedField()); +# endif + + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + const std::size_t nqp = this->n_quadrature_points; + const std::size_t nqp_d = nqp * dim; + Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out; + const VectorizedArrayType *gradients = + this->gradients_quad + q_point * dim; + + + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + // Cartesian cell + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + this->jacobian[1]; + const VectorizedArrayType inv_det = + (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] : + this->jacobian[0][0][0] * this->jacobian[0][1][1] * + this->jacobian[0][2][2]; + + // J * grad_quad * J^-1 * det(J^-1) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int comp = 0; comp < n_components; ++comp) + grad_out[comp][d] = gradients[comp * nqp_d + d] * + inv_t_jac[d][d] * + (jac[comp][comp] * inv_det); + } + else if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + // Affine cell + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into + // account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac); + + VectorizedArrayType tmp[dim][dim]; + // J * grad_quad * J^-1 * det(J^-1) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0]; + for (unsigned int f = 1; f < dim; ++f) + tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f]; + } + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType res = jac[comp][0] * tmp[d][0]; + for (unsigned int f = 1; f < dim; ++f) + res += jac[comp][f] * tmp[d][f]; + + grad_out[comp][d] = res * inv_det; + } + } + else + { + // General cell + + // This assert could be removed if we make sure that this is + // updated even though update_hessians or update_jacobian_grads is + // not passed, i.e make the necessary changes in + // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags + Assert(this->jacobian_gradients_non_inverse != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_hessians")); + + const auto jac_grad = + this->jacobian_gradients_non_inverse[q_point]; + const Tensor<2, dim, VectorizedArrayType> inv_t_jac = + this->jacobian[q_point]; + + // Derivatives are reordered for faces. Need to take this into + // account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac); + const Tensor<2, dim, VectorizedArrayType> t_jac = + invert(inv_t_jac); + + // (J * grad_quad) * J^-1 * det(J^-1), part in braces + VectorizedArrayType tmp[dim][dim]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e]; + for (unsigned int f = 1; f < dim; ++f) + tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e]; + } + + // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms + // outside braces with gradient part from above + for (unsigned int d = 0; d < dim; ++d) + { + for (unsigned int e = 0; e < dim; ++e) + tmp[e][d] += + jac_grad[e][d] * this->values_quad[e * nqp + q_point]; + for (unsigned int f = 0, r = dim; f < dim; ++f) + for (unsigned int k = f + 1; k < dim; ++k, ++r) + { + tmp[k][d] += + jac_grad[r][d] * this->values_quad[f * nqp + q_point]; + tmp[f][d] += + jac_grad[r][d] * this->values_quad[k * nqp + q_point]; + } + } + + // Apply J^{-1} appearing in both terms outside braces above + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0]; + for (unsigned int f = 1; f < dim; ++f) + res += tmp[f][d] * inv_t_jac[e][f]; + grad_out[d][e] = res; + } + + // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})), + // which can be expressed as a rank-1 update tmp[d] * tmp4[e], + // where tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1}) + VectorizedArrayType tmp3[dim], tmp4[dim]; + for (unsigned int d = 0; d < dim; ++d) + { + tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0]; + for (unsigned int e = 1; e < dim; ++e) + tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e]; + } + for (unsigned int e = 0, k = dim; e < dim; ++e) + for (unsigned int f = e + 1; f < dim; ++k, ++f) + for (unsigned int d = 0; d < dim; ++d) + { + tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d]; + tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d]; + } + for (unsigned int d = 0; d < dim; ++d) + { + tmp4[d] = tmp3[0] * inv_t_jac[d][0]; + for (unsigned int e = 1; e < dim; ++e) + tmp4[d] += tmp3[e] * inv_t_jac[d][e]; + } + + VectorizedArrayType tmp2[dim]; + for (unsigned int d = 0; d < dim; ++d) + { + tmp2[d] = t_jac[0][d] * this->values_quad[q_point]; + for (unsigned e = 1; e < dim; ++e) + tmp2[d] += + t_jac[e][d] * this->values_quad[e * nqp + q_point]; + } + + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + grad_out[d][e] -= tmp4[e] * tmp2[d]; + + // finally multiply by det(J^{-1}) necessary in all + // contributions above + grad_out[d][e] *= inv_det; + } + } + return grad_out; + } + } Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out; // Cartesian cell @@ -4477,7 +4310,10 @@ inline DEAL_II_ALWAYS_INLINE this->gradients_quad[(comp * nqp + q_point) * dim + e]; } } - return grad_out; + if constexpr (n_components == 1) + return grad_out[0]; + else + return grad_out; } @@ -4487,9 +4323,14 @@ template -inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType> -FEEvaluationBase:: - get_normal_derivative(const unsigned int q_point) const +inline DEAL_II_ALWAYS_INLINE + typename FEEvaluationBase::value_type + FEEvaluationBase:: + get_normal_derivative(const unsigned int q_point) const { AssertIndexRange(q_point, this->n_quadrature_points); # ifdef DEBUG @@ -4523,7 +4364,10 @@ FEEvaluationBase:: this->normal_x_jacobian[index][d]; } } - return grad_out; + if constexpr (n_components == 1) + return grad_out[0]; + else + return grad_out; } @@ -4586,7 +4430,11 @@ template -inline Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>> +inline typename FEEvaluationBase::hessian_type FEEvaluationBase:: get_hessian(const unsigned int q_point) const { @@ -4713,7 +4561,10 @@ FEEvaluationBase:: hessian_out[comp][e][d] = hessian_out[comp][d][e]; } } - return hessian_out; + if constexpr (n_components == 1) + return hessian_out[0]; + else + return hessian_out; } @@ -4723,7 +4574,11 @@ template -inline Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> +inline typename FEEvaluationBase::gradient_type FEEvaluationBase:: get_hessian_diagonal(const unsigned int q_point) const { @@ -4802,7 +4657,11 @@ FEEvaluationBase:: this->gradients_quad[(comp * nqp + q_point) * dim + e]; } } - return hessian_out; + + if constexpr (n_components == 1) + return hessian_out[0]; + else + return hessian_out; } @@ -4812,7 +4671,11 @@ template -inline Tensor<1, n_components_, VectorizedArrayType> +inline typename FEEvaluationBase::value_type FEEvaluationBase:: get_laplacian(const unsigned int q_point) const { @@ -4823,15 +4686,25 @@ FEEvaluationBase:: # endif AssertIndexRange(q_point, this->n_quadrature_points); - Tensor<1, n_components_, VectorizedArrayType> laplacian_out; - const auto hess_diag = get_hessian_diagonal(q_point); - for (unsigned int comp = 0; comp < n_components; ++comp) + const gradient_type hess_diag = get_hessian_diagonal(q_point); + if constexpr (n_components == 1) { - laplacian_out[comp] = hess_diag[comp][0]; + VectorizedArrayType sum = hess_diag[0]; for (unsigned int d = 1; d < dim; ++d) - laplacian_out[comp] += hess_diag[comp][d]; + sum += hess_diag[d]; + return sum; + } + else + { + Tensor<1, n_components_, VectorizedArrayType> laplacian_out; + for (unsigned int comp = 0; comp < n_components; ++comp) + { + laplacian_out[comp] = hess_diag[comp][0]; + for (unsigned int d = 1; d < dim; ++d) + laplacian_out[comp] += hess_diag[comp][d]; + } + return laplacian_out; } - return laplacian_out; } @@ -4843,8 +4716,7 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationBase:: - submit_dof_value(const Tensor<1, n_components_, VectorizedArrayType> val_in, - const unsigned int dof) + submit_dof_value(const value_type val_in, const unsigned int dof) { # ifdef DEBUG this->dof_values_initialized = true; @@ -4852,7 +4724,10 @@ FEEvaluationBase:: const std::size_t dofs = this->data->dofs_per_component_on_cell; AssertIndexRange(dof, this->data->dofs_per_component_on_cell); for (unsigned int comp = 0; comp < n_components; ++comp) - this->values_dofs[comp * dofs + dof] = val_in[comp]; + if constexpr (n_components == 1) + this->values_dofs[comp * dofs + dof] = val_in; + else + this->values_dofs[comp * dofs + dof] = val_in[comp]; } @@ -4864,8 +4739,7 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationBase:: - submit_value(const Tensor<1, n_components_, VectorizedArrayType> val_in, - const unsigned int q_point) + submit_value(const value_type val_in, const unsigned int q_point) { # ifdef DEBUG Assert(this->is_reinitialized, ExcNotInitialized()); @@ -4881,18 +4755,78 @@ FEEvaluationBase:: const std::size_t nqp = this->n_quadrature_points; VectorizedArrayType *values = this->values_quad + q_point; - if (this->cell_type <= internal::MatrixFreeFunctions::affine) - { - const VectorizedArrayType JxW = - this->J_value[0] * this->quadrature_weights[q_point]; - for (unsigned int comp = 0; comp < n_components; ++comp) - values[comp * nqp] = val_in[comp] * JxW; - } + const VectorizedArrayType JxW = + this->cell_type <= internal::MatrixFreeFunctions::affine ? + this->J_value[0] * this->quadrature_weights[q_point] : + this->J_value[q_point]; + if constexpr (n_components == 1) + values[0] = val_in * JxW; else { - const VectorizedArrayType JxW = this->J_value[q_point]; - for (unsigned int comp = 0; comp < n_components; ++comp) - values[comp * nqp] = val_in[comp] * JxW; + if (n_components == dim && + this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { + // Piola transform is required + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_value")); +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); + this->values_quad_submitted = true; +# endif + + VectorizedArrayType *values = this->values_quad + q_point; + const std::size_t nqp = this->n_quadrature_points; + + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; + const VectorizedArrayType weight = + this->quadrature_weights[q_point]; + + for (unsigned int comp = 0; comp < n_components; ++comp) + values[comp * nqp] = val_in[comp] * weight * jac[comp][comp]; + } + else + { + // Affine or general cell + const Tensor<2, dim, VectorizedArrayType> inv_t_jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->jacobian[q_point] : + this->jacobian[0]; + + // Derivatives are reordered for faces. Need to take this into + // account and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? + this->quadrature_weights[q_point] : + (((this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->J_value[q_point] : + this->J_value[0] * this->quadrature_weights[q_point]) * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac))); + const Tensor<2, dim, VectorizedArrayType> jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + transpose(invert(inv_t_jac)) : + this->jacobian[1]; + + // J^T * u * factor + for (unsigned int comp = 0; comp < n_components; ++comp) + { + values[comp * nqp] = val_in[0] * jac[0][comp]; + for (unsigned int e = 1; e < dim; ++e) + values[comp * nqp] += val_in[e] * jac[e][comp]; + values[comp * nqp] *= fac; + } + } + } + else + for (unsigned int comp = 0; comp < n_components; ++comp) + values[comp * nqp] = val_in[comp] * JxW; } } @@ -4905,9 +4839,7 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationBase:: - submit_gradient( - const Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_in, - const unsigned int q_point) + submit_gradient(const gradient_type grad_in, const unsigned int q_point) { # ifdef DEBUG Assert(this->is_reinitialized, ExcNotInitialized()); @@ -4923,41 +4855,242 @@ FEEvaluationBase:: this->gradients_quad_submitted = true; # endif - const std::size_t nqp_d = this->n_quadrature_points * dim; - VectorizedArrayType *gradients = this->gradients_quad + q_point * dim; - - if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) + if constexpr (dim > 1 && n_components == dim) { - const VectorizedArrayType JxW = - this->J_value[0] * this->quadrature_weights[q_point]; - std::array jac; - for (unsigned int d = 0; d < dim; ++d) - jac[d] = this->jacobian[0][d][d]; - for (unsigned int d = 0; d < dim; ++d) + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) { - const VectorizedArrayType factor = jac[d] * JxW; - for (unsigned int comp = 0; comp < n_components; ++comp) - gradients[comp * nqp_d + d] = grad_in[comp][d] * factor; - } - } - else - { - const Tensor<2, dim, VectorizedArrayType> jac = - this->cell_type > internal::MatrixFreeFunctions::affine ? - this->jacobian[q_point] : - this->jacobian[0]; - const VectorizedArrayType JxW = - this->cell_type > internal::MatrixFreeFunctions::affine ? - this->J_value[q_point] : - this->J_value[0] * this->quadrature_weights[q_point]; - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0]; - for (unsigned int e = 1; e < dim; ++e) - new_val += (jac[e][d] * grad_in[comp][e]); - gradients[comp * nqp_d + d] = new_val * JxW; - } + // Piola transform is required + +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); +# endif + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); +# ifdef DEBUG + this->gradients_quad_submitted = true; +# endif + + VectorizedArrayType *gradients = this->gradients_quad + q_point * dim; + VectorizedArrayType *values = + this->values_from_gradients_quad + q_point; + const std::size_t nqp = this->n_quadrature_points; + const std::size_t nqp_d = nqp * dim; + + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + // Cartesian cell + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + this->jacobian[1]; + const VectorizedArrayType weight = + this->quadrature_weights[q_point]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int comp = 0; comp < n_components; ++comp) + gradients[comp * nqp_d + d] = grad_in[comp][d] * + inv_t_jac[d][d] * + (jac[comp][comp] * weight); + } + else if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + // Affine cell + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into + // account and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? + this->quadrature_weights[q_point] : + this->J_value[0] * this->quadrature_weights[q_point] * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac)); + + // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor + VectorizedArrayType tmp[dim][dim]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0]; + for (unsigned int f = 1; f < dim; ++f) + tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f]; + } + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType res = jac[0][comp] * tmp[d][0]; + for (unsigned int f = 1; f < dim; ++f) + res += jac[f][comp] * tmp[d][f]; + + gradients[comp * nqp_d + d] = res * fac; + } + } + else + { + // General cell + + const auto jac_grad = + this->jacobian_gradients_non_inverse[q_point]; + const Tensor<2, dim, VectorizedArrayType> inv_t_jac = + this->jacobian[q_point]; + + // Derivatives are reordered for faces. Need to take this into + // account and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? this->quadrature_weights[q_point] : + this->J_value[q_point] * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac)); + const Tensor<2, dim, VectorizedArrayType> t_jac = + invert(inv_t_jac); + + // Start evaluation for values part below to enable the compiler + // to possibly re-use the same computation in get_gradient() + // without interfering with stores to 'gradients' + VectorizedArrayType tmp3[dim], tmp4[dim]; + for (unsigned int d = 0; d < dim; ++d) + { + tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0]; + for (unsigned int e = 1; e < dim; ++e) + tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e]; + } + for (unsigned int e = 0, k = dim; e < dim; ++e) + for (unsigned int f = e + 1; f < dim; ++k, ++f) + for (unsigned int d = 0; d < dim; ++d) + { + tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d]; + tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d]; + } + for (unsigned int d = 0; d < dim; ++d) + { + tmp4[d] = tmp3[0] * inv_t_jac[d][0]; + for (unsigned int e = 1; e < dim; ++e) + tmp4[d] += tmp3[e] * inv_t_jac[d][e]; + } + + const Tensor<2, dim, VectorizedArrayType> grad_in_scaled = + fac * grad_in; + + VectorizedArrayType tmp[dim][dim]; + + // J * (J^{-1} * (grad_in * factor)) + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0]; + for (unsigned int f = 1; f < dim; ++f) + tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f]; + } + + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + { + VectorizedArrayType res = t_jac[d][0] * tmp[e][0]; + for (unsigned int f = 1; f < dim; ++f) + res += t_jac[d][f] * tmp[e][f]; + + gradients[d * nqp_d + e] = res; + } + + // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces + // as 'tmp' from above + VectorizedArrayType value[dim]; + for (unsigned int d = 0; d < dim; ++d) + { + value[d] = tmp[d][0] * jac_grad[d][0]; + for (unsigned int e = 1; e < dim; ++e) + value[d] += tmp[d][e] * jac_grad[d][e]; + } + for (unsigned int e = 0, k = dim; e < dim; ++e) + for (unsigned int f = e + 1; f < dim; ++k, ++f) + for (unsigned int d = 0; d < dim; ++d) + { + value[e] += tmp[f][d] * jac_grad[k][d]; + value[f] += tmp[e][d] * jac_grad[k][d]; + } + + // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1}) + // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ ) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0]; + for (unsigned int e = 1; e < dim; ++e) + tmp2 += grad_in_scaled[d][e] * tmp4[e]; + for (unsigned int e = 0; e < dim; ++e) + value[e] -= t_jac[e][d] * tmp2; + } + + for (unsigned int d = 0; d < dim; ++d) + values[d * nqp] = value[d]; + } + return; + } + } + + const std::size_t nqp_d = this->n_quadrature_points * dim; + VectorizedArrayType *gradients = this->gradients_quad + q_point * dim; + + if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const VectorizedArrayType JxW = + this->J_value[0] * this->quadrature_weights[q_point]; + + // Load all entries before starting to write back to make sure the + // compiler sees opportunity of loads in a possibly nearby + // get_gradient() function (i.e., the compiler should not think that + // 'jacobian' could alias with 'gradients'). + std::array jac; + for (unsigned int d = 0; d < dim; ++d) + jac[d] = this->jacobian[0][d][d]; + + for (unsigned int d = 0; d < dim; ++d) + { + const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW; + if constexpr (n_components == 1) + gradients[d] = grad_in[d] * factor; + else + for (unsigned int comp = 0; comp < n_components; ++comp) + gradients[comp * nqp_d + d] = grad_in[comp][d] * factor; + } + } + else + { + const Tensor<2, dim, VectorizedArrayType> jac = + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->jacobian[q_point] : + this->jacobian[0]; + const VectorizedArrayType JxW = + this->cell_type > internal::MatrixFreeFunctions::affine ? + this->J_value[q_point] : + this->J_value[0] * this->quadrature_weights[q_point]; + if constexpr (n_components == 1) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType new_val = jac[0][d] * grad_in[0]; + for (unsigned int e = 1; e < dim; ++e) + new_val += (jac[e][d] * grad_in[e]); + gradients[d] = new_val * JxW; + } + else + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType new_val = jac[0][d] * grad_in[comp][0]; + for (unsigned int e = 1; e < dim; ++e) + new_val += (jac[e][d] * grad_in[comp][e]); + gradients[comp * nqp_d + d] = new_val * JxW; + } } } @@ -4970,9 +5103,7 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationBase:: - submit_normal_derivative( - const Tensor<1, n_components_, VectorizedArrayType> grad_in, - const unsigned int q_point) + submit_normal_derivative(const value_type grad_in, const unsigned int q_point) { AssertIndexRange(q_point, this->n_quadrature_points); Assert(this->normal_x_jacobian != nullptr, @@ -4994,7 +5125,10 @@ FEEvaluationBase:: { for (unsigned int d = 0; d < dim - 1; ++d) gradients[comp * nqp_d + d] = VectorizedArrayType(); - gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac; + if constexpr (n_components == 1) + gradients[dim - 1] = grad_in * JxW_jac; + else + gradients[comp * nqp_d + dim - 1] = grad_in[comp] * JxW_jac; } } else @@ -5008,10 +5142,11 @@ FEEvaluationBase:: this->J_value[index] * this->quadrature_weights[q_point] : this->J_value[index]; for (unsigned int comp = 0; comp < n_components; ++comp) - { - for (unsigned int d = 0; d < dim; ++d) + for (unsigned int d = 0; d < dim; ++d) + if constexpr (n_components == 1) + gradients[d] = (grad_in * JxW) * jac[d]; + else gradients[comp * nqp_d + d] = (grad_in[comp] * JxW) * jac[d]; - } } } @@ -5024,10 +5159,7 @@ template inline DEAL_II_ALWAYS_INLINE void FEEvaluationBase:: - submit_hessian( - const Tensor<1, n_components_, Tensor<2, dim, VectorizedArrayType>> - hessian_in, - const unsigned int q_point) + submit_hessian(const hessian_type hessian_in, const unsigned int q_point) { # ifdef DEBUG Assert(this->is_reinitialized, ExcNotInitialized()); @@ -5057,8 +5189,12 @@ FEEvaluationBase:: const auto jac_d = this->jacobian[0][d][d]; const VectorizedArrayType factor = jac_d * jac_d * JxW; for (unsigned int comp = 0; comp < n_components; ++comp) - this->hessians_quad[(comp * hdim + d) * nqp + q_point] = - hessian_in[comp][d][d] * factor; + if constexpr (n_components == 1) + this->hessians_quad[d * nqp + q_point] = + hessian_in[d][d] * factor; + else + this->hessians_quad[(comp * hdim + d) * nqp + q_point] = + hessian_in[comp][d][d] * factor; } // off diagonal part @@ -5069,8 +5205,12 @@ FEEvaluationBase:: const auto jac_e = this->jacobian[0][e][e]; const VectorizedArrayType factor = jac_d * jac_e * JxW; for (unsigned int comp = 0; comp < n_components; ++comp) - this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] = - (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor; + if constexpr (n_components == 1) + this->hessians_quad[off_dia * nqp + q_point] = + (hessian_in[d][e] + hessian_in[e][d]) * factor; + else + this->hessians_quad[(comp * hdim + off_dia) * nqp + q_point] = + (hessian_in[comp][d][e] + hessian_in[comp][e][d]) * factor; } } // cell with general Jacobian, but constant within the cell @@ -5081,14 +5221,20 @@ FEEvaluationBase:: this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int comp = 0; comp < n_components; ++comp) { - // 1. tmp = hessian_in(u) * J + Tensor<2, dim, VectorizedArrayType> hessian_c; + if constexpr (n_components == 1) + hessian_c = hessian_in; + else + hessian_c = hessian_in[comp]; + + // 1. tmp = hessian(u) * J VectorizedArrayType tmp[dim][dim]; for (unsigned int i = 0; i < dim; ++i) for (unsigned int j = 0; j < dim; ++j) { - tmp[i][j] = hessian_in[comp][i][0] * jac[0][j]; + tmp[i][j] = hessian_c[i][0] * jac[0][j]; for (unsigned int k = 1; k < dim; ++k) - tmp[i][j] += hessian_in[comp][i][k] * jac[k][j]; + tmp[i][j] += hessian_c[i][k] * jac[k][j]; } // 2. hessian_unit = J^T * tmp @@ -5120,14 +5266,20 @@ FEEvaluationBase:: const auto &jac_grad = this->jacobian_gradients[q_point]; for (unsigned int comp = 0; comp < n_components; ++comp) { - // 1. tmp = hessian_in(u) * J + Tensor<2, dim, VectorizedArrayType> hessian_c; + if constexpr (n_components == 1) + hessian_c = hessian_in; + else + hessian_c = hessian_in[comp]; + + // 1. tmp = hessian(u) * J VectorizedArrayType tmp[dim][dim]; for (unsigned int i = 0; i < dim; ++i) for (unsigned int j = 0; j < dim; ++j) { - tmp[i][j] = hessian_in[comp][i][0] * jac[0][j]; + tmp[i][j] = hessian_c[i][0] * jac[0][j]; for (unsigned int k = 1; k < dim; ++k) - tmp[i][j] += hessian_in[comp][i][k] * jac[k][j]; + tmp[i][j] += hessian_c[i][k] * jac[k][j]; } // 2. hessian_unit = J^T * tmp @@ -5136,907 +5288,81 @@ FEEvaluationBase:: for (unsigned int j = 0; j < dim; ++j) { tmp2[i][j] = jac[0][i] * tmp[0][j]; - for (unsigned int k = 1; k < dim; ++k) - tmp2[i][j] += jac[k][i] * tmp[k][j]; - } - - // diagonal part - for (unsigned int d = 0; d < dim; ++d) - this->hessians_quad[(comp * hdim + d) * nqp + q_point] = - tmp2[d][d] * JxW; - - // off diagonal part - for (unsigned int d = 0, off_diag = dim; d < dim; ++d) - for (unsigned int e = d + 1; e < dim; ++e, ++off_diag) - this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] = - (tmp2[d][e] + tmp2[e][d]) * JxW; - - // 3. gradient_unit = J' ** hessian_in - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType sum = 0; - for (unsigned int e = 0; e < dim; ++e) - sum += hessian_in[comp][e][e] * jac_grad[e][d]; - for (unsigned int e = 0, count = dim; e < dim; ++e) - for (unsigned int f = e + 1; f < dim; ++f, ++count) - sum += (hessian_in[comp][e][f] + hessian_in[comp][f][e]) * - jac_grad[count][d]; - this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim + - d] = sum * JxW; - } - } - } -} - - - -template -inline Tensor<1, n_components_, VectorizedArrayType> -FEEvaluationBase:: - integrate_value() const -{ -# ifdef DEBUG - Assert(this->is_reinitialized, ExcNotInitialized()); - Assert(this->values_quad_submitted == true, - internal::ExcAccessToUninitializedField()); -# endif - - Tensor<1, n_components_, VectorizedArrayType> return_value; - const std::size_t nqp = this->n_quadrature_points; - for (unsigned int q = 0; q < nqp; ++q) - for (unsigned int comp = 0; comp < n_components; ++comp) - return_value[comp] += this->values_quad[comp * nqp + q]; - return (return_value); -} - - - -/*----------------------- FEEvaluationAccess --------------------------------*/ - - -template -inline FEEvaluationAccess< - dim, - n_components_, - Number, - is_face, - VectorizedArrayType, - Enable>::FEEvaluationAccess(const MatrixFree - &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index, - const unsigned int active_quad_index, - const unsigned int face_type) - : FEEvaluationBase( - matrix_free, - dof_no, - first_selected_component, - quad_no, - fe_degree, - n_q_points, - is_interior_face, - active_fe_index, - active_quad_index, - face_type) -{} - - - -template -inline FEEvaluationAccess:: - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other) - : FEEvaluationBase( - mapping, - fe, - quadrature, - update_flags, - first_selected_component, - other) -{} - - - -template -inline FEEvaluationAccess< - dim, - n_components_, - Number, - is_face, - VectorizedArrayType, - Enable>::FEEvaluationAccess(const FEEvaluationAccess &other) - : FEEvaluationBase( - other) -{} - - - -template -inline FEEvaluationAccess & -FEEvaluationAccess< - dim, - n_components_, - Number, - is_face, - VectorizedArrayType, - Enable>::operator=(const FEEvaluationAccess &other) -{ - this->FEEvaluationBase::operator=(other); - return *this; -} - - - -/*-------------------- FEEvaluationAccess scalar ----------------------------*/ - - -template -inline FEEvaluationAccess:: - FEEvaluationAccess( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index, - const unsigned int active_quad_index, - const unsigned int face_type) - : FEEvaluationBase( - matrix_free, - dof_no, - first_selected_component, - quad_no, - fe_degree, - n_q_points, - is_interior_face, - active_fe_index, - active_quad_index, - face_type) -{} - - - -template -inline FEEvaluationAccess:: - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other) - : FEEvaluationBase( - mapping, - fe, - quadrature, - update_flags, - first_selected_component, - other) -{} - - - -template -inline FEEvaluationAccess:: - FEEvaluationAccess( - const FEEvaluationAccess - &other) - : FEEvaluationBase(other) -{} - - - -template -inline FEEvaluationAccess & -FEEvaluationAccess::operator=( - const FEEvaluationAccess &other) -{ - this - ->FEEvaluationBase::operator=( - other); - return *this; -} - - - -template -inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationAccess::get_dof_value( - const unsigned int dof) const -{ - AssertIndexRange(dof, this->data->dofs_per_component_on_cell); - return this->values_dofs[dof]; -} - - - -template -inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationAccess::get_value( - const unsigned int q_point) const -{ -# ifdef DEBUG - Assert(this->values_quad_initialized == true, - internal::ExcAccessToUninitializedField()); -# endif - AssertIndexRange(q_point, this->n_quadrature_points); - return this->values_quad[q_point]; -} - - - -template -inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationAccess:: - get_normal_derivative(const unsigned int q_point) const -{ - return BaseClass::get_normal_derivative(q_point)[0]; -} - - - -template -inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> -FEEvaluationAccess::get_gradient( - const unsigned int q_point_in) const -{ - // could use the base class gradient, but that involves too many expensive - // initialization operations on tensors - -# ifdef DEBUG - Assert(this->gradients_quad_initialized == true, - internal::ExcAccessToUninitializedField()); -# endif - AssertIndexRange(q_point_in, this->n_quadrature_points); - - Assert(this->jacobian != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); - - Tensor<1, dim, VectorizedArrayType> grad_out; - - const std::size_t q_point = q_point_in; - if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - for (unsigned int d = 0; d < dim; ++d) - grad_out[d] = - this->gradients_quad[dim * q_point + d] * this->jacobian[0][d][d]; - } - // cell with general/affine Jacobian - else - { - const Tensor<2, dim, VectorizedArrayType> &jac = - this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ? - q_point : - 0]; - for (unsigned int d = 0; d < dim; ++d) - { - grad_out[d] = jac[d][0] * this->gradients_quad[dim * q_point]; - for (unsigned int e = 1; e < dim; ++e) - grad_out[d] += jac[d][e] * this->gradients_quad[dim * q_point + e]; - } - } - return grad_out; -} - - - -template -inline Tensor<2, dim, VectorizedArrayType> -FEEvaluationAccess::get_hessian( - const unsigned int q_point) const -{ - return BaseClass::get_hessian(q_point)[0]; -} - - - -template -inline Tensor<1, dim, VectorizedArrayType> -FEEvaluationAccess:: - get_hessian_diagonal(const unsigned int q_point) const -{ - return BaseClass::get_hessian_diagonal(q_point)[0]; -} - - - -template -inline VectorizedArrayType -FEEvaluationAccess::get_laplacian( - const unsigned int q_point) const -{ - return BaseClass::get_laplacian(q_point)[0]; -} - - - -template -inline void DEAL_II_ALWAYS_INLINE -FEEvaluationAccess:: - submit_dof_value(const VectorizedArrayType val_in, const unsigned int dof) -{ -# ifdef DEBUG - this->dof_values_initialized = true; - AssertIndexRange(dof, this->data->dofs_per_component_on_cell); -# endif - this->values_dofs[dof] = val_in; -} - - - -template -inline void DEAL_II_ALWAYS_INLINE -FEEvaluationAccess::submit_value( - const VectorizedArrayType val_in, - const unsigned int q_point) -{ -# ifdef DEBUG - Assert(this->is_reinitialized, ExcNotInitialized()); -# endif - AssertIndexRange(q_point, this->n_quadrature_points); - Assert(this->J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_value")); -# ifdef DEBUG - this->values_quad_submitted = true; -# endif - - if (this->cell_type <= internal::MatrixFreeFunctions::affine) - { - const VectorizedArrayType JxW = - this->J_value[0] * this->quadrature_weights[q_point]; - this->values_quad[q_point] = val_in * JxW; - } - else // if (this->cell_type < internal::MatrixFreeFunctions::general) - { - this->values_quad[q_point] = val_in * this->J_value[q_point]; - } -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess::submit_value( - const Tensor<1, 1, VectorizedArrayType> val_in, - const unsigned int q_point) -{ - submit_value(val_in[0], q_point); -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess:: - submit_normal_derivative(const VectorizedArrayType grad_in, - const unsigned int q_point) -{ - Tensor<1, 1, VectorizedArrayType> grad; - grad[0] = grad_in; - BaseClass::submit_normal_derivative(grad, q_point); -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess:: - submit_gradient(const Tensor<1, dim, VectorizedArrayType> grad_in, - const unsigned int q_point_in) -{ -# ifdef DEBUG - Assert(this->is_reinitialized, ExcNotInitialized()); -# endif - AssertIndexRange(q_point_in, this->n_quadrature_points); - Assert(this->J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); - Assert(this->jacobian != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); -# ifdef DEBUG - this->gradients_quad_submitted = true; -# endif - - const std::size_t q_point = q_point_in; - VectorizedArrayType *grad_ptr = this->gradients_quad + dim * q_point; - if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - const VectorizedArrayType JxW = - this->J_value[0] * this->quadrature_weights[q_point]; - - // Make sure the compiler does not think 'jacobian' is aliased with - // 'gradients_quad' - std::array jac; - for (unsigned int d = 0; d < dim; ++d) - jac[d] = this->jacobian[0][d][d]; - - for (unsigned int d = 0; d < dim; ++d) - grad_ptr[d] = grad_in[d] * jac[d] * JxW; - } - // general/affine cell type - else - { - const Tensor<2, dim, VectorizedArrayType> jac = - this->cell_type > internal::MatrixFreeFunctions::affine ? - this->jacobian[q_point] : - this->jacobian[0]; - const VectorizedArrayType JxW = - this->cell_type > internal::MatrixFreeFunctions::affine ? - this->J_value[q_point] : - this->J_value[0] * this->quadrature_weights[q_point]; - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType new_val = jac[0][d] * grad_in[0]; - for (unsigned int e = 1; e < dim; ++e) - new_val += jac[e][d] * grad_in[e]; - grad_ptr[d] = new_val * JxW; - } - } -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess:: - submit_hessian(const Tensor<2, dim, VectorizedArrayType> hessian_in, - const unsigned int q_point) -{ - Tensor<1, 1, Tensor<2, dim, VectorizedArrayType>> hessian; - hessian[0] = hessian_in; - BaseClass::submit_hessian(hessian, q_point); -} - - - -template -inline VectorizedArrayType -FEEvaluationAccess:: - integrate_value() const -{ - return BaseClass::integrate_value()[0]; -} - - - -/*----------------- FEEvaluationAccess vector-valued ------------------------*/ - - -template -inline FEEvaluationAccess>:: - FEEvaluationAccess( - const MatrixFree &matrix_free, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face, - const unsigned int active_fe_index, - const unsigned int active_quad_index, - const unsigned int face_type) - : FEEvaluationBase( - matrix_free, - dof_no, - first_selected_component, - quad_no, - fe_degree, - n_q_points, - is_interior_face, - active_fe_index, - active_quad_index, - face_type) -{} - - - -template -inline FEEvaluationAccess>:: - FEEvaluationAccess( - const Mapping &mapping, - const FiniteElement &fe, - const Quadrature<1> &quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationData *other) - : FEEvaluationBase( - mapping, - fe, - quadrature, - update_flags, - first_selected_component, - other) -{} - - - -template -inline FEEvaluationAccess>:: - FEEvaluationAccess(const FEEvaluationAccess &other) - : FEEvaluationBase(other) -{} - - - -template -inline FEEvaluationAccess> & -FEEvaluationAccess< - dim, - dim, - Number, - is_face, - VectorizedArrayType, - std::enable_if_t>::operator=(const FEEvaluationAccess &other) -{ - this->FEEvaluationBase:: - operator=(other); - return *this; -} - - -template -inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> - FEEvaluationAccess>::get_value(const unsigned int - q_point) const -{ - if (this->data->element_type == - internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) - { - // Piola transform is required -# ifdef DEBUG - Assert(this->values_quad_initialized == true, - internal::ExcAccessToUninitializedField()); -# endif - - AssertIndexRange(q_point, this->n_quadrature_points); - Assert(this->J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_values")); - const std::size_t nqp = this->n_quadrature_points; - Tensor<1, dim, VectorizedArrayType> value_out; - - if (!is_face && - this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - // Cartesian cell - const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; - const VectorizedArrayType inv_det = - (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] : - this->jacobian[0][0][0] * this->jacobian[0][1][1] * - this->jacobian[0][2][2]; - - // J * u * det(J^-1) - for (unsigned int comp = 0; comp < n_components; ++comp) - value_out[comp] = this->values_quad[comp * nqp + q_point] * - jac[comp][comp] * inv_det; - } - else - { - // Affine or general cell - const Tensor<2, dim, VectorizedArrayType> inv_t_jac = - (this->cell_type > internal::MatrixFreeFunctions::affine) ? - this->jacobian[q_point] : - this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> jac = - (this->cell_type > internal::MatrixFreeFunctions::affine) ? - transpose(invert(inv_t_jac)) : - this->jacobian[1]; - - // Derivatives are reordered for faces. Need to take this into account - const VectorizedArrayType inv_det = - (is_face && dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac); - // J * u * det(J^-1) - for (unsigned int comp = 0; comp < n_components; ++comp) - { - value_out[comp] = this->values_quad[q_point] * jac[comp][0]; - for (unsigned int e = 1; e < dim; ++e) - value_out[comp] += - this->values_quad[e * nqp + q_point] * jac[comp][e]; - value_out[comp] *= inv_det; - } - } - return value_out; - } - else - { - // No Piola needed - return BaseClass::get_value(q_point); - } -} - - - -template -inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> - FEEvaluationAccess>::get_gradient(const unsigned int - q_point) const -{ - if (this->data->element_type == - internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) - { - // Piola transform is required -# ifdef DEBUG - Assert(this->gradients_quad_initialized == true, - internal::ExcAccessToUninitializedField()); -# endif - - AssertIndexRange(q_point, this->n_quadrature_points); - Assert(this->jacobian != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); - const std::size_t nqp = this->n_quadrature_points; - const std::size_t nqp_d = nqp * dim; - Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out; - const VectorizedArrayType *gradients = - this->gradients_quad + q_point * dim; - - - if (!is_face && - this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - // Cartesian cell - const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = - this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; - const VectorizedArrayType inv_det = - (dim == 2) ? this->jacobian[0][0][0] * this->jacobian[0][1][1] : - this->jacobian[0][0][0] * this->jacobian[0][1][1] * - this->jacobian[0][2][2]; - - // J * grad_quad * J^-1 * det(J^-1) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int comp = 0; comp < n_components; ++comp) - grad_out[comp][d] = gradients[comp * nqp_d + d] * - inv_t_jac[d][d] * (jac[comp][comp] * inv_det); - } - else if (this->cell_type <= internal::MatrixFreeFunctions::affine) - { - // Affine cell - const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = - this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; - - // Derivatives are reordered for faces. Need to take this into account - const VectorizedArrayType inv_det = - (is_face && dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac); - - VectorizedArrayType tmp[dim][dim]; - // J * grad_quad * J^-1 * det(J^-1) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - tmp[d][e] = inv_t_jac[d][0] * gradients[e * nqp_d + 0]; - for (unsigned int f = 1; f < dim; ++f) - tmp[d][e] += inv_t_jac[d][f] * gradients[e * nqp_d + f]; - } - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType res = jac[comp][0] * tmp[d][0]; - for (unsigned int f = 1; f < dim; ++f) - res += jac[comp][f] * tmp[d][f]; - - grad_out[comp][d] = res * inv_det; - } - } - else - { - // General cell - - // This assert could be removed if we make sure that this is updated - // even though update_hessians or update_jacobian_grads is not passed, - // i.e make the necessary changes in - // MatrixFreeFunctions::MappingInfoStorage::compute_update_flags - Assert(this->jacobian_gradients_non_inverse != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_hessians")); - - const auto jac_grad = this->jacobian_gradients_non_inverse[q_point]; - const Tensor<2, dim, VectorizedArrayType> inv_t_jac = - this->jacobian[q_point]; - - // Derivatives are reordered for faces. Need to take this into account - const VectorizedArrayType inv_det = - (is_face && dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac); - const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac); - - // (J * grad_quad) * J^-1 * det(J^-1), part in braces - VectorizedArrayType tmp[dim][dim]; - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - tmp[e][d] = t_jac[0][d] * gradients[0 * nqp_d + e]; - for (unsigned int f = 1; f < dim; ++f) - tmp[e][d] += t_jac[f][d] * gradients[f * nqp_d + e]; - } - - // Add (jac_grad * values) * J^{-1} * det(J^{-1}), combine terms - // outside braces with gradient part from above - for (unsigned int d = 0; d < dim; ++d) - { - for (unsigned int e = 0; e < dim; ++e) - tmp[e][d] += - jac_grad[e][d] * this->values_quad[e * nqp + q_point]; - for (unsigned int f = 0, r = dim; f < dim; ++f) - for (unsigned int k = f + 1; k < dim; ++k, ++r) - { - tmp[k][d] += - jac_grad[r][d] * this->values_quad[f * nqp + q_point]; - tmp[f][d] += - jac_grad[r][d] * this->values_quad[k * nqp + q_point]; - } - } - - // Apply J^{-1} appearing in both terms outside braces above - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - VectorizedArrayType res = tmp[0][d] * inv_t_jac[e][0]; - for (unsigned int f = 1; f < dim; ++f) - res += tmp[f][d] * inv_t_jac[e][f]; - grad_out[d][e] = res; - } - - // Add -(J^{-T} * jac_grad * J^{-1} * J * values * det(J^{-1})), - // which can be expressed as a rank-1 update tmp[d] * tmp4[e], where - // tmp = J * values and tmp4 = (J^{-T} * jac_grad * J^{-1}) - VectorizedArrayType tmp3[dim], tmp4[dim]; - for (unsigned int d = 0; d < dim; ++d) - { - tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0]; - for (unsigned int e = 1; e < dim; ++e) - tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e]; - } - for (unsigned int e = 0, k = dim; e < dim; ++e) - for (unsigned int f = e + 1; f < dim; ++k, ++f) - for (unsigned int d = 0; d < dim; ++d) - { - tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d]; - tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d]; - } + for (unsigned int k = 1; k < dim; ++k) + tmp2[i][j] += jac[k][i] * tmp[k][j]; + } + + // diagonal part for (unsigned int d = 0; d < dim; ++d) - { - tmp4[d] = tmp3[0] * inv_t_jac[d][0]; - for (unsigned int e = 1; e < dim; ++e) - tmp4[d] += tmp3[e] * inv_t_jac[d][e]; - } + this->hessians_quad[(comp * hdim + d) * nqp + q_point] = + tmp2[d][d] * JxW; + + // off diagonal part + for (unsigned int d = 0, off_diag = dim; d < dim; ++d) + for (unsigned int e = d + 1; e < dim; ++e, ++off_diag) + this->hessians_quad[(comp * hdim + off_diag) * nqp + q_point] = + (tmp2[d][e] + tmp2[e][d]) * JxW; - VectorizedArrayType tmp2[dim]; + // 3. gradient_unit = J' * hessian for (unsigned int d = 0; d < dim; ++d) { - tmp2[d] = t_jac[0][d] * this->values_quad[q_point]; - for (unsigned e = 1; e < dim; ++e) - tmp2[d] += t_jac[e][d] * this->values_quad[e * nqp + q_point]; + VectorizedArrayType sum = 0; + for (unsigned int e = 0; e < dim; ++e) + sum += hessian_c[e][e] * jac_grad[e][d]; + for (unsigned int e = 0, count = dim; e < dim; ++e) + for (unsigned int f = e + 1; f < dim; ++f, ++count) + sum += + (hessian_c[e][f] + hessian_c[f][e]) * jac_grad[count][d]; + this->gradients_from_hessians_quad[(comp * nqp + q_point) * dim + + d] = sum * JxW; } - - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - grad_out[d][e] -= tmp4[e] * tmp2[d]; - - // finally multiply by det(J^{-1}) necessary in all - // contributions above - grad_out[d][e] *= inv_det; - } } - return grad_out; } +} + + + +template +inline typename FEEvaluationBase::value_type +FEEvaluationBase:: + integrate_value() const +{ +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); + Assert(this->values_quad_submitted == true, + internal::ExcAccessToUninitializedField()); +# endif + + Tensor<1, n_components, VectorizedArrayType> return_value; + const std::size_t nqp = this->n_quadrature_points; + for (unsigned int q = 0; q < nqp; ++q) + for (unsigned int comp = 0; comp < n_components; ++comp) + return_value[comp] += this->values_quad[comp * nqp + q]; + if constexpr (n_components == 1) + return return_value[0]; else - { - return BaseClass::get_gradient(q_point); - } + return return_value; } -template +template +template inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -FEEvaluationAccess< - dim, - dim, - Number, - is_face, - VectorizedArrayType, - std::enable_if_t>::get_divergence(const unsigned int q_point) const +FEEvaluationBase:: + get_divergence(const unsigned int q_point) const { # ifdef DEBUG Assert(this->gradients_quad_initialized == true, @@ -6108,14 +5434,14 @@ FEEvaluationAccess< -template +template +template inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArrayType> - FEEvaluationAccess>:: +FEEvaluationBase:: get_symmetric_gradient(const unsigned int q_point) const { // copy from generic function into dim-specialization function @@ -6148,27 +5474,22 @@ inline DEAL_II_ALWAYS_INLINE SymmetricTensor<2, dim, VectorizedArrayType> -template +template +template inline DEAL_II_ALWAYS_INLINE Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> - FEEvaluationAccess>::get_curl(const unsigned int - q_point) const + FEEvaluationBase:: + get_curl(const unsigned int q_point) const { // copy from generic function into dim-specialization function const Tensor<2, dim, VectorizedArrayType> grad = get_gradient(q_point); Tensor<1, (dim == 2 ? 1 : dim), VectorizedArrayType> curl; switch (dim) { - case 1: - Assert(false, - ExcMessage( - "Computing the curl in 1d is not a useful operation")); - break; case 2: curl[0] = grad[1][0] - grad[0][1]; break; @@ -6185,344 +5506,14 @@ inline DEAL_II_ALWAYS_INLINE -template -inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> - FEEvaluationAccess>:: - get_hessian_diagonal(const unsigned int q_point) const -{ - return BaseClass::get_hessian_diagonal(q_point); -} - - - -template -inline DEAL_II_ALWAYS_INLINE Tensor<3, dim, VectorizedArrayType> - FEEvaluationAccess>::get_hessian(const unsigned int - q_point) const -{ -# ifdef DEBUG - Assert(this->hessians_quad_initialized == true, - internal::ExcAccessToUninitializedField()); -# endif - AssertIndexRange(q_point, this->n_quadrature_points); - return BaseClass::get_hessian(q_point); -} - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: - submit_value(const Tensor<1, dim, VectorizedArrayType> val_in, - const unsigned int q_point) -{ - if (this->data->element_type == - internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) - { - // Piola transform is required - AssertIndexRange(q_point, this->n_quadrature_points); - Assert(this->J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_value")); -# ifdef DEBUG - Assert(this->is_reinitialized, ExcNotInitialized()); - this->values_quad_submitted = true; -# endif - - VectorizedArrayType *values = this->values_quad + q_point; - const std::size_t nqp = this->n_quadrature_points; - - if (!is_face && - this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - const Tensor<2, dim, VectorizedArrayType> jac = this->jacobian[1]; - const VectorizedArrayType weight = this->quadrature_weights[q_point]; - - for (unsigned int comp = 0; comp < n_components; ++comp) - values[comp * nqp] = val_in[comp] * weight * jac[comp][comp]; - } - else - { - // Affine or general cell - const Tensor<2, dim, VectorizedArrayType> inv_t_jac = - (this->cell_type > internal::MatrixFreeFunctions::affine) ? - this->jacobian[q_point] : - this->jacobian[0]; - - // Derivatives are reordered for faces. Need to take this into account - // and 1/inv_det != J_value for faces - const VectorizedArrayType fac = - (!is_face) ? - this->quadrature_weights[q_point] : - (((this->cell_type > internal::MatrixFreeFunctions::affine) ? - this->J_value[q_point] : - this->J_value[0] * this->quadrature_weights[q_point]) * - ((dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac))); - const Tensor<2, dim, VectorizedArrayType> jac = - (this->cell_type > internal::MatrixFreeFunctions::affine) ? - transpose(invert(inv_t_jac)) : - this->jacobian[1]; - - // J^T * u * factor - for (unsigned int comp = 0; comp < n_components; ++comp) - { - values[comp * nqp] = val_in[0] * jac[0][comp]; - for (unsigned int e = 1; e < dim; ++e) - values[comp * nqp] += val_in[e] * jac[e][comp]; - values[comp * nqp] *= fac; - } - } - } - else - { - // No Piola transform - BaseClass::submit_value(val_in, q_point); - } -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: - submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in, - const unsigned int q_point) -{ - if (this->data->element_type == - internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) - { - // Piola transform is required - -# ifdef DEBUG - Assert(this->is_reinitialized, ExcNotInitialized()); -# endif - AssertIndexRange(q_point, this->n_quadrature_points); - Assert(this->J_value != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); - Assert(this->jacobian != nullptr, - internal::ExcMatrixFreeAccessToUninitializedMappingField( - "update_gradients")); -# ifdef DEBUG - this->gradients_quad_submitted = true; -# endif - - VectorizedArrayType *gradients = this->gradients_quad + q_point * dim; - VectorizedArrayType *values = this->values_from_gradients_quad + q_point; - const std::size_t nqp = this->n_quadrature_points; - const std::size_t nqp_d = nqp * dim; - - if (!is_face && - this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - // Cartesian cell - const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = - this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; - const VectorizedArrayType weight = this->quadrature_weights[q_point]; - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int comp = 0; comp < n_components; ++comp) - gradients[comp * nqp_d + d] = - grad_in[comp][d] * inv_t_jac[d][d] * (jac[comp][comp] * weight); - } - else if (this->cell_type <= internal::MatrixFreeFunctions::affine) - { - // Affine cell - const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = - this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; - - // Derivatives are reordered for faces. Need to take this into account - // and 1/inv_det != J_value for faces - const VectorizedArrayType fac = - (!is_face) ? this->quadrature_weights[q_point] : - this->J_value[0] * this->quadrature_weights[q_point] * - ((dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac)); - - // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor - VectorizedArrayType tmp[dim][dim]; - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - tmp[d][e] = inv_t_jac[0][d] * grad_in[e][0]; - for (unsigned int f = 1; f < dim; ++f) - tmp[d][e] += inv_t_jac[f][d] * grad_in[e][f]; - } - for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType res = jac[0][comp] * tmp[d][0]; - for (unsigned int f = 1; f < dim; ++f) - res += jac[f][comp] * tmp[d][f]; - - gradients[comp * nqp_d + d] = res * fac; - } - } - else - { - // General cell - - const auto jac_grad = this->jacobian_gradients_non_inverse[q_point]; - const Tensor<2, dim, VectorizedArrayType> inv_t_jac = - this->jacobian[q_point]; - - // Derivatives are reordered for faces. Need to take this into account - // and 1/inv_det != J_value for faces - const VectorizedArrayType fac = - (!is_face) ? - this->quadrature_weights[q_point] : - this->J_value[q_point] * ((dim == 2 && this->get_face_no() < 2) ? - -determinant(inv_t_jac) : - determinant(inv_t_jac)); - const Tensor<2, dim, VectorizedArrayType> t_jac = invert(inv_t_jac); - - // Start evaluation for values part below to enable the compiler to - // possibly re-use the same computation in get_gradient() without - // interfering with stores to 'gradients' - VectorizedArrayType tmp3[dim], tmp4[dim]; - for (unsigned int d = 0; d < dim; ++d) - { - tmp3[d] = inv_t_jac[0][d] * jac_grad[d][0]; - for (unsigned int e = 1; e < dim; ++e) - tmp3[d] += inv_t_jac[e][d] * jac_grad[d][e]; - } - for (unsigned int e = 0, k = dim; e < dim; ++e) - for (unsigned int f = e + 1; f < dim; ++k, ++f) - for (unsigned int d = 0; d < dim; ++d) - { - tmp3[f] += inv_t_jac[d][e] * jac_grad[k][d]; - tmp3[e] += inv_t_jac[d][f] * jac_grad[k][d]; - } - for (unsigned int d = 0; d < dim; ++d) - { - tmp4[d] = tmp3[0] * inv_t_jac[d][0]; - for (unsigned int e = 1; e < dim; ++e) - tmp4[d] += tmp3[e] * inv_t_jac[d][e]; - } - - const Tensor<2, dim, VectorizedArrayType> grad_in_scaled = - fac * grad_in; - - VectorizedArrayType tmp[dim][dim]; - - // J * (J^{-1} * (grad_in * factor)) - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - tmp[d][e] = inv_t_jac[0][d] * grad_in_scaled[e][0]; - for (unsigned int f = 1; f < dim; ++f) - tmp[d][e] += inv_t_jac[f][d] * grad_in_scaled[e][f]; - } - - for (unsigned int d = 0; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - { - VectorizedArrayType res = t_jac[d][0] * tmp[e][0]; - for (unsigned int f = 1; f < dim; ++f) - res += t_jac[d][f] * tmp[e][f]; - - gradients[d * nqp_d + e] = res; - } - - // jac_grad * (J^{-1} * (grad_in * factor)), re-use part in braces - // as 'tmp' from above - VectorizedArrayType value[dim]; - for (unsigned int d = 0; d < dim; ++d) - { - value[d] = tmp[d][0] * jac_grad[d][0]; - for (unsigned int e = 1; e < dim; ++e) - value[d] += tmp[d][e] * jac_grad[d][e]; - } - for (unsigned int e = 0, k = dim; e < dim; ++e) - for (unsigned int f = e + 1; f < dim; ++k, ++f) - for (unsigned int d = 0; d < dim; ++d) - { - value[e] += tmp[f][d] * jac_grad[k][d]; - value[f] += tmp[e][d] * jac_grad[k][d]; - } - - // -(grad_in * factor) * J * (J^{-T} * jac_grad * J^{-1}) - // = -(grad_in * factor) * J * ( \------- tmp4 ---------/ ) - for (unsigned int d = 0; d < dim; ++d) - { - VectorizedArrayType tmp2 = grad_in_scaled[d][0] * tmp4[0]; - for (unsigned int e = 1; e < dim; ++e) - tmp2 += grad_in_scaled[d][e] * tmp4[e]; - for (unsigned int e = 0; e < dim; ++e) - value[e] -= t_jac[e][d] * tmp2; - } - - for (unsigned int d = 0; d < dim; ++d) - values[d * nqp] = value[d]; - } - } - else - { - BaseClass::submit_gradient(grad_in, q_point); - } -} - - - -template -inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: - submit_gradient( - const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in, - const unsigned int q_point) -{ - if (this->data->element_type == - internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) - { - // Piola transform is required - const Tensor<2, dim, VectorizedArrayType> grad = grad_in; - FEEvaluationAccess:: - submit_gradient(grad, q_point); - } - else - { - BaseClass::submit_gradient(grad_in, q_point); - } -} - - - -template +template +template inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: +FEEvaluationBase:: submit_divergence(const VectorizedArrayType div_in, const unsigned int q_point) { @@ -6607,14 +5598,14 @@ FEEvaluationAccess +template +template inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: +FEEvaluationBase:: submit_symmetric_gradient( const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad, const unsigned int q_point) @@ -6695,25 +5686,20 @@ FEEvaluationAccess +template +template inline DEAL_II_ALWAYS_INLINE void -FEEvaluationAccess>:: +FEEvaluationBase:: submit_curl(const Tensor<1, dim == 2 ? 1 : dim, VectorizedArrayType> curl, const unsigned int q_point) { Tensor<2, dim, VectorizedArrayType> grad; switch (dim) { - case 1: - Assert(false, - ExcMessage( - "Testing by the curl in 1d is not a useful operation")); - break; case 2: grad[1][0] = curl[0]; grad[0][1] = -curl[0]; @@ -6763,7 +5749,8 @@ inline FEEvaluationdata->dofs_per_component_on_cell) , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_) , n_q_points(this->data->n_q_points) -- 2.39.5