From 6918f728e9368e9ff05c6968d60613bdea815074 Mon Sep 17 00:00:00 2001 From: Niklas Wik Date: Tue, 1 Mar 2022 18:30:34 +0100 Subject: [PATCH] Matrix-free evaluation of FE_RaviartThomasNodal Test added for matrix-free RT Split of evaluation_template_factory.templates.h Decrease FE_EVAL_FACTORY_DEGREE_MAX to 2 for windows vm --- .github/workflows/windows.yml | 2 +- include/deal.II/fe/fe_raviart_thomas.h | 8 + .../deal.II/matrix_free/evaluation_kernels.h | 1332 ++++++++++++++++- ...aluation_template_face_factory.templates.h | 125 ++ .../evaluation_template_factory.templates.h | 91 -- .../matrix_free/matrix_free.templates.h | 15 +- include/deal.II/matrix_free/shape_info.h | 13 +- .../matrix_free/shape_info.templates.h | 176 ++- .../matrix_free/tensor_product_kernels.h | 497 +++++- source/fe/fe_raviart_thomas_nodal.cc | 123 +- source/matrix_free/CMakeLists.txt | 11 + .../evaluation_template_face_factory.cc | 28 + .../evaluation_template_face_factory.inst.in | 32 + .../evaluation_template_face_factory_inst2.cc | 17 + .../evaluation_template_face_factory_inst3.cc | 17 + .../evaluation_template_face_factory_inst4.cc | 17 + .../evaluation_template_face_factory_inst5.cc | 17 + .../evaluation_template_face_factory_inst6.cc | 17 + .../evaluation_template_factory.cc | 2 +- .../evaluation_template_factory.inst.in | 13 - .../evaluation_template_factory_inst10.cc | 17 + .../evaluation_template_factory_inst7.cc | 17 + .../evaluation_template_factory_inst8.cc | 17 + .../evaluation_template_factory_inst9.cc | 17 + tests/matrix_free/matrix_vector_rt_01.cc | 44 + tests/matrix_free/matrix_vector_rt_01.output | 25 + tests/matrix_free/matrix_vector_rt_common.h | 241 +++ tests/matrix_free/matrix_vector_rt_face_01.cc | 44 + .../matrix_vector_rt_face_01.output | 19 + .../matrix_vector_rt_face_common.h | 311 ++++ 30 files changed, 3092 insertions(+), 213 deletions(-) create mode 100644 include/deal.II/matrix_free/evaluation_template_face_factory.templates.h create mode 100644 source/matrix_free/evaluation_template_face_factory.cc create mode 100644 source/matrix_free/evaluation_template_face_factory.inst.in create mode 100644 source/matrix_free/evaluation_template_face_factory_inst2.cc create mode 100644 source/matrix_free/evaluation_template_face_factory_inst3.cc create mode 100644 source/matrix_free/evaluation_template_face_factory_inst4.cc create mode 100644 source/matrix_free/evaluation_template_face_factory_inst5.cc create mode 100644 source/matrix_free/evaluation_template_face_factory_inst6.cc create mode 100644 source/matrix_free/evaluation_template_factory_inst10.cc create mode 100644 source/matrix_free/evaluation_template_factory_inst7.cc create mode 100644 source/matrix_free/evaluation_template_factory_inst8.cc create mode 100644 source/matrix_free/evaluation_template_factory_inst9.cc create mode 100644 tests/matrix_free/matrix_vector_rt_01.cc create mode 100644 tests/matrix_free/matrix_vector_rt_01.output create mode 100644 tests/matrix_free/matrix_vector_rt_common.h create mode 100644 tests/matrix_free/matrix_vector_rt_face_01.cc create mode 100644 tests/matrix_free/matrix_vector_rt_face_01.output create mode 100644 tests/matrix_free/matrix_vector_rt_face_common.h diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml index b46066f105..cd90aafcaf 100644 --- a/.github/workflows/windows.yml +++ b/.github/workflows/windows.yml @@ -29,7 +29,7 @@ jobs: mkdir build mkdir c:/project cd build - cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX" -T host=x64 -A x64 .. + cmake -DCMAKE_BUILD_TYPE=Debug -DCMAKE_INSTALL_PREFIX=c:/project -DDEAL_II_WITH_ZLIB=off -DDEAL_II_CXX_FLAGS="-WX /D FE_EVAL_FACTORY_DEGREE_MAX=2" -T host=x64 -A x64 .. - name: archive logs uses: actions/upload-artifact@v1 with: diff --git a/include/deal.II/fe/fe_raviart_thomas.h b/include/deal.II/fe/fe_raviart_thomas.h index ab12990c5e..fcd9e93e4d 100644 --- a/include/deal.II/fe/fe_raviart_thomas.h +++ b/include/deal.II/fe/fe_raviart_thomas.h @@ -327,6 +327,14 @@ private: * FiniteElementData::degree is higher by one than the constructor * argument! */ + +namespace internal +{ + template + std::vector + get_lexicographic_numbering_rt_nodal(const unsigned int degree); +} // namespace internal + template class FE_RaviartThomasNodal : public FE_PolyTensor { diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index e2dc5a5147..1faeb7a335 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -86,6 +86,12 @@ namespace internal static const EvaluatorVariant variant = evaluate_evenodd; }; + template + struct EvaluatorSelector + { + static const EvaluatorVariant variant = evaluate_raviart_thomas; + }; + /** @@ -180,6 +186,79 @@ namespace internal const bool add_into_values_array); }; + /** + * Specialization for MatrixFreeFunctions::tensor_raviart_thomas, which use + * specific sum-factorization kernels and with normal/tangential shape_data + */ + template + struct FEEvaluationImpl + { + using EvalNormal = + EvaluatorTensorProduct; + + using EvalTangent = + EvaluatorTensorProduct; + template + static void + evaluate_or_integrate( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array = false); + + private: + template + static EvalType + create_evaluator_tensor_product( + const MatrixFreeFunctions::ShapeInfo &shape_info) + { + // Shape_data is of size 2 with normal direction in [0] and tangential + // direction in [1]. See shape_info.template.h + const MatrixFreeFunctions::UnivariateShapeData *shape_data; + if (std::is_same::value) + { + shape_data = &shape_info.data[0]; + } + else + { + shape_data = &shape_info.data[1]; + } + return EvalType(shape_data->shape_values, + shape_data->shape_gradients, + shape_data->shape_hessians); + } + + template + static void + evaluate_tensor_product_per_component( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array, + std::integral_constant); + + template + static void + evaluate_tensor_product_per_component( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array, + std::integral_constant); + }; + template + template + inline void + FEEvaluationImpl:: + evaluate_or_integrate( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array) + { + if (evaluation_flag == EvaluationFlags::nothing) + return; + + AssertDimension(fe_eval.get_shape_info().data.size(), 2); + // First component: + evaluate_tensor_product_per_component<0>( + evaluation_flag, + values_dofs_actual, + fe_eval, + add_into_values_array, + std::integral_constant()); + // Second component : + evaluate_tensor_product_per_component<1>( + evaluation_flag, + values_dofs_actual, + fe_eval, + add_into_values_array, + std::integral_constant()); + if (dim == 3) + { + // Third component + evaluate_tensor_product_per_component<2>( + evaluation_flag, + values_dofs_actual, + fe_eval, + add_into_values_array, + std::integral_constant()); + } + } + + // Helper function that applies the 1d evaluation kernals. + // std::integral_constant is the interpolation path, and + // std::integral_constant bellow is the integration path. + template + template + inline void + FEEvaluationImpl:: + evaluate_tensor_product_per_component( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array, + std::integral_constant) + { + (void)add_into_values_array; + using Eval0 = + typename std::conditional::type; + using Eval1 = + typename std::conditional::type; + using Eval2 = + typename std::conditional::type; + + Eval0 eval0 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + Eval1 eval1 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + Eval2 eval2 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + + Number *temp1 = fe_eval.get_scratch_data().begin(); + Number *temp2; + + temp2 = temp1 + std::max(Utilities::fixed_power( + fe_eval.get_shape_info().data[0].fe_degree + 1), + Utilities::fixed_power( + fe_eval.get_shape_info().data[0].n_q_points_1d)); + + const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points; + const std::size_t dofs_per_comp = + fe_eval.get_shape_info().dofs_per_component_on_cell; + + // Initial shift depending on component (normal_dir) + Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; + Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir; + Number *gradients_quad = + fe_eval.begin_gradients() + dim * n_q_points * normal_dir; + Number *hessians_quad = + (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir : + fe_eval.begin_hessians() + 6 * n_q_points * normal_dir; + + switch (dim) + { + case 2: + if (evaluation_flag & EvaluationFlags::gradients) + { + eval0.template gradients<0, true, false, normal_dir>(values_dofs, + temp1); + eval1.template values<1, true, false, normal_dir>(temp1, + gradients_quad); + } + if (evaluation_flag & EvaluationFlags::hessians) + { + // The evaluation/integration here *should* work, however + // the piola transform is not implemented. + AssertThrow(false, ExcNotImplemented()); + // grad xy + if (!(evaluation_flag & EvaluationFlags::gradients)) + eval0.template gradients<0, true, false, normal_dir>( + values_dofs, temp1); + eval1.template gradients<1, true, false, normal_dir>( + temp1, hessians_quad + 2 * n_q_points); + + // grad xx + eval0.template hessians<0, true, false, normal_dir>(values_dofs, + temp1); + eval1.template values<1, true, false, normal_dir>(temp1, + hessians_quad); + } + + // grad y + eval0.template values<0, true, false, normal_dir>(values_dofs, temp1); + if (evaluation_flag & EvaluationFlags::gradients) + eval1.template gradients<1, true, false, normal_dir>( + temp1, gradients_quad + n_q_points); + + // grad yy + if (evaluation_flag & EvaluationFlags::hessians) + eval1.template hessians<1, true, false, normal_dir>(temp1, + hessians_quad + + n_q_points); + + // val: can use values applied in x + if (evaluation_flag & EvaluationFlags::values) + eval1.template values<1, true, false, normal_dir>(temp1, + values_quad); + break; + case 3: + if (evaluation_flag & EvaluationFlags::gradients) + { + // grad x + eval0.template gradients<0, true, false, normal_dir>(values_dofs, + temp1); + eval1.template values<1, true, false, normal_dir>(temp1, temp2); + eval2.template values<2, true, false, normal_dir>(temp2, + gradients_quad); + } + + if (evaluation_flag & EvaluationFlags::hessians) + { + // The evaluation/integration here *should* work, however + // the piola transform is not implemented. + AssertThrow(false, ExcNotImplemented()); + // grad xz + if (!(evaluation_flag & EvaluationFlags::gradients)) + { + eval0.template gradients<0, true, false, normal_dir>( + values_dofs, temp1); + eval1.template values<1, true, false, normal_dir>(temp1, + temp2); + } + eval2.template gradients<2, true, false, normal_dir>( + temp2, hessians_quad + 4 * n_q_points); + + // grad xy + eval1.template gradients<1, true, false, normal_dir>(temp1, + temp2); + eval2.template values<2, true, false, normal_dir>( + temp2, hessians_quad + 3 * n_q_points); + + // grad xx + eval0.template hessians<0, true, false, normal_dir>(values_dofs, + temp1); + eval1.template values<1, true, false, normal_dir>(temp1, temp2); + eval2.template values<2, true, false, normal_dir>(temp2, + hessians_quad); + } + + // grad y + eval0.template values<0, true, false, normal_dir>(values_dofs, temp1); + if (evaluation_flag & EvaluationFlags::gradients) + { + eval1.template gradients<1, true, false, normal_dir>(temp1, + temp2); + eval2.template values<2, true, false, normal_dir>(temp2, + gradients_quad + + n_q_points); + } + + if (evaluation_flag & EvaluationFlags::hessians) + { + // grad yz + if (!(evaluation_flag & EvaluationFlags::gradients)) + eval1.template gradients<1, true, false, normal_dir>(temp1, + temp2); + eval2.template gradients<2, true, false, normal_dir>( + temp2, hessians_quad + 5 * n_q_points); + + // grad yy + eval1.template hessians<1, true, false, normal_dir>(temp1, temp2); + eval2.template values<2, true, false, normal_dir>(temp2, + hessians_quad + + n_q_points); + } + + // grad z: can use the values applied in x direction stored in + // temp1 + eval1.template values<1, true, false, normal_dir>(temp1, temp2); + if (evaluation_flag & EvaluationFlags::gradients) + eval2.template gradients<2, true, false, normal_dir>( + temp2, gradients_quad + 2 * n_q_points); + + // grad zz: can use the values applied in x and y direction stored + // in temp2 + if (evaluation_flag & EvaluationFlags::hessians) + eval2.template hessians<2, true, false, normal_dir>( + temp2, hessians_quad + 2 * n_q_points); + + // val: can use the values applied in x & y direction stored in + // temp2 + if (evaluation_flag & EvaluationFlags::values) + eval2.template values<2, true, false, normal_dir>(temp2, + values_quad); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + + template + template + inline void + FEEvaluationImpl:: + evaluate_tensor_product_per_component( + const EvaluationFlags::EvaluationFlags evaluation_flag, + Number * values_dofs_actual, + FEEvaluationData & fe_eval, + const bool add_into_values_array, + std::integral_constant) + { + using Eval0 = + typename std::conditional::type; + using Eval1 = + typename std::conditional::type; + using Eval2 = + typename std::conditional::type; + + Eval0 eval0 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + Eval1 eval1 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + Eval2 eval2 = + create_evaluator_tensor_product(fe_eval.get_shape_info()); + + Number *temp1 = fe_eval.get_scratch_data().begin(); + Number *temp2; + + temp2 = temp1 + std::max(Utilities::fixed_power( + fe_eval.get_shape_info().data[0].fe_degree + 1), + Utilities::fixed_power( + fe_eval.get_shape_info().data[0].n_q_points_1d)); + + const std::size_t n_q_points = fe_eval.get_shape_info().n_q_points; + const std::size_t dofs_per_comp = + fe_eval.get_shape_info().dofs_per_component_on_cell; + + // Initial shift depending on component (normal_dir) + Number *values_dofs = values_dofs_actual + dofs_per_comp * normal_dir; + Number *values_quad = fe_eval.begin_values() + n_q_points * normal_dir; + Number *gradients_quad = + fe_eval.begin_gradients() + dim * n_q_points * normal_dir; + Number *hessians_quad = + (dim == 2) ? fe_eval.begin_hessians() + 3 * n_q_points * normal_dir : + fe_eval.begin_hessians() + 6 * n_q_points * normal_dir; + + // Integrate path + switch (dim) + { + case 2: + if ((evaluation_flag & EvaluationFlags::values) && + !(evaluation_flag & EvaluationFlags::gradients)) + { + eval1.template values<1, false, false, normal_dir>(values_quad, + temp1); + if (add_into_values_array == false) + eval0.template values<0, false, false, normal_dir>(temp1, + values_dofs); + else + eval0.template values<0, false, true, normal_dir>(temp1, + values_dofs); + } + if (evaluation_flag & EvaluationFlags::gradients) + { + eval1.template gradients<1, false, false, normal_dir>( + gradients_quad + n_q_points, temp1); + if ((evaluation_flag & EvaluationFlags::values)) + eval1.template values<1, false, true, normal_dir>(values_quad, + temp1); + if (add_into_values_array == false) + eval0.template values<0, false, false, normal_dir>(temp1, + values_dofs); + else + eval0.template values<0, false, true, normal_dir>(temp1, + values_dofs); + eval1.template values<1, false, false, normal_dir>(gradients_quad, + temp1); + eval0.template gradients<0, false, true, normal_dir>(temp1, + values_dofs); + } + if (evaluation_flag & EvaluationFlags::hessians) + { + // grad xx + eval1.template values<1, false, false, normal_dir>(hessians_quad, + temp1); + + if ((evaluation_flag & EvaluationFlags::values) || + (evaluation_flag & EvaluationFlags::gradients) || + add_into_values_array == true) + eval0.template hessians<0, false, true, normal_dir>( + temp1, values_dofs); + else + eval0.template hessians<0, false, false, normal_dir>( + temp1, values_dofs); + + // grad yy + eval1.template hessians<1, false, false, normal_dir>( + hessians_quad + n_q_points, temp1); + eval0.template values<0, false, true, normal_dir>(temp1, + values_dofs); + + // grad xy + eval1.template gradients<1, false, false, normal_dir>( + hessians_quad + 2 * n_q_points, temp1); + eval0.template gradients<0, false, true, normal_dir>(temp1, + values_dofs); + } + break; + + case 3: + if ((evaluation_flag & EvaluationFlags::values) && + !(evaluation_flag & EvaluationFlags::gradients)) + { + eval2.template values<2, false, false, normal_dir>(values_quad, + temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + if (add_into_values_array == false) + eval0.template values<0, false, false, normal_dir>(temp2, + values_dofs); + else + eval0.template values<0, false, true, normal_dir>(temp2, + values_dofs); + } + if (evaluation_flag & EvaluationFlags::gradients) + { + eval2.template gradients<2, false, false, normal_dir>( + gradients_quad + 2 * n_q_points, temp1); + if ((evaluation_flag & EvaluationFlags::values)) + eval2.template values<2, false, true, normal_dir>(values_quad, + temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval2.template values<2, false, false, normal_dir>( + gradients_quad + n_q_points, temp1); + eval1.template gradients<1, false, true, normal_dir>(temp1, + temp2); + if (add_into_values_array == false) + eval0.template values<0, false, false, normal_dir>(temp2, + values_dofs); + else + eval0.template values<0, false, true, normal_dir>(temp2, + values_dofs); + eval2.template values<2, false, false, normal_dir>(gradients_quad, + temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval0.template gradients<0, false, true, normal_dir>(temp2, + values_dofs); + } + if (evaluation_flag & EvaluationFlags::hessians) + { + // grad xx + eval2.template values<2, false, false, normal_dir>(hessians_quad, + temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + + if ((evaluation_flag & EvaluationFlags::values) || + (evaluation_flag & EvaluationFlags::gradients) || + add_into_values_array == true) + eval0.template hessians<0, false, true, normal_dir>( + temp2, values_dofs); + else + eval0.template hessians<0, false, false, normal_dir>( + temp2, values_dofs); + + // grad yy + eval2.template values<2, false, false, normal_dir>(hessians_quad + + n_q_points, + temp1); + eval1.template hessians<1, false, false, normal_dir>(temp1, + temp2); + eval0.template values<0, false, true, normal_dir>(temp2, + values_dofs); + + // grad zz + eval2.template hessians<2, false, false, normal_dir>( + hessians_quad + 2 * n_q_points, temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval0.template values<0, false, true, normal_dir>(temp2, + values_dofs); + + // grad xy + eval2.template values<2, false, false, normal_dir>( + hessians_quad + 3 * n_q_points, temp1); + eval1.template gradients<1, false, false, normal_dir>(temp1, + temp2); + eval0.template gradients<0, false, true, normal_dir>(temp2, + values_dofs); + + // grad xz + eval2.template gradients<2, false, false, normal_dir>( + hessians_quad + 4 * n_q_points, temp1); + eval1.template values<1, false, false, normal_dir>(temp1, temp2); + eval0.template gradients<0, false, true, normal_dir>(temp2, + values_dofs); + + // grad yz + eval2.template gradients<2, false, false, normal_dir>( + hessians_quad + 5 * n_q_points, temp1); + eval1.template gradients<1, false, false, normal_dir>(temp1, + temp2); + eval0.template values<0, false, true, normal_dir>(temp2, + values_dofs); + } + + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } /** * This struct implements the change between two different bases. This is an @@ -1779,7 +2306,8 @@ namespace internal Assert(fe_eval.get_shape_info().data.size() == 1 || (fe_eval.get_shape_info().data.size() == dim && - element_type == ElementType::tensor_general), + element_type == ElementType::tensor_general) || + element_type == ElementType::tensor_raviart_thomas, ExcNotImplemented()); if (fe_degree >= 0 && fe_degree + 1 == n_q_points_1d && @@ -1847,6 +2375,18 @@ namespace internal values_dofs, fe_eval); } + else if (element_type == ElementType::tensor_raviart_thomas) + { + FEEvaluationImpl< + ElementType::tensor_raviart_thomas, + dim, + (fe_degree == -1) ? 1 : fe_degree, + (n_q_points_1d < 1) ? 1 : n_q_points_1d, + Number>::template evaluate_or_integrate(evaluation_flag, + const_cast( + values_dofs), + fe_eval); + } else { FEEvaluationImpl= 0 && fe_degree + 1 == n_q_points_1d && @@ -1973,6 +2514,19 @@ namespace internal fe_eval, sum_into_values_array); } + else if (element_type == ElementType::tensor_raviart_thomas) + { + FEEvaluationImpl:: + template evaluate_or_integrate(integration_flag, + const_cast( + values_dofs), + fe_eval, + sum_into_values_array); + } else { FEEvaluationImpl + struct FEFaceEvaluationImplRaviartThomas + { + using EvalNormal = + EvaluatorTensorProduct; + using EvalTangent = + EvaluatorTensorProduct; + + using EvalGeneral = EvaluatorTensorProduct; + template + static EvalType + create_evaluator_tensor_product( + const MatrixFreeFunctions::UnivariateShapeData &data, + const unsigned int subface_index, + const unsigned int direction) + { + if (subface_index >= GeometryInfo::max_children_per_cell) + return EvalType(data.shape_values, + data.shape_gradients, + data.shape_hessians); + else + { + Assert(false, ExcNotImplemented()); + + const unsigned int index = + direction == 0 ? subface_index % 2 : subface_index / 2; + return EvalType(data.values_within_subface[index], + data.gradients_within_subface[index], + data.hessians_within_subface[index]); + } + } + + template + static void + evaluate_or_integrate_in_face( + const EvaluationFlags::EvaluationFlags evaluation_flag, + const MatrixFreeFunctions::ShapeInfo &shape_info, + Number * values_dofs, + FEEvaluationData & fe_eval, + Number * scratch_data, + const unsigned int subface_index, + const unsigned int face_no) + { + // TODO. Make sure hanging nodes also are supported. + // The following part probably needs a rethink. + EvalNormal eval_normal = + create_evaluator_tensor_product(shape_info.data.front(), + subface_index, + 0); + EvalTangent eval_tangent = + create_evaluator_tensor_product(shape_info.data.back(), + subface_index, + 1); + + // Used for normal faces which are isotropic + EvalGeneral eval_general = + create_evaluator_tensor_product(shape_info.data.back(), + subface_index, + 0); + + // Note, n_dofs on tangent face + const std::size_t n_dofs_tangent = shape_info.dofs_per_component_on_face; + const std::size_t n_dofs_normal = + n_dofs_tangent - Utilities::pow(fe_degree, dim - 2); + + const unsigned int face_direction = face_no / 2; + + if (face_direction == 0) + { + evaluate_in_face_apply<-1, 0>( + eval_general, + eval_general, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_normal, + std::integral_constant()); + + values_dofs += 3 * n_dofs_normal; + + evaluate_in_face_apply<0, 1>( + eval_normal, + eval_tangent, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + + values_dofs += 3 * n_dofs_tangent; + + if (dim == 3) + { + evaluate_in_face_apply<1, 2>( + eval_tangent, + eval_normal, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + } + } + else if (face_direction == 1) + { + // NOTE. Take zx-coordinates into account for dim == 3 + if (dim == 3) + evaluate_in_face_apply<1, 0>( + eval_tangent, + eval_normal, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + else + evaluate_in_face_apply<0, 0>( + eval_normal, + eval_tangent, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + + values_dofs += 3 * n_dofs_tangent; + + evaluate_in_face_apply<-1, 1>( + eval_general, + eval_general, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_normal, + std::integral_constant()); + + values_dofs += 3 * n_dofs_normal; + + if (dim == 3) + { + // NOTE. Take zx-coordinates into account + evaluate_in_face_apply<0, 2>( + eval_normal, + eval_tangent, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + } + } + else + { + evaluate_in_face_apply<0, 0>( + eval_normal, + eval_tangent, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + + values_dofs += 3 * n_dofs_tangent; + + evaluate_in_face_apply<1, 1>( + eval_tangent, + eval_normal, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_tangent, + std::integral_constant()); + + values_dofs += 3 * n_dofs_tangent; + + if (dim == 3) + { + evaluate_in_face_apply<-1, 2>( + eval_general, + eval_general, + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + n_dofs_normal, + std::integral_constant()); + } + } + } + + /* + * Helper function which applies the 1D kernels for on one + * component in a face. normal_dir indicates the direction of the continuous + * component of the RT space. std::integral_constant is the + * evaluation path, and std::integral_constant bellow is the + * integration path. + */ + template + static inline void + evaluate_in_face_apply( + const Eval0 & eval0, + const Eval1 & eval1, + Number * values_dofs, + FEEvaluationData & fe_eval, + Number * scratch_data, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const std::size_t dofs_stride, + std::integral_constant) + { + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + + Number *values_quad = fe_eval.begin_values() + n_q_points * component; + Number *gradients_quad = + fe_eval.begin_gradients() + dim * n_q_points * component; + Number *hessians_quad = + fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component; + + // Evaluation path + + if ((evaluation_flag & EvaluationFlags::values) && + !(evaluation_flag & EvaluationFlags::gradients)) + { + switch (dim) + { + case 3: + eval0.template values<0, true, false, normal_dir>(values_dofs, + values_quad); + eval1.template values<1, true, false, normal_dir>(values_quad, + values_quad); + break; + case 2: + eval0.template values<0, true, false, normal_dir>(values_dofs, + values_quad); + break; + default: + Assert(false, ExcNotImplemented()); + } + } + else if (evaluation_flag & EvaluationFlags::gradients) + { + switch (dim) + { + case 3: + // grad x + eval0.template gradients<0, true, false, normal_dir>( + values_dofs, scratch_data); + eval1.template values<1, true, false, normal_dir>( + scratch_data, gradients_quad); + + // grad y + eval0.template values<0, true, false, normal_dir>(values_dofs, + scratch_data); + eval1.template gradients<1, true, false, normal_dir>( + scratch_data, gradients_quad + n_q_points); + + if (evaluation_flag & EvaluationFlags::values) + eval1.template values<1, true, false, normal_dir>( + scratch_data, values_quad); + + // grad z + eval0.template values<0, true, false, normal_dir>(values_dofs + + dofs_stride, + scratch_data); + eval1.template values<1, true, false, normal_dir>( + scratch_data, gradients_quad + 2 * n_q_points); + + break; + case 2: + eval0.template values<0, true, false, normal_dir>( + values_dofs + dofs_stride, gradients_quad + n_q_points); + eval0.template gradients<0, true, false, normal_dir>( + values_dofs, gradients_quad); + if ((evaluation_flag & EvaluationFlags::values)) + eval0.template values<0, true, false, normal_dir>( + values_dofs, values_quad); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + + if (evaluation_flag & EvaluationFlags::hessians) + { + switch (dim) + { + case 3: + // grad xx + eval0.template hessians<0, true, false, normal_dir>( + values_dofs, scratch_data); + eval1.template values<1, true, false, normal_dir>( + scratch_data, hessians_quad); + + // grad yy + eval0.template values<0, true, false, normal_dir>(values_dofs, + scratch_data); + eval1.template hessians<1, true, false, normal_dir>( + scratch_data, hessians_quad + n_q_points); + + // grad zz + eval0.template values<0, true, false, normal_dir>( + values_dofs + 2 * dofs_stride, scratch_data); + eval1.template values<1, true, false, normal_dir>( + scratch_data, hessians_quad + 2 * n_q_points); + + // grad xy + eval0.template gradients<0, true, false, normal_dir>( + values_dofs, scratch_data); + eval1.template gradients<1, true, false, normal_dir>( + scratch_data, hessians_quad + 3 * n_q_points); + + // grad xz + eval0.template gradients<0, true, false, normal_dir>( + values_dofs + dofs_stride, scratch_data); + eval1.template values<1, true, false, normal_dir>( + scratch_data, hessians_quad + 4 * n_q_points); + + // grad yz + eval0.template values<0, true, false, normal_dir>(values_dofs + + dofs_stride, + scratch_data); + eval1.template gradients<1, true, false, normal_dir>( + scratch_data, hessians_quad + 5 * n_q_points); + + break; + case 2: + // grad xx + eval0.template hessians<0, true, false, normal_dir>( + values_dofs, hessians_quad); + // grad yy + eval0.template values<0, true, false, normal_dir>( + values_dofs + 2 * dofs_stride, hessians_quad + n_q_points); + // grad xy + eval0.template gradients<0, true, false, normal_dir>( + values_dofs + dofs_stride, hessians_quad + 2 * n_q_points); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + } + + template + static inline void + evaluate_in_face_apply( + const Eval0 & eval0, + const Eval1 & eval1, + Number * values_dofs, + FEEvaluationData & fe_eval, + Number * scratch_data, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const std::size_t dofs_stride, + std::integral_constant) + { + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + + Number *values_quad = fe_eval.begin_values() + n_q_points * component; + Number *gradients_quad = + fe_eval.begin_gradients() + dim * n_q_points * component; + Number *hessians_quad = + fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component; + + // Integration path + if ((evaluation_flag & EvaluationFlags::values) && + !(evaluation_flag & EvaluationFlags::gradients)) + { + switch (dim) + { + case 3: + eval1.template values<1, false, false, normal_dir>(values_quad, + values_quad); + eval0.template values<0, false, false, normal_dir>(values_quad, + values_dofs); + break; + case 2: + eval0.template values<0, false, false, normal_dir>(values_quad, + values_dofs); + break; + default: + Assert(false, ExcNotImplemented()); + } + } + else if (evaluation_flag & EvaluationFlags::gradients) + { + switch (dim) + { + case 3: + // grad z + eval1.template values<1, false, false, normal_dir>( + gradients_quad + 2 * n_q_points, + gradients_quad + 2 * n_q_points); + eval0.template values<0, false, false, normal_dir>( + gradients_quad + 2 * n_q_points, values_dofs + dofs_stride); + + if (evaluation_flag & EvaluationFlags::values) + { + eval1.template values<1, false, false, normal_dir>( + values_quad, scratch_data); + eval1.template gradients<1, false, true, normal_dir>( + gradients_quad + n_q_points, scratch_data); + } + else + eval1.template gradients<1, false, false, normal_dir>( + gradients_quad + n_q_points, scratch_data); + + // grad y + eval0.template values<0, false, false, normal_dir>(scratch_data, + values_dofs); + + // grad x + eval1.template values<1, false, false, normal_dir>( + gradients_quad, scratch_data); + eval0.template gradients<0, false, true, normal_dir>( + scratch_data, values_dofs); + + break; + case 2: + eval0.template values<0, false, false, normal_dir>( + gradients_quad + n_q_points, values_dofs + dofs_stride); + eval0.template gradients<0, false, false, normal_dir>( + gradients_quad, values_dofs); + if (evaluation_flag & EvaluationFlags::values) + eval0.template values<0, false, true, normal_dir>( + values_quad, values_dofs); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + + if (evaluation_flag & EvaluationFlags::hessians) + { + switch (dim) + { + case 3: + // grad xx + eval1.template values<1, false, false, normal_dir>( + hessians_quad, scratch_data); + if ((evaluation_flag & + (EvaluationFlags::values | EvaluationFlags::gradients))) + eval0.template hessians<0, false, true, normal_dir>( + scratch_data, values_dofs); + else + eval0.template hessians<0, false, false, normal_dir>( + scratch_data, values_dofs); + + // grad yy + eval1.template hessians<1, false, false, normal_dir>( + hessians_quad + n_q_points, scratch_data); + eval0.template values<0, false, true, normal_dir>(scratch_data, + values_dofs); + + // grad zz + eval1.template values<1, false, false, normal_dir>( + hessians_quad + 2 * n_q_points, scratch_data); + eval0.template values<0, false, false, normal_dir>( + scratch_data, values_dofs + 2 * dofs_stride); + + // grad xy + eval1.template gradients<1, false, false, normal_dir>( + hessians_quad + 3 * n_q_points, scratch_data); + eval0.template gradients<0, false, true, normal_dir>( + scratch_data, values_dofs); + + // grad xz + eval1.template values<1, false, false, normal_dir>( + hessians_quad + 4 * n_q_points, scratch_data); + if ((evaluation_flag & EvaluationFlags::gradients)) + eval0.template gradients<0, false, true, normal_dir>( + scratch_data, values_dofs + dofs_stride); + else + eval0.template gradients<0, false, false, normal_dir>( + scratch_data, values_dofs + dofs_stride); + + // grad yz + eval1.template gradients<1, false, false, normal_dir>( + hessians_quad + 5 * n_q_points, scratch_data); + eval0.template values<0, false, true, normal_dir>( + scratch_data, values_dofs + dofs_stride); + + break; + case 2: + // grad xx + if (evaluation_flag & + (EvaluationFlags::values | EvaluationFlags::gradients)) + eval0.template hessians<0, false, true, normal_dir>( + hessians_quad, values_dofs); + else + eval0.template hessians<0, false, false, normal_dir>( + hessians_quad, values_dofs); + + // grad yy + eval0.template values<0, false, false, normal_dir>( + hessians_quad + n_q_points, values_dofs + 2 * dofs_stride); + // grad xy + if ((evaluation_flag & EvaluationFlags::gradients)) + eval0.template gradients<0, false, true, normal_dir>( + hessians_quad + 2 * n_q_points, values_dofs + dofs_stride); + else + eval0.template gradients<0, false, false, normal_dir>( + hessians_quad + 2 * n_q_points, values_dofs + dofs_stride); + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + } + }; template struct FEFaceNormalEvaluationImpl { + using EvalNormal = + EvaluatorTensorProduct; + using EvalTangent = + EvaluatorTensorProduct; + template static void interpolate(const unsigned int n_components, @@ -2495,17 +3590,20 @@ namespace internal shape_info.data.front().fe_degree || fe_degree == -1, ExcInternalError()); - - interpolate_generic( - n_components, - input, - output, - flags, - face_no, - shape_info.data.front().fe_degree + 1, - shape_info.data.front().shape_data_on_face, - shape_info.dofs_per_component_on_cell, - 3 * shape_info.dofs_per_component_on_face); + if (shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + interpolate_generic_raviart_thomas( + n_components, input, output, flags, face_no, shape_info); + else + interpolate_generic( + n_components, + input, + output, + flags, + face_no, + shape_info.data.front().fe_degree + 1, + shape_info.data.front().shape_data_on_face, + shape_info.dofs_per_component_on_cell, + 3 * shape_info.dofs_per_component_on_face); } /** @@ -2611,6 +3709,168 @@ namespace internal dofs_per_component_on_face); } } + + private: + template + static EvalType + create_evaluator_tensor_product( + const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int face_no) + { + // Shape_data is of size 2 with normal direction in [0] and tangential + // direction in [1]. See shape_info.template.h + const MatrixFreeFunctions::UnivariateShapeData *shape_data; + if (std::is_same::value) + { + shape_data = &shape_info.data[0]; + } + else + { + shape_data = &shape_info.data[1]; + } + return EvalType(shape_data->shape_data_on_face[face_no % 2], + AlignedVector(), + AlignedVector()); + } + + template + static void + interpolate_generic_raviart_thomas( + const unsigned int n_components, + const Number * input, + Number * output, + const EvaluationFlags::EvaluationFlags flag, + const unsigned int face_no, + const MatrixFreeFunctions::ShapeInfo &shape_info) + { + bool increase_max_der = false; + if ((flag & EvaluationFlags::hessians && max_derivative < 2) || + (flag & EvaluationFlags::gradients && max_derivative < 1)) + increase_max_der = true; + + if (face_direction == face_no / 2 && !increase_max_der) + { + interpolate_generic_raviart_thomas_apply_face( + shape_info, face_no, input, output); + } + else if (face_direction == face_no / 2) + { + // Only increase max_derivative + interpolate_generic_raviart_thomas( + n_components, input, output, flag, face_no, shape_info); + } + else if (face_direction < dim) + { + if (increase_max_der) + { + interpolate_generic_raviart_thomas< + do_evaluate, + add_into_output, + std::min(face_direction + 1, dim - 1), + std::min(max_derivative + 1, 2)>( + n_components, input, output, flag, face_no, shape_info); + } + else + { + interpolate_generic_raviart_thomas( + n_components, input, output, flag, face_no, shape_info); + } + } + } + + /* Help function for interpolate_generic_raviart_thomas */ + template + static inline void + interpolate_generic_raviart_thomas_apply_face( + const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int face_no, + const Number * input, + Number * output) + { + using Evalf0 = typename std:: + conditional::type; + using Evalf1 = typename std:: + conditional::type; + using Evalf2 = typename std:: + conditional::type; + + Evalf0 evalf0 = + create_evaluator_tensor_product(shape_info, face_no); + Evalf1 evalf1 = + create_evaluator_tensor_product(shape_info, face_no); + Evalf2 evalf2 = + create_evaluator_tensor_product(shape_info, face_no); + + const unsigned int dofs_per_component_on_cell = + shape_info.dofs_per_component_on_cell; + const unsigned int dofs_per_component_on_face = + 3 * shape_info.dofs_per_component_on_face; + + // NOTE! dofs_per_component_on_face is in the tangent direction, + // i.e (fe.degree+1)*fe.degree. Normal faces are only + // fe.degree*fe.degree + const unsigned int in_stride = + do_evaluate ? dofs_per_component_on_cell : dofs_per_component_on_face; + const unsigned int out_stride = + do_evaluate ? dofs_per_component_on_face : dofs_per_component_on_cell; + + const unsigned int in_stride_after_normal = + do_evaluate ? + dofs_per_component_on_cell : + dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2); + const unsigned int out_stride_after_normal = + do_evaluate ? + dofs_per_component_on_face - 3 * Utilities::pow(fe_degree, dim - 2) : + dofs_per_component_on_cell; + + evalf0.template apply_face(input, output); + // stride to next component + input += (face_direction == 0) ? in_stride_after_normal : in_stride; + output += (face_direction == 0) ? out_stride_after_normal : out_stride; + + evalf1.template apply_face(input, output); + + if (dim == 3) + { + // stride to next component + input += (face_direction == 1) ? in_stride_after_normal : in_stride; + output += + (face_direction == 1) ? out_stride_after_normal : out_stride; + + evalf2.template apply_face(input, output); + } + } }; @@ -2949,9 +4209,25 @@ namespace internal constexpr unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0; - if (fe_degree > -1 && - subface_index >= GeometryInfo::max_children_per_cell && - shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) + if (fe_degree >= 1 && + shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + { + FEFaceEvaluationImplRaviartThomas:: + template evaluate_or_integrate_in_face(evaluation_flag, + shape_info, + temp, + fe_eval, + scratch_data, + subface_index, + face_no); + } + else if (fe_degree > -1 && + subface_index >= GeometryInfo::max_children_per_cell && + shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) FEFaceEvaluationImpl -1 ? n_q_points_1d : 0; const unsigned int subface_index = fe_eval.get_subface_index(); - if (fe_degree > -1 && - fe_eval.get_subface_index() >= - GeometryInfo::max_children_per_cell && - shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) + if (fe_degree >= 1 && + shape_info.element_type == MatrixFreeFunctions::tensor_raviart_thomas) + { + FEFaceEvaluationImplRaviartThomas:: + template evaluate_or_integrate_in_face(integration_flag, + shape_info, + temp, + fe_eval, + scratch_data, + subface_index, + face_no); + } + else if (fe_degree > -1 && + fe_eval.get_subface_index() >= + GeometryInfo::max_children_per_cell && + shape_info.element_type <= MatrixFreeFunctions::tensor_symmetric) FEFaceEvaluationImpl< true, dim, diff --git a/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h new file mode 100644 index 0000000000..fa068963ba --- /dev/null +++ b/include/deal.II/matrix_free/evaluation_template_face_factory.templates.h @@ -0,0 +1,125 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#ifndef dealii_matrix_free_evaluation_template_factory_templates_h +#define dealii_matrix_free_evaluation_template_factory_templates_h + + +#include + +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace internal +{ + template + void + FEFaceEvaluationFactory::evaluate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * values_dofs, + FEEvaluationData & fe_eval) + { + instantiation_helper_run<1, + FEFaceEvaluationImplEvaluateSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + evaluation_flag, + values_dofs, + fe_eval); + } + + + + template + void + FEFaceEvaluationFactory::integrate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * values_dofs, + FEEvaluationData & fe_eval) + { + instantiation_helper_run< + 1, + FEFaceEvaluationImplIntegrateSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + integration_flag, + values_dofs, + fe_eval); + } + + + + template + void + FEFaceEvaluationGatherFactory::evaluate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags evaluation_flag, + const Number * src_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &fe_eval) + { + instantiation_helper_run< + 1, + FEFaceEvaluationImplGatherEvaluateSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + evaluation_flag, + src_ptr, + sm_ptr, + fe_eval); + } + + + + template + void + FEFaceEvaluationGatherFactory::integrate( + const unsigned int n_components, + const EvaluationFlags::EvaluationFlags integration_flag, + Number * dst_ptr, + const std::vector> * sm_ptr, + FEEvaluationData &fe_eval) + { + instantiation_helper_run< + 1, + FEFaceEvaluationImplIntegrateScatterSelector>( + fe_eval.get_shape_info().data[0].fe_degree, + fe_eval.get_shape_info().data[0].n_q_points_1d, + n_components, + integration_flag, + dst_ptr, + sm_ptr, + fe_eval); + } + +} // end of namespace internal + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/include/deal.II/matrix_free/evaluation_template_factory.templates.h b/include/deal.II/matrix_free/evaluation_template_factory.templates.h index cdaca44458..d70ba5758d 100644 --- a/include/deal.II/matrix_free/evaluation_template_factory.templates.h +++ b/include/deal.II/matrix_free/evaluation_template_factory.templates.h @@ -93,97 +93,6 @@ namespace internal - template - void - FEFaceEvaluationFactory::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * values_dofs, - FEEvaluationData & fe_eval) - { - instantiation_helper_run<1, - FEFaceEvaluationImplEvaluateSelector>( - fe_eval.get_shape_info().data[0].fe_degree, - fe_eval.get_shape_info().data[0].n_q_points_1d, - n_components, - evaluation_flag, - values_dofs, - fe_eval); - } - - - - template - void - FEFaceEvaluationFactory::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * values_dofs, - FEEvaluationData & fe_eval) - { - instantiation_helper_run< - 1, - FEFaceEvaluationImplIntegrateSelector>( - fe_eval.get_shape_info().data[0].fe_degree, - fe_eval.get_shape_info().data[0].n_q_points_1d, - n_components, - integration_flag, - values_dofs, - fe_eval); - } - - - - template - void - FEFaceEvaluationGatherFactory::evaluate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags evaluation_flag, - const Number * src_ptr, - const std::vector> * sm_ptr, - FEEvaluationData &fe_eval) - { - instantiation_helper_run< - 1, - FEFaceEvaluationImplGatherEvaluateSelector>( - fe_eval.get_shape_info().data[0].fe_degree, - fe_eval.get_shape_info().data[0].n_q_points_1d, - n_components, - evaluation_flag, - src_ptr, - sm_ptr, - fe_eval); - } - - - - template - void - FEFaceEvaluationGatherFactory::integrate( - const unsigned int n_components, - const EvaluationFlags::EvaluationFlags integration_flag, - Number * dst_ptr, - const std::vector> * sm_ptr, - FEEvaluationData &fe_eval) - { - instantiation_helper_run< - 1, - FEFaceEvaluationImplIntegrateScatterSelector>( - fe_eval.get_shape_info().data[0].fe_degree, - fe_eval.get_shape_info().data[0].n_q_points_1d, - n_components, - integration_flag, - dst_ptr, - sm_ptr, - fe_eval); - } - - - template void CellwiseInverseMassFactory::apply( diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index f406690258..4f1ce43dd9 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1160,16 +1160,21 @@ namespace internal dof_info[no].fe_index_conversion[fe_index].clear(); for (unsigned int c = 0; c < dof_info[no].n_base_elements; ++c) { - dof_info[no].n_components[c] = fe.element_multiplicity(c); + dof_info[no].n_components[c] = + fe.element_multiplicity(c) * + fe.base_element(c).n_components(); for (unsigned int l = 0; l < dof_info[no].n_components[c]; ++l) { dof_info[no].component_to_base_index.push_back(c); dof_info[no] .component_dof_indices_offset[fe_index] - .push_back(dof_info[no] - .component_dof_indices_offset[fe_index] - .back() + - fe.base_element(c).n_dofs_per_cell()); + .push_back( + dof_info[no] + .component_dof_indices_offset[fe_index] + .back() + + shape_infos(dof_info[no].global_base_element_offset + c, + fe_index) + .dofs_per_component_on_cell); dof_info[no].fe_index_conversion[fe_index].push_back( fe.base_element(c).degree); } diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h index 1187d0708c..ee7394ee58 100644 --- a/include/deal.II/matrix_free/shape_info.h +++ b/include/deal.II/matrix_free/shape_info.h @@ -100,9 +100,18 @@ namespace internal tensor_symmetric_plus_dg0 = 5, /** - * Shape functions without an tensor product properties. + * Special case of the FE_RaviartThomasNodal element with anisotropic + * tensor product shape functions, i.e degree (k + 1) in normal direction, + * and k in tangential direction. */ - tensor_none = 6 + tensor_raviart_thomas = 6, + + /** + * Shape functions without a tensor product properties. + */ + tensor_none = 7 + + }; diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h index 9f1a084a81..63a2d6e92a 100644 --- a/include/deal.II/matrix_free/shape_info.templates.h +++ b/include/deal.II/matrix_free/shape_info.templates.h @@ -33,6 +33,7 @@ #include #include #include +#include #include #include #include @@ -170,8 +171,6 @@ namespace internal } } - - template std::unique_ptr> create_fe(const FiniteElement &fe) @@ -194,7 +193,6 @@ namespace internal } - // ----------------- actual ShapeInfo implementation -------------------- template @@ -238,15 +236,171 @@ namespace internal { static_assert(dim == spacedim, "Currently, only the case dim=spacedim is implemented"); - if (quad_in.is_tensor_product() == false || - dynamic_cast *>( - &fe_in.base_element(base_element_number)) || - dynamic_cast *>( - &fe_in.base_element(base_element_number)) || - dynamic_cast *>( - &fe_in.base_element(base_element_number)) || - dynamic_cast *>( + + // ShapeInfo for RT elements. Here, data is of size 2 instead of 1. + // data[0] is univariate_shape_data in normal direction and + // data[1] is univariate_shape_data in tangential direction + // + if (dynamic_cast *>( &fe_in.base_element(base_element_number))) + { + element_type = tensor_raviart_thomas; + + const auto quad = quad_in.get_tensor_basis()[0]; + + const FiniteElement &fe = + fe_in.base_element(base_element_number); + n_dimensions = dim; + n_components = fe_in.n_components(); + + data.resize(2); + const unsigned int n_q_points_1d = quad.size(); + + n_q_points = Utilities::fixed_power(n_q_points_1d); + n_q_points_face = Utilities::fixed_power(n_q_points_1d); + + dofs_per_component_on_cell = fe_in.n_dofs_per_cell() / n_components; + + // NOTE dofs_per_component_on_face is in tangential direction! + dofs_per_component_on_face = + fe_in.n_dofs_per_face() + Utilities::pow(fe_in.degree, dim - 2); + const unsigned int dofs_per_face_normal = fe_in.n_dofs_per_face(); + + lexicographic_numbering = + get_lexicographic_numbering_rt_nodal(fe_in.degree); + + // To get the right shape_values of the RT element + std::vector lex_normal, lex_tangent; + for (unsigned int i = 0; i < fe.degree; ++i) + lex_tangent.push_back(i); + + lex_normal.push_back(0); + for (unsigned int i = dofs_per_face_normal * 2 * dim; + i < dofs_per_face_normal * 2 * dim + fe.degree - 1; + ++i) + lex_normal.push_back(i); + lex_normal.push_back(dofs_per_face_normal); + + // 'direction' distingusishes between normal and tangential direction + for (unsigned int direction = 0; direction < 2; ++direction) + { + UnivariateShapeData &univariate_shape_data = + (direction == 0) ? data.front() : data.back(); + + univariate_shape_data.element_type = tensor_raviart_thomas; + univariate_shape_data.quadrature = quad; + univariate_shape_data.n_q_points_1d = n_q_points_1d; + univariate_shape_data.fe_degree = fe.degree - direction; + + // grant write access to common univariate shape data + auto &shape_values = univariate_shape_data.shape_values; + auto &shape_gradients = univariate_shape_data.shape_gradients; + auto &shape_hessians = univariate_shape_data.shape_hessians; + + auto &values_within_subface = + univariate_shape_data.values_within_subface; + auto &gradients_within_subface = + univariate_shape_data.gradients_within_subface; + auto &hessians_within_subface = + univariate_shape_data.hessians_within_subface; + + auto &shape_data_on_face = + univariate_shape_data.shape_data_on_face; + + const unsigned int n_dofs_1d = fe.degree + 1 - direction; + const unsigned int array_size = n_dofs_1d * n_q_points_1d; + + shape_gradients.resize_fast(array_size); + shape_values.resize_fast(array_size); + shape_hessians.resize_fast(array_size); + + values_within_subface[0].resize(array_size); + values_within_subface[1].resize(array_size); + gradients_within_subface[0].resize(array_size); + gradients_within_subface[1].resize(array_size); + hessians_within_subface[0].resize(array_size); + hessians_within_subface[1].resize(array_size); + + shape_data_on_face[0].resize(3 * n_dofs_1d); + shape_data_on_face[1].resize(3 * n_dofs_1d); + + Point unit_point; + for (unsigned int i = 0; i < n_dofs_1d; ++i) + { + // need to reorder from hierarchical to lexicographic to get + // the DoFs correct + const unsigned int my_i = + (direction == 0) ? lex_normal[i] : lex_tangent[i]; + for (unsigned int q = 0; q < n_q_points_1d; ++q) + { + Point q_point = unit_point; + q_point[direction] = quad.get_points()[q][0]; + + shape_values[i * n_q_points_1d + q] = + fe.shape_value_component(my_i, q_point, 0); + shape_gradients[i * n_q_points_1d + q] = + fe.shape_grad_component(my_i, q_point, 0)[direction]; + shape_hessians[i * n_q_points_1d + q] = + fe.shape_grad_grad_component(my_i, + q_point, + 0)[direction][direction]; + + // evaluate basis functions on the two 1D subfaces (i.e., + // at the positions divided by one half and shifted by one + // half, respectively) for hanging nodes + q_point[direction] *= 0.5; + values_within_subface[0][i * n_q_points_1d + q] = + fe.shape_value_component(my_i, q_point, 0); + gradients_within_subface[0][i * n_q_points_1d + q] = + fe.shape_grad_component(my_i, q_point, 0)[direction]; + hessians_within_subface[0][i * n_q_points_1d + q] = + fe.shape_grad_grad_component(my_i, + q_point, + 0)[direction][direction]; + q_point[direction] += 0.5; + values_within_subface[1][i * n_q_points_1d + q] = + fe.shape_value_component(my_i, q_point, 0); + gradients_within_subface[1][i * n_q_points_1d + q] = + fe.shape_grad_component(my_i, q_point, 0)[direction]; + hessians_within_subface[1][i * n_q_points_1d + q] = + fe.shape_grad_grad_component(my_i, + q_point, + 0)[direction][direction]; + } + // evaluate basis functions on the 1D faces, i.e., in zero and + // one + Point q_point = unit_point; + q_point[direction] = 0; + shape_data_on_face[0][i] = + fe.shape_value_component(my_i, q_point, 0); + shape_data_on_face[0][i + n_dofs_1d] = + fe.shape_grad_component(my_i, q_point, 0)[direction]; + shape_data_on_face[0][i + 2 * n_dofs_1d] = + fe.shape_grad_grad_component(my_i, + q_point, + 0)[direction][direction]; + q_point[direction] = 1; + shape_data_on_face[1][i] = + fe.shape_value_component(my_i, q_point, 0); + shape_data_on_face[1][i + n_dofs_1d] = + fe.shape_grad_component(my_i, q_point, 0)[direction]; + shape_data_on_face[1][i + 2 * n_dofs_1d] = + fe.shape_grad_grad_component(my_i, + q_point, + 0)[direction][direction]; + } + } + return; + } + else if (quad_in.is_tensor_product() == false || + dynamic_cast *>( + &fe_in.base_element(base_element_number)) || + dynamic_cast *>( + &fe_in.base_element(base_element_number)) || + dynamic_cast *>( + &fe_in.base_element(base_element_number)) || + dynamic_cast *>( + &fe_in.base_element(base_element_number))) { // specialization for arbitrary finite elements and quadrature rules // as needed in the context, e.g., of simplices diff --git a/include/deal.II/matrix_free/tensor_product_kernels.h b/include/deal.II/matrix_free/tensor_product_kernels.h index 608367b3b2..2e8885f01f 100644 --- a/include/deal.II/matrix_free/tensor_product_kernels.h +++ b/include/deal.II/matrix_free/tensor_product_kernels.h @@ -63,7 +63,11 @@ namespace internal * coefficient arrays. See the documentation of the EvaluatorTensorProduct * specialization for more information. */ - evaluate_symmetric_hierarchical + evaluate_symmetric_hierarchical, + /** + * Raviart-Thomas elements with anisotropic polynomials. + */ + evaluate_raviart_thomas }; @@ -194,21 +198,30 @@ namespace internal (void)dummy2; } - template + template void values(const Number in[], Number out[]) const { apply(shape_values, in, out); } - template + template void gradients(const Number in[], Number out[]) const { apply(shape_gradients, in, out); } - template + template void hessians(const Number in[], Number out[]) const { @@ -2328,6 +2341,482 @@ namespace internal + /** + * Internal evaluator for shape function in 2D and 3D using the + * tensor product form of the anisotropic basis functions of the + * raviart-thomas element, with degree k+1 in normal direction and + * k in tangential direction. + * + * @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 constexpr unsigned int n_rows_of_product = + numbers::invalid_unsigned_int; + static constexpr unsigned int n_columns_of_product = + numbers::invalid_unsigned_int; + + /** + * Empty constructor. Does nothing. Be careful when using 'values' and + * related methods because they need to be filled with the other pointer + */ + EvaluatorTensorProduct() + : shape_values(nullptr) + , shape_gradients(nullptr) + , shape_hessians(nullptr) + {} + + /** + * 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) + : 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 + void + values(const Number in[], Number out[]) const + { + apply(shape_values, + in, + out); + } + + template + void + gradients(const Number in[], Number out[]) const + { + apply(shape_gradients, + in, + out); + } + + template + void + hessians(const Number in[], Number out[]) const + { + 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 normal_dir Indicates the direction of the continuous component of the + * RT space in terms of the normal onto the face, e.g + * 0 if the is in x-direction, 1 if in y-direction + * etc. + * @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); + + 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; + }; + + template + template + inline void + 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); + constexpr int mm = contract_over_rows ? n_rows : n_columns, + nn = contract_over_rows ? n_columns : n_rows; + + constexpr int stride = Utilities::pow(n_columns, direction); + constexpr int n_blocks1 = one_line ? 1 : stride; + + // The number of blocks depend on both direction and dimension. + constexpr int n_blocks2 = + (dim - direction - 1 == 0) ? + 1 : + ((direction == normal_dir) ? + Utilities::pow((n_rows - 1), + (direction >= dim) ? 0 : dim - direction - 1) : + (((direction < normal_dir) ? (n_rows + 1) : n_rows) * + ((dim - direction == 3) ? n_rows : 1))); + + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + Number x[mm]; + for (int i = 0; i < mm; ++i) + x[i] = in[stride * i]; + + for (int col = 0; col < nn; ++col) + { + Number2 val0; + + if (contract_over_rows) + val0 = shape_data[col]; + else + val0 = shape_data[col * n_columns]; + + Number res0 = val0 * x[0]; + for (int i = 1; i < mm; ++i) + { + if (contract_over_rows) + val0 = shape_data[i * n_columns + col]; + else + val0 = shape_data[col * n_columns + i]; + + res0 += val0 * x[i]; + } + if (add) + out[stride * col] += res0; + + else + out[stride * col] = res0; + } + + if (one_line == false) + { + ++in; + ++out; + } + } + if (one_line == false) + { + in += stride * (mm - 1); + out += stride * (nn - 1); + } + } + } + + template + template + inline void + EvaluatorTensorProduct::apply_face(const Number *DEAL_II_RESTRICT in, + Number *DEAL_II_RESTRICT + out) const + { + Assert(dim > 1 && (lex_faces || dim < 4), + ExcMessage("Only dim=2,3 supported")); + static_assert(max_derivative >= 0 && max_derivative < 3, + "Only derivative orders 0-2 implemented"); + static_assert(!lex_faces, + "lex_faces = True is not implemented for Raviart-Thomas"); + Assert(shape_values != nullptr, + ExcMessage( + "The given array shape_values must not be the null pointer.")); + + // Determine the number of blocks depending on the face and normaldirection, + // as well as dimension. + constexpr int n_blocks1 = (face_direction == normal_direction) ? + (n_rows - 1) : + ((face_direction == 0 && normal_direction == 2) || + (face_direction == 1 && normal_direction == 2) || + (face_direction == 2 && normal_direction == 1)) ? + n_rows : + (n_rows + 1); + constexpr int n_blocks2 = + (dim == 2) ? 1 : + ((face_direction == normal_direction) ? + (n_rows - 1) : + (((face_direction == 0 && normal_direction == 1) || + (face_direction == 1 && normal_direction == 0) || + (face_direction == 2 && normal_direction == 0)) ? + n_rows : + (n_rows + 1))); + + AssertIndexRange(face_direction, dim); + + constexpr int stride = + (face_direction == normal_direction) ? + Utilities::pow(n_rows - 1, face_direction) : + ((face_direction == 0) ? + 1 : + ((face_direction == 2) ? + n_rows * (n_rows + 1) : + ((face_direction == 1 && normal_direction == 0) ? (n_rows + 1) : + n_rows))); + constexpr int out_stride = n_blocks1 * n_blocks2; + + const Number *DEAL_II_RESTRICT shape_values = this->shape_values; + + for (int i2 = 0; i2 < n_blocks2; ++i2) + { + for (int i1 = 0; i1 < n_blocks1; ++i1) + { + if (contract_onto_face == true) + { + Number res0 = shape_values[0] * in[0]; + Number res1, res2; + + if (max_derivative > 0) + res1 = shape_values[n_rows] * in[0]; + + if (max_derivative > 1) + res2 = shape_values[2 * n_rows] * in[0]; + + for (int ind = 1; ind < n_rows; ++ind) + { + res0 += shape_values[ind] * in[stride * ind]; + if (max_derivative > 0) + res1 += shape_values[ind + n_rows] * in[stride * ind]; + + if (max_derivative > 1) + res2 += shape_values[ind + 2 * n_rows] * in[stride * ind]; + } + if (add) + { + 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 (int col = 0; col < n_rows; ++col) + { + if (add) + out[col * stride] += shape_values[col] * in[0]; + else + out[col * stride] = shape_values[col] * in[0]; + + if (max_derivative > 0) + out[col * stride] += + shape_values[col + n_rows] * in[out_stride]; + + if (max_derivative > 1) + out[col * stride] += + shape_values[col + 2 * n_rows] * in[2 * out_stride]; + } + } + + // increment: in regular case, just go to the next point in + // x-direction. If we are at the end of one chunk in x-dir, need + // to jump over to the next layer in z-direction + switch (face_direction) + { + case 0: + 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 (normal_direction == 0) + { + if (contract_onto_face) + out += n_rows - 1; + else + in += n_rows - 1; + } + if (normal_direction == 1) + { + if (contract_onto_face) + out += n_rows - 2; + else + in += n_rows - 2; + } + if (normal_direction == 2) + { + if (contract_onto_face) + out += n_rows; + else + in += n_rows; + } + } + break; + + case 2: + ++in; + ++out; + break; + + default: + Assert(false, ExcNotImplemented()); + } + } + if (face_direction == 1 && dim == 3) + { + // adjust for local coordinate system zx + if (contract_onto_face) + { + if (normal_direction == 0) + { + in += (n_rows + 1) * (n_rows - 1); + out -= n_rows * (n_rows + 1) - 1; + } + if (normal_direction == 1) + { + in += (n_rows - 1) * (n_rows - 1); + out -= (n_rows - 1) * (n_rows - 1) - 1; + } + if (normal_direction == 2) + { + in += (n_rows - 1) * (n_rows); + out -= (n_rows) * (n_rows + 1) - 1; + } + } + else + { + if (normal_direction == 0) + { + out += (n_rows + 1) * (n_rows - 1); + in -= n_rows * (n_rows + 1) - 1; + } + if (normal_direction == 1) + { + out += (n_rows - 1) * (n_rows - 1); + in -= (n_rows - 1) * (n_rows - 1) - 1; + } + if (normal_direction == 2) + { + out += (n_rows - 1) * (n_rows); + in -= (n_rows) * (n_rows + 1) - 1; + } + } + } + } + } + + + /** * Struct to avoid using Tensor<1, dim, Point> in * evaluate_tensor_product_value_and_gradient because a Point cannot be used diff --git a/source/fe/fe_raviart_thomas_nodal.cc b/source/fe/fe_raviart_thomas_nodal.cc index 08ecfed20f..8397f71be1 100644 --- a/source/fe/fe_raviart_thomas_nodal.cc +++ b/source/fe/fe_raviart_thomas_nodal.cc @@ -147,67 +147,6 @@ namespace } - - // set up the numbering of the rt polynomials - std::vector - compute_rt_hierarchic_to_lexicographic(const unsigned int dim, - const unsigned int degree) - { - const unsigned int n_pols = - (degree + 2) * Utilities::pow(degree + 1, dim - 1); - - std::vector hierarchic_to_lexicographic; - - // dofs on faces - for (unsigned int face_no = 0; face_no < 2 * dim; ++face_no) - { - const unsigned int stride_x = face_no < 2 ? degree + 2 : 1; - const unsigned int stride_y = - face_no < 4 ? (degree + 2) * (degree + 1) : degree + 1; - const unsigned int offset = - (face_no % 2) * Utilities::pow(degree + 1, 1 + face_no / 2); - for (unsigned int j = 0; j < (dim > 2 ? degree + 1 : 1); ++j) - for (unsigned int i = 0; i < degree + 1; ++i) - hierarchic_to_lexicographic.push_back( - (face_no / 2) * n_pols + offset + i * stride_x + j * stride_y); - } - // dofs on cells, starting with x component... - for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k) - for (unsigned int j = 0; j < (dim > 1 ? degree + 1 : 1); ++j) - for (unsigned int i = 1; i < degree + 1; ++i) - hierarchic_to_lexicographic.push_back( - k * (degree + 1) * (degree + 2) + j * (degree + 2) + i); - // ... then y component ... - if (dim > 1) - for (unsigned int k = 0; k < (dim > 2 ? degree + 1 : 1); ++k) - for (unsigned int j = 1; j < degree + 1; ++j) - for (unsigned int i = 0; i < degree + 1; ++i) - hierarchic_to_lexicographic.push_back( - n_pols + k * (degree + 1) * (degree + 2) + j * (degree + 1) + i); - // ... and finally z component - if (dim > 2) - for (unsigned int k = 1; k < degree + 1; ++k) - for (unsigned int j = 0; j < degree + 1; ++j) - for (unsigned int i = 0; i < degree + 1; ++i) - hierarchic_to_lexicographic.push_back( - 2 * n_pols + k * (degree + 1) * (degree + 1) + j * (degree + 1) + - i); - - AssertDimension(hierarchic_to_lexicographic.size(), n_pols * dim); - -#ifdef DEBUG - // assert that we have a valid permutation - std::vector copy(hierarchic_to_lexicographic); - std::sort(copy.begin(), copy.end()); - for (unsigned int i = 0; i < copy.size(); ++i) - AssertDimension(i, copy[i]); -#endif - - return hierarchic_to_lexicographic; - } - - - template PolynomialsRaviartThomasNodal::PolynomialsRaviartThomasNodal( const unsigned int degree) @@ -219,11 +158,11 @@ namespace // actual order required by the finite element class with unknowns on // faces placed first const unsigned int n_pols = polynomial_space.n(); - hierarchic_to_lexicographic = - compute_rt_hierarchic_to_lexicographic(dim, degree); - lexicographic_to_hierarchic = - Utilities::invert_permutation(hierarchic_to_lexicographic); + internal::get_lexicographic_numbering_rt_nodal(degree + 1); + + hierarchic_to_lexicographic = + Utilities::invert_permutation(lexicographic_to_hierarchic); // since we only store an anisotropic polynomial for the first component, // we set up a second numbering to switch out the actual coordinate @@ -426,7 +365,61 @@ namespace } } // namespace +namespace internal +{ + template + std::vector + get_lexicographic_numbering_rt_nodal(const unsigned int degree) + { + const unsigned int n_dofs_face = Utilities::pow(degree, dim - 1); + std::vector lexicographic_numbering; + // component 1 + for (unsigned int j = 0; j < n_dofs_face; j++) + { + lexicographic_numbering.push_back(j); + for (unsigned int i = n_dofs_face * 2 * dim; + i < n_dofs_face * 2 * dim + degree - 1; + i++) + lexicographic_numbering.push_back(i + j * (degree - 1)); + lexicographic_numbering.push_back(n_dofs_face + j); + } + + // component 2 + unsigned int layers = (dim == 3) ? degree : 1; + for (unsigned int k = 0; k < layers; k++) + { + unsigned int k_add = k * degree; + for (unsigned int j = n_dofs_face * 2; j < n_dofs_face * 2 + degree; + j++) + lexicographic_numbering.push_back(j + k_add); + + for (unsigned int i = n_dofs_face * (2 * dim + (degree - 1)); + i < n_dofs_face * (2 * dim + (degree - 1)) + (degree - 1) * degree; + i++) + { + lexicographic_numbering.push_back(i + k_add * (degree - 1)); + } + for (unsigned int j = n_dofs_face * 3; j < n_dofs_face * 3 + degree; + j++) + lexicographic_numbering.push_back(j + k_add); + } + // component 3 + if (dim == 3) + { + for (unsigned int i = 4 * n_dofs_face; i < 5 * n_dofs_face; i++) + lexicographic_numbering.push_back(i); + for (unsigned int i = 6 * n_dofs_face + n_dofs_face * 2 * (degree - 1); + i < 6 * n_dofs_face + n_dofs_face * 3 * (degree - 1); + i++) + lexicographic_numbering.push_back(i); + for (unsigned int i = 5 * n_dofs_face; i < 6 * n_dofs_face; i++) + lexicographic_numbering.push_back(i); + } + + return lexicographic_numbering; + } +} // namespace internal // --------------------- actual implementation of element -------------------- diff --git a/source/matrix_free/CMakeLists.txt b/source/matrix_free/CMakeLists.txt index bdaa5a9220..ecc6b35ab3 100644 --- a/source/matrix_free/CMakeLists.txt +++ b/source/matrix_free/CMakeLists.txt @@ -23,7 +23,17 @@ SET(_src evaluation_template_factory_inst4.cc evaluation_template_factory_inst5.cc evaluation_template_factory_inst6.cc + evaluation_template_factory_inst7.cc + evaluation_template_factory_inst8.cc + evaluation_template_factory_inst9.cc + evaluation_template_factory_inst10.cc evaluation_template_factory_hanging_nodes.cc + evaluation_template_face_factory.cc + evaluation_template_face_factory_inst2.cc + evaluation_template_face_factory_inst3.cc + evaluation_template_face_factory_inst4.cc + evaluation_template_face_factory_inst5.cc + evaluation_template_face_factory_inst6.cc fe_point_evaluation.cc mapping_info.cc mapping_info_inst2.cc @@ -36,6 +46,7 @@ SET(_src SET(_inst evaluation_template_factory.inst.in + evaluation_template_face_factory.inst.in evaluation_template_factory_hanging_nodes.inst.in fe_point_evaluation.inst.in mapping_info.inst.in diff --git a/source/matrix_free/evaluation_template_face_factory.cc b/source/matrix_free/evaluation_template_face_factory.cc new file mode 100644 index 0000000000..9a2df3aeb5 --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory.cc @@ -0,0 +1,28 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + +#define SPLIT_INSTANTIATIONS_COUNT 6 +#ifndef SPLIT_INSTANTIATIONS_INDEX +# define SPLIT_INSTANTIATIONS_INDEX 0 +#endif + +#include "evaluation_template_face_factory.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/matrix_free/evaluation_template_face_factory.inst.in b/source/matrix_free/evaluation_template_face_factory.inst.in new file mode 100644 index 0000000000..d98cac1e9f --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory.inst.in @@ -0,0 +1,32 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; + deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED) + { + template struct dealii::internal:: + FEFaceEvaluationFactory; + + template struct dealii::internal::FEFaceEvaluationGatherFactory< + deal_II_dimension, + double, + deal_II_scalar_vectorized>; + + template struct dealii::internal::FEFaceEvaluationGatherFactory< + deal_II_dimension, + float, + deal_II_scalar_vectorized>; + } diff --git a/source/matrix_free/evaluation_template_face_factory_inst2.cc b/source/matrix_free/evaluation_template_face_factory_inst2.cc new file mode 100644 index 0000000000..5b3b84081d --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory_inst2.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 1 +#include "evaluation_template_face_factory.cc" diff --git a/source/matrix_free/evaluation_template_face_factory_inst3.cc b/source/matrix_free/evaluation_template_face_factory_inst3.cc new file mode 100644 index 0000000000..979c7e432b --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory_inst3.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 2 +#include "evaluation_template_face_factory.cc" diff --git a/source/matrix_free/evaluation_template_face_factory_inst4.cc b/source/matrix_free/evaluation_template_face_factory_inst4.cc new file mode 100644 index 0000000000..b444b1107f --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory_inst4.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 3 +#include "evaluation_template_face_factory.cc" diff --git a/source/matrix_free/evaluation_template_face_factory_inst5.cc b/source/matrix_free/evaluation_template_face_factory_inst5.cc new file mode 100644 index 0000000000..bcfd4482bb --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory_inst5.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 4 +#include "evaluation_template_face_factory.cc" diff --git a/source/matrix_free/evaluation_template_face_factory_inst6.cc b/source/matrix_free/evaluation_template_face_factory_inst6.cc new file mode 100644 index 0000000000..8aaeedec00 --- /dev/null +++ b/source/matrix_free/evaluation_template_face_factory_inst6.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 5 +#include "evaluation_template_face_factory.cc" diff --git a/source/matrix_free/evaluation_template_factory.cc b/source/matrix_free/evaluation_template_factory.cc index cab955b323..9d2bf69319 100644 --- a/source/matrix_free/evaluation_template_factory.cc +++ b/source/matrix_free/evaluation_template_factory.cc @@ -18,7 +18,7 @@ DEAL_II_NAMESPACE_OPEN -#define SPLIT_INSTANTIATIONS_COUNT 6 +#define SPLIT_INSTANTIATIONS_COUNT 10 #ifndef SPLIT_INSTANTIATIONS_INDEX # define SPLIT_INSTANTIATIONS_INDEX 0 #endif diff --git a/source/matrix_free/evaluation_template_factory.inst.in b/source/matrix_free/evaluation_template_factory.inst.in index 03b1de0a40..441855d054 100644 --- a/source/matrix_free/evaluation_template_factory.inst.in +++ b/source/matrix_free/evaluation_template_factory.inst.in @@ -20,19 +20,6 @@ for (deal_II_dimension : DIMENSIONS; template struct dealii::internal:: FEEvaluationFactory; - template struct dealii::internal:: - FEFaceEvaluationFactory; - - template struct dealii::internal::FEFaceEvaluationGatherFactory< - deal_II_dimension, - double, - deal_II_scalar_vectorized>; - - template struct dealii::internal::FEFaceEvaluationGatherFactory< - deal_II_dimension, - float, - deal_II_scalar_vectorized>; - // inverse mass template struct dealii::internal:: CellwiseInverseMassFactory; diff --git a/source/matrix_free/evaluation_template_factory_inst10.cc b/source/matrix_free/evaluation_template_factory_inst10.cc new file mode 100644 index 0000000000..b93c12d5a0 --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_inst10.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 9 +#include "evaluation_template_factory.cc" diff --git a/source/matrix_free/evaluation_template_factory_inst7.cc b/source/matrix_free/evaluation_template_factory_inst7.cc new file mode 100644 index 0000000000..d71730bdc3 --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_inst7.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 6 +#include "evaluation_template_factory.cc" diff --git a/source/matrix_free/evaluation_template_factory_inst8.cc b/source/matrix_free/evaluation_template_factory_inst8.cc new file mode 100644 index 0000000000..47f1f7a3fe --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_inst8.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 7 +#include "evaluation_template_factory.cc" diff --git a/source/matrix_free/evaluation_template_factory_inst9.cc b/source/matrix_free/evaluation_template_factory_inst9.cc new file mode 100644 index 0000000000..94dc699ee6 --- /dev/null +++ b/source/matrix_free/evaluation_template_factory_inst9.cc @@ -0,0 +1,17 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#define SPLIT_INSTANTIATIONS_INDEX 8 +#include "evaluation_template_factory.cc" diff --git a/tests/matrix_free/matrix_vector_rt_01.cc b/tests/matrix_free/matrix_vector_rt_01.cc new file mode 100644 index 0000000000..d358717dc3 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_01.cc @@ -0,0 +1,44 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// This function tests the correctness of the matrix-free implementation +// of the FE_RaviartThomasNodal element by evaluating a simple fe operator +// and comparing the result with FEVaules which is considered the +// reference. The mesh is a hypercube mesh with no hanging nodes and no other +// constraints + +#include "../tests.h" + +#include "matrix_vector_rt_common.h" + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(2); + + FE_RaviartThomasNodal fe(fe_degree - 1); + + DoFHandler dof(tria); + + dof.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + do_test(dof, constraints); +} diff --git a/tests/matrix_free/matrix_vector_rt_01.output b/tests/matrix_free/matrix_vector_rt_01.output new file mode 100644 index 0000000000..0e53c41bb8 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_01.output @@ -0,0 +1,25 @@ + +DEAL:2d::Testing FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 16 +DEAL:2d::Number of degrees of freedom: 144 +DEAL:2d:: +DEAL:2d::Norm of difference: 6.20418e-16 +DEAL:2d:: +DEAL:2d::Testing FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 16 +DEAL:2d::Number of degrees of freedom: 312 +DEAL:2d:: +DEAL:2d::Norm of difference: 7.03559e-16 +DEAL:2d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 64 +DEAL:3d::Number of degrees of freedom: 1728 +DEAL:3d:: +DEAL:3d::Norm of difference: 1.04798e-15 +DEAL:3d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(2) +DEAL:3d::Number of cells: 64 +DEAL:3d::Number of degrees of freedom: 5616 +DEAL:3d:: +DEAL:3d::Norm of difference: 1.49180e-15 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_common.h b/tests/matrix_free/matrix_vector_rt_common.h new file mode 100644 index 0000000000..e2d0ae2324 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_common.h @@ -0,0 +1,241 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// This test template evaluates a simple operator on FE_PolyTensor +// elements using FEEvaluation and compares the result with the output +// of FEValues (which is considered to be the reference) on cartesian +// meshes without hanging nodes. (It will be extended to also handle general +// meshes and hanging nodes in the future.) The test do not include +// multithreading because FEValues is not thread-safe. +// See matrix_vector_rt_01.cc for an example that uses this template. + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +// forward declare this function. will be implemented in .cc files +template +void +test(); + + + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + virtual ~MatrixFreeTest(){}; + + void + operator()(const MatrixFree & data, + Vector & dst, + const Vector & src, + const std::pair &cell_range) const + { + FEEvaluation fe_eval(data); + + // OBS! This will need to be modified once the Piola transform is + // implemented + unsigned int n_cells = + data.get_dof_handler().get_triangulation().n_active_cells(); + Number piola = + (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4); + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + fe_eval.submit_value(Number(10 * piola) * fe_eval.get_value(q), q); + fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q); + } + + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }; + + + void + test_functions(Vector &dst, const Vector &src) const + { + data.cell_loop(&MatrixFreeTest::operator(), this, dst, src); + }; + +protected: + const MatrixFree &data; +}; + + +template +void +do_test(const DoFHandler & dof, + const AffineConstraints &constraints) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + deallog << "Number of cells: " << dof.get_triangulation().n_active_cells() + << std::endl; + deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl + << std::endl; + + + // constraints.distribute(solution); + MatrixFree mf_data; + { + const QGaussLobatto<1> quad(fe_degree + 2); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_hessians; + mf_data.reinit(dof, constraints, quad, data); + } + + // create vector with random entries + Vector solution, initial_condition; + + mf_data.initialize_dof_vector(solution); + mf_data.initialize_dof_vector(initial_condition); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + initial_condition[i] = random_value(); + } + + // MatrixFree solution + MatrixFreeTest mf(mf_data); + mf.test_functions(solution, initial_condition); + + + SparsityPattern sp; + SparseMatrix system_matrix; + DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern(dof, dsp); + sp.copy_from(dsp); + system_matrix.reinit(sp); + + FEValues fe_val(dof.get_fe(), + Quadrature(mf_data.get_quadrature(0)), + update_values | update_gradients | update_JxW_values); + + + const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell; + std::vector local_dof_indices(dofs_per_cell); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + + const FEValuesExtractors::Vector velocities(0); + + // Assemble matrix + for (const auto &cell : dof.active_cell_iterators()) + { + fe_val.reinit(cell); + local_matrix = 0; + + for (const auto q : fe_val.quadrature_point_indices()) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.; + const Tensor<2, dim> grad_phi_i = fe_val[velocities].gradient(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q); + const Tensor<2, dim> grad_phi_j = + fe_val[velocities].gradient(j, q); + + local_matrix(i, j) += + (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) * + fe_val.JxW(q); + } + } + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + local_matrix(i, j)); + } + + Vector ref(solution.size()); + + // Compute reference + system_matrix.vmult(ref, initial_condition); + + ref -= solution; + + const double diff_norm = ref.linfty_norm() / solution.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int +main() +{ + initlog(); + + { + deallog.push("2d"); + test<2, 2>(); + test<2, 3>(); + deallog.pop(); + deallog.push("3d"); + test<3, 2>(); + test<3, 3>(); + deallog.pop(); + } +} diff --git a/tests/matrix_free/matrix_vector_rt_face_01.cc b/tests/matrix_free/matrix_vector_rt_face_01.cc new file mode 100644 index 0000000000..efb1e82ee6 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_01.cc @@ -0,0 +1,44 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// This function tests the correctness of the matrix-free implementation +// of the FE_RaviartThomasNodal element by evaluating a face operator +// and comparing the result with FEVaules which is considered the +// reference. The mesh is a hypercube mesh with no hanging nodes and no other +// constraints + +#include "../tests.h" + +#include "matrix_vector_rt_face_common.h" + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + + FE_RaviartThomasNodal fe(fe_degree - 1); + + DoFHandler dof(tria); + + dof.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + do_test(dof, constraints); +} diff --git a/tests/matrix_free/matrix_vector_rt_face_01.output b/tests/matrix_free/matrix_vector_rt_face_01.output new file mode 100644 index 0000000000..99a9ae91f1 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_01.output @@ -0,0 +1,19 @@ + +DEAL:2d::Testing FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 40 +DEAL:2d:: +DEAL:2d::Norm of difference: 3.44596e-16 +DEAL:2d:: +DEAL:2d::Testing FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 84 +DEAL:2d:: +DEAL:2d::Norm of difference: 3.47432e-15 +DEAL:2d:: +DEAL:3d::Testing FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 8 +DEAL:3d::Number of degrees of freedom: 240 +DEAL:3d:: +DEAL:3d::Norm of difference: 2.28531e-15 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_face_common.h b/tests/matrix_free/matrix_vector_rt_face_common.h new file mode 100644 index 0000000000..922ab008dc --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_common.h @@ -0,0 +1,311 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// This test template evaluates a simple operator on FE_PolyTensor +// elements using FEFaceEvaluation and compares the result with the output +// of FEFaceValues (which is considered to be the reference) on cartesian +// meshes without hanging nodes. (It will be extended to also handle general +// meshes and hanging nodes in the future.) The test do not include +// multithreading because FEFaceValues is not thread-safe. +// See matrix_vector_rt_face_01.cc for an example that uses this template. + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "../tests.h" + + +// forward declare this function. will be implemented in .cc files +template +void +test(); + + + +template +class MatrixFreeTest +{ +public: + MatrixFreeTest(const MatrixFree &data_in) + : data(data_in){}; + + virtual ~MatrixFreeTest(){}; + + void + dummy(const MatrixFree &, + Vector &, + const Vector &, + const std::pair &) const + {} + + void + operator_face(const MatrixFree & data, + Vector & dst, + const Vector & src, + const std::pair &face_range) const + { + FEFaceEvaluation fe_eval(data, + true); + FEFaceEvaluation fe_eval_n( + data, false); + + // Note that this will need to be modified once the Piola transform is + // implemented + const unsigned int n_cells = + data.get_dof_handler().get_triangulation().n_active_cells(); + const Number piola = + (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4); + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_eval.reinit(face); + fe_eval_n.reinit(face); + + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + fe_eval_n.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q); + fe_eval_n.submit_value(Number(10. * piola) * fe_eval_n.get_value(q), + q); + + fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q); + fe_eval_n.submit_gradient(Number(piola) * fe_eval_n.get_gradient(q), + q); + } + + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + fe_eval_n.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }; + + void + operator_boundary( + const MatrixFree & data, + Vector & dst, + const Vector & src, + const std::pair &face_range) const + { + FEFaceEvaluation fe_eval(data, + true); + + // Note that this will need to be modified once the Piola transform is + // implemented + const unsigned int n_cells = + data.get_dof_handler().get_triangulation().n_active_cells(); + const Number piola = + (dim == 2) ? n_cells : Utilities::pow((int)std::cbrt(n_cells), 4); + + for (unsigned int face = face_range.first; face < face_range.second; ++face) + { + fe_eval.reinit(face); + fe_eval.gather_evaluate(src, + EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < fe_eval.n_q_points; ++q) + { + fe_eval.submit_value(Number(10. * piola) * fe_eval.get_value(q), q); + fe_eval.submit_gradient(Number(piola) * fe_eval.get_gradient(q), q); + } + + fe_eval.integrate_scatter(EvaluationFlags::values | + EvaluationFlags::gradients, + dst); + } + }; + + + void + test_functions(Vector &dst, const Vector &src) const + { + data.loop(&MatrixFreeTest::dummy, + &MatrixFreeTest::operator_face, + &MatrixFreeTest::operator_boundary, + this, + dst, + src); + }; + +protected: + const MatrixFree &data; +}; + + + +template +void +do_test(const DoFHandler & dof, + const AffineConstraints &constraints) +{ + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + deallog << "Number of cells: " << dof.get_triangulation().n_active_cells() + << std::endl; + deallog << "Number of degrees of freedom: " << dof.n_dofs() << std::endl + << std::endl; + + + // constraints.distribute(solution); + MatrixFree mf_data; + { + const QGaussLobatto<1> quad(fe_degree + 2); + const MappingQ mapping(fe_degree); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients | update_JxW_values; + data.mapping_update_flags_inner_faces = + (update_gradients | update_JxW_values); + data.mapping_update_flags_boundary_faces = + (update_gradients | update_JxW_values); + mf_data.reinit(mapping, dof, constraints, quad, data); + } + + Vector solution, initial_condition; + + // create vector with random entries + mf_data.initialize_dof_vector(solution); + mf_data.initialize_dof_vector(initial_condition); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i)) + continue; + initial_condition[i] = random_value(); + } + + MatrixFreeTest mf(mf_data); + mf.test_functions(solution, initial_condition); + + + SparsityPattern sp; + SparseMatrix system_matrix; + DynamicSparsityPattern dsp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern(dof, dsp); + sp.copy_from(dsp); + system_matrix.reinit(sp); + + FEFaceValues fe_val( + dof.get_fe(), + QGaussLobatto(mf_data.get_quadrature(0).size()), + update_values | update_gradients | update_JxW_values); + + + const unsigned int dofs_per_cell = fe_val.get_fe().dofs_per_cell; + std::vector local_dof_indices(dofs_per_cell); + FullMatrix local_matrix(dofs_per_cell, dofs_per_cell); + + + const FEValuesExtractors::Vector velocities(0); + // Assemble matrix + for (const auto &cell : dof.active_cell_iterators()) + for (const auto &face : cell->face_iterators()) + { + fe_val.reinit(cell, face); + local_matrix = 0; + + for (const auto q : fe_val.quadrature_point_indices()) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const Tensor<1, dim> phi_i = fe_val[velocities].value(i, q) * 10.; + const Tensor<2, dim> grad_phi_i = + fe_val[velocities].gradient(i, q); + + for (unsigned int j = 0; j < dofs_per_cell; ++j) + { + const Tensor<1, dim> phi_j = fe_val[velocities].value(j, q); + const Tensor<2, dim> grad_phi_j = + fe_val[velocities].gradient(j, q); + local_matrix(i, j) += + (phi_j * phi_i + scalar_product(grad_phi_i, grad_phi_j)) * + fe_val.JxW(q); + } + } + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add(local_dof_indices[i], + local_dof_indices[j], + local_matrix(i, j)); + } + + Vector ref(solution.size()); + + // Compute reference + system_matrix.vmult(ref, initial_condition); + + ref -= solution; + + const double diff_norm = ref.linfty_norm() / solution.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int +main() +{ + initlog(); + + { + deallog.push("2d"); + test<2, 2>(); + test<2, 3>(); + deallog.pop(); + deallog.push("3d"); + test<3, 2>(); + // test<3, 3>(); //Takes way too long + deallog.pop(); + } +} -- 2.39.5