/**
- * 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<Number>
*
* @ingroup matrixfree
*/
template <int dim,
- int n_components_,
typename Number,
bool is_face = false,
typename VectorizedArrayType = VectorizedArray<Number>>
-class FEEvaluationBase
+class FEEvaluationBaseData
{
static_assert(
std::is_same<Number, typename VectorizedArrayType::value_type>::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
const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
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<VectorizedArrayType> &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<VectorizedArrayType> &array,
+ const VectorizedArrayType & value) const;
+
+ /**
+ * The same as above, just for std::array of length of VectorizedArrayType for
+ * arbitrary data type.
+ */
+ template <typename T>
+ std::array<T, VectorizedArrayType::size()>
+ read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
+ &array) const;
+
+ /**
+ * The same as above, just for std::array of length of VectorizedArrayType for
+ * arbitrary data type.
+ */
+ template <typename T>
+ void
+ set_cell_data(
+ AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+ const std::array<T, VectorizedArrayType::size()> & value) const;
+
+ /**
+ * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
+ * associated with.
+ */
+ std::array<unsigned int, VectorizedArrayType::size()>
+ get_cell_ids() const;
+
+ /**
+ * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
+ * associated with.
+ */
+ std::array<unsigned int, VectorizedArrayType::size()>
+ 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<unsigned int> &
+ 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<VectorizedArrayType>
+ 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<dim, Number, VectorizedArrayType> &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<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *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<VectorizedArrayType> *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<dim, Number, VectorizedArrayType> *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<VectorizedArrayType> *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<internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
+ mapped_geometry;
+
+ // Make FEEvaluation objects friends for access to protected member
+ // mapped_geometry.
+ template <int, int, int, int, typename, typename>
+ 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<Number>
+ *
+ * @note Currently only VectorizedArray<Number, width> is supported as
+ * VectorizedArrayType.
+ *
+ *
+ * @ingroup matrixfree
+ */
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face = false,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+class FEEvaluationBase
+ : public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+{
+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
*/
//@{
/**
//@}
/**
- * @name 3: Data access
+ * @name 2: Access to data at quadrature points or the gather vector data
*/
//@{
/**
#ifdef DOXYGEN
// doxygen does not anyhow mention functions coming from partial template
// specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
- // 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
#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
*
* @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<VectorizedArrayType> &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<VectorizedArrayType> &array,
- const VectorizedArrayType & value) const;
-
- /**
- * The same as above, just for std::array of length of VectorizedArrayType for
- * arbitrary data type.
- */
- template <typename T>
- std::array<T, VectorizedArrayType::size()>
- read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
- &array) const;
-
- /**
- * The same as above, just for std::array of length of VectorizedArrayType for
- * arbitrary data type.
- */
- template <typename T>
- void
- set_cell_data(
- AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
- const std::array<T, VectorizedArrayType::size()> & value) const;
-
- /**
- * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
- * associated with.
- */
- std::array<unsigned int, VectorizedArrayType::size()>
- 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<unsigned int, VectorizedArrayType::size()>
- get_cell_or_face_ids() const;
+ value_type
+ integrate_value() const;
//@}
/**
- * @name 4: Access to internal data
+ * @name 3: Access to internal data
*/
//@{
/**
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<unsigned int> &
- 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<VectorizedArrayType>
- get_scratch_data() const;
-
//@}
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 <int n_components_other>
- FEEvaluationBase(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *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<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other);
/**
* Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
read_write_operation_global(const VectorOperation &operation,
VectorType * vectors[]) const;
- /**
- * This is the general array for all data fields.
- */
- AlignedVector<VectorizedArrayType> *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
*/
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<dim, Number, VectorizedArrayType> *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<VectorizedArrayType> *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
*/
bool gradients_quad_submitted;
- /**
- * Geometry data that can be generated FEValues on the fly with the
- * respective constructor.
- */
- std::shared_ptr<internal::MatrixFreeFunctions::
- MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
- mapped_geometry;
-
/**
* For a FiniteElement with more than one base element, select at which
* component this data structure should start.
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 <int, int, typename, bool, typename>
- friend class FEEvaluationBase;
- template <int, int, int, int, typename, typename>
- friend class FEEvaluation;
};
* Constructor with reduced functionality for similar usage of FEEvaluation
* as FEValues, including matrix assembly.
*/
- template <int n_components_other>
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other);
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other);
/**
* Copy constructor
* Constructor with reduced functionality for similar usage of FEEvaluation
* as FEValues, including matrix assembly.
*/
- template <int n_components_other>
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other);
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other);
/**
* Copy constructor
* Constructor with reduced functionality for similar usage of FEEvaluation
* as FEValues, including matrix assembly.
*/
- template <int n_components_other>
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other);
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other);
/**
* Copy constructor
* Constructor with reduced functionality for similar usage of FEEvaluation
* as FEValues, including matrix assembly.
*/
- template <int n_components_other>
- FEEvaluationAccess(const Mapping<1> & mapping,
- const FiniteElement<1> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<1,
- n_components_other,
- Number,
- 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
* use the FEEvaluation object in %parallel to the given one because
* otherwise the intended sharing may create race conditions.
*/
- template <int n_components_other>
- FEEvaluation(const FiniteElement<dim> & fe,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- false,
- VectorizedArrayType> &other,
- const unsigned int first_selected_component = 0);
+ FEEvaluation(
+ const FiniteElement<dim> & fe,
+ const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &other,
+ const unsigned int first_selected_component = 0);
/**
* Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
#ifndef DOXYGEN
+/*----------------------- FEEvaluationBaseData ------------------------*/
-/*----------------------- FEEvaluationBase ----------------------------------*/
-
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType>::
- FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &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 <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(
+ const MatrixFree<dim, Number, VectorizedArrayType> &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,
internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
, cell_type(internal::MatrixFreeFunctions::general)
- , 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(matrix_info->mapping_initialized() == true, ExcNotInitialized());
AssertDimension(matrix_info->get_task_info().vectorization_length,
VectorizedArrayType::size());
n_quadrature_points);
AssertDimension(n_quadrature_points,
mapping_data->descriptor[active_quad_index].n_q_points);
- Assert(
- dof_info->start_components.back() == 1 ||
- static_cast<int>(n_components_) <=
- static_cast<int>(
- 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 <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-template <int n_components_other>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType>::
- FEEvaluationBase(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other)
: scratch_data_array(new AlignedVector<VectorizedArrayType>())
, 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(
, J_value(nullptr)
, normal_vectors(nullptr)
, normal_x_jacobian(nullptr)
- , quadrature_weights(nullptr)
- , cell(0)
- , cell_type(internal::MatrixFreeFunctions::general)
- , is_interior_face(true)
- , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
- , 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)
+ , quadrature_weights(nullptr)
+ , cell(0)
+ , cell_type(internal::MatrixFreeFunctions::general)
+ , is_interior_face(true)
+ , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
{
- set_data_pointers();
-
Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
ExcInternalError());
if (other != nullptr &&
mapping_data = &mapped_geometry->get_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 <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType>::
- FEEvaluationBase(const FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> &other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ &other)
: scratch_data_array(other.matrix_info == nullptr ?
new AlignedVector<VectorizedArrayType>() :
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)
, cell_type(internal::MatrixFreeFunctions::general)
, is_interior_face(other.is_interior_face)
, dof_access_index(other.dof_access_index)
- , 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)
{
- set_data_pointers();
-
// Create deep copy of mapped geometry for use in parallel...
if (other.mapped_geometry.get() != nullptr)
{
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> &
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-operator=(const FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> &other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::operator=(
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &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)
data = other.data;
scratch_data_array = matrix_info->acquire_scratch_data();
}
- set_data_pointers();
quadrature_weights =
(mapping_data != nullptr ?
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType>::~FEEvaluationBase()
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ ~FEEvaluationBaseData()
{
if (matrix_info != nullptr)
{
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- set_data_pointers()
-{
- Assert(scratch_data_array != nullptr, ExcInternalError());
-
- const unsigned int tensor_dofs_per_component =
- Utilities::fixed_power<dim>(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 <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline unsigned int
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
get_mapping_data_index_offset() const
{
if (matrix_info == nullptr)
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline internal::MatrixFreeFunctions::GeometryType
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- get_cell_type() const
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
+ const
{
Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
return cell_type;
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
get_shape_info() const
{
Assert(data != nullptr, ExcInternalError());
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
- FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
get_normal_vector(const unsigned int q_point) const
{
AssertIndexRange(q_point, n_quadrature_points);
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
- FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::JxW(
+ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
const unsigned int q_point) const
{
AssertIndexRange(q_point, n_quadrature_points);
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
- FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
inverse_jacobian(const unsigned int q_point) const
{
AssertIndexRange(q_point, n_quadrature_points);
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
inline std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- get_cell_ids() const
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
+ const
{
+ Assert(this->matrix_info != nullptr, ExcNotInitialized());
+
const unsigned int n_lanes = VectorizedArrayType::size();
std::array<unsigned int, n_lanes> cells;
namespace internal
{
template <int dim,
- int n_components_,
typename Number,
bool is_face,
typename VectorizedArrayType,
typename FU>
inline void
process_cell_data(
- const FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> & phi,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &phi,
const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info,
GlobalVectorType & array,
VectorizedArrayType2 & out,
-template <int dim,
- int n_components_,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
get_cell_or_face_ids() const
{
const unsigned int v_len = VectorizedArrayType::size();
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 <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline VectorizedArrayType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::read_cell_data(
+ const AlignedVector<VectorizedArrayType> &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 <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::set_cell_data(
+ AlignedVector<VectorizedArrayType> &array,
+ const VectorizedArrayType & in) const
+{
+ internal::process_cell_data(
+ *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+ global = local;
+ });
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <typename T>
+inline std::array<T, VectorizedArrayType::size()>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::read_cell_data(
+ const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
+{
+ std::array<T, VectorizedArrayType::size()> out;
+ internal::process_cell_data(
+ *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
+ local = global;
+ });
+ return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <typename T>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::set_cell_data(
+ AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+ const std::array<T, VectorizedArrayType::size()> & in) const
+{
+ internal::process_cell_data(
+ *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+ global = local;
+ });
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const std::vector<unsigned int> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_internal_dof_numbering() const
+{
+ return data->lexicographic_numbering;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline ArrayView<VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_scratch_data() const
+{
+ return ArrayView<VectorizedArrayType>(
+ const_cast<VectorizedArrayType *>(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 <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::
+ FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &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<dim, Number, is_face, VectorizedArrayType>(
+ 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<int>(n_components_) <=
+ static_cast<int>(
+ 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
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-inline VectorizedArrayType
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::
+ FEEvaluationBase(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other)
+ : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+ 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;
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- set_cell_data(AlignedVector<VectorizedArrayType> &array,
- const VectorizedArrayType & in) const
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::
+ FEEvaluationBase(const FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType> &other)
+ : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(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();
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-template <typename T>
-inline std::array<T, VectorizedArrayType::size()>
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType> &
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
- &array) const
+operator=(const FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType> &other)
{
- std::array<T, VectorizedArrayType::size()> out;
- internal::process_cell_data(
- *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
- local = global;
- });
- return out;
+ this->FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ operator=(other);
+ AssertDimension(n_fe_components, other.n_fe_components);
+ AssertDimension(first_selected_component, other.first_selected_component);
+
+ return *this;
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-template <typename T>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- set_cell_data(
- AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
- const std::array<T, VectorizedArrayType::size()> & 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<dim>(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);
}
// 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)
{
"FEEvaluation. In that case, you must pass an "
"std::vector<VectorType> 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);
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;
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<unsigned int, unsigned int> *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
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<unsigned int, unsigned int> *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;
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;
}
for (; index_indicators != next_index_indicators; ++index_indicators)
{
const std::pair<unsigned short, unsigned short> 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)
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],
++index_indicators)
{
const std::pair<unsigned short, unsigned short> 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)
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],
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;
}
}
}
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]);
}
internal::is_vectorizable<VectorType, Number>::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<unsigned int> &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;
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,
}
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);
}
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,
}
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,
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)
{
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],
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]);
}
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)
{
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],
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]);
/*------------------------------ access to data fields ----------------------*/
-template <int dim,
- int n_components,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline const std::vector<unsigned int> &
-FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
- get_internal_dof_numbering() const
-{
- return data->lexicographic_numbering;
-}
-
-
-
-template <int dim,
- int n_components,
- typename Number,
- bool is_face,
- typename VectorizedArrayType>
-inline ArrayView<VectorizedArrayType>
-FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
- get_scratch_data() const
-{
- return ArrayView<VectorizedArrayType>(
- const_cast<VectorizedArrayType *>(scratch_data),
- scratch_data_array->end() - scratch_data);
-}
-
template <int dim,
# endif
AssertIndexRange(q_point, this->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;
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)
{
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;
# 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;
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++)
# 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;
{
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++)
{
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;
}
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;
{
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)
{
typename Number,
bool is_face,
typename VectorizedArrayType>
-template <int n_components_other>
inline FEEvaluationAccess<dim,
n_components_,
Number,
is_face,
VectorizedArrayType>::
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other)
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other)
: FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
mapping,
fe,
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other)
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other)
: FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
mapping,
fe,
template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
- FEEvaluationAccess(const Mapping<dim> & mapping,
- const FiniteElement<dim> &fe,
- const Quadrature<1> & quadrature,
- const UpdateFlags update_flags,
- const unsigned int first_selected_component,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- is_face,
- VectorizedArrayType> *other)
+ FEEvaluationAccess(
+ const Mapping<dim> & mapping,
+ const FiniteElement<dim> &fe,
+ const Quadrature<1> & quadrature,
+ const UpdateFlags update_flags,
+ const unsigned int first_selected_component,
+ const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+ *other)
: FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
mapping,
fe,
template <typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
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,
quadrature,
update_flags,
first_selected_component,
- static_cast<
- FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
- 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)
quadrature,
update_flags,
first_selected_component,
- static_cast<
- FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
- 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)
int n_components_,
typename Number,
typename VectorizedArrayType>
-template <int n_components_other>
inline FEEvaluation<dim,
fe_degree,
n_q_points_1d,
n_components_,
Number,
VectorizedArrayType>::
- FEEvaluation(const FiniteElement<dim> & fe,
- const FEEvaluationBase<dim,
- n_components_other,
- Number,
- false,
- VectorizedArrayType> &other,
- const unsigned int first_selected_component)
+ FEEvaluation(
+ const FiniteElement<dim> & fe,
+ const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &other,
+ const unsigned int first_selected_component)
: BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
fe,
other.mapped_geometry->get_quadrature(),
projects / dealii.git / commitdiff