#include <deal.II/matrix_free/evaluation_kernels.h>
#include <deal.II/matrix_free/evaluation_selector.h>
#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
#include <deal.II/matrix_free/mapping_data_on_the_fly.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/shape_info.h>
-namespace internal
-{
- DeclException0(ExcAccessToUninitializedField);
-
- DeclException1(
- ExcMatrixFreeAccessToUninitializedMappingField,
- std::string,
- << "You are requesting information from an FEEvaluation/FEFaceEvaluation "
- << "object for which this kind of information has not been computed. What "
- << "information these objects compute is determined by the update_* flags you "
- << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, "
- << "the operation you are attempting requires the <" << arg1
- << "> flag to be set, but it was apparently not specified "
- << "upon initialization.");
-} // namespace internal
-
-template <int dim,
- int fe_degree,
- int n_q_points_1d = fe_degree + 1,
- int n_components_ = 1,
- typename Number = double,
- typename VectorizedArrayType = VectorizedArray<Number>>
-class FEEvaluation;
-
-
-
-/**
- * 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 four template arguments:
- *
- * @tparam dim Dimension in which this class is to be used
- *
- * @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, typename Number, bool is_face, typename VectorizedArrayType>
-class FEEvaluationBaseData
-{
- static_assert(
- std::is_same<Number, typename VectorizedArrayType::value_type>::value,
- "Type of Number and of VectorizedArrayType do not match.");
-
-public:
- static constexpr unsigned int dimension = dim;
-
- /**
- * Destructor.
- */
- ~FEEvaluationBaseData();
-
- /**
- * Return the index offset within the geometry fields for the cell the @p
- * reinit() function has been called for. This index can be used to access
- * an index into a field that has the same compression behavior as the
- * Jacobian of the geometry, e.g., to store an effective coefficient tensors
- * that combines a coefficient with the geometry for lower memory transfer
- * as the available data fields.
- */
- unsigned int
- get_mapping_data_index_offset() const;
-
- /**
- * Return the type of the cell the @p reinit() function has been called for.
- * Valid values are @p cartesian for Cartesian cells (which allows for
- * considerable data compression), @p affine for cells with affine mappings,
- * and @p general for general cells without any compressed storage applied.
- */
- internal::MatrixFreeFunctions::GeometryType
- get_cell_type() const;
-
- /**
- * Return a reference to the ShapeInfo object currently in use.
- */
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
- get_shape_info() const;
-
- /**
- * Return a reference to the DoFInfo object currently in use.
- */
- const internal::MatrixFreeFunctions::DoFInfo &
- get_dof_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 $J^{-\mathrm T}$ of the
- * Jacobian of the mapping between the unit to the real cell defined as
- * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor
- * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space
- * coordinates and rows to real cell coordinates. Thus, the returned tensor
- * represents a covariant transformation, which is used in the
- * FEEvaluationBase::get_gradient() function to transform the unit cell
- * gradients to gradients on the real cell by a multiplication $J^{-\mathrm
- * T} \hat{\nabla} u_h$.
- */
- 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_cell_batches() +
- * MatrixFree::n_ghost_cell_batches() 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_cell_batches() +
- * MatrixFree::n_ghost_cell_batches() 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;
-
- /**
- * Return the number of the quadrature formula of the present cell.
- */
- unsigned int
- get_quadrature_index() const;
-
- /**
- * Return index of the current cell or face.
- */
- unsigned int
- get_current_cell_index() const;
-
- /**
- * Return the active FE index for this class for efficient indexing in the
- * hp-case.
- */
- unsigned int
- get_active_fe_index() const;
-
- /**
- * Return the active quadrature index for this class for efficient indexing in
- * the hp-case.
- */
- unsigned int
- get_active_quadrature_index() const;
-
- /**
- * Return the underlying MatrixFree object.
- */
- const MatrixFree<dim, Number, VectorizedArrayType> &
- get_matrix_free() 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,
- const unsigned int active_fe_index,
- const unsigned int active_quad_index,
- const unsigned int face_type);
-
- /**
- * 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;
-
- /**
- * 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;
-
- /**
- * 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;
-
- /**
- * A pointer to the underlying quadrature formula specified at construction.
- * Also contained in matrix_info, but it simplifies code if we store a
- * reference to it.
- */
- const typename internal::MatrixFreeFunctions::MappingInfoStorage<
- (is_face ? dim - 1 : dim),
- dim,
- Number,
- VectorizedArrayType>::QuadratureDescriptor *descriptor;
-
- /**
- * The number of quadrature points in the current evaluation context.
- */
- const unsigned int n_quadrature_points;
-
- /**
- * 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. For cells it defines if
- * the dof values should be read from the actual cell corresponding to the
- * interior face or the neighboring cell corresponding to the exterior face.
- */
- 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
template <int dim,
int n_components_,
typename Number,
- bool is_face = false,
- typename VectorizedArrayType = VectorizedArray<Number>>
+ bool is_face,
+ typename VectorizedArrayType>
class FEEvaluationBase
: public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
{
//@}
+ /**
+ * Provides a unified interface to access data in a vector of
+ * VectorizedArray fields of length MatrixFree::n_cell_batches() +
+ * MatrixFree::n_ghost_cell_batches() 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_cell_batches() +
+ * MatrixFree::n_ghost_cell_batches() 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 first selected component.
*/
unsigned int
get_first_selected_component() const;
+ /**
+ * Return the underlying MatrixFree object.
+ */
+ const MatrixFree<dim, Number, VectorizedArrayType> &
+ get_matrix_free() const;
+
protected:
/**
* Constructor. Made protected to prevent users from directly using this
FEEvaluationBase &
operator=(const FEEvaluationBase &other);
+ /**
+ * Destructor.
+ */
+ ~FEEvaluationBase();
+
/**
* A unified function to read from and write into vectors based on the given
* template operation. It can perform the operation for @p read_dof_values,
void
apply_hanging_node_constraints() const;
+ /**
+ * This is the general array for all data fields.
+ */
+ AlignedVector<VectorizedArrayType> *scratch_data_array;
+
+ /**
+ * A pointer to the underlying data.
+ */
+ const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
+
/**
* This field stores the values for local degrees of freedom (e.g. after
* reading out from a vector but before applying unit cell transformations
#ifndef DOXYGEN
-/*----------------------- FEEvaluationBaseData ------------------------*/
+namespace internal
+{
+ template <typename MatrixFreeType>
+ inline unsigned int
+ get_active_fe_index(const MatrixFreeType &data,
+ const unsigned int dof_no,
+ const unsigned int fe_degree,
+ const unsigned int first_selected_component,
+ const unsigned int active_fe_index)
+ {
+ return fe_degree != numbers::invalid_unsigned_int ?
+ data.get_dof_info(dof_no).fe_index_from_degree(
+ first_selected_component, fe_degree) :
+ (active_fe_index != numbers::invalid_unsigned_int ?
+ active_fe_index :
+ 0);
+ }
-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,
- const unsigned int active_fe_index_in,
- const unsigned int active_quad_index_in,
- const unsigned int face_type)
- : scratch_data_array(data_in.acquire_scratch_data())
- , quad_no(quad_no_in)
- , matrix_info(&data_in)
- , dof_info(&data_in.get_dof_info(dof_no))
- , mapping_data(
+
+
+ template <bool is_face,
+ int dim,
+ typename Number,
+ typename VectorizedArrayType>
+ inline unsigned int
+ get_active_quad_index(const MatrixFree<dim, Number, VectorizedArrayType> &mf,
+ const unsigned int quad_no,
+ const unsigned int fe_degree,
+ const unsigned int n_q_points,
+ const unsigned int active_fe_index,
+ const unsigned int active_quad_index)
+ {
+ const auto &mapping_storage = internal::MatrixFreeFunctions::
+ MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
+ mf.get_mapping_info(), quad_no);
+ if (fe_degree == numbers::invalid_unsigned_int)
+ return active_quad_index != numbers::invalid_unsigned_int ?
+ active_quad_index :
+ std::min<unsigned int>(active_fe_index,
+ mapping_storage.descriptor.size() - 1);
+ else
+ return mapping_storage.quad_index_from_n_q_points(n_q_points);
+ }
+} // namespace internal
+
+
+/*----------------------- 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,
+ const unsigned int fe_degree,
+ const unsigned int n_q_points,
+ const bool is_interior_face,
+ const unsigned int active_fe_index,
+ const unsigned int active_quad_index,
+ const unsigned int face_type)
+ : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+ data_in.get_shape_info(
+ dof_no,
+ quad_no,
+ data_in.get_dof_info(dof_no)
+ .component_to_base_index[first_selected_component],
+ internal::get_active_fe_index(data_in,
+ dof_no,
+ fe_degree,
+ first_selected_component,
+ active_fe_index),
+ internal::get_active_quad_index<is_face>(
+ data_in,
+ quad_no,
+ fe_degree,
+ n_q_points,
+ internal::get_active_fe_index(data_in,
+ dof_no,
+ fe_degree,
+ first_selected_component,
+ active_fe_index),
+ active_quad_index)),
+ data_in.get_dof_info(dof_no),
internal::MatrixFreeFunctions::
MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
data_in.get_mapping_info(),
- quad_no))
- , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
- data_in.get_dof_info(dof_no).fe_index_from_degree(
- first_selected_component,
- fe_degree) :
- (active_fe_index_in != numbers::invalid_unsigned_int ?
- active_fe_index_in :
- 0))
- , active_quad_index(
- fe_degree != numbers::invalid_unsigned_int ?
- (mapping_data->quad_index_from_n_q_points(n_q_points)) :
- (active_quad_index_in != numbers::invalid_unsigned_int ?
- active_quad_index_in :
- std::min<unsigned int>(active_fe_index,
- mapping_data->descriptor.size() - 1)))
- , descriptor(
- &mapping_data->descriptor
- [is_face ?
- (active_quad_index * std::max<unsigned int>(1, dim - 1) +
- (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
- active_quad_index])
- , n_quadrature_points(descriptor->n_q_points)
- , data(&data_in.get_shape_info(
- dof_no,
- quad_no_in,
- dof_info->component_to_base_index[first_selected_component],
- active_fe_index,
- active_quad_index))
- , jacobian(nullptr)
- , J_value(nullptr)
- , normal_vectors(nullptr)
- , normal_x_jacobian(nullptr)
- , quadrature_weights(descriptor->quadrature_weights.begin())
- , cell(numbers::invalid_unsigned_int)
- , is_interior_face(is_interior_face)
- , dof_access_index(
- is_face ?
- (is_interior_face ?
- internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
- internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
- internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
- , cell_type(internal::MatrixFreeFunctions::general)
-{
- Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
- AssertDimension(matrix_info->get_task_info().vectorization_length,
- VectorizedArrayType::size());
- AssertDimension(n_quadrature_points, descriptor->n_q_points);
+ quad_no),
+ quad_no,
+ is_interior_face,
+ internal::get_active_fe_index(data_in,
+ dof_no,
+ fe_degree,
+ first_selected_component,
+ active_fe_index),
+ internal::get_active_quad_index<is_face>(
+ data_in,
+ quad_no,
+ fe_degree,
+ n_q_points,
+ internal::get_active_fe_index(data_in,
+ dof_no,
+ fe_degree,
+ first_selected_component,
+ active_fe_index),
+ active_quad_index),
+ face_type)
+ , scratch_data_array(data_in.acquire_scratch_data())
+ , matrix_info(&data_in)
+ , 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) + ")"));
+
+ // do not check for correct dimensions of data fields here, should be done
+ // in derived classes
}
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- FEEvaluationBaseData(
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::
+ FEEvaluationBase(
const Mapping<dim> & mapping,
const FiniteElement<dim> &fe,
const Quadrature<1> & quadrature,
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)
- , active_fe_index(numbers::invalid_unsigned_int)
- , active_quad_index(numbers::invalid_unsigned_int)
- , descriptor(nullptr)
- , n_quadrature_points(
- Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
+ : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+ other != nullptr &&
+ other->mapped_geometry->get_quadrature() == quadrature ?
+ other->mapped_geometry :
+ std::make_shared<
+ internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
+ mapping,
+ quadrature,
+ update_flags))
+ , scratch_data_array(new AlignedVector<VectorizedArrayType>())
, matrix_info(nullptr)
- , dof_info(nullptr)
- , mapping_data(nullptr)
- ,
- // select the correct base element from the given FE component
- data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- Quadrature<dim - is_face>(quadrature),
- fe,
- fe.component_to_base_index(first_selected_component).first))
- , jacobian(nullptr)
- , 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)
-{
- Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
- ExcInternalError());
- if (other != nullptr &&
- other->mapped_geometry->get_quadrature() == quadrature)
- mapped_geometry = other->mapped_geometry;
- else
- mapped_geometry =
- std::make_shared<internal::MatrixFreeFunctions::
- MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
- mapping, quadrature, update_flags);
- cell = 0;
+ , 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)
+{
+ 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;
+
+ Assert(this->data == nullptr, ExcInternalError());
+ this->data =
+ new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+ Quadrature<dim - is_face>(quadrature),
+ fe,
+ fe.component_to_base_index(first_selected_component).first);
- 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();
+ set_data_pointers();
}
-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 ?
+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)
+ : 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)
- , active_fe_index(other.active_fe_index)
- , active_quad_index(other.active_quad_index)
- , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
- , n_quadrature_points(other.n_quadrature_points)
, matrix_info(other.matrix_info)
- , dof_info(other.dof_info)
- , mapping_data(other.mapping_data)
- , data(other.matrix_info == nullptr ?
- new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- *other.data) :
- other.data)
- , jacobian(nullptr)
- , J_value(nullptr)
- , normal_vectors(nullptr)
- , normal_x_jacobian(nullptr)
- , quadrature_weights(other.descriptor == nullptr ?
- nullptr :
- descriptor->quadrature_weights.begin())
- , cell(numbers::invalid_unsigned_int)
- , cell_type(internal::MatrixFreeFunctions::general)
- , is_interior_face(other.is_interior_face)
- , dof_access_index(other.dof_access_index)
-{
- // Create deep copy of mapped geometry for use in parallel...
- if (other.mapped_geometry.get() != nullptr)
+ , 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)
+{
+ if (other.matrix_info == nullptr)
{
- mapped_geometry = std::make_shared<
+ Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
+ this->data =
+ new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+ *other.data);
+
+ // Create deep copy of mapped geometry for use in parallel
+ this->mapped_geometry = std::make_shared<
internal::MatrixFreeFunctions::
MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
other.mapped_geometry->get_fe_values().get_mapping(),
other.mapped_geometry->get_quadrature(),
other.mapped_geometry->get_fe_values().get_update_flags());
- mapping_data = &mapped_geometry->get_data_storage();
- cell = 0;
+ this->mapping_data = &this->mapped_geometry->get_data_storage();
+ this->cell = 0;
- jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
- J_value = mapped_geometry->get_data_storage().JxW_values.begin();
+ this->jacobian =
+ this->mapped_geometry->get_data_storage().jacobians[0].begin();
+ this->J_value =
+ this->mapped_geometry->get_data_storage().JxW_values.begin();
}
+
+ set_data_pointers();
}
-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)
+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)
{
- AssertDimension(quad_no, other.quad_no);
- AssertDimension(active_fe_index, other.active_fe_index);
- AssertDimension(active_quad_index, other.active_quad_index);
-
// release old memory
if (matrix_info == nullptr)
{
- delete data;
+ delete this->data;
delete scratch_data_array;
}
else
matrix_info->release_scratch_data(scratch_data_array);
}
- matrix_info = other.matrix_info;
- dof_info = other.dof_info;
- descriptor = other.descriptor;
- mapping_data = other.mapping_data;
+ this->FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ operator=(other);
+
+ matrix_info = other.matrix_info;
+
if (other.matrix_info == nullptr)
{
- data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
- *other.data);
+ Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
+ this->data =
+ new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+ *other.data);
scratch_data_array = new AlignedVector<VectorizedArrayType>();
- }
- else
- {
- data = other.data;
- scratch_data_array = matrix_info->acquire_scratch_data();
- }
- quadrature_weights =
- (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr);
- cell = numbers::invalid_unsigned_int;
- cell_type = internal::MatrixFreeFunctions::general;
- is_interior_face = other.is_interior_face;
- dof_access_index = other.dof_access_index;
-
- // Create deep copy of mapped geometry for use in parallel...
- if (other.mapped_geometry.get() != nullptr)
- {
- mapped_geometry = std::make_shared<
+ // Create deep copy of mapped geometry for use in parallel
+ this->mapped_geometry = std::make_shared<
internal::MatrixFreeFunctions::
MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
other.mapped_geometry->get_fe_values().get_mapping(),
other.mapped_geometry->get_quadrature(),
other.mapped_geometry->get_fe_values().get_update_flags());
- cell = 0;
- 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();
- }
-
- return *this;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- ~FEEvaluationBaseData()
-{
- if (matrix_info != nullptr)
- {
- try
- {
- matrix_info->release_scratch_data(scratch_data_array);
- }
- catch (...)
- {}
- }
- else
- {
- delete scratch_data_array;
- delete data;
- data = nullptr;
+ this->cell = 0;
+ this->mapping_data = &this->mapped_geometry->get_data_storage();
+ this->jacobian =
+ this->mapped_geometry->get_data_storage().jacobians[0].begin();
+ this->J_value =
+ this->mapped_geometry->get_data_storage().JxW_values.begin();
}
- scratch_data_array = nullptr;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_mapping_data_index_offset() const
-{
- if (matrix_info == nullptr)
- return 0;
else
{
- AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
- return this->mapping_data->data_index_offsets[cell];
+ scratch_data_array = matrix_info->acquire_scratch_data();
}
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline internal::MatrixFreeFunctions::GeometryType
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
- const
-{
- Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- return cell_type;
-}
+ AssertDimension(n_fe_components, other.n_fe_components);
+ AssertDimension(first_selected_component, other.first_selected_component);
+ set_data_pointers();
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_shape_info() const
-{
- Assert(data != nullptr, ExcInternalError());
- return *data;
+ return *this;
}
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const internal::MatrixFreeFunctions::DoFInfo &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_dof_info()
- const
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+ n_components_,
+ Number,
+ is_face,
+ VectorizedArrayType>::~FEEvaluationBase()
{
- Assert(dof_info != nullptr,
- ExcMessage(
- "FEEvaluation was not initialized with a MatrixFree object!"));
- return *dof_info;
+ if (matrix_info != nullptr)
+ {
+ try
+ {
+ matrix_info->release_scratch_data(scratch_data_array);
+ }
+ catch (...)
+ {}
+ }
+ else
+ {
+ delete scratch_data_array;
+ delete this->data;
+ }
}
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_normal_vector(const unsigned int q_point) const
+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()
{
- AssertIndexRange(q_point, n_quadrature_points);
- Assert(normal_vectors != nullptr,
- internal::ExcMatrixFreeAccessToUninitializedMappingField(
- "update_normal_vectors"));
- if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
- return normal_vectors[0];
- else
- return normal_vectors[q_point];
-}
+ 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 = this->n_quadrature_points;
+ const unsigned int size_scratch_data =
+ std::max(tensor_dofs_per_component + 1, dofs_per_component) *
+ n_components_ * 3 +
+ 2 * n_quadrature_points;
+ const unsigned int allocated_size =
+ size_scratch_data + n_components_ * dofs_per_component +
+ (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+ n_quadrature_points);
+ this->scratch_data_array->resize_fast(allocated_size);
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
- const unsigned int q_point) const
-{
- AssertIndexRange(q_point, n_quadrature_points);
- Assert(J_value != nullptr,
- internal::ExcMatrixFreeAccessToUninitializedMappingField(
- "update_values|update_gradients"));
- if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ // set the pointers to the correct position in the data array
+ for (unsigned int c = 0; c < n_components_; ++c)
{
- Assert(this->quadrature_weights != nullptr, ExcInternalError());
- return J_value[0] * this->quadrature_weights[q_point];
+ values_dofs[c] =
+ this->scratch_data_array->begin() + c * dofs_per_component;
}
- else
- return J_value[q_point];
+ 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);
+ gradients_from_hessians_quad =
+ this->scratch_data_array->begin() +
+ n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
+ hessians_quad =
+ this->scratch_data_array->begin() +
+ n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+
+ this->scratch_data.reinit(this->scratch_data_array->begin() +
+ n_components_ * dofs_per_component +
+ (n_components_ *
+ ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+ n_quadrature_points),
+ size_scratch_data);
}
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- inverse_jacobian(const unsigned int q_point) const
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
+inline const MatrixFree<dim, Number, VectorizedArrayType> &
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ get_matrix_free() const
{
- AssertIndexRange(q_point, n_quadrature_points);
- Assert(this->jacobian != nullptr,
- internal::ExcMatrixFreeAccessToUninitializedMappingField(
- "update_gradients"));
- if (this->cell_type <= internal::MatrixFreeFunctions::affine)
- return jacobian[0];
- else
- return jacobian[q_point];
+ Assert(matrix_info != nullptr,
+ ExcMessage(
+ "FEEvaluation was not initialized with a MatrixFree object!"));
+ return *matrix_info;
}
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool is_face,
+ typename VectorizedArrayType>
inline std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
- const
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ get_cell_ids() const
{
Assert(this->matrix_info != nullptr, ExcNotInitialized());
{
// cell or interior face face (element-centric loop)
for (unsigned int i = 0; i < n_lanes; ++i)
- cells[i] = cell * n_lanes + i;
+ cells[i] = this->cell * n_lanes + i;
}
else if (is_face &&
this->dof_access_index ==
{
// face-centric faces
const unsigned int *cells_ =
- is_interior_face ?
- &this->matrix_info->get_face_info(cell).cells_interior[0] :
- &this->matrix_info->get_face_info(cell).cells_exterior[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 i = 0; i < VectorizedArrayType::size(); ++i)
if (cells_[i] != numbers::invalid_unsigned_int)
cells[i] = cells_[i];
namespace internal
{
template <int dim,
+ int n_components,
typename Number,
bool is_face,
typename VectorizedArrayType,
typename FU>
inline void
process_cell_data(
- const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &phi,
+ const FEEvaluationBase<dim,
+ n_components,
+ Number,
+ is_face,
+ VectorizedArrayType> & phi,
const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info,
GlobalVectorType & array,
VectorizedArrayType2 & out,
for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
if (cells[i] != numbers::invalid_unsigned_int)
fu(out[i],
- array[cells[i] / VectorizedArrayType::size()]
- [cells[i] % VectorizedArrayType::size()]);
- }
-} // namespace internal
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_cell_or_face_ids() const
-{
- const unsigned int v_len = VectorizedArrayType::size();
- std::array<unsigned int, VectorizedArrayType::size()> cells;
-
- // initialize array
- for (unsigned int i = 0; i < v_len; ++i)
- cells[i] = numbers::invalid_unsigned_int;
-
- if ((is_face &&
- this->dof_access_index ==
- internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
- this->is_interior_face == false) ||
- (!is_face && !this->is_interior_face))
- {
- // cell-based face-loop: plus face
- for (unsigned int i = 0; i < v_len; ++i)
- {
- // compute actual (non vectorized) cell ID
- const unsigned int cell_this = this->cell * v_len + i;
- // compute face ID
- unsigned int fn =
- this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
- this->face_no,
- i);
-
- 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);
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_quadrature_index() const
-{
- return this->quad_no;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_current_cell_index() const
-{
- if (is_face && this->dof_access_index ==
- internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
- return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
- else
- return this->cell;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_active_fe_index() const
-{
- return active_fe_index;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_active_quadrature_index() const
-{
- return active_quad_index;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const MatrixFree<dim, Number, VectorizedArrayType> &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
- get_matrix_free() const
-{
- Assert(matrix_info != nullptr,
- ExcMessage(
- "FEEvaluation was not initialized with a MatrixFree object!"));
- return *matrix_info;
-}
+ array[cells[i] / VectorizedArrayType::size()]
+ [cells[i] % VectorizedArrayType::size()]);
+ }
+} // namespace internal
-/*----------------------- 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,
- const unsigned int active_fe_index,
- const unsigned int active_quad_index,
- const unsigned int face_type)
- : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
- data_in,
- dof_no,
- first_selected_component,
- quad_no_in,
- fe_degree,
- n_q_points,
- is_interior_face,
- active_fe_index,
- active_quad_index,
- face_type)
- , 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)
+std::array<unsigned int, VectorizedArrayType::size()>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ get_cell_or_face_ids() const
{
- 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) + ")"));
+ const unsigned int v_len = VectorizedArrayType::size();
+ std::array<unsigned int, VectorizedArrayType::size()> cells;
- // do not check for correct dimensions of data fields here, should be done
- // in derived classes
+ // initialize array
+ for (unsigned int i = 0; i < v_len; ++i)
+ cells[i] = numbers::invalid_unsigned_int;
+
+ if ((is_face &&
+ this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+ this->is_interior_face == false) ||
+ (!is_face && !this->is_interior_face))
+ {
+ // cell-based face-loop: plus face
+ for (unsigned int i = 0; i < v_len; ++i)
+ {
+ // compute actual (non vectorized) cell ID
+ const unsigned int cell_this = this->cell * v_len + i;
+ // compute face ID
+ unsigned int fn =
+ this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+ this->face_no,
+ i);
+
+ 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] = this->cell * v_len + i;
+ }
+
+ return cells;
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-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)
+inline VectorizedArrayType
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
{
- 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;
+ VectorizedArrayType out = Number(1.);
+ internal::process_cell_data(
+ *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
+ local = global;
+ });
+ return out;
}
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)
- : 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)
+inline void
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+ set_cell_data(AlignedVector<VectorizedArrayType> &array,
+ const VectorizedArrayType & in) const
{
- set_data_pointers();
+ internal::process_cell_data(
+ *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+ global = local;
+ });
}
typename Number,
bool is_face,
typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> &
+template <typename T>
+inline std::array<T, VectorizedArrayType::size()>
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-operator=(const FEEvaluationBase<dim,
- n_components_,
- Number,
- is_face,
- VectorizedArrayType> &other)
+ read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
+ &array) const
{
- 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;
+ 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;
}
typename Number,
bool is_face,
typename VectorizedArrayType>
+template <typename T>
inline void
FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
- set_data_pointers()
+ set_cell_data(
+ AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+ const std::array<T, VectorizedArrayType::size()> & in) const
{
- 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 = this->n_quadrature_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 + 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);
- gradients_from_hessians_quad =
- this->scratch_data_array->begin() +
- n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
- hessians_quad =
- this->scratch_data_array->begin() +
- n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
- this->scratch_data = this->scratch_data_array->begin() +
- n_components_ * dofs_per_component +
- (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
- n_quadrature_points);
+ internal::process_cell_data(
+ *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+ global = local;
+ });
}
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 quad_no,
const unsigned int fe_degree,
const unsigned int n_q_points,
const bool is_interior_face,
data_in,
dof_no,
first_selected_component,
- quad_no_in,
+ quad_no,
fe_degree,
n_q_points,
is_interior_face,
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 quad_no,
const unsigned int fe_degree,
const unsigned int n_q_points,
const bool is_interior_face,
data_in,
dof_no,
first_selected_component,
- quad_no_in,
+ quad_no,
fe_degree,
n_q_points,
is_interior_face,
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 quad_no,
const unsigned int fe_degree,
const unsigned int n_q_points,
const bool is_interior_face,
data_in,
dof_no,
first_selected_component,
- quad_no_in,
+ quad_no,
fe_degree,
n_q_points,
is_interior_face,
FEEvaluationAccess(const MatrixFree<1, Number, VectorizedArrayType> &data_in,
const unsigned int dof_no,
const unsigned int first_selected_component,
- const unsigned int quad_no_in,
+ const unsigned int quad_no,
const unsigned int fe_degree,
const unsigned int n_q_points,
const bool is_interior_face,
data_in,
dof_no,
first_selected_component,
- quad_no_in,
+ quad_no,
fe_degree,
n_q_points,
is_interior_face,
this->values_quad,
this->gradients_quad,
this->hessians_quad,
- this->scratch_data);
+ this->scratch_data.begin());
else
internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
n_components,
this->values_quad,
this->gradients_quad,
this->hessians_quad,
- this->scratch_data);
+ this->scratch_data.begin());
if (hessians_on_general_cells)
for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i)
this->values_quad,
this->gradients_quad,
this->hessians_quad,
- this->scratch_data,
+ this->scratch_data.begin(),
false);
else
internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
this->values_quad,
this->gradients_quad,
this->hessians_quad,
- this->scratch_data,
+ this->scratch_data.begin(),
false);
# ifdef DEBUG
this->values_quad,
this->gradients_quad,
this->hessians_quad,
- this->scratch_data,
+ this->scratch_data.begin(),
this->data,
integration_flag) == false)
{
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag_actual & EvaluationFlags::values,
evaluation_flag_actual & EvaluationFlags::gradients,
evaluation_flag_actual & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag_actual & EvaluationFlags::values,
evaluation_flag_actual & EvaluationFlags::gradients,
evaluation_flag_actual & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag_actual & EvaluationFlags::values,
evaluation_flag_actual & EvaluationFlags::gradients,
evaluation_flag_actual & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag_actual & EvaluationFlags::values,
evaluation_flag_actual & EvaluationFlags::gradients,
evaluation_flag_actual & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag_actual & EvaluationFlags::values,
evaluation_flag_actual & EvaluationFlags::gradients,
evaluation_flag_actual & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag & EvaluationFlags::values,
evaluation_flag & EvaluationFlags::gradients,
evaluation_flag & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag & EvaluationFlags::values,
evaluation_flag & EvaluationFlags::gradients,
evaluation_flag & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag & EvaluationFlags::values,
evaluation_flag & EvaluationFlags::gradients,
evaluation_flag & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag & EvaluationFlags::values,
evaluation_flag & EvaluationFlags::gradients,
evaluation_flag & EvaluationFlags::hessians,
this->begin_values(),
this->begin_gradients(),
this->begin_hessians(),
- this->scratch_data,
+ this->scratch_data.begin(),
evaluation_flag & EvaluationFlags::values,
evaluation_flag & EvaluationFlags::gradients,
evaluation_flag & EvaluationFlags::hessians,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_fe_evaluation_base_data_h
+#define dealii_matrix_free_fe_evaluation_base_data_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/dof_info.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
+#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/type_traits.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace internal
+{
+ DeclException0(ExcAccessToUninitializedField);
+
+ DeclException1(
+ ExcMatrixFreeAccessToUninitializedMappingField,
+ std::string,
+ << "You are requesting information from an FEEvaluation/FEFaceEvaluation "
+ << "object for which this kind of information has not been computed. What "
+ << "information these objects compute is determined by the update_* flags you "
+ << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, "
+ << "the operation you are attempting requires the <" << arg1
+ << "> flag to be set, but it was apparently not specified "
+ << "upon initialization.");
+} // namespace internal
+
+// forward declarations
+template <int dim,
+ int fe_degree,
+ int n_q_points_1d = fe_degree + 1,
+ int n_components_ = 1,
+ typename Number = double,
+ typename VectorizedArrayType = VectorizedArray<Number>>
+class FEEvaluation;
+
+template <int dim,
+ int n_components_,
+ typename Number,
+ bool = false,
+ typename = VectorizedArray<Number>>
+class FEEvaluationBase;
+
+namespace internal
+{
+ namespace MatrixFreeFunctions
+ {
+ template <int, typename, typename>
+ class MappingDataOnTheFly;
+ }
+} // namespace internal
+
+
+/**
+ * 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 four template arguments:
+ *
+ * @tparam dim Dimension in which this class is to be used
+ *
+ * @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, typename Number, bool is_face, typename VectorizedArrayType>
+class FEEvaluationBaseData
+{
+ static_assert(
+ std::is_same<Number, typename VectorizedArrayType::value_type>::value,
+ "Type of Number and of VectorizedArrayType do not match.");
+
+ using ShapeInfoType =
+ internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>;
+ using MappingInfoStorageType =
+ internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 : dim),
+ dim,
+ Number,
+ VectorizedArrayType>;
+ using DoFInfo = internal::MatrixFreeFunctions::DoFInfo;
+
+public:
+ static constexpr unsigned int dimension = dim;
+
+ /**
+ * Return the index offset within the geometry fields for the cell the @p
+ * reinit() function has been called for. This index can be used to access
+ * an index into a field that has the same compression behavior as the
+ * Jacobian of the geometry, e.g., to store an effective coefficient tensors
+ * that combines a coefficient with the geometry for lower memory transfer
+ * as the available data fields.
+ */
+ unsigned int
+ get_mapping_data_index_offset() const;
+
+ /**
+ * Return the type of the cell the @p reinit() function has been called for.
+ * Valid values are @p cartesian for Cartesian cells (which allows for
+ * considerable data compression), @p affine for cells with affine mappings,
+ * and @p general for general cells without any compressed storage applied.
+ */
+ internal::MatrixFreeFunctions::GeometryType
+ get_cell_type() const;
+
+ /**
+ * Return a reference to the ShapeInfo object currently in use.
+ */
+ const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+ get_shape_info() const;
+
+ /**
+ * Return a reference to the DoFInfo object currently in use.
+ */
+ const internal::MatrixFreeFunctions::DoFInfo &
+ get_dof_info() 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;
+
+ /**
+ * 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 $J^{-\mathrm T}$ of the
+ * Jacobian of the mapping between the unit to the real cell defined as
+ * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor
+ * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space
+ * coordinates and rows to real cell coordinates. Thus, the returned tensor
+ * represents a covariant transformation, which is used in the
+ * FEEvaluationBase::get_gradient() function to transform the unit cell
+ * gradients to gradients on the real cell by a multiplication $J^{-\mathrm
+ * T} \hat{\nabla} u_h$.
+ */
+ 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;
+
+ /**
+ * 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 the number of the quadrature formula of the present cell.
+ */
+ unsigned int
+ get_quadrature_index() const;
+
+ /**
+ * Return index of the current cell or face.
+ */
+ unsigned int
+ get_current_cell_index() const;
+
+ /**
+ * Return the active FE index for this class for efficient indexing in the
+ * hp-case.
+ */
+ unsigned int
+ get_active_fe_index() const;
+
+ /**
+ * Return the active quadrature index for this class for efficient indexing in
+ * the hp-case.
+ */
+ unsigned int
+ get_active_quadrature_index() 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 ShapeInfoType & shape_info,
+ const DoFInfo & dof_info,
+ const MappingInfoStorageType &mapping_data,
+ const unsigned int quad_no,
+ const bool is_interior_face,
+ const unsigned int active_fe_index,
+ const unsigned int active_quad_index,
+ const unsigned int face_type);
+
+ /**
+ * Constructor that comes with reduced functionality and works similar as
+ * FEValues.
+ */
+ FEEvaluationBaseData(
+ const std::shared_ptr<
+ internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapping_data);
+
+ /**
+ * Copy constructor.
+ */
+ 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);
+
+ /**
+ * 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 ShapeInfoType *data;
+
+ /**
+ * 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 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 MappingInfoStorageType *mapping_data;
+
+ /**
+ * The number of the quadrature formula of the present cell among all
+ * quadrature formulas available in the MatrixFree objects pass to derived
+ * classes.
+ */
+ const unsigned int quad_no;
+
+ /**
+ * 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;
+
+ /**
+ * A pointer to the underlying quadrature formula specified at construction.
+ * Also contained in mapping_data, but it simplifies code if we store a
+ * reference to it.
+ */
+ const typename MappingInfoStorageType::QuadratureDescriptor *descriptor;
+
+ /**
+ * The number of quadrature points in the current evaluation context.
+ */
+ const unsigned int n_quadrature_points;
+
+ /**
+ * 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;
+
+ /**
+ * This is the user-visible part of FEEvaluationBase::scratch_data_array,
+ * only showing the part that can be consumed by various users. It is set
+ * during the allocation of the internal data structures in
+ * FEEvaluationBase.
+ */
+ mutable ArrayView<VectorizedArrayType> scratch_data;
+
+ /**
+ * 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. For cells it defines if
+ * the dof values should be read from the actual cell corresponding to the
+ * interior face or the neighboring cell corresponding to the exterior face.
+ */
+ 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, as an alternative to the entry point with
+ * MatrixFree.
+ */
+ std::shared_ptr<internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
+ mapped_geometry;
+
+ // Make FEEvaluation and FEEvaluationBase objects friends for access to
+ // protected member mapped_geometry.
+ template <int, int, typename, bool, typename>
+ friend class FEEvaluationBase;
+
+ template <int, int, int, int, typename, typename>
+ friend class FEEvaluation;
+};
+
+
+
+/*----------------------- Inline functions ----------------------------------*/
+
+#ifndef DOXYGEN
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(const ShapeInfoType & shape_info,
+ const DoFInfo & dof_info,
+ const MappingInfoStorageType &mapping_data,
+ const unsigned int quad_no,
+ const bool is_interior_face,
+ const unsigned int active_fe_index,
+ const unsigned int active_quad_index,
+ const unsigned int face_type)
+ : data(&shape_info)
+ , dof_info(&dof_info)
+ , mapping_data(&mapping_data)
+ , quad_no(quad_no)
+ , active_fe_index(active_fe_index)
+ , active_quad_index(active_quad_index)
+ , descriptor(
+ &mapping_data.descriptor
+ [is_face ?
+ (active_quad_index * std::max<unsigned int>(1, dim - 1) +
+ (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
+ active_quad_index])
+ , n_quadrature_points(descriptor->n_q_points)
+ , jacobian(nullptr)
+ , J_value(nullptr)
+ , normal_vectors(nullptr)
+ , normal_x_jacobian(nullptr)
+ , quadrature_weights(descriptor->quadrature_weights.begin())
+ , cell(numbers::invalid_unsigned_int)
+ , is_interior_face(is_interior_face)
+ , dof_access_index(
+ is_face ?
+ (is_interior_face ?
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+ , face_no(0)
+ , face_orientation(0)
+ , subface_index(0)
+ , cell_type(internal::MatrixFreeFunctions::general)
+{}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(
+ const std::shared_ptr<
+ internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapped_geometry)
+ : data(nullptr)
+ , dof_info(nullptr)
+ , mapping_data(nullptr)
+ , quad_no(numbers::invalid_unsigned_int)
+ , active_fe_index(numbers::invalid_unsigned_int)
+ , active_quad_index(numbers::invalid_unsigned_int)
+ , descriptor(nullptr)
+ , n_quadrature_points(
+ mapped_geometry->get_data_storage().descriptor[0].n_q_points)
+ , jacobian(nullptr)
+ , 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)
+ , mapped_geometry(mapped_geometry)
+{
+ 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();
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ FEEvaluationBaseData(const FEEvaluationBaseData &other)
+ : data(other.data)
+ , dof_info(other.dof_info)
+ , mapping_data(other.mapping_data)
+ , quad_no(other.quad_no)
+ , active_fe_index(other.active_fe_index)
+ , active_quad_index(other.active_quad_index)
+ , descriptor(other.descriptor)
+ , n_quadrature_points(other.n_quadrature_points)
+ , jacobian(nullptr)
+ , J_value(nullptr)
+ , normal_vectors(nullptr)
+ , normal_x_jacobian(nullptr)
+ , quadrature_weights(other.quadrature_weights)
+ , cell(numbers::invalid_unsigned_int)
+ , is_interior_face(other.is_interior_face)
+ , dof_access_index(
+ is_face ?
+ (is_interior_face ?
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+ , face_no(0)
+ , face_orientation(0)
+ , subface_index(0)
+ , cell_type(internal::MatrixFreeFunctions::general)
+{}
+
+
+
+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 &other)
+{
+ AssertDimension(quad_no, other.quad_no);
+ AssertDimension(active_fe_index, other.active_fe_index);
+ AssertDimension(active_quad_index, other.active_quad_index);
+ AssertDimension(n_quadrature_points, descriptor->n_q_points);
+
+ data = other.data;
+ dof_info = other.dof_info;
+ mapping_data = other.mapping_data;
+ descriptor = other.descriptor;
+ jacobian = nullptr;
+ J_value = nullptr;
+ normal_vectors = nullptr;
+ normal_x_jacobian = nullptr;
+ quadrature_weights = other.quadrature_weights;
+ cell = numbers::invalid_unsigned_int;
+ is_interior_face = other.is_interior_face;
+ dof_access_index =
+ is_face ?
+ (is_interior_face ?
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
+ face_no = 0;
+ face_orientation = 0;
+ subface_index = 0;
+ cell_type = internal::MatrixFreeFunctions::general;
+
+ return *this;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_mapping_data_index_offset() const
+{
+ Assert(mapping_data != nullptr, ExcInternalError());
+
+ if (dof_info == nullptr)
+ return 0;
+ else
+ {
+ AssertIndexRange(cell, mapping_data->data_index_offsets.size());
+ return mapping_data->data_index_offsets[cell];
+ }
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline internal::MatrixFreeFunctions::GeometryType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
+ const
+{
+ Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+ return cell_type;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_shape_info() const
+{
+ Assert(data != nullptr, ExcInternalError());
+ return *data;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const internal::MatrixFreeFunctions::DoFInfo &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_dof_info()
+ const
+{
+ Assert(dof_info != nullptr,
+ ExcMessage(
+ "FEEvaluation was not initialized with a MatrixFree object!"));
+ return *dof_info;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_normal_vector(const unsigned int q_point) const
+{
+ AssertIndexRange(q_point, n_quadrature_points);
+ Assert(normal_vectors != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_normal_vectors"));
+ if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
+ return normal_vectors[0];
+ else
+ return normal_vectors[q_point];
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
+ const unsigned int q_point) const
+{
+ AssertIndexRange(q_point, n_quadrature_points);
+ Assert(J_value != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_values|update_gradients"));
+ if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ {
+ Assert(this->quadrature_weights != nullptr, ExcInternalError());
+ return J_value[0] * this->quadrature_weights[q_point];
+ }
+ else
+ return J_value[q_point];
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ inverse_jacobian(const unsigned int q_point) const
+{
+ AssertIndexRange(q_point, n_quadrature_points);
+ Assert(this->jacobian != nullptr,
+ internal::ExcMatrixFreeAccessToUninitializedMappingField(
+ "update_gradients"));
+ if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ return jacobian[0];
+ else
+ return jacobian[q_point];
+}
+
+
+
+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 unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_quadrature_index() const
+{
+ return this->quad_no;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_current_cell_index() const
+{
+ if (is_face && this->dof_access_index ==
+ internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+ return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
+ else
+ return this->cell;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_active_fe_index() const
+{
+ return active_fe_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_active_quadrature_index() const
+{
+ return active_quad_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline ArrayView<VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+ get_scratch_data() const
+{
+ return scratch_data;
+}
+
+
+
+#endif // ifndef DOXYGEN
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
projects / dealii.git / commitdiff