From: Katharina Kormann Date: Mon, 5 Feb 2018 20:53:43 +0000 (+0100) Subject: Integrators for DG face integrals with tests. X-Git-Tag: v9.0.0-rc1~406^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5870%2Fhead;p=dealii.git Integrators for DG face integrals with tests. --- 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 Threads::Mutex::ScopedLock lock(this->mutex); const unsigned int n = Utilities::fixed_power(size > 0 ? size : eigenvalues[0].size()); tmp_array.resize_fast(n*2); - constexpr int kernel_size = size > 0 ? size-1 : -1; - internal::EvaluatorTensorProduct + constexpr int kernel_size = size > 0 ? size : 0; + internal::EvaluatorTensorProduct eval(AlignedVector {}, AlignedVector {}, - AlignedVector {}, mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows()); + AlignedVector {}, 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 Threads::Mutex::ScopedLock lock(this->mutex); const unsigned int n = size > 0 ? size : eigenvalues[0].size(); tmp_array.resize_fast (Utilities::fixed_power(n)); - constexpr int kernel_size = size > 0 ? size-1 : -1; - internal::EvaluatorTensorProduct + constexpr int kernel_size = size > 0 ? size : 0; + internal::EvaluatorTensorProduct eval(AlignedVector(), AlignedVector(), - AlignedVector(), mass_matrix[0].n_rows()-1, mass_matrix[0].n_rows()); + AlignedVector(), 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> &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + void evaluate (const MatrixFreeFunctions::ShapeInfo &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> &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool evaluate_values, - const bool evaluate_gradients); + void integrate (const MatrixFreeFunctions::ShapeInfo &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 inline void FEEvaluationImpl - ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + ::evaluate (const MatrixFreeFunctions::ShapeInfo &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 = EvaluatorSelector4)>::variant; - typedef EvaluatorTensorProduct > Eval; + typedef EvaluatorTensorProduct 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 *temp1; - VectorizedArray *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 **values_dofs = values_dofs_actual; - VectorizedArray *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(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(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(); + values_dofs_tmp[c*dofs_per_comp+count_q] = Number(); } for (int j=degree+1-i; j(); + values_dofs_tmp[c*dofs_per_comp+count_q] = Number(); } - AssertDimension(count_q, Utilities::fixed_power(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 (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 - ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + ::integrate (const MatrixFreeFunctions::ShapeInfo &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 = EvaluatorSelector4)>::variant; - typedef EvaluatorTensorProduct > Eval; + typedef EvaluatorTensorProduct 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 *temp1; - VectorizedArray *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(shape_info.fe_degree+1) : shape_info.dofs_per_component_on_cell; // expand dof_values to tensor product for truncated tensor products - VectorizedArray **values_dofs = values_dofs_actual; - VectorizedArray *expanded_dof_values[n_components]; - if (type == MatrixFreeFunctions::truncated_tensor) - { - values_dofs = expanded_dof_values; - for (unsigned int c=0; c(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 (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 (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 (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(); + { + for (unsigned int c=0; c2?degree+1:1); ++i) @@ -512,7 +541,7 @@ namespace internal for (int k=0; k + 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 &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 eval_val (transformation_matrix, + AlignedVector(), + AlignedVector(), + 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(np_1); + values_out = values_out + n_components*Utilities::fixed_power(np_2); + for (unsigned int c=n_components; c!=0; --c) + { + values_in -= Utilities::fixed_power(np_1); + values_out -= Utilities::fixed_power(np_2); + if (next_dim < dim) + for (unsigned int q=np_1; q!=0; --q) + FEEvaluationImplBasisChange + ::do_forward(transformation_matrix, + values_in + (q-1)*Utilities::fixed_power(np_1), + values_out + (q-1)*Utilities::fixed_power(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(values_in, values_out); + else + eval_val.template values(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 &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 eval_val (transformation_matrix, + AlignedVector(), + AlignedVector(), + 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 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(values_in, values_in); + } + if (next_dim < dim) + for (unsigned int q=0; q + ::do_backward(transformation_matrix, + add_into_result, + values_in + q*Utilities::fixed_power(np_2), + values_out + q*Utilities::fixed_power(np_1), + basis_size_1_variable, basis_size_2_variable); + + values_in += Utilities::fixed_power(np_2); + values_out += Utilities::fixed_power(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 &transformation_matrix, + const AlignedVector &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 + ::do_forward(transformation_matrix, + values_in + (q-1)*Utilities::fixed_int_power::value, + my_scratch + (q-1)*Utilities::fixed_int_power::value); + EvaluatorTensorProduct 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::value; + for (unsigned int ii=0; ii (my_scratch+i, my_scratch+i); + for (unsigned int q=0; q(my_scratch+i, my_scratch+i); + } + for (unsigned int q=0; q + ::do_backward(transformation_matrix, false, + my_scratch + q*Utilities::fixed_int_power::value, + values_out + q*Utilities::fixed_int_power::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 - struct FEEvaluationImplTransformToCollocation + struct FEEvaluationImplCollocation { static - void evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians); + void evaluate (const MatrixFreeFunctions::ShapeInfo &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> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool integrate_values, - const bool integrate_gradients); + void integrate (const MatrixFreeFunctions::ShapeInfo &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 inline void - FEEvaluationImplTransformToCollocation - ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *, - const bool, - const bool evaluate_gradients, - const bool evaluate_hessians) + FEEvaluationImplCollocation + ::evaluate (const MatrixFreeFunctions::ShapeInfo &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 > Eval; - Eval eval_val (shape_info.shape_values_eo, - AlignedVector >(), - AlignedVector >(), - shape_info.fe_degree, - shape_info.n_q_points_1d); - Eval eval(AlignedVector >(), - 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 + eval(AlignedVector(), + shape_info.shape_gradients_collocation_eo, + shape_info.shape_hessians_collocation_eo); + constexpr unsigned int n_q_points = Utilities::fixed_int_power::value; for (unsigned int c=0; c(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(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 inline void - FEEvaluationImplTransformToCollocation - ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *, - const bool integrate_values, - const bool integrate_gradients) + FEEvaluationImplCollocation + ::integrate (const MatrixFreeFunctions::ShapeInfo &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 > Eval; - Eval eval_val (shape_info.shape_values_eo, - AlignedVector >(), - AlignedVector >(), - shape_info.fe_degree, - shape_info.n_q_points_1d); - Eval eval(AlignedVector >(), - 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 + eval(AlignedVector(), + shape_info.shape_gradients_collocation_eo, + shape_info.shape_hessians_collocation_eo); + constexpr unsigned int n_q_points = Utilities::fixed_int_power::value; for (unsigned int c=0; c(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 - struct FEEvaluationImplCollocation + * @author Katharina Kormann, Martin Kronbichler, 2017 + */ + template + struct FEEvaluationImplTransformToCollocation { static - void evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians); + void evaluate (const MatrixFreeFunctions::ShapeInfo &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> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, - const bool integrate_values, - const bool integrate_gradients); + void integrate (const MatrixFreeFunctions::ShapeInfo &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 + + + template inline void - FEEvaluationImplCollocation - ::evaluate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *, - const bool evaluate_values, - const bool evaluate_gradients, - const bool evaluate_hessians) + FEEvaluationImplTransformToCollocation + ::evaluate (const MatrixFreeFunctions::ShapeInfo &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 > Eval; - Eval eval(AlignedVector >(), - 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::value; for (unsigned int c=0; c=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:: + 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 + + + template inline void - FEEvaluationImplCollocation - ::integrate (const MatrixFreeFunctions::ShapeInfo> &shape_info, - VectorizedArray *values_dofs[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *, - const bool integrate_values, - const bool integrate_gradients) + FEEvaluationImplTransformToCollocation + ::integrate (const MatrixFreeFunctions::ShapeInfo &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 > Eval; - Eval eval(AlignedVector >(), - 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::value; for (unsigned int c=0; c(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:: + 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=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 + struct FEFaceEvaluationImpl + { + static + void evaluate_in_face (const MatrixFreeFunctions::ShapeInfo &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 &val1 + = symmetric_evaluate ? data.shape_values_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_values : data.values_within_subface[subface_index%2]); + const AlignedVector &val2 + = symmetric_evaluate ? data.shape_values_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_values : data.values_within_subface[subface_index/2]); + + const AlignedVector &grad1 + = symmetric_evaluate ? data.shape_gradients_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_gradients : data.gradients_within_subface[subface_index%2]); + const AlignedVector &grad2 + = symmetric_evaluate ? data.shape_gradients_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_gradients : data.gradients_within_subface[subface_index/2]); + + typedef internal::EvaluatorTensorProduct + Eval; + typedef internal::EvaluatorTensorProduct + EvalGeneric; + Eval eval1(val1,grad1,AlignedVector(), + data.fe_degree+1, data.n_q_points_1d); + Eval eval2(val2,grad2,AlignedVector(), + data.fe_degree+1, data.n_q_points_1d); + + const unsigned int size_deg = fe_degree > -1 ? + Utilities::fixed_int_power::value : + (dim > 1 ? Utilities::fixed_power(data.fe_degree+1) : 1); + + const unsigned int n_q_points = fe_degree > -1 ? + Utilities::fixed_int_power::value : data.n_q_points_face; + + if (evaluate_grad == false) + for (unsigned int c=0; c(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 fe_degree) + { + eval1.template values<0,true,false>(values_dofs, values_quad); + eval1.template values<1,true,false>(values_quad, values_quad); + internal::EvaluatorTensorProduct + eval_grad + (AlignedVector(), + data.shape_gradients_collocation_eo, + AlignedVector()); + 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 &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 &val1 + = symmetric_evaluate ? data.shape_values_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_values : data.values_within_subface[subface_index%2]); + const AlignedVector &val2 + = symmetric_evaluate ? data.shape_values_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_values : data.values_within_subface[subface_index/2]); + + const AlignedVector &grad1 + = symmetric_evaluate ? data.shape_gradients_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_gradients : data.gradients_within_subface[subface_index%2]); + const AlignedVector &grad2 + = symmetric_evaluate ? data.shape_gradients_eo : + (subface_index >= GeometryInfo::max_children_per_cell ? + data.shape_gradients : data.gradients_within_subface[subface_index/2]); + + typedef internal::EvaluatorTensorProduct + Eval; + typedef internal::EvaluatorTensorProduct + 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::value : + (dim > 1 ? Utilities::fixed_power(data.fe_degree+1) : 1); + + const unsigned int n_q_points = fe_degree > -1 ? + Utilities::fixed_int_power::value : data.n_q_points_face; + + if (integrate_grad == false) + for (unsigned int c=0; c(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 (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 eval_grad + (AlignedVector(), + data.shape_gradients_collocation_eo, + AlignedVector()); + 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 + struct FEFaceNormalEvaluationImpl + { + template + static void interpolate(const MatrixFreeFunctions::ShapeInfo &data, + const Number *input, + Number *output, + const bool do_gradients, + const unsigned int face_no) + { + internal::EvaluatorTensorProduct + evalf(data.shape_data_on_face[face_no%2], + AlignedVector(), + AlignedVector(), + 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(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 struct Default { - static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo &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 > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo &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 { static inline void evaluate ( - const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + const internal::MatrixFreeFunctions::ShapeInfo &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 > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + const internal::MatrixFreeFunctions::ShapeInfo &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 struct Factory::type> { - static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + static inline void evaluate (const internal::MatrixFreeFunctions::ShapeInfo &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 - ::evaluate(shape_info, values_dofs_actual, values_quad, - gradients_quad, hessians_quad, scratch_data, - evaluate_values, evaluate_gradients, evaluate_hessians); - else - internal::FEEvaluationImplTransformToCollocation - ::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 + ::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 + ::evaluate(shape_info, values_dofs_actual, values_quad, + gradients_quad, hessians_quad, scratch_data, + evaluate_values, evaluate_gradients, evaluate_hessians); else internal::FEEvaluationImpl ::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 > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + static inline void integrate (const internal::MatrixFreeFunctions::ShapeInfo &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 - ::integrate(shape_info, values_dofs_actual, values_quad, - gradients_quad, scratch_data, - integrate_values, integrate_gradients); - else - internal::FEEvaluationImplTransformToCollocation - ::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 + ::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 + ::integrate(shape_info, values_dofs_actual, values_quad, + gradients_quad, scratch_data, + integrate_values, integrate_gradients, false); else internal::FEEvaluationImpl ::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 @@ -239,12 +235,12 @@ namespace * for the 'evaluate' function. */ template - void symmetric_selector_evaluate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + void symmetric_selector_evaluate (const internal::MatrixFreeFunctions::ShapeInfo &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 - void symmetric_selector_integrate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + void symmetric_selector_integrate (const internal::MatrixFreeFunctions::ShapeInfo &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 > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo &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 > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + static void integrate(const internal::MatrixFreeFunctions::ShapeInfo &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 * internal::FEEvaluationImplTransformToCollocation::evaluate() with appropriate * template parameters. */ - static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, + static void evaluate(const internal::MatrixFreeFunctions::ShapeInfo &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 * internal::FEEvaluationImplTransformToCollocation::integrate() with appropriate * template parameters. */ - static void integrate(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, + static void integrate(const internal::MatrixFreeFunctions::ShapeInfo &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 ::evaluate -(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, +(const internal::MatrixFreeFunctions::ShapeInfo &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::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 + internal::FEEvaluationImplTransformToCollocation ::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 ::integrate -(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, +(const internal::MatrixFreeFunctions::ShapeInfo &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::integrate internal::FEEvaluationImplCollocation ::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 + internal::FEEvaluationImplTransformToCollocation ::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::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::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::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::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 inline void SelectEvaluator::evaluate -(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *hessians_quad[][(dim*(dim+1))/2], - VectorizedArray *scratch_data, +(const internal::MatrixFreeFunctions::ShapeInfo &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 inline void SelectEvaluator::integrate -(const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *values_dofs_actual[], - VectorizedArray *values_quad[], - VectorizedArray *gradients_quad[][dim], - VectorizedArray *scratch_data, +(const internal::MatrixFreeFunctions::ShapeInfo &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::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::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 ::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 (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 Assert(this->matrix_info != nullptr || this->mapped_geometry->is_initialized(), ExcNotInitialized()); - SelectEvaluator - ::evaluate (*this->data, &this->values_dofs[0], this->values_quad, - this->gradients_quad, this->hessians_quad, this->scratch_data, + SelectEvaluator > + ::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 Assert(this->matrix_info != nullptr || this->mapped_geometry->is_initialized(), ExcNotInitialized()); - SelectEvaluator - ::integrate (*this->data, &this->values_dofs[0], this->values_quad, - this->gradients_quad, this->scratch_data, + SelectEvaluator > + ::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 > 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 + template 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 - struct EvaluatorTensorProduct + template + struct EvaluatorTensorProduct { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::value; + static const unsigned int n_rows_of_product = Utilities::fixed_int_power::value; + static const unsigned int n_columns_of_product = Utilities::fixed_int_power::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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + template void values (const Number in [], Number out[]) const { - apply(shape_values, in, out); + apply(shape_values, in, out); } - template + template void gradients (const Number in [], Number out[]) const { - apply(shape_gradients, in, out); + apply(shape_gradients, in, out); } - template + template void hessians (const Number in [], Number out[]) const { - apply(shape_hessians, in, out); + apply(shape_hessians, in, out); + } + + template + void + values_one_line (const Number in [], + Number out[]) const + { + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } + + template + void + gradients_one_line (const Number in [], + Number out[]) const + { + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } + + template + void + hessians_one_line (const Number in [], + Number out[]) const + { + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, in, out); } - template - 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 + 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 + 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 - template + + + template + template inline void - EvaluatorTensorProduct - ::apply (const Number *shape_data, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::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::value; + constexpr int stride = Utilities::fixed_int_power::value; + constexpr int n_blocks1 = one_line ? 1 : stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=dim)?0:(dim-direction-1)>::value; for (int i2=0; i2 + template + template inline void - apply_tensor_product_face (const Number *shape_data, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::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::value; + AssertIndexRange (face_direction, dim); + constexpr int stride = Utilities::fixed_int_power::value; + constexpr int out_stride = Utilities::fixed_int_power::value; + const Number *DEAL_II_RESTRICT shape_values = this->shape_values; for (int i2=0; i2 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 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 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 - struct EvaluatorTensorProduct + template + struct EvaluatorTensorProduct { - 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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int fe_degree, - const unsigned int n_q_points_1d) + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + template void values (const Number *in, Number *out) const { - apply(shape_values, in, out); + apply(shape_values, in, out); } - template + template void gradients (const Number *in, Number *out) const { - apply(shape_gradients, in, out); + apply(shape_gradients, in, out); } - template + template void hessians (const Number *in, Number *out) const { - apply(shape_hessians, in, out); + apply(shape_hessians, in, out); + } + + template + void + values_one_line (const Number in [], + Number out[]) const + { + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } + + template + void + gradients_one_line (const Number in [], + Number out[]) const + { + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } + + template + void + hessians_one_line (const Number in [], + Number out[]) const + { + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, in, out); } - template - void apply (const Number *shape_data, - const Number *in, - Number *out) const; + template + void apply (const Number2 *DEAL_II_RESTRICT shape_data, + const Number *in, + Number *out) const; + + template + 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 - template + + + template + template inline void - EvaluatorTensorProduct - ::apply (const Number *shape_data, - const Number *in, - Number *out) const + EvaluatorTensorProduct + ::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(nn); + const int stride = direction==0 ? 1 : Utilities::fixed_power(n_columns); + const int n_blocks1 = one_line ? 1 : stride; + const int n_blocks2 = direction >= dim-1 ? 1 : Utilities::fixed_power(n_rows); + Assert(n_rows <= 128, ExcNotImplemented()); for (int i2=0; i2 + template + inline + void + EvaluatorTensorProduct + ::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(n_rows) : 1; + const int out_stride = dim > 1 ? Utilities::fixed_power(n_rows) : 1; + + for (int i2=0; i2 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 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 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 - struct EvaluatorTensorProduct + template + struct EvaluatorTensorProduct { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::value; + static const unsigned int n_rows_of_product = Utilities::fixed_int_power::value; + static const unsigned int n_columns_of_product = Utilities::fixed_int_power::value; /** * Constructor, taking the data from ShapeInfo */ - EvaluatorTensorProduct (const AlignedVector &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) + EvaluatorTensorProduct (const AlignedVector &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + template void values (const Number in [], Number out[]) const; - template + template void gradients (const Number in [], Number out[]) const; - template + template 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 - template + template + template inline void - EvaluatorTensorProduct + EvaluatorTensorProduct ::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::value; + constexpr int stride = Utilities::fixed_int_power::value; + constexpr int n_blocks1 = stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=dim)?0:(dim-direction-1)>::value; for (int i2=0; i2 0) { @@ -565,15 +932,15 @@ namespace internal res1 += val0 * in1; for (int ind=1; ind 0) { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; + res0 = val0 * (in[0] + in[stride*(mm-1)]); for (int ind=1; ind 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 - template + template + template inline void - EvaluatorTensorProduct + EvaluatorTensorProduct ::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::value; + constexpr int stride = Utilities::fixed_int_power::value; + constexpr int n_blocks1 = stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=dim)?0:(dim-direction-1)>::value; for (int i2=0; i2 0) { @@ -764,15 +1110,15 @@ namespace internal res1 -= val0 * in1; for (int ind=1; ind - template + template + template inline void - EvaluatorTensorProduct + EvaluatorTensorProduct ::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::value; + constexpr int stride = Utilities::fixed_int_power::value; + constexpr int n_blocks1 = stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=dim)?0:(dim-direction-1)>::value; for (int i2=0; i2 0) { @@ -908,15 +1235,15 @@ namespace internal res1 += val0 * in1; for (int ind=1; ind 0) { - res0 = in[0] + in[stride*(mm-1)]; - res0 *= val0; + res0 = val0 * (in[0] + in[stride*(mm-1)]); for (int ind=1; ind - struct EvaluatorTensorProduct + template + struct EvaluatorTensorProduct { - static const unsigned int dofs_per_cell = Utilities::fixed_int_power::value; - static const unsigned int n_q_points = Utilities::fixed_int_power::value; + static const unsigned int n_rows_of_product = Utilities::fixed_int_power::value; + static const unsigned int n_columns_of_product = Utilities::fixed_int_power::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 &shape_values, - const AlignedVector &shape_gradients, - const AlignedVector &shape_hessians, - const unsigned int dummy1 = 0, - const unsigned int dummy2 = 0) + EvaluatorTensorProduct (const AlignedVector &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 &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + template void values (const Number in [], Number out[]) const { - apply(shape_values, in, out); + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); } - template + template void gradients (const Number in [], Number out[]) const { - apply(shape_gradients, in, out); + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); } - template + template void hessians (const Number in [], Number out[]) const { - apply(shape_hessians, in, out); + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, in, out); } - template - static void apply (const Number *shape_data, - const Number in [], - Number out []); + template + void + values_one_line (const Number in [], + Number out[]) const + { + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } + + template + void + gradients_one_line (const Number in [], + Number out[]) const + { + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } - const Number *shape_values; - const Number *shape_gradients; - const Number *shape_hessians; + template + void + hessians_one_line (const Number in [], + Number out[]) const + { + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(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 + 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 - template + template + template inline void - EvaluatorTensorProduct - ::apply (const Number *shapes, - const Number in [], - Number out []) + EvaluatorTensorProduct + ::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::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::value; + constexpr int n_blocks1 = one_line ? 1 : stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=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 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 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 + struct EvaluatorTensorProduct + { + static const unsigned int n_rows_of_product = Utilities::fixed_int_power::value; + static const unsigned int n_columns_of_product = Utilities::fixed_int_power::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 &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 &shape_values, + const AlignedVector &shape_gradients, + const AlignedVector &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 + void + values (const Number in [], + Number out[]) const + { + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } + + template + void + gradients (const Number in [], + Number out[]) const + { + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } + + template + void + hessians (const Number in [], + Number out[]) const + { + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(shape_hessians, in, out); + } + + template + void + values_one_line (const Number in [], + Number out[]) const + { + Assert (shape_values != nullptr, ExcNotInitialized()); + apply(shape_values, in, out); + } + + template + void + gradients_one_line (const Number in [], + Number out[]) const + { + Assert (shape_gradients != nullptr, ExcNotInitialized()); + apply(shape_gradients, in, out); + } + + template + void + hessians_one_line (const Number in [], + Number out[]) const + { + Assert (shape_hessians != nullptr, ExcNotInitialized()); + apply(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 + 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 + template + inline + void + EvaluatorTensorProduct + ::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::value; + constexpr int n_blocks1 = one_line ? 1 : stride; + constexpr int n_blocks2 = Utilities::fixed_int_power=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 0) + { + r0 = shapes[col] * x[0]; + r1 = shapes[col+n_columns] * x[1]; + for (unsigned int ind=1; ind0) + { + r0 = shapes[n_cols + shift*n_columns] * x[shift]; + for (unsigned int ind=1; ind0?mid:1]; + for (int i=0; i 0) + { + r0 = shapes[2*col*n_columns] * xp[0]; + r1 = shapes[(2*col+1)*n_columns] * xm[0]; + for (unsigned int ind=1; ind 0) + { + r0 = shapes[(nn-1)*n_columns] * xp[0]; + for (unsigned int ind=1; ind *values_dofs_ptr[n_comp]; - VectorizedArray *values_quad_ptr[n_comp]; - VectorizedArray *gradients_quad_ptr[n_comp][dim]; - VectorizedArray *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 &renumber_to_lexicographic = data.shape_info.lexicographic_numbering; for (unsigned int i=0; i::evaluate - (data.shape_info, &(values_dofs_ptr[0]), &(values_quad_ptr[0]), - &(gradients_quad_ptr[0]), &(hessians_quad_ptr[0]), &(data.scratch[0]), + SelectEvaluator >::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::integrate + SelectEvaluator >::integrate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *[], VectorizedArray *[], - VectorizedArray *[][deal_II_dimension], VectorizedArray *, + VectorizedArray *, VectorizedArray *, + VectorizedArray *, VectorizedArray *, const bool, const bool); template void - SelectEvaluator::evaluate + SelectEvaluator >::evaluate (const internal::MatrixFreeFunctions::ShapeInfo > &shape_info, - VectorizedArray *[], VectorizedArray *[], - VectorizedArray *[][deal_II_dimension], - VectorizedArray *[][(deal_II_dimension*(deal_II_dimension+1))/2], + VectorizedArray *, VectorizedArray *, + VectorizedArray *, + VectorizedArray *, VectorizedArray *, 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 #include -#include +#include #include @@ -138,7 +138,7 @@ void MGTransferMatrixFree::build } } - evaluation_data.resize(3*n_child_cell_dofs); + evaluation_data.resize(n_child_cell_dofs); } @@ -271,55 +271,6 @@ void MGTransferMatrixFree namespace { - template - 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 > &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 *t0 = &evaluation_data[0]; - VectorizedArray *t1 = &evaluation_data[n_child_cell_dofs]; - VectorizedArray *t2 = &evaluation_data[2*n_child_cell_dofs]; - - for (unsigned int c=0; c(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 void weight_dofs_on_child (const VectorizedArray *weights, const unsigned int n_components, @@ -368,11 +319,11 @@ void MGTransferMatrixFree 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 (parent_child_connect[to_level-1][cell+v].second, @@ -399,43 +350,37 @@ void MGTransferMatrixFree // perform tensorized operation if (element_is_continuous) { - typedef internal::EvaluatorTensorProduct > Evaluator; - Evaluator evaluator(prolongation_matrix_1d, - prolongation_matrix_1d, - prolongation_matrix_1d, - fe_degree, - 2*fe_degree+1); - perform_tensorized_op(evaluator, - Utilities::fixed_power(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, VectorizedArray > + ::do_forward(prolongation_matrix_1d, + evaluation_data.begin()+ c*Utilities::fixed_power(degree_size), + evaluation_data.begin()+ c*n_scalar_cell_dofs, + fe_degree+1, 2*fe_degree+1); weight_dofs_on_child(&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 > Evaluator; - Evaluator evaluator(prolongation_matrix_1d, - prolongation_matrix_1d, - prolongation_matrix_1d, - fe_degree, - 2*(fe_degree+1)); - perform_tensorized_op(evaluator, - Utilities::fixed_power(fe_degree+1), - n_child_cell_dofs, - n_components, - evaluation_data); + for (int c=n_components-1; c>=0; --c) + internal::FEEvaluationImplBasisChange, VectorizedArray > + ::do_forward(prolongation_matrix_1d, + evaluation_data.begin() + c*Utilities::fixed_power(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 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 // perform tensorized operation if (element_is_continuous) { - typedef internal::EvaluatorTensorProduct > Evaluator; - Evaluator evaluator(prolongation_matrix_1d, - prolongation_matrix_1d, - prolongation_matrix_1d, - fe_degree, - 2*fe_degree+1); weight_dofs_on_child(&weights_on_refined[from_level-1][(cell/vec_size)*three_to_dim], n_components, fe_degree, &evaluation_data[0]); - perform_tensorized_op(evaluator, - Utilities::fixed_power(fe_degree+1), - n_child_cell_dofs, - n_components, - evaluation_data); + for (unsigned int c=0; c,VectorizedArray > + ::do_backward(prolongation_matrix_1d, false, + evaluation_data.begin() + c*n_scalar_cell_dofs, + evaluation_data.begin() + c*Utilities::fixed_power(degree_size), + fe_degree+1, 2*fe_degree+1); } else { - typedef internal::EvaluatorTensorProduct > Evaluator; - Evaluator evaluator(prolongation_matrix_1d, - prolongation_matrix_1d, - prolongation_matrix_1d, - fe_degree, - 2*(fe_degree+1)); - perform_tensorized_op(evaluator, - Utilities::fixed_power(fe_degree+1), - n_child_cell_dofs, - n_components, - evaluation_data); + for (unsigned int c=0; c,VectorizedArray > + ::do_backward(prolongation_matrix_1d, false, + evaluation_data.begin() + c*n_scalar_cell_dofs, + evaluation_data.begin() + c*Utilities::fixed_power(degree_size), + fe_degree+1, 2*fe_degree+2); } // write into dst vector - for (unsigned int v=0; v (parent_child_connect[from_level-1][cell+v].second, @@ -523,7 +460,7 @@ void MGTransferMatrixFree { // apply Dirichlet boundary conditions on parent cell for (std::vector::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 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 @@ -63,7 +61,7 @@ void test() } // apply function for tensor product - internal::EvaluatorTensorProduct evaluator(shape, shape, shape); + internal::EvaluatorTensorProduct 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 -#include +#include template @@ -77,7 +74,7 @@ void test() } // apply function for tensor product - internal::EvaluatorTensorProduct evaluator(shape_sym, shape_sym, shape_sym); + internal::EvaluatorTensorProduct 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 + +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + AlignedVector shape(M*N); + for (unsigned int i=0; i 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 (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(); + 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 + +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + AlignedVector shape(M*N); + for (unsigned int i=0; i(); + + double x[N+M], x_ref[N], y_ref[M]; + for (unsigned int i=0; i(); + + // compute reference + for (unsigned int i=0; i + 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(); + + // compute reference + for (unsigned int i=0; i (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(); + 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 + +#include + + +template +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(); + 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 shape_sym(M*((N+1)/2)); + for (unsigned int i=0; i(); + + // compute reference + for (unsigned int i=0; i 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(); + + // compute reference + for (unsigned int i=0; i (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(); + 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 + +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + AlignedVector shape(M*N); + for (unsigned int i=0; i(); + + // compute reference + for (unsigned int i=0; i + 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(); + + // compute reference + for (unsigned int i=0; i (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(); + 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 + +#include +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + AlignedVector shape(M*N); + for (unsigned int i=0; i(); + + VectorizedArray x[N], x_ref[N], y[M], y_ref[M]; + for (unsigned int i=0; i::n_array_elements; ++v) + x[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j,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::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::n_array_elements; ++v) + y[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j (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_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 + +#include +#include + + +template +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(); + 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 shape_sym(M*((N+1)/2)); + for (unsigned int i=0; i x[N], x_ref[N], y[M], y_ref[M]; + for (unsigned int i=0; i::n_array_elements; ++v) + x[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j,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::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::n_array_elements; ++v) + y[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j (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_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 + +#include +#include + + +template +void test() +{ + deallog << "Test " << M << " x " << N << std::endl; + AlignedVector shape(M*N); + for (unsigned int i=0; i x[N], x_ref[N], y[M], y_ref[M]; + for (unsigned int i=0; i::n_array_elements; ++v) + x[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j,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::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::n_array_elements; ++v) + y[i][v] = random_value(); + + // compute reference + for (unsigned int i=0; i(); + for (unsigned int j=0; j (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_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