cudaError_t cuda_error =
cudaMemcpyToSymbol(internal::get_global_shape_values<Number>(),
- shape_info.shape_values.data(),
+ shape_info.data.front().shape_values.data(),
size_shape_values,
0,
cudaMemcpyHostToDevice);
{
cuda_error =
cudaMemcpyToSymbol(internal::get_global_shape_gradients<Number>(),
- shape_info.shape_gradients.data(),
+ shape_info.data.front().shape_gradients.data(),
size_shape_values,
0,
cudaMemcpyHostToDevice);
fe_degree + 1,
n_q_points_1d,
Number>;
- Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
- shape_info.shape_values,
- variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
- shape_info.shape_gradients,
- variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
- shape_info.shape_hessians,
- shape_info.fe_degree + 1,
- shape_info.n_q_points_1d);
+ Eval eval(variant == evaluate_evenodd ?
+ shape_info.data.front().shape_values_eo :
+ shape_info.data.front().shape_values,
+ variant == evaluate_evenodd ?
+ shape_info.data.front().shape_gradients_eo :
+ shape_info.data.front().shape_gradients,
+ variant == evaluate_evenodd ?
+ shape_info.data.front().shape_hessians_eo :
+ shape_info.data.front().shape_hessians,
+ shape_info.data.front().fe_degree + 1,
+ shape_info.data.front().n_q_points_1d);
const unsigned int temp_size =
Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
if (temp_size == 0)
{
temp1 = scratch_data;
- temp2 = temp1 +
- std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
- Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+ shape_info.data.front().fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ shape_info.data.front().n_q_points_1d));
}
else
{
temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
const unsigned int dofs_per_comp =
(type == MatrixFreeFunctions::truncated_tensor) ?
- Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
+ Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
shape_info.dofs_per_component_on_cell;
const Number *values_dofs = values_dofs_actual;
if (type == MatrixFreeFunctions::truncated_tensor)
Number *values_dofs_tmp =
scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
shape_info.n_q_points));
- const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+ const int degree =
+ fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
unsigned int count_p = 0, count_q = 0;
for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
{
fe_degree + 1,
n_q_points_1d,
Number>;
- Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
- shape_info.shape_values,
- variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
- shape_info.shape_gradients,
- variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
- shape_info.shape_hessians,
- shape_info.fe_degree + 1,
- shape_info.n_q_points_1d);
+ Eval eval(variant == evaluate_evenodd ?
+ shape_info.data.front().shape_values_eo :
+ shape_info.data.front().shape_values,
+ variant == evaluate_evenodd ?
+ shape_info.data.front().shape_gradients_eo :
+ shape_info.data.front().shape_gradients,
+ variant == evaluate_evenodd ?
+ shape_info.data.front().shape_hessians_eo :
+ shape_info.data.front().shape_hessians,
+ shape_info.data.front().fe_degree + 1,
+ shape_info.data.front().n_q_points_1d);
const unsigned int temp_size =
Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
if (temp_size == 0)
{
temp1 = scratch_data;
- temp2 = temp1 +
- std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
- Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+ temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+ shape_info.data.front().fe_degree + 1),
+ Utilities::fixed_power<dim>(
+ shape_info.data.front().n_q_points_1d));
}
else
{
temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
const unsigned int dofs_per_comp =
(type == MatrixFreeFunctions::truncated_tensor) ?
- Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
+ Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
shape_info.dofs_per_component_on_cell;
// expand dof_values to tensor product for truncated tensor products
Number *values_dofs =
{
values_dofs -= dofs_per_comp * n_components;
unsigned int count_p = 0, count_q = 0;
- const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+ const int degree =
+ fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
{
for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
count_q += i * (degree + 1);
}
AssertDimension(count_q,
- Utilities::fixed_power<dim>(shape_info.fe_degree + 1));
+ Utilities::fixed_power<dim>(
+ shape_info.data.front().fe_degree + 1));
}
}
const bool evaluate_gradients,
const bool evaluate_hessians)
{
- AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
- (fe_degree + 2) / 2 * (fe_degree + 1));
+ AssertDimension(
+ shape_info.data.front().shape_gradients_collocation_eo.size(),
+ (fe_degree + 2) / 2 * (fe_degree + 1));
EvaluatorTensorProduct<evaluate_evenodd,
dim,
fe_degree + 1,
Number>
eval(AlignedVector<Number>(),
- shape_info.shape_gradients_collocation_eo,
- shape_info.shape_hessians_collocation_eo);
+ shape_info.data.front().shape_gradients_collocation_eo,
+ shape_info.data.front().shape_hessians_collocation_eo);
constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
for (unsigned int c = 0; c < n_components; c++)
const bool integrate_gradients,
const bool add_into_values_array)
{
- AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
- (fe_degree + 2) / 2 * (fe_degree + 1));
+ AssertDimension(
+ shape_info.data.front().shape_gradients_collocation_eo.size(),
+ (fe_degree + 2) / 2 * (fe_degree + 1));
EvaluatorTensorProduct<evaluate_evenodd,
dim,
fe_degree + 1,
Number>
eval(AlignedVector<Number>(),
- shape_info.shape_gradients_collocation_eo,
- shape_info.shape_hessians_collocation_eo);
+ shape_info.data.front().shape_gradients_collocation_eo,
+ shape_info.data.front().shape_hessians_collocation_eo);
constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
for (unsigned int c = 0; c < n_components; c++)
n_q_points_1d,
1,
Number,
- Number>::do_forward(shape_info.shape_values_eo,
+ Number>::do_forward(shape_info.data.front().shape_values_eo,
values_dofs,
values_quad);
"of lower degree, so the evaluation results would be "
"wrong. Thus, this class does not permit the desired "
"operation."));
- AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
- (n_q_points_1d + 1) / 2 * n_q_points_1d);
+ AssertDimension(
+ shape_info.data.front().shape_gradients_collocation_eo.size(),
+ (n_q_points_1d + 1) / 2 * n_q_points_1d);
constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
for (unsigned int c = 0; c < n_components; c++)
n_q_points_1d,
1,
Number,
- Number>::do_backward(shape_info.shape_values_eo,
+ Number>::do_backward(shape_info.data.front().shape_values_eo,
add_into_values_array,
values_quad,
values_dofs);
{
const AlignedVector<Number> &val1 =
symmetric_evaluate ?
- data.shape_values_eo :
+ data.data.front().shape_values_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_values :
- data.values_within_subface[subface_index % 2]);
+ data.data.front().shape_values :
+ data.data.front().values_within_subface[subface_index % 2]);
const AlignedVector<Number> &val2 =
symmetric_evaluate ?
- data.shape_values_eo :
+ data.data.front().shape_values_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_values :
- data.values_within_subface[subface_index / 2]);
+ data.data.front().shape_values :
+ data.data.front().values_within_subface[subface_index / 2]);
const AlignedVector<Number> &grad1 =
symmetric_evaluate ?
- data.shape_gradients_eo :
+ data.data.front().shape_gradients_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_gradients :
- data.gradients_within_subface[subface_index % 2]);
+ data.data.front().shape_gradients :
+ data.data.front().gradients_within_subface[subface_index % 2]);
const AlignedVector<Number> &grad2 =
symmetric_evaluate ?
- data.shape_gradients_eo :
+ data.data.front().shape_gradients_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_gradients :
- data.gradients_within_subface[subface_index / 2]);
+ data.data.front().shape_gradients :
+ data.data.front().gradients_within_subface[subface_index / 2]);
using Eval =
internal::EvaluatorTensorProduct<symmetric_evaluate ?
Eval eval1(val1,
grad1,
AlignedVector<Number>(),
- data.fe_degree + 1,
- data.n_q_points_1d);
+ data.data.front().fe_degree + 1,
+ data.data.front().n_q_points_1d);
Eval eval2(val2,
grad2,
AlignedVector<Number>(),
- data.fe_degree + 1,
- data.n_q_points_1d);
+ data.data.front().fe_degree + 1,
+ data.data.front().n_q_points_1d);
const unsigned int size_deg =
fe_degree > -1 ?
Utilities::pow(fe_degree + 1, dim - 1) :
- (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
+ (dim > 1 ?
+ Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
+ 1);
const unsigned int n_q_points = fe_degree > -1 ?
Utilities::pow(n_q_points_1d, dim - 1) :
n_q_points_1d,
n_q_points_1d,
Number>
- eval_grad(AlignedVector<Number>(),
- data.shape_gradients_collocation_eo,
- AlignedVector<Number>());
+ eval_grad(
+ AlignedVector<Number>(),
+ data.data.front().shape_gradients_collocation_eo,
+ AlignedVector<Number>());
eval_grad.template gradients<0, true, false>(
values_quad, gradients_quad);
eval_grad.template gradients<1, true, false>(
{
const AlignedVector<Number> &val1 =
symmetric_evaluate ?
- data.shape_values_eo :
+ data.data.front().shape_values_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_values :
- data.values_within_subface[subface_index % 2]);
+ data.data.front().shape_values :
+ data.data.front().values_within_subface[subface_index % 2]);
const AlignedVector<Number> &val2 =
symmetric_evaluate ?
- data.shape_values_eo :
+ data.data.front().shape_values_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_values :
- data.values_within_subface[subface_index / 2]);
+ data.data.front().shape_values :
+ data.data.front().values_within_subface[subface_index / 2]);
const AlignedVector<Number> &grad1 =
symmetric_evaluate ?
- data.shape_gradients_eo :
+ data.data.front().shape_gradients_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_gradients :
- data.gradients_within_subface[subface_index % 2]);
+ data.data.front().shape_gradients :
+ data.data.front().gradients_within_subface[subface_index % 2]);
const AlignedVector<Number> &grad2 =
symmetric_evaluate ?
- data.shape_gradients_eo :
+ data.data.front().shape_gradients_eo :
(subface_index >= GeometryInfo<dim>::max_children_per_cell ?
- data.shape_gradients :
- data.gradients_within_subface[subface_index / 2]);
+ data.data.front().shape_gradients :
+ data.data.front().gradients_within_subface[subface_index / 2]);
using Eval =
internal::EvaluatorTensorProduct<symmetric_evaluate ?
fe_degree + 1,
n_q_points_1d,
Number>;
- Eval eval1(val1, grad1, val1, data.fe_degree + 1, data.n_q_points_1d);
- Eval eval2(val2, grad2, val1, data.fe_degree + 1, data.n_q_points_1d);
+ Eval eval1(val1,
+ grad1,
+ val1,
+ data.data.front().fe_degree + 1,
+ data.data.front().n_q_points_1d);
+ Eval eval2(val2,
+ grad2,
+ val1,
+ data.data.front().fe_degree + 1,
+ data.data.front().n_q_points_1d);
const unsigned int size_deg =
fe_degree > -1 ?
Utilities::pow(fe_degree + 1, dim - 1) :
- (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
+ (dim > 1 ?
+ Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
+ 1);
const unsigned int n_q_points = fe_degree > -1 ?
Utilities::pow(n_q_points_1d, dim - 1) :
n_q_points_1d,
n_q_points_1d,
Number>
- eval_grad(AlignedVector<Number>(),
- data.shape_gradients_collocation_eo,
- AlignedVector<Number>());
+ eval_grad(
+ AlignedVector<Number>(),
+ data.data.front().shape_gradients_collocation_eo,
+ AlignedVector<Number>());
if (integrate_val)
eval_grad.template gradients<1, false, true>(
gradients_quad + n_q_points, values_quad);
fe_degree + 1,
0,
Number>
- evalf(data.shape_data_on_face[face_no % 2],
+ evalf(data.data.front().shape_data_on_face[face_no % 2],
AlignedVector<Number>(),
AlignedVector<Number>(),
- data.fe_degree + 1,
+ data.data.front().fe_degree + 1,
0);
const unsigned int in_stride = do_evaluate ?
fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
numbers::invalid_unsigned_int;
const unsigned int dofs_per_face =
- fe_degree > -1 ? static_dofs_per_face :
- Utilities::pow(data.fe_degree + 1, dim - 1);
+ fe_degree > -1 ?
+ static_dofs_per_face :
+ Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
// we allocate small amounts of data on the stack to signal the compiler
// that this temporary data is only needed for the calculations but the
fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
numbers::invalid_unsigned_int;
const unsigned int dofs_per_face =
- fe_degree > -1 ? static_dofs_per_face :
- Utilities::pow(data.fe_degree + 1, dim - 1);
+ fe_degree > -1 ?
+ static_dofs_per_face :
+ Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
constexpr unsigned int stack_array_size_threshold = 100;
fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
numbers::invalid_unsigned_int;
const unsigned int dofs_per_face =
- fe_degree > -1 ? static_dofs_per_face :
- Utilities::pow(data.fe_degree + 1, dim - 1);
+ fe_degree > -1 ?
+ static_dofs_per_face :
+ Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
constexpr unsigned int stack_array_size_threshold = 100;
// case 1: contiguous and interleaved indices
if (((evaluate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 1 + side];
+ data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
const unsigned int *index_array =
// case 2: contiguous and interleaved indices with fixed stride
else if (((evaluate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 1 + side];
+ data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
// case 3: contiguous and interleaved indices with mixed stride
else if (((evaluate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 1 + side];
+ data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
// case 4: contiguous indices without interleaving
else if (((evaluate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 1 + side];
+ data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
numbers::invalid_unsigned_int;
const unsigned int dofs_per_face =
- fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
- Utilities::pow(data.fe_degree + 1, dim - 1);
+ fe_degree > -1 ?
+ Utilities::pow(fe_degree + 1, dim - 1) :
+ Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
constexpr unsigned int stack_array_size_threshold = 100;
// case 1: contiguous and interleaved indices
if (((integrate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 2 - side];
+ data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
const unsigned int *index_array =
// case 2: contiguous and interleaved indices with fixed stride
else if (((integrate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 2 - side];
+ data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
// case 3: contiguous and interleaved indices with mixed stride
else if (((integrate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 2 - side];
+ data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
// case 4: contiguous indices without interleaving
else if (((integrate_gradients == false &&
- data.nodal_at_cell_boundaries == true) ||
+ data.data.front().nodal_at_cell_boundaries == true) ||
(data.element_type ==
internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
fe_degree > 1)) &&
// right (side==1) are the negative from the value at the left
// (side==0), so we only read out one of them.
const VectorizedArrayType grad_weight =
- data.shape_data_on_face[0][fe_degree + 2 - side];
+ data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
AssertDimension(data.face_to_cell_index_hermite.size(1),
2 * dofs_per_face);
const unsigned int *index_array =
fe_degree + 1,
fe_degree + 1,
Number>
- evaluator(AlignedVector<Number>(),
- AlignedVector<Number>(),
- fe_eval.get_shape_info().inverse_shape_values_eo);
+ evaluator(
+ AlignedVector<Number>(),
+ AlignedVector<Number>(),
+ fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
for (unsigned int d = 0; d < n_components; ++d)
{
const bool evaluate_gradients,
const bool evaluate_hessians)
{
- const unsigned int runtime_degree = shape_info.fe_degree;
+ const unsigned int runtime_degree = shape_info.data.front().fe_degree;
constexpr unsigned int start_n_q_points = degree + 1;
if (runtime_degree == degree)
Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
const bool integrate_gradients,
const bool sum_into_values_array = false)
{
- const int runtime_degree = shape_info.fe_degree;
+ const int runtime_degree = shape_info.data.front().fe_degree;
constexpr unsigned int start_n_q_points = degree + 1;
if (runtime_degree == degree)
Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
const bool evaluate_gradients,
const bool evaluate_hessians)
{
- const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+ const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d;
if (runtime_n_q_points_1d == n_q_points_1d)
{
if (n_q_points_1d == degree + 1 &&
const bool integrate_gradients,
const bool sum_into_values_array)
{
- const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+ const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d;
if (runtime_n_q_points_1d == n_q_points_1d)
{
if (n_q_points_1d == degree + 1 &&
Assert(scratch_data_array != nullptr, ExcInternalError());
const unsigned int tensor_dofs_per_component =
- Utilities::fixed_power<dim>(this->data->fe_degree + 1);
+ Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
const unsigned int dofs_per_component =
this->data->dofs_per_component_on_cell;
const unsigned int n_quadrature_points =
// print error message when the dimensions do not match. Propose a possible
// fix
if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
- static_cast<unsigned int>(fe_degree) != this->data->fe_degree) ||
+ static_cast<unsigned int>(fe_degree) !=
+ this->data->data.front().fe_degree) ||
n_q_points != this->n_quadrature_points)
{
std::string message =
proposed_quad_comp = numbers::invalid_unsigned_int;
if (dof_no != numbers::invalid_unsigned_int)
{
- if (static_cast<unsigned int>(fe_degree) == this->data->fe_degree)
+ if (static_cast<unsigned int>(fe_degree) ==
+ this->data->data.front().fe_degree)
{
proposed_dof_comp = dof_no;
proposed_fe_comp = first_selected_component;
++nf)
if (this->matrix_info
->get_shape_info(no, 0, nf, this->active_fe_index, 0)
+ .data.front()
.fe_degree == static_cast<unsigned int>(fe_degree))
{
proposed_dof_comp = no;
std::pow(1.001 * this->n_quadrature_points, 1. / dim));
message += "Wrong template arguments:\n";
message += " Did you mean FEEvaluation<dim,";
- message += Utilities::int_to_string(this->data->fe_degree) + ",";
+ message +=
+ Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
message += Utilities::int_to_string(proposed_n_q_points_1d);
message += "," + Utilities::int_to_string(n_components);
message += ",Number>(data";
}
message += ")?\n";
std::string correct_pos;
- if (this->data->fe_degree != static_cast<unsigned int>(fe_degree))
+ if (this->data->data.front().fe_degree !=
+ static_cast<unsigned int>(fe_degree))
correct_pos = " ^";
else
correct_pos = " ";
correct_pos += " \n";
message += " " + correct_pos;
- Assert(static_cast<unsigned int>(fe_degree) == this->data->fe_degree &&
+ Assert(static_cast<unsigned int>(fe_degree) ==
+ this->data->data.front().fe_degree &&
n_q_points == this->n_quadrature_points,
ExcMessage(message));
}
AssertIndexRange(q, n_q_points);
const unsigned int n_q_points_1d_actual =
- fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;
+ fe_degree == -1 ? this->data->data.front().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
bool all_nodal = true;
for (unsigned int c = 0; c < di.n_base_elements; ++c)
if (!shape_info(di.global_base_element_offset + c, 0, 0, 0)
+ .data.front()
.nodal_at_cell_boundaries)
all_nodal = false;
if (all_nodal == false)
fe_degree,
n_components,
VectorizedArrayType>::
- apply(fe_eval.get_shape_info().inverse_shape_values_eo,
+ apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
inverse_coefficients,
n_actual_components,
in_array,
n_components,
VectorizedArrayType>::
transform_from_q_points_to_basis(
- fe_eval.get_shape_info().inverse_shape_values_eo,
+ fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
n_actual_components,
in_array,
out_array);
tensor_symmetric_plus_dg0 = 5
};
- /**
- * The class that stores the shape functions, gradients and Hessians
- * evaluated for a tensor product finite element and tensor product
- * quadrature formula on the unit cell. Because of this structure, only
- * one-dimensional data is stored.
- *
- * @ingroup matrixfree
- *
- * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
- */
template <typename Number>
- struct ShapeInfo
+ struct UnivariateShapeData
{
- /**
- * Empty constructor. Does nothing.
- */
- ShapeInfo();
-
- /**
- * Constructor that initializes the data fields using the reinit method.
- */
- template <int dim>
- ShapeInfo(const Quadrature<1> & quad,
- const FiniteElement<dim> &fe,
- const unsigned int base_element = 0);
-
- /**
- * Initializes the data fields. Takes a one-dimensional quadrature
- * formula and a finite element as arguments and evaluates the shape
- * functions, gradients and Hessians on the one-dimensional unit cell.
- * This function assumes that the finite element is derived from a one-
- * dimensional element by a tensor product and that the zeroth shape
- * function in zero evaluates to one.
- */
- template <int dim>
- void
- reinit(const Quadrature<1> & quad,
- const FiniteElement<dim> &fe_dim,
- const unsigned int base_element = 0);
-
/**
* Return the memory consumption of this class in bytes.
*/
* will select the most efficient algorithm based on the given element
* type.
*/
- ElementType element_type;
+ ElementType element_type = ElementType::tensor_general;
/**
* Stores the shape values of the 1D finite element evaluated on all 1D
*/
Quadrature<1> quadrature;
+ /**
+ * Stores the degree of the element.
+ */
+ unsigned int fe_degree = 0;
+
+ /**
+ * Stores the number of quadrature points per dimension.
+ */
+ unsigned int n_q_points_1d = 0;
+
+ /**
+ * Indicates whether the basis functions are nodal in 0 and 1, i.e., the
+ * end points of the unit cell.
+ */
+ bool nodal_at_cell_boundaries = false;
+ };
+
+
+
+ /**
+ * The class that stores the shape functions, gradients and Hessians
+ * evaluated for a tensor product finite element and tensor product
+ * quadrature formula on the unit cell. Because of this structure, only
+ * one-dimensional data is stored.
+ *
+ * @ingroup matrixfree
+ *
+ * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
+ */
+ template <typename Number>
+ struct ShapeInfo
+ {
+ /**
+ * Encodes the type of element detected at construction. FEEvaluation
+ * will select the most efficient algorithm based on the given element
+ * type.
+ */
+ ElementType element_type;
+
+ /**
+ * Empty constructor. Does nothing.
+ */
+ ShapeInfo();
+
+ /**
+ * Constructor that initializes the data fields using the reinit method.
+ */
+ template <int dim>
+ ShapeInfo(const Quadrature<1> & quad,
+ const FiniteElement<dim> &fe,
+ const unsigned int base_element = 0);
+
+ /**
+ * Initializes the data fields. Takes a one-dimensional quadrature
+ * formula and a finite element as arguments and evaluates the shape
+ * functions, gradients and Hessians on the one-dimensional unit cell.
+ * This function assumes that the finite element is derived from a one-
+ * dimensional element by a tensor product and that the zeroth shape
+ * function in zero evaluates to one.
+ */
+ template <int dim>
+ void
+ reinit(const Quadrature<1> & quad,
+ const FiniteElement<dim> &fe_dim,
+ const unsigned int base_element = 0);
+
+ /**
+ * Return data of univariate shape functions which defines the
+ * dimension @p dimension of tensor product shape functions
+ * regarding vector component @p component of the underlying
+ * finite element.
+ */
+ const UnivariateShapeData<Number> &
+ get_shape_data(const unsigned int dimension = 0,
+ const unsigned int component = 0) const;
+
+ /**
+ * Return the memory consumption of this class in bytes.
+ */
+ std::size_t
+ memory_consumption() const;
+
/**
* Renumbering from deal.II's numbering of cell degrees of freedom to
* lexicographic numbering used inside the FEEvaluation schemes of the
std::vector<unsigned int> lexicographic_numbering;
/**
- * Stores the degree of the element.
+ * Stores data of univariate shape functions defining the
+ * underlying tensor product finite element.
*/
- unsigned int fe_degree;
+ std::vector<UnivariateShapeData<Number>> data;
/**
- * Stores the number of quadrature points per dimension.
+ * Grants access to univariate shape function data of given
+ * dimension and vector component. Rows identify dimensions and
+ * columns identify vector components.
+ */
+ dealii::Table<2, UnivariateShapeData<Number> *> data_access;
+
+ /**
+ * Stores the number of space dimensions.
*/
- unsigned int n_q_points_1d;
+ unsigned int n_dimensions;
+
+ /**
+ * Stores the number of vector components of the underlying
+ * vector-valued finite element.
+ */
+ unsigned int n_components;
/**
* Stores the number of quadrature points in @p dim dimensions for a
*/
unsigned int dofs_per_component_on_face;
- /**
- * Indicates whether the basis functions are nodal in 0 and 1, i.e., the
- * end points of the unit cell.
- */
- bool nodal_at_cell_boundaries;
-
/**
* For nodal basis functions with nodes located at the boundary of the
* unit cell, face integrals that involve only the values of the shape
* also fill the shape_???_eo fields.
*/
bool
- check_1d_shapes_symmetric(const unsigned int n_q_points_1d);
+ check_1d_shapes_symmetric(
+ UnivariateShapeData<Number> &univariate_shape_data);
/**
* Check whether symmetric 1D basis functions are such that the shape
* that save some operations in the evaluation.
*/
bool
- check_1d_shapes_collocation() const;
+ check_1d_shapes_collocation(
+ const UnivariateShapeData<Number> &univariate_shape_data) const;
};
const FiniteElement<dim> &fe_in,
const unsigned int base_element_number)
: element_type(tensor_general)
- , fe_degree(0)
- , n_q_points_1d(0)
+ , n_dimensions(0)
+ , n_components(0)
, n_q_points(0)
, dofs_per_component_on_cell(0)
, n_q_points_face(0)
, dofs_per_component_on_face(0)
- , nodal_at_cell_boundaries(false)
{
reinit(quad, fe_in, base_element_number);
}
+ template <typename Number>
+ inline const UnivariateShapeData<Number> &
+ ShapeInfo<Number>::get_shape_data(const unsigned int dimension,
+ const unsigned int component) const
+ {
+ AssertDimension(n_dimensions, data_access.size(0));
+ AssertDimension(n_components, data_access.size(1));
+ AssertIndexRange(dimension, n_dimensions);
+ AssertIndexRange(component, n_components);
+ return *(data_access(dimension, component));
+ }
+
} // end of namespace MatrixFreeFunctions
} // end of namespace internal
template <typename Number>
ShapeInfo<Number>::ShapeInfo()
: element_type(tensor_general)
- , fe_degree(numbers::invalid_unsigned_int)
- , n_q_points_1d(0)
+ , n_dimensions(0)
+ , n_components(0)
, n_q_points(0)
, dofs_per_component_on_cell(0)
, n_q_points_face(0)
, dofs_per_component_on_face(0)
- , nodal_at_cell_boundaries(false)
{}
const FiniteElement<dim> &fe_in,
const unsigned int base_element_number)
{
- this->quadrature = quad;
const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
+ n_dimensions = dim;
+ n_components = fe_in.n_components();
Assert(fe->n_components() == 1,
ExcMessage("FEEvaluation only works for scalar finite elements."));
- fe_degree = fe->degree;
- n_q_points_1d = quad.size();
-
+ // assuming isotropy of dimensions and components
+ data.resize(1);
+ UnivariateShapeData<Number> *univariate_shape_data = &data.front();
+ data_access.reinit(n_dimensions, n_components);
+ data_access.fill(univariate_shape_data);
+ univariate_shape_data->quadrature = quad;
+ univariate_shape_data->fe_degree = fe->degree;
+ univariate_shape_data->n_q_points_1d = quad.size();
+
+ // grant write access to common univariate shape data
+ auto &shape_values = univariate_shape_data->shape_values;
+ auto &shape_gradients = univariate_shape_data->shape_gradients;
+ auto &shape_hessians = univariate_shape_data->shape_hessians;
+ auto &shape_gradients_collocation =
+ univariate_shape_data->shape_gradients_collocation;
+ auto &shape_hessians_collocation =
+ univariate_shape_data->shape_hessians_collocation;
+ auto &inverse_shape_values = univariate_shape_data->inverse_shape_values;
+ auto &shape_data_on_face = univariate_shape_data->shape_data_on_face;
+ auto &values_within_subface =
+ univariate_shape_data->values_within_subface;
+ auto &gradients_within_subface =
+ univariate_shape_data->gradients_within_subface;
+ auto &hessians_within_subface =
+ univariate_shape_data->hessians_within_subface;
+ auto &nodal_at_cell_boundaries =
+ univariate_shape_data->nodal_at_cell_boundaries;
+
+ const unsigned int fe_degree = fe->degree;
+ const unsigned int n_q_points_1d = quad.size();
const unsigned int n_dofs_1d = std::min(fe->dofs_per_cell, fe_degree + 1);
// renumber (this is necessary for FE_Q, for example, since there the
dim > 1 ? Utilities::fixed_power<dim - 1>(fe_degree + 1) : 1;
const unsigned int array_size = n_dofs_1d * n_q_points_1d;
- this->shape_gradients.resize_fast(array_size);
- this->shape_values.resize_fast(array_size);
- this->shape_hessians.resize_fast(array_size);
-
- this->shape_data_on_face[0].resize(3 * n_dofs_1d);
- this->shape_data_on_face[1].resize(3 * n_dofs_1d);
- this->values_within_subface[0].resize(array_size);
- this->values_within_subface[1].resize(array_size);
- this->gradients_within_subface[0].resize(array_size);
- this->gradients_within_subface[1].resize(array_size);
- this->hessians_within_subface[0].resize(array_size);
- this->hessians_within_subface[1].resize(array_size);
+ shape_gradients.resize_fast(array_size);
+ shape_values.resize_fast(array_size);
+ shape_hessians.resize_fast(array_size);
+
+ shape_data_on_face[0].resize(3 * n_dofs_1d);
+ shape_data_on_face[1].resize(3 * n_dofs_1d);
+ values_within_subface[0].resize(array_size);
+ values_within_subface[1].resize(array_size);
+ gradients_within_subface[0].resize(array_size);
+ gradients_within_subface[1].resize(array_size);
+ hessians_within_subface[0].resize(array_size);
+ hessians_within_subface[1].resize(array_size);
for (unsigned int i = 0; i < n_dofs_1d; ++i)
{
}
// evaluate basis functions on the 1D faces, i.e., in zero and one
- Point<dim> q_point = unit_point;
- q_point[0] = 0;
- this->shape_data_on_face[0][i] = fe->shape_value(my_i, q_point);
- this->shape_data_on_face[0][i + n_dofs_1d] =
+ Point<dim> q_point = unit_point;
+ q_point[0] = 0;
+ shape_data_on_face[0][i] = fe->shape_value(my_i, q_point);
+ shape_data_on_face[0][i + n_dofs_1d] =
fe->shape_grad(my_i, q_point)[0];
- this->shape_data_on_face[0][i + 2 * n_dofs_1d] =
+ shape_data_on_face[0][i + 2 * n_dofs_1d] =
fe->shape_grad_grad(my_i, q_point)[0][0];
- q_point[0] = 1;
- this->shape_data_on_face[1][i] = fe->shape_value(my_i, q_point);
- this->shape_data_on_face[1][i + n_dofs_1d] =
+ q_point[0] = 1;
+ shape_data_on_face[1][i] = fe->shape_value(my_i, q_point);
+ shape_data_on_face[1][i + n_dofs_1d] =
fe->shape_grad(my_i, q_point)[0];
- this->shape_data_on_face[1][i + 2 * n_dofs_1d] =
+ shape_data_on_face[1][i + 2 * n_dofs_1d] =
fe->shape_grad_grad(my_i, q_point)[0][0];
}
}
if (element_type == tensor_general &&
- check_1d_shapes_symmetric(n_q_points_1d))
+ check_1d_shapes_symmetric(*univariate_shape_data))
{
- if (check_1d_shapes_collocation())
+ if (check_1d_shapes_collocation(*univariate_shape_data))
element_type = tensor_symmetric_collocation;
else
element_type = tensor_symmetric;
}
}
else if (element_type == tensor_symmetric_plus_dg0)
- check_1d_shapes_symmetric(n_q_points_1d);
+ check_1d_shapes_symmetric(*univariate_shape_data);
nodal_at_cell_boundaries = true;
for (unsigned int i = 1; i < n_dofs_1d; ++i)
}
}
}
+
+ // TODO !!!
+ univariate_shape_data->element_type = this->element_type;
}
template <typename Number>
bool
ShapeInfo<Number>::check_1d_shapes_symmetric(
- const unsigned int n_q_points_1d)
+ UnivariateShapeData<Number> &univariate_shape_data)
{
if (dofs_per_component_on_cell == 0)
return false;
+ const auto n_q_points_1d = univariate_shape_data.n_q_points_1d;
+ const auto fe_degree = univariate_shape_data.fe_degree;
+ auto & shape_values = univariate_shape_data.shape_values;
+ auto & shape_gradients = univariate_shape_data.shape_gradients;
+ auto & shape_hessians = univariate_shape_data.shape_hessians;
+ auto & shape_gradients_collocation =
+ univariate_shape_data.shape_gradients_collocation;
+ auto &shape_hessians_collocation =
+ univariate_shape_data.shape_hessians_collocation;
+ auto &shape_values_eo = univariate_shape_data.shape_values_eo;
+ auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo;
+ auto &shape_hessians_eo = univariate_shape_data.shape_hessians_eo;
+ auto &shape_gradients_collocation_eo =
+ univariate_shape_data.shape_gradients_collocation_eo;
+ auto &shape_hessians_collocation_eo =
+ univariate_shape_data.shape_hessians_collocation_eo;
+ auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
+ auto &inverse_shape_values_eo =
+ univariate_shape_data.inverse_shape_values_eo;
+
const double zero_tol =
std::is_same<Number, double>::value == true ? 1e-12 : 1e-7;
// symmetry for values
template <typename Number>
bool
- ShapeInfo<Number>::check_1d_shapes_collocation() const
+ ShapeInfo<Number>::check_1d_shapes_collocation(
+ const UnivariateShapeData<Number> &univariate_shape_data) const
{
if (dofs_per_component_on_cell != n_q_points)
return false;
+ const auto fe_degree = univariate_shape_data.fe_degree;
+ auto & shape_values = univariate_shape_data.shape_values;
+
const double zero_tol =
std::is_same<Number, double>::value == true ? 1e-12 : 1e-7;
// check: identity operation for shape values
template <typename Number>
std::size_t
ShapeInfo<Number>::memory_consumption() const
+ {
+ std::size_t memory = sizeof(*this);
+ for (const auto &univariate_shape_data : data)
+ memory += univariate_shape_data.memory_consumption();
+ return memory;
+ }
+
+ template <typename Number>
+ std::size_t
+ UnivariateShapeData<Number>::memory_consumption() const
{
std::size_t memory = sizeof(*this);
memory += MemoryConsumption::memory_consumption(shape_values);
deallog << "Testing " << fe.get_name() << " with " << quadrature_name << "("
<< quad.size() << ")" << std::endl;
deallog << "shape values: " << std::endl;
- for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+ const auto &univariate_shape_data = shape_info.get_shape_data(0, 0);
+ for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
{
for (unsigned int q = 0; q < quad.size(); ++q)
- deallog << std::setw(15) << shape_info.shape_values[i * quad.size() + q]
+ deallog << std::setw(15)
+ << univariate_shape_data.shape_values[i * quad.size() + q]
<< " ";
deallog << std::endl;
}
deallog << "inverse shape values: " << std::endl;
- for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+ for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
{
for (unsigned int q = 0; q < quad.size(); ++q)
- deallog << std::setw(15)
- << shape_info.inverse_shape_values[i * quad.size() + q] << " ";
+ deallog
+ << std::setw(15)
+ << univariate_shape_data.inverse_shape_values[i * quad.size() + q]
+ << " ";
deallog << std::endl;
}
deallog << "inverse shapes' * shapes: " << std::endl;
for (unsigned int j = 0; j < quad.size(); ++j)
{
double sum = 0;
- for (unsigned int k = 0; k < shape_info.fe_degree + 1; ++k)
- sum += shape_info.inverse_shape_values[k * quad.size() + i] *
- shape_info.shape_values[k * quad.size() + j];
+ for (unsigned int k = 0; k < univariate_shape_data.fe_degree + 1; ++k)
+ sum +=
+ univariate_shape_data.inverse_shape_values[k * quad.size() + i] *
+ univariate_shape_data.shape_values[k * quad.size() + j];
deallog << std::setw(15) << sum << " ";
}
deallog << std::endl;
}
deallog << "inverse shapes * shapes': " << std::endl;
- for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+ for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
{
- for (unsigned int j = 0; j < shape_info.fe_degree + 1; ++j)
+ for (unsigned int j = 0; j < univariate_shape_data.fe_degree + 1; ++j)
{
double sum = 0;
for (unsigned int k = 0; k < quad.size(); ++k)
- sum += shape_info.inverse_shape_values[i * quad.size() + k] *
- shape_info.shape_values[j * quad.size() + k];
+ sum +=
+ univariate_shape_data.inverse_shape_values[i * quad.size() + k] *
+ univariate_shape_data.shape_values[j * quad.size() + k];
deallog << std::setw(15) << sum << " ";
}
deallog << std::endl;