From 13fc76dc4a97ed584032e86493d74df85bbbcedb Mon Sep 17 00:00:00 2001 From: Katharina Kormann <katharina.kormann@tum.de> Date: Mon, 5 Feb 2018 21:53:43 +0100 Subject: [PATCH] Integrators for DG face integrals with tests. --- include/deal.II/lac/tensor_product_matrix.h | 12 +- .../deal.II/matrix_free/evaluation_kernels.h | 1352 ++++++++++---- .../deal.II/matrix_free/evaluation_selector.h | 254 ++- include/deal.II/matrix_free/fe_evaluation.h | 12 +- include/deal.II/matrix_free/operators.h | 2 +- .../matrix_free/tensor_product_kernels.h | 1659 ++++++++++++----- source/fe/mapping_q_generic.cc | 30 +- .../matrix_free/evaluation_selector.inst.in | 14 +- source/multigrid/mg_transfer_matrix_free.cc | 151 +- ...te_1d_shape.cc => evaluate_1d_shape_01.cc} | 9 +- ...ape.output => evaluate_1d_shape_01.output} | 0 ...ape_evenodd.cc => evaluate_1d_shape_02.cc} | 12 +- ...odd.output => evaluate_1d_shape_02.output} | 0 tests/matrix_free/evaluate_1d_shape_03.cc | 156 ++ tests/matrix_free/evaluate_1d_shape_03.output | 91 + tests/matrix_free/evaluate_1d_shape_04.cc | 118 ++ tests/matrix_free/evaluate_1d_shape_04.output | 46 + tests/matrix_free/evaluate_1d_shape_05.cc | 145 ++ tests/matrix_free/evaluate_1d_shape_05.output | 46 + tests/matrix_free/evaluate_1d_shape_06.cc | 126 ++ tests/matrix_free/evaluate_1d_shape_06.output | 46 + tests/matrix_free/evaluate_1d_shape_07.cc | 160 ++ tests/matrix_free/evaluate_1d_shape_07.output | 91 + tests/matrix_free/evaluate_1d_shape_08.cc | 186 ++ tests/matrix_free/evaluate_1d_shape_08.output | 91 + tests/matrix_free/evaluate_1d_shape_09.cc | 168 ++ tests/matrix_free/evaluate_1d_shape_09.output | 91 + 27 files changed, 3912 insertions(+), 1156 deletions(-) rename tests/matrix_free/{evaluate_1d_shape.cc => evaluate_1d_shape_01.cc} (92%) rename tests/matrix_free/{evaluate_1d_shape.output => evaluate_1d_shape_01.output} (100%) rename tests/matrix_free/{evaluate_1d_shape_evenodd.cc => evaluate_1d_shape_02.cc} (90%) rename tests/matrix_free/{evaluate_1d_shape_evenodd.output => evaluate_1d_shape_02.output} (100%) create mode 100644 tests/matrix_free/evaluate_1d_shape_03.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_03.output create mode 100644 tests/matrix_free/evaluate_1d_shape_04.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_04.output create mode 100644 tests/matrix_free/evaluate_1d_shape_05.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_05.output create mode 100644 tests/matrix_free/evaluate_1d_shape_06.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_06.output create mode 100644 tests/matrix_free/evaluate_1d_shape_07.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_07.output create mode 100644 tests/matrix_free/evaluate_1d_shape_08.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_08.output create mode 100644 tests/matrix_free/evaluate_1d_shape_09.cc create mode 100644 tests/matrix_free/evaluate_1d_shape_09.output diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 5642f9439c..80c6481ce0 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -454,10 +454,10 @@ TensorProductMatrixSymmetricSumBase<dim,Number,size> 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]); @@ -515,10 +515,10 @@ TensorProductMatrixSymmetricSumBase<dim,Number,size> 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]); diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index f5bc0cd9e1..7f69067eb7 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -100,38 +100,40 @@ namespace internal 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) @@ -141,22 +143,22 @@ namespace internal 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; @@ -169,15 +171,13 @@ namespace internal 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) @@ -186,39 +186,37 @@ namespace internal { 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; @@ -228,33 +226,39 @@ namespace internal // 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; @@ -264,9 +268,9 @@ namespace internal 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) @@ -274,27 +278,27 @@ namespace internal // 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) @@ -302,26 +306,32 @@ namespace internal // 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; @@ -332,9 +342,14 @@ namespace internal // 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]; + } } @@ -344,32 +359,33 @@ namespace internal 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; @@ -382,23 +398,14 @@ namespace internal 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) { @@ -406,79 +413,93 @@ namespace internal 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; @@ -489,20 +510,28 @@ namespace internal // 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) @@ -512,7 +541,7 @@ namespace internal 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; } @@ -525,337 +554,846 @@ namespace internal /** - * 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 diff --git a/include/deal.II/matrix_free/evaluation_selector.h b/include/deal.II/matrix_free/evaluation_selector.h index 8757c4119a..db65bbcbd6 100644 --- a/include/deal.II/matrix_free/evaluation_selector.h +++ b/include/deal.II/matrix_free/evaluation_selector.h @@ -45,12 +45,12 @@ namespace 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) @@ -62,11 +62,11 @@ namespace 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) { @@ -74,7 +74,7 @@ namespace 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); } }; @@ -108,12 +108,12 @@ namespace 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) @@ -131,11 +131,11 @@ namespace } 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) { @@ -158,12 +158,12 @@ namespace 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) @@ -171,19 +171,17 @@ namespace 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, @@ -196,34 +194,32 @@ namespace 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> @@ -239,12 +235,12 @@ namespace * 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) @@ -264,11 +260,11 @@ namespace * 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) { @@ -303,12 +299,12 @@ struct SelectEvaluator * 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); @@ -320,11 +316,11 @@ struct SelectEvaluator * 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); }; @@ -350,12 +346,12 @@ struct SelectEvaluator<dim, -1, n_q_points_1d, n_components, Number> * 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); @@ -368,11 +364,11 @@ struct SelectEvaluator<dim, -1, n_q_points_1d, n_components, Number> * 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); }; @@ -384,12 +380,12 @@ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename 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) @@ -404,10 +400,10 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::evaluate 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); @@ -454,11 +450,11 @@ template <int dim, int fe_degree, int n_q_points_1d, int n_components, typename 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) { @@ -470,15 +466,15 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate 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) { @@ -486,7 +482,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate 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) { @@ -494,7 +490,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate 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) { @@ -502,7 +498,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate 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) { @@ -510,7 +506,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate 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()); @@ -522,12 +518,12 @@ template <int dim, int dummy, int n_components, typename Number> 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) @@ -567,11 +563,11 @@ template <int dim, int dummy, int n_components, typename Number> 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) { @@ -581,7 +577,7 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate 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) { @@ -589,14 +585,14 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate 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, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 1ca1a606ee..5a286358c4 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -5393,9 +5393,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number> 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 @@ -5427,9 +5427,9 @@ FEEvaluation<dim,fe_degree,n_q_points_1d,n_components_,Number> 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 diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index bb1e5f1bc3..88e0bbd0fb 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -892,7 +892,7 @@ namespace MatrixFreeOperators 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); diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 8fe79420e1..91b3f21bc0 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// 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. // @@ -50,26 +50,71 @@ namespace internal * 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 @@ -85,94 +130,200 @@ namespace internal /** * 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; @@ -180,77 +331,94 @@ namespace internal 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 @@ -259,8 +427,8 @@ namespace internal 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; @@ -270,10 +438,10 @@ namespace internal // 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: @@ -286,13 +454,17 @@ namespace internal } 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; + } } } } @@ -300,15 +472,23 @@ namespace internal /** - * 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 @@ -319,102 +499,157 @@ namespace internal 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; @@ -422,16 +657,115 @@ namespace internal 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; @@ -440,10 +774,19 @@ namespace internal 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; + } } } } @@ -456,48 +799,70 @@ namespace internal * 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; }; @@ -520,23 +885,24 @@ namespace internal // 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) { @@ -544,16 +910,17 @@ namespace internal { 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) { @@ -565,15 +932,15 @@ namespace internal 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)]; @@ -585,24 +952,24 @@ namespace internal } 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) @@ -616,27 +983,24 @@ namespace internal 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 @@ -646,20 +1010,18 @@ namespace internal 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]; @@ -672,29 +1034,11 @@ namespace internal 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); } } @@ -719,23 +1063,24 @@ namespace internal // 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) { @@ -743,16 +1088,17 @@ namespace internal { 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) { @@ -764,15 +1110,15 @@ namespace internal 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)]; @@ -786,13 +1132,13 @@ namespace internal 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) { @@ -807,22 +1153,21 @@ namespace internal } 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; @@ -830,31 +1175,11 @@ namespace internal 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); } } @@ -863,23 +1188,24 @@ namespace internal // 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) { @@ -887,16 +1213,17 @@ namespace internal { 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) { @@ -908,15 +1235,15 @@ namespace internal 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)]; @@ -930,13 +1257,13 @@ namespace internal 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) { @@ -951,34 +1278,33 @@ namespace internal } 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) @@ -987,29 +1313,11 @@ namespace internal 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); } } @@ -1032,99 +1340,202 @@ namespace internal * 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 @@ -1137,7 +1548,7 @@ namespace internal 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)]; @@ -1148,51 +1559,52 @@ namespace internal 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; @@ -1200,20 +1612,20 @@ namespace internal 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) @@ -1225,14 +1637,14 @@ namespace internal 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) @@ -1254,38 +1666,375 @@ namespace internal 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); } } } diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 9c2647fd0d..3d0fb42da5 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -1451,11 +1451,6 @@ namespace internal 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); @@ -1478,6 +1473,9 @@ namespace internal 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) @@ -1489,26 +1487,10 @@ namespace internal = 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); } diff --git a/source/matrix_free/evaluation_selector.inst.in b/source/matrix_free/evaluation_selector.inst.in index 6313bb8980..968a2a22ab 100644 --- a/source/matrix_free/evaluation_selector.inst.in +++ b/source/matrix_free/evaluation_selector.inst.in @@ -18,18 +18,18 @@ for (deal_II_dimension : DIMENSIONS; components : SPACE_DIMENSIONS; scalar_type { 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); } diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc index 57de8ca64b..2952db3dec 100644 --- a/source/multigrid/mg_transfer_matrix_free.cc +++ b/source/multigrid/mg_transfer_matrix_free.cc @@ -27,7 +27,7 @@ #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> @@ -138,7 +138,7 @@ void MGTransferMatrixFree<dim,Number>::build } } - evaluation_data.resize(3*n_child_cell_dofs); + evaluation_data.resize(n_child_cell_dofs); } @@ -271,55 +271,6 @@ void MGTransferMatrixFree<dim,Number> 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, @@ -368,11 +319,11 @@ void MGTransferMatrixFree<dim,Number> 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, @@ -399,43 +350,37 @@ void MGTransferMatrixFree<dim,Number> // 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; } } @@ -459,14 +404,14 @@ void MGTransferMatrixFree<dim,Number> 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]); @@ -479,38 +424,30 @@ void MGTransferMatrixFree<dim,Number> // 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, @@ -523,7 +460,7 @@ void MGTransferMatrixFree<dim,Number> { // 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) @@ -531,7 +468,7 @@ void MGTransferMatrixFree<dim,Number> 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]; } } } diff --git a/tests/matrix_free/evaluate_1d_shape.cc b/tests/matrix_free/evaluate_1d_shape_01.cc similarity index 92% rename from tests/matrix_free/evaluate_1d_shape.cc rename to tests/matrix_free/evaluate_1d_shape_01.cc index 4304c1a4c0..0a27f9082e 100644 --- a/tests/matrix_free/evaluate_1d_shape.cc +++ b/tests/matrix_free/evaluate_1d_shape_01.cc @@ -15,10 +15,8 @@ -// 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> @@ -63,7 +61,7 @@ void test() } // 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) @@ -161,4 +159,3 @@ int main () return 0; } - diff --git a/tests/matrix_free/evaluate_1d_shape.output b/tests/matrix_free/evaluate_1d_shape_01.output similarity index 100% rename from tests/matrix_free/evaluate_1d_shape.output rename to tests/matrix_free/evaluate_1d_shape_01.output diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.cc b/tests/matrix_free/evaluate_1d_shape_02.cc similarity index 90% rename from tests/matrix_free/evaluate_1d_shape_evenodd.cc rename to tests/matrix_free/evaluate_1d_shape_02.cc index 80363a2bcf..d0c5e11eee 100644 --- a/tests/matrix_free/evaluate_1d_shape_evenodd.cc +++ b/tests/matrix_free/evaluate_1d_shape_02.cc @@ -15,16 +15,13 @@ -// 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> @@ -77,7 +74,7 @@ void test() } // 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) @@ -176,4 +173,3 @@ int main () return 0; } - diff --git a/tests/matrix_free/evaluate_1d_shape_evenodd.output b/tests/matrix_free/evaluate_1d_shape_02.output similarity index 100% rename from tests/matrix_free/evaluate_1d_shape_evenodd.output rename to tests/matrix_free/evaluate_1d_shape_02.output diff --git a/tests/matrix_free/evaluate_1d_shape_03.cc b/tests/matrix_free/evaluate_1d_shape_03.cc new file mode 100644 index 0000000000..556f4186e9 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_03.cc @@ -0,0 +1,156 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_03.output b/tests/matrix_free/evaluate_1d_shape_03.output new file mode 100644 index 0000000000..2bd9530d6c --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_03.output @@ -0,0 +1,91 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_04.cc b/tests/matrix_free/evaluate_1d_shape_04.cc new file mode 100644 index 0000000000..e8d90956c2 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_04.cc @@ -0,0 +1,118 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_04.output b/tests/matrix_free/evaluate_1d_shape_04.output new file mode 100644 index 0000000000..ad134a113d --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_04.output @@ -0,0 +1,46 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_05.cc b/tests/matrix_free/evaluate_1d_shape_05.cc new file mode 100644 index 0000000000..43c715b8aa --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_05.cc @@ -0,0 +1,145 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_05.output b/tests/matrix_free/evaluate_1d_shape_05.output new file mode 100644 index 0000000000..ad134a113d --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_05.output @@ -0,0 +1,46 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_06.cc b/tests/matrix_free/evaluate_1d_shape_06.cc new file mode 100644 index 0000000000..003e8ae299 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_06.cc @@ -0,0 +1,126 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_06.output b/tests/matrix_free/evaluate_1d_shape_06.output new file mode 100644 index 0000000000..ad134a113d --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_06.output @@ -0,0 +1,46 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_07.cc b/tests/matrix_free/evaluate_1d_shape_07.cc new file mode 100644 index 0000000000..463e9bc2ae --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_07.cc @@ -0,0 +1,160 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_07.output b/tests/matrix_free/evaluate_1d_shape_07.output new file mode 100644 index 0000000000..bce88a2f41 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_07.output @@ -0,0 +1,91 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_08.cc b/tests/matrix_free/evaluate_1d_shape_08.cc new file mode 100644 index 0000000000..4871a6768d --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_08.cc @@ -0,0 +1,186 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_08.output b/tests/matrix_free/evaluate_1d_shape_08.output new file mode 100644 index 0000000000..bce88a2f41 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_08.output @@ -0,0 +1,91 @@ + +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 diff --git a/tests/matrix_free/evaluate_1d_shape_09.cc b/tests/matrix_free/evaluate_1d_shape_09.cc new file mode 100644 index 0000000000..62734a91ea --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_09.cc @@ -0,0 +1,168 @@ +// --------------------------------------------------------------------- +// +// 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; +} diff --git a/tests/matrix_free/evaluate_1d_shape_09.output b/tests/matrix_free/evaluate_1d_shape_09.output new file mode 100644 index 0000000000..bce88a2f41 --- /dev/null +++ b/tests/matrix_free/evaluate_1d_shape_09.output @@ -0,0 +1,91 @@ + +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 -- 2.39.5