Threads::Mutex::ScopedLock lock(this->mutex);
const unsigned int n = Utilities::fixed_power<dim>(size > 0 ? size : eigenvalues[0].size());
tmp_array.resize_fast(n*2);
- constexpr int kernel_size = size > 0 ? size-1 : -1;
- internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size+1,Number>
+ constexpr int kernel_size = size > 0 ? size : 0;
+ internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size,Number>
eval(AlignedVector<Number> {}, AlignedVector<Number> {},
- AlignedVector<Number> {}, mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows());
+ AlignedVector<Number> {}, mass_matrix[0].n_rows(), mass_matrix[0].n_rows());
Number *t = tmp_array.begin();
const Number *src = src_view.begin();
Number *dst = &(dst_view[0]);
Threads::Mutex::ScopedLock lock(this->mutex);
const unsigned int n = size > 0 ? size : eigenvalues[0].size();
tmp_array.resize_fast (Utilities::fixed_power<dim>(n));
- constexpr int kernel_size = size > 0 ? size-1 : -1;
- internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size+1,Number>
+ constexpr int kernel_size = size > 0 ? size : 0;
+ internal::EvaluatorTensorProduct<internal::evaluate_general,dim,kernel_size,kernel_size,Number>
eval(AlignedVector<Number>(), AlignedVector<Number>(),
- AlignedVector<Number>(), mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows());
+ AlignedVector<Number>(), mass_matrix[0].n_rows(), mass_matrix[0].n_rows());
Number *t = tmp_array.begin();
const Number *src = src_view.data();
Number *dst = &(dst_view[0]);
struct FEEvaluationImpl
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_values,
- const bool evaluate_gradients);
+ void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const bool add_into_values_array);
};
+
template <MatrixFreeFunctions::ElementType type, int dim, int fe_degree,
int n_q_points_1d, int n_components, typename Number>
inline
void
FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
const EvaluatorVariant variant =
EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
- typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
- VectorizedArray<Number> > Eval;
+ typedef EvaluatorTensorProduct<variant, dim, fe_degree+1, n_q_points_1d,
+ Number> Eval;
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,
+ shape_info.fe_degree+1,
shape_info.n_q_points_1d);
- const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
- : (Eval::dofs_per_cell > Eval::n_q_points ?
- Eval::dofs_per_cell : Eval::n_q_points);
- VectorizedArray<Number> *temp1;
- VectorizedArray<Number> *temp2;
+ const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? 0
+ : (Eval::n_rows_of_product > Eval::n_columns_of_product ?
+ Eval::n_rows_of_product : Eval::n_columns_of_product);
+ Number *temp1;
+ Number *temp2;
if (temp_size == 0)
{
temp1 = scratch_data;
temp2 = temp1 + temp_size;
}
- VectorizedArray<Number> **values_dofs = values_dofs_actual;
- VectorizedArray<Number> *expanded_dof_values[n_components];
+ const unsigned int n_q_points = 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) : shape_info.dofs_per_component_on_cell;
+ const Number *values_dofs = values_dofs_actual;
if (type == MatrixFreeFunctions::truncated_tensor)
{
- values_dofs = expanded_dof_values;
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell,
- shape_info.n_q_points)) +
- c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
+ 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;
unsigned int count_p = 0, count_q = 0;
for (int i=0; i<(dim>2?degree+1:1); ++i)
{
for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = values_dofs_actual[c][count_p];
+ values_dofs_tmp[c*dofs_per_comp+count_q] = values_dofs_actual[c*shape_info.dofs_per_component_on_cell+count_p];
for (int k=degree+1-j-i; k<degree+1; ++k, ++count_q)
for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+ values_dofs_tmp[c*dofs_per_comp+count_q] = Number();
}
for (int j=degree+1-i; j<degree+1; ++j)
for (int k=0; k<degree+1; ++k, ++count_q)
for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c][count_q] = VectorizedArray<Number>();
+ values_dofs_tmp[c*dofs_per_comp+count_q] = Number();
}
- AssertDimension(count_q, Utilities::fixed_power<dim>(shape_info.fe_degree+1));
+ AssertDimension(count_q, dofs_per_comp);
+ values_dofs = values_dofs_tmp;
}
- // These avoid compiler warnings; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
- const unsigned int d3 = dim>2?3:0;
- const unsigned int d4 = dim>2?4:0;
- const unsigned int d5 = dim>2?5:0;
-
switch (dim)
{
case 1:
for (unsigned int c=0; c<n_components; c++)
{
if (evaluate_values == true)
- eval.template values<0,true,false> (values_dofs[c], values_quad[c]);
+ eval.template values<0,true,false> (values_dofs, values_quad);
if (evaluate_gradients == true)
- eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
+ eval.template gradients<0,true,false>(values_dofs, gradients_quad);
if (evaluate_hessians == true)
- eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
+ eval.template hessians<0,true,false> (values_dofs, hessians_quad);
+
+ // advance the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += n_q_points;
+ hessians_quad += n_q_points;
}
break;
// grad x
if (evaluate_gradients == true)
{
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, gradients_quad[c][0]);
+ eval.template gradients<0,true,false> (values_dofs, temp1);
+ eval.template values<1,true,false> (temp1, gradients_quad);
}
if (evaluate_hessians == true)
{
// grad xy
if (evaluate_gradients == false)
- eval.template gradients<0,true,false>(values_dofs[c], temp1);
- eval.template gradients<1,true,false> (temp1, hessians_quad[c][d1+d1]);
+ eval.template gradients<0,true,false>(values_dofs, temp1);
+ eval.template gradients<1,true,false> (temp1, hessians_quad+2*n_q_points);
// grad xx
- eval.template hessians<0,true,false>(values_dofs[c], temp1);
- eval.template values<1,true,false> (temp1, hessians_quad[c][0]);
+ eval.template hessians<0,true,false>(values_dofs, temp1);
+ eval.template values<1,true,false> (temp1, hessians_quad);
}
// grad y
- eval.template values<0,true,false> (values_dofs[c], temp1);
+ eval.template values<0,true,false> (values_dofs, temp1);
if (evaluate_gradients == true)
- eval.template gradients<1,true,false> (temp1, gradients_quad[c][d1]);
+ eval.template gradients<1,true,false> (temp1, gradients_quad+n_q_points);
// grad yy
if (evaluate_hessians == true)
- eval.template hessians<1,true,false> (temp1, hessians_quad[c][d1]);
+ eval.template hessians<1,true,false> (temp1, hessians_quad+n_q_points);
// val: can use values applied in x
if (evaluate_values == true)
- eval.template values<1,true,false> (temp1, values_quad[c]);
+ eval.template values<1,true,false> (temp1, values_quad);
+
+ // advance to the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += 2*n_q_points;
+ hessians_quad += 3*n_q_points;
}
break;
if (evaluate_gradients == true)
{
// grad x
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
+ eval.template gradients<0,true,false> (values_dofs, temp1);
eval.template values<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, gradients_quad[c][0]);
+ eval.template values<2,true,false> (temp2, gradients_quad);
}
if (evaluate_hessians == true)
// grad xz
if (evaluate_gradients == false)
{
- eval.template gradients<0,true,false> (values_dofs[c], temp1);
+ eval.template gradients<0,true,false> (values_dofs, temp1);
eval.template values<1,true,false> (temp1, temp2);
}
- eval.template gradients<2,true,false> (temp2, hessians_quad[c][d4]);
+ eval.template gradients<2,true,false> (temp2, hessians_quad+4*n_q_points);
// grad xy
eval.template gradients<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][d3]);
+ eval.template values<2,true,false> (temp2, hessians_quad+3*n_q_points);
// grad xx
- eval.template hessians<0,true,false>(values_dofs[c], temp1);
+ eval.template hessians<0,true,false>(values_dofs, temp1);
eval.template values<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][0]);
+ eval.template values<2,true,false> (temp2, hessians_quad);
}
// grad y
- eval.template values<0,true,false> (values_dofs[c], temp1);
+ eval.template values<0,true,false> (values_dofs, temp1);
if (evaluate_gradients == true)
{
eval.template gradients<1,true,false>(temp1, temp2);
- eval.template values<2,true,false> (temp2, gradients_quad[c][d1]);
+ eval.template values<2,true,false> (temp2, gradients_quad+n_q_points);
}
if (evaluate_hessians == true)
// grad yz
if (evaluate_gradients == false)
eval.template gradients<1,true,false>(temp1, temp2);
- eval.template gradients<2,true,false> (temp2, hessians_quad[c][d5]);
+ eval.template gradients<2,true,false> (temp2, hessians_quad+5*n_q_points);
// grad yy
eval.template hessians<1,true,false> (temp1, temp2);
- eval.template values<2,true,false> (temp2, hessians_quad[c][d1]);
+ eval.template values<2,true,false> (temp2, hessians_quad+n_q_points);
}
// grad z: can use the values applied in x direction stored in temp1
eval.template values<1,true,false> (temp1, temp2);
if (evaluate_gradients == true)
- eval.template gradients<2,true,false> (temp2, gradients_quad[c][d2]);
+ eval.template gradients<2,true,false> (temp2, gradients_quad+2*n_q_points);
// grad zz: can use the values applied in x and y direction stored
// in temp2
if (evaluate_hessians == true)
- eval.template hessians<2,true,false>(temp2, hessians_quad[c][d2]);
+ eval.template hessians<2,true,false>(temp2, hessians_quad+2*n_q_points);
// val: can use the values applied in x & y direction stored in temp2
if (evaluate_values == true)
- eval.template values<2,true,false> (temp2, values_quad[c]);
+ eval.template values<2,true,false> (temp2, values_quad);
+
+ // advance to the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += 3*n_q_points;
+ hessians_quad += 6*n_q_points;
}
break;
// case additional dof for FE_Q_DG0: add values; gradients and second
// derivatives evaluate to zero
if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0 && evaluate_values)
- for (unsigned int c=0; c<n_components; ++c)
- for (unsigned int q=0; q<shape_info.n_q_points; ++q)
- values_quad[c][q] += values_dofs[c][shape_info.dofs_per_component_on_cell-1];
+ {
+ values_quad -= n_components*n_q_points;
+ values_dofs -= n_components*dofs_per_comp;
+ for (unsigned int c=0; c<n_components; ++c)
+ for (unsigned int q=0; q<shape_info.n_q_points; ++q)
+ values_quad[c*shape_info.n_q_points+q] +=
+ values_dofs[(c+1)*shape_info.dofs_per_component_on_cell-1];
+ }
}
inline
void
FEEvaluationImpl<type,dim,fe_degree,n_q_points_1d,n_components,Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
- const bool integrate_gradients)
+ const bool integrate_gradients,
+ const bool add_into_values_array)
{
const EvaluatorVariant variant =
EvaluatorSelector<type,(fe_degree+n_q_points_1d>4)>::variant;
- typedef EvaluatorTensorProduct<variant, dim, fe_degree, n_q_points_1d,
- VectorizedArray<Number> > Eval;
+ typedef EvaluatorTensorProduct<variant, dim, fe_degree+1, n_q_points_1d,
+ Number> Eval;
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,
+ shape_info.fe_degree+1,
shape_info.n_q_points_1d);
- const unsigned int temp_size = Eval::dofs_per_cell == numbers::invalid_unsigned_int ? 0
- : (Eval::dofs_per_cell > Eval::n_q_points ?
- Eval::dofs_per_cell : Eval::n_q_points);
- VectorizedArray<Number> *temp1;
- VectorizedArray<Number> *temp2;
+ const unsigned int temp_size = Eval::n_rows_of_product == numbers::invalid_unsigned_int ? 0
+ : (Eval::n_rows_of_product > Eval::n_columns_of_product ?
+ Eval::n_rows_of_product : Eval::n_columns_of_product);
+ Number *temp1;
+ Number *temp2;
if (temp_size == 0)
{
temp1 = scratch_data;
temp2 = temp1 + temp_size;
}
+ const unsigned int n_q_points = 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) : shape_info.dofs_per_component_on_cell;
// expand dof_values to tensor product for truncated tensor products
- VectorizedArray<Number> **values_dofs = values_dofs_actual;
- VectorizedArray<Number> *expanded_dof_values[n_components];
- if (type == MatrixFreeFunctions::truncated_tensor)
- {
- values_dofs = expanded_dof_values;
- for (unsigned int c=0; c<n_components; ++c)
- expanded_dof_values[c] = scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell,
- shape_info.n_q_points)) +
- c*Utilities::fixed_power<dim>(shape_info.fe_degree+1);
- }
-
- // These avoid compiler warnings; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
+ Number *values_dofs = (type == MatrixFreeFunctions::truncated_tensor) ?
+ scratch_data+2*(std::max(shape_info.dofs_per_component_on_cell,
+ shape_info.n_q_points)) :
+ values_dofs_actual;
switch (dim)
{
for (unsigned int c=0; c<n_components; c++)
{
if (integrate_values == true)
- eval.template values<0,false,false> (values_quad[c], values_dofs[c]);
+ {
+ if (add_into_values_array == false)
+ eval.template values<0,false,false> (values_quad, values_dofs);
+ else
+ eval.template values<0,false,true> (values_quad, values_dofs);
+ }
if (integrate_gradients == true)
{
- if (integrate_values == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], values_dofs[c]);
+ if (integrate_values == true || add_into_values_array == true)
+ eval.template gradients<0,false,true> (gradients_quad, values_dofs);
else
- eval.template gradients<0,false,false> (gradients_quad[c][0], values_dofs[c]);
+ eval.template gradients<0,false,false> (gradients_quad, values_dofs);
}
+
+ // advance to the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += n_q_points;
}
break;
case 2:
for (unsigned int c=0; c<n_components; c++)
{
- if (integrate_values == true)
+ if (integrate_values == true &&
+ integrate_gradients == false)
{
- // val
- eval.template values<0,false,false> (values_quad[c], temp1);
- //grad x
- if (integrate_gradients == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
- eval.template values<1,false,false>(temp1, values_dofs[c]);
+ eval.template values<1,false,false> (values_quad, temp1);
+ if (add_into_values_array == false)
+ eval.template values<0,false,false>(temp1, values_dofs);
+ else
+ eval.template values<0,false,true>(temp1, values_dofs);
}
if (integrate_gradients == true)
{
- // grad y
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- if (integrate_values == false)
- {
- eval.template gradients<1,false,false>(temp1, values_dofs[c]);
- //grad x
- eval.template gradients<0,false,false> (gradients_quad[c][0], temp1);
- eval.template values<1,false,true> (temp1, values_dofs[c]);
- }
+ eval.template gradients<1,false,false> (gradients_quad+n_q_points, temp1);
+ if (integrate_values)
+ eval.template values<1,false,true> (values_quad, temp1);
+ if (add_into_values_array == false)
+ eval.template values<0,false,false>(temp1, values_dofs);
else
- eval.template gradients<1,false,true>(temp1, values_dofs[c]);
+ eval.template values<0,false,true>(temp1, values_dofs);
+ eval.template values<1,false,false> (gradients_quad, temp1);
+ eval.template gradients<0,false,true> (temp1, values_dofs);
}
+
+ // advance to the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += 2*n_q_points;
}
break;
case 3:
for (unsigned int c=0; c<n_components; c++)
{
- if (integrate_values == true)
- {
- // val
- eval.template values<0,false,false> (values_quad[c], temp1);
- //grad x: can sum to temporary value in temp1
- if (integrate_gradients == true)
- eval.template gradients<0,false,true> (gradients_quad[c][0], temp1);
- eval.template values<1,false,false>(temp1, temp2);
- if (integrate_gradients == true)
- {
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- eval.template gradients<1,false,true>(temp1, temp2);
- }
- eval.template values<2,false,false> (temp2, values_dofs[c]);
- }
- else if (integrate_gradients == true)
+ if (integrate_values == true &&
+ integrate_gradients == false)
{
- eval.template gradients<0,false,false>(gradients_quad[c][0], temp1);
- eval.template values<1,false,false> (temp1, temp2);
- eval.template values<0,false,false> (gradients_quad[c][d1], temp1);
- eval.template gradients<1,false,true>(temp1, temp2);
- eval.template values<2,false,false> (temp2, values_dofs[c]);
+ eval.template values<2,false,false> (values_quad, temp1);
+ eval.template values<1,false,false> (temp1, temp2);
+ if (add_into_values_array == false)
+ eval.template values<0,false,false>(temp2, values_dofs);
+ else
+ eval.template values<0,false,true> (temp2, values_dofs);
}
if (integrate_gradients == true)
{
- // grad z: can sum to temporary x and y value in output
- eval.template values<0,false,false> (gradients_quad[c][d2], temp1);
- eval.template values<1,false,false> (temp1, temp2);
- eval.template gradients<2,false,true> (temp2, values_dofs[c]);
+ eval.template gradients<2,false,false>(gradients_quad+2*n_q_points, temp1);
+ if (integrate_values)
+ eval.template values<2,false,true> (values_quad, temp1);
+ eval.template values<1,false,false> (temp1, temp2);
+ eval.template values<2,false,false> (gradients_quad+n_q_points, temp1);
+ eval.template gradients<1,false,true> (temp1, temp2);
+ if (add_into_values_array == false)
+ eval.template values<0,false,false> (temp2, values_dofs);
+ else
+ eval.template values<0,false,true> (temp2, values_dofs);
+ eval.template values<2,false,false> (gradients_quad, temp1);
+ eval.template values<1,false,false> (temp1, temp2);
+ eval.template gradients<0,false,true> (temp2, values_dofs);
}
+
+ // advance to the next component in 1D array
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
+ gradients_quad += 3*n_q_points;
}
break;
// case FE_Q_DG0: add values, gradients and second derivatives are zero
if (type == MatrixFreeFunctions::tensor_symmetric_plus_dg0)
{
+ values_dofs -= n_components * dofs_per_comp - shape_info.dofs_per_component_on_cell + 1;
+ values_quad -= n_components * n_q_points;
if (integrate_values)
for (unsigned int c=0; c<n_components; ++c)
{
- values_dofs[c][shape_info.dofs_per_component_on_cell-1] = values_quad[c][0];
+ values_dofs[0] = values_quad[0];
for (unsigned int q=1; q<shape_info.n_q_points; ++q)
- values_dofs[c][shape_info.dofs_per_component_on_cell-1] += values_quad[c][q];
+ values_dofs[0] += values_quad[q];
+ values_dofs += dofs_per_comp;
+ values_quad += n_q_points;
}
else
- for (unsigned int c=0; c<n_components; ++c)
- values_dofs[c][shape_info.dofs_per_component_on_cell-1] = VectorizedArray<Number>();
+ {
+ for (unsigned int c=0; c<n_components; ++c)
+ values_dofs[c*shape_info.dofs_per_component_on_cell] = Number();
+ values_dofs += n_components*shape_info.dofs_per_component_on_cell;
+ }
}
if (type == MatrixFreeFunctions::truncated_tensor)
{
+ 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;
for (int i=0; i<(dim>2?degree+1:1); ++i)
for (int k=0; k<degree+1-j-i; ++k, ++count_p, ++count_q)
{
for (unsigned int c=0; c<n_components; ++c)
- values_dofs_actual[c][count_p] = expanded_dof_values[c][count_q];
+ values_dofs_actual[c*shape_info.dofs_per_component_on_cell+count_p] = values_dofs[c*dofs_per_comp+count_q];
}
count_q += j+i;
}
/**
- * This struct performs the evaluation of function values, gradients and
- * Hessians for tensor-product finite elements. This a specialization for
- * symmetric basis functions about the mid point 0.5 of the unit interval
- * with the same number of quadrature points as degrees of freedom. In that
- * case, we can first transform the basis to one that has the nodal points
- * in the quadrature points (i.e., the collocation space) and then perform
- * the evaluation of the first and second derivatives in this transformed
- * space, using the identity operation for the shape values.
+ * This struct implements the change between two different bases. This is an
+ * ingredient in the FEEvaluationImplTransformToCollocation class where we
+ * first transform to the appropriate basis where we can compute the
+ * derivative through collocation techniques.
+ *
+ * This class allows for dimension-independent application of the operation,
+ * implemented by template recursion. It has been tested up to 6D.
*
* @author Katharina Kormann, Martin Kronbichler, 2017
*/
+ template <EvaluatorVariant variant, int dim, int basis_size_1, int basis_size_2, int n_components,
+ typename Number, typename Number2>
+ struct FEEvaluationImplBasisChange
+ {
+ static_assert(basis_size_1 == 0 || basis_size_1 <= basis_size_2,
+ "The second dimension must not be smaller than the first");
+
+ /**
+ * This applies the transformation that contracts over the rows of the
+ * coefficient array, generating values along the columns of the
+ * coefficient array.
+ *
+ * @param transformation_matrix The coefficient matrix handed in as a
+ * vector, using @p basis_size_1 rows and @p basis_size_2
+ * columns if interpreted as a matrix.
+ * @param values_in The array of the input of size basis_size_1^dim. It
+ * may alias with values_out
+ * @param values_out The array of size basis_size_2^dim where the results
+ * of the transformation are stored. It may alias with
+ * the values_in array.
+ * @param basis_size_1_variable In case the template argument basis_size_1 is
+ * zero, the size of the first basis can alternatively be
+ * passed in as a run time argument. The template
+ * argument takes precedence in case it is nonzero
+ * for efficiency reasons.
+ * @param basis_size_2_variable In case the template argument basis_size_1 is
+ * zero, the size of the second basis can alternatively be
+ * passed in as a run time argument.
+ */
+#ifndef DEBUG
+ DEAL_II_ALWAYS_INLINE
+#endif
+ static void do_forward (const AlignedVector<Number2> &transformation_matrix,
+ const Number *values_in,
+ Number *values_out,
+ const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
+ const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
+ {
+ Assert(basis_size_1 != 0 ||
+ basis_size_1_variable <= basis_size_2_variable,
+ ExcMessage("The second dimension must not be smaller than the first"));
+
+ // we do recursion until dim==1 or dim==2 and we have
+ // basis_size_1==basis_size_2. The latter optimization increases
+ // optimization possibilities for the compiler but does only work for
+ // aliased pointers if the sizes are equal.
+ constexpr int next_dim = (dim > 2 || ((basis_size_1 == 0 || basis_size_2>basis_size_1)
+ && dim>1)) ? dim-1 : dim;
+
+ EvaluatorTensorProduct<variant, dim, basis_size_1, (basis_size_1==0 ? 0 : basis_size_2),
+ Number,Number2> eval_val (transformation_matrix,
+ AlignedVector<Number2>(),
+ AlignedVector<Number2>(),
+ basis_size_1_variable,
+ basis_size_2_variable);
+ const unsigned int np_1 = basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
+ const unsigned int np_2 = basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
+ Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
+ ExcMessage("Cannot transform with 0-point basis"));
+ Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
+ ExcMessage("Cannot transform with 0-point basis"));
+
+ // run loop backwards to ensure correctness if values_in aliases with
+ // values_out in case with basis_size_1 < basis_size_2
+ values_in = values_in + n_components*Utilities::fixed_power<dim>(np_1);
+ values_out = values_out + n_components*Utilities::fixed_power<dim>(np_2);
+ for (unsigned int c=n_components; c!=0; --c)
+ {
+ values_in -= Utilities::fixed_power<dim>(np_1);
+ values_out -= Utilities::fixed_power<dim>(np_2);
+ if (next_dim < dim)
+ for (unsigned int q=np_1; q!=0; --q)
+ FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,1,Number,Number2>
+ ::do_forward(transformation_matrix,
+ values_in + (q-1)*Utilities::fixed_power<next_dim>(np_1),
+ values_out + (q-1)*Utilities::fixed_power<next_dim>(np_2),
+ basis_size_1_variable,
+ basis_size_2_variable);
+
+ // the recursion stops if dim==1 or if dim==2 and
+ // basis_size_1==basis_size_2 (the latter is used because the
+ // compiler generates nicer code)
+ if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
+ {
+ eval_val.template values<0,true,false>(values_in, values_out);
+ eval_val.template values<1,true,false>(values_out, values_out);
+ }
+ else if (dim==1)
+ eval_val.template values<dim-1,true,false>(values_in, values_out);
+ else
+ eval_val.template values<dim-1,true,false>(values_out, values_out);
+ }
+ }
+
+ /**
+ * This applies the transformation that contracts over the columns of the
+ * coefficient array, generating values along the rows of the coefficient
+ * array.
+ *
+ * @param transformation_matrix The coefficient matrix handed in as a
+ * vector, using @p basis_size_1 rows and @p basis_size_2
+ * columns if interpreted as a matrix.
+ * @param add_into_result Define whether the result should be added into the
+ * array @p values_out (if true) or overwrite the
+ * previous content. The result is undefined in case
+ * values_in and values_out point to the same array and
+ * @p add_into_result is true, in which case an
+ * exception is thrown.
+ * @param values_in The array of the input of size basis_size_2^dim. It
+ * may alias with values_out. Note that the previous
+ * content of @p values_in is overwritten within the
+ * function.
+ * @param values_out The array of size basis_size_1^dim where the results
+ * of the transformation are stored. It may alias with
+ * the @p values_in array.
+ * @param basis_size_1_variable In case the template argument basis_size_1 is
+ * zero, the size of the first basis can alternatively be
+ * passed in as a run time argument. The template
+ * argument takes precedence in case it is nonzero
+ * for efficiency reasons.
+ * @param basis_size_2_variable In case the template argument basis_size_1 is
+ * zero, the size of the second basis can alternatively be
+ * passed in as a run time argument.
+ */
+#ifndef DEBUG
+ DEAL_II_ALWAYS_INLINE
+#endif
+ static void do_backward (const AlignedVector<Number2> &transformation_matrix,
+ const bool add_into_result,
+ Number *values_in,
+ Number *values_out,
+ const unsigned int basis_size_1_variable = numbers::invalid_unsigned_int,
+ const unsigned int basis_size_2_variable = numbers::invalid_unsigned_int)
+ {
+ Assert(basis_size_1 != 0 ||
+ basis_size_1_variable <= basis_size_2_variable,
+ ExcMessage("The second dimension must not be smaller than the first"));
+ Assert(add_into_result == false || values_in != values_out,
+ ExcMessage("Input and output cannot alias with each other when "
+ "adding the result of the basis change to existing data"));
+
+ constexpr int next_dim = (dim > 2 || ((basis_size_1 == 0 || basis_size_2>basis_size_1)
+ && dim>1)) ? dim-1 : dim;
+ EvaluatorTensorProduct<variant, dim, basis_size_1, (basis_size_1==0 ? 0 : basis_size_2),
+ Number,Number2> eval_val (transformation_matrix,
+ AlignedVector<Number2>(),
+ AlignedVector<Number2>(),
+ basis_size_1_variable,
+ basis_size_2_variable);
+ const unsigned int np_1 = basis_size_1 > 0 ? basis_size_1 : basis_size_1_variable;
+ const unsigned int np_2 = basis_size_1 > 0 ? basis_size_2 : basis_size_2_variable;
+ Assert(np_1 > 0 && np_1 != numbers::invalid_unsigned_int,
+ ExcMessage("Cannot transform with 0-point basis"));
+ Assert(np_2 > 0 && np_2 != numbers::invalid_unsigned_int,
+ ExcMessage("Cannot transform with 0-point basis"));
+
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ if (basis_size_1 > 0 && basis_size_2 == basis_size_1 && dim == 2)
+ {
+ eval_val.template values<1,false,false>(values_in, values_in);
+ if (add_into_result)
+ eval_val.template values<0,false,true>(values_in, values_out);
+ else
+ eval_val.template values<0,false,false>(values_in, values_out);
+ }
+ else
+ {
+ if (dim==1 && add_into_result)
+ eval_val.template values<0,false,true>(values_in, values_out);
+ else if (dim==1)
+ eval_val.template values<0,false,false>(values_in, values_out);
+ else
+ eval_val.template values<dim-1,false,false>(values_in, values_in);
+ }
+ if (next_dim < dim)
+ for (unsigned int q=0; q<np_1; ++q)
+ FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,1,Number,Number2>
+ ::do_backward(transformation_matrix,
+ add_into_result,
+ values_in + q*Utilities::fixed_power<next_dim>(np_2),
+ values_out + q*Utilities::fixed_power<next_dim>(np_1),
+ basis_size_1_variable, basis_size_2_variable);
+
+ values_in += Utilities::fixed_power<dim>(np_2);
+ values_out += Utilities::fixed_power<dim>(np_1);
+ }
+ }
+
+ /**
+ * This operation applies a mass-matrix-like operation, consisting of a
+ * do_forward() operation, multiplication by the coefficients in the
+ * quadrature points, and the do_backward() operation.
+ *
+ * @param transformation_matrix The coefficient matrix handed in as a
+ * vector, using @p basis_size_1 rows and @p basis_size_2
+ * columns if interpreted as a matrix.
+ * @param coefficients The array of coefficients by which the result is
+ * multiplied. Its length must be either
+ * basis_size_2^dim or n_components*basis_size_2^dim
+ * @param values_in The array of the input of size basis_size_2^dim. It
+ * may alias with values_out
+ * @param scratch_data Array to hold temporary data during the operation.
+ * Must be of length basis_size_2^dim
+ * @param values_out The array of size basis_size_1^dim where the results
+ * of the transformation are stored. It may alias with
+ * the values_in array.
+ */
+ static void do_mass (const AlignedVector<Number2> &transformation_matrix,
+ const AlignedVector<Number> &coefficients,
+ const Number *values_in,
+ Number *scratch_data,
+ Number *values_out)
+ {
+ constexpr int next_dim = dim > 1 ? dim-1 : dim;
+ Number *my_scratch = basis_size_1 != basis_size_2 ? scratch_data : values_out;
+ for (unsigned int q=basis_size_1; q!=0; --q)
+ FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,n_components,Number,Number2>
+ ::do_forward(transformation_matrix,
+ values_in + (q-1)*Utilities::fixed_int_power<basis_size_1,dim-1>::value,
+ my_scratch + (q-1)*Utilities::fixed_int_power<basis_size_2,dim-1>::value);
+ EvaluatorTensorProduct<variant, dim, basis_size_1, basis_size_2,
+ Number,Number2> eval_val (transformation_matrix);
+ const unsigned int n_inner_blocks = (dim > 1 && basis_size_2 < 10) ? basis_size_2 : 1;
+ const unsigned int n_blocks = Utilities::fixed_int_power<basis_size_2,dim-1>::value;
+ for (unsigned int ii=0; ii<n_blocks; ii+=n_inner_blocks)
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
+ eval_val.template values_one_line<dim-1,true,false> (my_scratch+i, my_scratch+i);
+ for (unsigned int q=0; q<basis_size_2; ++q)
+ for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
+ my_scratch[i+q*n_blocks] *= coefficients[i+q*n_blocks];
+ for (unsigned int i=ii; i<ii+n_inner_blocks; ++i)
+ eval_val.template values_one_line<dim-1,false,false>(my_scratch+i, my_scratch+i);
+ }
+ for (unsigned int q=0; q<basis_size_1; ++q)
+ FEEvaluationImplBasisChange<variant,next_dim,basis_size_1,basis_size_2,n_components,Number,Number2>
+ ::do_backward(transformation_matrix, false,
+ my_scratch + q*Utilities::fixed_int_power<basis_size_2,dim-1>::value,
+ values_out + q*Utilities::fixed_int_power<basis_size_1,dim-1>::value);
+ }
+ };
+
+
+
+ /**
+ * This struct performs the evaluation of function values, gradients and
+ * Hessians for tensor-product finite elements. This a specialization for
+ * elements where the nodal points coincide with the quadrature points like
+ * FE_Q shape functions on Gauss-Lobatto elements integrated with
+ * Gauss-Lobatto quadrature. The assumption of this class is that the shape
+ * 'values' operation is identity, which allows us to write shorter code.
+ *
+ * In literature, this form of evaluation is often called spectral
+ * evaluation, spectral collocation or simply collocation, meaning the same
+ * location for shape functions and evaluation space (quadrature points).
+ *
+ * @author Katharina Kormann, 2012
+ */
template <int dim, int fe_degree, int n_components, typename Number>
- struct FEEvaluationImplTransformToCollocation
+ struct FEEvaluationImplCollocation
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_values,
- const bool evaluate_gradients,
- const bool evaluate_hessians);
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
+ const bool evaluate_values,
+ const bool evaluate_gradients,
+ const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool integrate_values,
- const bool integrate_gradients);
+ void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const bool add_into_values_array);
};
+
+
template <int dim, int fe_degree, int n_components, typename Number>
inline
void
- FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *,
- const bool,
- const bool evaluate_gradients,
- const bool evaluate_hessians)
+ FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *,
+ const bool evaluate_values,
+ const bool evaluate_gradients,
+ const bool evaluate_hessians)
{
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval_val (shape_info.shape_values_eo,
- AlignedVector<VectorizedArray<Number> >(),
- AlignedVector<VectorizedArray<Number> >(),
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
- Eval eval(AlignedVector<VectorizedArray<Number> >(),
- shape_info.shape_gradients_collocation_eo,
- shape_info.shape_hessians_collocation_eo,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler warnings; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:d1;
- const unsigned int d3 = d1+d2;
- const unsigned int d4 = dim>2?4:d3;
- const unsigned int d5 = dim>2?5:d4;
+ AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
+ (fe_degree+2)/2*(fe_degree+1));
+
+ EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree+1, fe_degree+1, Number>
+ eval(AlignedVector<Number>(),
+ shape_info.shape_gradients_collocation_eo,
+ shape_info.shape_hessians_collocation_eo);
+ constexpr unsigned int n_q_points = Utilities::fixed_int_power<fe_degree+1,dim>::value;
for (unsigned int c=0; c<n_components; c++)
{
- // transform to the basis functions of the collocation space. use
- // gradients_quad[c][0] as a temporary array (it gets overwritten by
- // the gradient contributions later)
- if (dim == 1)
- eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]);
- else if (dim == 2)
- {
- eval_val.template values<0,true,false>(values_dofs[c], gradients_quad[c][0]);
- eval_val.template values<1,true,false>(gradients_quad[c][0], values_quad[c]);
- }
- else if (dim == 3)
- {
- eval_val.template values<0,true,false>(values_dofs[c], values_quad[c]);
- eval_val.template values<1,true,false>(values_quad[c], gradients_quad[c][0]);
- eval_val.template values<2,true,false>(gradients_quad[c][0], values_quad[c]);
- }
-
- // apply derivatives in the collocation space
+ if (evaluate_values == true)
+ for (unsigned int i=0; i<n_q_points; ++i)
+ values_quad[i] = values_dofs[i];
if (evaluate_gradients == true || evaluate_hessians == true)
{
- eval.template gradients<0,true,false>(values_quad[c], gradients_quad[c][0]);
+ eval.template gradients<0,true,false>(values_dofs, gradients_quad);
if (dim > 1)
- eval.template gradients<1,true,false>(values_quad[c], gradients_quad[c][d1]);
+ eval.template gradients<1,true,false>(values_dofs, gradients_quad+n_q_points);
if (dim > 2)
- eval.template gradients<2,true,false>(values_quad[c], gradients_quad[c][d2]);
+ eval.template gradients<2,true,false>(values_dofs, gradients_quad+2*n_q_points);
}
if (evaluate_hessians == true)
{
- eval.template hessians<0,true,false> (values_quad[c], hessians_quad[c][0]);
+ eval.template hessians<0,true,false> (values_dofs, hessians_quad);
if (dim > 1)
{
- // re-use grad_x already in gradients
- eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]);
- eval.template hessians<1,true,false> (values_quad[c], hessians_quad[c][d1]);
+ eval.template gradients<1,true,false> (gradients_quad, hessians_quad+dim*n_q_points);
+ eval.template hessians<1,true,false> (values_dofs, hessians_quad+n_q_points);
}
if (dim > 2)
{
- // re-use grad_x and grad_y already in gradients
- eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]);
- eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]);
- eval.template hessians<2,true,false> (values_quad[c], hessians_quad[c][d2]);
+ eval.template gradients<2,true,false> (gradients_quad, hessians_quad+4*n_q_points);
+ eval.template gradients<2,true,false> (gradients_quad+n_q_points, hessians_quad+5*n_q_points);
+ eval.template hessians<2,true,false> (values_dofs, hessians_quad+2*n_q_points);
}
+ hessians_quad += (dim*(dim+1))/2*n_q_points;
}
+ gradients_quad += dim*n_q_points;
+ values_quad += n_q_points;
+ values_dofs += n_q_points;
}
}
+
+
template <int dim, int fe_degree, int n_components, typename Number>
inline
void
- FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *,
- const bool integrate_values,
- const bool integrate_gradients)
+ FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const bool add_into_values_array)
{
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval_val (shape_info.shape_values_eo,
- AlignedVector<VectorizedArray<Number> >(),
- AlignedVector<VectorizedArray<Number> >(),
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
- Eval eval(AlignedVector<VectorizedArray<Number> >(),
- shape_info.shape_gradients_collocation_eo,
- shape_info.shape_hessians_collocation_eo,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler warnings; they are only used in sensible context but
- // compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
+ AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
+ (fe_degree+2)/2*(fe_degree+1));
+
+ EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree+1, fe_degree+1, Number>
+ eval(AlignedVector<Number>(),
+ shape_info.shape_gradients_collocation_eo,
+ shape_info.shape_hessians_collocation_eo);
+ constexpr unsigned int n_q_points = Utilities::fixed_int_power<fe_degree+1,dim>::value;
for (unsigned int c=0; c<n_components; c++)
{
- // apply derivatives in collocation space
+ if (integrate_values == true && add_into_values_array == false)
+ for (unsigned int i=0; i<n_q_points; ++i)
+ values_dofs[i] = values_quad[i];
+ else if (integrate_values == true)
+ for (unsigned int i=0; i<n_q_points; ++i)
+ values_dofs[i] += values_quad[i];
if (integrate_gradients == true)
{
- if (integrate_values)
- eval.template gradients<0,false,true>(gradients_quad[c][0], values_quad[c]);
+ if (integrate_values == true || add_into_values_array == true)
+ eval.template gradients<0,false,true>(gradients_quad, values_dofs);
else
- eval.template gradients<0,false,false>(gradients_quad[c][0], values_quad[c]);
+ eval.template gradients<0,false,false>(gradients_quad, values_dofs);
if (dim > 1)
- eval.template gradients<1,false,true>(gradients_quad[c][d1], values_quad[c]);
+ eval.template gradients<1,false,true>(gradients_quad+n_q_points, values_dofs);
if (dim > 2)
- eval.template gradients<2,false,true>(gradients_quad[c][d2], values_quad[c]);
- }
-
- // transform back to the original space
- if (dim == 1)
- eval_val.template values<0,false,false>(values_quad[c], values_dofs[c]);
- else if (dim == 2)
- {
- eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
- eval_val.template values<1,false,false>(gradients_quad[c][0], values_dofs[c]);
- }
- else if (dim == 3)
- {
- eval_val.template values<0,false,false>(values_quad[c], gradients_quad[c][0]);
- eval_val.template values<1,false,false>(gradients_quad[c][0], values_quad[c]);
- eval_val.template values<2,false,false>(values_quad[c], values_dofs[c]);
+ eval.template gradients<2,false,true>(gradients_quad+2*n_q_points, values_dofs);
}
+ gradients_quad += dim*n_q_points;
+ values_quad += n_q_points;
+ values_dofs += n_q_points;
}
}
+
/**
* This struct performs the evaluation of function values, gradients and
* Hessians for tensor-product finite elements. This a specialization for
- * elements where the nodal points coincide with the quadrature points like
- * FE_Q shape functions on Gauss-Lobatto elements integrated with
- * Gauss-Lobatto quadrature. The assumption of this class is that the shape
- * 'values' operation is identity, which allows us to write shorter code.
- *
- * In literature, this form of evaluation is often called spectral
- * evaluation, spectral collocation or simply collocation, meaning the same
- * location for shape functions and evaluation space (quadrature points).
+ * symmetric basis functions about the mid point 0.5 of the unit interval
+ * with the same number of quadrature points as degrees of freedom. In that
+ * case, we can first transform the basis to one that has the nodal points
+ * in the quadrature points (i.e., the collocation space) and then perform
+ * the evaluation of the first and second derivatives in this transformed
+ * space, using the identity operation for the shape values.
*
- * @author Katharina Kormann, 2012
- */
- template <int dim, int fe_degree, int n_components, typename Number>
- struct FEEvaluationImplCollocation
+ * @author Katharina Kormann, Martin Kronbichler, 2017
+ */
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+ struct FEEvaluationImplTransformToCollocation
{
static
- void evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
- const bool evaluate_values,
- const bool evaluate_gradients,
- const bool evaluate_hessians);
+ void evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
+ const bool evaluate_values,
+ const bool evaluate_gradients,
+ const bool evaluate_hessians);
static
- void integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
- const bool integrate_values,
- const bool integrate_gradients);
+ void integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const bool add_into_values_array);
};
- template <int dim, int fe_degree, int n_components, typename Number>
+
+
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
inline
void
- FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
- ::evaluate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *,
- const bool evaluate_values,
- const bool evaluate_gradients,
- const bool evaluate_hessians)
+ FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
+ ::evaluate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ const Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *,
+ const bool ,
+ const bool evaluate_gradients,
+ const bool evaluate_hessians)
{
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval(AlignedVector<VectorizedArray<Number> >(),
- shape_info.shape_gradients_eo,
- shape_info.shape_hessians_eo,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler warnings; they are only used in sensible context
- // but compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:d1;
- const unsigned int d3 = d1+d2;
- const unsigned int d4 = dim>2?4:d3;
- const unsigned int d5 = dim>2?5:d4;
+ Assert(n_q_points_1d > fe_degree,
+ ExcMessage("You lose information when going to a collocation space "
+ "of lower degree, so the evaluation results would be "
+ "wrong. Thus, this class does not permit the desired "
+ "operation."));
+ constexpr unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
for (unsigned int c=0; c<n_components; c++)
{
- if (evaluate_values == true)
- for (unsigned int i=0; i<Eval::dofs_per_cell; ++i)
- values_quad[c][i] = values_dofs[c][i];
+ FEEvaluationImplBasisChange<evaluate_evenodd, dim,
+ (fe_degree>=n_q_points_1d?n_q_points_1d:fe_degree+1),
+ n_q_points_1d,1,Number,Number>
+ ::do_forward(shape_info.shape_values_eo,
+ values_dofs, values_quad);
+
+ // apply derivatives in the collocation space
if (evaluate_gradients == true || evaluate_hessians == true)
- {
- eval.template gradients<0,true,false>(values_dofs[c], gradients_quad[c][0]);
- if (dim > 1)
- eval.template gradients<1,true,false>(values_dofs[c], gradients_quad[c][d1]);
- if (dim > 2)
- eval.template gradients<2,true,false>(values_dofs[c], gradients_quad[c][d2]);
- }
- if (evaluate_hessians == true)
- {
- eval.template hessians<0,true,false> (values_dofs[c], hessians_quad[c][0]);
- if (dim > 1)
- {
- // re-use grad_x already in gradients
- eval.template gradients<1,true,false> (gradients_quad[c][0], hessians_quad[c][d3]);
- eval.template hessians<1,true,false> (values_dofs[c], hessians_quad[c][d1]);
- }
- if (dim > 2)
- {
- // re-use grad_x already in gradients
- eval.template gradients<2,true,false> (gradients_quad[c][0], hessians_quad[c][d4]);
- eval.template gradients<2,true,false> (gradients_quad[c][d1], hessians_quad[c][d5]);
- eval.template hessians<2,true,false> (values_dofs[c], hessians_quad[c][d2]);
- }
- }
+ FEEvaluationImplCollocation<dim,n_q_points_1d-1,1,Number>::
+ evaluate(shape_info, values_quad, nullptr, gradients_quad, hessians_quad,
+ nullptr, false, evaluate_gradients, evaluate_hessians);
+
+ values_dofs += shape_info.dofs_per_component_on_cell;
+ values_quad += n_q_points;
+ gradients_quad += dim*n_q_points;
+ hessians_quad += (dim*(dim+1))/2*n_q_points;
}
}
- template <int dim, int fe_degree, int n_components, typename Number>
+
+
+ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
inline
void
- FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
- ::integrate (const MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number>> &shape_info,
- VectorizedArray<Number> *values_dofs[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *,
- const bool integrate_values,
- const bool integrate_gradients)
+ FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
+ ::integrate (const MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *,
+ const bool integrate_values,
+ const bool integrate_gradients,
+ const bool add_into_values_array)
{
- typedef EvaluatorTensorProduct<evaluate_evenodd, dim, fe_degree, fe_degree+1,
- VectorizedArray<Number> > Eval;
- Eval eval(AlignedVector<VectorizedArray<Number> >(),
- shape_info.shape_gradients_eo,
- shape_info.shape_hessians_eo,
- shape_info.fe_degree,
- shape_info.n_q_points_1d);
-
- // These avoid compiler warnings; they are only used in sensible context
- // but compilers typically cannot detect when we access something like
- // gradients_quad[2] only for dim==3.
- const unsigned int d1 = dim>1?1:0;
- const unsigned int d2 = dim>2?2:0;
+ Assert(n_q_points_1d > fe_degree,
+ ExcMessage("You lose information when going to a collocation space "
+ "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);
+ constexpr unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
for (unsigned int c=0; c<n_components; c++)
{
- if (integrate_values == true)
- for (unsigned int i=0; i<Eval::dofs_per_cell; ++i)
- values_dofs[c][i] = values_quad[c][i];
+
+ // apply derivatives in collocation space
if (integrate_gradients == true)
- {
- if (integrate_values == true)
- eval.template gradients<0,false,true>(gradients_quad[c][0], values_dofs[c]);
- else
- eval.template gradients<0,false,false>(gradients_quad[c][0], values_dofs[c]);
- if (dim > 1)
- eval.template gradients<1,false,true>(gradients_quad[c][d1], values_dofs[c]);
- if (dim > 2)
- eval.template gradients<2,false,true>(gradients_quad[c][d2], values_dofs[c]);
- }
+ FEEvaluationImplCollocation<dim,n_q_points_1d-1,1,Number>::
+ integrate(shape_info, values_quad, nullptr, gradients_quad, nullptr, false,
+ integrate_gradients,/*add_into_values_array=*/integrate_values);
+
+ // transform back to the original space
+ FEEvaluationImplBasisChange<evaluate_evenodd, dim,
+ (fe_degree>=n_q_points_1d?n_q_points_1d:fe_degree+1),
+ n_q_points_1d,1,Number,Number>
+ ::do_backward(shape_info.shape_values_eo,
+ add_into_values_array,
+ values_quad,
+ values_dofs);
+ gradients_quad += dim*n_q_points;
+ values_quad += n_q_points;
+ values_dofs += shape_info.dofs_per_component_on_cell;
}
}
+
+
+ template <bool symmetric_evaluate, int dim, int fe_degree, int n_q_points_1d, int n_components, typename Number>
+ struct FEFaceEvaluationImpl
+ {
+ static
+ void evaluate_in_face (const MatrixFreeFunctions::ShapeInfo<Number> &data,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
+ const bool evaluate_val,
+ const bool evaluate_grad,
+ const unsigned int subface_index)
+ {
+ const AlignedVector<Number> &val1
+ = symmetric_evaluate ? data.shape_values_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_values : data.values_within_subface[subface_index%2]);
+ const AlignedVector<Number> &val2
+ = symmetric_evaluate ? data.shape_values_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_values : data.values_within_subface[subface_index/2]);
+
+ const AlignedVector<Number> &grad1
+ = symmetric_evaluate ? data.shape_gradients_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_gradients : data.gradients_within_subface[subface_index%2]);
+ const AlignedVector<Number> &grad2
+ = symmetric_evaluate ? data.shape_gradients_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_gradients : data.gradients_within_subface[subface_index/2]);
+
+ typedef internal::EvaluatorTensorProduct
+ <symmetric_evaluate ? internal::evaluate_evenodd :internal::evaluate_general,
+ dim-1,fe_degree+1,n_q_points_1d,Number> Eval;
+ typedef internal::EvaluatorTensorProduct
+ <internal::evaluate_general,dim-1,fe_degree+1,n_q_points_1d,
+ Number> EvalGeneric;
+ Eval eval1(val1,grad1,AlignedVector<Number>(),
+ data.fe_degree+1, data.n_q_points_1d);
+ Eval eval2(val2,grad2,AlignedVector<Number>(),
+ data.fe_degree+1, data.n_q_points_1d);
+
+ const unsigned int size_deg = fe_degree > -1 ?
+ Utilities::fixed_int_power<fe_degree+1,dim-1>::value :
+ (dim > 1 ? Utilities::fixed_power<dim-1>(data.fe_degree+1) : 1);
+
+ const unsigned int n_q_points = fe_degree > -1 ?
+ Utilities::fixed_int_power<n_q_points_1d,dim-1>::value : data.n_q_points_face;
+
+ if (evaluate_grad == false)
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ switch (dim)
+ {
+ case 3:
+ eval1.template values<0,true,false>(values_dofs, values_quad);
+ eval2.template values<1,true,false>(values_quad, values_quad);
+ break;
+ case 2:
+ eval1.template values<0,true,false>(values_dofs, values_quad);
+ break;
+ case 1:
+ values_quad[c] = values_dofs[2*c];
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ values_dofs += 2*size_deg;
+ values_quad += n_q_points;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ switch (dim)
+ {
+ case 3:
+ if (symmetric_evaluate && n_q_points_1d > fe_degree)
+ {
+ eval1.template values<0,true,false>(values_dofs, values_quad);
+ eval1.template values<1,true,false>(values_quad, values_quad);
+ internal::EvaluatorTensorProduct
+ <internal::evaluate_evenodd,dim-1,n_q_points_1d,n_q_points_1d,Number> eval_grad
+ (AlignedVector<Number>(),
+ data.shape_gradients_collocation_eo,
+ AlignedVector<Number>());
+ eval_grad.template gradients<0,true,false>(values_quad, gradients_quad);
+ eval_grad.template gradients<1,true,false>(values_quad,
+ gradients_quad+n_q_points);
+ }
+ else
+ {
+ eval1.template gradients<0,true,false>(values_dofs, scratch_data);
+ eval2.template values<1,true,false>(scratch_data, gradients_quad);
+
+ eval1.template values<0,true,false>(values_dofs, scratch_data);
+ eval2.template gradients<1,true,false>(scratch_data, gradients_quad+n_q_points);
+
+ if (evaluate_val == true)
+ eval2.template values<1,true,false>(scratch_data, values_quad);
+ }
+ eval1.template values<0,true,false>(values_dofs+size_deg, scratch_data);
+ eval2.template values<1,true,false>(scratch_data,
+ gradients_quad+(dim-1)*n_q_points);
+
+ break;
+ case 2:
+ eval1.template values<0,true,false>(values_dofs+size_deg,
+ gradients_quad+(dim-1)*n_q_points);
+ eval1.template gradients<0,true,false>(values_dofs, gradients_quad);
+ if (evaluate_val == true)
+ eval1.template values<0,true,false>(values_dofs, values_quad);
+ break;
+ case 1:
+ values_quad[0] = values_dofs[0];
+ gradients_quad[0] = values_dofs[1];
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ values_dofs += 2*size_deg;
+ values_quad += n_q_points;
+ gradients_quad += dim*n_q_points;
+ }
+ }
+
+ static
+ void integrate_in_face (const MatrixFreeFunctions::ShapeInfo<Number> &data,
+ Number *values_dofs,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
+ const bool integrate_val,
+ const bool integrate_grad,
+ const unsigned int subface_index)
+ {
+ const AlignedVector<Number> &val1
+ = symmetric_evaluate ? data.shape_values_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_values : data.values_within_subface[subface_index%2]);
+ const AlignedVector<Number> &val2
+ = symmetric_evaluate ? data.shape_values_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_values : data.values_within_subface[subface_index/2]);
+
+ const AlignedVector<Number> &grad1
+ = symmetric_evaluate ? data.shape_gradients_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_gradients : data.gradients_within_subface[subface_index%2]);
+ const AlignedVector<Number> &grad2
+ = symmetric_evaluate ? data.shape_gradients_eo :
+ (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
+ data.shape_gradients : data.gradients_within_subface[subface_index/2]);
+
+ typedef internal::EvaluatorTensorProduct
+ <symmetric_evaluate ? internal::evaluate_evenodd :internal::evaluate_general,
+ dim-1,fe_degree+1,n_q_points_1d,Number> Eval;
+ typedef internal::EvaluatorTensorProduct
+ <internal::evaluate_general,dim-1,fe_degree+1,n_q_points_1d,
+ Number> EvalGeneric;
+ 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);
+
+ const unsigned int size_deg = fe_degree > -1 ?
+ Utilities::fixed_int_power<fe_degree+1,dim-1>::value :
+ (dim > 1 ? Utilities::fixed_power<dim-1>(data.fe_degree+1) : 1);
+
+ const unsigned int n_q_points = fe_degree > -1 ?
+ Utilities::fixed_int_power<n_q_points_1d,dim-1>::value : data.n_q_points_face;
+
+ if (integrate_grad == false)
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ switch (dim)
+ {
+ case 3:
+ eval2.template values<1,false,false>(values_quad, values_quad);
+ eval1.template values<0,false,false>(values_quad, values_dofs);
+ break;
+ case 2:
+ eval1.template values<0,false,false>(values_quad, values_dofs);
+ break;
+ case 1:
+ values_dofs[2*c] = values_quad[c][0];
+ break;
+ default:
+ Assert(false, ExcNotImplemented());
+ }
+ values_dofs += 2*size_deg;
+ values_quad += n_q_points;
+ }
+ else
+ for (unsigned int c=0; c<n_components; ++c)
+ {
+ switch (dim)
+ {
+ case 3:
+ eval2.template values<1,false,false> (gradients_quad+2*n_q_points,
+ gradients_quad+2*n_q_points);
+ eval1.template values<0,false,false> (gradients_quad+2*n_q_points,
+ values_dofs+size_deg);
+ if (symmetric_evaluate && n_q_points_1d > fe_degree)
+ {
+ internal::EvaluatorTensorProduct <internal::evaluate_evenodd,
+ dim-1,n_q_points_1d,n_q_points_1d,Number> eval_grad
+ (AlignedVector<Number>(),
+ data.shape_gradients_collocation_eo,
+ AlignedVector<Number>());
+ if (integrate_val)
+ eval_grad.template gradients<1,false,true>(gradients_quad+n_q_points,
+ values_quad);
+ else
+ eval_grad.template gradients<1,false,false>(gradients_quad+n_q_points,
+ values_quad);
+ eval_grad.template gradients<0,false,true>(gradients_quad,
+ values_quad);
+ eval1.template values<1,false,false>(values_quad, values_quad);
+ eval1.template values<0,false,false>(values_quad, values_dofs);
+ }
+ else
+ {
+ if (integrate_val)
+ {
+ eval2.template values<1,false,false> (values_quad, scratch_data);
+ eval2.template gradients<1,false,true> (gradients_quad+n_q_points,
+ scratch_data);
+ }
+ else
+ eval2.template gradients<1,false,false> (gradients_quad+n_q_points,
+ scratch_data);
+
+ eval1.template values<0,false,false> (scratch_data, values_dofs);
+ eval2.template values<1,false,false> (gradients_quad, scratch_data);
+ eval1.template gradients<0,false,true> (scratch_data, values_dofs);
+ }
+ break;
+ case 2:
+ eval1.template values<0,false,false>(gradients_quad+n_q_points,
+ values_dofs+size_deg);
+ eval1.template gradients<0,false,false>(gradients_quad, values_dofs);
+ if (integrate_val == true)
+ eval1.template values<0,false,true>(values_quad, values_dofs);
+ break;
+ case 1:
+ values_dofs[0] = values_quad[0];
+ values_dofs[1] = gradients_quad[0];
+ break;
+ default:
+ AssertThrow(false, ExcNotImplemented());
+ }
+ values_dofs += 2*size_deg;
+ values_quad += n_q_points;
+ gradients_quad += dim*n_q_points;
+ }
+ }
+ };
+
+
+
+ template <int dim, int fe_degree, int n_components, typename Number>
+ struct FEFaceNormalEvaluationImpl
+ {
+ template <bool do_evaluate, bool add_into_output>
+ static void interpolate(const MatrixFreeFunctions::ShapeInfo<Number> &data,
+ const Number *input,
+ Number *output,
+ const bool do_gradients,
+ const unsigned int face_no)
+ {
+ internal::EvaluatorTensorProduct<internal::evaluate_general,dim,
+ fe_degree+1,0,Number>
+ evalf(data.shape_data_on_face[face_no%2],
+ AlignedVector<Number>(),
+ AlignedVector<Number>(),
+ data.fe_degree+1, 0);
+
+ const unsigned int in_stride = do_evaluate ? data.dofs_per_component_on_cell : 2*data.dofs_per_component_on_face;
+ const unsigned int out_stride = do_evaluate ? 2*data.dofs_per_component_on_face : data.dofs_per_component_on_cell;
+ const unsigned int face_direction = face_no / 2;
+ for (unsigned int c=0; c<n_components; c++)
+ {
+ if (do_gradients)
+ {
+ if (face_direction == 0)
+ evalf.template apply_face<0,do_evaluate,add_into_output,1>(input, output);
+ else if (face_direction == 1)
+ evalf.template apply_face<1,do_evaluate,add_into_output,1>(input, output);
+ else
+ evalf.template apply_face<2,do_evaluate,add_into_output,1>(input, output);
+ }
+ else
+ {
+ if (face_direction == 0)
+ evalf.template apply_face<0,do_evaluate,add_into_output,0>(input, output);
+ else if (face_direction == 1)
+ evalf.template apply_face<1,do_evaluate,add_into_output,0>(input, output);
+ else
+ evalf.template apply_face<2,do_evaluate,add_into_output,0>(input, output);
+ }
+ input += in_stride;
+ output += out_stride;
+ }
+ }
+ };
} // end of namespace internal
template <int dim, int n_components, typename Number>
struct Default
{
- static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
evaluate_values, evaluate_gradients, evaluate_hessians);
}
- static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
dim, -1, 0, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
};
struct Factory<dim, n_components, Number, 0, degree, n_q_points_1d>
{
static inline void evaluate (
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
}
static inline void integrate (
- const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
template<int degree, int n_q_points_1d, int dim, int n_components, typename Number>
struct Factory<dim, n_components, Number, 1, degree, n_q_points_1d, typename std::enable_if<(n_q_points_1d<degree+3)>::type>
{
- static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
if (runtime_n_q_points_1d == n_q_points_1d)
{
- if (n_q_points_1d == degree+1)
- {
- if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
- internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
- ::evaluate(shape_info, values_dofs_actual, values_quad,
- gradients_quad, hessians_quad, scratch_data,
- evaluate_values, evaluate_gradients, evaluate_hessians);
- else
- internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
- ::evaluate(shape_info, values_dofs_actual, values_quad,
- gradients_quad, hessians_quad, scratch_data,
- evaluate_values, evaluate_gradients, evaluate_hessians);
- }
+ if (n_q_points_1d == degree+1 &&
+ shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+ internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
+ ::evaluate(shape_info, values_dofs_actual, values_quad,
+ gradients_quad, hessians_quad, scratch_data,
+ evaluate_values, evaluate_gradients, evaluate_hessians);
+ else if (degree < n_q_points_1d)
+ internal::FEEvaluationImplTransformToCollocation<dim, degree, n_q_points_1d, n_components, Number>
+ ::evaluate(shape_info, values_dofs_actual, values_quad,
+ gradients_quad, hessians_quad, scratch_data,
+ evaluate_values, evaluate_gradients, evaluate_hessians);
else
internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
::evaluate(shape_info, values_dofs_actual, values_quad,
evaluate_values, evaluate_gradients, evaluate_hessians);
}
- static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
if (runtime_n_q_points_1d == n_q_points_1d)
{
- if (n_q_points_1d == degree+1)
- {
- if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
- internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
- ::integrate(shape_info, values_dofs_actual, values_quad,
- gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
- else
- internal::FEEvaluationImplTransformToCollocation<dim, degree, n_components, Number>
- ::integrate(shape_info, values_dofs_actual, values_quad,
- gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
- }
+ if (n_q_points_1d == degree+1 &&
+ shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_collocation)
+ internal::FEEvaluationImplCollocation<dim, degree, n_components, Number>
+ ::integrate(shape_info, values_dofs_actual, values_quad,
+ gradients_quad, scratch_data,
+ integrate_values, integrate_gradients, false);
+ else if (degree < n_q_points_1d)
+ internal::FEEvaluationImplTransformToCollocation<dim, degree, n_q_points_1d, n_components, Number>
+ ::integrate(shape_info, values_dofs_actual, values_quad,
+ gradients_quad, scratch_data,
+ integrate_values, integrate_gradients, false);
else
internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_symmetric, dim, degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad, gradients_quad,
- scratch_data, integrate_values, integrate_gradients);
+ scratch_data, integrate_values, integrate_gradients, false);
}
else
Factory<dim, n_components, Number, 1, degree, n_q_points_1d+1>
* for the 'evaluate' function.
*/
template<int dim, int n_components, typename Number>
- void symmetric_selector_evaluate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ void symmetric_selector_evaluate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
* for the 'integrate' function.
*/
template<int dim, int n_components, typename Number>
- void symmetric_selector_integrate (const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ void symmetric_selector_integrate (const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
* internal::FEEvaluationImplTransformToCollocation::evaluate() with appropriate
* template parameters.
*/
- static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians);
* internal::FEEvaluationImplTransformToCollocation::integrate() with appropriate
* template parameters.
*/
- static void integrate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ static void integrate(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients);
};
* internal::FEEvaluationImplTransformToCollocation::evaluate() with appropriate
* template parameters.
*/
- static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+ static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians);
* internal::FEEvaluationImplTransformToCollocation::integrate() with appropriate
* template parameters.
*/
- static void integrate(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+ static void integrate(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients);
};
inline
void
SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate
-(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
gradients_quad, hessians_quad, scratch_data,
evaluate_values, evaluate_gradients, evaluate_hessians);
}
- else if (fe_degree+1 == n_q_points_1d &&
+ else if (fe_degree < n_q_points_1d &&
shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
{
- internal::FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
+ internal::FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
::evaluate(shape_info, values_dofs_actual, values_quad,
gradients_quad, hessians_quad, scratch_data,
evaluate_values, evaluate_gradients, evaluate_hessians);
inline
void
SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate
-(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
internal::FEEvaluationImplCollocation<dim, fe_degree, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
- else if (fe_degree+1 == n_q_points_1d &&
+ else if (fe_degree < n_q_points_1d &&
shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
{
- internal::FEEvaluationImplTransformToCollocation<dim, fe_degree, n_components, Number>
+ internal::FEEvaluationImplTransformToCollocation<dim, fe_degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric)
{
dim, fe_degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
{
dim, fe_degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
{
dim, fe_degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
{
dim, fe_degree, n_q_points_1d, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else
AssertThrow(false, ExcNotImplemented());
inline
void
SelectEvaluator<dim, -1, dummy, n_components, Number>::evaluate
-(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *hessians_quad[][(dim*(dim+1))/2],
- VectorizedArray<Number> *scratch_data,
+(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *hessians_quad,
+ Number *scratch_data,
const bool evaluate_values,
const bool evaluate_gradients,
const bool evaluate_hessians)
inline
void
SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate
-(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<Number> > &shape_info,
- VectorizedArray<Number> *values_dofs_actual[],
- VectorizedArray<Number> *values_quad[],
- VectorizedArray<Number> *gradients_quad[][dim],
- VectorizedArray<Number> *scratch_data,
+(const internal::MatrixFreeFunctions::ShapeInfo<Number> &shape_info,
+ Number *values_dofs_actual,
+ Number *values_quad,
+ Number *gradients_quad,
+ Number *scratch_data,
const bool integrate_values,
const bool integrate_gradients)
{
dim, -1, 0, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::truncated_tensor)
{
dim, -1, 0, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
}
else if (shape_info.element_type == internal::MatrixFreeFunctions::tensor_general)
internal::FEEvaluationImpl<internal::MatrixFreeFunctions::tensor_general,
dim, -1, 0, n_components, Number>
::integrate(shape_info, values_dofs_actual, values_quad,
gradients_quad, scratch_data,
- integrate_values, integrate_gradients);
+ integrate_values, integrate_gradients, false);
else
symmetric_selector_integrate<dim, n_components, Number>
(shape_info, values_dofs_actual, values_quad,
Assert(this->matrix_info != nullptr ||
this->mapped_geometry->is_initialized(), ExcNotInitialized());
- SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>
- ::evaluate (*this->data, &this->values_dofs[0], this->values_quad,
- this->gradients_quad, this->hessians_quad, this->scratch_data,
+ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, VectorizedArray<Number> >
+ ::evaluate (*this->data, this->values_dofs[0], this->values_quad[0],
+ this->gradients_quad[0][0], this->hessians_quad[0][0], this->scratch_data,
evaluate_values, evaluate_gradients, evaluate_hessians);
#ifdef DEBUG
Assert(this->matrix_info != nullptr ||
this->mapped_geometry->is_initialized(), ExcNotInitialized());
- SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>
- ::integrate (*this->data, &this->values_dofs[0], this->values_quad,
- this->gradients_quad, this->scratch_data,
+ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, VectorizedArray<Number> >
+ ::integrate (*this->data, this->values_dofs[0], this->values_quad[0],
+ this->gradients_quad[0][0], this->scratch_data,
integrate_values, integrate_gradients);
#ifdef DEBUG
Assert(dim == 2 || dim == 3, ExcNotImplemented());
- internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,fe_degree,
+ internal::EvaluatorTensorProduct<internal::evaluate_evenodd,dim,fe_degree+1,
fe_degree+1, VectorizedArray<Number> >
evaluator(inverse_shape, inverse_shape, inverse_shape);
// ---------------------------------------------------------------------
//
-// Copyright (C) 2017 by the deal.II authors
+// Copyright (C) 2017-2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
* vector separately: see the documentation of the EvaluatorTensorProduct
* specialization for more information.
*/
- evaluate_evenodd
+ evaluate_evenodd,
+ /**
+ * Use symmetry in Legendre and similar polynomial spaces where the shape
+ * functions with even number are symmetric about the center of the
+ * quadrature points (think about even polynomial degrees) and the shape
+ * functions with odd number are anti-symmetric about the center of the
+ * quadrature points (think about odd polynomial degrees). This allows to
+ * use a strategy similar to the even-odd technique but without separate
+ * coefficient arrays. See the documentation of the EvaluatorTensorProduct
+ * specialization for more information.
+ */
+ evaluate_symmetric_hierarchical
};
+
+
/**
- * Generic evaluator framework
+ * Generic evaluator framework that valuates the given shape data in general
+ * dimensions using the tensor product form. Depending on the particular
+ * layout in the matrix entries, this corresponds to a usual matrix-matrix
+ * product or a matrix-matrix product including some symmetries.
+ *
+ * @tparam variant Variant of evaluation used for creating template
+ * specializations
+ * @tparam dim Dimension of the function
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number to be valid
*/
- template <EvaluatorVariant variant, int dim, int fe_degree, int n_q_points_1d,
- typename Number>
+ template <EvaluatorVariant variant, int dim, int n_rows, int n_columns,
+ typename Number, typename Number2=Number>
struct EvaluatorTensorProduct
{};
+
+
/**
- * Internal evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions
+ * Internal evaluator for shape function in arbitrary dimension using the
+ * tensor product form of the basis functions.
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
*/
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ struct EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
{
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+ static const unsigned int n_rows_of_product = Utilities::fixed_int_power<n_rows,dim>::value;
+ static const unsigned int n_columns_of_product = Utilities::fixed_int_power<n_columns,dim>::value;
/**
* Empty constructor. Does nothing. Be careful when using 'values' and
/**
* Constructor, taking the data from ShapeInfo
*/
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
:
shape_values (shape_values.begin()),
shape_gradients (shape_gradients.begin()),
shape_hessians (shape_hessians.begin())
{
+ // We can enter this function either for the apply() path that has
+ // n_rows * n_columns entries or for the apply_face() path that only has
+ // n_rows * 3 entries in the array. Since we cannot decide about the use
+ // we must allow for both here.
+ Assert(shape_values.size() == 0 ||
+ shape_values.size() == n_rows*n_columns ||
+ shape_values.size() == 3*n_rows,
+ ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
+ Assert(shape_gradients.size() == 0 ||
+ shape_gradients.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
+ Assert(shape_hessians.size() == 0 ||
+ shape_hessians.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
(void)dummy1;
(void)dummy2;
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
values (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add>(shape_values, in, out);
+ apply<direction,contract_over_rows,add>(shape_values, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
gradients (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add>(shape_gradients, in, out);
+ apply<direction,contract_over_rows,add>(shape_gradients, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
hessians (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add>(shape_hessians, in, out);
+ apply<direction,contract_over_rows,add>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ values_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_values, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ gradients_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ hessians_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_hessians, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
- static void apply (const Number *shape_data,
- const Number in [],
- Number out []);
+ /**
+ * This function applies the tensor product kernel, corresponding to a
+ * multiplication of 1D stripes, along the given @p direction of the tensor
+ * data in the input array. This function allows the @p in and @p out
+ * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
+ * perform the contraction in place where @p in and @p out point to the
+ * same address. For the case n_rows != n_columns, the output is in general
+ * not correct.
+ *
+ * @tparam direction Direction that is evaluated
+ * @tparam contract_over_rows If true, the tensor contraction sums
+ * over the rows in the given @p shape_data
+ * array, otherwise it sums over the columns
+ * @tparam add If true, the result is added to the output vector, else
+ * the computed values overwrite the content in the output
+ * @tparam one_line If true, the kernel is only applied along a single 1D
+ * stripe within a dim-dimensional tensor, not the full
+ * n_rows^dim points as in the @p false case.
+ *
+ * @param shape_data Transformation matrix with @p n_rows rows and
+ * @p n_columns columns, stored in row-major format
+ * @param in Pointer to the start of the input data vector
+ * @param out Pointer to the start of the output data vector
+ */
+ template <int direction, bool contract_over_rows, bool add, bool one_line=false>
+ static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out);
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
+ /**
+ * This function applies the tensor product operation to produce face values
+ * from cell values. As opposed to the apply method, this method assumes
+ * that the directions orthogonal to the face have n_rows degrees of
+ * freedom per direction and not n_columns for those directions lower than
+ * the one currently applied. In other words, apply_face() must be called
+ * before calling any interpolation within the face.
+ *
+ * @tparam face_direction Direction of the normal vector (0=x, 1=y, etc)
+ * @tparam contract_onto_face If true, the input vector is of size n_rows^dim
+ * and interpolation into n_rows^(dim-1) points
+ * is performed. This is a typical scenario in
+ * FEFaceEvaluation::evaluate() calls. If false,
+ * data from n_rows^(dim-1) points is expanded
+ * into the n_rows^dim points of the higher-
+ * dimensional data array. Derivatives in the
+ * case contract_onto_face==false are summed
+ * together
+ * @tparam add If true, the result is added to the output vector, else
+ * the computed values overwrite the content in the output
+ * @tparam max_derivative Sets the number of derivatives that should be
+ * computed. 0 means only values, 1 means values and first
+ * derivatives, 2 second derivates. Note that all the
+ * derivatives access the data in @p shape_values passed to
+ * the constructor of the class
+ *
+ * @param in address of the input data vector
+ * @param out address of the output data vector
+ */
+ template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
+ void apply_face (const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT out) const;
+
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
};
- // evaluates the given shape data in 1d-3d using the tensor product
- // form. does not use a particular layout of entries in the matrices
- // like the functions below and corresponds to a usual matrix-matrix
- // product
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add>
+
+
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add, bool one_line>
inline
void
- EvaluatorTensorProduct<evaluate_general,dim,fe_degree,n_q_points_1d,Number>
- ::apply (const Number *shape_data,
- const Number in [],
- Number out [])
+ EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
+ ::apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out)
{
+ static_assert (one_line == false || direction==dim-1,
+ "Single-line evaluation only works for direction=dim-1.");
+ Assert(shape_data != nullptr,
+ ExcMessage("The given array shape_data must not be the null pointer!"));
+ Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
+ ExcMessage("In-place operation only supported for "
+ "n_rows==n_columns or single-line interpolation"));
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+ constexpr int mm = contract_over_rows ? n_rows : n_columns,
+ nn = contract_over_rows ? n_columns : n_rows;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = Utilities::fixed_int_power<nn,direction>::value;
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = one_line ? 1 : stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
+ Number x[mm];
+ for (int i=0; i<mm; ++i)
+ x[i] = in[stride*i];
for (int col=0; col<nn; ++col)
{
- Number val0;
- if (dof_to_quad == true)
+ Number2 val0;
+ if (contract_over_rows == true)
val0 = shape_data[col];
else
- val0 = shape_data[col*n_q_points_1d];
- Number res0 = val0 * in[0];
- for (int ind=1; ind<mm; ++ind)
+ val0 = shape_data[col*n_columns];
+ Number res0 = val0 * x[0];
+ for (int i=1; i<mm; ++i)
{
- if (dof_to_quad == true)
- val0 = shape_data[ind*n_q_points_1d+col];
+ if (contract_over_rows == true)
+ val0 = shape_data[i*n_columns+col];
else
- val0 = shape_data[col*n_q_points_1d+ind];
- res0 += val0 * in[stride*ind];
+ val0 = shape_data[col*n_columns+i];
+ res0 += val0 * x[i];
}
if (add == false)
out[stride*col] = res0;
out[stride*col] += res0;
}
- // increment: in regular case, just go to the next point in
- // x-direction. If we are at the end of one chunk in x-dir, need
- // to jump over to the next layer in z-direction
- switch (direction)
+ if (one_line == false)
{
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
++in;
++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
}
}
- if (direction == 1)
+ if (one_line == false)
{
- in += nn*(mm-1);
- out += nn*(nn-1);
+ in += stride*(mm-1);
+ out += stride*(nn-1);
}
}
}
- // This method applies the tensor product operation to produce face values
- // out from cell values. As opposed to the apply_tensor_product method, this
- // method assumes that the directions orthogonal to the face have
- // fe_degree+1 degrees of freedom per direction and not n_q_points_1d for
- // those directions lower than the one currently applied
- template <int dim, int fe_degree, typename Number, int face_direction,
- bool dof_to_quad, bool add>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
inline
void
- apply_tensor_product_face (const Number *shape_data,
- const Number in [],
- Number out [])
+ EvaluatorTensorProduct<evaluate_general,dim,n_rows,n_columns,Number,Number2>
+ ::apply_face (const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT out) const
{
- const int n_blocks1 = dim > 1 ? (fe_degree+1) : 1;
- const int n_blocks2 = dim > 2 ? (fe_degree+1) : 1;
+ static_assert(dim > 0 && dim<4, "Only dim=1,2,3 supported");
+ static_assert(max_derivative >= 0 && max_derivative<3,
+ "Only derivative orders 0-2 implemented");
+ Assert(shape_values != nullptr,
+ ExcMessage("The given array shape_values must not be the null pointer."));
- AssertIndexRange (face_direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : 1,
- nn = dof_to_quad ? 1 : (fe_degree+1);
+ constexpr int n_blocks1 = dim > 1 ? n_rows : 1;
+ constexpr int n_blocks2 = dim > 2 ? n_rows : 1;
- const int stride = Utilities::fixed_int_power<fe_degree+1,face_direction>::value;
+ AssertIndexRange (face_direction, dim);
+ constexpr int stride = Utilities::fixed_int_power<n_rows,face_direction>::value;
+ constexpr int out_stride = Utilities::fixed_int_power<n_rows,dim-1>::value;
+ const Number *DEAL_II_RESTRICT shape_values = this->shape_values;
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
- if (dof_to_quad == true)
+ if (contract_onto_face == true)
{
- Number res0 = shape_data[0] * in[0];
- for (int ind=1; ind<mm; ++ind)
- res0 += shape_data[ind] * in[stride*ind];
+ Number res0 = shape_values[0] * in[0];
+ Number res1, res2;
+ if (max_derivative > 0)
+ res1 = shape_values[n_rows] * in[0];
+ if (max_derivative > 1)
+ res2 = shape_values[2*n_rows] * in[0];
+ for (int ind=1; ind<n_rows; ++ind)
+ {
+ res0 += shape_values[ind] * in[stride*ind];
+ if (max_derivative > 0)
+ res1 += shape_values[ind+n_rows] * in[stride*ind];
+ if (max_derivative > 1)
+ res2 += shape_values[ind+2*n_rows] * in[stride*ind];
+ }
if (add == false)
- out[0] = res0;
+ {
+ out[0] = res0;
+ if (max_derivative > 0)
+ out[out_stride] = res1;
+ if (max_derivative > 1)
+ out[2*out_stride] = res2;
+ }
else
- out[0] += res0;
+ {
+ out[0] += res0;
+ if (max_derivative > 0)
+ out[out_stride] += res1;
+ if (max_derivative > 1)
+ out[2*out_stride] += res2;
+ }
}
else
{
- for (int col=0; col<nn; ++col)
- if (add == false)
- out[col*stride] = shape_data[col] * in[0];
- else
- out[col*stride] += shape_data[col] * in[0];
+ for (int col=0; col<n_rows; ++col)
+ {
+ if (add == false)
+ out[col*stride] = shape_values[col] * in[0];
+ else
+ out[col*stride] += shape_values[col] * in[0];
+ if (max_derivative > 0)
+ out[col*stride] += shape_values[col+n_rows] * in[out_stride];
+ if (max_derivative > 1)
+ out[col*stride] += shape_values[col+2*n_rows] * in[2*out_stride];
+ }
}
// increment: in regular case, just go to the next point in
switch (face_direction)
{
case 0:
- in += mm;
- out += nn;
+ in += contract_onto_face ? n_rows : 1;
+ out += contract_onto_face ? 1 : n_rows;
break;
case 1:
++in;
// product. Need to take that into account.
if (dim == 3)
{
- if (dof_to_quad)
- out += fe_degree;
+ if (contract_onto_face)
+ out += n_rows-1;
else
- in += fe_degree;
+ in += n_rows-1;
}
break;
case 2:
}
if (face_direction == 1 && dim == 3)
{
- in += mm*(mm-1);
- out += nn*(nn-1);
// adjust for local coordinate system zx
- if (dof_to_quad)
- out -= (fe_degree+1)*(fe_degree+1)-1;
+ if (contract_onto_face)
+ {
+ in += n_rows*(n_rows-1);
+ out -= n_rows*n_rows-1;
+ }
else
- in -= (fe_degree+1)*(fe_degree+1)-1;
+ {
+ out += n_rows*(n_rows-1);
+ in -= n_rows*n_rows-1;
+ }
}
}
}
/**
- * Internal evaluator for 1d-3d shape function using the tensor product form
- * of the basis functions. The same as above but without making use of
- * template arguments and rather variable loop bounds.
+ * Internal evaluator for shape function using the tensor product form
+ * of the basis functions. The same as the other templated class but
+ * without making use of template arguments and variable loop bounds
+ * instead.
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
*/
- template <int dim, typename Number>
- struct EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
+ template <int dim, typename Number, typename Number2>
+ struct EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
{
- static const unsigned int dofs_per_cell = numbers::invalid_unsigned_int;
- static const unsigned int n_q_points = numbers::invalid_unsigned_int;
+ static const unsigned int n_rows_of_product = numbers::invalid_unsigned_int;
+ static const unsigned int n_columns_of_product = numbers::invalid_unsigned_int;
/**
* Empty constructor. Does nothing. Be careful when using 'values' and
shape_values (nullptr),
shape_gradients (nullptr),
shape_hessians (nullptr),
- fe_degree (numbers::invalid_unsigned_int),
- n_q_points_1d (numbers::invalid_unsigned_int)
+ n_rows (numbers::invalid_unsigned_int),
+ n_columns (numbers::invalid_unsigned_int)
{}
/**
* Constructor, taking the data from ShapeInfo
*/
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int fe_degree,
- const unsigned int n_q_points_1d)
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int n_rows,
+ const unsigned int n_columns)
:
shape_values (shape_values.begin()),
shape_gradients (shape_gradients.begin()),
shape_hessians (shape_hessians.begin()),
- fe_degree (fe_degree),
- n_q_points_1d (n_q_points_1d)
- {}
+ n_rows (n_rows),
+ n_columns (n_columns)
+ {
+ // We can enter this function either for the apply() path that has
+ // n_rows * n_columns entries or for the apply_face() path that only has
+ // n_rows * 3 entries in the array. Since we cannot decide about the use
+ // we must allow for both here.
+ Assert(shape_values.size() == 0 ||
+ shape_values.size() == n_rows*n_columns ||
+ shape_values.size() == n_rows*3,
+ ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
+ Assert(shape_gradients.size() == 0 ||
+ shape_gradients.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
+ Assert(shape_hessians.size() == 0 ||
+ shape_hessians.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
+ }
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
values (const Number *in,
Number *out) const
{
- apply<direction,dof_to_quad,add>(shape_values, in, out);
+ apply<direction,contract_over_rows,add>(shape_values, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
gradients (const Number *in,
Number *out) const
{
- apply<direction,dof_to_quad,add>(shape_gradients, in, out);
+ apply<direction,contract_over_rows,add>(shape_gradients, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
hessians (const Number *in,
Number *out) const
{
- apply<direction,dof_to_quad,add>(shape_hessians, in, out);
+ apply<direction,contract_over_rows,add>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ values_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_values, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ gradients_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ hessians_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,true>(shape_hessians, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
- void apply (const Number *shape_data,
- const Number *in,
- Number *out) const;
+ template <int direction, bool contract_over_rows, bool add, bool one_line=false>
+ void apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out) const;
+
+ template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
+ void apply_face (const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT out) const;
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
- const unsigned int fe_degree;
- const unsigned int n_q_points_1d;
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
+ const unsigned int n_rows;
+ const unsigned int n_columns;
};
- // evaluates the given shape data in 1d-3d using the tensor product
- // form. does not use a particular layout of entries in the matrices
- // like the functions below and corresponds to a usual matrix-matrix
- // product
- template <int dim, typename Number>
- template <int direction, bool dof_to_quad, bool add>
+
+
+ template <int dim, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add, bool one_line>
inline
void
- EvaluatorTensorProduct<evaluate_general,dim,-1,0,Number>
- ::apply (const Number *shape_data,
- const Number *in,
- Number *out) const
+ EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
+ ::apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out) const
{
+ static_assert (one_line == false || direction==dim-1,
+ "Single-line evaluation only works for direction=dim-1.");
+ Assert(shape_data != nullptr,
+ ExcMessage("The given array shape_data must not be the null pointer!"));
+ Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
+ ExcMessage("In-place operation only supported for "
+ "n_rows==n_columns or single-line interpolation"));
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
+ const int mm = contract_over_rows ? n_rows : n_columns,
+ nn = contract_over_rows ? n_columns : n_rows;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = direction==0 ? 1 : Utilities::fixed_power<direction>(nn);
+ const int stride = direction==0 ? 1 : Utilities::fixed_power<direction>(n_columns);
+ const int n_blocks1 = one_line ? 1 : stride;
+ const int n_blocks2 = direction >= dim-1 ? 1 : Utilities::fixed_power<dim-direction-1>(n_rows);
+ Assert(n_rows <= 128, ExcNotImplemented());
for (int i2=0; i2<n_blocks2; ++i2)
{
for (int i1=0; i1<n_blocks1; ++i1)
{
+ Number x[129];
+ for (int i=0; i<mm; ++i)
+ x[i] = in[stride*i];
for (int col=0; col<nn; ++col)
{
- Number val0;
- if (dof_to_quad == true)
+ Number2 val0;
+ if (contract_over_rows == true)
val0 = shape_data[col];
else
- val0 = shape_data[col*n_q_points_1d];
- Number res0 = val0 * in[0];
- for (int ind=1; ind<mm; ++ind)
+ val0 = shape_data[col*n_columns];
+ Number res0 = val0 * x[0];
+ for (int i=1; i<mm; ++i)
{
- if (dof_to_quad == true)
- val0 = shape_data[ind*n_q_points_1d+col];
+ if (contract_over_rows == true)
+ val0 = shape_data[i*n_columns+col];
else
- val0 = shape_data[col*n_q_points_1d+ind];
- res0 += val0 * in[stride*ind];
+ val0 = shape_data[col*n_columns+i];
+ res0 += val0 * x[i];
}
if (add == false)
out[stride*col] = res0;
out[stride*col] += res0;
}
+ if (one_line == false)
+ {
+ ++in;
+ ++out;
+ }
+ }
+ if (one_line == false)
+ {
+ in += stride*(mm-1);
+ out += stride*(nn-1);
+ }
+ }
+ }
+
+
+
+ template <int dim, typename Number, typename Number2>
+ template <int face_direction, bool contract_onto_face, bool add, int max_derivative>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_general,dim,0,0,Number,Number2>
+ ::apply_face (const Number *DEAL_II_RESTRICT in,
+ Number *DEAL_II_RESTRICT out) const
+ {
+ Assert(shape_values != nullptr,
+ ExcMessage("The given array shape_data must not be the null pointer!"));
+ static_assert(dim > 0 && dim<4, "Only dim=1,2,3 supported");
+ const int n_blocks1 = dim > 1 ? n_rows : 1;
+ const int n_blocks2 = dim > 2 ? n_rows : 1;
+
+ AssertIndexRange (face_direction, dim);
+ const int stride = face_direction > 0 ? Utilities::fixed_power<face_direction>(n_rows) : 1;
+ const int out_stride = dim > 1 ? Utilities::fixed_power<dim-1>(n_rows) : 1;
+
+ for (int i2=0; i2<n_blocks2; ++i2)
+ {
+ for (int i1=0; i1<n_blocks1; ++i1)
+ {
+ if (contract_onto_face == true)
+ {
+ Number res0 = shape_values[0] * in[0];
+ Number res1, res2;
+ if (max_derivative > 0)
+ res1 = shape_values[n_rows] * in[0];
+ if (max_derivative > 1)
+ res2 = shape_values[2*n_rows] * in[0];
+ for (unsigned int ind=1; ind<n_rows; ++ind)
+ {
+ res0 += shape_values[ind] * in[stride*ind];
+ if (max_derivative > 0)
+ res1 += shape_values[ind+n_rows] * in[stride*ind];
+ if (max_derivative > 1)
+ res2 += shape_values[ind+2*n_rows] * in[stride*ind];
+ }
+ if (add == false)
+ {
+ out[0] = res0;
+ if (max_derivative > 0)
+ out[out_stride] = res1;
+ if (max_derivative > 1)
+ out[2*out_stride] = res2;
+ }
+ else
+ {
+ out[0] += res0;
+ if (max_derivative > 0)
+ out[out_stride] += res1;
+ if (max_derivative > 1)
+ out[2*out_stride] += res2;
+ }
+ }
+ else
+ {
+ for (unsigned int col=0; col<n_rows; ++col)
+ {
+ if (add == false)
+ out[col*stride] = shape_values[col] * in[0];
+ else
+ out[col*stride] += shape_values[col] * in[0];
+ if (max_derivative > 0)
+ out[col*stride] += shape_values[col+n_rows] * in[out_stride];
+ if (max_derivative > 1)
+ out[col*stride] += shape_values[col+2*n_rows] * in[2*out_stride];
+ }
+ }
+
// increment: in regular case, just go to the next point in
// x-direction. If we are at the end of one chunk in x-dir, need
// to jump over to the next layer in z-direction
- switch (direction)
+ switch (face_direction)
{
case 0:
- in += mm;
- out += nn;
+ in += contract_onto_face ? n_rows : 1;
+ out += contract_onto_face ? 1 : n_rows;
break;
case 1:
+ ++in;
+ ++out;
+ // faces 2 and 3 in 3D use local coordinate system zx, which
+ // is the other way around compared to the tensor
+ // product. Need to take that into account.
+ if (dim == 3)
+ {
+ if (contract_onto_face)
+ out += n_rows-1;
+ else
+ in += n_rows-1;
+ }
+ break;
case 2:
++in;
++out;
Assert (false, ExcNotImplemented());
}
}
- if (direction == 1)
+ if (face_direction == 1 && dim == 3)
{
- in += nn*(mm-1);
- out += nn*(nn-1);
+ // adjust for local coordinate system zx
+ if (contract_onto_face)
+ {
+ in += n_rows*(n_rows-1);
+ out -= n_rows*n_rows-1;
+ }
+ else
+ {
+ out += n_rows*(n_rows-1);
+ in -= n_rows*n_rows-1;
+ }
}
}
}
* tensor-product based elements for "symmetric" finite elements, i.e., when
* the shape functions are symmetric about 0.5 and the quadrature points
* are, too.
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
*/
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ struct EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
{
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+ static const unsigned int n_rows_of_product = Utilities::fixed_int_power<n_rows,dim>::value;
+ static const unsigned int n_columns_of_product = Utilities::fixed_int_power<n_columns,dim>::value;
/**
* Constructor, taking the data from ShapeInfo
*/
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
:
shape_values (shape_values.begin()),
shape_gradients (shape_gradients.begin()),
shape_hessians (shape_hessians.begin())
{
+ Assert(shape_values.size() == 0 ||
+ shape_values.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_values.size(), n_rows*n_columns));
+ Assert(shape_gradients.size() == 0 ||
+ shape_gradients.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_gradients.size(), n_rows*n_columns));
+ Assert(shape_hessians.size() == 0 ||
+ shape_hessians.size() == n_rows*n_columns,
+ ExcDimensionMismatch(shape_hessians.size(), n_rows*n_columns));
(void)dummy1;
(void)dummy2;
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
values (const Number in [],
Number out[]) const;
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
gradients (const Number in [],
Number out[]) const;
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
hessians (const Number in [],
Number out[]) const;
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
};
// In these matrices, we want to use avoid computations involving zeros and
// ones and in addition use the symmetry in entries to reduce the number of
// read operations.
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add>
inline
void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
::values (const Number in [],
Number out []) const
{
+ Assert (shape_values != nullptr, ExcNotInitialized());
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
- const int n_cols = nn / 2;
- const int mid = mm / 2;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns,
+ nn = contract_over_rows ? n_columns : n_rows;
+ constexpr int n_cols = nn / 2;
+ constexpr int mid = mm / 2;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = Utilities::fixed_int_power<nn,direction>::value;
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
for (int i2=0; i2<n_blocks2; ++i2)
{
{
for (int col=0; col<n_cols; ++col)
{
- Number val0, val1, in0, in1, res0, res1;
- if (dof_to_quad == true)
+ Number2 val0, val1;
+ Number in0, in1, res0, res1;
+ if (contract_over_rows == true)
{
val0 = shape_values[col];
val1 = shape_values[nn-1-col];
}
else
{
- val0 = shape_values[col*n_q_points_1d];
- val1 = shape_values[(col+1)*n_q_points_1d-1];
+ val0 = shape_values[col*n_columns];
+ val1 = shape_values[(col+1)*n_columns-1];
}
if (mid > 0)
{
res1 += val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
- val0 = shape_values[ind*n_q_points_1d+col];
- val1 = shape_values[ind*n_q_points_1d+nn-1-col];
+ val0 = shape_values[ind*n_columns+col];
+ val1 = shape_values[ind*n_columns+nn-1-col];
}
else
{
- val0 = shape_values[col*n_q_points_1d+ind];
- val1 = shape_values[(col+1)*n_q_points_1d-1-ind];
+ val0 = shape_values[col*n_columns+ind];
+ val1 = shape_values[(col+1)*n_columns-1-ind];
}
in0 = in[stride*ind];
in1 = in[stride*(mm-1-ind)];
}
else
res0 = res1 = Number();
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
if (mm % 2 == 1)
{
- val0 = shape_values[mid*n_q_points_1d+col];
- val1 = val0 * in[stride*mid];
- res0 += val1;
- res1 += val1;
+ val0 = shape_values[mid*n_columns+col];
+ in1 = val0 * in[stride*mid];
+ res0 += in1;
+ res1 += in1;
}
}
else
{
if (mm % 2 == 1 && nn % 2 == 0)
{
- val0 = shape_values[col*n_q_points_1d+mid];
- val1 = val0 * in[stride*mid];
- res0 += val1;
- res1 += val1;
+ val0 = shape_values[col*n_columns+mid];
+ in1 = val0 * in[stride*mid];
+ res0 += in1;
+ res1 += in1;
}
}
if (add == false)
out[stride*(nn-1-col)] += res1;
}
}
- if ( dof_to_quad == true && nn%2==1 && mm%2==1 )
+ if ( contract_over_rows == true && nn%2==1 && mm%2==1 )
{
if (add==false)
out[stride*n_cols] = in[stride*mid];
else
out[stride*n_cols] += in[stride*mid];
}
- else if (dof_to_quad == true && nn%2==1)
+ else if (contract_over_rows == true && nn%2==1)
{
Number res0;
- Number val0 = shape_values[n_cols];
+ Number2 val0 = shape_values[n_cols];
if (mid > 0)
{
- res0 = in[0] + in[stride*(mm-1)];
- res0 *= val0;
+ res0 = val0 * (in[0] + in[stride*(mm-1)]);
for (int ind=1; ind<mid; ++ind)
{
- val0 = shape_values[ind*n_q_points_1d+n_cols];
- Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
+ val0 = shape_values[ind*n_columns+n_cols];
+ res0 += val0 * (in[stride*ind] + in[stride*(mm-1-ind)]);
}
}
else
else
out[stride*n_cols] += res0;
}
- else if (dof_to_quad == false && nn%2 == 1)
+ else if (contract_over_rows == false && nn%2 == 1)
{
Number res0;
if (mid > 0)
{
- Number val0 = shape_values[n_cols*n_q_points_1d];
- res0 = in[0] + in[stride*(mm-1)];
- res0 *= val0;
+ Number2 val0 = shape_values[n_cols*n_columns];
+ res0 = val0 * (in[0] + in[stride*(mm-1)]);
for (int ind=1; ind<mid; ++ind)
{
- val0 = shape_values[n_cols*n_q_points_1d+ind];
- Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
+ val0 = shape_values[n_cols*n_columns+ind];
+ Number in1 = val0 * (in[stride*ind] + in[stride*(mm-1-ind)]);
+ res0 += in1;
}
if (mm % 2)
res0 += in[stride*mid];
out[stride*n_cols] += res0;
}
- // increment: in regular case, just go to the next point in
- // x-direction. If we are at the end of one chunk in x-dir, need to
- // jump over to the next layer in z-direction
- switch (direction)
- {
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
- }
- if (direction == 1)
- {
- in += nn*(mm-1);
- out += nn*(nn-1);
+ ++in;
+ ++out;
}
+ in += stride*(mm-1);
+ out += stride*(nn-1);
}
}
// In these matrices, we want to use avoid computations involving
// zeros and ones and in addition use the symmetry in entries to
// reduce the number of read operations.
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add>
inline
void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
::gradients (const Number in [],
Number out []) const
{
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
- const int n_cols = nn / 2;
- const int mid = mm / 2;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns,
+ nn = contract_over_rows ? n_columns : n_rows;
+ constexpr int n_cols = nn / 2;
+ constexpr int mid = mm / 2;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = Utilities::fixed_int_power<nn,direction>::value;
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
for (int i2=0; i2<n_blocks2; ++i2)
{
{
for (int col=0; col<n_cols; ++col)
{
- Number val0, val1, in0, in1, res0, res1;
- if (dof_to_quad == true)
+ Number2 val0, val1;
+ Number in0, in1, res0, res1;
+ if (contract_over_rows == true)
{
val0 = shape_gradients[col];
val1 = shape_gradients[nn-1-col];
}
else
{
- val0 = shape_gradients[col*n_q_points_1d];
- val1 = shape_gradients[(nn-col-1)*n_q_points_1d];
+ val0 = shape_gradients[col*n_columns];
+ val1 = shape_gradients[(nn-col-1)*n_columns];
}
if (mid > 0)
{
res1 -= val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
- val0 = shape_gradients[ind*n_q_points_1d+col];
- val1 = shape_gradients[ind*n_q_points_1d+nn-1-col];
+ val0 = shape_gradients[ind*n_columns+col];
+ val1 = shape_gradients[ind*n_columns+nn-1-col];
}
else
{
- val0 = shape_gradients[col*n_q_points_1d+ind];
- val1 = shape_gradients[(nn-col-1)*n_q_points_1d+ind];
+ val0 = shape_gradients[col*n_columns+ind];
+ val1 = shape_gradients[(nn-col-1)*n_columns+ind];
}
in0 = in[stride*ind];
in1 = in[stride*(mm-1-ind)];
res0 = res1 = Number();
if (mm % 2 == 1)
{
- if (dof_to_quad == true)
- val0 = shape_gradients[mid*n_q_points_1d+col];
+ if (contract_over_rows == true)
+ val0 = shape_gradients[mid*n_columns+col];
else
- val0 = shape_gradients[col*n_q_points_1d+mid];
- val1 = val0 * in[stride*mid];
- res0 += val1;
- res1 -= val1;
+ val0 = shape_gradients[col*n_columns+mid];
+ in1 = val0 * in[stride*mid];
+ res0 += in1;
+ res1 -= in1;
}
if (add == false)
{
}
if ( nn%2 == 1 )
{
- Number val0, res0;
- if (dof_to_quad == true)
+ Number2 val0;
+ Number res0;
+ if (contract_over_rows == true)
val0 = shape_gradients[n_cols];
else
- val0 = shape_gradients[n_cols*n_q_points_1d];
- res0 = in[0] - in[stride*(mm-1)];
- res0 *= val0;
+ val0 = shape_gradients[n_cols*n_columns];
+ res0 = val0 * (in[0] - in[stride*(mm-1)]);
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
- val0 = shape_gradients[ind*n_q_points_1d+n_cols];
+ if (contract_over_rows == true)
+ val0 = shape_gradients[ind*n_columns+n_cols];
else
- val0 = shape_gradients[n_cols*n_q_points_1d+ind];
- Number val1 = in[stride*ind] - in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
+ val0 = shape_gradients[n_cols*n_columns+ind];
+ Number in1 = val0 * (in[stride*ind] - in[stride*(mm-1-ind)]);
+ res0 += in1;
}
if (add == false)
out[stride*n_cols] = res0;
out[stride*n_cols] += res0;
}
- // increment: in regular case, just go to the next point in
- // x-direction. for y-part in 3D and if we are at the end of one
- // chunk in x-dir, need to jump over to the next layer in
- // z-direction
- switch (direction)
- {
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
- }
-
- if (direction == 1)
- {
- in += nn * (mm-1);
- out += nn * (nn-1);
+ ++in;
+ ++out;
}
+ in += stride*(mm-1);
+ out += stride*(nn-1);
}
}
// evaluates the given shape data in 1d-3d using the tensor product
// form assuming the symmetries of unit cell shape hessians for
// finite elements in FEEvaluation
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add>
inline
void
- EvaluatorTensorProduct<evaluate_symmetric,dim,fe_degree,n_q_points_1d,Number>
+ EvaluatorTensorProduct<evaluate_symmetric,dim,n_rows,n_columns,Number,Number2>
::hessians (const Number in [],
Number out []) const
{
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
- const int n_cols = nn / 2;
- const int mid = mm / 2;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns;
+ constexpr int nn = contract_over_rows ? n_columns : n_rows;
+ constexpr int n_cols = nn / 2;
+ constexpr int mid = mm / 2;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = Utilities::fixed_int_power<nn,direction>::value;
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
for (int i2=0; i2<n_blocks2; ++i2)
{
{
for (int col=0; col<n_cols; ++col)
{
- Number val0, val1, in0, in1, res0, res1;
- if (dof_to_quad == true)
+ Number2 val0, val1;
+ Number in0, in1, res0, res1;
+ if (contract_over_rows == true)
{
val0 = shape_hessians[col];
val1 = shape_hessians[nn-1-col];
}
else
{
- val0 = shape_hessians[col*n_q_points_1d];
- val1 = shape_hessians[(col+1)*n_q_points_1d-1];
+ val0 = shape_hessians[col*n_columns];
+ val1 = shape_hessians[(col+1)*n_columns-1];
}
if (mid > 0)
{
res1 += val0 * in1;
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
- val0 = shape_hessians[ind*n_q_points_1d+col];
- val1 = shape_hessians[ind*n_q_points_1d+nn-1-col];
+ val0 = shape_hessians[ind*n_columns+col];
+ val1 = shape_hessians[ind*n_columns+nn-1-col];
}
else
{
- val0 = shape_hessians[col*n_q_points_1d+ind];
- val1 = shape_hessians[(col+1)*n_q_points_1d-1-ind];
+ val0 = shape_hessians[col*n_columns+ind];
+ val1 = shape_hessians[(col+1)*n_columns-1-ind];
}
in0 = in[stride*ind];
in1 = in[stride*(mm-1-ind)];
res0 = res1 = Number();
if (mm % 2 == 1)
{
- if (dof_to_quad == true)
- val0 = shape_hessians[mid*n_q_points_1d+col];
+ if (contract_over_rows == true)
+ val0 = shape_hessians[mid*n_columns+col];
else
- val0 = shape_hessians[col*n_q_points_1d+mid];
- val1 = val0 * in[stride*mid];
- res0 += val1;
- res1 += val1;
+ val0 = shape_hessians[col*n_columns+mid];
+ in1 = val0 * in[stride*mid];
+ res0 += in1;
+ res1 += in1;
}
if (add == false)
{
}
if ( nn%2 == 1 )
{
- Number val0, res0;
- if (dof_to_quad == true)
+ Number2 val0;
+ Number res0;
+ if (contract_over_rows == true)
val0 = shape_hessians[n_cols];
else
- val0 = shape_hessians[n_cols*n_q_points_1d];
+ val0 = shape_hessians[n_cols*n_columns];
if (mid > 0)
{
- res0 = in[0] + in[stride*(mm-1)];
- res0 *= val0;
+ res0 = val0 * (in[0] + in[stride*(mm-1)]);
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
- val0 = shape_hessians[ind*n_q_points_1d+n_cols];
+ if (contract_over_rows == true)
+ val0 = shape_hessians[ind*n_columns+n_cols];
else
- val0 = shape_hessians[n_cols*n_q_points_1d+ind];
- Number val1 = in[stride*ind] + in[stride*(mm-1-ind)];
- val1 *= val0;
- res0 += val1;
+ val0 = shape_hessians[n_cols*n_columns+ind];
+ Number in1 = val0*(in[stride*ind] + in[stride*(mm-1-ind)]);
+ res0 += in1;
}
}
else
res0 = Number();
if (mm % 2 == 1)
{
- if (dof_to_quad == true)
- val0 = shape_hessians[mid*n_q_points_1d+n_cols];
+ if (contract_over_rows == true)
+ val0 = shape_hessians[mid*n_columns+n_cols];
else
- val0 = shape_hessians[n_cols*n_q_points_1d+mid];
+ val0 = shape_hessians[n_cols*n_columns+mid];
res0 += val0 * in[stride*mid];
}
if (add == false)
out[stride*n_cols] += res0;
}
- // increment: in regular case, just go to the next point in
- // x-direction. If we are at the end of one chunk in x-dir, need to
- // jump over to the next layer in z-direction
- switch (direction)
- {
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
- }
- }
- if (direction == 1)
- {
- in += nn*(mm-1);
- out += nn*(nn-1);
+ ++in;
+ ++out;
}
+ in += stride*(mm-1);
+ out += stride*(nn-1);
}
}
* experiments in the book say that the method is not efficient for N<20, it
* is more efficient in the context where the loop bounds are compile-time
* constants (templates).
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
*/
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- struct EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ struct EvaluatorTensorProduct<evaluate_evenodd,dim,n_rows,n_columns,Number,Number2>
{
- static const unsigned int dofs_per_cell = Utilities::fixed_int_power<fe_degree+1,dim>::value;
- static const unsigned int n_q_points = Utilities::fixed_int_power<n_q_points_1d,dim>::value;
+ static const unsigned int n_rows_of_product = Utilities::fixed_int_power<n_rows,dim>::value;
+ static const unsigned int n_columns_of_product = Utilities::fixed_int_power<n_columns,dim>::value;
/**
* Empty constructor. Does nothing. Be careful when using 'values' and
- * related methods because they need to be filled with the other pointer
+ * related methods because they need to be filled with the other
+ * constructor passing in at least an array for the values.
*/
EvaluatorTensorProduct ()
:
- shape_values (0),
- shape_gradients (0),
- shape_hessians (0)
+ shape_values (nullptr),
+ shape_gradients (nullptr),
+ shape_hessians (nullptr)
{}
/**
* Constructor, taking the data from ShapeInfo (using the even-odd
* variants stored there)
*/
- EvaluatorTensorProduct (const AlignedVector<Number> &shape_values,
- const AlignedVector<Number> &shape_gradients,
- const AlignedVector<Number> &shape_hessians,
- const unsigned int dummy1 = 0,
- const unsigned int dummy2 = 0)
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (nullptr),
+ shape_hessians (nullptr)
+ {
+ AssertDimension(shape_values.size(), n_rows*((n_columns+1)/2));
+ }
+
+ /**
+ * Constructor, taking the data from ShapeInfo (using the even-odd
+ * variants stored there)
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
:
shape_values (shape_values.begin()),
shape_gradients (shape_gradients.begin()),
shape_hessians (shape_hessians.begin())
{
+ // In this function, we allow for dummy pointers if some of values,
+ // gradients or hessians should not be computed
+ if (!shape_values.empty())
+ AssertDimension(shape_values.size(), n_rows*((n_columns+1)/2));
+ if (!shape_gradients.empty())
+ AssertDimension(shape_gradients.size(), n_rows*((n_columns+1)/2));
+ if (!shape_hessians.empty())
+ AssertDimension(shape_hessians.size(), n_rows*((n_columns+1)/2));
(void)dummy1;
(void)dummy2;
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
values (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add,0>(shape_values, in, out);
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0>(shape_values, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
gradients (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add,1>(shape_gradients, in, out);
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,1>(shape_gradients, in, out);
}
- template <int direction, bool dof_to_quad, bool add>
+ template <int direction, bool contract_over_rows, bool add>
void
hessians (const Number in [],
Number out[]) const
{
- apply<direction,dof_to_quad,add,2>(shape_hessians, in, out);
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,2>(shape_hessians, in, out);
}
- template <int direction, bool dof_to_quad, bool add, int type>
- static void apply (const Number *shape_data,
- const Number in [],
- Number out []);
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ values_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0,true>(shape_values, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ gradients_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,1,true>(shape_gradients, in, out);
+ }
- const Number *shape_values;
- const Number *shape_gradients;
- const Number *shape_hessians;
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ hessians_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,2,true>(shape_hessians, in, out);
+ }
+
+ /**
+ * This function applies the tensor product kernel, corresponding to a
+ * multiplication of 1D stripes, along the given @p direction of the tensor
+ * data in the input array. This function allows the @p in and @p out
+ * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
+ * perform the contraction in place where @p in and @p out point to the
+ * same address. For the case n_rows != n_columns, the output is only
+ * correct if @p one_line is set to true.
+ *
+ * @tparam direction Direction that is evaluated
+ * @tparam contract_over_rows If true, the tensor contraction sums
+ * over the rows in the given @p shape_data
+ * array, otherwise it sums over the columns
+ * @tparam add If true, the result is added to the output vector, else
+ * the computed values overwrite the content in the output
+ * @tparam type Determines whether to use the symmetries appearing in
+ * shape values (type=0), shape gradients (type=1) or
+ * second derivatives (type=2, similar to type 0 but
+ * without two additional zero entries)
+ * @tparam one_line If true, the kernel is only applied along a single 1D
+ * stripe within a dim-dimensional tensor, not the full
+ * n_rows^dim points as in the @p false case.
+ *
+ * @param shape_data Transformation matrix with @p n_rows rows and
+ * @p n_columns columns, stored in row-major format
+ * @param in Pointer to the start of the input data vector
+ * @param out Pointer to the start of the output data vector
+ */
+ template <int direction, bool contract_over_rows, bool add, int type,
+ bool one_line=false>
+ static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out);
+
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
};
- template <int dim, int fe_degree, int n_q_points_1d, typename Number>
- template <int direction, bool dof_to_quad, bool add, int type>
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add, int type, bool one_line>
inline
void
- EvaluatorTensorProduct<evaluate_evenodd,dim,fe_degree,n_q_points_1d,Number>
- ::apply (const Number *shapes,
- const Number in [],
- Number out [])
+ EvaluatorTensorProduct<evaluate_evenodd,dim,n_rows,n_columns,Number,Number2>
+ ::apply (const Number2 *DEAL_II_RESTRICT shapes,
+ const Number *in,
+ Number *out)
{
- AssertIndexRange (type, 3);
+ static_assert (type < 3, "Only three variants type=0,1,2 implemented");
+ static_assert (one_line == false || direction==dim-1,
+ "Single-line evaluation only works for direction=dim-1.");
+ Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
+ ExcMessage("In-place operation only supported for "
+ "n_rows==n_columns or single-line interpolation"));
+
+ // We cannot statically assert that direction is less than dim, so must do
+ // an additional dynamic check
AssertIndexRange (direction, dim);
- const int mm = dof_to_quad ? (fe_degree+1) : n_q_points_1d,
- nn = dof_to_quad ? n_q_points_1d : (fe_degree+1);
- const int n_cols = nn / 2;
- const int mid = mm / 2;
- const int n_blocks1 = (dim > 1 ? (direction > 0 ? nn : mm) : 1);
- const int n_blocks2 = (dim > 2 ? (direction > 1 ? nn : mm) : 1);
- const int stride = Utilities::fixed_int_power<nn,direction>::value;
+ constexpr int nn = contract_over_rows ? n_columns : n_rows;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns;
+ constexpr int n_cols = nn / 2;
+ constexpr int mid = mm / 2;
- const int offset = (n_q_points_1d+1)/2;
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = one_line ? 1 : stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
+
+ constexpr int offset = (n_columns+1)/2;
// this code may look very inefficient at first sight due to the many
// different cases with if's at the innermost loop part, but all of the
Number xp[mid>0?mid:1], xm[mid>0?mid:1];
for (int i=0; i<mid; ++i)
{
- if (dof_to_quad == true && type == 1)
+ if (contract_over_rows == true && type == 1)
{
xp[i] = in[stride*i] - in[stride*(mm-1-i)];
xm[i] = in[stride*i] + in[stride*(mm-1-i)];
xm[i] = in[stride*i] - in[stride*(mm-1-i)];
}
}
+ Number xmid = in[stride*mid];
for (int col=0; col<n_cols; ++col)
{
Number r0, r1;
if (mid > 0)
{
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
- r0 = shapes[col] * xp[0];
- r1 = shapes[fe_degree*offset + col] * xm[0];
+ r0 = shapes[col] * xp[0];
+ r1 = shapes[(n_rows-1)*offset + col] * xm[0];
}
else
{
- r0 = shapes[col*offset] * xp[0];
- r1 = shapes[(fe_degree-col)*offset] * xm[0];
+ r0 = shapes[col*offset] * xp[0];
+ r1 = shapes[(n_rows-1-col)*offset] * xm[0];
}
for (int ind=1; ind<mid; ++ind)
{
- if (dof_to_quad == true)
+ if (contract_over_rows == true)
{
- r0 += shapes[ind*offset+col] * xp[ind];
- r1 += shapes[(fe_degree-ind)*offset+col] * xm[ind];
+ r0 += shapes[ind*offset+col] * xp[ind];
+ r1 += shapes[(n_rows-1-ind)*offset+col] * xm[ind];
}
else
{
- r0 += shapes[col*offset+ind] * xp[ind];
- r1 += shapes[(fe_degree-col)*offset+ind] * xm[ind];
+ r0 += shapes[col*offset+ind] * xp[ind];
+ r1 += shapes[(n_rows-1-col)*offset+ind] * xm[ind];
}
}
}
else
r0 = r1 = Number();
- if (mm % 2 == 1 && dof_to_quad == true)
+ if (mm % 2 == 1 && contract_over_rows == true)
{
if (type == 1)
- r1 += shapes[mid*offset+col] * in[stride*mid];
+ r1 += shapes[mid*offset+col] * xmid;
else
- r0 += shapes[mid*offset+col] * in[stride*mid];
+ r0 += shapes[mid*offset+col] * xmid;
}
else if (mm % 2 == 1 && (nn % 2 == 0 || type > 0))
- r0 += shapes[col*offset+mid] * in[stride*mid];
+ r0 += shapes[col*offset+mid] * xmid;
if (add == false)
{
out[stride*col] = r0 + r1;
- if (type == 1 && dof_to_quad == false)
+ if (type == 1 && contract_over_rows == false)
out[stride*(nn-1-col)] = r1 - r0;
else
out[stride*(nn-1-col)] = r0 - r1;
else
{
out[stride*col] += r0 + r1;
- if (type == 1 && dof_to_quad == false)
+ if (type == 1 && contract_over_rows == false)
out[stride*(nn-1-col)] += r1 - r0;
else
out[stride*(nn-1-col)] += r0 - r1;
}
}
- if ( type == 0 && dof_to_quad == true && nn%2==1 && mm%2==1 )
+ if ( type == 0 && contract_over_rows == true && nn%2==1 && mm%2==1 )
{
if (add==false)
- out[stride*n_cols] = in[stride*mid];
+ out[stride*n_cols] = xmid;
else
- out[stride*n_cols] += in[stride*mid];
+ out[stride*n_cols] += xmid;
}
- else if (dof_to_quad == true && nn%2==1)
+ else if (contract_over_rows == true && nn%2==1)
{
Number r0;
if (mid > 0)
else
r0 = Number();
if (type != 1 && mm % 2 == 1)
- r0 += shapes[mid*offset+n_cols] * in[stride*mid];
+ r0 += shapes[mid*offset+n_cols] * xmid;
if (add == false)
out[stride*n_cols] = r0;
else
out[stride*n_cols] += r0;
}
- else if (dof_to_quad == false && nn%2 == 1)
+ else if (contract_over_rows == false && nn%2 == 1)
{
Number r0;
if (mid > 0)
r0 = Number();
if (type == 0 && mm % 2 == 1)
- r0 += in[stride*mid];
+ r0 += xmid;
else if (type == 2 && mm % 2 == 1)
- r0 += shapes[n_cols*offset+mid] * in[stride*mid];
+ r0 += shapes[n_cols*offset+mid] * xmid;
if (add == false)
out[stride*n_cols] = r0;
else
out[stride*n_cols] += r0;
}
+ if (one_line == false)
+ {
+ in += 1;
+ out += 1;
+ }
+ }
+ if (one_line == false)
+ {
+ in += stride * (mm-1);
+ out += stride * (nn-1);
+ }
+ }
+ }
- // increment: in regular case, just go to the next point in
- // x-direction. If we are at the end of one chunk in x-dir, need to
- // jump over to the next layer in z-direction
- switch (direction)
+
+
+ /**
+ * Internal evaluator for 1d-3d shape function using the tensor product form
+ * of the basis functions.
+ *
+ * This class implements an approach similar to the even-odd decomposition
+ * but with a different type of symmetry. In this case, we assume that a
+ * single shape function already shows the symmetry over the quadrature
+ * points, rather than the complete basis that is considered in the even-odd
+ * case. In particular, we assume that the shape functions are ordered as in
+ * the Legendre basis, with symmetric shape functions in the even slots
+ * (rows of the values array) and point-symmetric in the odd slots. Like the
+ * even-odd decomposition, the number of operations are N^2/2 rather than
+ * N^2 FMAs (fused multiply-add), where N is the number of 1D dofs. The
+ * difference is in the way the input and output quantities are symmetrized.
+ *
+ * @tparam dim Space dimension in which this class is applied
+ * @tparam n_rows Number of rows in the transformation matrix, which corresponds
+ * to the number of 1d shape functions in the usual tensor
+ * contraction setting
+ * @tparam n_columns Number of columns in the transformation matrix, which
+ * corresponds to the number of 1d shape functions in the
+ * usual tensor contraction setting
+ * @tparam Number Abstract number type for input and output arrays
+ * @tparam Number2 Abstract number type for coefficient arrays (defaults to
+ * same type as the input/output arrays); must implement
+ * operator* with Number and produce Number as an output to
+ * be a valid type
+ */
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ struct EvaluatorTensorProduct<evaluate_symmetric_hierarchical,dim,n_rows,n_columns,Number,Number2>
+ {
+ static const unsigned int n_rows_of_product = Utilities::fixed_int_power<n_rows,dim>::value;
+ static const unsigned int n_columns_of_product = Utilities::fixed_int_power<n_columns,dim>::value;
+
+ /**
+ * Empty constructor. Does nothing. Be careful when using 'values' and
+ * related methods because they need to be filled with the other
+ * constructor passing in at least an array for the values.
+ */
+ EvaluatorTensorProduct ()
+ :
+ shape_values (nullptr),
+ shape_gradients (nullptr),
+ shape_hessians (nullptr)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo (using the even-odd
+ * variants stored there)
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number> &shape_values)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (nullptr),
+ shape_hessians (nullptr)
+ {}
+
+ /**
+ * Constructor, taking the data from ShapeInfo (using the even-odd
+ * variants stored there)
+ */
+ EvaluatorTensorProduct (const AlignedVector<Number2> &shape_values,
+ const AlignedVector<Number2> &shape_gradients,
+ const AlignedVector<Number2> &shape_hessians,
+ const unsigned int dummy1 = 0,
+ const unsigned int dummy2 = 0)
+ :
+ shape_values (shape_values.begin()),
+ shape_gradients (shape_gradients.begin()),
+ shape_hessians (shape_hessians.begin())
+ {
+ (void)dummy1;
+ (void)dummy2;
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ values (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0>(shape_values, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ gradients (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,1>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ hessians (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0>(shape_hessians, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ values_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_values != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0,true>(shape_values, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ gradients_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_gradients != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,1,true>(shape_gradients, in, out);
+ }
+
+ template <int direction, bool contract_over_rows, bool add>
+ void
+ hessians_one_line (const Number in [],
+ Number out[]) const
+ {
+ Assert (shape_hessians != nullptr, ExcNotInitialized());
+ apply<direction,contract_over_rows,add,0,true>(shape_hessians, in, out);
+ }
+
+ /**
+ * This function applies the tensor product kernel, corresponding to a
+ * multiplication of 1D stripes, along the given @p direction of the tensor
+ * data in the input array. This function allows the @p in and @p out
+ * arrays to alias for the case n_rows == n_columns, i.e., it is safe to
+ * perform the contraction in place where @p in and @p out point to the
+ * same address. For the case n_rows != n_columns, the output is only
+ * correct if @p one_line is set to true.
+ *
+ * @tparam direction Direction that is evaluated
+ * @tparam contract_over_rows If true, the tensor contraction sums
+ * over the rows in the given @p shape_data
+ * array, otherwise it sums over the columns
+ * @tparam add If true, the result is added to the output vector, else
+ * the computed values overwrite the content in the output
+ * @tparam type Determines whether the evaluation is symmetric in even
+ * rows (type=0) or odd rows (type=1) of @p shape_data and
+ * skew-symmetric in odd rows (type=0) or even rows (type=1)
+ * @tparam one_line If true, the kernel is only applied along a single 1D
+ * stripe within a dim-dimensional tensor, not the full
+ * n_rows^dim points as in the @p false case.
+ *
+ * @param shape_data Transformation matrix with @p n_rows rows and
+ * @p n_columns columns, stored in row-major format
+ * @param in Pointer to the start of the input data vector
+ * @param out Pointer to the start of the output data vector
+ */
+ template <int direction, bool contract_over_rows, bool add, int type,
+ bool one_line=false>
+ static void apply (const Number2 *DEAL_II_RESTRICT shape_data,
+ const Number *in,
+ Number *out);
+
+ const Number2 *shape_values;
+ const Number2 *shape_gradients;
+ const Number2 *shape_hessians;
+ };
+
+
+
+ template <int dim, int n_rows, int n_columns, typename Number, typename Number2>
+ template <int direction, bool contract_over_rows, bool add, int type, bool one_line>
+ inline
+ void
+ EvaluatorTensorProduct<evaluate_symmetric_hierarchical,dim,n_rows,n_columns,Number,Number2>
+ ::apply (const Number2 *DEAL_II_RESTRICT shapes,
+ const Number *in,
+ Number *out)
+ {
+ static_assert (one_line == false || direction==dim-1,
+ "Single-line evaluation only works for direction=dim-1.");
+ static_assert (type == 0 || type == 1,
+ "Only types 0 and 1 implemented for evaluate_symmetric_hierarchical.");
+ Assert (dim == direction+1 || one_line == true || n_rows == n_columns || in != out,
+ ExcMessage("In-place operation only supported for "
+ "n_rows==n_columns or single-line interpolation"));
+
+ // We cannot statically assert that direction is less than dim, so must do
+ // an additional dynamic check
+ AssertIndexRange (direction, dim);
+
+ constexpr int nn = contract_over_rows ? n_columns : n_rows;
+ constexpr int mm = contract_over_rows ? n_rows : n_columns;
+ constexpr int n_cols = nn / 2;
+ constexpr int mid = mm / 2;
+
+ constexpr int stride = Utilities::fixed_int_power<n_columns,direction>::value;
+ constexpr int n_blocks1 = one_line ? 1 : stride;
+ constexpr int n_blocks2 = Utilities::fixed_int_power<n_rows,(direction>=dim)?0:(dim-direction-1)>::value;
+
+ // this code may look very inefficient at first sight due to the many
+ // different cases with if's at the innermost loop part, but all of the
+ // conditionals can be evaluated at compile time because they are
+ // templates, so the compiler should optimize everything away
+ for (int i2=0; i2<n_blocks2; ++i2)
+ {
+ for (int i1=0; i1<n_blocks1; ++i1)
+ {
+ if (contract_over_rows)
{
- case 0:
- in += mm;
- out += nn;
- break;
- case 1:
- case 2:
- ++in;
- ++out;
- break;
- default:
- Assert (false, ExcNotImplemented());
+ Number x[mm];
+ for (unsigned int i=0; i<mm; ++i)
+ x[i] = in[stride*i];
+ for (unsigned int col=0; col<n_cols; ++col)
+ {
+ Number r0, r1;
+ if (mid > 0)
+ {
+ r0 = shapes[col] * x[0];
+ r1 = shapes[col+n_columns] * x[1];
+ for (unsigned int ind=1; ind<mid; ++ind)
+ {
+ r0 += shapes[col+2*ind*n_columns] * x[2*ind];
+ r1 += shapes[col+(2*ind+1)*n_columns] * x[2*ind+1];
+ }
+ }
+ else
+ r0 = r1 = Number();
+ if (mm%2 == 1)
+ r0 += shapes[col+(mm-1)*n_columns] * x[mm-1];
+ if (add == false)
+ {
+ out[stride*col] = r0 + r1;
+ if (type == 1)
+ out[stride*(nn-1-col)] = r1 - r0;
+ else
+ out[stride*(nn-1-col)] = r0 - r1;
+ }
+ else
+ {
+ out[stride*col] += r0 + r1;
+ if (type == 1)
+ out[stride*(nn-1-col)] += r1 - r0;
+ else
+ out[stride*(nn-1-col)] += r0 - r1;
+ }
+ }
+ if (nn%2 == 1)
+ {
+ Number r0;
+ const unsigned int shift = type==1 ? 1 : 0;
+ if (mid>0)
+ {
+ r0 = shapes[n_cols + shift*n_columns] * x[shift];
+ for (unsigned int ind=1; ind<mid; ++ind)
+ r0 += shapes[n_cols + (2*ind+shift)*n_columns] * x[2*ind+shift];
+ }
+ else
+ r0 = 0;
+ if (type != 1 && mm%2 == 1)
+ r0 += shapes[n_cols + (mm-1)*n_columns] * x[mm-1];
+ if (add == false)
+ out[stride*n_cols] = r0;
+ else
+ out[stride*n_cols] += r0;
+ }
+ }
+ else
+ {
+ Number xp[mid+1], xm[mid>0?mid:1];
+ for (int i=0; i<mid; ++i)
+ if (type == 0)
+ {
+ xp[i] = in[stride*i] + in[stride*(mm-1-i)];
+ xm[i] = in[stride*i] - in[stride*(mm-1-i)];
+ }
+ else
+ {
+ xp[i] = in[stride*i] - in[stride*(mm-1-i)];
+ xm[i] = in[stride*i] + in[stride*(mm-1-i)];
+ }
+ if (mm%2 == 1)
+ xp[mid] = in[stride*mid];
+ for (unsigned int col=0; col<n_cols; ++col)
+ {
+ Number r0, r1;
+ if (mid > 0)
+ {
+ r0 = shapes[2*col*n_columns] * xp[0];
+ r1 = shapes[(2*col+1)*n_columns] * xm[0];
+ for (unsigned int ind=1; ind<mid; ++ind)
+ {
+ r0 += shapes[2*col*n_columns+ind] * xp[ind];
+ r1 += shapes[(2*col+1)*n_columns+ind] * xm[ind];
+ }
+ }
+ else
+ r0 = r1 = Number();
+ if (mm%2 == 1)
+ {
+ if (type == 1)
+ r1 += shapes[(2*col+1)*n_columns+mid] * xp[mid];
+ else
+ r0 += shapes[2*col*n_columns+mid] * xp[mid];
+ }
+ if (add == false)
+ {
+ out[stride*(2*col)] = r0;
+ out[stride*(2*col+1)] = r1;
+ }
+ else
+ {
+ out[stride*(2*col)] += r0;
+ out[stride*(2*col+1)] += r1;
+ }
+ }
+ if (nn%2 == 1)
+ {
+ Number r0;
+ if (mid > 0)
+ {
+ r0 = shapes[(nn-1)*n_columns] * xp[0];
+ for (unsigned int ind=1; ind<mid; ++ind)
+ r0 += shapes[(nn-1)*n_columns+ind] * xp[ind];
+ }
+ else
+ r0 = Number();
+ if (mm%2 == 1 && type == 0)
+ r0 += shapes[(nn-1)*n_columns+mid] * xp[mid];
+ if (add == false)
+ out[stride*(nn-1)] = r0;
+ else
+ out[stride*(nn-1)] += r0;
+ }
+ }
+ if (one_line == false)
+ {
+ in += 1;
+ out += 1;
}
}
- if (direction == 1)
+ if (one_line == false)
{
- in += nn*(mm-1);
- out += nn*(nn-1);
+ in += stride * (mm-1);
+ out += stride * (nn-1);
}
}
}
const unsigned int n_comp = 1+ (spacedim-1)/vec_length;
const unsigned int n_hessians = (dim*(dim+1))/2;
- VectorizedArray<double> *values_dofs_ptr[n_comp];
- VectorizedArray<double> *values_quad_ptr[n_comp];
- VectorizedArray<double> *gradients_quad_ptr[n_comp][dim];
- VectorizedArray<double> *hessians_quad_ptr[n_comp][n_hessians];
-
const bool evaluate_values = update_flags & update_quadrature_points;
const bool evaluate_gradients= (cell_similarity != CellSimilarity::translation)
&&(update_flags & update_contravariant_transformation);
data.values_quad.resize(n_comp*n_q_points);
data.gradients_quad.resize (n_comp*n_q_points*dim);
+ if (evaluate_hessians)
+ data.hessians_quad.resize(n_comp*n_q_points*n_hessians);
+
const std::vector<unsigned int> &renumber_to_lexicographic
= data.shape_info.lexicographic_numbering;
for (unsigned int i=0; i<n_shape_values; ++i)
= data.mapping_support_points[renumber_to_lexicographic[i]][d];
}
- for (unsigned int c=0; c<n_comp; ++c)
- {
- values_dofs_ptr[c] = &(data.values_dofs[c*n_shape_values]);
- values_quad_ptr[c] = &(data.values_quad[c*n_q_points]);
- for (unsigned int j=0; j<dim; ++j)
- gradients_quad_ptr[c][j] = &(data.gradients_quad[(c*dim+j)*n_q_points]);
- }
-
- if (evaluate_hessians)
- {
- data.hessians_quad.resize(n_comp*n_q_points*n_hessians);
- for (unsigned int c=0; c<n_comp; ++c)
- for (unsigned int j=0; j<n_hessians; ++j)
- hessians_quad_ptr[c][j] = &(data.hessians_quad[(c*n_hessians+j)*n_q_points]);
- }
-
// do the actual tensorized evaluation
- SelectEvaluator<dim, -1, 0, n_comp, double>::evaluate
- (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]),
- &(gradients_quad_ptr[0]), &(hessians_quad_ptr[0]), &(data.scratch[0]),
+ SelectEvaluator<dim, -1, 0, n_comp, VectorizedArray<double> >::evaluate
+ (data.shape_info, data.values_dofs.begin(), data.values_quad.begin(),
+ data.gradients_quad.begin(), data.hessians_quad.begin(), data.scratch.begin(),
evaluate_values, evaluate_gradients, evaluate_hessians);
}
{
template
void
- SelectEvaluator<deal_II_dimension, -1, 0, components, scalar_type>::integrate
+ SelectEvaluator<deal_II_dimension, -1, 0, components, VectorizedArray<scalar_type> >::integrate
(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<scalar_type> > &shape_info,
- VectorizedArray<scalar_type> *[], VectorizedArray<scalar_type> *[],
- VectorizedArray<scalar_type> *[][deal_II_dimension], VectorizedArray<scalar_type> *,
+ VectorizedArray<scalar_type> *, VectorizedArray<scalar_type> *,
+ VectorizedArray<scalar_type> *, VectorizedArray<scalar_type> *,
const bool, const bool);
template
void
- SelectEvaluator<deal_II_dimension, -1, 0, components, scalar_type>::evaluate
+ SelectEvaluator<deal_II_dimension, -1, 0, components, VectorizedArray<scalar_type> >::evaluate
(const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArray<scalar_type> > &shape_info,
- VectorizedArray<scalar_type> *[], VectorizedArray<scalar_type> *[],
- VectorizedArray<scalar_type> *[][deal_II_dimension],
- VectorizedArray<scalar_type> *[][(deal_II_dimension*(deal_II_dimension+1))/2],
+ VectorizedArray<scalar_type> *, VectorizedArray<scalar_type> *,
+ VectorizedArray<scalar_type> *,
+ VectorizedArray<scalar_type> *,
VectorizedArray<scalar_type> *, const bool, const bool, const bool);
}
#include <deal.II/multigrid/mg_transfer_matrix_free.h>
#include <deal.II/multigrid/mg_transfer_internal.h>
-#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/matrix_free/evaluation_kernels.h>
#include <algorithm>
}
}
- evaluation_data.resize(3*n_child_cell_dofs);
+ evaluation_data.resize(n_child_cell_dofs);
}
namespace
{
- template <int dim, typename Eval, typename Number, bool prolongate>
- void
- perform_tensorized_op(const Eval &evaluator,
- const unsigned int n_points_in,
- const unsigned int n_child_cell_dofs,
- const unsigned int n_components,
- AlignedVector<VectorizedArray<Number> > &evaluation_data)
- {
- if (Eval::n_q_points != numbers::invalid_unsigned_int)
- AssertDimension(n_components * Eval::n_q_points, n_child_cell_dofs);
- if (Eval::dofs_per_cell != numbers::invalid_unsigned_int)
- AssertDimension(Eval::dofs_per_cell, n_points_in);
- const unsigned int n_points_out = n_child_cell_dofs / n_components;
- VectorizedArray<Number> *t0 = &evaluation_data[0];
- VectorizedArray<Number> *t1 = &evaluation_data[n_child_cell_dofs];
- VectorizedArray<Number> *t2 = &evaluation_data[2*n_child_cell_dofs];
-
- for (unsigned int c=0; c<n_components; ++c)
- {
- // for the prolongate case, we go from dofs (living on the parent cell) to
- // quads (living on all children) in the FEEvaluation terminology
- if (dim == 1)
- evaluator.template values<0,prolongate,false>(t0, t2);
- else if (dim == 2)
- {
- evaluator.template values<0,prolongate,false>(t0, t1);
- evaluator.template values<1,prolongate,false>(t1, t2);
- }
- else if (dim == 3)
- {
- evaluator.template values<0,prolongate,false>(t0, t2);
- evaluator.template values<1,prolongate,false>(t2, t1);
- evaluator.template values<2,prolongate,false>(t1, t2);
- }
- else
- Assert(false, ExcNotImplemented());
- if (prolongate)
- {
- t0 += n_points_in;
- t2 += n_points_out;
- }
- else
- {
- t0 += n_points_out;
- t2 += n_points_in;
- }
- }
- }
-
template <int dim, int degree, typename Number>
void weight_dofs_on_child (const VectorizedArray<Number> *weights,
const unsigned int n_components,
for (unsigned int cell=0; cell < n_owned_level_cells[to_level-1];
cell += vec_size)
{
- const unsigned int n_chunks = cell+vec_size > n_owned_level_cells[to_level-1] ?
- n_owned_level_cells[to_level-1] - cell : vec_size;
+ const unsigned int n_lanes = cell+vec_size > n_owned_level_cells[to_level-1] ?
+ n_owned_level_cells[to_level-1] - cell : vec_size;
// read from source vector
- for (unsigned int v=0; v<n_chunks; ++v)
+ for (unsigned int v=0; v<n_lanes; ++v)
{
const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
(parent_child_connect[to_level-1][cell+v].second,
// perform tensorized operation
if (element_is_continuous)
{
- typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,degree!=-1 ? 2*degree+1 : 0,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(prolongation_matrix_1d,
- prolongation_matrix_1d,
- prolongation_matrix_1d,
- fe_degree,
- 2*fe_degree+1);
- perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
- Utilities::fixed_power<dim>(fe_degree+1),
- n_child_cell_dofs,
- n_components,
- evaluation_data);
+ // must go through the components backwards because we want to write
+ // the output to the same array as the input
+ for (int c=n_components-1; c>=0; --c)
+ internal::FEEvaluationImplBasisChange<internal::evaluate_general,dim,degree+1,2*degree+1,
+ 1,VectorizedArray<Number>, VectorizedArray<Number> >
+ ::do_forward(prolongation_matrix_1d,
+ evaluation_data.begin()+ c*Utilities::fixed_power<dim>(degree_size),
+ evaluation_data.begin()+ c*n_scalar_cell_dofs,
+ fe_degree+1, 2*fe_degree+1);
weight_dofs_on_child<dim,degree,Number>(&weights_on_refined[to_level-1][(cell/vec_size)*three_to_dim],
n_components, fe_degree,
- &evaluation_data[2*n_child_cell_dofs]);
+ evaluation_data.begin());
}
else
{
- typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*(degree+1),VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(prolongation_matrix_1d,
- prolongation_matrix_1d,
- prolongation_matrix_1d,
- fe_degree,
- 2*(fe_degree+1));
- perform_tensorized_op<dim,Evaluator,Number,true>(evaluator,
- Utilities::fixed_power<dim>(fe_degree+1),
- n_child_cell_dofs,
- n_components,
- evaluation_data);
+ for (int c=n_components-1; c>=0; --c)
+ internal::FEEvaluationImplBasisChange<internal::evaluate_general,dim,degree+1,2*degree+2,
+ 1,VectorizedArray<Number>, VectorizedArray<Number> >
+ ::do_forward(prolongation_matrix_1d,
+ evaluation_data.begin() + c*Utilities::fixed_power<dim>(degree_size),
+ evaluation_data.begin() + c*n_scalar_cell_dofs,
+ fe_degree+1, 2*fe_degree+2);
}
// write into dst vector
const unsigned int *indices = &level_dof_indices[to_level][cell*
n_child_cell_dofs];
- for (unsigned int v=0; v<n_chunks; ++v)
+ for (unsigned int v=0; v<n_lanes; ++v)
{
for (unsigned int i=0; i<n_child_cell_dofs; ++i)
- dst.local_element(indices[i]) += evaluation_data[2*n_child_cell_dofs+i][v];
+ dst.local_element(indices[i]) += evaluation_data[i][v];
indices += n_child_cell_dofs;
}
}
for (unsigned int cell=0; cell < n_owned_level_cells[from_level-1];
cell += vec_size)
{
- const unsigned int n_chunks = cell+vec_size > n_owned_level_cells[from_level-1] ?
- n_owned_level_cells[from_level-1] - cell : vec_size;
+ const unsigned int n_lanes = cell+vec_size > n_owned_level_cells[from_level-1] ?
+ n_owned_level_cells[from_level-1] - cell : vec_size;
// read from source vector
{
const unsigned int *indices = &level_dof_indices[from_level][cell*
n_child_cell_dofs];
- for (unsigned int v=0; v<n_chunks; ++v)
+ for (unsigned int v=0; v<n_lanes; ++v)
{
for (unsigned int i=0; i<n_child_cell_dofs; ++i)
evaluation_data[i][v] = src.local_element(indices[i]);
// perform tensorized operation
if (element_is_continuous)
{
- typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,degree!=-1 ? 2*degree+1 : 0,VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(prolongation_matrix_1d,
- prolongation_matrix_1d,
- prolongation_matrix_1d,
- fe_degree,
- 2*fe_degree+1);
weight_dofs_on_child<dim,degree,Number>(&weights_on_refined[from_level-1][(cell/vec_size)*three_to_dim],
n_components, fe_degree,
&evaluation_data[0]);
- perform_tensorized_op<dim,Evaluator,Number,false>(evaluator,
- Utilities::fixed_power<dim>(fe_degree+1),
- n_child_cell_dofs,
- n_components,
- evaluation_data);
+ for (unsigned int c=0; c<n_components; ++c)
+ internal::FEEvaluationImplBasisChange<internal::evaluate_general,dim,degree+1,2*degree+1,
+ 1,VectorizedArray<Number>,VectorizedArray<Number> >
+ ::do_backward(prolongation_matrix_1d, false,
+ evaluation_data.begin() + c*n_scalar_cell_dofs,
+ evaluation_data.begin() + c*Utilities::fixed_power<dim>(degree_size),
+ fe_degree+1, 2*fe_degree+1);
}
else
{
- typedef internal::EvaluatorTensorProduct<internal::evaluate_general,dim,degree,2*(degree+1),VectorizedArray<Number> > Evaluator;
- Evaluator evaluator(prolongation_matrix_1d,
- prolongation_matrix_1d,
- prolongation_matrix_1d,
- fe_degree,
- 2*(fe_degree+1));
- perform_tensorized_op<dim,Evaluator,Number,false>(evaluator,
- Utilities::fixed_power<dim>(fe_degree+1),
- n_child_cell_dofs,
- n_components,
- evaluation_data);
+ for (unsigned int c=0; c<n_components; ++c)
+ internal::FEEvaluationImplBasisChange<internal::evaluate_general,dim,degree+1,2*degree+2,
+ 1,VectorizedArray<Number>,VectorizedArray<Number> >
+ ::do_backward(prolongation_matrix_1d, false,
+ evaluation_data.begin() + c*n_scalar_cell_dofs,
+ evaluation_data.begin() + c*Utilities::fixed_power<dim>(degree_size),
+ fe_degree+1, 2*fe_degree+2);
}
// write into dst vector
- for (unsigned int v=0; v<n_chunks; ++v)
+ for (unsigned int v=0; v<n_lanes; ++v)
{
const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
(parent_child_connect[from_level-1][cell+v].second,
{
// apply Dirichlet boundary conditions on parent cell
for (std::vector<unsigned short>::const_iterator i=dirichlet_indices[from_level-1][cell+v].begin(); i!=dirichlet_indices[from_level-1][cell+v].end(); ++i)
- evaluation_data[2*n_child_cell_dofs+(*i)][v] = 0.;
+ evaluation_data[*i][v] = 0.;
for (unsigned int k=0; k<(dim>2 ? degree_size : 1); ++k)
for (unsigned int j=0; j<(dim>1 ? degree_size : 1); ++j)
dst.local_element(indices[c*n_scalar_cell_dofs +
k*n_child_dofs_1d*n_child_dofs_1d+
j*n_child_dofs_1d+i])
- += evaluation_data[2*n_child_cell_dofs+m][v];
+ += evaluation_data[m][v];
}
}
}
-// this function tests the correctness of the 1d evaluation functions used in
-// FEEvaluation. These functions are marked 'internal' but it is much easier
-// to check their correctness directly rather than from the results in
-// dependent functions
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_symmetric
#include "../tests.h"
#include <iostream>
}
// apply function for tensor product
- internal::EvaluatorTensorProduct<internal::evaluate_symmetric,1,M-1,N,double> evaluator(shape, shape, shape);
+ internal::EvaluatorTensorProduct<internal::evaluate_symmetric,1,M,N,double> evaluator(shape, shape, shape);
if (type == 0)
evaluator.template values<0,false,add> (x,y);
if (type == 1)
return 0;
}
-
-// this function tests the correctness of the 1d evaluation functions used in
-// FEEvaluation. These functions are marked 'internal' but it is much easier
-// to check their correctness directly rather than from the results in
-// dependent functions. this function tests the even-odd path of the
-// evaluation functions
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_evenodd
#include "../tests.h"
#include <iostream>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/tensor_product_kernels.h>
template <int M, int N, int type, bool add>
}
// apply function for tensor product
- internal::EvaluatorTensorProduct<internal::evaluate_evenodd,1,M-1,N,double> evaluator(shape_sym, shape_sym, shape_sym);
+ internal::EvaluatorTensorProduct<internal::evaluate_evenodd,1,M,N,double> evaluator(shape_sym, shape_sym, shape_sym);
if (type == 0)
evaluator.template values<0,false,add> (x,y);
if (type == 1)
return 0;
}
-
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2013 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_symmetric_hierarchical
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ AlignedVector<double> shape(M*N);
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int j=0; j<(N+1)/2; ++j)
+ {
+ shape[i*N+j] = -1. + 2. * (double)Testing::rand()/RAND_MAX;
+ if (((i+type)%2)==1)
+ shape[i*N+N-1-j] = -shape[i*N+j];
+ else
+ shape[i*N+N-1-j] = shape[i*N+j];
+ if (j==N/2 && ((i+type)%2)==1)
+ shape[i*N+j] = 0.;
+ }
+
+ double x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = (double)Testing::rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i*N+j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_symmetric_hierarchical,1,M,N,double> evaluator(shape, shape, shape);
+ if (type == 0)
+ evaluator.template values<0,false,add> (x,y);
+ if (type == 1)
+ evaluator.template gradients<0,false,add> (x,y);
+ if (type == 2)
+ evaluator.template hessians<0,false,add> (x,y);
+
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ deallog << y[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ y[i] = (double)Testing::rand()/RAND_MAX;
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j*N+i] * y[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,add> (y,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,add> (y,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,add> (y,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ test<3,3,0,false>();
+ test<4,3,0,false>();
+ test<3,4,0,false>();
+ test<3,5,0,false>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,false>();
+ test<3,3,1,false>();
+ test<4,3,1,false>();
+ test<3,4,1,false>();
+ test<3,5,1,false>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,false>();
+ test<3,3,2,false>();
+ test<4,3,2,false>();
+ test<3,4,2,false>();
+ test<3,5,2,false>();
+ deallog.pop();
+
+ deallog.push("add");
+
+ deallog.push("values");
+ test<4,4,0,true>();
+ test<3,3,0,true>();
+ test<4,3,0,true>();
+ test<3,4,0,true>();
+ test<3,5,0,true>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,true>();
+ test<3,3,1,true>();
+ test<4,3,1,true>();
+ test<3,4,1,true>();
+ test<3,5,1,true>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,true>();
+ test<3,3,2,true>();
+ test<4,3,2,true>();
+ test<3,4,2,true>();
+ test<3,5,2,true>();
+ deallog.pop();
+
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0 0 0 0
+DEAL:values::Errors transpose: 0 0 0 0
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0 0 0
+DEAL:values::Errors transpose: 0 0 0
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0 0 0 0
+DEAL:values::Errors transpose: 0 0 0
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0 0 0
+DEAL:values::Errors transpose: 0 0 0 0
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0 0 0
+DEAL:values::Errors transpose: 0 0 0 0 0
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0 0 0 0
+DEAL:gradients::Errors transpose: 0 0 0 0
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0 0 0
+DEAL:gradients::Errors transpose: 0 0 0
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0 0 0 0
+DEAL:gradients::Errors transpose: 0 0 0
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0 0 0
+DEAL:gradients::Errors transpose: 0 0 0 0
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0 0 0
+DEAL:gradients::Errors transpose: 0 0 0 0 0
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0 0 0 0
+DEAL:hessians::Errors transpose: 0 0 0 0
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0 0 0
+DEAL:hessians::Errors transpose: 0 0 0
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0 0 0 0
+DEAL:hessians::Errors transpose: 0 0 0
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0 0 0
+DEAL:hessians::Errors transpose: 0 0 0 0
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0 0 0
+DEAL:hessians::Errors transpose: 0 0 0 0 0
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 0 0 0 0
+DEAL:add:values::Errors transpose: 0 0 0 0
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0 0 0
+DEAL:add:values::Errors transpose: 0 0 0
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 0 0 0 0
+DEAL:add:values::Errors transpose: 0 0 0
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 0 0 0
+DEAL:add:values::Errors transpose: 0 0 0 0
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 0 0 0
+DEAL:add:values::Errors transpose: 0 0 0 0 0
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0 0
+DEAL:add:gradients::Errors transpose: 0 0 0 0
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0
+DEAL:add:gradients::Errors transpose: 0 0 0
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0 0 0 0
+DEAL:add:gradients::Errors transpose: 0 0 0
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: 0 0 0
+DEAL:add:gradients::Errors transpose: 0 0 0 0
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 0 0 0
+DEAL:add:gradients::Errors transpose: 0 0 0 0 0
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0 0
+DEAL:add:hessians::Errors transpose: 0 0 0 0
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0
+DEAL:add:hessians::Errors transpose: 0 0 0
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0 0 0 0
+DEAL:add:hessians::Errors transpose: 0 0 0
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 0 0 0
+DEAL:add:hessians::Errors transpose: 0 0 0 0
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 0 0 0
+DEAL:add:hessians::Errors transpose: 0 0 0 0 0
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_general, when using same array for in and out
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+
+template <int M, int N, int type>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ AlignedVector<double> shape(M*N);
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ shape[i*N+j] = -1. + 2. * random_value<double>();
+
+ double x[N+M], x_ref[N], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y_ref[i] = 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i*N+j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_general,1,M,N,double>
+ evaluator(shape, shape, shape);
+ if (type == 0)
+ evaluator.template values<0,false,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,false,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,false,false> (x,x);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ deallog << x[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x_ref[i] = 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j*N+i] * x[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,false> (x,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0>();
+ test<3,3,0>();
+ test<4,3,0>();
+ test<3,4,0>();
+ test<3,5,0>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1>();
+ test<3,3,1>();
+ test<4,3,1>();
+ test<3,4,1>();
+ test<3,5,1>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2>();
+ test<3,3,2>();
+ test<4,3,2>();
+ test<3,4,2>();
+ test<3,5,2>();
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_evenodd, when using same array for in and out
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+
+template <int M, int N, int type>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ double shape[M][N];
+ for (unsigned int i=0; i<(M+1)/2; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ {
+ shape[i][j] = -1. + 2. * random_value<double>();
+ if (type == 1)
+ shape[M-1-i][N-1-j] = -shape[i][j];
+ else
+ shape[M-1-i][N-1-j] = shape[i][j];
+ }
+ if (type == 0 && M%2 == 1 && N%2 == 1)
+ {
+ for (unsigned int i=0; i<M; ++i)
+ shape[i][N/2] = 0.;
+ shape[M/2][N/2] = 1;
+ }
+ if (type == 1 && M%2 == 1 && N%2 == 1)
+ shape[M/2][N/2] = 0.;
+
+
+ // create symmetrized shape array exactly as expected by the evenodd
+ // function
+ AlignedVector<double> shape_sym(M*((N+1)/2));
+ for (unsigned int i=0; i<M/2; ++i)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ {
+ shape_sym[i*((N+1)/2)+q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
+ shape_sym[(M-1-i)*((N+1)/2)+q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
+ }
+ if (M % 2 == 1)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ shape_sym[(M-1)/2*((N+1)/2)+q] = shape[(M-1)/2][q];
+
+ double x[N+M], x_ref[N], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y_ref[i] = 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i][j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_evenodd,1,M,N,double> evaluator(shape_sym, shape_sym, shape_sym);
+ if (type == 0)
+ evaluator.template values<0,false,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,false,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,false,false> (x,x);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ deallog << x[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x_ref[i] = 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j][i] * x[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,false> (x,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0>();
+ test<3,3,0>();
+ test<4,3,0>();
+ test<3,4,0>();
+ test<3,5,0>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1>();
+ test<3,3,1>();
+ test<4,3,1>();
+ test<3,4,1>();
+ test<3,5,1>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2>();
+ test<3,3,2>();
+ test<4,3,2>();
+ test<3,4,2>();
+ test<3,5,2>();
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_symmetric_hierarchical, when using same array for in and out
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+
+
+template <int M, int N, int type>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ AlignedVector<double> shape(M*N);
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int j=0; j<(N+1)/2; ++j)
+ {
+ shape[i*N+j] = -1. + 2. * (double)Testing::rand()/RAND_MAX;
+ if (((i+type)%2)==1)
+ shape[i*N+N-1-j] = -shape[i*N+j];
+ else
+ shape[i*N+N-1-j] = shape[i*N+j];
+ if (j==N/2 && ((i+type)%2)==1)
+ shape[i*N+j] = 0.;
+ }
+
+ double x[N+M], x_ref[N], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y_ref[i] = 0.;
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i*N+j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_symmetric_hierarchical,1,M,N,double>
+ evaluator(shape, shape, shape);
+ if (type == 0)
+ evaluator.template values<0,false,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,false,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,false,false> (x,x);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ deallog << x[i] - y_ref[i] << " ";
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ x[i] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x_ref[i] = 0.;
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j*N+i] * x[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,false> (x,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,false> (x,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,false> (x,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ deallog << x[i] - x_ref[i] << " ";
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0>();
+ test<3,3,0>();
+ test<4,3,0>();
+ test<3,4,0>();
+ test<3,5,0>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1>();
+ test<3,3,1>();
+ test<4,3,1>();
+ test<3,4,1>();
+ test<3,5,1>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2>();
+ test<3,3,2>();
+ test<4,3,2>();
+ test<3,4,2>();
+ test<3,5,2>();
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_general, when using a double array for coefficients but
+// VectorizedArray for the input and output vector
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/base/vectorization.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ AlignedVector<double> shape(M*N);
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ shape[i*N+j] = -1. + 2. * random_value<double>();
+
+ VectorizedArray<double> x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ x[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i*N+j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_general,1,M,N,VectorizedArray<double>,double>
+ evaluator(shape, shape, shape);
+ if (type == 0)
+ evaluator.template values<0,false,add> (x,y);
+ if (type == 1)
+ evaluator.template gradients<0,false,add> (x,y);
+ if (type == 2)
+ evaluator.template hessians<0,false,add> (x,y);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ {
+ deallog << y[i][0] - y_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(y[i][v] - y_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ y[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j*N+i] * y[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,add> (y,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,add> (y,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,add> (y,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ {
+ deallog << x[i][0] - x_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(x[i][v] - x_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ test<3,3,0,false>();
+ test<4,3,0,false>();
+ test<3,4,0,false>();
+ test<3,5,0,false>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,false>();
+ test<3,3,1,false>();
+ test<4,3,1,false>();
+ test<3,4,1,false>();
+ test<3,5,1,false>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,false>();
+ test<3,3,2,false>();
+ test<4,3,2,false>();
+ test<3,4,2,false>();
+ test<3,5,2,false>();
+ deallog.pop();
+
+ deallog.push("add");
+
+ deallog.push("values");
+ test<4,4,0,true>();
+ test<3,3,0,true>();
+ test<4,3,0,true>();
+ test<3,4,0,true>();
+ test<3,5,0,true>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,true>();
+ test<3,3,1,true>();
+ test<4,3,1,true>();
+ test<3,4,1,true>();
+ test<3,5,1,true>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,true>();
+ test<3,3,2,true>();
+ test<4,3,2,true>();
+ test<3,4,2,true>();
+ test<3,5,2,true>();
+ deallog.pop();
+
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 1.11022e-16 2.77556e-17 5.55112e-17 1.11022e-16
+DEAL:values::Errors transpose: 2.22045e-16 0.00000 0.00000 2.22045e-16
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 5.55112e-17 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 -1.11022e-16
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 -1.11022e-16 1.38778e-17
+DEAL:values::Errors transpose: 0.00000 0.00000 -2.22045e-16 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 5.55112e-17 0.00000
+DEAL:values::Errors transpose: -5.55112e-17 0.00000 0.00000 0.00000 1.11022e-16
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 2.77556e-17 -2.77556e-17 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 -2.77556e-17 0.00000
+DEAL:gradients::Errors transpose: -5.55112e-17 1.38778e-17 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: -2.77556e-17 0.00000 1.11022e-16 0.00000
+DEAL:gradients::Errors transpose: 0.00000 -5.55112e-17 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 -7.63278e-17 0.00000
+DEAL:gradients::Errors transpose: -1.38778e-17 0.00000 1.11022e-16 -1.38778e-17
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: -1.11022e-16 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: -1.11022e-16 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 -5.55112e-17 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 5.55112e-17 0.00000 0.00000
+DEAL:hessians::Errors transpose: -1.11022e-16 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: -4.44089e-16 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 2.22045e-16
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 5.55112e-17 -5.55112e-17 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 -1.11022e-16 0.00000 -1.11022e-16
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 4.44089e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 4.44089e-16
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 -1.11022e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 -2.22045e-16
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 1.66533e-16 2.22045e-16 -1.11022e-16
+DEAL:add:values::Errors transpose: 0.00000 4.44089e-16 0.00000 0.00000
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 1.66533e-16 4.44089e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 0.00000 -5.55112e-17 0.00000
+DEAL:add:gradients::Errors transpose: 0.00000 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 1.11022e-16 -2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -2.22045e-16 0.00000
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 -2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: -2.22045e-16 1.11022e-16 2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -4.44089e-16 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 -2.22045e-16 -5.55112e-17
+DEAL:add:gradients::Errors transpose: -4.44089e-16 4.44089e-16 -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0.00000 -2.22045e-16 1.11022e-16 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 -2.22045e-16 1.11022e-16
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0.00000 0.00000 2.22045e-16 -1.11022e-16
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: -2.22045e-16 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 2.22045e-16
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_general, when using a double array for coefficients but
+// VectorizedArray for the input and output vector
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/base/vectorization.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ double shape[M][N];
+ for (unsigned int i=0; i<(M+1)/2; ++i)
+ for (unsigned int j=0; j<N; ++j)
+ {
+ shape[i][j] = -1. + 2. * random_value<double>();
+ if (type == 1)
+ shape[M-1-i][N-1-j] = -shape[i][j];
+ else
+ shape[M-1-i][N-1-j] = shape[i][j];
+ }
+ if (type == 0 && M%2 == 1 && N%2 == 1)
+ {
+ for (unsigned int i=0; i<M; ++i)
+ shape[i][N/2] = 0.;
+ shape[M/2][N/2] = 1;
+ }
+ if (type == 1 && M%2 == 1 && N%2 == 1)
+ shape[M/2][N/2] = 0.;
+
+ // create symmetrized shape array exactly as expected by the evenodd
+ // function
+ AlignedVector<double> shape_sym(M*((N+1)/2));
+ for (unsigned int i=0; i<M/2; ++i)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ {
+ shape_sym[i*((N+1)/2)+q] = 0.5 * (shape[i][q] + shape[i][N-1-q]);
+ shape_sym[(M-1-i)*((N+1)/2)+q] = 0.5 * (shape[i][q] - shape[i][N-1-q]);
+ }
+ if (M % 2 == 1)
+ for (unsigned int q=0; q<(N+1)/2; ++q)
+ shape_sym[(M-1)/2*((N+1)/2)+q] = shape[(M-1)/2][q];
+
+ VectorizedArray<double> x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ x[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i][j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_evenodd,1,M,N,VectorizedArray<double>,double>
+ evaluator(shape_sym, shape_sym, shape_sym);
+ if (type == 0)
+ evaluator.template values<0,false,add> (x,y);
+ if (type == 1)
+ evaluator.template gradients<0,false,add> (x,y);
+ if (type == 2)
+ evaluator.template hessians<0,false,add> (x,y);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ {
+ deallog << y[i][0] - y_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(y[i][v] - y_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ y[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j][i] * y[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,add> (y,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,add> (y,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,add> (y,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ {
+ deallog << x[i][0] - x_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(x[i][v] - x_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ test<3,3,0,false>();
+ test<4,3,0,false>();
+ test<3,4,0,false>();
+ test<3,5,0,false>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,false>();
+ test<3,3,1,false>();
+ test<4,3,1,false>();
+ test<3,4,1,false>();
+ test<3,5,1,false>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,false>();
+ test<3,3,2,false>();
+ test<4,3,2,false>();
+ test<3,4,2,false>();
+ test<3,5,2,false>();
+ deallog.pop();
+
+ deallog.push("add");
+
+ deallog.push("values");
+ test<4,4,0,true>();
+ test<3,3,0,true>();
+ test<4,3,0,true>();
+ test<3,4,0,true>();
+ test<3,5,0,true>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,true>();
+ test<3,3,1,true>();
+ test<4,3,1,true>();
+ test<3,4,1,true>();
+ test<3,5,1,true>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,true>();
+ test<3,3,2,true>();
+ test<4,3,2,true>();
+ test<3,4,2,true>();
+ test<3,5,2,true>();
+ deallog.pop();
+
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 1.11022e-16 2.77556e-17 5.55112e-17 1.11022e-16
+DEAL:values::Errors transpose: 2.22045e-16 0.00000 0.00000 2.22045e-16
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 5.55112e-17 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 -1.11022e-16
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 -1.11022e-16 1.38778e-17
+DEAL:values::Errors transpose: 0.00000 0.00000 -2.22045e-16 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 5.55112e-17 0.00000
+DEAL:values::Errors transpose: -5.55112e-17 0.00000 0.00000 0.00000 1.11022e-16
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 2.77556e-17 -2.77556e-17 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 -2.77556e-17 0.00000
+DEAL:gradients::Errors transpose: -5.55112e-17 1.38778e-17 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: -2.77556e-17 0.00000 1.11022e-16 0.00000
+DEAL:gradients::Errors transpose: 0.00000 -5.55112e-17 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 -7.63278e-17 0.00000
+DEAL:gradients::Errors transpose: -1.38778e-17 0.00000 1.11022e-16 -1.38778e-17
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: -1.11022e-16 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: -1.11022e-16 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 -5.55112e-17 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 5.55112e-17 0.00000 0.00000
+DEAL:hessians::Errors transpose: -1.11022e-16 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: -4.44089e-16 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 2.22045e-16
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 5.55112e-17 -5.55112e-17 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 -1.11022e-16 0.00000 -1.11022e-16
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 4.44089e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 4.44089e-16
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 -1.11022e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 -2.22045e-16
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 1.66533e-16 2.22045e-16 -1.11022e-16
+DEAL:add:values::Errors transpose: 0.00000 4.44089e-16 0.00000 0.00000
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 1.66533e-16 4.44089e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 0.00000 -5.55112e-17 0.00000
+DEAL:add:gradients::Errors transpose: 0.00000 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 1.11022e-16 -2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -2.22045e-16 0.00000
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 -2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: -2.22045e-16 1.11022e-16 2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -4.44089e-16 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 -2.22045e-16 -5.55112e-17
+DEAL:add:gradients::Errors transpose: -4.44089e-16 4.44089e-16 -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0.00000 -2.22045e-16 1.11022e-16 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 -2.22045e-16 1.11022e-16
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0.00000 0.00000 2.22045e-16 -1.11022e-16
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: -2.22045e-16 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 2.22045e-16
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the correctness of the 1d evaluation functions used in FEEvaluation,
+// path evaluate_symmetric_hierarchical, when using a double array for
+// coefficients but VectorizedArray for the input and output vector
+
+#include "../tests.h"
+#include <iostream>
+
+#include <deal.II/matrix_free/tensor_product_kernels.h>
+#include <deal.II/base/vectorization.h>
+
+
+template <int M, int N, int type, bool add>
+void test()
+{
+ deallog << "Test " << M << " x " << N << std::endl;
+ AlignedVector<double> shape(M*N);
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int j=0; j<(N+1)/2; ++j)
+ {
+ shape[i*N+j] = -1. + 2. * (double)Testing::rand()/RAND_MAX;
+ if (((i+type)%2)==1)
+ shape[i*N+N-1-j] = -shape[i*N+j];
+ else
+ shape[i*N+N-1-j] = shape[i*N+j];
+ if (j==N/2 && ((i+type)%2)==1)
+ shape[i*N+j] = 0.;
+ }
+
+ VectorizedArray<double> x[N], x_ref[N], y[M], y_ref[M];
+ for (unsigned int i=0; i<N; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ x[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<M; ++i)
+ {
+ y[i] = 1.;
+ y_ref[i] = add ? y[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<N; ++j)
+ y_ref[i] += shape[i*N+j] * x[j];
+ }
+
+ // apply function for tensor product
+ internal::EvaluatorTensorProduct<internal::evaluate_symmetric_hierarchical,1,M,N,VectorizedArray<double>,double>
+ evaluator(shape, shape, shape);
+ if (type == 0)
+ evaluator.template values<0,false,add> (x,y);
+ if (type == 1)
+ evaluator.template gradients<0,false,add> (x,y);
+ if (type == 2)
+ evaluator.template hessians<0,false,add> (x,y);
+
+ deallog << "Errors no transpose: ";
+ for (unsigned int i=0; i<M; ++i)
+ {
+ deallog << y[i][0] - y_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(y[i][v] - y_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+
+
+ for (unsigned int i=0; i<M; ++i)
+ for (unsigned int v=0; v<VectorizedArray<double>::n_array_elements; ++v)
+ y[i][v] = random_value<double>();
+
+ // compute reference
+ for (unsigned int i=0; i<N; ++i)
+ {
+ x[i] = 2.;
+ x_ref[i] = add ? x[i] : VectorizedArray<double>();
+ for (unsigned int j=0; j<M; ++j)
+ x_ref[i] += shape[j*N+i] * y[j];
+ }
+
+ // apply function for tensor product
+ if (type == 0)
+ evaluator.template values<0,true,add> (y,x);
+ if (type == 1)
+ evaluator.template gradients<0,true,add> (y,x);
+ if (type == 2)
+ evaluator.template hessians<0,true,add> (y,x);
+
+ deallog << "Errors transpose: ";
+ for (unsigned int i=0; i<N; ++i)
+ {
+ deallog << x[i][0] - x_ref[i][0] << " ";
+ for (unsigned int v=1; v<VectorizedArray<double>::n_array_elements; ++v)
+ AssertThrow(std::abs(x[i][v] - x_ref[i][v])<1e-12, ExcInternalError());
+ }
+ deallog << std::endl;
+}
+
+int main ()
+{
+ initlog();
+
+ deallog.push("values");
+ test<4,4,0,false>();
+ test<3,3,0,false>();
+ test<4,3,0,false>();
+ test<3,4,0,false>();
+ test<3,5,0,false>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,false>();
+ test<3,3,1,false>();
+ test<4,3,1,false>();
+ test<3,4,1,false>();
+ test<3,5,1,false>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,false>();
+ test<3,3,2,false>();
+ test<4,3,2,false>();
+ test<3,4,2,false>();
+ test<3,5,2,false>();
+ deallog.pop();
+
+ deallog.push("add");
+
+ deallog.push("values");
+ test<4,4,0,true>();
+ test<3,3,0,true>();
+ test<4,3,0,true>();
+ test<3,4,0,true>();
+ test<3,5,0,true>();
+ deallog.pop();
+
+ deallog.push("gradients");
+ test<4,4,1,true>();
+ test<3,3,1,true>();
+ test<4,3,1,true>();
+ test<3,4,1,true>();
+ test<3,5,1,true>();
+ deallog.pop();
+
+ deallog.push("hessians");
+ test<4,4,2,true>();
+ test<3,3,2,true>();
+ test<4,3,2,true>();
+ test<3,4,2,true>();
+ test<3,5,2,true>();
+ deallog.pop();
+
+ deallog.pop();
+
+ return 0;
+}
--- /dev/null
+
+DEAL:values::Test 4 x 4
+DEAL:values::Errors no transpose: 1.11022e-16 2.77556e-17 5.55112e-17 1.11022e-16
+DEAL:values::Errors transpose: 2.22045e-16 0.00000 0.00000 2.22045e-16
+DEAL:values::Test 3 x 3
+DEAL:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:values::Test 4 x 3
+DEAL:values::Errors no transpose: 5.55112e-17 0.00000 0.00000 0.00000
+DEAL:values::Errors transpose: 0.00000 0.00000 -1.11022e-16
+DEAL:values::Test 3 x 4
+DEAL:values::Errors no transpose: 0.00000 -1.11022e-16 1.38778e-17
+DEAL:values::Errors transpose: 0.00000 0.00000 -2.22045e-16 0.00000
+DEAL:values::Test 3 x 5
+DEAL:values::Errors no transpose: 0.00000 5.55112e-17 0.00000
+DEAL:values::Errors transpose: -5.55112e-17 0.00000 0.00000 0.00000 1.11022e-16
+DEAL:gradients::Test 4 x 4
+DEAL:gradients::Errors no transpose: 0.00000 0.00000 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 2.77556e-17 -2.77556e-17 0.00000
+DEAL:gradients::Test 3 x 3
+DEAL:gradients::Errors no transpose: 0.00000 -2.77556e-17 0.00000
+DEAL:gradients::Errors transpose: -5.55112e-17 1.38778e-17 0.00000
+DEAL:gradients::Test 4 x 3
+DEAL:gradients::Errors no transpose: -2.77556e-17 0.00000 1.11022e-16 0.00000
+DEAL:gradients::Errors transpose: 0.00000 -5.55112e-17 0.00000
+DEAL:gradients::Test 3 x 4
+DEAL:gradients::Errors no transpose: 0.00000 -7.63278e-17 0.00000
+DEAL:gradients::Errors transpose: -1.38778e-17 0.00000 1.11022e-16 -1.38778e-17
+DEAL:gradients::Test 3 x 5
+DEAL:gradients::Errors no transpose: -1.11022e-16 0.00000 0.00000
+DEAL:gradients::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:hessians::Test 4 x 4
+DEAL:hessians::Errors no transpose: -1.11022e-16 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 -5.55112e-17 0.00000
+DEAL:hessians::Test 3 x 3
+DEAL:hessians::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 -1.38778e-17 0.00000
+DEAL:hessians::Test 4 x 3
+DEAL:hessians::Errors no transpose: 0.00000 5.55112e-17 0.00000 0.00000
+DEAL:hessians::Errors transpose: -1.11022e-16 0.00000 0.00000
+DEAL:hessians::Test 3 x 4
+DEAL:hessians::Errors no transpose: -4.44089e-16 0.00000 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 0.00000 2.22045e-16
+DEAL:hessians::Test 3 x 5
+DEAL:hessians::Errors no transpose: 5.55112e-17 -5.55112e-17 0.00000
+DEAL:hessians::Errors transpose: 0.00000 0.00000 -1.11022e-16 0.00000 -1.11022e-16
+DEAL:add:values::Test 4 x 4
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 4.44089e-16 0.00000 0.00000 0.00000
+DEAL:add:values::Test 3 x 3
+DEAL:add:values::Errors no transpose: 0.00000 0.00000 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 4.44089e-16
+DEAL:add:values::Test 4 x 3
+DEAL:add:values::Errors no transpose: 2.22045e-16 0.00000 -1.11022e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 -2.22045e-16
+DEAL:add:values::Test 3 x 4
+DEAL:add:values::Errors no transpose: 1.66533e-16 2.22045e-16 -1.11022e-16
+DEAL:add:values::Errors transpose: 0.00000 4.44089e-16 0.00000 0.00000
+DEAL:add:values::Test 3 x 5
+DEAL:add:values::Errors no transpose: 1.66533e-16 4.44089e-16 0.00000
+DEAL:add:values::Errors transpose: 0.00000 0.00000 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 4 x 4
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 0.00000 -5.55112e-17 0.00000
+DEAL:add:gradients::Errors transpose: 0.00000 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Test 3 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 1.11022e-16 -2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -2.22045e-16 0.00000
+DEAL:add:gradients::Test 4 x 3
+DEAL:add:gradients::Errors no transpose: 0.00000 -2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:gradients::Errors transpose: 0.00000 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 4
+DEAL:add:gradients::Errors no transpose: -2.22045e-16 1.11022e-16 2.22045e-16
+DEAL:add:gradients::Errors transpose: 0.00000 -4.44089e-16 0.00000 0.00000
+DEAL:add:gradients::Test 3 x 5
+DEAL:add:gradients::Errors no transpose: 2.22045e-16 -2.22045e-16 -5.55112e-17
+DEAL:add:gradients::Errors transpose: -4.44089e-16 4.44089e-16 -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Test 4 x 4
+DEAL:add:hessians::Errors no transpose: 0.00000 -2.22045e-16 1.11022e-16 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 -2.22045e-16 1.11022e-16
+DEAL:add:hessians::Test 3 x 3
+DEAL:add:hessians::Errors no transpose: -2.22045e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 4 x 3
+DEAL:add:hessians::Errors no transpose: 0.00000 0.00000 2.22045e-16 -1.11022e-16
+DEAL:add:hessians::Errors transpose: 0.00000 2.22045e-16 0.00000
+DEAL:add:hessians::Test 3 x 4
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: -2.22045e-16 2.22045e-16 0.00000 -4.44089e-16
+DEAL:add:hessians::Test 3 x 5
+DEAL:add:hessians::Errors no transpose: 1.11022e-16 0.00000 0.00000
+DEAL:add:hessians::Errors transpose: 0.00000 0.00000 0.00000 0.00000 2.22045e-16
projects / dealii.git / commitdiff