#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/helper_functions.h>
+#include <deal.II/matrix_free/mapping_info.h>
DEAL_II_NAMESPACE_OPEN
void reinit (const typename Triangulation<dim>::cell_iterator &cell);
/**
- * For the transformation information stored in MappingInfo, this function
- * returns the index which belongs to the current cell as specified in @p
- * reinit. Note that MappingInfo has different fields for Cartesian cells,
- * cells with affine mapping and with general mappings, so in order to
- * access the correct data, this interface must be used together with
- * get_cell_type.
+ * @deprecated Use get_mapping_data_index_offset() instead.
*/
+ DEAL_II_DEPRECATED
unsigned int get_cell_data_number() const;
/**
- * Return the type of the cell the @p reinit function has been called for.
+ * 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.
*/
VectorizedArray<Number> JxW(const unsigned int q_point) const;
+ /**
+ * Gets the inverse and transposed version of Jacobian of the mapping
+ * between the unit to the real cell (representing the covariant
+ * transformation). This is exactly the matrix used internally to transform
+ * the unit cell gradients to gradients on the real cell.
+ */
+ Tensor<2,dim,VectorizedArray<Number> >
+ inverse_jacobian(const unsigned int q_index) const;
+
//@}
/**
const internal::MatrixFreeFunctions::DoFInfo *dof_info;
/**
- * Stores a pointer to the underlying transformation data from unit to real
- * cells for the given quadrature formula specified at construction. Also
- * contained in matrix_info, but it simplifies code if we store a reference
- * to it.
+ * Stores a pointer to the underlying transformation data from unit to
+ * real cells for the given quadrature formula specified at construction.
+ * Also contained in matrix_info, but it simplifies code if we store a
+ * reference to it.
*/
- const internal::MatrixFreeFunctions::MappingInfo<dim,Number> *mapping_info;
+ const internal::MatrixFreeFunctions::MappingInfoStorage<dim,dim,Number> *mapping_data;
/**
* Stores a pointer to the unit cell shape data, i.e., values, gradients and
*/
const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> *data;
- /**
- * A pointer to the Cartesian Jacobian information of the present cell. Only
- * set to a useful value if on a Cartesian cell, otherwise zero.
- */
- const Tensor<1,dim,VectorizedArray<Number> > *cartesian_data;
-
/**
* A pointer to the Jacobian information of the present cell. Only set to a
* useful value if on a non-Cartesian cell.
/**
* A pointer to the quadrature weights of the underlying quadrature formula.
*/
- const VectorizedArray<Number> *quadrature_weights;
-
- /**
- * A pointer to the quadrature points on the present cell.
- */
- const Point<dim,VectorizedArray<Number> > *quadrature_points;
-
- /**
- * A pointer to the diagonal part of the Jacobian gradient on the present
- * cell. Only set to a useful value if on a general cell with non-constant
- * Jacobian.
- */
- const Tensor<2,dim,VectorizedArray<Number> > *jacobian_grad;
-
- /**
- * A pointer to the upper diagonal part of the Jacobian gradient on the
- * present cell. Only set to a useful value if on a general cell with non-
- * constant Jacobian.
- */
- const Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > * jacobian_grad_upper;
+ const Number *quadrature_weights;
/**
* After a call to reinit(), stores the number of the cell we are currently
*/
internal::MatrixFreeFunctions::CellType cell_type;
- /**
- * The stride to access the correct data in MappingInfo.
- */
- unsigned int cell_data_number;
-
/**
* Debug information to track whether dof values have been initialized
* before accessed. Used to control exceptions when uninitialized data is
*/
std::shared_ptr<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> > mapped_geometry;
- /**
- * For use with on-the-fly evaluation, provide a data structure to store the
- * global dof indices on the current cell from a reinit call.
- */
- std::vector<types::global_dof_index> old_style_dof_indices;
-
/**
* For a FiniteElement with more than one finite element, select at which
* component this data structure should start.
:
0),
active_quad_index (fe_degree != numbers::invalid_unsigned_int ?
- data_in.get_mapping_info().
- mapping_data_gen[quad_no_in].
- quad_index_from_n_q_points(n_q_points)
+ (data_in.get_mapping_info().cell_data[quad_no_in].
+ quad_index_from_n_q_points(n_q_points))
:
0),
matrix_info (&data_in),
dof_info (&data_in.get_dof_info(fe_no_in)),
- mapping_info (&data_in.get_mapping_info()),
+ mapping_data (internal::MatrixFreeFunctions::MappingInfoCellsOrFaces<dim,Number,false>::get
+ (data_in.get_mapping_info(), quad_no)),
data (&data_in.get_shape_info
(fe_no_in, quad_no_in, active_fe_index,
active_quad_index)),
- cartesian_data (nullptr),
jacobian (nullptr),
J_value (nullptr),
- quadrature_weights (mapping_info->mapping_data_gen[quad_no].
- quadrature_weights[active_quad_index].begin()),
- quadrature_points (nullptr),
- jacobian_grad (nullptr),
- jacobian_grad_upper(nullptr),
+ quadrature_weights (mapping_data->descriptor[active_quad_index].quadrature_weights.begin()),
cell (numbers::invalid_unsigned_int),
- cell_type (internal::MatrixFreeFunctions::undefined),
- cell_data_number (numbers::invalid_unsigned_int),
+ cell_type (internal::MatrixFreeFunctions::general),
dof_values_initialized (false),
values_quad_initialized (false),
gradients_quad_initialized(false),
AssertDimension (data->dofs_per_component_on_cell*n_fe_components,
dof_info->dofs_per_cell[active_fe_index]);
AssertDimension (data->n_q_points,
- mapping_info->mapping_data_gen[quad_no].n_q_points[active_quad_index]);
+ mapping_data->descriptor[active_quad_index].n_q_points);
Assert (n_fe_components == 1 ||
n_components == 1 ||
n_components == n_fe_components,
active_quad_index (numbers::invalid_unsigned_int),
matrix_info (nullptr),
dof_info (nullptr),
- mapping_info (nullptr),
+ mapping_data (nullptr),
// select the correct base element from the given FE component
data (new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(quadrature, fe, fe.component_to_base_index(first_selected_component).first)),
- cartesian_data (nullptr),
jacobian (nullptr),
J_value (nullptr),
quadrature_weights (nullptr),
- quadrature_points (nullptr),
- jacobian_grad (nullptr),
- jacobian_grad_upper(nullptr),
cell (0),
cell_type (internal::MatrixFreeFunctions::general),
- cell_data_number (numbers::invalid_unsigned_int),
dof_values_initialized (false),
values_quad_initialized (false),
gradients_quad_initialized(false),
mapped_geometry
= std::make_shared<internal::MatrixFreeFunctions::MappingDataOnTheFly<dim,Number> >
(mapping, quadrature, update_flags);
- jacobian = mapped_geometry->get_inverse_jacobians().begin();
- J_value = mapped_geometry->get_JxW_values().begin();
- quadrature_points = mapped_geometry->get_quadrature_points().begin();
+
+ 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();
Assert(fe.element_multiplicity(base_element_number) == 1 ||
fe.element_multiplicity(base_element_number)-first_selected_component >= n_components_,
active_quad_index (other.active_quad_index),
matrix_info (other.matrix_info),
dof_info (other.dof_info),
- mapping_info (other.mapping_info),
+ mapping_data (other.mapping_data),
data (other.matrix_info == nullptr ?
new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data) :
other.data),
- cartesian_data (nullptr),
jacobian (nullptr),
J_value (nullptr),
- quadrature_weights (mapping_info != nullptr ?
- mapping_info->mapping_data_gen[quad_no].
- quadrature_weights[active_quad_index].begin()
- :
- nullptr),
- quadrature_points (nullptr),
- jacobian_grad (nullptr),
- jacobian_grad_upper(nullptr),
+ quadrature_weights (other.matrix_info == nullptr ? nullptr :
+ mapping_data->descriptor[active_quad_index].quadrature_weights.begin()),
cell (numbers::invalid_unsigned_int),
cell_type (internal::MatrixFreeFunctions::general),
- cell_data_number (numbers::invalid_unsigned_int),
dof_values_initialized (false),
values_quad_initialized (false),
gradients_quad_initialized(false),
MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
other.mapped_geometry->get_quadrature(),
other.mapped_geometry->get_fe_values().get_update_flags()));
- jacobian = mapped_geometry->get_inverse_jacobians().begin();
- J_value = mapped_geometry->get_JxW_values().begin();
- quadrature_points = mapped_geometry->get_quadrature_points().begin();
+ mapping_data = &mapped_geometry->get_data_storage();
cell = 0;
+
+ jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
+ J_value = mapped_geometry->get_data_storage().JxW_values.begin();
}
}
matrix_info = other.matrix_info;
dof_info = other.dof_info;
- mapping_info = other.mapping_info;
+ mapping_data = other.mapping_data;
if (other.matrix_info == nullptr)
{
data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>>(*other.data);
}
set_data_pointers();
- cartesian_data = nullptr;
jacobian = nullptr;
J_value = nullptr;
- quadrature_weights = mapping_info != nullptr ?
- mapping_info->mapping_data_gen[quad_no].
- quadrature_weights[active_quad_index].begin()
- :
- nullptr;
- quadrature_points = nullptr;
- jacobian_grad = nullptr;
- jacobian_grad_upper = nullptr;
+ quadrature_weights = (mapping_data != nullptr ?
+ mapping_data->descriptor[active_quad_index].quadrature_weights.begin()
+ :
+ nullptr);
cell = numbers::invalid_unsigned_int;
cell_type = internal::MatrixFreeFunctions::general;
- cell_data_number = numbers::invalid_unsigned_int;
// Create deep copy of mapped geometry for use in parallel...
if (other.mapped_geometry.get() != nullptr)
MappingDataOnTheFly<dim,Number>(other.mapped_geometry->get_fe_values().get_mapping(),
other.mapped_geometry->get_quadrature(),
other.mapped_geometry->get_fe_values().get_update_flags()));
- jacobian = mapped_geometry->get_inverse_jacobians().begin();
- J_value = mapped_geometry->get_JxW_values().begin();
- quadrature_points = mapped_geometry->get_quadrature_points().begin();
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, int n_components_, typename Number>
inline
void
-FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_in)
+FEEvaluationBase<dim,n_components_,Number>::reinit (const unsigned int cell_index)
{
Assert (mapped_geometry == nullptr,
ExcMessage("FEEvaluation was initialized without a matrix-free object."
" Integer indexing is not possible"));
if (mapped_geometry != nullptr)
return;
- Assert (dof_info != nullptr, ExcNotInitialized());
- Assert (mapping_info != nullptr, ExcNotInitialized());
- AssertIndexRange (cell_in, dof_info->row_starts.size()-1);
- AssertDimension (((dof_info->cell_active_fe_index.size() > 0) ?
- dof_info->cell_active_fe_index[cell_in] : 0),
- active_fe_index);
- cell = cell_in;
- cell_type = mapping_info->get_cell_type(cell);
- cell_data_number = mapping_info->get_cell_data_index(cell);
-
- if (mapping_info->quadrature_points_initialized == true)
- {
- AssertIndexRange (cell_data_number, mapping_info->
- mapping_data_gen[quad_no].rowstart_q_points.size());
- const unsigned int index = mapping_info->mapping_data_gen[quad_no].
- rowstart_q_points[cell];
- AssertIndexRange (index, mapping_info->mapping_data_gen[quad_no].
- quadrature_points.size());
- quadrature_points =
- &mapping_info->mapping_data_gen[quad_no].quadrature_points[index];
- }
- if (cell_type == internal::MatrixFreeFunctions::cartesian)
- {
- cartesian_data = &mapping_info->cartesian_data[cell_data_number].first;
- J_value = &mapping_info->cartesian_data[cell_data_number].second;
- }
- else if (cell_type == internal::MatrixFreeFunctions::affine)
- {
- jacobian = &mapping_info->affine_data[cell_data_number].first;
- J_value = &mapping_info->affine_data[cell_data_number].second;
- }
- else
- {
- const unsigned int rowstart = mapping_info->
- mapping_data_gen[quad_no].rowstart_jacobians[cell_data_number];
- AssertIndexRange (rowstart, mapping_info->
- mapping_data_gen[quad_no].jacobians.size());
- jacobian =
- &mapping_info->mapping_data_gen[quad_no].jacobians[rowstart];
- if (mapping_info->JxW_values_initialized == true)
- {
- AssertIndexRange (rowstart, mapping_info->
- mapping_data_gen[quad_no].JxW_values.size());
- J_value = &(mapping_info->mapping_data_gen[quad_no].
- JxW_values[rowstart]);
- }
- if (mapping_info->second_derivatives_initialized == true)
- {
- AssertIndexRange(rowstart, mapping_info->
- mapping_data_gen[quad_no].jacobians_grad_diag.size());
- jacobian_grad = &mapping_info->mapping_data_gen[quad_no].
- jacobians_grad_diag[rowstart];
- AssertIndexRange(rowstart, mapping_info->
- mapping_data_gen[quad_no].jacobians_grad_upper.size());
- jacobian_grad_upper = &mapping_info->mapping_data_gen[quad_no].
- jacobians_grad_upper[rowstart];
- }
- }
+ Assert (this->dof_info != nullptr, ExcNotInitialized());
+ Assert (this->mapping_data != nullptr, ExcNotInitialized());
+ this->cell = cell_index;
+ this->cell_type = this->matrix_info->get_mapping_info().get_cell_type(cell_index);
+
+ const unsigned int offsets = this->mapping_data->data_index_offsets[cell_index];
+ this->jacobian = &this->mapping_data->jacobians[0][offsets];
+ this->J_value = &this->mapping_data->JxW_values[offsets];
#ifdef DEBUG
dof_values_initialized = false;
FEEvaluationBase<dim,n_components_,Number>
::get_cell_data_number () const
{
- Assert (cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- return cell_data_number;
+ return get_mapping_data_index_offset();
+}
+
+
+
+template <int dim, int n_components_, typename Number>
+inline
+unsigned int
+FEEvaluationBase<dim,n_components_,Number>
+::get_mapping_data_index_offset () const
+{
+ if (matrix_info == 0)
+ return 0;
+ else
+ {
+ AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
+ return this->mapping_data->data_index_offsets[cell];
+ }
}
if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
this->cell_type == internal::MatrixFreeFunctions::affine)
{
- Assert (this->mapping_info != nullptr, ExcNotImplemented());
+ Assert (this->mapping_data != nullptr, ExcNotImplemented());
VectorizedArray<Number> J = this->J_value[0];
for (unsigned int q=0; q<this->data->n_q_points; ++q)
JxW_values[q] = J * this->quadrature_weights[q];
VectorizedArray<Number>
FEEvaluationBase<dim,n_components_,Number>::JxW(const unsigned int q_point) const
{
+ AssertIndexRange(q_point, data->n_q_points);
Assert (this->J_value != nullptr, ExcNotInitialized());
if (this->cell_type == internal::MatrixFreeFunctions::cartesian ||
this->cell_type == internal::MatrixFreeFunctions::affine)
{
- Assert (this->mapping_info != nullptr, ExcInternalError());
+ Assert (this->quadrature_weights != nullptr, ExcInternalError());
return this->J_value[0] * this->quadrature_weights[q_point];
}
else
+template <int dim, int n_components_, typename Number>
+inline
+Tensor<2,dim,VectorizedArray<Number> >
+FEEvaluationBase<dim,n_components_,Number>
+::inverse_jacobian(const unsigned int q_index) const
+{
+ AssertIndexRange(q_index, data->n_q_points);
+ Assert (this->jacobian != nullptr, ExcNotImplemented());
+ if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+ return jacobian[0];
+ else
+ return jacobian[q_index];
+}
+
+
+
namespace internal
{
// write access to generic vectors that have operator ().
for (unsigned int comp=0; comp<n_components; comp++)
for (unsigned int d=0; d<dim; ++d)
grad_out[comp][d] = (this->gradients_quad[comp][d][q_point] *
- cartesian_data[0][d]);
+ jacobian[0][d][d]);
}
// cell with general/affine Jacobian
else
{
const Tensor<2,dim,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
jacobian[q_point] : jacobian[0];
for (unsigned int comp=0; comp<n_components; comp++)
{
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, this->data->n_q_points);
+ Assert(jacobian != nullptr, ExcNotImplemented());
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point];
+
Tensor<2,dim,VectorizedArray<Number> > hessian_out [n_components];
// Cartesian cell
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
- const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
for (unsigned int comp=0; comp<n_components; comp++)
for (unsigned int d=0; d<dim; ++d)
{
hessian_out[comp][d][d] = (this->hessians_quad[comp][d][q_point] *
- jac[d] * jac[d]);
+ jac[d][d] * jac[d][d]);
switch (dim)
{
case 1:
break;
case 2:
hessian_out[comp][0][1] = (this->hessians_quad[comp][2][q_point] *
- jac[0] * jac[1]);
+ jac[0][0] * jac[1][1]);
break;
case 3:
hessian_out[comp][0][1] = (this->hessians_quad[comp][3][q_point] *
- jac[0] * jac[1]);
+ jac[0][0] * jac[1][1]);
hessian_out[comp][0][2] = (this->hessians_quad[comp][4][q_point] *
- jac[0] * jac[2]);
+ jac[0][0] * jac[2][2]);
hessian_out[comp][1][2] = (this->hessians_quad[comp][5][q_point] *
- jac[1] * jac[2]);
+ jac[1][1] * jac[2][2]);
break;
default:
Assert (false, ExcNotImplemented());
hessian_out[comp][e][d] = hessian_out[comp][d][e];
}
}
- // cell with general Jacobian
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
+ // cell with general Jacobian, but constant within the cell
+ else if (this->cell_type == internal::MatrixFreeFunctions::affine)
{
- Assert (this->mapping_info->second_derivatives_initialized == true,
- ExcNotInitialized());
- const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
- const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
- const Tensor<1,(dim>1?dim*(dim-1)/2:1),
- Tensor<1,dim,VectorizedArray<Number> > >
- & jac_grad_UT = jacobian_grad_upper[q_point];
for (unsigned int comp=0; comp<n_components; comp++)
{
// compute laplacian before the gradient because it needs to access
hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
}
- // add diagonal part of J' * grad(u)
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- hessian_out[comp][d][d] += (jac_grad[d][e] *
- this->gradients_quad[comp][e][q_point]);
-
- // add off-diagonal part of J' * grad(u)
- for (unsigned int d=0, count=0; d<dim; ++d)
- for (unsigned int e=d+1; e<dim; ++e, ++count)
- for (unsigned int f=0; f<dim; ++f)
- hessian_out[comp][d][e] += (jac_grad_UT[count][f] *
- this->gradients_quad[comp][f][q_point]);
+ // no J' * grad(u) part here because the Jacobian is constant
+ // throughout the cell and hence, its derivative is zero
// take symmetric part
for (unsigned int d=0; d<dim; ++d)
hessian_out[comp][e][d] = hessian_out[comp][d][e];
}
}
- // cell with general Jacobian, but constant within the cell
- else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // cell with general Jacobian
+ else
{
- const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
+ const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > &jac_grad =
+ mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point];
for (unsigned int comp=0; comp<n_components; comp++)
{
// compute laplacian before the gradient because it needs to access
hessian_out[comp][d][e] += jac[d][f] * tmp[f][e];
}
- // no J' * grad(u) part here because the Jacobian is constant
- // throughout the cell and hence, its derivative is zero
+ // add diagonal part of J' * grad(u)
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ hessian_out[comp][d][d] += (jac_grad[d][e] *
+ this->gradients_quad[comp][e][q_point]);
+
+ // add off-diagonal part of J' * grad(u)
+ for (unsigned int d=0, count=dim; d<dim; ++d)
+ for (unsigned int e=d+1; e<dim; ++e, ++count)
+ for (unsigned int f=0; f<dim; ++f)
+ hessian_out[comp][d][e] += (jac_grad[count][f] *
+ this->gradients_quad[comp][f][q_point]);
// take symmetric part
for (unsigned int d=0; d<dim; ++d)
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, this->data->n_q_points);
+ Assert(jacobian != nullptr, ExcNotImplemented());
+ const Tensor<2,dim,VectorizedArray<Number> > &jac =
+ jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ? 0 : q_point];
+
Tensor<1,n_components_,Tensor<1,dim,VectorizedArray<Number> > > hessian_out;
// Cartesian cell
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
- const Tensor<1,dim,VectorizedArray<Number> > &jac = cartesian_data[0];
for (unsigned int comp=0; comp<n_components; comp++)
for (unsigned int d=0; d<dim; ++d)
hessian_out[comp][d] = (this->hessians_quad[comp][d][q_point] *
- jac[d] * jac[d]);
+ jac[d][d] * jac[d][d]);
}
- // cell with general Jacobian
- else if (this->cell_type == internal::MatrixFreeFunctions::general)
+ // cell with general Jacobian, but constant within the cell
+ else if (this->cell_type == internal::MatrixFreeFunctions::affine)
{
- Assert (this->mapping_info->second_derivatives_initialized == true,
- ExcNotInitialized());
- const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[q_point];
- const Tensor<2,dim,VectorizedArray<Number> > &jac_grad = jacobian_grad[q_point];
for (unsigned int comp=0; comp<n_components; comp++)
{
// compute laplacian before the gradient because it needs to access
for (unsigned int f=1; f<dim; ++f)
hessian_out[comp][d] += jac[d][f] * tmp[f][d];
}
-
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- hessian_out[comp][d] += (jac_grad[d][e] *
- this->gradients_quad[comp][e][q_point]);
}
}
- // cell with general Jacobian, but constant within the cell
- else // if (this->cell_type == internal::MatrixFreeFunctions::affine)
+ // cell with general Jacobian
+ else
{
- const Tensor<2,dim,VectorizedArray<Number> > &jac = jacobian[0];
+ const Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > &jac_grad =
+ mapping_data->jacobian_gradients[0][this->get_mapping_data_index_offset()+q_point];
for (unsigned int comp=0; comp<n_components; comp++)
{
// compute laplacian before the gradient because it needs to access
for (unsigned int f=1; f<dim; ++f)
hessian_out[comp][d] += jac[d][f] * tmp[f][d];
}
+
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ hessian_out[comp][d] += (jac_grad[d][e] *
+ this->gradients_quad[comp][e][q_point]);
}
}
return hessian_out;
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
AssertIndexRange (q_point, this->data->n_q_points);
this->gradients_quad_submitted = true;
+ Assert (this->J_value != nullptr, ExcNotInitialized());
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
#endif
+
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
const VectorizedArray<Number> JxW = J_value[0] * quadrature_weights[q_point];
for (unsigned int comp=0; comp<n_components; comp++)
for (unsigned int d=0; d<dim; ++d)
this->gradients_quad[comp][d][q_point] = (grad_in[comp][d] *
- cartesian_data[0][d] * JxW);
+ jacobian[0][d][d] * JxW);
}
else
{
const Tensor<2,dim,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
jacobian[q_point] : jacobian[0];
const VectorizedArray<Number> JxW =
- this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
J_value[q_point] : J_value[0] * quadrature_weights[q_point];
for (unsigned int comp=0; comp<n_components; ++comp)
for (unsigned int d=0; d<dim; ++d)
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, this->data->n_q_points);
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
+
Tensor<1,dim,VectorizedArray<Number> > grad_out;
- // Cartesian cell
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
for (unsigned int d=0; d<dim; ++d)
grad_out[d] = (this->gradients_quad[0][d][q_point] *
- this->cartesian_data[0][d]);
+ this->jacobian[0][d][d]);
}
- // cell with general/constant Jacobian
+ // cell with general/affine Jacobian
else
{
const Tensor<2,dim,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] : this->jacobian[0];
+ this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
+ q_point : 0];
for (unsigned int d=0; d<dim; ++d)
{
- grad_out[d] = (jac[d][0] * this->gradients_quad[0][0][q_point]);
+ grad_out[d] = jac[d][0] * this->gradients_quad[0][0][q_point];
for (unsigned int e=1; e<dim; ++e)
- grad_out[d] += (jac[d][e] * this->gradients_quad[0][e][q_point]);
+ grad_out[d] += jac[d][e] * this->gradients_quad[0][e][q_point];
}
}
return grad_out;
void
FEEvaluationAccess<dim,1,Number>
::submit_gradient (const Tensor<1,dim,VectorizedArray<Number> > grad_in,
- const unsigned int q_point)
+ const unsigned int q_index)
{
#ifdef DEBUG
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
- AssertIndexRange (q_point, this->data->n_q_points);
+ AssertIndexRange (q_index, this->data->n_q_points);
this->gradients_quad_submitted = true;
+ Assert (this->J_value != nullptr, ExcNotInitialized());
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
#endif
+
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
- const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
+ const VectorizedArray<Number> JxW = this->J_value[0] *
+ this->quadrature_weights[q_index];
for (unsigned int d=0; d<dim; ++d)
- this->gradients_quad[0][d][q_point] = (grad_in[d] *
- this->cartesian_data[0][d] *
+ this->gradients_quad[0][d][q_index] = (grad_in[d] *
+ this->jacobian[0][d][d] *
JxW);
}
// general/affine cell type
else
{
const Tensor<2,dim,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] : this->jacobian[0];
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
+ this->jacobian[q_index] : this->jacobian[0];
const VectorizedArray<Number> JxW =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
+ this->J_value[q_index] : this->J_value[0] * this->quadrature_weights[q_index];
for (unsigned int d=0; d<dim; ++d)
{
VectorizedArray<Number> new_val = jac[0][d] * grad_in[0];
for (unsigned int e=1; e<dim; ++e)
new_val += jac[e][d] * grad_in[e];
- this->gradients_quad[0][d][q_point] = new_val * JxW;
+ this->gradients_quad[0][d][q_index] = new_val * JxW;
}
}
}
Assert (this->gradients_quad_initialized==true,
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, this->data->n_q_points);
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
VectorizedArray<Number> divergence;
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
divergence = (this->gradients_quad[0][0][q_point] *
- this->cartesian_data[0][0]);
+ this->jacobian[0][0][0]);
for (unsigned int d=1; d<dim; ++d)
divergence += (this->gradients_quad[d][d][q_point] *
- this->cartesian_data[0][d]);
+ this->jacobian[0][d][d]);
}
// cell with general/constant Jacobian
else
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
AssertIndexRange (q_point, this->data->n_q_points);
this->gradients_quad_submitted = true;
+ Assert (this->J_value != nullptr, ExcNotInitialized());
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
#endif
+
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
const VectorizedArray<Number> fac = this->J_value[0] *
for (unsigned int d=0; d<dim; ++d)
{
this->gradients_quad[d][d][q_point] = (fac *
- this->cartesian_data[0][d]);
+ this->jacobian[0][d][d]);
for (unsigned int e=d+1; e<dim; ++e)
{
this->gradients_quad[d][e][q_point] = VectorizedArray<Number>();
Assert (this->cell != numbers::invalid_unsigned_int, ExcNotInitialized());
AssertIndexRange (q_point, this->data->n_q_points);
this->gradients_quad_submitted = true;
+ Assert (this->J_value != nullptr, ExcNotInitialized());
+ Assert (this->jacobian != nullptr, ExcNotInitialized());
#endif
+
if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
{
const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
for (unsigned int d=0; d<dim; ++d)
this->gradients_quad[d][d][q_point] = (sym_grad.access_raw_entry(d) *
JxW *
- this->cartesian_data[0][d]);
+ this->jacobian[0][d][d]);
for (unsigned int e=0, counter=dim; e<dim; ++e)
for (unsigned int d=e+1; d<dim; ++d, ++counter)
{
const VectorizedArray<Number> value = sym_grad.access_raw_entry(counter) * JxW;
this->gradients_quad[e][d][q_point] = (value *
- this->cartesian_data[0][d]);
+ this->jacobian[0][d][d]);
this->gradients_quad[d][e][q_point] = (value *
- this->cartesian_data[0][e]);
+ this->jacobian[0][e][e]);
}
}
// general/affine cell type
internal::ExcAccessToUninitializedField());
AssertIndexRange (q_point, this->data->n_q_points);
- Tensor<1,1,VectorizedArray<Number> > grad_out;
+ const Tensor<2,1,VectorizedArray<Number> > &jac =
+ this->cell_type == internal::MatrixFreeFunctions::general ?
+ this->jacobian[q_point] : this->jacobian[0];
- // Cartesian cell
- if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
- {
- grad_out[0] = (this->gradients_quad[0][0][q_point] *
- this->cartesian_data[0][0]);
- }
- // cell with general/constant Jacobian
- else
- {
- const Tensor<2,1,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] : this->jacobian[0];
+ Tensor<1,1,VectorizedArray<Number> > grad_out;
+ grad_out[0] = jac[0][0] * this->gradients_quad[0][0][q_point];
- grad_out[0] = (jac[0][0] * this->gradients_quad[0][0][q_point]);
- }
return grad_out;
}
AssertIndexRange (q_point, this->data->n_q_points);
this->gradients_quad_submitted = true;
#endif
- if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
- {
- const VectorizedArray<Number> JxW = this->J_value[0] * this->quadrature_weights[q_point];
- this->gradients_quad[0][0][q_point] = (grad_in[0] *
- this->cartesian_data[0][0] *
- JxW);
- }
- // general/affine cell type
- else
- {
- const Tensor<2,1,VectorizedArray<Number> > &jac =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->jacobian[q_point] : this->jacobian[0];
- const VectorizedArray<Number> JxW =
- this->cell_type == internal::MatrixFreeFunctions::general ?
- this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
- this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW;
- }
+ const Tensor<2,1,VectorizedArray<Number> > &jac =
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
+ this->jacobian[q_point] : this->jacobian[0];
+ const VectorizedArray<Number> JxW =
+ this->cell_type > internal::MatrixFreeFunctions::affine ?
+ this->J_value[q_point] : this->J_value[0] * this->quadrature_weights[q_point];
+
+ this->gradients_quad[0][0][q_point] = jac[0][0] * grad_in[0] * JxW;
}
proposed_dof_comp = no;
break;
}
- if (static_n_q_points ==
- this->mapping_info->mapping_data_gen[this->quad_no].n_q_points[this->active_quad_index])
+ if (n_q_points ==
+ this->mapping_data->descriptor[this->active_quad_index].n_q_points)
proposed_quad_comp = this->quad_no;
else
- for (unsigned int no=0; no<this->mapping_info->mapping_data_gen.size(); ++no)
- if (this->mapping_info->mapping_data_gen[no].n_q_points[this->active_quad_index]
- == static_n_q_points)
+ for (unsigned int no=0; no<this->matrix_info->get_mapping_info().cell_data.size(); ++no)
+ if (this->matrix_info->get_mapping_info().cell_data[no].descriptor[this->active_quad_index].n_q_points
+ == n_q_points)
{
proposed_quad_comp = no;
break;
if (fe_no != numbers::invalid_unsigned_int)
{
AssertDimension (n_q_points,
- this->mapping_info->mapping_data_gen[this->quad_no].
- n_q_points[this->active_quad_index]);
+ this->mapping_data->descriptor[this->active_quad_index].n_q_points);
AssertDimension (this->data->dofs_per_component_on_cell * this->n_fe_components,
this->dof_info->dofs_per_cell[this->active_fe_index]);
}
FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number>
::quadrature_point (const unsigned int q) const
{
- Assert (this->mapping_info->quadrature_points_initialized == true,
- ExcNotInitialized());
- Assert (this->quadrature_points != nullptr, ExcNotInitialized());
+ if (this->matrix_info == nullptr)
+ {
+ Assert ((this->mapped_geometry->get_fe_values().get_update_flags() |
+ update_quadrature_points),
+ ExcNotInitialized());
+ }
+ else
+ {
+ Assert(this->mapping_data->quadrature_point_offsets.empty() == false,
+ ExcNotInitialized());
+ }
+
AssertIndexRange (q, n_q_points);
+ const Point<dim,VectorizedArray<Number> > *quadrature_points = &this->mapping_data->
+ quadrature_points[this->mapping_data->quadrature_point_offsets[this->cell]];
+
+ const unsigned int n_q_points_1d_actual =
+ fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;
// Cartesian mesh: not all quadrature points are stored, only the
// diagonal. Hence, need to find the tensor product index and retrieve the
switch (dim)
{
case 1:
- return this->quadrature_points[q];
+ return quadrature_points[q];
case 2:
- point[0] = this->quadrature_points[q%n_q_points_1d][0];
- point[1] = this->quadrature_points[q/n_q_points_1d][1];
+ point[0] = quadrature_points[q%n_q_points_1d_actual][0];
+ point[1] = quadrature_points[q/n_q_points_1d_actual][1];
return point;
case 3:
- point[0] = this->quadrature_points[q%n_q_points_1d][0];
- point[1] = this->quadrature_points[(q/n_q_points_1d)%n_q_points_1d][1];
- point[2] = this->quadrature_points[q/(n_q_points_1d*n_q_points_1d)][2];
+ point[0] = quadrature_points[q%n_q_points_1d_actual][0];
+ point[1] = quadrature_points[(q/n_q_points_1d_actual)%n_q_points_1d_actual][1];
+ point[2] = quadrature_points[q/(n_q_points_1d_actual*n_q_points_1d_actual)][2];
return point;
default:
Assert (false, ExcNotImplemented());
}
// all other cases: just return the respective data as it is fully stored
else
- return this->quadrature_points[q];
+ return quadrature_points[q];
}
unsigned int n_procs;
};
- /**
- * Data type to identify cell type.
- */
- enum CellType
- {
- /**
- * The cell is Cartesian.
- */
- cartesian = 0,
- /**
- * The cell may be described with an affine mapping.
- */
- affine = 1,
- /**
- * There is no special information available for compressing the
- * representation of the cell.
- */
- general = 2,
- /**
- * The cell type is undefined.
- */
- undefined = 3
- };
-
-
- /**
- * A class that is used to compare floating point arrays (e.g.
- * std::vectors, Tensor<1,dim>, etc.). The idea of this class is to
- * consider two arrays as equal if they are the same within a given
- * tolerance. We use this comparator class within an std::map<> of the
- * given arrays. Note that this comparison operator does not satisfy all
- * the mathematical properties one usually wants to have (consider e.g.
- * the numbers a=0, b=0.1, c=0.2 with tolerance 0.15; the operator gives
- * a<c, but neither of a<b? or b<c? is satisfied). This is not a problem
- * in the use cases for this class, but be careful when using it in other
- * contexts.
- */
- template <typename Number>
- struct FPArrayComparator
- {
- FPArrayComparator (const Number scaling);
-
- bool operator() (const std::vector<Number> &v1,
- const std::vector<Number> &v2) const;
-
- template <int dim>
- bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
- const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
-
- template <int dim>
- bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
- const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
-
- Number tolerance;
- };
-
- // Note: Implementation in matrix_free.templates.h
-
} // end of namespace MatrixFreeFunctions
} // end of namespace internal
#include <deal.II/matrix_free/mapping_info.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/mapping_q1.h>
DEAL_II_NAMESPACE_OPEN
const dealii::FEValues<dim> &get_fe_values () const;
/**
- * Return a vector of inverse transpose Jacobians. For compatibility
- * with FEEvaluation, it returns tensors of vectorized arrays, even
- * though all components are equal.
+ * Return a reference to the underlying storage field of type
+ * MappingInfoStorage of the same format as the data fields in
+ * MappingInfo. This ensures compatibility with the precomputed data
+ * fields in the MappingInfo class.
*/
- const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
- get_inverse_jacobians() const;
-
- /**
- * Return a vector of quadrature weights times the Jacobian determinant
- * (JxW). For compatibility with FEEvaluation, it returns tensors of
- * vectorized arrays, even though all components are equal.
- */
- const AlignedVector<VectorizedArray<Number> > &
- get_JxW_values() const;
-
- /**
- * Return a vector of quadrature points in real space on the given cell.
- * For compatibility with FEEvaluation, it returns tensors of vectorized
- * arrays, even though all components are equal.
- */
- const AlignedVector<Point<dim,VectorizedArray<Number> > > &
- get_quadrature_points() const;
-
- /**
- * Return a vector of normal vectors in real space on the given cell.
- * For compatibility with FEEvaluation, it returns tensors of vectorized
- * arrays, even though all components are equal.
- */
- const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
- get_normal_vectors() const;
+ const MappingInfoStorage<dim,dim,Number> &
+ get_data_storage() const;
/**
* Return a reference to 1D quadrature underlying this object.
const Quadrature<1> quadrature_1d;
/**
- * Inverse Jacobians, stored in vectorized array form.
- */
- AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > inverse_jacobians;
-
- /**
- * Stored Jacobian determinants and quadrature weights
- */
- AlignedVector<VectorizedArray<Number> > jxw_values;
-
- /**
- * Stored quadrature points
+ * The storage part created for a single cell and held in analogy to
+ * MappingInfo.
*/
- AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
-
- /**
- * Stored normal vectors (for face integration)
- */
- AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
+ MappingInfoStorage<dim,dim,Number> mapping_info_storage;
};
const UpdateFlags update_flags)
:
fe_values(mapping, fe_dummy, Quadrature<dim>(quadrature),
- internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
- quadrature_1d(quadrature),
- inverse_jacobians(fe_values.get_quadrature().size()),
- jxw_values(fe_values.get_quadrature().size()),
- quadrature_points(fe_values.get_quadrature().size()),
- normal_vectors(fe_values.get_quadrature().size())
+ MappingInfo<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);
+ if (update_flags & update_quadrature_points)
+ {
+ mapping_info_storage.quadrature_point_offsets.resize(1, 0);
+ mapping_info_storage.quadrature_points.resize(fe_values.n_quadrature_points);
+ }
+ if (fe_values.get_update_flags() & update_normal_vectors)
+ {
+ mapping_info_storage.normal_vectors.resize(fe_values.n_quadrature_points);
+ mapping_info_storage.normals_times_jacobians[0].resize(fe_values.n_quadrature_points);
+ }
Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
ExcNotImplemented());
}
MappingDataOnTheFly<dim,Number>::MappingDataOnTheFly (const Quadrature<1> &quadrature,
const UpdateFlags update_flags)
:
- fe_values(fe_dummy, Quadrature<dim>(quadrature),
- internal::MatrixFreeFunctions::MappingInfo<dim,Number>::compute_update_flags(update_flags)),
- quadrature_1d(quadrature),
- inverse_jacobians(fe_values.get_quadrature().size()),
- jxw_values(fe_values.get_quadrature().size()),
- quadrature_points(fe_values.get_quadrature().size()),
- normal_vectors(fe_values.get_quadrature().size())
- {
- Assert(!(fe_values.get_update_flags() & update_jacobian_grads),
- ExcNotImplemented());
- }
+ MappingDataOnTheFly(::dealii::StaticMappingQ1<dim,dim>::mapping,
+ quadrature, update_flags)
+ {}
fe_values.reinit(present_cell);
for (unsigned int q=0; q<fe_values.get_quadrature().size(); ++q)
{
- if (fe_values.get_update_flags() & update_inverse_jacobians)
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- inverse_jacobians[q][d][e] = fe_values.inverse_jacobian(q)[e][d];
+ 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);
+ 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)
for (unsigned int d=0; d<dim; ++d)
- quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
+ mapping_info_storage.quadrature_points[q][d] = fe_values.quadrature_point(q)[d];
if (fe_values.get_update_flags() & update_normal_vectors)
- for (unsigned int d=0; d<dim; ++d)
- normal_vectors[q][d] = fe_values.normal_vector(q)[d];
- if (fe_values.get_update_flags() & update_JxW_values)
- jxw_values[q] = fe_values.JxW(q);
+ {
+ for (unsigned int d=0; d<dim; ++d)
+ mapping_info_storage.normal_vectors[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];
+ }
}
}
template <int dim, typename Number>
inline
- const AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > &
- MappingDataOnTheFly<dim,Number>::get_inverse_jacobians() const
- {
- return inverse_jacobians;
- }
-
-
-
- template <int dim, typename Number>
- inline
- const AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > &
- MappingDataOnTheFly<dim,Number>::get_normal_vectors() const
- {
- return normal_vectors;
- }
-
-
-
- template <int dim, typename Number>
- inline
- const AlignedVector<Point<dim,VectorizedArray<Number> > > &
- MappingDataOnTheFly<dim,Number>::get_quadrature_points() const
- {
- return quadrature_points;
- }
-
-
-
- template <int dim, typename Number>
- inline
- const AlignedVector<VectorizedArray<Number> > &
- MappingDataOnTheFly<dim,Number>::get_JxW_values() const
+ const MappingInfoStorage<dim,dim,Number> &
+ MappingDataOnTheFly<dim,Number>::get_data_storage() const
{
- return jxw_values;
+ return mapping_info_storage;
}
namespace MatrixFreeFunctions
{
/**
- * The class that stores all geometry-dependent data related with cell
- * interiors for use in the matrix-free class.
+ * An enum to identify various types of cells. The most general type is
+ * what we typically compute in the FEValues context but for many
+ * geometries we can save significant storage.
+ */
+ enum CellType : unsigned char
+ {
+ /**
+ * The cell is Cartesian.
+ */
+ cartesian = 0,
+ /**
+ * The cell may be described with an affine mapping.
+ */
+ affine = 1,
+ /**
+ * The current face is flat, i.e., the normal factor on a face is the
+ * same on all quadrature points.
+ */
+ flat_faces = 2,
+ /**
+ * There is no special information available for compressing the
+ * representation of the cell.
+ */
+ general = 3
+ };
+
+
+
+ /**
+ * Definition of a structure that stores all cached data related to the
+ * evaluated geometry from the mapping. In order to support hp-adaptivity
+ * and compressed storage (in particular for Jacobians, JxW values, and
+ * normals), storage length can be different for different rows. Thus, it
+ * allows to jump at the data of individual rows similar to compressed row
+ * storage in sparse matrices. We have two different start indices for
+ * fields with different sizes. The first category of offsets are the
+ * indices for Jacobians of the transformation from unit to real cell (we
+ * store the inverse Jacobian), second derivatives, JxW values, and normal
+ * vectors. We keep separate arrays for all these data structures because
+ * a user code might access only some of them. In such a case, one array
+ * will be gone through in a contiguous order with access to all entries,
+ * which makes it easy for the processor to prefetch data. Having all data
+ * in a single array would require some strides in the access pattern,
+ * which is much more complicated for the processor to predict (and indeed
+ * leads to prefetching of data that does not get used on Intel processors
+ * such as BroadwellEP).
*
- * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
+ * The second category of indices are the offsets for the quadrature
+ * points. Quadrature points can be compressed less than the other fields
+ * and thus need longer fields. Quadrature point indices are often used in
+ * other contexts such as evaluation of right hand sides.
+ *
+ * The third component is a descriptor of data from the unit cells, called
+ * QuadratureDescriptor, which contains the quadrature weights and
+ * permutations of how to go through quadrature points in case of face
+ * data. The latter comes in a vector for the support of hp adaptivity,
+ * with several data fields for the individual quadrature formulas.
+ *
+ * @author Katharina Kormann, Martin Kronbichler, 2018
*/
- template <int dim, typename Number>
- struct MappingInfo
+ template <int structdim, int spacedim, typename Number>
+ struct MappingInfoStorage
{
+ struct QuadratureDescriptor
+ {
+ /**
+ * Constructor. Does nothing.
+ */
+ QuadratureDescriptor();
+
+ /**
+ * Set up the lengths in the various members of this struct.
+ */
+ void initialize(const Quadrature<1> &quadrature_1d,
+ const UpdateFlags update_flags_inner_faces = update_default);
+
+ /**
+ * Returns the memory consumption in bytes.
+ */
+ std::size_t memory_consumption() const;
+
+ /**
+ * Number of quadrature points applied on the given cell or face.
+ */
+ unsigned int n_q_points;
+
+ /**
+ * Quadrature formula applied on the given cell or face.
+ */
+ Quadrature<structdim> quadrature;
+
+ /**
+ * Quadrature weights separated by dimension for use in specific
+ * situations.
+ */
+ std::array<AlignedVector<Number>, structdim> tensor_quadrature_weights;
+
+ /**
+ * A cached vector of quadrature weights in the given number format
+ * (non-vectorized, as it is cheap to broadcast the value to all lanes
+ * when it is used in a vectorized context).
+ */
+ AlignedVector<Number> quadrature_weights;
+
+ /**
+ * For quadrature on faces, the evaluation of basis functions is not
+ * in the correct order if a face is not in the standard orientation
+ * to a given element. This data structure is used to re-order the
+ * data evaluated on quadrature points to represent the correct order.
+ */
+ dealii::Table<2,unsigned int> face_orientations;
+ };
+
/**
- * Determines how many bits of an unsigned int are used to distinguish
- * the cell types (Cartesian, with constant Jacobian, or general)
+ * A class describing the layout of the sections in the @p data_storage
+ * field and also includes some data that depends on the number of
+ * quadrature points in the hp context such as the inner quadrature
+ * formula and re-indexing for faces that are not in the standard
+ * orientation.
*/
- static const std::size_t n_cell_type_bits = 2;
+ std::vector<QuadratureDescriptor> descriptor;
/**
- * Determines how many types of different cells can be detected at most.
- * Corresponds to the number of bits we reserved for it.
+ * Stores the index offset into the arrays @p jxw_values, @p jacobians,
+ * @p normal_vectors and the second derivatives. Note that affine cells
+ * have shorter fields of length 1, where the others have lengths equal
+ * to the numer of quadrature points of the given cell.
*/
- static const unsigned int n_cell_types = 1U<<n_cell_type_bits;
+ AlignedVector<unsigned int> data_index_offsets;
/**
- * An abbreviation for the length of vector lines of the current data
- * type.
+ * The storage of the Jacobian determinant (times the quadrature weight
+ * in case the transformation is non-affine) on quadrature
+ * points.
+ *
+ * Indexed by @p data_index_offsets.
*/
- static const unsigned int n_vector_elements = VectorizedArray<Number>::n_array_elements;
+ AlignedVector<VectorizedArray<Number> > JxW_values;
+ /**
+ * Stores the normal vectors.
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normal_vectors;
+
+ /**
+ * The storage of covariant transformation on quadrature points, i.e.,
+ * the inverse and transposed Jacobians of the transformation from the
+ * unit to the real cell.
+ *
+ * Indexed by @p data_index_offsets.
+ *
+ * Contains two fields for access from both sides for interior faces,
+ * but the default case (cell integrals or boundary integrals) only
+ * fills the zeroth component and ignores the first one.
+ */
+ AlignedVector<Tensor<2,spacedim,VectorizedArray<Number> > > jacobians[2];
+
+ /**
+ * The storage of the gradients of the inverse Jacobian
+ * transformation. Because of symmetry, only the upper diagonal and
+ * diagonal part are needed. The first index runs through the
+ * derivatives, starting with the diagonal and then continuing row-wise,
+ * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then
+ * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index
+ * is the spatial coordinate.
+ *
+ * Indexed by @p data_index_offsets.
+ *
+ * Contains two fields for access from both sides for interior faces,
+ * but the default case (cell integrals or boundary integrals) only
+ * fills the zeroth component and ignores the first one.
+ */
+ AlignedVector<Tensor<1,spacedim *(spacedim+1)/2,
+ Tensor<1,spacedim,VectorizedArray<Number> > > > jacobian_gradients[2];
+
+ /**
+ * Stores the Jacobian transformations times the normal vector (this
+ * represents a shortcut that is accessed often and can thus get higher
+ * performance).
+ *
+ * Indexed by @p data_index_offsets.
+ */
+ AlignedVector<Tensor<1,spacedim,VectorizedArray<Number> > > normals_times_jacobians [2];
+
+ /**
+ * Stores the index offset of a particular cell into the quadrature
+ * points array in real coordinates. Note that Cartesian cells have
+ * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells
+ * (length is @p n_q_points) or faces.
+ */
+ AlignedVector<unsigned int> quadrature_point_offsets;
+
+ /**
+ * Stores the quadrature points in real coordinates, including a
+ * compression scheme for Cartesian cells where we do not need to store
+ * the full data on all points.
+ *
+ * Indexed by @p quadrature_point_offsets.
+ */
+ AlignedVector<Point<spacedim,VectorizedArray<Number> > > quadrature_points;
+
+ /**
+ * Returns the quadrature index for a given number of quadrature
+ * points. If not in hp mode or if the index is not found, this
+ * function always returns index 0. Hence, this function does not
+ * check whether the given degree is actually present.
+ */
+ unsigned int
+ quad_index_from_n_q_points (const unsigned int n_q_points) const;
+
+ /**
+ * Prints a detailed summary of memory consumption in the different
+ * structures of this class to the given output stream.
+ */
+ template <typename StreamType>
+ void print_memory_consumption(StreamType &out,
+ const SizeInfo &task_info) const;
+
+ /**
+ * Returns the memory consumption in bytes.
+ */
+ std::size_t memory_consumption () const;
+ };
+
+
+
+ /**
+ * The class that stores all geometry-dependent data related with cell
+ * interiors for use in the matrix-free class.
+ *
+ * @author Katharina Kormann and Martin Kronbichler, 2010, 2011, 2017
+ */
+ template <int dim, typename Number>
+ struct MappingInfo
+ {
/**
* Empty constructor.
*/
MappingInfo();
/**
- * Compute the information in the given cells. The cells are specified
- * by the level and the index within the level (as given by
+ * Compute the information in the given cells and faces. The cells are
+ * specified by the level and the index within the level (as given by
* CellIterator::level() and CellIterator::index(), in order to allow
* for different kinds of iterators, e.g. standard DoFHandler,
* multigrid, etc.) on a fixed Triangulation. In addition, a mapping
* and several quadrature formulas are given.
*/
void initialize (const dealii::Triangulation<dim> &tria,
- const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+ const std::vector<std::pair<unsigned int,unsigned int> > &cells,
const std::vector<unsigned int> &active_fe_index,
const Mapping<dim> &mapping,
const std::vector<dealii::hp::QCollection<1> > &quad,
- const UpdateFlags update_flags);
-
- /**
- * Helper function to determine which update flags must be set in the
- * internal functions to initialize all data as requested by the user.
- */
- static UpdateFlags
- compute_update_flags (const UpdateFlags update_flags,
- const std::vector<dealii::hp::QCollection<1> > &quad =
- std::vector<dealii::hp::QCollection<1> >());
+ const UpdateFlags update_flags_cells);
/**
* Return the type of a given cell as detected during initialization.
*/
CellType get_cell_type (const unsigned int cell_chunk_no) const;
- /**
- * Return the type of a given cell as detected during initialization.
- */
- unsigned int get_cell_data_index (const unsigned int cell_chunk_no) const;
-
/**
* Clear all data fields in this class.
*/
*/
template <typename StreamType>
void print_memory_consumption(StreamType &out,
- const SizeInfo &size_info) const;
+ const SizeInfo &task_info) const;
/**
- * Stores whether a cell is Cartesian, has constant transform data
- * (Jacobians) or is general. cell_type % 4 gives this information (0:
- * Cartesian, 1: constant Jacobian throughout cell, 2: general cell),
- * and cell_type / 4 gives the index in the data field of where to find
- * the information in the fields Jacobian and JxW values (except for
- * quadrature points, for which the index runs as usual).
+ * Stores whether a cell is Cartesian (cell type 0), has constant
+ * transform data (Jacobians) (cell type 1), or is general (cell type
+ * 3). Type 2 is only used for faces and no cells are assigned this
+ * value.
*/
- std::vector<unsigned int> cell_type;
+ std::vector<CellType> cell_type;
/**
- * The first field stores the inverse Jacobian for Cartesian cells:
- * There, it is a diagonal rank-2 tensor, so we actually just store a
- * rank-1 tensor. It is the same on all cells, therefore we only store
- * it once per cell, and use similarities from one cell to another, too
- * (on structured meshes, there are usually many cells with the same
- * Jacobian).
- *
- * The second field stores the Jacobian determinant for Cartesian cells
- * (without the quadrature weight, which depends on the quadrature
- * point, whereas the determinant is the same on each quadrature point).
+ * Stores whether a face (and both cells adjacent to the face) is
+ * Cartesian (face type 0), whether it represents an affine situation
+ * (face type 1), whether it is a flat face where the normal vector is
+ * the same throughout the face (face type 2), or is general (face type
+ * 3).
*/
- AlignedVector<std::pair<Tensor<1,dim,VectorizedArray<Number> >,
- VectorizedArray<Number> > > cartesian_data;
+ std::vector<CellType> face_type;
/**
- * The first field stores the Jacobian for non-Cartesian cells where all
- * the Jacobians on the cell are the same (i.e., constant, which comes
- * from a linear transformation from unit to real cell). Also use
- * similarities from one cell to another (on structured meshes, there
- * are usually many cells with the same Jacobian).
- *
- * The second field stores the Jacobian determinant for non-Cartesian
- * cells with constant Jacobian throughout the cell (without the
- * quadrature weight, which depends on the quadrature point, whereas the
- * determinant is the same on each quadrature point).
+ * The data cache for the cells.
*/
- AlignedVector<std::pair<Tensor<2,dim,VectorizedArray<Number> >,
- VectorizedArray<Number> > > affine_data;
+ std::vector<MappingInfoStorage<dim,dim,Number> > cell_data;
/**
- * Definition of a structure that stores data that depends on the
- * quadrature formula (if we have more than one quadrature formula on a
- * given problem, these fields will be different)
+ * The data cache for the faces.
*/
- struct MappingInfoDependent
- {
- /**
- * This field stores the row starts for the inverse Jacobian
- * transformations, quadrature weights and second derivatives.
- */
- std::vector<unsigned int> rowstart_jacobians;
-
- /**
- * This field stores the inverse Jacobian transformation from unit to
- * real cell, which is needed for most gradient transformations
- * (corresponds to FEValues::inverse_jacobian) for general cells.
- */
- AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians;
-
- /**
- * This field stores the Jacobian determinant times the quadrature
- * weights (JxW in deal.II speak) for general cells.
- */
- AlignedVector<VectorizedArray<Number> > JxW_values;
-
- /**
- * Stores the diagonal part of the gradient of the inverse Jacobian
- * transformation. The first index runs over the derivatives
- * $\partial^2/\partial x_i^2$, the second over the space coordinate.
- * Needed for computing the Laplacian of FE functions on the real
- * cell. Uses a separate storage from the off-diagonal part
- * $\partial^2/\partial x_i \partial x_j, i\neq j$ because that is
- * only needed for computing a full Hessian.
- */
- AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > jacobians_grad_diag;
-
- /**
- * Stores the off-diagonal part of the gradient of the inverse
- * Jacobian transformation. Because of symmetry, only the upper
- * diagonal part is needed. The first index runs through the
- * derivatives row-wise, i.e., $\partial^2/\partial x_1 \partial x_2$
- * first, then $\partial^2/\partial x_1 \partial x_3$, and so on. The
- * second index is the spatial coordinate. Not filled currently.
- */
- AlignedVector<Tensor<1,(dim>1?dim*(dim-1)/2:1),
- Tensor<1,dim,VectorizedArray<Number> > > > jacobians_grad_upper;
+ std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data;
- /**
- * Stores the row start for quadrature points in real coordinates for
- * both types of cells. Note that Cartesian cells will have shorter
- * fields (length is @p n_q_points_1d) than non-Cartesian cells
- * (length is @p n_q_points).
- */
- std::vector<unsigned int> rowstart_q_points;
-
- /**
- * Stores the quadrature points in real coordinates for Cartesian
- * cells (does not need to store the full data on all points)
- */
- AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
-
- /**
- * The dim-dimensional quadrature formula underlying the problem
- * (constructed from a 1D tensor product quadrature formula).
- */
- dealii::hp::QCollection<dim> quadrature;
+ /**
+ * The data cache for the face-associated-with-cell topology, following
+ * the @p cell_type variable for the cell types.
+ */
+ std::vector<MappingInfoStorage<dim-1,dim,Number> > face_data_by_cells;
- /**
- * The (dim-1)-dimensional quadrature formula corresponding to face
- * evaluation (constructed from a 1D tensor product quadrature
- * formula).
- */
- dealii::hp::QCollection<dim-1> face_quadrature;
+ /**
+ * Computes the information in the given cells, called within
+ * initialize.
+ */
+ void initialize_cells (const dealii::Triangulation<dim> &tria,
+ const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+ const std::vector<unsigned int> &active_fe_index,
+ const Mapping<dim> &mapping,
+ const std::vector<dealii::hp::QCollection<1> > &quad,
+ const UpdateFlags update_flags_cells);
- /**
- * The number of quadrature points for the current quadrature formula.
- */
- std::vector<unsigned int> n_q_points;
+ /**
+ * Helper function to determine which update flags must be set in the
+ * internal functions to initialize all data as requested by the user.
+ */
+ static UpdateFlags
+ compute_update_flags (const UpdateFlags update_flags,
+ const std::vector<dealii::hp::QCollection<1> > &quad =
+ std::vector<dealii::hp::QCollection<1> >());
+ };
- /**
- * The number of quadrature points for the current quadrature formula
- * when applied to a face. Only set if the quadrature formula is
- * derived from a tensor product, since it is not defined from the
- * full quadrature formula otherwise.
- */
- std::vector<unsigned int> n_q_points_face;
- /**
- * The quadrature weights (vectorized data format) on the unit cell.
- */
- std::vector<AlignedVector<VectorizedArray<Number> > > quadrature_weights;
- /**
- * This variable stores the number of quadrature points for all
- * quadrature indices in the underlying element for easier access to
- * data in the hp case.
- */
- std::vector<unsigned int> quad_index_conversion;
+ /**
+ * A helper class to extract either cell or face data from mapping info
+ * for use in FEEvaluationBase.
+ *
+ * @author Katharina Kormann, Martin Kronbichler, 2018
+ */
+ template <int, typename, bool> struct MappingInfoCellsOrFaces;
- /**
- * Return the quadrature index for a given number of quadrature
- * points. If not in hp mode or if the index is not found, this
- * function always returns index 0. Hence, this function does not
- * check whether the given degree is actually present.
- */
- unsigned int
- quad_index_from_n_q_points (const unsigned int n_q_points) const;
+ template <int dim, typename Number>
+ struct MappingInfoCellsOrFaces<dim,Number,false>
+ {
+ static const MappingInfoStorage<dim,dim,Number> *
+ get(const MappingInfo<dim,Number> &mapping_info,
+ const unsigned int quad_no)
+ {
+ AssertIndexRange(quad_no, mapping_info.cell_data.size());
+ return &mapping_info.cell_data[quad_no];
+ }
+ };
+ template <int dim, typename Number>
+ struct MappingInfoCellsOrFaces<dim,Number,true>
+ {
+ static const MappingInfoStorage<dim-1,dim,Number> *
+ get(const MappingInfo<dim,Number> &mapping_info,
+ const unsigned int quad_no)
+ {
+ AssertIndexRange(quad_no, mapping_info.face_data.size());
+ return &mapping_info.face_data[quad_no];
+ }
+ };
- /**
- * Prints a detailed summary of memory consumption in the different
- * structures of this class to the given output stream.
- */
- template <typename StreamType>
- void print_memory_consumption(StreamType &out,
- const SizeInfo &size_info) const;
- /**
- * Return the memory consumption in bytes.
- */
- std::size_t memory_consumption () const;
- };
- /**
- * Contains all the stuff that depends on the quadrature formula
- */
- std::vector<MappingInfoDependent> mapping_data_gen;
+ /**
+ * A class that is used to compare floating point arrays (e.g. std::vectors,
+ * Tensor<1,dim>, etc.). The idea of this class is to consider two arrays as
+ * equal if they are the same within a given tolerance. We use this
+ * comparator class within an std::map<> of the given arrays. Note that this
+ * comparison operator does not satisfy all the mathematical properties one
+ * usually wants to have (consider e.g. the numbers a=0, b=0.1, c=0.2 with
+ * tolerance 0.15; the operator gives a<c, but neither of a<b? or b<c? is
+ * satisfied). This is not a problem in the use cases for this class, but be
+ * careful when using it in other contexts.
+ */
+ template <typename Number>
+ struct FPArrayComparator
+ {
+ FPArrayComparator (const Number scaling);
- /**
- * Stores whether JxW values have been initialized
- */
- bool JxW_values_initialized;
+ bool operator() (const std::vector<Number> &v1,
+ const std::vector<Number> &v2) const;
- /**
- * Stores whether we computed second derivatives.
- */
- bool second_derivatives_initialized;
+ bool operator ()(const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t1,
+ const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t2) const;
- /**
- * Stores whether we computed quadrature points.
- */
- bool quadrature_points_initialized;
+ template <int dim>
+ bool operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+ const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
- /**
- * Internal temporary data used for the initialization.
- */
- struct CellData
- {
- CellData (const double jac_size);
- void resize (const unsigned int size);
-
- AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > quadrature_points;
- AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > general_jac;
- AlignedVector<Tensor<3,dim,VectorizedArray<Number> > > general_jac_grad;
- Tensor<2,dim,VectorizedArray<Number> > const_jac;
- const double jac_size;
- };
+ template <int dim>
+ bool operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+ const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const;
- /**
- * Helper function called internally during the initialize function.
- */
- void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
- const std::pair<unsigned int,unsigned int> *cells,
- const unsigned int cell,
- const unsigned int my_q,
- CellType (&cell_t_prev)[n_vector_elements],
- CellType (&cell_t)[n_vector_elements],
- dealii::FEValues<dim,dim> &fe_values,
- CellData &cell_data) const;
+ Number tolerance;
};
/* ------------------- inline functions ----------------------------- */
- template <int dim, typename Number>
+ template <int structdim, int spacedim, typename Number>
inline
unsigned int
- MappingInfo<dim,Number>::MappingInfoDependent::
- quad_index_from_n_q_points (const unsigned int n_q_points) const
+ MappingInfoStorage<structdim,spacedim,Number>
+ ::quad_index_from_n_q_points (const unsigned int n_q_points) const
{
- for (unsigned int i=0; i<quad_index_conversion.size(); ++i)
- if (n_q_points == quad_index_conversion[i])
+ for (unsigned int i=0; i<descriptor.size(); ++i)
+ if (n_q_points == descriptor[i].n_q_points)
return i;
return 0;
}
MappingInfo<dim,Number>::get_cell_type (const unsigned int cell_no) const
{
AssertIndexRange (cell_no, cell_type.size());
- CellType enum_cell_type = (CellType)(cell_type[cell_no] % n_cell_types);
- Assert(enum_cell_type != undefined, ExcInternalError());
- return enum_cell_type;
- }
-
-
-
- template <int dim, typename Number>
- inline
- unsigned int
- MappingInfo<dim,Number>::get_cell_data_index (const unsigned int cell_no) const
- {
- AssertIndexRange (cell_no, cell_type.size());
- return cell_type[cell_no] >> n_cell_type_bits;
+ return cell_type[cell_no];
}
} // end of namespace MatrixFreeFunctions
#include <deal.II/base/utilities.h>
#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/thread_management.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
{
namespace MatrixFreeFunctions
{
- // ----------------- actual MappingInfo functions -------------------------
+
+ /* ------------------------ MappingInfoStorage implementation ---------- */
+
+ template <int structdim, int spacedim, typename Number>
+ MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+ ::QuadratureDescriptor()
+ :
+ n_q_points (numbers::invalid_unsigned_int)
+ {
+ }
+
+
+
+ template <int structdim, int spacedim, typename Number>
+ void
+ MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+ ::initialize(const Quadrature<1> &quadrature_1d,
+ const UpdateFlags update_flags_inner_faces)
+ {
+ Assert(structdim+1 <= spacedim ||
+ update_flags_inner_faces == update_default,
+ ExcMessage("Volume cells do not allow for setting inner faces"));
+ quadrature = Quadrature<structdim>(quadrature_1d);
+ n_q_points = quadrature.size();
+ quadrature_weights.resize(n_q_points);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ quadrature_weights[i] = quadrature.weight(i);
+
+ for (unsigned int d=0; d<structdim; ++d)
+ {
+ tensor_quadrature_weights[d].resize(quadrature_1d.size());
+ for (unsigned int i=0; i<quadrature_1d.size(); ++i)
+ tensor_quadrature_weights[d][i] = quadrature_1d.weight(i);
+ }
+
+ // face orientation for faces in 3D
+ if (structdim == spacedim-1 && spacedim == 3 &&
+ update_flags_inner_faces != update_default)
+ {
+ const unsigned int n=quadrature_1d.size();
+ face_orientations.reinit(8, n*n);
+ for (unsigned int j=0, i=0; j<n; ++j)
+ for (unsigned int k=0; k<n; ++k, ++i)
+ {
+ // face_orientation=true, face_flip=false, face_rotation=false
+ face_orientations[0][i] = i;
+ // face_orientation=false, face_flip=false, face_rotation=false
+ face_orientations[1][i] = j + k *n;
+ // face_orientation=true, face_flip=true, face_rotation=false
+ face_orientations[2][i] = (n-1-k) + (n-1-j)*n;
+ // face_orientation=false, face_flip=true, face_rotation=false
+ face_orientations[3][i] = (n-1-j) + (n-1-k)*n;
+ // face_orientation=true, face_flip=false, face_rotation=true
+ face_orientations[4][i] = j + (n-1-k)*n;
+ // face_orientation=false, face_flip=false, face_rotation=true
+ face_orientations[5][i] = k + (n-1-j)*n;
+ // face_orientation=true, face_flip=true, face_rotation=true
+ face_orientations[6][i] = (n-1-j) + k *n;
+ // face_orientation=false, face_flip=true, face_rotation=true
+ face_orientations[7][i] = (n-1-k) + j *n;
+ }
+ }
+ }
+
+
+
+ template <int structdim, int spacedim, typename Number>
+ std::size_t
+ MappingInfoStorage<structdim, spacedim, Number>::QuadratureDescriptor
+ ::memory_consumption() const
+ {
+ std::size_t memory =
+ sizeof (this) +
+ quadrature.memory_consumption() +
+ quadrature_weights.memory_consumption() +
+ face_orientations.memory_consumption();
+ for (unsigned int d=0; d<structdim; ++d)
+ memory += tensor_quadrature_weights[d].memory_consumption();
+ return memory;
+ }
+
+
+
+ template <int structdim, int spacedim, typename Number>
+ std::size_t
+ MappingInfoStorage<structdim,spacedim,Number>::memory_consumption() const
+ {
+ return
+ MemoryConsumption::memory_consumption (descriptor) +
+ MemoryConsumption::memory_consumption (data_index_offsets) +
+ MemoryConsumption::memory_consumption (JxW_values) +
+ MemoryConsumption::memory_consumption (normal_vectors) +
+ MemoryConsumption::memory_consumption (jacobians[0]) +
+ MemoryConsumption::memory_consumption (jacobians[1]) +
+ MemoryConsumption::memory_consumption (jacobian_gradients[0]) +
+ MemoryConsumption::memory_consumption (jacobian_gradients[1]) +
+ MemoryConsumption::memory_consumption (normals_times_jacobians[0]) +
+ MemoryConsumption::memory_consumption (normals_times_jacobians[1]) +
+ MemoryConsumption::memory_consumption (quadrature_point_offsets) +
+ MemoryConsumption::memory_consumption (quadrature_points);
+ }
+
+
+
+ template <int structdim, int spacedim, typename Number>
+ template <typename StreamType>
+ void
+ MappingInfoStorage<structdim,spacedim,Number>::print_memory_consumption
+ (StreamType &out,
+ const SizeInfo &task_info) const
+ {
+ // print_memory_statistics involves global communication, so we can
+ // disable the check here only if no processor has any such data
+ const std::size_t size = Utilities::MPI::sum(jacobians[0].size(),
+ task_info.communicator);
+ if (size > 0)
+ {
+ out << " Memory JxW data: ";
+ task_info.print_memory_statistics
+ (out, MemoryConsumption::memory_consumption (data_index_offsets) +
+ MemoryConsumption::memory_consumption (JxW_values));
+ out << " Memory Jacobian data: ";
+ task_info.print_memory_statistics
+ (out, MemoryConsumption::memory_consumption (jacobians[0]) +
+ MemoryConsumption::memory_consumption (jacobians[1]));
+ out << " Memory second derivative data: ";
+ task_info.print_memory_statistics
+ (out, MemoryConsumption::memory_consumption (jacobian_gradients[0]) +
+ MemoryConsumption::memory_consumption (jacobian_gradients[1]));
+ }
+ const std::size_t normal_size = Utilities::MPI::sum(normal_vectors.size(),
+ task_info.communicator);
+ if (normal_size > 0)
+ {
+ out << " Memory normal vectors data: ";
+ task_info.print_memory_statistics
+ (out, MemoryConsumption::memory_consumption (normal_vectors) +
+ MemoryConsumption::memory_consumption (normals_times_jacobians[0]) +
+ MemoryConsumption::memory_consumption (normals_times_jacobians[1]));
+ }
+
+ const std::size_t quad_size =
+ Utilities::MPI::sum(quadrature_points.size(),
+ task_info.communicator);
+ if (quad_size > 0)
+ {
+ out << " Memory quadrature points: ";
+ task_info.print_memory_statistics
+ (out, MemoryConsumption::memory_consumption (quadrature_point_offsets) +
+ MemoryConsumption::memory_consumption (quadrature_points));
+ }
+ }
+
+
+
+ /* ------------------------ MappingInfo implementation ----------------- */
template <int dim, typename Number>
MappingInfo<dim,Number>::MappingInfo()
- :
- JxW_values_initialized (false),
- second_derivatives_initialized (false),
- quadrature_points_initialized (false)
{}
void
MappingInfo<dim,Number>::clear ()
{
- JxW_values_initialized = false;
- quadrature_points_initialized = false;
- second_derivatives_initialized = false;
- mapping_data_gen.clear();
+ cell_data.clear();
+ face_data.clear();
+ face_data_by_cells.clear();
cell_type.clear();
- cartesian_data.clear();
- affine_data.clear();
+ face_type.clear();
}
compute_update_flags (const UpdateFlags update_flags,
const std::vector<dealii::hp::QCollection<1> > &quad)
{
- // this class is build around the evaluation this class is build around
- // the evaluation of inverse gradients, so compute them in any case
- UpdateFlags new_flags = update_inverse_jacobians;
-
- // if the user requested gradients, need inverse Jacobians
- if (update_flags & update_gradients || update_flags & update_inverse_jacobians)
- new_flags |= update_inverse_jacobians;
-
- // for JxW, would only need JxW values.
- if (update_flags & update_JxW_values)
- new_flags |= update_JxW_values;
+ // this class is build around the evaluation of jacobians, so compute
+ // them in any case. The Jacobians will be inverted manually. Since we
+ // always do support integration, we also include the JxW values
+ UpdateFlags new_flags = update_jacobians | update_JxW_values;
// for Hessian information, need inverse Jacobians and the derivative of
// Jacobians (these two together will give use the gradients of the
- namespace internal
+ template <int dim, typename Number>
+ void
+ MappingInfo<dim,Number>::initialize
+ (const dealii::Triangulation<dim> &tria,
+ const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+ const std::vector<unsigned int> &active_fe_index,
+ const Mapping<dim> &mapping,
+ const std::vector<dealii::hp::QCollection<1> > &quad,
+ const UpdateFlags update_flags_cells)
+ {
+ clear();
+
+ // Could call these functions in parallel, but not useful because the
+ // work inside is nicely split up already
+ initialize_cells(tria, cells, active_fe_index, mapping, quad, update_flags_cells);
+ }
+
+
+
+ /* ------------------------- initialization of cells ------------------- */
+
+ // Anonymous namespace with implementation of extraction of values on cell
+ // range
+ namespace
{
template <int dim>
double get_jacobian_size (const dealii::Triangulation<dim> &tria)
return 1;
else return tria.begin()->diameter();
}
- }
- template <int dim, typename Number>
- void
- MappingInfo<dim,Number>::initialize
- (const dealii::Triangulation<dim> &tria,
- const std::vector<std::pair<unsigned int,unsigned int> > &cells,
- const std::vector<unsigned int> &active_fe_index,
- const Mapping<dim> &mapping,
- const std::vector<dealii::hp::QCollection<1> > &quad,
- const UpdateFlags update_flags_input)
- {
- clear();
- const unsigned int n_quads = quad.size();
- const unsigned int n_cells = cells.size();
- const unsigned int vectorization_length =
- VectorizedArray<Number>::n_array_elements;
- Assert (n_cells%vectorization_length == 0, ExcInternalError());
- const unsigned int n_macro_cells = n_cells/vectorization_length;
- mapping_data_gen.resize (n_quads);
- cell_type.resize (n_macro_cells);
+ template <int dim, typename Number>
+ struct CompressedCellData
+ {
+ CompressedCellData(const double expected_size)
+ :
+ data(FPArrayComparator<Number>(expected_size))
+ {}
+
+ std::map<Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> >
+ , unsigned int, FPArrayComparator<Number> > data;
+ };
+
+ /**
+ * Internal temporary data used for the initialization.
+ */
+ template <int dim, typename Number>
+ struct LocalData
+ {
+ LocalData (const double jac_size);
+ void resize (const unsigned int size);
- // dummy FE that is used to set up an FEValues object. Do not need the
- // actual finite element because we will only evaluate quantities for
- // the mapping that are independent of the FE
- FE_Nothing<dim> dummy_fe;
- UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
+ AlignedVector<Point<dim,VectorizedArray<Number> > > quadrature_points;
+ AlignedVector<Tensor<2,dim,VectorizedArray<Number> > > general_jac;
+ AlignedVector<VectorizedArray<Number> > JxW_values;
+ AlignedVector<Tensor<3,dim,VectorizedArray<Number> > > general_jac_grad;
+ AlignedVector<Tensor<1,dim,VectorizedArray<Number> > > normal_vectors;
+ Tensor<2,dim,VectorizedArray<Number> > const_jac;
+ const double jac_size;
+ };
- if (update_flags & update_JxW_values)
- JxW_values_initialized = true;
- if (update_flags & update_jacobian_grads)
- second_derivatives_initialized = true;
- if (update_flags & update_quadrature_points)
- quadrature_points_initialized = true;
-
- // when we make comparisons about the size of Jacobians we need to know
- // the approximate size of typical entries in Jacobians. We need to fix
- // the Jacobian size once and for all. We choose the diameter of the
- // first cell (on level zero, which is the best accuracy we can hope
- // for, since diameters on finer levels are computed by differences of
- // nearby cells). If the mesh extends over a certain domain, the
- // precision of double values is essentially limited by this precision.
- const double jacobian_size = internal::get_jacobian_size(tria);
-
- // objects that hold the data for up to vectorization_length cells while
- // we fill them up. Only after all vectorization_length cells have been
- // processed, we can insert the data into the data structures of this
- // class
- CellData data (jacobian_size);
- for (unsigned int my_q=0; my_q<n_quads; ++my_q)
- {
- MappingInfoDependent ¤t_data = mapping_data_gen[my_q];
- const unsigned int n_hp_quads = quad[my_q].size();
- AssertIndexRange (0, n_hp_quads);
- current_data.n_q_points.reserve (n_hp_quads);
- current_data.n_q_points_face.reserve (n_hp_quads);
- current_data.quadrature_weights.resize (n_hp_quads);
- std::vector<unsigned int> n_q_points_1d (n_hp_quads),
- step_size_cartesian (n_hp_quads);
- if (n_hp_quads > 1)
- current_data.quad_index_conversion.resize(n_hp_quads);
- for (unsigned int q=0; q<n_hp_quads; ++q)
- {
- n_q_points_1d[q] = quad[my_q][q].size();
- const unsigned int n_q_points =
- Utilities::fixed_power<dim>(n_q_points_1d[q]);
- current_data.n_q_points.push_back (n_q_points);
-
- current_data.n_q_points_face.push_back
- (dim>1 ? Utilities::fixed_power<dim-1>(n_q_points_1d[q]) : 1);
- current_data.quadrature.push_back
- (Quadrature<dim>(quad[my_q][q]));
- current_data.face_quadrature.push_back
- (Quadrature<dim-1>(quad[my_q][q]));
-
- // set quadrature weights in vectorized form
- current_data.quadrature_weights[q].resize(n_q_points);
- for (unsigned int i=0; i<n_q_points; ++i)
- current_data.quadrature_weights[q][i] =
- current_data.quadrature[q].get_weights()[i];
-
- if (n_hp_quads > 1)
- current_data.quad_index_conversion[q] = n_q_points;
-
- // To walk on the diagonal for lexicographic ordering, we have
- // to jump one index ahead in each direction. For direction 0,
- // this is just the next point, for direction 1, it means adding
- // n_q_points_1d, and so on.
- step_size_cartesian[q] = 0;
- unsigned int factor = 1;
- for (unsigned int d=0; d<dim; ++d)
+
+ template <int dim, typename Number>
+ LocalData<dim,Number>::LocalData (const double jac_size_in)
+ :
+ jac_size (jac_size_in)
+ {}
+
+
+
+ template <int dim, typename Number>
+ void
+ LocalData<dim,Number>::resize (const unsigned int size)
+ {
+ if (JxW_values.size() != size)
+ {
+ quadrature_points.resize_fast(size);
+ general_jac.resize_fast(size*2);
+ JxW_values.resize_fast(size);
+ general_jac_grad.resize_fast(size*2);
+ normal_vectors.resize_fast(size);
+ }
+ }
+
+ /**
+ * Helper function called internally during the initialize function.
+ */
+ template <int dim, typename Number>
+ void evaluate_on_cell (const dealii::Triangulation<dim> &tria,
+ const std::pair<unsigned int,unsigned int> *cells,
+ const unsigned int my_q,
+ CellType &cell_t_prev,
+ CellType (&cell_t)[VectorizedArray<Number>::n_array_elements],
+ dealii::FEValues<dim,dim> &fe_val,
+ LocalData<dim,Number> &cell_data)
+ {
+ const unsigned int n_q_points = fe_val.n_quadrature_points;
+ const UpdateFlags update_flags = fe_val.get_update_flags();
+
+ cell_data.const_jac = Tensor<2,dim,VectorizedArray<Number> >();
+
+ // this should be the same value as used in HashValue::scaling (but we
+ // not have that field here)
+ const double zero_tolerance_double = cell_data.jac_size *
+ std::numeric_limits<double>::epsilon() * 1024.;
+ for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
+ {
+ typename dealii::Triangulation<dim>::cell_iterator
+ cell_it (&tria, cells[j].first, cells[j].second);
+ fe_val.reinit(cell_it);
+ cell_t[j] = general;
+
+ // extract quadrature points and store them temporarily. if we have
+ // Cartesian cells, we can compress the indices
+ if (update_flags & update_quadrature_points)
+ for (unsigned int q=0; q<n_q_points; ++q)
{
- step_size_cartesian[q] += factor;
- factor *= n_q_points_1d[q];
+ const Point<dim> &point = fe_val.quadrature_point(q);
+ for (unsigned int d=0; d<dim; ++d)
+ cell_data.quadrature_points[q][d][j] = point[d];
}
- }
- // if there are no cells, there is nothing to do
- if (cells.size() == 0)
- continue;
-
- Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
- Tensor<1,dim,VectorizedArray<Number> > tmp;
-
- // encodes the cell types of the current cell. Since several cells
- // must be considered together, this variable holds the individual
- // info of the last chunk of cells
- CellType cell_t [vectorization_length],
- cell_t_prev [vectorization_length];
- for (unsigned int j=0; j<vectorization_length; ++j)
- cell_t_prev[j] = undefined;
-
- // fe_values object that is used to compute the mapping data. for
- // the hp case there might be more than one finite element. since we
- // manually select the active FE index and not via a
- // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
- // select the correct finite element, so just hold a vector of
- // FEValues
- std::vector<std::shared_ptr<dealii::FEValues<dim> > >
- fe_values (current_data.quadrature.size());
- UpdateFlags update_flags_feval =
- ((update_flags & update_inverse_jacobians) ? update_jacobians : update_default) |
- ((update_flags & update_jacobian_grads) ? update_jacobian_grads : update_default) |
- ((update_flags & update_quadrature_points) ? update_quadrature_points : update_default);
-
- // resize the fields that have fixed size or for which we know
- // something from an earlier loop
- current_data.rowstart_q_points.resize (n_macro_cells+1);
- if (my_q > 0)
- {
- const unsigned int n_cells_var =
- mapping_data_gen[0].rowstart_jacobians.size()-1;
- current_data.rowstart_jacobians.reserve (n_cells_var+1);
- const unsigned int reserve_size = n_cells_var *
- current_data.n_q_points[0];
- if (mapping_data_gen[0].jacobians.size() > 0)
- current_data.jacobians.reserve (reserve_size);
- if (mapping_data_gen[0].JxW_values.size() > 0)
- current_data.jacobians.reserve (reserve_size);
- if (mapping_data_gen[0].jacobians_grad_diag.size() > 0)
- current_data.jacobians_grad_diag.reserve (reserve_size);
- if (mapping_data_gen[0].jacobians_grad_upper.size() > 0)
- current_data.jacobians_grad_upper.reserve (reserve_size);
- }
+ // if this is not the first quadrature formula and we already have
+ // determined that this cell is either Cartesian or with constant
+ // Jacobian, we have nothing more to do.
+ if (my_q > 0 && cell_t_prev <= affine)
+ continue;
+
+ // first round: if the transformation is detected to be the same as
+ // on the old cell, we only need to copy over the data.
+ if (fe_val.get_cell_similarity() == CellSimilarity::translation
+ &&
+ my_q == 0)
+ {
+ if (j==0)
+ cell_t[j] = cell_t_prev;
+ else
+ cell_t[j] = cell_t[j-1];
+ }
+
+ const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0);
+
+ if (my_q == 0)
+ {
+ // check whether the Jacobian is constant on this cell the first
+ // time we come around here
+ if (cell_t[j] == general)
+ {
+ bool jacobian_constant = true;
+ for (unsigned int q=1; q<n_q_points; ++q)
+ {
+ const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ if (std::fabs(jac_0[d][e]-jac[d][e]) >
+ zero_tolerance_double)
+ jacobian_constant = false;
+ if (jacobian_constant == false)
+ break;
+ }
+
+ // check whether the Jacobian is diagonal to machine
+ // accuracy
+ bool cell_cartesian = jacobian_constant;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ if (d!=e)
+ if (std::fabs(jac_0[d][e]) >
+ zero_tolerance_double)
+ {
+ cell_cartesian=false;
+ break;
+ }
+
+ // in case we have only one quadrature point, we can have
+ // non-constant Jacobians, but we cannot detect it by
+ // comparison from one quadrature point to the next: in that
+ // case, need to look at second derivatives and see whether
+ // there are some non-zero entries (this is necessary since
+ // we determine the constness of the Jacobian for the first
+ // quadrature formula and might not look at them any more
+ // for the second, third quadrature formula). in any case,
+ // the flag update_jacobian_grads will be set in that case
+ if (cell_cartesian == false && n_q_points == 1 &&
+ update_flags & update_jacobian_grads)
+ {
+ const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0);
+ const DerivativeForm<2,dim,dim> &jacobian_grad =
+ fe_val.jacobian_grad(0);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ for (unsigned int f=0; f<dim; ++f)
+ {
+ double jac_grad_comp = (jac[f][0] *
+ jacobian_grad[d][e][0]);
+ for (unsigned int g=1; g<dim; ++g)
+ jac_grad_comp += (jac[f][g] *
+ jacobian_grad[d][e][g]);
+ if (std::fabs(jac_grad_comp) >
+ zero_tolerance_double)
+ jacobian_constant = false;
+ }
+ }
+ // set cell type
+ if (cell_cartesian == true)
+ cell_t[j] = cartesian;
+ else if (jacobian_constant == true)
+ cell_t[j] = affine;
+ else
+ cell_t[j] = general;
+ }
+
+ // Cartesian cell
+ if (cell_t[j] == cartesian)
+ {
+ // set Jacobian into diagonal (off-diagonal part is already
+ // zeroed out)
+ for (unsigned int d=0; d<dim; ++d)
+ cell_data.const_jac[d][d][j] = jac_0[d][d];
+ continue;
+ }
- // we would like to put a Tensor<1,dim,VectorizedArray<Number> > as
- // key into the std::map, but std::map allocation does not align the
- // allocated memory correctly, so put it into a tensor of the
- // correct length instead
- FPArrayComparator<Number> comparator(jacobian_size);
- typedef Tensor<1,VectorizedArray<Number>::n_array_elements,Number> VEC_ARRAY;
- std::map<Tensor<1,dim,VEC_ARRAY>, unsigned int,
- FPArrayComparator<Number> > cartesians(comparator);
- std::map<Tensor<2,dim,VEC_ARRAY>, unsigned int,
- FPArrayComparator<Number> > affines(comparator);
-
- // loop over all cells
- for (unsigned int cell=0; cell<n_macro_cells; ++cell)
+ // cell with affine mapping
+ else if (cell_t[j] == affine)
+ {
+ // compress out very small values
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ if (std::fabs(jac_0[d][e]))
+ cell_data.const_jac[d][e][j] = jac_0[d][e];
+ continue;
+ }
+ }
+
+ // general cell case
+
+ // go through all quadrature points and fill in the data into the
+ // temporary data structures with slots for the vectorized data
+ // types
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ // compress out very small numbers which are only noise. Then it
+ // is cleaner to use zero straight away (though it does not save
+ // any memory)
+ const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ cell_data.general_jac[q][d][e][j] =
+ std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
+
+ // need to do some calculus based on the gradient of the
+ // Jacobian, in order to find the gradient of the inverse
+ // Jacobian which is needed in user code. however, we would like
+ // to perform that on vectorized data types instead of doubles
+ // or floats. to this end, copy the gradients first
+ if (update_flags & update_jacobian_grads)
+ {
+ const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ for (unsigned int f=0; f<dim; ++f)
+ cell_data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
+ }
+ }
+ } // end loop over entries of vectorization (n_array_elements cells)
+
+ // set information for next cell
+ cell_t_prev = cell_t[VectorizedArray<Number>::n_array_elements-1];
+ }
+
+
+
+ template <int dim, typename Number>
+ void
+ initialize_cell_range
+ (const std::pair<unsigned int,unsigned int> cell_range,
+ const dealii::Triangulation<dim> &tria,
+ const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+ const std::vector<unsigned int> &active_fe_index,
+ const Mapping<dim> &mapping,
+ const std::vector<dealii::hp::QCollection<1> > &quad,
+ const UpdateFlags update_flags,
+ MappingInfo<dim,Number> &mapping_info,
+ std::pair<std::vector<MappingInfoStorage<dim,dim,Number> >,
+ CompressedCellData<dim,Number> > &data)
+ {
+ FE_Nothing<dim> dummy_fe;
+
+ Tensor<3,dim,VectorizedArray<Number> > jac_grad, grad_jac_inv;
+ Tensor<1,dim,VectorizedArray<Number> > tmp;
+
+ // when we make comparisons about the size of Jacobians we need to
+ // know the approximate size of typical entries in Jacobians. We need
+ // to fix the Jacobian size once and for all. We choose the diameter
+ // of the first cell (on level zero, which is the best accuracy we can
+ // hope for, since diameters on finer levels are computed by
+ // differences of nearby cells) as the order of magnitude by which we
+ // make comparisons "relative."
+ const double jacobian_size = get_jacobian_size(tria);
+
+ // objects that hold the data for up to vectorization_width cells while
+ // we fill them up. Only after all vectorization_width cells have been
+ // processed, we can insert the data into the data structures of this
+ // class
+ LocalData<dim,Number> cell_data (jacobian_size);
+
+ // encodes the cell types of the current cell. Since several cells
+ // must be considered together, this variable holds the individual
+ // info of the last chunk of cells
+ CellType cell_t [VectorizedArray<Number>::n_array_elements];
+ CellType cell_t_prev = general;
+
+ // fe_values object that is used to compute the mapping data. for
+ // the hp case there might be more than one finite element. since we
+ // manually select the active FE index and not via a
+ // hp::DoFHandler<dim>::active_cell_iterator, we need to manually
+ // select the correct finite element, so just hold a vector of
+ // FEValues
+ std::vector<std::vector<std::shared_ptr<dealii::FEValues<dim> > > >
+ fe_values (mapping_info.cell_data.size());
+ for (unsigned int i=0; i<fe_values.size(); ++i)
+ fe_values[i].resize(mapping_info.cell_data[i].descriptor.size());
+ UpdateFlags update_flags_feval =
+ (update_flags & update_jacobians ? update_jacobians : update_default) |
+ (update_flags & update_jacobian_grads ? update_jacobian_grads : update_default) |
+ (update_flags & update_quadrature_points ? update_quadrature_points : update_default);
+
+ std::vector<std::vector<unsigned int> > n_q_points_1d (quad.size()),
+ step_size_cartesian (quad.size());
+ for (unsigned int my_q=0; my_q<quad.size(); ++my_q)
+ {
+ n_q_points_1d[my_q].resize(quad[my_q].size());
+ step_size_cartesian[my_q].resize(quad[my_q].size());
+ for (unsigned int hpq=0; hpq<quad[my_q].size(); ++hpq)
+ {
+ n_q_points_1d[my_q][hpq] = quad[my_q][hpq].size();
+
+ // To walk on the diagonal for lexicographic ordering, we have
+ // to jump one index ahead in each direction. For direction 0,
+ // this is just the next point, for direction 1, it means adding
+ // n_q_points_1d, and so on.
+ step_size_cartesian[my_q][hpq] = 0;
+ unsigned int factor = 1;
+ for (unsigned int d=0; d<dim; ++d)
+ {
+ step_size_cartesian[my_q][hpq] += factor;
+ factor *= n_q_points_1d[my_q][hpq];
+ }
+ }
+ }
+
+ const unsigned int end_cell = std::min(mapping_info.cell_type.size(),
+ std::size_t(cell_range.second));
+ // loop over given cells
+ for (unsigned int cell=cell_range.first; cell<end_cell; ++cell)
+ for (unsigned int my_q=0; my_q<mapping_info.cell_data.size(); ++my_q)
{
// GENERAL OUTLINE: First generate the data in format "number"
- // for vectorization_length cells, and then find the most
+ // for vectorization_width cells, and then find the most
// general type of cell for appropriate vectorized formats. then
// fill this data in
const unsigned int fe_index = active_fe_index.size() > 0 ?
active_fe_index[cell] : 0;
- const unsigned int n_q_points = current_data.n_q_points[fe_index];
- if (fe_values[fe_index].get() == nullptr)
- fe_values[fe_index].reset
- (new dealii::FEValues<dim> (mapping, dummy_fe,
- current_data.quadrature[fe_index],
+ const unsigned int n_q_points =
+ mapping_info.cell_data[my_q].descriptor[fe_index].n_q_points;
+ if (fe_values[my_q][fe_index].get() == nullptr)
+ fe_values[my_q][fe_index].reset
+ (new dealii::FEValues<dim> (mapping, dummy_fe, mapping_info.cell_data[my_q].
+ descriptor[fe_index].quadrature,
update_flags_feval));
- dealii::FEValues<dim> &fe_val = *fe_values[fe_index];
- data.resize (n_q_points);
+ dealii::FEValues<dim> &fe_val = *fe_values[my_q][fe_index];
+ cell_data.resize (n_q_points);
// if the fe index has changed from the previous cell, set the
// old cell type to invalid (otherwise, we might detect
// similarity due to some cells further ahead)
- if (cell > 0 && active_fe_index.size() > 0 &&
- active_fe_index[cell] != active_fe_index[cell-1])
- cell_t_prev[vectorization_length-1] = undefined;
- evaluate_on_cell (tria, &cells[cell*vectorization_length],
- cell, my_q, cell_t_prev, cell_t, fe_val, data);
+ if (my_q > 0)
+ cell_t_prev = CellType(mapping_info.cell_type[cell]);
+ else if (cell > cell_range.first && active_fe_index.size() > 0 &&
+ active_fe_index[cell] != active_fe_index[cell-1])
+ cell_t_prev = general;
+
+ evaluate_on_cell (tria, &cells[cell*VectorizedArray<Number>::n_array_elements],
+ my_q, cell_t_prev, cell_t, fe_val,
+ cell_data);
// now reorder the data into vectorized types. if we are here
// for the first time, we need to find out whether the Jacobian
// allows for some simplification (Cartesian, affine) taking
- // vectorization_length cell together and we have to insert that
- // data into the respective fields. Also, we have to compress
- // different cell indicators into one structure.
+ // vectorization_width cell together
if (my_q == 0)
{
- // find the most general cell type (most general type is 2
+ // find the most general cell type (most general type is 3
// (general cell))
CellType most_general_type = cartesian;
- for (unsigned int j=0; j<vectorization_length; ++j)
+ for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
if (cell_t[j] > most_general_type)
most_general_type = cell_t[j];
- AssertIndexRange (most_general_type, 3);
- unsigned int insert_position = numbers::invalid_unsigned_int;
-
- // Cartesian cell with diagonal Jacobian: only insert the
- // diagonal of the inverse and the Jacobian determinant. We
- // do this by using an std::map that collects pointers to
- // all Cartesian Jacobians. We need a pointer in the
- // std::map because it cannot store data based on
- // VectorizedArray (alignment issue). We circumvent the
- // problem by temporarily filling the next value into the
- // cartesian_data field and, in case we did an insertion,
- // the data is already in the correct place.
- if (most_general_type == cartesian)
- {
- std::pair<Tensor<1,dim,VEC_ARRAY>,unsigned int> new_entry;
- new_entry.second = cartesians.size();
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
- new_entry.first[d][v] = data.const_jac[d][d][v];
-
- std::pair<typename std::map<Tensor<1,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator,
- bool> it = cartesians.insert(new_entry);
- if (it.second == false)
- insert_position = it.first->second;
- else
- insert_position = new_entry.second;
- }
+ AssertIndexRange ((unsigned int)most_general_type, 4U);
+ mapping_info.cell_type[cell] = most_general_type;
+ }
+
+ AssertThrow(data.first[my_q].JxW_values.size() <
+ static_cast<std::size_t>(std::numeric_limits<unsigned int>::max()),
+ ExcMessage("Index overflow. Cannot fit data in 32 bit integers"));
- // Constant Jacobian case. same strategy as before, but with
- // other data fields
- else if (most_general_type == affine)
+ unsigned int insert_position = data.first[my_q].JxW_values.size();
+ // Cartesian/affine cell with constant Jacobians throughout the
+ // cell. We need to store the data in another data field because
+ // std::map cannot store data based on VectorizedArray directly
+ // (alignment issue).
+ if (mapping_info.cell_type[cell] <= affine)
+ {
+ if (my_q == 0)
{
- std::pair<Tensor<2,dim,VEC_ARRAY>,unsigned int> new_entry;
- new_entry.second = affines.size();
+ std::pair<Tensor<2,dim,Tensor<1,VectorizedArray<Number>::
+ n_array_elements,Number> >,unsigned int> new_entry;
+ // This number overlaps with the general data but we
+ // take care of that when we merge data from different
+ // threads
+ new_entry.second = data.second.data.size();
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
- new_entry.first[d][e][v] = data.const_jac[d][e][v];
-
- std::pair<typename std::map<Tensor<2,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator,
- bool> it = affines.insert(new_entry);
- if (it.second == false)
- insert_position = it.first->second;
- else
- insert_position = new_entry.second;
- }
+ new_entry.first[d][e][v] = cell_data.const_jac[d][e][v];
- // general cell case: first resize the data field to fit the
- // new data. if we are here the first time, assume that
- // there are many general cells to come, so reserve some
- // memory in order to not have too many reallocations and
- // memcpy's. The scheme used here involves at most one
- // reallocation.
- else
- {
- Assert (most_general_type == general, ExcInternalError());
- insert_position = current_data.rowstart_jacobians.size();
- if (current_data.rowstart_jacobians.size() == 0)
- {
- unsigned int reserve_size = (n_macro_cells-cell+1)/2;
- current_data.rowstart_jacobians.reserve
- (reserve_size);
- reserve_size *= n_q_points;
- current_data.jacobians.reserve (reserve_size);
- if (update_flags & update_JxW_values)
- current_data.JxW_values.reserve (reserve_size);
- if (update_flags & update_jacobian_grads)
- {
- current_data.jacobians_grad_diag.reserve (reserve_size);
- current_data.jacobians_grad_upper.reserve (reserve_size);
- }
- }
+ insert_position = data.second.data.insert(new_entry).first->second;
}
-
- cell_type[cell] = ((insert_position << n_cell_type_bits) +
- (unsigned int)most_general_type);
-
- } // end if (my_q == 0)
+ else
+ insert_position = data.first[0].data_index_offsets[cell-cell_range.first];
+ }
// general cell case: now go through all quadrature points and
// collect the data. done for all different quadrature formulas,
// so do it outside the above loop.
- if (get_cell_type(cell) == general)
+ data.first[my_q].data_index_offsets.push_back(insert_position);
+ if (mapping_info.get_cell_type(cell) == general)
{
- const unsigned int previous_size =
- current_data.jacobians.size();
- current_data.rowstart_jacobians.push_back (previous_size);
- if (update_flags & update_JxW_values)
- {
- AssertDimension (previous_size,
- current_data.JxW_values.size());
- }
- if (update_flags & update_jacobian_grads)
- {
- AssertDimension (previous_size,
- current_data.jacobians_grad_diag.size());
- AssertDimension (previous_size,
- current_data.jacobians_grad_upper.size());
- }
for (unsigned int q=0; q<n_q_points; ++q)
{
- Tensor<2,dim,VectorizedArray<Number> > &jac = data.general_jac[q];
- Tensor<3,dim,VectorizedArray<Number> > &jacobian_grad = data.general_jac_grad[q];
- for (unsigned int j=0; j<vectorization_length; ++j)
- if (cell_t[j] == cartesian || cell_t[j] == affine)
+ Tensor<2,dim,VectorizedArray<Number> > &jac = cell_data.general_jac[q];
+ Tensor<3,dim,VectorizedArray<Number> > &jacobian_grad = cell_data.general_jac_grad[q];
+ for (unsigned int j=0; j<VectorizedArray<Number>::n_array_elements; ++j)
+ if (cell_t[j] < general)
{
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
{
- jac[d][e][j] = data.const_jac[d][e][j];
+ jac[d][e][j] = cell_data.const_jac[d][e][j];
for (unsigned int f=0; f<dim; ++f)
jacobian_grad[d][e][f][j] = 0.;
}
}
- const VectorizedArray<Number> det = determinant (jac);
- current_data.jacobians.push_back (transpose(invert(jac)));
- const Tensor<2,dim,VectorizedArray<Number> > &inv_jac = current_data.jacobians.back();
-
- if (update_flags & update_JxW_values)
- current_data.JxW_values.push_back
- (det * current_data.quadrature_weights[fe_index][q]);
+ data.first[my_q].JxW_values.push_back(determinant(jac)*
+ fe_val.get_quadrature().weight(q));
+ Tensor<2,dim,VectorizedArray<Number> > inv_jac = transpose(invert(jac));
+ data.first[my_q].jacobians[0].push_back(inv_jac);
if (update_flags & update_jacobian_grads)
{
}
}
- {
- VectorizedArray<Number> grad_diag[dim][dim];
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- grad_diag[d][e] = grad_jac_inv[d][d][e];
- current_data.jacobians_grad_diag.push_back
- (Tensor<2,dim,VectorizedArray<Number> >(grad_diag));
- }
+ // the diagonal part of Jacobian gradient comes first
+ Tensor<1,dim*(dim+1)/2,Tensor<1,dim,VectorizedArray<Number> > > final_grad;
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ final_grad[d][e] = grad_jac_inv[d][d][e];
- // sets upper-diagonal part of Jacobian
- Tensor<1,(dim>1?dim*(dim-1)/2:1),Tensor<1,dim,VectorizedArray<Number> > > grad_upper;
+ // then the upper-diagonal part
for (unsigned int d=0, count=0; d<dim; ++d)
for (unsigned int e=d+1; e<dim; ++e, ++count)
for (unsigned int f=0; f<dim; ++f)
- grad_upper[count][f] = grad_jac_inv[d][e][f];
- current_data.jacobians_grad_upper.push_back(grad_upper);
+ final_grad[dim+count][f] = grad_jac_inv[d][e][f];
+ data.first[my_q].jacobian_gradients[0].push_back(final_grad);
}
}
}
// (with some little indexing) the location of all
// quadrature points.
const unsigned int old_size =
- current_data.quadrature_points.size();
- current_data.rowstart_q_points[cell] = old_size;
+ data.first[my_q].quadrature_points.size();
+ data.first[my_q].quadrature_point_offsets.push_back(old_size);
- Tensor<1,dim,VectorizedArray<Number> > quad_point;
-
- if (get_cell_type(cell) == cartesian)
+ if (mapping_info.get_cell_type(cell) == cartesian)
{
- current_data.quadrature_points.resize (old_size+
- n_q_points_1d[fe_index]);
- for (unsigned int q=0; q<n_q_points_1d[fe_index]; ++q)
- for (unsigned int d=0; d<dim; ++d)
- current_data.quadrature_points[old_size+q][d] =
- data.quadrature_points[q*step_size_cartesian[fe_index]][d];
+ for (unsigned int q=0; q<n_q_points_1d[my_q][fe_index]; ++q)
+ {
+ Point<dim,VectorizedArray<Number> > quad_point;
+ for (unsigned int d=0; d<dim; ++d)
+ quad_point[d] =
+ cell_data.quadrature_points[q*step_size_cartesian[my_q][fe_index]][d];
+ data.first[my_q].quadrature_points.push_back(quad_point);
+ }
+
}
else
{
- current_data.quadrature_points.resize (old_size + n_q_points);
for (unsigned int q=0; q<n_q_points; ++q)
- for (unsigned int d=0; d<dim; ++d)
- current_data.quadrature_points[old_size+q][d] =
- data.quadrature_points[q][d];
+ data.first[my_q].quadrature_points.push_back
+ (cell_data.quadrature_points[q]);
}
}
- } // end for ( cell < n_macro_cells )
- current_data.rowstart_jacobians.push_back
- (current_data.jacobians.size());
- current_data.rowstart_q_points[n_macro_cells] =
- current_data.quadrature_points.size();
-
- // finally, fill the accumulated data for Cartesian and affine cells
- // into cartesian_data and affine_data, invert and transpose the
- // Jacobians, and compute the JxW value.
- if (my_q == 0)
- {
- cartesian_data.resize(cartesians.size());
- for (typename std::map<Tensor<1,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator
- it = cartesians.begin(); it != cartesians.end(); ++it)
- {
- VectorizedArray<Number> det = make_vectorized_array<Number>(1.);
- for (unsigned int d=0; d<dim; ++d)
- {
- VectorizedArray<Number> jac_d;
- for (unsigned int v=0;
- v<VectorizedArray<Number>::n_array_elements; ++v)
- jac_d[v] = it->first[d][v];
- cartesian_data[it->second].first[d] = 1./jac_d;
- det *= jac_d;
- }
- cartesian_data[it->second].second = det;
- }
- affine_data.resize(affines.size());
- for (typename std::map<Tensor<2,dim,VEC_ARRAY>,
- unsigned int, FPArrayComparator<Number> >::iterator
- it = affines.begin(); it != affines.end(); ++it)
- {
- Tensor<2,dim,VectorizedArray<Number> > jac;
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- for (unsigned int v=0;
- v<VectorizedArray<Number>::n_array_elements; ++v)
- jac[d][e][v] = it->first[d][e][v];
+ } // end for ( cell < end_cells )
+ }
- affine_data[it->second].second = determinant(jac);
- affine_data[it->second].first = transpose(invert(jac));
- }
- }
- }
- }
+
+
+ template <typename CONTAINER>
+ void
+ merge_compressed_data(const CONTAINER &source,
+ CONTAINER &destination,
+ std::vector<unsigned int> &indices)
+ {
+ indices.resize(source.size());
+ typename CONTAINER::iterator lookup = destination.begin();
+ for (typename CONTAINER::const_iterator it = source.begin();
+ it != source.end(); ++it)
+ {
+ typename CONTAINER::value_type entry = *it;
+ entry.second = destination.size();
+ lookup = destination.insert(lookup, entry);
+ AssertIndexRange(it->second, indices.size());
+ indices[it->second] = lookup->second;
+ // best guess for insert position of next item
+ ++lookup;
+ }
+ }
+
+
+
+ template <int structdim, int dim, typename Number>
+ void
+ copy_data (const unsigned int first_cell,
+ const std::array<std::size_t,2> &data_shift,
+ const std::vector<unsigned int> &indices_compressed,
+ const std::vector<CellType> &cell_type,
+ MappingInfoStorage<structdim,dim,Number> &data_cells_local,
+ MappingInfoStorage<structdim,dim,Number> &data_cells)
+ {
+ // Copy the index offsets and shift by the appropriate value
+ for (unsigned int lcell=0;
+ lcell<data_cells_local.data_index_offsets.size(); ++lcell)
+ {
+ const unsigned int cell = lcell + first_cell;
+ data_cells.data_index_offsets[cell]
+ = cell_type[cell] <= static_cast<unsigned int>(affine) ?
+ indices_compressed[data_cells_local.data_index_offsets[lcell]]
+ :
+ data_cells_local.data_index_offsets[lcell] + data_shift[0];
+ if (data_cells_local.quadrature_point_offsets.size()>lcell)
+ data_cells.quadrature_point_offsets[cell] =
+ data_cells_local.quadrature_point_offsets[lcell] +
+ data_shift[1];
+ }
+
+ // Copy quadrature points
+ if (data_cells.quadrature_point_offsets.empty() == false)
+ {
+ Point<dim,VectorizedArray<Number> > *out_point =
+ &data_cells.quadrature_points[data_shift[1]];
+ for (const Point<dim,VectorizedArray<Number> > *point =
+ data_cells_local.quadrature_points.begin(); point !=
+ data_cells_local.quadrature_points.end(); ++point, ++out_point)
+ *out_point = *point;
+ data_cells_local.quadrature_points.clear();
+ }
+
+ // If we have collected Jacobian data, copy Jacobians, JxW values,
+ // Jacobian gradients
+ if (data_cells_local.JxW_values.empty())
+ return;
+
+ std::copy(data_cells_local.JxW_values.begin(),
+ data_cells_local.JxW_values.end(),
+ data_cells.JxW_values.begin()+data_shift[0]);
+ data_cells_local.JxW_values.clear();
+ std::copy(data_cells_local.normal_vectors.begin(),
+ data_cells_local.normal_vectors.end(),
+ data_cells.normal_vectors.begin()+data_shift[0]);
+ data_cells_local.normal_vectors.clear();
+ for (unsigned int i=0; i<2; ++i)
+ {
+ std::copy(data_cells_local.jacobians[i].begin(),
+ data_cells_local.jacobians[i].end(),
+ data_cells.jacobians[i].begin()+data_shift[0]);
+ data_cells_local.jacobians[i].clear();
+ std::copy(data_cells_local.jacobian_gradients[i].begin(),
+ data_cells_local.jacobian_gradients[i].end(),
+ data_cells.jacobian_gradients[i].begin()+data_shift[0]);
+ data_cells_local.jacobian_gradients[i].clear();
+ std::copy(data_cells_local.normals_times_jacobians[i].begin(),
+ data_cells_local.normals_times_jacobians[i].end(),
+ data_cells.normals_times_jacobians[i].begin()+data_shift[0]);
+ data_cells_local.normals_times_jacobians[i].clear();
+ }
+ }
+
+ } // end of anonymous namespace
template <int dim, typename Number>
void
- MappingInfo<dim,Number>::evaluate_on_cell (const dealii::Triangulation<dim> &tria,
- const std::pair<unsigned int,unsigned int> *cells,
- const unsigned int cell,
- const unsigned int my_q,
- CellType (&cell_t_prev)[n_vector_elements],
- CellType (&cell_t)[n_vector_elements],
- dealii::FEValues<dim,dim> &fe_val,
- CellData &data) const
+ MappingInfo<dim,Number>::initialize_cells
+ (const dealii::Triangulation<dim> &tria,
+ const std::vector<std::pair<unsigned int,unsigned int> > &cells,
+ const std::vector<unsigned int> &active_fe_index,
+ const Mapping<dim> &mapping,
+ const std::vector<dealii::hp::QCollection<1> > &quad,
+ const UpdateFlags update_flags_input)
{
- const unsigned int n_q_points = fe_val.n_quadrature_points;
- const UpdateFlags update_flags = fe_val.get_update_flags();
-
- // this should be the same value as used in HashValue::scaling (but we
- // not have that field here)
- const double zero_tolerance_double = data.jac_size *
- std::numeric_limits<double>::epsilon() * 1024.;
- for (unsigned int j=0; j<n_vector_elements; ++j)
- {
- typename dealii::Triangulation<dim>::cell_iterator
- cell_it (&tria, cells[j].first, cells[j].second);
- fe_val.reinit(cell_it);
- cell_t[j] = undefined;
+ const unsigned int n_quads = quad.size();
+ const unsigned int n_cells = cells.size();
+ const unsigned int vectorization_width =
+ VectorizedArray<Number>::n_array_elements;
+ Assert (n_cells%vectorization_width == 0, ExcInternalError());
+ const unsigned int n_macro_cells = n_cells/vectorization_width;
+ cell_data.resize (n_quads);
+ cell_type.resize (n_macro_cells);
- // extract quadrature points and store them temporarily. if we have
- // Cartesian cells, we can compress the indices
- if (update_flags & update_quadrature_points)
- for (unsigned int q=0; q<n_q_points; ++q)
- {
- const Point<dim> &point = fe_val.quadrature_point(q);
- for (unsigned int d=0; d<dim; ++d)
- data.quadrature_points[q][d][j] = point[d];
- }
+ // dummy FE that is used to set up an FEValues object. Do not need the
+ // actual finite element because we will only evaluate quantities for
+ // the mapping that are independent of the FE
+ UpdateFlags update_flags = compute_update_flags (update_flags_input, quad);
- // if this is not the first quadrature formula and we already have
- // determined that this cell is either Cartesian or with constant
- // Jacobian, we have nothing more to do.
- if (my_q > 0 && (get_cell_type(cell) == cartesian
- || get_cell_type(cell) == affine) )
- continue;
-
- // first round: if the transformation is detected to be the same as
- // on the old cell, we only need to copy over the data.
- if (fe_val.get_cell_similarity() == CellSimilarity::translation
- &&
- my_q == 0)
- {
- if (j==0)
- {
- Assert (cell>0, ExcInternalError());
- cell_t[j] = cell_t_prev[n_vector_elements-1];
- }
- else
- cell_t[j] = cell_t[j-1];
- }
+ for (unsigned int my_q=0; my_q<n_quads; ++my_q)
+ {
+ const unsigned int n_hp_quads = quad[my_q].size();
+ AssertIndexRange (0, n_hp_quads);
+ cell_data[my_q].descriptor.resize(n_hp_quads);
+ for (unsigned int q=0; q<n_hp_quads; ++q)
+ cell_data[my_q].descriptor[q].initialize(quad[my_q][q],
+ update_default);
+ }
- const DerivativeForm<1,dim,dim> &jac_0 = fe_val.jacobian(0);
+ if (n_macro_cells == 0)
+ return;
- if (my_q == 0)
- {
- // check whether the Jacobian is constant on this cell the first
- // time we come around here
- if (cell_t[j] == undefined)
- {
- bool jacobian_constant = true;
- for (unsigned int q=1; q<n_q_points; ++q)
- {
- const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- if (std::fabs(jac_0[d][e]-jac[d][e]) >
- zero_tolerance_double)
- jacobian_constant = false;
- if (jacobian_constant == false)
- break;
- }
+ // Create as many chunks of cells as we have threads and spawn the work
+ unsigned int work_per_chunk =
+ std::max(8U, (n_macro_cells + MultithreadInfo::n_threads() - 1) /
+ MultithreadInfo::n_threads());
- // check whether the Jacobian is diagonal to machine
- // accuracy
- bool cell_cartesian = jacobian_constant;
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- if (d!=e)
- if (std::fabs(jac_0[d][e]) >
- zero_tolerance_double)
- {
- cell_cartesian=false;
- break;
- }
+ std::vector<std::pair<std::vector<MappingInfoStorage<dim,dim,Number> >,
+ CompressedCellData<dim,Number> > > data_cells_local;
+ // Reserve enough space to avoid re-allocation (which would break the
+ // references to the data fields passed to the tasks!)
+ data_cells_local.reserve(MultithreadInfo::n_threads());
- // in case we have only one quadrature point, we can have
- // non-constant Jacobians, but we cannot detect it by
- // comparison from one quadrature point to the next: in that
- // case, need to look at second derivatives and see whether
- // there are some non-zero entries (this is necessary since
- // we determine the constness of the Jacobian for the first
- // quadrature formula and might not look at them any more
- // for the second, third quadrature formula). in any case,
- // the flag update_jacobian_grads will be set in that case
- if (cell_cartesian == false && n_q_points == 1 &&
- update_flags & update_jacobian_grads)
- {
- const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(0);
- const DerivativeForm<2,dim,dim> &jacobian_grad =
- fe_val.jacobian_grad(0);
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- for (unsigned int f=0; f<dim; ++f)
- {
- double jac_grad_comp = (jac[f][0] *
- jacobian_grad[d][e][0]);
- for (unsigned int g=1; g<dim; ++g)
- jac_grad_comp += (jac[f][g] *
- jacobian_grad[d][e][g]);
- if (std::fabs(jac_grad_comp) >
- zero_tolerance_double)
- jacobian_constant = false;
- }
- }
- // set cell type
- if (cell_cartesian == true)
- cell_t[j] = cartesian;
- else if (jacobian_constant == true)
- cell_t[j] = affine;
- else
- cell_t[j] = general;
- }
+ {
+ Threads::TaskGroup<> tasks;
+ std::pair<unsigned int,unsigned int> cell_range(0U, work_per_chunk);
+ while (cell_range.first < n_macro_cells)
+ {
+ data_cells_local.push_back
+ (std::make_pair (std::vector<MappingInfoStorage<dim,dim,Number> >(n_quads),
+ CompressedCellData<dim,Number>(get_jacobian_size(tria))));
+ tasks += Threads::new_task(&initialize_cell_range<dim,Number>,
+ cell_range, tria,
+ cells, active_fe_index, mapping, quad,
+ update_flags, *this,
+ data_cells_local.back());
+ cell_range.first = cell_range.second;
+ cell_range.second += work_per_chunk;
+ }
+ tasks.join_all();
+ }
- // Cartesian cell
- if (cell_t[j] == cartesian)
- {
- // set Jacobian into diagonal and clear off-diagonal part
- for (unsigned int d=0; d<dim; ++d)
- {
- data.const_jac[d][d][j] = jac_0[d][d];
- for (unsigned int e=d+1; e<dim; ++e)
- {
- data.const_jac[d][e][j] = 0.;
- data.const_jac[e][d][j] = 0.;
- }
- }
- continue;
- }
+ // Fill in each thread's constant Jacobians into the data of the zeroth
+ // chunk in serial
+ std::vector<std::vector<unsigned int> > indices_compressed(data_cells_local.size());
+ for (unsigned int i=0; i<data_cells_local.size(); ++i)
+ merge_compressed_data(data_cells_local[i].second.data,
+ data_cells_local[0].second.data,
+ indices_compressed[i]);
+
+ // Collect all data in the final data fields.
+ // First allocate the memory
+ const unsigned int n_constant_jacobians = data_cells_local[0].second.data.size();
+ for (unsigned int my_q=0; my_q<cell_data.size(); ++my_q)
+ {
+ cell_data[my_q].data_index_offsets.resize(cell_type.size());
+ std::vector<std::array<std::size_t,2> > shift(data_cells_local.size());
+ shift[0][0] = n_constant_jacobians;
+ shift[0][1] = 0;
+ for (unsigned int i=1; i<data_cells_local.size(); ++i)
+ {
+ shift[i][0] = shift[i-1][0] + data_cells_local[i-1].first[my_q].JxW_values.size();
+ shift[i][1] = shift[i-1][1] + data_cells_local[i-1].first[my_q].quadrature_points.size();
+ }
+ cell_data[my_q].JxW_values.
+ resize_fast(shift.back()[0] + data_cells_local.back().first[my_q].
+ JxW_values.size());
+ cell_data[my_q].jacobians[0].resize_fast(cell_data[my_q].JxW_values.size());
+ if (update_flags & update_jacobian_grads)
+ cell_data[my_q].jacobian_gradients[0].resize_fast(cell_data[my_q].JxW_values.size());
+ if (update_flags & update_quadrature_points)
+ {
+ cell_data[my_q].quadrature_point_offsets.resize(cell_type.size());
+ cell_data[my_q].quadrature_points.
+ resize_fast(shift.back()[1] + data_cells_local.back().first[my_q].
+ quadrature_points.size());
+ }
- // cell with affine mapping
- else if (cell_t[j] == affine)
+ // Start tasks that copy the local data
+ Threads::TaskGroup<> tasks;
+ for (unsigned int i=0; i<data_cells_local.size(); ++i)
+ tasks += Threads::new_task(©_data<dim,dim,Number>,
+ work_per_chunk * i, shift[i],
+ indices_compressed[i], cell_type,
+ data_cells_local[i].first[my_q],
+ cell_data[my_q]);
+
+ // finally, insert the constant cell data at the beginning (the
+ // other tasks can already start copying the non-constant data)
+ if (my_q == 0)
+ {
+ for (auto &it : data_cells_local[0].second.data)
{
- // compress out very small values
+ Tensor<2,dim,VectorizedArray<Number> > jac;
for (unsigned int d=0; d<dim; ++d)
for (unsigned int e=0; e<dim; ++e)
- data.const_jac[d][e][j] =
- std::fabs(jac_0[d][e]) < zero_tolerance_double ?
- 0 : jac_0[d][e];
- continue;
+ for (unsigned int v=0; v<VectorizedArray<Number>::n_array_elements; ++v)
+ jac[d][e][v] = it.first[d][e][v];
+ AssertIndexRange(it.second, n_constant_jacobians);
+ const std::size_t index = it.second;
+ cell_data[my_q].JxW_values[index] = determinant(jac);
+ // invert and transpose jac
+ cell_data[my_q].jacobians[0][index] = transpose(invert(jac));
+ // second derivative of transformation is zero on affine cells
}
}
-
- // general cell case
-
- // go through all quadrature points and fill in the data into the
- // temporary data structures with slots for the vectorized data
- // types
- for (unsigned int q=0; q<n_q_points; ++q)
+ else
{
-
- // compress out very small numbers which are only noise. Then it
- // is cleaner to use zero straight away (though it does not save
- // any memory)
- const DerivativeForm<1,dim,dim> &jac = fe_val.jacobian(q);
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- data.general_jac[q][d][e][j] =
- std::fabs(jac[d][e]) < zero_tolerance_double ? 0. : jac[d][e];
-
- // need to do some calculus based on the gradient of the
- // Jacobian, in order to find the gradient of the inverse
- // Jacobian which is needed in user code. however, we would like
- // to perform that on vectorized data types instead of doubles
- // or floats. to this end, copy the gradients first
- if (update_flags & update_jacobian_grads)
+ for (unsigned int i=0; i<n_constant_jacobians; ++i)
{
- const DerivativeForm<2,dim,dim> &jacobian_grad = fe_val.jacobian_grad(q);
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- for (unsigned int f=0; f<dim; ++f)
- data.general_jac_grad[q][d][e][f][j] = jacobian_grad[d][e][f];
+ cell_data[my_q].JxW_values[i] = cell_data[0].JxW_values[i];
+ cell_data[my_q].jacobians[0][i] = cell_data[0].jacobians[0][i];
}
}
- } // end loop over all entries in vectorization (n_vector_elements
- // cells)
- // set information for next cell
- for (unsigned int j=0; j<n_vector_elements; ++j)
- cell_t_prev[j] = cell_t[j];
+ // ... wait for the parallel work to finish
+ tasks.join_all();
+ }
}
+
template <int dim, typename Number>
- MappingInfo<dim,Number>::CellData::CellData (const double jac_size_in)
- :
- jac_size (jac_size_in)
- {}
+ std::size_t MappingInfo<dim,Number>::memory_consumption() const
+ {
+ std::size_t
+ memory = MemoryConsumption::memory_consumption (cell_data);
+ memory += MemoryConsumption::memory_consumption (face_data);
+ memory += cell_type.capacity()*sizeof(CellType);
+ memory += face_type.capacity()*sizeof(CellType);
+ memory += sizeof (*this);
+ return memory;
+ }
template <int dim, typename Number>
- void
- MappingInfo<dim,Number>::CellData::resize (const unsigned int size)
+ template <typename StreamType>
+ void MappingInfo<dim,Number>::print_memory_consumption(StreamType &out,
+ const SizeInfo &task_info) const
{
- if (general_jac.size() != size)
+ out << " Cell types: ";
+ task_info.print_memory_statistics
+ (out, cell_type.capacity()*sizeof(CellType));
+ out << " Face types: ";
+ task_info.print_memory_statistics
+ (out, face_type.capacity()*sizeof(CellType));
+ for (unsigned int j=0; j<cell_data.size(); ++j)
{
- quadrature_points.resize(size);
- general_jac.resize(size);
- general_jac_grad.resize(size);
+ out << " Data component " << j << std::endl;
+ cell_data[j].print_memory_consumption(out, task_info);
+ face_data[j].print_memory_consumption(out, task_info);
}
}
- template <int dim, typename Number>
- std::size_t MappingInfo<dim,Number>::MappingInfoDependent::memory_consumption() const
+ /* ------------------------------------------------------------------ */
+
+ template <typename Number>
+ FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
+ :
+ tolerance (scaling * std::numeric_limits<double>::epsilon() * 1024.)
+ {}
+
+
+
+ template <typename Number>
+ bool
+ FPArrayComparator<Number>::operator() (const std::vector<Number> &v1,
+ const std::vector<Number> &v2) const
{
- std::size_t
- memory = MemoryConsumption::memory_consumption (jacobians);
- memory += MemoryConsumption::memory_consumption (JxW_values);
- memory += MemoryConsumption::memory_consumption (jacobians_grad_diag);
- memory += MemoryConsumption::memory_consumption (jacobians_grad_upper);
- memory += MemoryConsumption::memory_consumption (rowstart_q_points);
- memory += MemoryConsumption::memory_consumption (quadrature_points);
- memory += MemoryConsumption::memory_consumption (quadrature);
- memory += MemoryConsumption::memory_consumption (face_quadrature);
- memory += MemoryConsumption::memory_consumption (quadrature_weights);
- memory += MemoryConsumption::memory_consumption (n_q_points);
- memory += MemoryConsumption::memory_consumption (n_q_points_face);
- memory += MemoryConsumption::memory_consumption (quad_index_conversion);
- return memory;
+ const unsigned int s1 = v1.size(), s2 = v2.size();
+ if (s1 < s2)
+ return true;
+ else if (s1 > s2)
+ return false;
+ else
+ for (unsigned int i=0; i<s1; ++i)
+ if (v1[i] < v2[i] - tolerance)
+ return true;
+ else if (v1[i] > v2[i] + tolerance)
+ return false;
+ return false;
}
- template <int dim, typename Number>
- std::size_t MappingInfo<dim,Number>::memory_consumption() const
+ template <typename Number>
+ bool
+ FPArrayComparator<Number>::
+ operator ()(const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t1,
+ const Tensor<1,VectorizedArray<Number>::n_array_elements,Number> &t2) const
{
- std::size_t
- memory= MemoryConsumption::memory_consumption (mapping_data_gen);
- memory += MemoryConsumption::memory_consumption (affine_data);
- memory += MemoryConsumption::memory_consumption (cartesian_data);
- memory += MemoryConsumption::memory_consumption (cell_type);
- memory += sizeof (*this);
- return memory;
+ for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+ if (t1[k] < t2[k] - tolerance)
+ return true;
+ else if (t1[k] > t2[k] + tolerance)
+ return false;
+ return false;
}
- template <int dim, typename Number>
- template <typename StreamType>
- void MappingInfo<dim,Number>::MappingInfoDependent::print_memory_consumption
- (StreamType &out,
- const SizeInfo &size_info) const
+ template <typename Number>
+ template <int dim>
+ bool
+ FPArrayComparator<Number>::
+ operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+ const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
{
- // print_memory_statistics involves global communication, so we can
- // disable the check here only if no processor has any such data
-#ifdef DEAL_II_WITH_MPI
- unsigned int general_size_glob = 0, general_size_loc = jacobians.size();
- int ierr = MPI_Allreduce (&general_size_loc, &general_size_glob, 1,
- MPI_UNSIGNED, MPI_MAX, size_info.communicator);
- AssertThrowMPI (ierr);
-#else
- unsigned int general_size_glob = jacobians.size();
-#endif
- if (general_size_glob > 0)
- {
- out << " Memory Jacobian data: ";
- size_info.print_memory_statistics
- (out, MemoryConsumption::memory_consumption (jacobians) +
- MemoryConsumption::memory_consumption (JxW_values));
- out << " Memory second derivative data: ";
- size_info.print_memory_statistics
- (out,MemoryConsumption::memory_consumption (jacobians_grad_diag) +
- MemoryConsumption::memory_consumption (jacobians_grad_upper));
- }
-
-#ifdef DEAL_II_WITH_MPI
- unsigned int quad_size_glob = 0, quad_size_loc = quadrature_points.size();
- ierr = MPI_Allreduce (&quad_size_loc, &quad_size_glob, 1, MPI_UNSIGNED,
- MPI_MAX, size_info.communicator);
- AssertThrowMPI (ierr);
-#else
- unsigned int quad_size_glob = quadrature_points.size();
-#endif
- if (quad_size_glob > 0)
- {
- out << " Memory quadrature points: ";
- size_info.print_memory_statistics
- (out, MemoryConsumption::memory_consumption (rowstart_q_points) +
- MemoryConsumption::memory_consumption (quadrature_points));
- }
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+ if (t1[d][k] < t2[d][k] - tolerance)
+ return true;
+ else if (t1[d][k] > t2[d][k] + tolerance)
+ return false;
+ return false;
}
- template <int dim, typename Number>
- template <typename StreamType>
- void MappingInfo<dim,Number>::print_memory_consumption(StreamType &out,
- const SizeInfo &size_info) const
+ template <typename Number>
+ template <int dim>
+ bool
+ FPArrayComparator<Number>::
+ operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
+ const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
{
- out << " Cell types: ";
- size_info.print_memory_statistics
- (out, MemoryConsumption::memory_consumption (cell_type));
- out << " Memory transformations compr: ";
- size_info.print_memory_statistics
- (out, MemoryConsumption::memory_consumption (affine_data) +
- MemoryConsumption::memory_consumption (cartesian_data));
- for (unsigned int j=0; j<mapping_data_gen.size(); ++j)
- {
- out << " Data component " << j << std::endl;
- mapping_data_gen[j].print_memory_consumption(out, size_info);
- }
+ for (unsigned int d=0; d<dim; ++d)
+ for (unsigned int e=0; e<dim; ++e)
+ for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
+ if (t1[d][e][k] < t2[d][e][k] - tolerance)
+ return true;
+ else if (t1[d][e][k] > t2[d][e][k] + tolerance)
+ return false;
+ return false;
}
} // end of namespace MatrixFreeFunctions
MatrixFree<dim,Number>::get_n_q_points(const unsigned int quad_index,
const unsigned int active_fe_index) const
{
- AssertIndexRange (quad_index,
- mapping_info.mapping_data_gen.size());
- return mapping_info.mapping_data_gen[quad_index].n_q_points[active_fe_index];
+ AssertIndexRange (quad_index, mapping_info.cell_data.size());
+ return mapping_info.cell_data[quad_index].descriptor[active_fe_index].n_q_points;
}
MatrixFree<dim,Number>::get_n_q_points_face(const unsigned int quad_index,
const unsigned int active_fe_index) const
{
- AssertIndexRange (quad_index,
- mapping_info.mapping_data_gen.size());
- return mapping_info.mapping_data_gen[quad_index].n_q_points_face[active_fe_index];
+ AssertIndexRange (quad_index, mapping_info.face_data.size());
+ return mapping_info.face_data[quad_index].descriptor[active_fe_index].n_q_points;
}
MatrixFree<dim,Number>::get_quadrature (const unsigned int quad_index,
const unsigned int active_fe_index) const
{
- AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
- return mapping_info.mapping_data_gen[quad_index].
- quadrature[active_fe_index];
+ AssertIndexRange (quad_index, mapping_info.cell_data.size());
+ return mapping_info.cell_data[quad_index].descriptor[active_fe_index].quadrature;
}
MatrixFree<dim,Number>::get_face_quadrature (const unsigned int quad_index,
const unsigned int active_fe_index) const
{
- AssertIndexRange (quad_index, mapping_info.mapping_data_gen.size());
- return mapping_info.mapping_data_gen[quad_index].
- face_quadrature[active_fe_index];
+ AssertIndexRange (quad_index, mapping_info.face_data.size());
+ return mapping_info.face_data[quad_index].descriptor[active_fe_index].quadrature;
}
boundary_cells_start = boundary_cells_end = n_macro_cells;
}
-
-
- /* ------------------------------------------------------------------ */
-
- template <typename Number>
- FPArrayComparator<Number>::FPArrayComparator (const Number scaling)
- :
- tolerance (scaling * std::numeric_limits<double>::epsilon() * 1024.)
- {}
-
-
-
- template <typename Number>
- bool
- FPArrayComparator<Number>::operator() (const std::vector<Number> &v1,
- const std::vector<Number> &v2) const
- {
- const unsigned int s1 = v1.size(), s2 = v2.size();
- if (s1 < s2)
- return true;
- else if (s1 > s2)
- return false;
- else
- for (unsigned int i=0; i<s1; ++i)
- if (v1[i] < v2[i] - tolerance)
- return true;
- else if (v1[i] > v2[i] + tolerance)
- return false;
- return false;
- }
-
-
-
- template <typename Number>
- template <int dim>
- bool
- FPArrayComparator<Number>::
- operator ()(const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
- const Tensor<1,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
- {
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
- if ((t1)[d][k] < (t2)[d][k] - tolerance)
- return true;
- else if ((t1)[d][k] > (t2)[d][k] + tolerance)
- return false;
- return false;
- }
-
-
-
- template <typename Number>
- template <int dim>
- bool
- FPArrayComparator<Number>::
- operator ()(const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t1,
- const Tensor<2,dim,Tensor<1,VectorizedArray<Number>::n_array_elements,Number> > &t2) const
- {
- for (unsigned int d=0; d<dim; ++d)
- for (unsigned int e=0; e<dim; ++e)
- for (unsigned int k=0; k<VectorizedArray<Number>::n_array_elements; ++k)
- if ((t1)[d][e][k] < (t2)[d][e][k] - tolerance)
- return true;
- else if ((t1)[d][e][k] > (t2)[d][e][k] + tolerance)
- return false;
- return false;
- }
}
}
template struct internal::MatrixFreeFunctions::MappingInfo<deal_II_dimension,double>;
template struct internal::MatrixFreeFunctions::MappingInfo<deal_II_dimension,float>;
+ template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension,deal_II_dimension,double>;
+ template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension,deal_II_dimension,float>;
+ template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension-1,deal_II_dimension,double>;
+ template struct internal::MatrixFreeFunctions::MappingInfoStorage<deal_II_dimension-1,deal_II_dimension,float>;
+
#ifndef DEAL_II_MSVC
template
void
mf.reinit (dof, constraints, quad, data);
const unsigned int n_macro_cells = mf.n_macro_cells();
- const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
- const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
- const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+ std::vector<unsigned int> n_cell_types(4, 0);
+ for (unsigned int i=0; i<n_macro_cells; ++i)
+ n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
// should do at least some compression
- Assert(n_cartesian+n_affine+n_general < n_macro_cells, ExcInternalError());
- Assert(n_cartesian * 5 < n_macro_cells, ExcInternalError());
- Assert(n_affine * 10 < n_macro_cells, ExcInternalError());
+ Assert(n_cell_types[0]+n_cell_types[1] > 0, ExcInternalError());
+ Assert(mf.get_mapping_info().cell_data[0].jacobians[0].size() <
+ (n_cell_types[3]*quad.size()+n_macro_cells-n_cell_types[3]),
+ ExcInternalError());
deallog << "OK" << std::endl;
}
mf.reinit (dof, constraints, quad, data);
const unsigned int n_macro_cells = mf.n_macro_cells();
- const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
- const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
- const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+ std::vector<unsigned int> n_cell_types(4, 0);
+ for (unsigned int i=0; i<n_macro_cells; ++i)
+ n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
// should have one Cartesian cell and no other
// cell type
- AssertDimension(n_cartesian, 1);
- AssertDimension(n_affine, 0);
- AssertDimension(n_general, 0);
+ AssertDimension(n_cell_types[0], n_macro_cells);
+ AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1);
Assert(n_macro_cells > 1, ExcInternalError());
deallog << "OK" << std::endl;
}
mf.reinit (dof, constraints, quad, data);
const unsigned int n_macro_cells = mf.n_macro_cells();
- const unsigned int n_cartesian = mf.get_mapping_info().cartesian_data.size();
- const unsigned int n_affine = mf.get_mapping_info().affine_data.size();
- const unsigned int n_general = mf.get_mapping_info().mapping_data_gen[0].rowstart_jacobians.size()-1;
+ std::vector<unsigned int> n_cell_types(4, 0);
+ for (unsigned int i=0; i<n_macro_cells; ++i)
+ n_cell_types[mf.get_mapping_info().get_cell_type(i)]++;
// should have one affine cell and no other
// cell type
- AssertDimension(n_cartesian, 0);
- AssertDimension(n_affine, 1);
- AssertDimension(n_general, 0);
+ AssertDimension(n_cell_types[1], n_macro_cells);
+ AssertDimension(mf.get_mapping_info().cell_data[0].jacobians[0].size(), 1);
Assert(n_macro_cells > 1, ExcInternalError());
deallog << "OK" << std::endl;
}
projects / dealii.git / commitdiff