From 09aa10961063fc8f12ed35542e96aa39ec533d21 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 10 Sep 2020 16:57:31 +0200 Subject: [PATCH] Split out mapping-related part from FEEvaluationBase into base class --- .../deal.II/matrix_free/evaluation_kernels.h | 17 +- include/deal.II/matrix_free/fe_evaluation.h | 1849 +++++++++-------- 2 files changed, 966 insertions(+), 900 deletions(-) diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 19055df97d..57513e8977 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -32,6 +32,11 @@ DEAL_II_NAMESPACE_OPEN +// forward declaration +template +class FEEvaluationBaseData; + + namespace internal { @@ -3165,12 +3170,14 @@ namespace internal template struct CellwiseInverseMassMatrixImpl { - template static void - apply(const unsigned int n_components, - const FEEvaluationType &fe_eval, - const Number * in_array, - Number * out_array) + apply(const unsigned int n_components, + const FEEvaluationBaseData &fe_eval, + const Number * in_array, + Number * out_array) { constexpr unsigned int dofs_per_component = Utilities::pow(fe_degree + 1, dim); diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 4e4d7db727..38bb982967 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -60,26 +60,22 @@ class FEEvaluation; /** - * This is the base class for the FEEvaluation classes. This class is a base - * class and needs usually not be called in user code. It does not have any - * public constructor. The usage is through the class FEEvaluation instead. It - * implements a reinit method that is used to set pointers so that operations - * on quadrature points can be performed quickly, access functions to vectors - * for the @p read_dof_values, @p set_dof_values, and @p - * distribute_local_to_global functions, as well as methods to access values - * and gradients of finite element functions. + * This base class of the FEEvaluation and FEFaceEvaluation classes handles + * mapping-related information independent of the degrees of freedom and + * finite element in use. This class provides access functionality for user + * code but is otherwise invisible without any public constructor. The usage + * is through the class FEEvaluation instead. * - * This class has three template arguments: + * This class has four template arguments: * * @tparam dim Dimension in which this class is to be used * - * @tparam n_components Number of vector components when solving a system of - * PDEs. If the same operation is applied to several components of a PDE (e.g. - * a vector Laplace equation), they can be applied simultaneously with one - * call (and often more efficiently) - * * @tparam Number Number format, usually @p double or @p float * + * @tparam is_face Whether the class is used for a cell integrator (with + * quadrature dimension the same as the space dimension) or for a face + * integrator (with quadrature dimension one less) + * * @tparam VectorizedArrayType Type of array to be woked on in a vectorized * fashion, defaults to VectorizedArray * @@ -90,32 +86,22 @@ class FEEvaluation; * @ingroup matrixfree */ template > -class FEEvaluationBase +class FEEvaluationBaseData { static_assert( std::is_same::value, "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_; + static constexpr unsigned int dimension = dim; - /** - * @name 1: General operations - */ - //@{ /** * Destructor. */ - ~FEEvaluationBase(); + ~FEEvaluationBaseData(); /** * Return the index offset within the geometry fields for the cell the @p @@ -143,10 +129,365 @@ public: const internal::MatrixFreeFunctions::ShapeInfo & get_shape_info() const; - //@} + /** + * Return the determinant of the Jacobian from the unit to the real cell + * times the quadrature weight. + */ + VectorizedArrayType + JxW(const unsigned int q_point) const; + + /** + * Return the inverse and transposed version of Jacobian of the mapping + * between the unit to the real cell (representing the covariant + * transformation). This is exactly the matrix used internally to transform + * the unit cell gradients to gradients on the real cell. + */ + Tensor<2, dim, VectorizedArrayType> + inverse_jacobian(const unsigned int q_point) const; + + /** + * Return the unit normal vector on a face. Note that both sides of a face + * use the same orientation of the normal vector: For the faces enumerated + * as `interior` in FaceToCellTopology and selected with the + * `is_interior_face=true` flag of the constructor, this corresponds to the + * outer normal vector, whereas for faces enumerated as `exterior` in + * FaceToCellTopology and selected with the `is_interior_face=false` flag of + * the constructor, the normal points into the element as a consequence of + * the single normal vector. + * + * @note Only implemented in case `is_face == true`. + */ + Tensor<1, dim, VectorizedArrayType> + get_normal_vector(const unsigned int q_point) const; + + /** + * Provides a unified interface to access data in a vector of + * VectorizedArray fields of length MatrixFree::n_macro_cells() + + * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces + * (indirect addressing). + */ + VectorizedArrayType + read_cell_data(const AlignedVector &array) const; + + /** + * Provides a unified interface to set data in a vector of + * VectorizedArray fields of length MatrixFree::n_macro_cells() + + * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces + * (indirect addressing). + */ + void + set_cell_data(AlignedVector &array, + const VectorizedArrayType & value) const; + + /** + * The same as above, just for std::array of length of VectorizedArrayType for + * arbitrary data type. + */ + template + std::array + read_cell_data(const AlignedVector> + &array) const; + + /** + * The same as above, just for std::array of length of VectorizedArrayType for + * arbitrary data type. + */ + template + void + set_cell_data( + AlignedVector> &array, + const std::array & value) const; + + /** + * Return the id of the cells this FEEvaluation or FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_ids() const; + + /** + * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is + * associated with. + */ + std::array + get_cell_or_face_ids() const; + + + /** + * Return the numbering of local degrees of freedom within the evaluation + * routines of FEEvaluation in terms of the standard numbering on finite + * elements. + */ + const std::vector & + get_internal_dof_numbering() const; + + /** + * Return an ArrayView to internal memory for temporary use. Note that some + * of this memory is overwritten during evaluate() and integrate() calls so + * do not assume it to be stable over those calls. The maximum size you can + * write into is 3*dofs_per_cell+2*n_q_points. + */ + ArrayView + get_scratch_data() const; + +protected: + /** + * Constructor. Made protected to prevent users from directly using this + * class. Takes all data stored in MatrixFree. If applied to problems with + * more than one quadrature formula selected during construction of + * `matrix_free`, `quad_no` allows to select the appropriate formula. + */ + FEEvaluationBaseData( + 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); + + /** + * Constructor that comes with reduced functionality and works similar as + * FEValues. + */ + FEEvaluationBaseData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other); + + /** + * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, + * quadrature, and update flags, the underlying geometry evaluation based on + * FEValues will be deep-copied in order to allow for using in parallel with + * threads. + */ + FEEvaluationBaseData(const FEEvaluationBaseData &other); + + /** + * Copy assignment operator. If FEEvaluationBase was constructed from a + * mapping, fe, quadrature, and update flags, the underlying geometry + * evaluation based on FEValues will be deep-copied in order to allow for + * using in parallel with threads. + */ + FEEvaluationBaseData & + operator=(const FEEvaluationBaseData &other); + + /** + * This is the general array for all data fields. + */ + AlignedVector *scratch_data_array; + + /** + * This is the user-visible part of scratch_data_array, only showing the + * last part of scratch_data_array. The first part is consumed by + * values_dofs, values_quad, etc. + */ + VectorizedArrayType *scratch_data; + + /** + * The number of the quadrature formula of the present cell. + */ + const unsigned int quad_no; /** - * @name 2: Reading from and writing to vectors + * The active fe index for this class for efficient indexing in the hp case. + */ + const unsigned int active_fe_index; + + /** + * The active quadrature index for this class for efficient indexing in the + * hp case. + */ + const unsigned int active_quad_index; + + /** + * The number of quadrature points in the current evaluation context. + */ + const unsigned int n_quadrature_points; + + /** + * A pointer to the underlying data. + */ + const MatrixFree *matrix_info; + + /** + * A pointer to the underlying DoF indices and constraint description + * for the component specified at construction. Also contained in + * matrix_info, but it simplifies code if we store a reference to it. + */ + const internal::MatrixFreeFunctions::DoFInfo *dof_info; + + /** + * A pointer to the underlying transformation data from unit to real cells + * for the given quadrature formula specified at construction. Also + * contained in matrix_info, but it simplifies code if we store a reference + * to it. + */ + const internal::MatrixFreeFunctions::MappingInfoStorage< + (is_face ? dim - 1 : dim), + dim, + Number, + VectorizedArrayType> *mapping_data; + + /** + * A pointer to the unit cell shape data, i.e., values, gradients and + * Hessians in 1D at the quadrature points that constitute the tensor + * product. Also contained in matrix_info, but it simplifies code if we + * store a reference to it. + */ + const internal::MatrixFreeFunctions::ShapeInfo *data; + + /** + * A pointer to the Jacobian information of the present cell. Only set to a + * useful value if on a non-Cartesian cell. + */ + const Tensor<2, dim, VectorizedArrayType> *jacobian; + + /** + * A pointer to the Jacobian determinant of the present cell. If on a + * Cartesian cell or on a cell with constant Jacobian, this is just the + * Jacobian determinant, otherwise the Jacobian determinant times the + * quadrature weight. + */ + const VectorizedArrayType *J_value; + + /** + * A pointer to the normal vectors at faces. + */ + const Tensor<1, dim, VectorizedArrayType> *normal_vectors; + + /** + * A pointer to the normal vectors times the jacobian at faces. + */ + const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian; + + /** + * A pointer to the quadrature weights of the underlying quadrature formula. + */ + const Number *quadrature_weights; + + /** + * After a call to reinit(), stores the number of the cell we are currently + * working with. + */ + unsigned int cell; + + /** + * Flag holding information whether a face is an interior or exterior face + * according to the defined direction of the normal. Not used for cells. + */ + bool is_interior_face; + + /** + * Stores the index an FEFaceEvaluation object is currently pointing into + * (interior face, exterior face, data associated with cell). + */ + internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; + + /** + * Stores the current number of a face within the given cell in case + * `is_face==true`, using values between `0` and `2*dim`. + */ + unsigned int face_no; + + /** + * Stores the orientation of the given face with respect to the standard + * orientation, 0 if in standard orientation. + */ + unsigned int face_orientation; + + /** + * Stores the subface index of the given face. Usually, this variable takes + * the value numbers::invalid_unsigned_int to indicate integration over the + * full face, but in case the current physical face has a neighbor that is + * more refined, it is a subface and must scale the entries in ShapeInfo + * appropriately. + */ + unsigned int subface_index; + + /** + * Stores the type of the cell we are currently working with after a call to + * reinit(). Valid values are @p cartesian, @p affine and @p general, which + * have different implications on how the Jacobian transformations are + * stored internally in MappingInfo. + */ + internal::MatrixFreeFunctions::GeometryType cell_type; + + /** + * Geometry data that can be generated FEValues on the fly with the + * respective constructor. + */ + std::shared_ptr> + mapped_geometry; + + // Make FEEvaluation objects friends for access to protected member + // mapped_geometry. + template + friend class FEEvaluation; +}; + + + +/** + * This is the base class for the FEEvaluation classes. This class needs + * usually not be called in user code and does not have any public + * constructor. The usage is through the class FEEvaluation instead. It + * implements a reinit method that is used to set pointers so that operations + * on quadrature points can be performed quickly, access functions to vectors + * for the FEEvaluationBase::read_dof_values(), + * FEEvaluationBase::set_dof_values(), and + * FEEvaluationBase::distribute_local_to_global() functions, as well as + * methods to access values and gradients of finite element functions. It also + * inherits the geometry access functions provided by the class + * FEEvaluationBaseData. + * + * This class has five template arguments: + * + * @tparam dim Dimension in which this class is to be used + * + * @tparam n_components Number of vector components when solving a system of + * PDEs. If the same operation is applied to several components of a PDE (e.g. + * a vector Laplace equation), they can be applied simultaneously with one + * call (and often more efficiently) + * + * @tparam Number Number format, usually @p double or @p float + * + * @tparam is_face Whether the class is used for a cell integrator (with + * quadrature dimension the same as the space dimension) or for a face + * integrator (with quadrature dimension one less) + * + * @tparam VectorizedArrayType Type of array to be woked on in a vectorized + * fashion, defaults to VectorizedArray + * + * @note Currently only VectorizedArray is supported as + * VectorizedArrayType. + * + * + * @ingroup matrixfree + */ +template > +class FEEvaluationBase + : public FEEvaluationBaseData +{ +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_; + + /** + * @name 1: Reading from and writing to vectors */ //@{ /** @@ -305,7 +646,7 @@ public: //@} /** - * @name 3: Data access + * @name 2: Access to data at quadrature points or the gather vector data */ //@{ /** @@ -480,7 +821,7 @@ public: #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: + // For now, hack in those functions manually only to fix documentation: /** * Return the divergence of a vector-valued finite element at quadrature @@ -570,7 +911,7 @@ public: #endif /** - * Takes values on quadrature points, multiplies by the Jacobian determinant + * Takes values at quadrature points, multiplies by the Jacobian determinant * and quadrature weights (JxW) and sums the values for all quadrature * points on the cell. The result is a scalar, representing the integral * over the function over the cell. If a vector-element is used, the @@ -580,102 +921,18 @@ public: * * @note In case the FEEvaluation object is initialized with a batch of * cells where not all lanes in the SIMD vector VectorizedArray are - * representing actual data, this method performs computations on dummy data - * (that is copied from the last valid lane) that will not make sense. Thus, - * the user needs to make sure that it is not used in any computation - * explicitly, like when summing the results of several cells. - */ - value_type - integrate_value() const; - - /** - * Return the determinant of the Jacobian from the unit to the real cell - * times the quadrature weight. - */ - VectorizedArrayType - JxW(const unsigned int q_point) const; - - /** - * Return the inverse and transposed version of Jacobian of the mapping - * between the unit to the real cell (representing the covariant - * transformation). This is exactly the matrix used internally to transform - * the unit cell gradients to gradients on the real cell. - */ - Tensor<2, dim, VectorizedArrayType> - inverse_jacobian(const unsigned int q_point) const; - - /** - * Return the unit normal vector on a face. Note that both sides of a face - * use the same orientation of the normal vector: For the faces enumerated - * as `interior` in FaceToCellTopology and selected with the - * `is_interior_face=true` flag of the constructor, this corresponds to the - * outer normal vector, whereas for faces enumerated as `exterior` in - * FaceToCellTopology and selected with the `is_interior_face=false` flag of - * the constructor, the normal points into the element as a consequence of - * the single normal vector. - * - * @note Only implemented in case `is_face == true`. - */ - Tensor<1, dim, VectorizedArrayType> - get_normal_vector(const unsigned int q_point) const; - - /** - * Provides a unified interface to access data in a vector of - * VectorizedArray fields of length MatrixFree::n_macro_cells() + - * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces - * (indirect addressing). - */ - VectorizedArrayType - read_cell_data(const AlignedVector &array) const; - - /** - * Provides a unified interface to set data in a vector of - * VectorizedArray fields of length MatrixFree::n_macro_cells() + - * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces - * (indirect addressing). - */ - void - set_cell_data(AlignedVector &array, - const VectorizedArrayType & value) const; - - /** - * The same as above, just for std::array of length of VectorizedArrayType for - * arbitrary data type. - */ - template - std::array - read_cell_data(const AlignedVector> - &array) const; - - /** - * The same as above, just for std::array of length of VectorizedArrayType for - * arbitrary data type. - */ - template - void - set_cell_data( - AlignedVector> &array, - const std::array & value) const; - - /** - * Return the id of the cells this FEEvaluation or FEFaceEvaluation is - * associated with. - */ - std::array - get_cell_ids() const; - - - /** - * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is - * associated with. + * representing actual data, this method performs computations on dummy data + * (that is copied from the last valid lane) that will not make sense. Thus, + * the user needs to make sure that it is not used in any computation + * explicitly, like when summing the results of several cells. */ - std::array - get_cell_or_face_ids() const; + value_type + integrate_value() const; //@} /** - * @name 4: Access to internal data + * @name 3: Access to internal data */ //@{ /** @@ -784,23 +1041,6 @@ public: VectorizedArrayType * begin_hessians(); - /** - * Return the numbering of local degrees of freedom within the evaluation - * routines of FEEvaluation in terms of the standard numbering on finite - * elements. - */ - const std::vector & - get_internal_dof_numbering() const; - - /** - * Return an ArrayView to internal memory for temporary use. Note that some - * of this memory is overwritten during evaluate() and integrate() calls so - * do not assume it to be stable over those calls. The maximum size you can - * write into is 3*dofs_per_cell+2*n_q_points. - */ - ArrayView - get_scratch_data() const; - //@} protected: @@ -847,25 +1087,23 @@ protected: * possibly within the element if the evaluate/integrate routines are * combined inside user code (e.g. for computing cell matrices). * - * The optional FEEvaluationBase object allows several FEEvaluation objects - * to share the geometry evaluation, i.e., the underlying mapping and - * quadrature points do only need to be evaluated once. This only works if - * the quadrature formulas are the same. Otherwise, a new evaluation object - * is created. Make sure to not pass an optional object around when you - * intend to use the FEEvaluation object in %parallel with another one - * because otherwise the intended sharing may create race conditions. - */ - template - FEEvaluationBase(const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase *other); + * The optional FEEvaluationBaseData object allows several + * FEEvaluation objects to share the geometry evaluation, i.e., the + * underlying mapping and quadrature points do only need to be evaluated + * once. This only works if the quadrature formulas are the same. Otherwise, + * a new evaluation object is created. Make sure to not pass an optional + * object around when you intend to use the FEEvaluation object in %parallel + * with another one because otherwise the intended sharing may create race + * conditions. + */ + FEEvaluationBase( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other); /** * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, @@ -923,18 +1161,6 @@ protected: read_write_operation_global(const VectorOperation &operation, VectorType * vectors[]) const; - /** - * This is the general array for all data fields. - */ - AlignedVector *scratch_data_array; - - /** - * This is the user-visible part of scratch_data_array, only showing the - * last part of scratch_data_array. The first part is consumed by - * values_dofs, values_quad, etc. - */ - VectorizedArrayType *scratch_data; - /** * This field stores the values for local degrees of freedom (e.g. after * reading out from a vector but before applying unit cell transformations @@ -990,142 +1216,12 @@ protected: */ VectorizedArrayType *hessians_quad; - /** - * Stores the number of the quadrature formula of the present cell. - */ - const unsigned int quad_no; - /** * Stores the number of components in the finite element as detected in the * MatrixFree storage class for comparison with the template argument. */ const unsigned int n_fe_components; - /** - * Stores the active fe index for this class for efficient indexing in the - * hp case. - */ - const unsigned int active_fe_index; - - /** - * Stores the active quadrature index for this class for efficient indexing - * in the hp case. - */ - const unsigned int active_quad_index; - - /** - * Stores the number of quadrature points in the current evaluation context. - */ - const unsigned int n_quadrature_points; - - /** - * Stores a pointer to the underlying data. - */ - const MatrixFree *matrix_info; - - /** - * Stores a pointer to the underlying DoF indices and constraint description - * for the component specified at construction. Also contained in - * matrix_info, but it simplifies code if we store a reference to it. - */ - const internal::MatrixFreeFunctions::DoFInfo *dof_info; - - /** - * Stores a pointer to the underlying transformation data from unit to - * real cells for the given quadrature formula specified at construction. - * Also contained in matrix_info, but it simplifies code if we store a - * reference to it. - */ - const internal::MatrixFreeFunctions::MappingInfoStorage< - (is_face ? dim - 1 : dim), - dim, - Number, - VectorizedArrayType> *mapping_data; - - /** - * Stores a pointer to the unit cell shape data, i.e., values, gradients and - * Hessians in 1D at the quadrature points that constitute the tensor - * product. Also contained in matrix_info, but it simplifies code if we - * store a reference to it. - */ - const internal::MatrixFreeFunctions::ShapeInfo *data; - - /** - * A pointer to the Jacobian information of the present cell. Only set to a - * useful value if on a non-Cartesian cell. - */ - const Tensor<2, dim, VectorizedArrayType> *jacobian; - - /** - * A pointer to the Jacobian determinant of the present cell. If on a - * Cartesian cell or on a cell with constant Jacobian, this is just the - * Jacobian determinant, otherwise the Jacobian determinant times the - * quadrature weight. - */ - const VectorizedArrayType *J_value; - - /** - * A pointer to the normal vectors at faces. - */ - const Tensor<1, dim, VectorizedArrayType> *normal_vectors; - - /** - * A pointer to the normal vectors times the jacobian at faces. - */ - const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian; - - /** - * A pointer to the quadrature weights of the underlying quadrature formula. - */ - const Number *quadrature_weights; - - /** - * After a call to reinit(), stores the number of the cell we are currently - * working with. - */ - unsigned int cell; - - /** - * Flag holding information whether a face is an interior or exterior face - * according to the defined direction of the normal. Not used for cells. - */ - bool is_interior_face; - - /** - * Stores the index an FEFaceEvaluation object is currently pointing into - * (interior face, exterior face, data associated with cell). - */ - internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index; - - /** - * Stores the current number of a face within the given cell in case - * `is_face==true`, using values between `0` and `2*dim`. - */ - unsigned int face_no; - - /** - * Stores the orientation of the given face with respect to the standard - * orientation, 0 if in standard orientation. - */ - unsigned int face_orientation; - - /** - * Stores the subface index of the given face. Usually, this variable takes - * the value numbers::invalid_unsigned_int to indicate integration over the - * full face, but in case the current physical face has a neighbor that is - * more refined, it is a subface and must scale the entries in ShapeInfo - * appropriately. - */ - unsigned int subface_index; - - /** - * Stores the type of the cell we are currently working with after a call to - * reinit(). Valid values are @p cartesian, @p affine and @p general, which - * have different implications on how the Jacobian transformations are - * stored internally in MappingInfo. - */ - internal::MatrixFreeFunctions::GeometryType cell_type; - /** * Debug information to track whether dof values have been initialized * before accessed. Used to control exceptions when uninitialized data is @@ -1168,14 +1264,6 @@ protected: */ bool gradients_quad_submitted; - /** - * Geometry data that can be generated FEValues on the fly with the - * respective constructor. - */ - std::shared_ptr> - mapped_geometry; - /** * For a FiniteElement with more than one base element, select at which * component this data structure should start. @@ -1191,16 +1279,10 @@ protected: private: /** * Sets the pointers for values, gradients, hessians to the central - * scratch_data_array. + * scratch_data_array of the base class. */ void set_data_pointers(); - - // Make other FEEvaluationBase as well as FEEvaluation objects friends. - template - friend class FEEvaluationBase; - template - friend class FEEvaluation; }; @@ -1258,17 +1340,14 @@ protected: * Constructor with reduced functionality for similar usage of FEEvaluation * as FEValues, including matrix assembly. */ - template - FEEvaluationAccess(const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase *other); + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other); /** * Copy constructor @@ -1396,17 +1475,14 @@ protected: * Constructor with reduced functionality for similar usage of FEEvaluation * as FEValues, including matrix assembly. */ - template - FEEvaluationAccess(const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase *other); + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other); /** * Copy constructor @@ -1558,17 +1634,14 @@ protected: * Constructor with reduced functionality for similar usage of FEEvaluation * as FEValues, including matrix assembly. */ - template - FEEvaluationAccess(const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase *other); + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other); /** * Copy constructor @@ -1707,17 +1780,13 @@ protected: * Constructor with reduced functionality for similar usage of FEEvaluation * as FEValues, including matrix assembly. */ - template - FEEvaluationAccess(const Mapping<1> & mapping, - const FiniteElement<1> &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase<1, - n_components_other, - Number, - is_face, - VectorizedArrayType> *other); + FEEvaluationAccess( + const Mapping<1> & mapping, + const FiniteElement<1> &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other); /** * Copy constructor @@ -2450,14 +2519,10 @@ public: * use the FEEvaluation object in %parallel to the given one because * otherwise the intended sharing may create race conditions. */ - template - FEEvaluation(const FiniteElement & fe, - const FEEvaluationBase &other, - const unsigned int first_selected_component = 0); + FEEvaluation( + const FiniteElement & fe, + const FEEvaluationBaseData &other, + const unsigned int first_selected_component = 0); /** * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe, @@ -3107,29 +3172,20 @@ namespace internal #ifndef DOXYGEN +/*----------------------- FEEvaluationBaseData ------------------------*/ -/*----------------------- FEEvaluationBase ----------------------------------*/ - -template -inline FEEvaluationBase:: - FEEvaluationBase(const MatrixFree &data_in, - const unsigned int dof_no, - const unsigned int first_selected_component, - const unsigned int quad_no_in, - const unsigned int fe_degree, - const unsigned int n_q_points, - const bool is_interior_face) +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData( + const MatrixFree &data_in, + const unsigned int dof_no, + const unsigned int first_selected_component, + const unsigned int quad_no_in, + const unsigned int fe_degree, + const unsigned int n_q_points, + const bool is_interior_face) : scratch_data_array(data_in.acquire_scratch_data()) , quad_no(quad_no_in) - , n_fe_components(data_in.get_dof_info(dof_no).start_components.back()) , active_fe_index(fe_degree != numbers::invalid_unsigned_int ? data_in.get_dof_info(dof_no).fe_index_from_degree( first_selected_component, @@ -3185,15 +3241,7 @@ inline FEEvaluationBasemapping_initialized() == true, ExcNotInitialized()); AssertDimension(matrix_info->get_task_info().vectorization_length, VectorizedArrayType::size()); @@ -3201,59 +3249,22 @@ inline FEEvaluationBasedescriptor[active_quad_index].n_q_points); - Assert( - dof_info->start_components.back() == 1 || - static_cast(n_components_) <= - static_cast( - dof_info->start_components - [dof_info->component_to_base_index[first_selected_component] + 1]) - - first_selected_component, - ExcMessage( - "You tried to construct a vector-valued evaluator with " + - std::to_string(n_components) + - " components. However, " - "the current base element has only " + - std::to_string( - dof_info->start_components - [dof_info->component_to_base_index[first_selected_component] + 1] - - first_selected_component) + - " components left when starting from local element index " + - std::to_string( - first_selected_component - - dof_info->start_components - [dof_info->component_to_base_index[first_selected_component]]) + - " (global index " + std::to_string(first_selected_component) + ")")); - - // do not check for correct dimensions of data fields here, should be done - // in derived classes } -template -template -inline FEEvaluationBase:: - FEEvaluationBase(const Mapping & mapping, - const FiniteElement &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase *other) +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other) : scratch_data_array(new AlignedVector()) , quad_no(numbers::invalid_unsigned_int) - , n_fe_components(n_components_) , active_fe_index(numbers::invalid_unsigned_int) , active_quad_index(numbers::invalid_unsigned_int) , n_quadrature_points( @@ -3271,24 +3282,12 @@ inline FEEvaluationBasemapped_geometry.get() != nullptr, ExcInternalError()); if (other != nullptr && @@ -3304,40 +3303,19 @@ inline FEEvaluationBaseget_data_storage(); jacobian = mapped_geometry->get_data_storage().jacobians[0].begin(); J_value = mapped_geometry->get_data_storage().JxW_values.begin(); - - const unsigned int base_element_number = - fe.component_to_base_index(first_selected_component).first; - Assert(fe.element_multiplicity(base_element_number) == 1 || - fe.element_multiplicity(base_element_number) - - first_selected_component >= - n_components_, - ExcMessage("The underlying element must at least contain as many " - "components as requested by this class")); - (void)base_element_number; } -template -inline FEEvaluationBase:: - FEEvaluationBase(const FEEvaluationBase &other) +template +inline FEEvaluationBaseData:: + FEEvaluationBaseData( + const FEEvaluationBaseData + &other) : scratch_data_array(other.matrix_info == nullptr ? new AlignedVector() : other.matrix_info->acquire_scratch_data()) , quad_no(other.quad_no) - , n_fe_components(other.n_fe_components) , active_fe_index(other.active_fe_index) , active_quad_index(other.active_quad_index) , n_quadrature_points(other.n_quadrature_points) @@ -3360,16 +3338,7 @@ inline FEEvaluationBase -inline FEEvaluationBase & -FEEvaluationBase:: -operator=(const FEEvaluationBase &other) +template +inline FEEvaluationBaseData & +FEEvaluationBaseData::operator=( + const FEEvaluationBaseData &other) { AssertDimension(quad_no, other.quad_no); - AssertDimension(n_fe_components, other.n_fe_components); AssertDimension(active_fe_index, other.active_fe_index); AssertDimension(active_quad_index, other.active_quad_index); - AssertDimension(first_selected_component, other.first_selected_component); // release old memory if (matrix_info == nullptr) @@ -3437,7 +3392,6 @@ operator=(const FEEvaluationBaseacquire_scratch_data(); } - set_data_pointers(); quadrature_weights = (mapping_data != nullptr ? @@ -3468,16 +3422,9 @@ operator=(const FEEvaluationBase -inline FEEvaluationBase::~FEEvaluationBase() +template +inline FEEvaluationBaseData:: + ~FEEvaluationBaseData() { if (matrix_info != nullptr) { @@ -3499,61 +3446,9 @@ inline FEEvaluationBase -inline void -FEEvaluationBase:: - set_data_pointers() -{ - Assert(scratch_data_array != nullptr, ExcInternalError()); - - const unsigned int tensor_dofs_per_component = - Utilities::fixed_power(this->data->data.front().fe_degree + 1); - const unsigned int dofs_per_component = - this->data->dofs_per_component_on_cell; - const unsigned int n_quadrature_points = - is_face ? this->data->n_q_points_face : this->data->n_q_points; - - const unsigned int shift = - std::max(tensor_dofs_per_component + 1, dofs_per_component) * - n_components_ * 3 + - 2 * n_quadrature_points; - const unsigned int allocated_size = - shift + n_components_ * dofs_per_component + - (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points); - scratch_data_array->resize_fast(allocated_size); - - // set the pointers to the correct position in the data array - for (unsigned int c = 0; c < n_components_; ++c) - { - this->values_dofs[c] = - scratch_data_array->begin() + c * dofs_per_component; - } - this->values_quad = - scratch_data_array->begin() + n_components * dofs_per_component; - this->gradients_quad = - scratch_data_array->begin() + - n_components * (dofs_per_component + n_quadrature_points); - this->hessians_quad = - scratch_data_array->begin() + - n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); - scratch_data = - scratch_data_array->begin() + n_components_ * dofs_per_component + - (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points); -} - - - -template +template inline unsigned int -FEEvaluationBase:: +FEEvaluationBaseData:: get_mapping_data_index_offset() const { if (matrix_info == nullptr) @@ -3567,14 +3462,10 @@ FEEvaluationBase:: -template +template inline internal::MatrixFreeFunctions::GeometryType -FEEvaluationBase:: - get_cell_type() const +FEEvaluationBaseData::get_cell_type() + const { Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized()); return cell_type; @@ -3582,13 +3473,9 @@ FEEvaluationBase:: -template +template inline const internal::MatrixFreeFunctions::ShapeInfo & -FEEvaluationBase:: +FEEvaluationBaseData:: get_shape_info() const { Assert(data != nullptr, ExcInternalError()); @@ -3597,13 +3484,9 @@ FEEvaluationBase:: -template +template inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> - FEEvaluationBase:: + FEEvaluationBaseData:: get_normal_vector(const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); @@ -3616,13 +3499,9 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> -template +template inline DEAL_II_ALWAYS_INLINE VectorizedArrayType - FEEvaluationBase::JxW( + FEEvaluationBaseData::JxW( const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); @@ -3638,13 +3517,9 @@ inline DEAL_II_ALWAYS_INLINE VectorizedArrayType -template +template inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> - FEEvaluationBase:: + FEEvaluationBaseData:: inverse_jacobian(const unsigned int q_point) const { AssertIndexRange(q_point, n_quadrature_points); @@ -3657,15 +3532,13 @@ inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> -template +template inline std::array -FEEvaluationBase:: - get_cell_ids() const +FEEvaluationBaseData::get_cell_ids() + const { + Assert(this->matrix_info != nullptr, ExcNotInitialized()); + const unsigned int n_lanes = VectorizedArrayType::size(); std::array cells; @@ -3737,7 +3610,6 @@ FEEvaluationBase:: namespace internal { template inline void process_cell_data( - const FEEvaluationBase & phi, + const FEEvaluationBaseData &phi, const MatrixFree *matrix_info, GlobalVectorType & array, VectorizedArrayType2 & out, @@ -3775,13 +3643,9 @@ namespace internal -template +template std::array -FEEvaluationBase:: +FEEvaluationBaseData:: get_cell_or_face_ids() const { const unsigned int v_len = VectorizedArrayType::size(); @@ -3807,29 +3671,178 @@ FEEvaluationBase:: this->face_no, i); - if (fn == numbers::invalid_unsigned_int) - continue; // invalid face ID: no neighbor on boundary + if (fn == numbers::invalid_unsigned_int) + continue; // invalid face ID: no neighbor on boundary + + // get cell ID on both sides of face + auto cell_m = this->matrix_info->get_face_info(fn / v_len) + .cells_interior[fn % v_len]; + auto cell_p = this->matrix_info->get_face_info(fn / v_len) + .cells_exterior[fn % v_len]; + + // compare the IDs with the given cell ID + if (cell_m == cell_this) + cells[i] = cell_p; // neighbor has the other ID + else if (cell_p == cell_this) + cells[i] = cell_m; + } + } + else + { + for (unsigned int i = 0; i < v_len; ++i) + cells[i] = cell * v_len + i; + } + + return cells; +} + + + +template +inline VectorizedArrayType +FEEvaluationBaseData::read_cell_data( + const AlignedVector &array) const +{ + VectorizedArrayType out = Number(1.); + internal::process_cell_data( + *this, this->matrix_info, array, out, [](auto &local, const auto &global) { + local = global; + }); + return out; +} + + + +template +inline void +FEEvaluationBaseData::set_cell_data( + AlignedVector &array, + const VectorizedArrayType & in) const +{ + internal::process_cell_data( + *this, this->matrix_info, array, in, [](const auto &local, auto &global) { + global = local; + }); +} + + + +template +template +inline std::array +FEEvaluationBaseData::read_cell_data( + const AlignedVector> &array) const +{ + std::array out; + internal::process_cell_data( + *this, this->matrix_info, array, out, [](auto &local, const auto &global) { + local = global; + }); + return out; +} + + + +template +template +inline void +FEEvaluationBaseData::set_cell_data( + AlignedVector> &array, + const std::array & in) const +{ + internal::process_cell_data( + *this, this->matrix_info, array, in, [](const auto &local, auto &global) { + global = local; + }); +} + + + +template +inline const std::vector & +FEEvaluationBaseData:: + get_internal_dof_numbering() const +{ + return data->lexicographic_numbering; +} + + + +template +inline ArrayView +FEEvaluationBaseData:: + get_scratch_data() const +{ + return ArrayView( + const_cast(scratch_data), + scratch_data_array->end() - scratch_data); +} + - // get cell ID on both sides of face - auto cell_m = this->matrix_info->get_face_info(fn / v_len) - .cells_interior[fn % v_len]; - auto cell_p = this->matrix_info->get_face_info(fn / v_len) - .cells_exterior[fn % v_len]; +/*----------------------- FEEvaluationBase ----------------------------------*/ - // compare the IDs with the given cell ID - if (cell_m == cell_this) - cells[i] = cell_p; // neighbor has the other ID - else if (cell_p == cell_this) - cells[i] = cell_m; - } - } - else - { - for (unsigned int i = 0; i < v_len; ++i) - cells[i] = cell * v_len + i; - } +template +inline FEEvaluationBase:: + FEEvaluationBase(const MatrixFree &data_in, + const unsigned int dof_no, + const unsigned int first_selected_component, + const unsigned int quad_no_in, + const unsigned int fe_degree, + const unsigned int n_q_points, + const bool is_interior_face) + : FEEvaluationBaseData( + data_in, + dof_no, + first_selected_component, + quad_no_in, + fe_degree, + n_q_points, + is_interior_face) + , n_fe_components(data_in.get_dof_info(dof_no).start_components.back()) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + , first_selected_component(first_selected_component) +{ + set_data_pointers(); + Assert( + this->dof_info->start_components.back() == 1 || + static_cast(n_components_) <= + static_cast( + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component] + + 1]) - + first_selected_component, + ExcMessage( + "You tried to construct a vector-valued evaluator with " + + std::to_string(n_components) + + " components. However, " + "the current base element has only " + + std::to_string( + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component] + + 1] - + first_selected_component) + + " components left when starting from local element index " + + std::to_string( + first_selected_component - + this->dof_info->start_components + [this->dof_info->component_to_base_index[first_selected_component]]) + + " (global index " + std::to_string(first_selected_component) + ")")); - return cells; + // do not check for correct dimensions of data fields here, should be done + // in derived classes } @@ -3839,16 +3852,48 @@ template -inline VectorizedArrayType -FEEvaluationBase:: - read_cell_data(const AlignedVector &array) const +inline FEEvaluationBase:: + FEEvaluationBase( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other) + : FEEvaluationBaseData( + mapping, + fe, + quadrature, + update_flags, + first_selected_component, + other) + , n_fe_components(n_components_) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + // keep the number of the selected component within the current base element + // for reading dof values + , first_selected_component(first_selected_component) { - VectorizedArrayType out = Number(1.); - internal::process_cell_data( - *this, this->matrix_info, array, out, [](auto &local, const auto &global) { - local = global; - }); - return out; + set_data_pointers(); + + const unsigned int base_element_number = + fe.component_to_base_index(first_selected_component).first; + Assert(fe.element_multiplicity(base_element_number) == 1 || + fe.element_multiplicity(base_element_number) - + first_selected_component >= + n_components_, + ExcMessage("The underlying element must at least contain as many " + "components as requested by this class")); + (void)base_element_number; } @@ -3858,15 +3903,27 @@ template -inline void -FEEvaluationBase:: - set_cell_data(AlignedVector &array, - const VectorizedArrayType & in) const +inline FEEvaluationBase:: + FEEvaluationBase(const FEEvaluationBase &other) + : FEEvaluationBaseData(other) + , n_fe_components(other.n_fe_components) + , dof_values_initialized(false) + , values_quad_initialized(false) + , gradients_quad_initialized(false) + , hessians_quad_initialized(false) + , values_quad_submitted(false) + , gradients_quad_submitted(false) + , first_selected_component(other.first_selected_component) { - internal::process_cell_data( - *this, this->matrix_info, array, in, [](const auto &local, auto &global) { - global = local; - }); + set_data_pointers(); } @@ -3876,18 +3933,24 @@ template -template -inline std::array +inline FEEvaluationBase & FEEvaluationBase:: - read_cell_data(const AlignedVector> - &array) const +operator=(const FEEvaluationBase &other) { - std::array out; - internal::process_cell_data( - *this, this->matrix_info, array, out, [](auto &local, const auto &global) { - local = global; - }); - return out; + this->FEEvaluationBaseData:: + operator=(other); + AssertDimension(n_fe_components, other.n_fe_components); + AssertDimension(first_selected_component, other.first_selected_component); + + return *this; } @@ -3897,17 +3960,44 @@ template -template inline void FEEvaluationBase:: - set_cell_data( - AlignedVector> &array, - const std::array & in) const + set_data_pointers() { - internal::process_cell_data( - *this, this->matrix_info, array, in, [](const auto &local, auto &global) { - global = local; - }); + Assert(this->scratch_data_array != nullptr, ExcInternalError()); + + const unsigned int tensor_dofs_per_component = + Utilities::fixed_power(this->data->data.front().fe_degree + 1); + const unsigned int dofs_per_component = + this->data->dofs_per_component_on_cell; + const unsigned int n_quadrature_points = + is_face ? this->data->n_q_points_face : this->data->n_q_points; + + const unsigned int shift = + std::max(tensor_dofs_per_component + 1, dofs_per_component) * + n_components_ * 3 + + 2 * n_quadrature_points; + const unsigned int allocated_size = + shift + n_components_ * dofs_per_component + + (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points); + this->scratch_data_array->resize_fast(allocated_size); + + // set the pointers to the correct position in the data array + for (unsigned int c = 0; c < n_components_; ++c) + { + values_dofs[c] = + this->scratch_data_array->begin() + c * dofs_per_component; + } + values_quad = + this->scratch_data_array->begin() + n_components * dofs_per_component; + gradients_quad = this->scratch_data_array->begin() + + n_components * (dofs_per_component + n_quadrature_points); + hessians_quad = + this->scratch_data_array->begin() + + n_components * (dofs_per_component + (dim + 1) * n_quadrature_points); + this->scratch_data = + this->scratch_data_array->begin() + n_components_ * dofs_per_component + + (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points); } @@ -4006,14 +4096,14 @@ FEEvaluationBase:: // Case 1: No MatrixFree object given, simple case because we do not need to // process constraints and need not care about vectorization -> go to // separate function - if (matrix_info == nullptr) + if (this->matrix_info == nullptr) { read_write_operation_global(operation, src); return; } - Assert(dof_info != nullptr, ExcNotInitialized()); - Assert(matrix_info->indices_initialized() == true, ExcNotInitialized()); + Assert(this->dof_info != nullptr, ExcNotInitialized()); + Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized()); if (n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) { @@ -4025,21 +4115,22 @@ FEEvaluationBase:: "FEEvaluation. In that case, you must pass an " "std::vector or a BlockVector to " + "read_dof_values and distribute_local_to_global.")); - internal::check_vector_compatibility(*src[comp], *dof_info); + internal::check_vector_compatibility(*src[comp], *this->dof_info); } else { - internal::check_vector_compatibility(*src[0], *dof_info); + internal::check_vector_compatibility(*src[0], *this->dof_info); } // Case 2: contiguous indices which use reduced storage of indices and can // use vectorized load/store operations -> go to separate function - AssertIndexRange(cell, - dof_info->index_storage_variants[dof_access_index].size()); - if (dof_info->index_storage_variants - [is_face ? dof_access_index : + AssertIndexRange( + this->cell, + this->dof_info->index_storage_variants[this->dof_access_index].size()); + if (this->dof_info->index_storage_variants + [is_face ? this->dof_access_index : internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [cell] >= + [this->cell] >= internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous) { read_write_operation_contiguous(operation, src, mask); @@ -4059,17 +4150,19 @@ FEEvaluationBase:: const unsigned int dofs_per_component = this->data->dofs_per_component_on_cell; - if (dof_info->index_storage_variants - [is_face ? dof_access_index : + if (this->dof_info->index_storage_variants + [is_face ? this->dof_access_index : internal::MatrixFreeFunctions::DoFInfo::dof_access_cell] - [cell] == + [this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved) { const unsigned int *dof_indices = - dof_info->dof_indices_interleaved.data() + - dof_info->row_starts[cell * n_fe_components * n_lanes].first + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] * + this->dof_info->dof_indices_interleaved.data() + + this->dof_info->row_starts[this->cell * n_fe_components * n_lanes] + .first + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * n_lanes; if (n_components == 1 || n_fe_components == 1) for (unsigned int i = 0; i < dofs_per_component; @@ -4098,27 +4191,29 @@ FEEvaluationBase:: unsigned int cells_copied[n_lanes]; const unsigned int *cells; unsigned int n_vectorization_actual = - dof_info->n_vectorization_lanes_filled[dof_access_index][cell]; + this->dof_info + ->n_vectorization_lanes_filled[this->dof_access_index][this->cell]; bool has_constraints = false; if (is_face) { - if (dof_access_index == + if (this->dof_access_index == internal::MatrixFreeFunctions::DoFInfo::dof_access_cell) for (unsigned int v = 0; v < n_vectorization_actual; ++v) - cells_copied[v] = cell * VectorizedArrayType::size() + v; - cells = dof_access_index == - internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ? - &cells_copied[0] : - (is_interior_face ? - &this->matrix_info->get_face_info(cell).cells_interior[0] : - &this->matrix_info->get_face_info(cell).cells_exterior[0]); + cells_copied[v] = this->cell * VectorizedArrayType::size() + v; + cells = + this->dof_access_index == + internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ? + &cells_copied[0] : + (this->is_interior_face ? + &this->matrix_info->get_face_info(this->cell).cells_interior[0] : + &this->matrix_info->get_face_info(this->cell).cells_exterior[0]); for (unsigned int v = 0; v < n_vectorization_actual; ++v) { - Assert(cells[v] < dof_info->row_starts.size() - 1, + Assert(cells[v] < this->dof_info->row_starts.size() - 1, ExcInternalError()); const std::pair *my_index_start = - &dof_info->row_starts[cells[v] * n_fe_components + - first_selected_component]; + &this->dof_info->row_starts[cells[v] * n_fe_components + + first_selected_component]; // check whether any of the SIMD lanes has constraints, i.e., the // constraint indicator which is the second entry of row_starts @@ -4127,33 +4222,34 @@ FEEvaluationBase:: has_constraints = true; dof_indices[v] = - dof_info->dof_indices.data() + my_index_start[0].first; + this->dof_info->dof_indices.data() + my_index_start[0].first; } for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v) dof_indices[v] = nullptr; } else { - AssertIndexRange((cell + 1) * n_lanes * n_fe_components, - dof_info->row_starts.size()); + AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components, + this->dof_info->row_starts.size()); const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1; for (unsigned int v = 0; v < n_vectorization_actual; ++v) { const std::pair *my_index_start = - &dof_info->row_starts[(cell * n_lanes + v) * n_fe_components + - first_selected_component]; + &this->dof_info + ->row_starts[(this->cell * n_lanes + v) * n_fe_components + + first_selected_component]; if (my_index_start[n_components_read].second != my_index_start[0].second) has_constraints = true; Assert(my_index_start[n_components_read].first == my_index_start[0].first || - my_index_start[0].first < dof_info->dof_indices.size(), + my_index_start[0].first < this->dof_info->dof_indices.size(), ExcIndexRange(0, my_index_start[0].first, - dof_info->dof_indices.size())); + this->dof_info->dof_indices.size())); dof_indices[v] = - dof_info->dof_indices.data() + my_index_start[0].first; + this->dof_info->dof_indices.data() + my_index_start[0].first; } for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v) dof_indices[v] = nullptr; @@ -4200,30 +4296,33 @@ FEEvaluationBase:: operation.process_empty(values_dofs[comp][i]); for (unsigned int v = 0; v < n_vectorization_actual; ++v) { - const unsigned int cell_index = is_face ? cells[v] : cell * n_lanes + v; + const unsigned int cell_index = + is_face ? cells[v] : this->cell * n_lanes + v; const unsigned int cell_dof_index = cell_index * n_fe_components + first_selected_component; const unsigned int n_components_read = n_fe_components > 1 ? n_components : 1; unsigned int index_indicators = - dof_info->row_starts[cell_dof_index].second; + this->dof_info->row_starts[cell_dof_index].second; unsigned int next_index_indicators = - dof_info->row_starts[cell_dof_index + 1].second; + this->dof_info->row_starts[cell_dof_index + 1].second; // For read_dof_values_plain, redirect the dof_indices field to the // unconstrained indices if (apply_constraints == false && - dof_info->row_starts[cell_dof_index].second != - dof_info->row_starts[cell_dof_index + n_components_read].second) + this->dof_info->row_starts[cell_dof_index].second != + this->dof_info->row_starts[cell_dof_index + n_components_read] + .second) { - Assert(dof_info->row_starts_plain_indices[cell_index] != + Assert(this->dof_info->row_starts_plain_indices[cell_index] != numbers::invalid_unsigned_int, ExcNotInitialized()); dof_indices[v] = - dof_info->plain_dof_indices.data() + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] + - dof_info->row_starts_plain_indices[cell_index]; + this->dof_info->plain_dof_indices.data() + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] + + this->dof_info->row_starts_plain_indices[cell_index]; next_index_indicators = index_indicators; } @@ -4233,7 +4332,7 @@ FEEvaluationBase:: for (; index_indicators != next_index_indicators; ++index_indicators) { const std::pair indicator = - dof_info->constraint_indicator[index_indicators]; + this->dof_info->constraint_indicator[index_indicators]; // run through values up to next constraint for (unsigned int j = 0; j < indicator.first; ++j) for (unsigned int comp = 0; comp < n_components; ++comp) @@ -4252,9 +4351,9 @@ FEEvaluationBase:: value[comp]); const Number *data_val = - matrix_info->constraint_pool_begin(indicator.second); + this->matrix_info->constraint_pool_begin(indicator.second); const Number *end_pool = - matrix_info->constraint_pool_end(indicator.second); + this->matrix_info->constraint_pool_end(indicator.second); for (; data_val != end_pool; ++data_val, ++dof_indices[v]) for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_constraint(*dof_indices[v], @@ -4292,7 +4391,7 @@ FEEvaluationBase:: ++index_indicators) { const std::pair indicator = - dof_info->constraint_indicator[index_indicators]; + this->dof_info->constraint_indicator[index_indicators]; // run through values up to next constraint for (unsigned int j = 0; j < indicator.first; ++j) @@ -4309,9 +4408,9 @@ FEEvaluationBase:: value); const Number *data_val = - matrix_info->constraint_pool_begin(indicator.second); + this->matrix_info->constraint_pool_begin(indicator.second); const Number *end_pool = - matrix_info->constraint_pool_end(indicator.second); + this->matrix_info->constraint_pool_end(indicator.second); for (; data_val != end_pool; ++data_val, ++dof_indices[v]) operation.process_constraint(*dof_indices[v], @@ -4338,7 +4437,7 @@ FEEvaluationBase:: if (apply_constraints == true && comp + 1 < n_components) next_index_indicators = - dof_info->row_starts[cell_dof_index + comp + 2].second; + this->dof_info->row_starts[cell_dof_index + comp + 2].second; } } } @@ -4360,15 +4459,15 @@ FEEvaluationBase:: Assert(!local_dof_indices.empty(), ExcNotInitialized()); unsigned int index = - first_selected_component * data->dofs_per_component_on_cell; + first_selected_component * this->data->dofs_per_component_on_cell; for (unsigned int comp = 0; comp < n_components; ++comp) { - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell; ++i, ++index) { operation.process_empty(values_dofs[comp][i]); operation.process_dof_global( - local_dof_indices[data->lexicographic_numbering[index]], + local_dof_indices[this->data->lexicographic_numbering[index]], *src[0], values_dofs[comp][i][0]); } @@ -4401,55 +4500,57 @@ FEEvaluationBase:: internal::is_vectorizable::value> vector_selector; const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind = - is_face ? dof_access_index : + is_face ? this->dof_access_index : internal::MatrixFreeFunctions::DoFInfo::dof_access_cell; const unsigned int n_lanes = mask.count(); const std::vector &dof_indices_cont = - dof_info->dof_indices_contiguous[ind]; + this->dof_info->dof_indices_contiguous[ind]; // Simple case: We have contiguous storage, so we can simply copy out the // data - if (dof_info->index_storage_variants[ind][cell] == + if (this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous && n_lanes == VectorizedArrayType::size()) { const unsigned int dof_index = - dof_indices_cont[cell * VectorizedArrayType::size()] + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] * + dof_indices_cont[this->cell * VectorizedArrayType::size()] + + this->dof_info->component_dof_indices_offset[this->active_fe_index] + [first_selected_component] * VectorizedArrayType::size(); if (n_components == 1 || n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) - operation.process_dofs_vectorized(data->dofs_per_component_on_cell, - dof_index, - *src[comp], - values_dofs[comp], - vector_selector); + operation.process_dofs_vectorized( + this->data->dofs_per_component_on_cell, + dof_index, + *src[comp], + values_dofs[comp], + vector_selector); else - operation.process_dofs_vectorized(data->dofs_per_component_on_cell * - n_components, - dof_index, - *src[0], - values_dofs[0], - vector_selector); + operation.process_dofs_vectorized( + this->data->dofs_per_component_on_cell * n_components, + dof_index, + *src[0], + values_dofs[0], + vector_selector); return; } // More general case: Must go through the components one by one and apply // some transformations const unsigned int n_filled_lanes = - dof_info->n_vectorization_lanes_filled[ind][this->cell]; + this->dof_info->n_vectorization_lanes_filled[ind][this->cell]; unsigned int dof_indices[VectorizedArrayType::size()]; for (unsigned int v = 0; v < n_filled_lanes; ++v) dof_indices[v] = - dof_indices_cont[cell * VectorizedArrayType::size() + v] + - dof_info->component_dof_indices_offset[active_fe_index] - [first_selected_component] * - dof_info->dof_indices_interleave_strides - [ind][cell * VectorizedArrayType::size() + v]; + dof_indices_cont[this->cell * VectorizedArrayType::size() + v] + + this->dof_info + ->component_dof_indices_offset[this->active_fe_index] + [this->first_selected_component] * + this->dof_info->dof_indices_interleave_strides + [ind][this->cell * VectorizedArrayType::size() + v]; for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v) dof_indices[v] = numbers::invalid_unsigned_int; @@ -4459,32 +4560,33 @@ FEEvaluationBase:: if (n_filled_lanes == VectorizedArrayType::size() && n_lanes == VectorizedArrayType::size()) { - if (dof_info->index_storage_variants[ind][cell] == + if (this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: contiguous) { if (n_components == 1 || n_fe_components == 1) for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_dofs_vectorized_transpose( - data->dofs_per_component_on_cell, + this->data->dofs_per_component_on_cell, dof_indices, *src[comp], values_dofs[comp], vector_selector); else operation.process_dofs_vectorized_transpose( - data->dofs_per_component_on_cell * n_components, + this->data->dofs_per_component_on_cell * n_components, dof_indices, *src[0], &values_dofs[0][0], vector_selector); } - else if (dof_info->index_storage_variants[ind][cell] == + else if (this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: interleaved_contiguous_strided) { if (n_components == 1 || n_fe_components == 1) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) + for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell; + ++i) { for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_dof_gather(dof_indices, @@ -4495,13 +4597,14 @@ FEEvaluationBase:: } else for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + for (unsigned int i = 0; + i < this->data->dofs_per_component_on_cell; ++i) { operation.process_dof_gather( dof_indices, *src[0], - (comp * data->dofs_per_component_on_cell + i) * + (comp * this->data->dofs_per_component_on_cell + i) * VectorizedArrayType::size(), values_dofs[comp][i], vector_selector); @@ -4509,15 +4612,16 @@ FEEvaluationBase:: } else { - Assert(dof_info->index_storage_variants[ind][cell] == + Assert(this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo:: IndexStorageVariants::interleaved_contiguous_mixed_strides, ExcNotImplemented()); const unsigned int *offsets = - &dof_info->dof_indices_interleave_strides - [ind][VectorizedArrayType::size() * cell]; + &this->dof_info->dof_indices_interleave_strides + [ind][VectorizedArrayType::size() * this->cell]; if (n_components == 1 || n_fe_components == 1) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) + for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell; + ++i) { for (unsigned int comp = 0; comp < n_components; ++comp) operation.process_dof_gather(dof_indices, @@ -4531,7 +4635,8 @@ FEEvaluationBase:: } else for (unsigned int comp = 0; comp < n_components; ++comp) - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; + for (unsigned int i = 0; + i < this->data->dofs_per_component_on_cell; ++i) { operation.process_dof_gather(dof_indices, @@ -4548,9 +4653,10 @@ FEEvaluationBase:: else for (unsigned int comp = 0; comp < n_components; ++comp) { - for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i) + for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell; + ++i) operation.process_empty(values_dofs[comp][i]); - if (dof_info->index_storage_variants[ind][cell] == + if (this->dof_info->index_storage_variants[ind][this->cell] == internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants:: contiguous) { @@ -4559,7 +4665,7 @@ FEEvaluationBase:: for (unsigned int v = 0; v < n_filled_lanes; ++v) if (mask[v] == true) for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; + i < this->data->dofs_per_component_on_cell; ++i) operation.process_dof(dof_indices[v] + i, *src[comp], @@ -4570,11 +4676,11 @@ FEEvaluationBase:: for (unsigned int v = 0; v < n_filled_lanes; ++v) if (mask[v] == true) for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; + i < this->data->dofs_per_component_on_cell; ++i) operation.process_dof( dof_indices[v] + i + - comp * data->dofs_per_component_on_cell, + comp * this->data->dofs_per_component_on_cell, *src[0], values_dofs[comp][i][v]); } @@ -4582,8 +4688,8 @@ FEEvaluationBase:: else { const unsigned int *offsets = - &dof_info->dof_indices_interleave_strides - [ind][VectorizedArrayType::size() * cell]; + &this->dof_info->dof_indices_interleave_strides + [ind][VectorizedArrayType::size() * this->cell]; for (unsigned int v = 0; v < n_filled_lanes; ++v) AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1); if (n_components == 1 || n_fe_components == 1) @@ -4591,7 +4697,7 @@ FEEvaluationBase:: { if (mask[v] == true) for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; + i < this->data->dofs_per_component_on_cell; ++i) operation.process_dof(dof_indices[v] + i * offsets[v], *src[comp], @@ -4602,11 +4708,11 @@ FEEvaluationBase:: for (unsigned int v = 0; v < n_filled_lanes; ++v) if (mask[v] == true) for (unsigned int i = 0; - i < data->dofs_per_component_on_cell; + i < this->data->dofs_per_component_on_cell; ++i) operation.process_dof( dof_indices[v] + - (i + comp * data->dofs_per_component_on_cell) * + (i + comp * this->data->dofs_per_component_on_cell) * offsets[v], *src[0], values_dofs[comp][i][v]); @@ -4757,34 +4863,6 @@ FEEvaluationBase:: /*------------------------------ access to data fields ----------------------*/ -template -inline const std::vector & -FEEvaluationBase:: - get_internal_dof_numbering() const -{ - return data->lexicographic_numbering; -} - - - -template -inline ArrayView -FEEvaluationBase:: - get_scratch_data() const -{ - return ArrayView( - const_cast(scratch_data), - scratch_data_array->end() - scratch_data); -} - template n_quadrature_points); - Assert(jacobian != nullptr, ExcNotInitialized()); + Assert(this->jacobian != nullptr, ExcNotInitialized()); const std::size_t nqp = this->n_quadrature_points; Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out; @@ -4991,15 +5069,15 @@ inline DEAL_II_ALWAYS_INLINE for (unsigned int d = 0; d < dim; ++d) for (unsigned int comp = 0; comp < n_components; comp++) grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] * - jacobian[0][d][d]; + this->jacobian[0][d][d]; } // cell with general/affine Jacobian else { const Tensor<2, dim, VectorizedArrayType> &jac = - jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ? - q_point : - 0]; + this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ? + q_point : + 0]; for (unsigned int comp = 0; comp < n_components; comp++) for (unsigned int d = 0; d < dim; ++d) { @@ -5030,7 +5108,7 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType> internal::ExcAccessToUninitializedField()); # endif - Assert(normal_x_jacobian != nullptr, ExcNotInitialized()); + Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized()); const std::size_t nqp = this->n_quadrature_points; Tensor<1, n_components, VectorizedArrayType> grad_out; @@ -5126,11 +5204,11 @@ FEEvaluationBase:: # endif AssertIndexRange(q_point, this->n_quadrature_points); - Assert(jacobian != nullptr, ExcNotImplemented()); + Assert(this->jacobian != nullptr, ExcNotImplemented()); const Tensor<2, dim, VectorizedArrayType> &jac = - jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? - 0 : - q_point]; + this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? + 0 : + q_point]; Tensor<1, n_components, Tensor<2, dim, VectorizedArrayType>> hessian_out; @@ -5205,7 +5283,7 @@ FEEvaluationBase:: else { const auto &jac_grad = - mapping_data->jacobian_gradients + this->mapping_data->jacobian_gradients [1 - this->is_interior_face] [this->mapping_data->data_index_offsets[this->cell] + q_point]; for (unsigned int comp = 0; comp < n_components; comp++) @@ -5267,11 +5345,11 @@ FEEvaluationBase:: # endif AssertIndexRange(q_point, this->n_quadrature_points); - Assert(jacobian != nullptr, ExcNotImplemented()); + Assert(this->jacobian != nullptr, ExcNotImplemented()); const Tensor<2, dim, VectorizedArrayType> &jac = - jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? - 0 : - q_point]; + this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? + 0 : + q_point]; const std::size_t nqp = this->n_quadrature_points; constexpr unsigned int hdim = (dim * (dim + 1)) / 2; @@ -5312,7 +5390,7 @@ FEEvaluationBase:: { const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>> &jac_grad = - mapping_data->jacobian_gradients + this->mapping_data->jacobian_gradients [0][this->mapping_data->data_index_offsets[this->cell] + q_point]; for (unsigned int comp = 0; comp < n_components; comp++) { @@ -5412,13 +5490,14 @@ FEEvaluationBase:: const std::size_t nqp = this->n_quadrature_points; if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point]; + const VectorizedArrayType JxW = + this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int comp = 0; comp < n_components; ++comp) values_quad[comp * nqp + q_point] = val_in[comp] * JxW; } else { - const VectorizedArrayType JxW = J_value[q_point]; + const VectorizedArrayType JxW = this->J_value[q_point]; for (unsigned int comp = 0; comp < n_components; ++comp) values_quad[comp * nqp + q_point] = val_in[comp] * JxW; } @@ -5448,10 +5527,11 @@ FEEvaluationBase:: const std::size_t nqp = this->n_quadrature_points; if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) { - const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point]; + const VectorizedArrayType JxW = + this->J_value[0] * this->quadrature_weights[q_point]; for (unsigned int d = 0; d < dim; ++d) { - const VectorizedArrayType factor = jacobian[0][d][d] * JxW; + const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW; for (unsigned int comp = 0; comp < n_components; comp++) gradients_quad[(comp * dim + d) * nqp + q_point] = grad_in[comp][d] * factor; @@ -5461,12 +5541,12 @@ FEEvaluationBase:: { const Tensor<2, dim, VectorizedArrayType> jac = this->cell_type > internal::MatrixFreeFunctions::affine ? - jacobian[q_point] : - jacobian[0]; + this->jacobian[q_point] : + this->jacobian[0]; const VectorizedArrayType JxW = this->cell_type > internal::MatrixFreeFunctions::affine ? - J_value[q_point] : - J_value[0] * quadrature_weights[q_point]; + 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) { @@ -5591,22 +5671,19 @@ template -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 FEEvaluationBase *other) + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other) : FEEvaluationBase( mapping, fe, @@ -5692,18 +5769,15 @@ inline FEEvaluationAccess:: template -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 FEEvaluationBase *other) + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other) : FEEvaluationBase( mapping, fe, @@ -5994,18 +6068,15 @@ inline FEEvaluationAccess:: template -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 FEEvaluationBase *other) + FEEvaluationAccess( + const Mapping & mapping, + const FiniteElement &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData + *other) : FEEvaluationBase( mapping, fe, @@ -6392,18 +6463,14 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: template -template inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>:: - FEEvaluationAccess(const Mapping<1> & mapping, - const FiniteElement<1> &fe, - const Quadrature<1> & quadrature, - const UpdateFlags update_flags, - const unsigned int first_selected_component, - const FEEvaluationBase<1, - n_components_other, - Number, - is_face, - VectorizedArrayType> *other) + FEEvaluationAccess( + const Mapping<1> & mapping, + const FiniteElement<1> &fe, + const Quadrature<1> & quadrature, + const UpdateFlags update_flags, + const unsigned int first_selected_component, + const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other) : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>( mapping, fe, @@ -6708,9 +6775,7 @@ inline FEEvaluation *>( - nullptr)) + nullptr) , dofs_per_component(this->data->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) @@ -6741,9 +6806,7 @@ inline FEEvaluation *>( - nullptr)) + nullptr) , dofs_per_component(this->data->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) @@ -6759,20 +6822,16 @@ template -template inline FEEvaluation:: - FEEvaluation(const FiniteElement & fe, - const FEEvaluationBase &other, - const unsigned int first_selected_component) + FEEvaluation( + const FiniteElement & fe, + const FEEvaluationBaseData &other, + const unsigned int first_selected_component) : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(), fe, other.mapped_geometry->get_quadrature(), -- 2.39.5