template <int dim, int spacedim>
class ComputeMappingDataHelper;
}
+ template <int dim, int spacedim, typename Number>
+ class MappingInfo;
} // namespace NonMatching
*/
InternalDataBase(const InternalDataBase &) = delete;
+ /**
+ * This function initializes the data fields related to evaluation of the
+ * mapping on cells, implemented by (derived) classes. This function is
+ * used both when setting up a field of this class for the first time or
+ * when a new Quadrature formula should be considered without creating an
+ * entirely new object. This is used when the number of evaluation points
+ * is different on each cell, e.g. when using FEPointEvaluation for
+ * handling particles or with certain non-matching problem settings.
+ */
+ virtual void
+ reinit(const UpdateFlags update_flags, const Quadrature<dim> &quadrature);
+
/**
* Virtual destructor for derived classes
*/
memory_consumption() const;
};
-
protected:
/**
* Given a set of update flags, compute which other quantities <i>also</i>
friend class FESubfaceValues<dim, spacedim>;
friend class NonMatching::FEImmersedSurfaceValues<dim>;
friend class NonMatching::internal::ComputeMappingDataHelper<dim, spacedim>;
+ template <int, int, typename>
+ friend class NonMatching::MappingInfo;
};
*/
InternalData(const Quadrature<dim> &quadrature);
+ // Documentation see Mapping::InternalDataBase.
+ virtual void
+ reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature) override;
+
/**
* Return an estimate (in bytes) for the memory consumption of this object.
*/
mutable double volume_element;
/**
- * Vector of all quadrature points. Especially, all points on all faces.
+ * Location of quadrature points of faces or subfaces in 3d with all
+ * possible orientations. Can be accessed with the correct offset provided
+ * via QProjector::DataSetDescriptor. Not needed/used for cells.
*/
std::vector<Point<dim>> quadrature_points;
};
maybe_update_cell_quadrature_points(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const InternalData &data,
- std::vector<Point<dim>> &quadrature_points) const;
+ const ArrayView<const Point<dim>> &unit_quadrature_points,
+ std::vector<Point<dim>> &quadrature_points) const;
/**
* Compute the quadrature points if the UpdateFlags of the incoming
const InternalData &data,
std::vector<Point<dim>> &quadrature_points) const;
- /**
- * Transform quadrature points in InternalData to real space by scaling unit
- * coordinates with cell_extents in each direction.
- *
- * Called from the various maybe_update_*_quadrature_points functions.
- */
- void
- transform_quadrature_points(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const InternalData &data,
- const typename QProjector<dim>::DataSetDescriptor &offset,
- std::vector<Point<dim>> &quadrature_points) const;
-
/**
* Compute the normal vectors if the UpdateFlags of the incoming InternalData
* object say that they should be updated.
*/
InternalData(const FiniteElement<dim, spacedim> &fe);
- /**
- * Initialize the object's member variables related to cell data based on
- * the given arguments.
- *
- * The function also calls compute_shape_function_values() to actually set
- * the member variables related to the values and derivatives of the
- * mapping shape functions.
- */
- void
- initialize(const UpdateFlags update_flags,
- const Quadrature<dim> &quadrature,
- const unsigned int n_original_q_points);
+ // Documentation see Mapping::InternalDataBase.
+ virtual void
+ reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature) override;
/**
* Initialize the object's member variables related to cell and face data
* based on the given arguments. In order to initialize cell data, this
- * function calls initialize().
+ * function calls reinit().
*/
void
initialize_face(const UpdateFlags update_flags,
InternalData(const FiniteElement<dim, spacedim> &fe,
const ComponentMask &mask);
+ // Documentation see Mapping::InternalDataBase.
+ virtual void
+ reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature) override;
+
/**
* Shape function at quadrature point. Shape functions are in tensor
* product order, so vertices must be reordered to obtain transformation.
virtual std::size_t
memory_consumption() const override;
+ /**
+ * A pointer to the underlying finite element.
+ */
+ SmartPointer<const FiniteElement<dim, spacedim>> fe;
+
/**
* Values of shape functions. Access by function @p shape.
*
Point<spacedim>
do_transform_unit_to_real_cell(const InternalData &mdata) const;
-
/**
* Transform the point @p p on the real cell to the corresponding point on
* the unit cell @p cell by a Newton iteration.
* Takes a reference to an @p InternalData that is assumed to be previously
* created by the @p get_data function with @p UpdateFlags including @p
* update_transformation_values and @p update_transformation_gradients and a
- * one point Quadrature that includes the given initial guess for the
- * transformation @p initial_p_unit. Hence this function assumes that @p
- * mdata already includes the transformation shape values and gradients
- * computed at @p initial_p_unit.
+ * one point Quadrature that includes the given initial guess specified
+ * through the given @p point_quadrature object.
*
* @p mdata will be changed by this function.
*/
do_transform_real_to_unit_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<spacedim> &p,
- const Point<dim> &initial_p_unit,
- InternalData &mdata) const;
+ Quadrature<dim> &point_quadrature,
+ InternalData &mdata) const;
/**
* Update internal degrees of freedom.
const typename MappingFEField<dim, spacedim, VectorType>::InternalData
&data) const;
- /**
- * See the documentation of the base class for detailed information.
- */
- virtual void
- compute_shapes_virtual(
- const std::vector<Point<dim>> &unit_points,
- typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
- const;
-
/*
* Which components to use for the mapping.
*/
mutable Threads::Mutex fe_values_mutex;
void
- compute_data(const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points,
- InternalData &data) const;
-
- void
- compute_face_data(const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points,
- InternalData &data) const;
-
+ compute_face_data(const unsigned int n_original_q_points,
+ InternalData &data) const;
// Declare other MappingFEField classes friends.
template <int, int, class>
*/
InternalData() = default;
- /**
- * Initialize the object's member variables related to cell data based on
- * the given arguments.
- *
- * The function also calls compute_shape_function_values() to actually set
- * the member variables related to the values and derivatives of the
- * mapping shape functions.
- */
- void
- initialize(const UpdateFlags update_flags,
- const Quadrature<dim> &quadrature,
- const unsigned int n_original_q_points);
+ // Documentation see Mapping::InternalDataBase.
+ virtual void
+ reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature) override;
/**
* Initialize the object's member variables related to cell and face data
#ifndef DOXYGEN
-template <int dim, int spacedim>
-inline void
-MappingManifold<dim, spacedim>::InternalData::store_vertices(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
-{
- vertices.resize(GeometryInfo<dim>::vertices_per_cell);
- for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
- vertices[i] = cell->vertex(i);
- this->cell = cell;
-}
-
-
-template <int dim, int spacedim>
-inline void
-MappingManifold<dim, spacedim>::InternalData::
- compute_manifold_quadrature_weights(const Quadrature<dim> &quad)
-{
- cell_manifold_quadrature_weights.resize(
- quad.size(), std::vector<double>(GeometryInfo<dim>::vertices_per_cell));
- for (unsigned int q = 0; q < quad.size(); ++q)
- {
- for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
- {
- cell_manifold_quadrature_weights[q][i] =
- GeometryInfo<dim>::d_linear_shape_function(quad.point(q), i);
- }
- }
-}
-
-
-
template <int dim, int spacedim>
inline bool
MappingManifold<dim, spacedim>::preserves_vertex_locations() const
*/
InternalData(const unsigned int polynomial_degree);
- /**
- * Initialize the object's member variables related to cell data based on
- * the given arguments.
- */
- void
- initialize(const UpdateFlags update_flags,
- const Quadrature<dim> &quadrature,
- const unsigned int n_original_q_points);
+ // Documentation see Mapping::InternalDataBase.
+ virtual void
+ reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature) override;
/**
* Initialize the object's member variables related to cell and face data
MappingData &mapping_data)
{
mapping_data.initialize(quadrature.size(), update_flags_mapping);
-
- // reuse internal_mapping_data for MappingQ to avoid memory allocations
- if (const MappingQ<dim, spacedim> *mapping_q =
- dynamic_cast<const MappingQ<dim, spacedim> *>(mapping.get()))
- {
- (void)mapping_q;
- auto &data =
- static_cast<typename MappingQ<dim, spacedim>::InternalData &>(
- *internal_mapping_data);
- data.initialize(update_flags_mapping,
- quadrature,
- quadrature.size());
- }
- else
- {
- internal_mapping_data =
- mapping->get_data(update_flags_mapping, quadrature);
- }
+ internal_mapping_data->reinit(update_flags_mapping, quadrature);
cell_similarity = mapping->fill_fe_values(cell,
cell_similarity,
{
mapping_data.initialize(quadrature.size(), update_flags_mapping);
- // reuse internal_mapping_data for MappingQ to avoid memory allocations
- if (dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
- {
- auto &data =
- dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
- *internal_mapping_data);
- data.initialize(update_flags_mapping,
- quadrature,
- quadrature.size());
- }
- else
- {
- internal_mapping_data =
- mapping->get_data(update_flags_mapping, quadrature);
- }
+ internal_mapping_data->reinit(update_flags_mapping, quadrature);
mapping->fill_fe_immersed_surface_values(cell,
quadrature,
* template argument, e.g. VectorizedArray<double>) provides both mapping data
* and unit points in vectorized format. The Number template parameter of
* MappingInfo and FEPointEvaluation has to be identical.
+ *
+ * @note This class cannot be copied or copy-constructed. The intention of
+ * this class is to be available as a single use, e.g. inside
+ * FEPointEvaluation or as a storage container for the whole mesh. Use a
+ * shared or unique pointer of this class in case several objects should be
+ * held.
*/
template <int dim, int spacedim = dim, typename Number = double>
class MappingInfo : public Subscriptor
const UpdateFlags update_flags,
const AdditionalData additional_data = AdditionalData());
+ /**
+ * Do not allow making copies.
+ */
+ MappingInfo(const MappingInfo &) = delete;
+
+ /**
+ * Do not allow copy assignment of this class.
+ */
+ MappingInfo &
+ operator=(const MappingInfo &) = delete;
+
/**
* Compute the mapping information for the incoming cell and unit
* points. This overload is needed to resolve ambiguity.
Quadrature<dim> quadrature;
/**
- * An array containing the data fields of this class.
+ * Helper class that temporarily holds the data requested for one cell by
+ * Mapping::fill_fe_values() before it is filled into appropriate data
+ * structures for consumption by FEPointEvaluation.
*/
MappingData mapping_data;
internal::ComputeMappingDataHelper<dim, spacedim>::required_update_flags(
this->mapping, update_flags_mapping);
- // construct internal_mapping_data for MappingQ to be able to reuse it in
- // reinit() calls to avoid memory allocations
- if (const MappingQ<dim, spacedim> *mapping_q =
- dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
- {
- internal_mapping_data =
- std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
- mapping_q->get_degree());
- }
+ // construct internal_mapping_data for mappings for reuse in reinit()
+ // calls to avoid frequent memory allocations
+ internal_mapping_data = mapping.get_data(update_flags, Quadrature<dim>());
}
const typename Mapping<1, spacedim>::InternalDataBase &mapping_data,
Table<2, Number> &value_list)
{
- unsigned int n_q_points;
- double cell_extent = 1.0;
+ double cell_extent = 1.0;
// Check mapping_data is associated with a compatible mapping class
if (dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
const typename MappingCartesian<1>::InternalData *data =
dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
&mapping_data);
- n_q_points = data->quadrature_points.size();
cell_extent = data->cell_extents[0];
}
else
AssertDimension(value_list.size(0), n_dofs_per_cell);
AssertDimension(n_dofs_per_cell, 2 * regularity + 2);
- AssertIndexRange(n_q_points_out, n_q_points + 1);
- (void)n_q_points;
std::vector<unsigned int> l2h =
hermite_lexicographic_to_hierarchic_numbering<1>(regularity);
const typename Mapping<2, spacedim>::InternalDataBase &mapping_data,
Table<2, Number> &value_list)
{
- unsigned int n_q_points;
Tensor<1, 2> cell_extents;
// Check mapping_data is associated with a compatible mapping class
const typename MappingCartesian<2>::InternalData *data =
dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
&mapping_data);
- n_q_points = data->quadrature_points.size();
cell_extents = data->cell_extents;
}
else
AssertDimension(value_list.size(0), n_dofs_per_cell);
AssertDimension(n_dofs_per_dim * n_dofs_per_dim, n_dofs_per_cell);
- AssertIndexRange(n_q_points_out, n_q_points + 1);
- (void)n_q_points;
std::vector<unsigned int> l2h =
hermite_lexicographic_to_hierarchic_numbering<2>(regularity);
const typename Mapping<3, spacedim>::InternalDataBase &mapping_data,
Table<2, Number> &value_list)
{
- unsigned int n_q_points;
Tensor<1, 3> cell_extents;
// Check mapping_data is associated with a compatible mapping class
const typename MappingCartesian<3>::InternalData *data =
dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
&mapping_data);
- n_q_points = data->quadrature_points.size();
cell_extents = data->cell_extents;
}
else
AssertDimension(value_list.size(0), n_dofs_per_cell);
AssertDimension(Utilities::pow(n_dofs_per_dim, 3), n_dofs_per_cell);
- AssertIndexRange(n_q_points_out, n_q_points + 1);
- (void)n_q_points;
std::vector<unsigned int> l2h =
hermite_lexicographic_to_hierarchic_numbering<3>(regularity);
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::InternalDataBase::reinit(const UpdateFlags,
+ const Quadrature<dim> &)
+{
+ DEAL_II_ASSERT_UNREACHABLE();
+}
+
+
+
template <int dim, int spacedim>
std::size_t
Mapping<dim, spacedim>::InternalDataBase::memory_consumption() const
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::InternalData::reinit(
+ const UpdateFlags update_flags,
+ const Quadrature<dim> &)
+{
+ // store the flags in the internal data object so we can access them
+ // in fill_fe_*_values(). use the transitive hull of the required
+ // flags
+ this->update_each = update_flags;
+}
+
+
+
template <int dim, int spacedim>
std::size_t
MappingCartesian<dim, spacedim>::InternalData::memory_consumption() const
return (Mapping<dim, spacedim>::InternalDataBase::memory_consumption() +
MemoryConsumption::memory_consumption(cell_extents) +
MemoryConsumption::memory_consumption(inverse_cell_extents) +
- MemoryConsumption::memory_consumption(volume_element) +
- MemoryConsumption::memory_consumption(quadrature_points));
+ MemoryConsumption::memory_consumption(volume_element));
}
const Quadrature<dim> &q) const
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
- std::make_unique<InternalData>(q);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
-
- // store the flags in the internal data object so we can access them
- // in fill_fe_*_values(). use the transitive hull of the required
- // flags
- data.update_each = requires_update_flags(update_flags);
+ std::make_unique<InternalData>();
+ data_ptr->reinit(requires_update_flags(update_flags), q);
return data_ptr;
}
+namespace
+{
+ template <int dim>
+ void
+ transform_quadrature_points(
+ const Tensor<1, dim> first_vertex,
+ const Tensor<1, dim> cell_extents,
+ const ArrayView<const Point<dim>> &unit_quadrature_points,
+ const typename QProjector<dim>::DataSetDescriptor &offset,
+ std::vector<Point<dim>> &quadrature_points)
+ {
+ for (unsigned int i = 0; i < quadrature_points.size(); ++i)
+ {
+ quadrature_points[i] = first_vertex;
+ for (unsigned int d = 0; d < dim; ++d)
+ quadrature_points[i][d] +=
+ cell_extents[d] * unit_quadrature_points[i + offset][d];
+ }
+ }
+} // namespace
+
+
+
template <int dim, int spacedim>
void
MappingCartesian<dim, spacedim>::maybe_update_cell_quadrature_points(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const InternalData &data,
- std::vector<Point<dim>> &quadrature_points) const
+ const ArrayView<const Point<dim>> &unit_quadrature_points,
+ std::vector<Point<dim>> &quadrature_points) const
{
if (data.update_each & update_quadrature_points)
{
const auto offset = QProjector<dim>::DataSetDescriptor::cell();
- transform_quadrature_points(cell, data, offset, quadrature_points);
+ transform_quadrature_points(cell->vertex(0),
+ data.cell_extents,
+ unit_quadrature_points,
+ offset,
+ quadrature_points);
}
}
quadrature_points.size());
- transform_quadrature_points(cell, data, offset, quadrature_points);
+ transform_quadrature_points(cell->vertex(0),
+ data.cell_extents,
+ make_array_view(data.quadrature_points),
+ offset,
+ quadrature_points);
}
}
quadrature_points.size(),
cell->subface_case(face_no));
- transform_quadrature_points(cell, data, offset, quadrature_points);
- }
-}
-
-
-
-template <int dim, int spacedim>
-void
-MappingCartesian<dim, spacedim>::transform_quadrature_points(
- const typename Triangulation<dim, spacedim>::cell_iterator &cell,
- const InternalData &data,
- const typename QProjector<dim>::DataSetDescriptor &offset,
- std::vector<Point<dim>> &quadrature_points) const
-{
- // let @p{start} be the origin of a local coordinate system. it is chosen
- // as the (lower) left vertex
- const Point<dim> start = cell->vertex(0);
-
- for (unsigned int i = 0; i < quadrature_points.size(); ++i)
- {
- quadrature_points[i] = start;
- for (unsigned int d = 0; d < dim; ++d)
- quadrature_points[i][d] +=
- data.cell_extents[d] * data.quadrature_points[i + offset][d];
+ transform_quadrature_points(cell->vertex(0),
+ data.cell_extents,
+ make_array_view(data.quadrature_points),
+ offset,
+ quadrature_points);
}
}
maybe_update_cell_quadrature_points(cell,
data,
+ quadrature.get_points(),
output_data.quadrature_points);
// compute Jacobian determinant. all values are equal and are the
InternalData data;
data.update_each = update_flags;
- data.quadrature_points =
- std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
update_cell_extents(cell, CellSimilarity::none, data);
maybe_update_cell_quadrature_points(cell,
data,
+ unit_points,
output_data.quadrature_points);
maybe_update_jacobians(data, CellSimilarity::none, output_data);
maybe_update_cell_quadrature_points(cell,
data,
+ quadrature.get_points(),
output_data.quadrature_points);
if (data.update_each & update_normal_vectors)
}
+
template <int dim, int spacedim>
void
-MappingFE<dim, spacedim>::InternalData::initialize(
- const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points)
+MappingFE<dim, spacedim>::InternalData::reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &q)
{
// store the flags in the internal data object so we can access them
// in fill_fe_*_values()
const unsigned int n_q_points = q.size();
if (this->update_each & update_covariant_transformation)
- covariant.resize(n_original_q_points);
+ covariant.resize(n_q_points);
if (this->update_each & update_contravariant_transformation)
- contravariant.resize(n_original_q_points);
+ contravariant.resize(n_q_points);
if (this->update_each & update_volume_elements)
- volume_elements.resize(n_original_q_points);
+ volume_elements.resize(n_q_points);
// see if we need the (transformation) shape function values
// and/or gradients and resize the necessary arrays
const Quadrature<dim> &q,
const unsigned int n_original_q_points)
{
- initialize(update_flags, q, n_original_q_points);
+ reinit(update_flags, q);
if (this->update_each &
(update_boundary_forms | update_normal_vectors | update_jacobians |
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(*this->fe);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
- data.initialize(this->requires_update_flags(update_flags), q, q.size());
-
+ data_ptr->reinit(this->requires_update_flags(update_flags), q);
return data_ptr;
}
MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
const FiniteElement<dim, spacedim> &fe,
const ComponentMask &mask)
- : unit_tangentials()
+ : fe(&fe)
+ , unit_tangentials()
, n_shape_functions(fe.n_dofs_per_cell())
, mask(mask)
, local_dof_indices(fe.n_dofs_per_cell())
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::InternalData::reinit(
+ const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature)
+{
+ // store the flags in the internal data object so we can access them
+ // in fill_fe_*_values(). use the transitive hull of the required
+ // flags
+ this->update_each = update_flags;
+
+ const unsigned int n_q_points = quadrature.size();
+ const std::vector<Point<dim>> &points = quadrature.get_points();
+
+ // see if we need the (transformation) shape function values
+ // and/or gradients and resize the necessary arrays
+ if (update_flags & update_quadrature_points)
+ {
+ shape_values.resize(n_shape_functions * n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < n_shape_functions; ++i)
+ shape(point, i) = fe->shape_value(i, points[point]);
+ }
+
+ if (update_flags &
+ (update_covariant_transformation | update_contravariant_transformation |
+ update_JxW_values | update_boundary_forms | update_normal_vectors |
+ update_jacobians | update_jacobian_grads | update_inverse_jacobians))
+ {
+ shape_derivatives.resize(n_shape_functions * n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < n_shape_functions; ++i)
+ derivative(point, i) = fe->shape_grad(i, points[point]);
+ }
+
+ if (update_flags & update_covariant_transformation)
+ covariant.resize(n_q_points);
+
+ if (update_flags & update_contravariant_transformation)
+ contravariant.resize(n_q_points);
+
+ if (update_flags & update_volume_elements)
+ volume_elements.resize(n_q_points);
+
+ if (update_flags &
+ (update_jacobian_grads | update_jacobian_pushed_forward_grads))
+ {
+ shape_second_derivatives.resize(n_shape_functions * n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < n_shape_functions; ++i)
+ second_derivative(point, i) = fe->shape_grad_grad(i, points[point]);
+ }
+
+ if (update_flags & (update_jacobian_2nd_derivatives |
+ update_jacobian_pushed_forward_2nd_derivatives))
+ {
+ shape_third_derivatives.resize(n_shape_functions * n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < n_shape_functions; ++i)
+ third_derivative(point, i) =
+ fe->shape_3rd_derivative(i, points[point]);
+ }
+
+ if (update_flags & (update_jacobian_3rd_derivatives |
+ update_jacobian_pushed_forward_3rd_derivatives))
+ {
+ shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
+ for (unsigned int point = 0; point < n_q_points; ++point)
+ for (unsigned int i = 0; i < n_shape_functions; ++i)
+ fourth_derivative(point, i) =
+ fe->shape_4th_derivative(i, points[point]);
+ }
+
+ // This (for face values and simplices) can be different for different
+ // calls, so always copy
+ quadrature_weights = quadrature.get_weights();
+}
+
+
+
template <int dim, int spacedim, typename VectorType>
std::size_t
MappingFEField<dim, spacedim, VectorType>::InternalData::memory_consumption()
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::compute_shapes_virtual(
- const std::vector<Point<dim>> &unit_points,
- typename MappingFEField<dim, spacedim, VectorType>::InternalData &data) const
-{
- const auto fe = &euler_dof_handler->get_fe();
- const unsigned int n_points = unit_points.size();
-
- for (unsigned int point = 0; point < n_points; ++point)
- {
- if (data.shape_values.size() != 0)
- for (unsigned int i = 0; i < data.n_shape_functions; ++i)
- data.shape(point, i) = fe->shape_value(i, unit_points[point]);
-
- if (data.shape_derivatives.size() != 0)
- for (unsigned int i = 0; i < data.n_shape_functions; ++i)
- data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
-
- if (data.shape_second_derivatives.size() != 0)
- for (unsigned int i = 0; i < data.n_shape_functions; ++i)
- data.second_derivative(point, i) =
- fe->shape_grad_grad(i, unit_points[point]);
-
- if (data.shape_third_derivatives.size() != 0)
- for (unsigned int i = 0; i < data.n_shape_functions; ++i)
- data.third_derivative(point, i) =
- fe->shape_3rd_derivative(i, unit_points[point]);
-
- if (data.shape_fourth_derivatives.size() != 0)
- for (unsigned int i = 0; i < data.n_shape_functions; ++i)
- data.fourth_derivative(point, i) =
- fe->shape_4th_derivative(i, unit_points[point]);
- }
-}
-
-
template <int dim, int spacedim, typename VectorType>
UpdateFlags
MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
}
-
template <int dim, int spacedim, typename VectorType>
void
-MappingFEField<dim, spacedim, VectorType>::compute_data(
- const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points,
- InternalData &data) const
+MappingFEField<dim, spacedim, VectorType>::compute_face_data(
+ const unsigned int n_original_q_points,
+ InternalData &data) const
{
- // store the flags in the internal data object so we can access them
- // in fill_fe_*_values(). use the transitive hull of the required
- // flags
- data.update_each = requires_update_flags(update_flags);
-
- const unsigned int n_q_points = q.size();
-
- // see if we need the (transformation) shape function values
- // and/or gradients and resize the necessary arrays
- if (data.update_each & update_quadrature_points)
- data.shape_values.resize(data.n_shape_functions * n_q_points);
-
- if (data.update_each &
- (update_covariant_transformation | update_contravariant_transformation |
- update_JxW_values | update_boundary_forms | update_normal_vectors |
- update_jacobians | update_jacobian_grads | update_inverse_jacobians))
- data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
-
+ // Set to the size of a single quadrature object for faces, as the size set
+ // in in reinit() is for all points
if (data.update_each & update_covariant_transformation)
data.covariant.resize(n_original_q_points);
if (data.update_each & update_volume_elements)
data.volume_elements.resize(n_original_q_points);
- if (data.update_each &
- (update_jacobian_grads | update_jacobian_pushed_forward_grads))
- data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
-
- if (data.update_each & (update_jacobian_2nd_derivatives |
- update_jacobian_pushed_forward_2nd_derivatives))
- data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
-
- if (data.update_each & (update_jacobian_3rd_derivatives |
- update_jacobian_pushed_forward_3rd_derivatives))
- data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
-
- compute_shapes_virtual(q.get_points(), data);
-
- // This (for face values and simplices) can be different for different calls,
- // so always copy
- data.quadrature_weights = q.get_weights();
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::compute_face_data(
- const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points,
- InternalData &data) const
-{
- compute_data(update_flags, q, n_original_q_points, data);
-
if (dim > 1)
{
if (data.update_each & update_boundary_forms)
}
+
template <int dim, int spacedim, typename VectorType>
typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
MappingFEField<dim, spacedim, VectorType>::get_data(
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
- this->compute_data(update_flags, quadrature, quadrature.size(), data);
+ data_ptr->reinit(requires_update_flags(update_flags), quadrature);
return data_ptr;
}
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
+ auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
const Quadrature<dim> q(
QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
- this->compute_face_data(update_flags, q, quadrature[0].size(), data);
+ data.reinit(requires_update_flags(update_flags), q);
+ this->compute_face_data(quadrature[0].size(), data);
return data_ptr;
}
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
+ auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
const Quadrature<dim> q(
QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
- this->compute_face_data(update_flags, q, quadrature.size(), data);
+ data.reinit(requires_update_flags(update_flags), q);
+ this->compute_face_data(quadrature.size(), data);
return data_ptr;
}
Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
ExcInternalError());
- update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
+ update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
- return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
+ return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
}
initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
- const Quadrature<dim> point_quadrature(initial_p_unit);
+ Quadrature<dim> point_quadrature(initial_p_unit);
UpdateFlags update_flags = update_quadrature_points | update_jacobians;
if (spacedim > dim)
Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
ExcInternalError());
- update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
+ update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
return do_transform_real_to_unit_cell(cell,
p,
- initial_p_unit,
- dynamic_cast<InternalData &>(*mdata));
+ point_quadrature,
+ static_cast<InternalData &>(*mdata));
}
MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
const typename Triangulation<dim, spacedim>::cell_iterator &cell,
const Point<spacedim> &p,
- const Point<dim> &initial_p_unit,
+ Quadrature<dim> &point_quadrature,
InternalData &mdata) const
{
const unsigned int n_shapes = mdata.shape_values.size();
// The shape values and derivatives
// of the mapping at this point are
// previously computed.
- // f(x)
- Point<dim> p_unit = initial_p_unit;
+
+ AssertDimension(point_quadrature.size(), 1);
+ Point<dim> p_unit = point_quadrature.point(0);
Point<dim> f;
- compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
+ mdata.reinit(mdata.update_each, point_quadrature);
+
Point<spacedim> p_real(do_transform_unit_to_real_cell(mdata));
Tensor<1, spacedim> p_minus_F = p - p_real;
const double eps = 1.e-12 * cell->diameter();
p_unit_trial[i] -= step_length * delta[i];
// shape values and derivatives
// at new p_unit point
- compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
- mdata);
+ point_quadrature.initialize(
+ ArrayView<const Point<dim>>(&p_unit_trial, 1));
+ mdata.reinit(mdata.update_each, point_quadrature);
// f(x)
Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
const Tensor<1, spacedim> f_trial = p - p_real_trial;
}
+
template <int dim, int spacedim>
void
-MappingManifold<dim, spacedim>::InternalData::initialize(
+MappingManifold<dim, spacedim>::InternalData::reinit(
const UpdateFlags update_flags,
- const Quadrature<dim> &q,
- const unsigned int n_original_q_points)
+ const Quadrature<dim> &q)
{
// store the flags in the internal data object so we can access them
// in fill_fe_*_values()
this->update_each = update_flags;
+ const unsigned int n_q_points = q.size();
+
// Store the quadrature
- this->quad = q;
+ this->quad.initialize(q.get_points(), q.get_weights());
// Resize the weights
this->vertex_weights.resize(GeometryInfo<dim>::vertices_per_cell);
compute_manifold_quadrature_weights(q);
if (this->update_each & update_covariant_transformation)
- covariant.resize(n_original_q_points);
+ covariant.resize(n_q_points);
if (this->update_each & update_contravariant_transformation)
- contravariant.resize(n_original_q_points);
+ contravariant.resize(n_q_points);
if (this->update_each & update_volume_elements)
- volume_elements.resize(n_original_q_points);
+ volume_elements.resize(n_q_points);
}
+
template <int dim, int spacedim>
void
MappingManifold<dim, spacedim>::InternalData::initialize_face(
const Quadrature<dim> &q,
const unsigned int n_original_q_points)
{
- initialize(update_flags, q, n_original_q_points);
+ reinit(update_flags, q);
+
+ // Set to the size of a single quadrature object for faces, as the size set
+ // in in reinit() is for all points
+ if (this->update_each & update_covariant_transformation)
+ covariant.resize(n_original_q_points);
+
+ if (this->update_each & update_contravariant_transformation)
+ contravariant.resize(n_original_q_points);
+
+ if (this->update_each & update_volume_elements)
+ volume_elements.resize(n_original_q_points);
if (dim > 1)
{
+template <int dim, int spacedim>
+void
+MappingManifold<dim, spacedim>::InternalData::store_vertices(
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
+{
+ vertices.resize(GeometryInfo<dim>::vertices_per_cell);
+ for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+ vertices[i] = cell->vertex(i);
+ this->cell = cell;
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingManifold<dim, spacedim>::InternalData::
+ compute_manifold_quadrature_weights(const Quadrature<dim> &quad)
+{
+ cell_manifold_quadrature_weights.resize(
+ quad.size(), std::vector<double>(GeometryInfo<dim>::vertices_per_cell));
+ for (unsigned int q = 0; q < quad.size(); ++q)
+ {
+ for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+ {
+ cell_manifold_quadrature_weights[q][i] =
+ GeometryInfo<dim>::d_linear_shape_function(quad.point(q), i);
+ }
+ }
+}
+
+
+
template <int dim, int spacedim>
MappingManifold<dim, spacedim>::MappingManifold(
const MappingManifold<dim, spacedim> &)
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>();
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
- data.initialize(this->requires_update_flags(update_flags), q, q.size());
+ data_ptr->reinit(this->requires_update_flags(update_flags), q);
return data_ptr;
}
template <int dim, int spacedim>
void
-MappingQ<dim, spacedim>::InternalData::initialize(
- const UpdateFlags update_flags,
- const Quadrature<dim> &quadrature,
- const unsigned int n_original_q_points)
+MappingQ<dim, spacedim>::InternalData::reinit(const UpdateFlags update_flags,
+ const Quadrature<dim> &quadrature)
{
// store the flags in the internal data object so we can access them
// in fill_fe_*_values()
const unsigned int n_q_points = quadrature.size();
if (this->update_each & update_volume_elements)
- volume_elements.resize(n_original_q_points);
+ volume_elements.resize(n_q_points);
tensor_product_quadrature = quadrature.is_tensor_product();
const Quadrature<dim> &quadrature,
const unsigned int n_original_q_points)
{
- initialize(update_flags, quadrature, n_original_q_points);
+ reinit(update_flags, quadrature);
quadrature_points = quadrature.get_points();
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(polynomial_degree);
- auto &data = dynamic_cast<InternalData &>(*data_ptr);
- data.initialize(this->requires_update_flags(update_flags), q, q.size());
-
+ data_ptr->reinit(update_flags, q);
return data_ptr;
}
projects / dealii.git / commitdiff