* MatrixFree object was given at initialization.
*/
mutable std::vector<types::global_dof_index> local_dof_indices;
-
- /**
- * A temporary data structure for the Jacobian information necessary if
- * the reinit() function is used that allows to construct user batches.
- */
- AlignedVector<Tensor<2, dim, VectorizedArrayType>> jacobian_data;
-
- /**
- * A temporary data structure for the Jacobian determinant necessary if
- * the reinit() function is used that allows to construct user batches.
- */
- AlignedVector<VectorizedArrayType> J_value_data;
-
- /**
- * A temporary data structure for the gradients of the inverse Jacobian
- * transformation necessary if the reinit() function is used that allows
- * to construct user batches.
- */
- AlignedVector<
- Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>>
- jacobian_gradients_data;
-
- /**
- * A temporary data structure for the quadrature-point locations necessary if
- * the reinit() function is used that allows to construct user batches.
- */
- AlignedVector<Point<dim, VectorizedArrayType>> quadrature_points_data;
};
VectorizedArrayType>::
reinit(const std::array<unsigned int, VectorizedArrayType::size()> &cell_ids)
{
- Assert(this->mapped_geometry == nullptr,
- ExcMessage("FEEvaluation was initialized without a matrix-free object."
- " Integer indexing is not possible"));
- if (this->mapped_geometry != nullptr)
- return;
-
Assert(this->dof_info != nullptr, ExcNotInitialized());
Assert(this->mapping_data != nullptr, ExcNotInitialized());
}
// allocate memory for internal data storage
+ if (this->mapped_geometry == nullptr)
+ this->mapped_geometry =
+ std::make_shared<internal::MatrixFreeFunctions::
+ MappingDataOnTheFly<dim, VectorizedArrayType>>();
+
+ auto &mapping_storage = this->mapped_geometry->get_data_storage();
+
+ auto &this_jacobian_data = mapping_storage.jacobians[0];
+ auto &this_J_value_data = mapping_storage.JxW_values;
+ auto &this_jacobian_gradients_data = mapping_storage.jacobian_gradients[0];
+ auto &this_quadrature_points_data = mapping_storage.quadrature_points;
if (this->cell_type <= internal::MatrixFreeFunctions::GeometryType::affine)
{
if (this->mapping_data->jacobians[0].size() > 0)
- this->jacobian_data.resize_fast(2);
+ this_jacobian_data.resize_fast(2);
if (this->mapping_data->JxW_values.size() > 0)
- this->J_value_data.resize_fast(1);
+ this_J_value_data.resize_fast(1);
if (this->mapping_data->jacobian_gradients[0].size() > 0)
- this->jacobian_gradients_data.resize_fast(1);
+ this_jacobian_gradients_data.resize_fast(1);
if (this->mapping_data->quadrature_points.size() > 0)
- this->quadrature_points_data.resize_fast(1);
+ this_quadrature_points_data.resize_fast(1);
}
else
{
if (this->mapping_data->jacobians[0].size() > 0)
- this->jacobian_data.resize_fast(this->n_quadrature_points);
+ this_jacobian_data.resize_fast(this->n_quadrature_points);
if (this->mapping_data->JxW_values.size() > 0)
- this->J_value_data.resize_fast(this->n_quadrature_points);
+ this_J_value_data.resize_fast(this->n_quadrature_points);
if (this->mapping_data->jacobian_gradients[0].size() > 0)
- this->jacobian_gradients_data.resize_fast(this->n_quadrature_points);
+ this_jacobian_gradients_data.resize_fast(this->n_quadrature_points);
if (this->mapping_data->quadrature_points.size() > 0)
- this->quadrature_points_data.resize_fast(this->n_quadrature_points);
+ this_quadrature_points_data.resize_fast(this->n_quadrature_points);
}
// set pointers to internal data storage
- this->jacobian = this->jacobian_data.data();
- this->J_value = this->J_value_data.data();
- this->jacobian_gradients = this->jacobian_gradients_data.data();
- this->quadrature_points = this->quadrature_points_data.data();
+ this->jacobian = this_jacobian_data.data();
+ this->J_value = this_J_value_data.data();
+ this->jacobian_gradients = this_jacobian_gradients_data.data();
+ this->quadrature_points = this_quadrature_points_data.data();
// fill internal data storage lane by lane
for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
const unsigned int q = 0;
if (this->mapping_data->JxW_values.size() > 0)
- this->J_value_data[q][v] =
+ this_J_value_data[q][v] =
this->mapping_data->JxW_values[offsets + q][lane];
if (this->mapping_data->jacobians[0].size() > 0)
for (unsigned int q = 0; q < 2; ++q)
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = 0; j < dim; ++j)
- this->jacobian_data[q][i][j][v] =
+ this_jacobian_data[q][i][j][v] =
this->mapping_data->jacobians[0][offsets + q][i][j][lane];
if (this->mapping_data->jacobian_gradients[0].size() > 0)
for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
for (unsigned int j = 0; j < dim; ++j)
- this->jacobian_gradients_data[q][i][j][v] =
+ this_jacobian_gradients_data[q][i][j][v] =
this->mapping_data
->jacobian_gradients[0][offsets + q][i][j][lane];
if (this->mapping_data->quadrature_points.size() > 0)
for (unsigned int i = 0; i < dim; ++i)
- this->quadrature_points_data[q][i][v] =
+ this_quadrature_points_data[q][i][v] =
this->mapping_data->quadrature_points
[this->mapping_data
->quadrature_point_offsets[cell_batch_index] +
q;
if (this->mapping_data->JxW_values.size() > 0)
- this->J_value_data[q][v] =
+ this_J_value_data[q][v] =
this->mapping_data->JxW_values[offsets + q_src][lane];
if (this->mapping_data->jacobians[0].size() > 0)
for (unsigned int i = 0; i < dim; ++i)
for (unsigned int j = 0; j < dim; ++j)
- this->jacobian_data[q][i][j][v] =
+ this_jacobian_data[q][i][j][v] =
this->mapping_data
->jacobians[0][offsets + q_src][i][j][lane];
if (this->mapping_data->jacobian_gradients[0].size() > 0)
for (unsigned int i = 0; i < dim * (dim + 1) / 2; ++i)
for (unsigned int j = 0; j < dim; ++j)
- this->jacobian_gradients_data[q][i][j][v] =
+ this_jacobian_gradients_data[q][i][j][v] =
this->mapping_data
->jacobian_gradients[0][offsets + q_src][i][j][lane];
this->descriptor->quadrature.point(q)[e]);
for (unsigned int i = 0; i < dim; ++i)
- this->quadrature_points_data[q][i][v] = point[i][lane];
+ this_quadrature_points_data[q][i][v] = point[i][lane];
}
else
{
// general case: quadrature points are available
for (unsigned int i = 0; i < dim; ++i)
- this->quadrature_points_data[q][i][v] =
+ this_quadrature_points_data[q][i][v] =
this->mapping_data->quadrature_points
[this->mapping_data
->quadrature_point_offsets[cell_batch_index] +
class MappingDataOnTheFly
{
public:
+ /**
+ * Default constructor.
+ */
+ MappingDataOnTheFly() = default;
+
/**
* Constructor, similar to FEValues. Since this class only evaluates the
* geometry, no finite element has to be specified and the simplest
const MappingInfoStorage<dim, dim, Number> &
get_data_storage() const;
+ /**
+ * Return a non-const reference to the underlying storage field
+ * of type MappingInfoStorage of the same format as the data fields
+ * in MappingInfo. This function can be used to manually fill the
+ * data if the default constructor has been called and reinit() was
+ * not called for a given cell.
+ */
+ MappingInfoStorage<dim, dim, Number> &
+ get_data_storage();
+
/**
* Return a reference to 1D quadrature underlying this object.
*/
/**
* An underlying FEValues object that performs the (scalar) evaluation.
*/
- dealii::FEValues<dim> fe_values;
+ std::unique_ptr<dealii::FEValues<dim>> fe_values;
/**
* Get 1D quadrature formula to be used for reinitializing shape info.
const Mapping<dim> & mapping,
const Quadrature<1> &quadrature,
const UpdateFlags update_flags)
- : fe_values(mapping,
- fe_dummy,
- Quadrature<dim>(quadrature),
- MappingInfoStorage<dim, dim, Number>::compute_update_flags(
- update_flags))
+ : fe_values(std::make_unique<dealii::FEValues<dim>>(
+ mapping,
+ fe_dummy,
+ Quadrature<dim>(quadrature),
+ MappingInfoStorage<dim, dim, Number>::compute_update_flags(
+ update_flags)))
, quadrature_1d(quadrature)
{
mapping_info_storage.descriptor.resize(1);
mapping_info_storage.descriptor[0].initialize(quadrature);
mapping_info_storage.data_index_offsets.resize(1);
- mapping_info_storage.JxW_values.resize(fe_values.n_quadrature_points);
- mapping_info_storage.jacobians[0].resize(fe_values.n_quadrature_points);
+ mapping_info_storage.JxW_values.resize(fe_values->n_quadrature_points);
+ mapping_info_storage.jacobians[0].resize(fe_values->n_quadrature_points);
if ((update_flags & update_quadrature_points) != 0u)
{
mapping_info_storage.quadrature_point_offsets.resize(1, 0);
mapping_info_storage.quadrature_points.resize(
- fe_values.n_quadrature_points);
+ fe_values->n_quadrature_points);
}
- if (fe_values.get_update_flags() & update_normal_vectors)
+ if (fe_values->get_update_flags() & update_normal_vectors)
{
mapping_info_storage.normal_vectors.resize(
- fe_values.n_quadrature_points);
+ fe_values->n_quadrature_points);
mapping_info_storage.normals_times_jacobians[0].resize(
- fe_values.n_quadrature_points);
+ fe_values->n_quadrature_points);
}
- Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
+ Assert(!(fe_values->get_update_flags() & update_jacobian_grads),
ExcNotImplemented());
}
if (present_cell == cell)
return;
present_cell = cell;
- fe_values.reinit(present_cell);
- for (unsigned int q = 0; q < fe_values.get_quadrature().size(); ++q)
+ fe_values->reinit(present_cell);
+ for (unsigned int q = 0; q < fe_values->get_quadrature().size(); ++q)
{
- if (fe_values.get_update_flags() & update_JxW_values)
- mapping_info_storage.JxW_values[q] = fe_values.JxW(q);
- if (fe_values.get_update_flags() & update_jacobians)
+ if (fe_values->get_update_flags() & update_JxW_values)
+ mapping_info_storage.JxW_values[q] = fe_values->JxW(q);
+ if (fe_values->get_update_flags() & update_jacobians)
{
- Tensor<2, dim> jac = fe_values.jacobian(q);
+ Tensor<2, dim> jac = fe_values->jacobian(q);
jac = invert(transpose(jac));
for (unsigned int d = 0; d < dim; ++d)
for (unsigned int e = 0; e < dim; ++e)
mapping_info_storage.jacobians[0][q][d][e] = jac[d][e];
}
- if (fe_values.get_update_flags() & update_quadrature_points)
+ if (fe_values->get_update_flags() & update_quadrature_points)
for (unsigned int d = 0; d < dim; ++d)
mapping_info_storage.quadrature_points[q][d] =
- fe_values.quadrature_point(q)[d];
- if (fe_values.get_update_flags() & update_normal_vectors)
+ fe_values->quadrature_point(q)[d];
+ if (fe_values->get_update_flags() & update_normal_vectors)
{
for (unsigned int d = 0; d < dim; ++d)
mapping_info_storage.normal_vectors[q][d] =
- fe_values.normal_vector(q)[d];
+ fe_values->normal_vector(q)[d];
mapping_info_storage.normals_times_jacobians[0][q] =
mapping_info_storage.normal_vectors[q] *
mapping_info_storage.jacobians[0][q];
inline typename dealii::Triangulation<dim>::cell_iterator
MappingDataOnTheFly<dim, Number>::get_cell() const
{
- return fe_values.get_cell();
+ return fe_values->get_cell();
}
inline const dealii::FEValues<dim> &
MappingDataOnTheFly<dim, Number>::get_fe_values() const
{
- return fe_values;
+ return *fe_values;
}
+ template <int dim, typename Number>
+ inline MappingInfoStorage<dim, dim, Number> &
+ MappingDataOnTheFly<dim, Number>::get_data_storage()
+ {
+ return mapping_info_storage;
+ }
+
+
+
template <int dim, typename Number>
inline const Quadrature<1> &
MappingDataOnTheFly<dim, Number>::get_quadrature() const
projects / dealii.git / commitdiff