From: Niklas Wik Date: Mon, 2 May 2022 08:55:47 +0000 (+0200) Subject: Piola transform for affine cells X-Git-Tag: v9.4.0-rc1~170^2~4 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5ac08a75fe14ab11c12a3f06e324dd9cff8a58c8;p=dealii.git Piola transform for affine cells Updated tests, reordering of components in face evaluation, and storing jacobian for affine face evaluation --- diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index 90c53097c9..00c9d4c7e5 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3055,10 +3055,13 @@ namespace internal template static EvalType create_evaluator_tensor_product( - const MatrixFreeFunctions::UnivariateShapeData &data, - const unsigned int subface_index, - const unsigned int direction) + const MatrixFreeFunctions::ShapeInfo &shape_info, + const unsigned int subface_index, + const unsigned int direction) { + const MatrixFreeFunctions::UnivariateShapeData &data = + (std::is_same::value) ? shape_info.data.front() : + shape_info.data.back(); if (subface_index >= GeometryInfo::max_children_per_cell) return EvalType(data.shape_values, data.shape_gradients, @@ -3078,199 +3081,107 @@ namespace internal 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) + const EvaluationFlags::EvaluationFlags evaluation_flag, + 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<0, EvalNormal, EvalTangent>( + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + std::integral_constant()); - 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, EvalTangent, EvalNormal>( + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + std::integral_constant()); - 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()); - } - } + evaluate_in_face_apply<2, EvalGeneral, EvalGeneral>( + values_dofs, + fe_eval, + scratch_data, + evaluation_flag, + face_direction, + subface_index, + 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. + * evaluation path, and std::integral_constant below is the + * integration path. These two functions can be fused together since all + * offsets and pointers are the exact same. */ - template + 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, + const unsigned int face_direction, + const unsigned int subface_index, std::integral_constant) { - constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + // TODO. Need to make sure hanging nodes work, i.e call + // create_eval_tensor_product with correct direction. + Eval0 eval0 = + create_evaluator_tensor_product(fe_eval.get_shape_info(), + subface_index, + 0); + Eval1 eval1 = + create_evaluator_tensor_product(fe_eval.get_shape_info(), + subface_index, + 0); - Number *values_quad = fe_eval.begin_values() + n_q_points * component; + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + const std::size_t n_dofs_tangent = + fe_eval.get_shape_info().dofs_per_component_on_face; + const std::size_t n_dofs_normal = + n_dofs_tangent - Utilities::pow(fe_degree, dim - 2); + const std::size_t dofs_stride = + (std::is_same::value) ? n_dofs_normal : + n_dofs_tangent; + + const unsigned int component_table[3][3] = {{1, 2, 0}, + {2, 0, 1}, + {0, 1, 2}}; + const unsigned int component = + (dim == 2 && normal_dir == 0 && face_direction == 1) ? + 0 : + component_table[face_direction][normal_dir]; + + // Initial offsets + values_dofs += + 3 * ((component == 0) ? + 0 : + ((component == 1) ? + ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) : + ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent : + n_dofs_normal + n_dofs_tangent))); + const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir; + Number *values_quad = fe_eval.begin_values() + n_q_points * shift; Number *gradients_quad = - fe_eval.begin_gradients() + dim * n_q_points * component; + fe_eval.begin_gradients() + dim * n_q_points * shift; Number *hessians_quad = - fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component; + fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift; // Evaluation path - if ((evaluation_flag & EvaluationFlags::values) && !(evaluation_flag & EvaluationFlags::gradients)) { @@ -3393,25 +3304,59 @@ namespace internal } } - template + 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, + const unsigned int face_direction, + const unsigned int subface_index, std::integral_constant) { - constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + // TODO. Need to make sure hanging nodes work, i.e call + // create_eval_tensor_product with correct direction. + Eval0 eval0 = + create_evaluator_tensor_product(fe_eval.get_shape_info(), + subface_index, + 0); + Eval1 eval1 = + create_evaluator_tensor_product(fe_eval.get_shape_info(), + subface_index, + 0); - Number *values_quad = fe_eval.begin_values() + n_q_points * component; + constexpr std::size_t n_q_points = Utilities::pow(n_q_points_1d, dim - 1); + const std::size_t n_dofs_tangent = + fe_eval.get_shape_info().dofs_per_component_on_face; + const std::size_t n_dofs_normal = + n_dofs_tangent - Utilities::pow(fe_degree, dim - 2); + const std::size_t dofs_stride = + (std::is_same::value) ? n_dofs_normal : + n_dofs_tangent; + + const unsigned int component_table[3][3] = {{1, 2, 0}, + {2, 0, 1}, + {0, 1, 2}}; + const unsigned int component = + (dim == 2 && normal_dir == 0 && face_direction == 1) ? + 0 : + component_table[face_direction][normal_dir]; + + // Initial offsets + values_dofs += + 3 * ((component == 0) ? + 0 : + ((component == 1) ? + ((face_direction == 0) ? n_dofs_normal : n_dofs_tangent) : + ((face_direction == 2) ? n_dofs_tangent + n_dofs_tangent : + n_dofs_normal + n_dofs_tangent))); + const unsigned int shift = (dim == 2) ? normal_dir / 2 : normal_dir; + Number *values_quad = fe_eval.begin_values() + n_q_points * shift; Number *gradients_quad = - fe_eval.begin_gradients() + dim * n_q_points * component; + fe_eval.begin_gradients() + dim * n_q_points * shift; Number *hessians_quad = - fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * component; + fe_eval.begin_hessians() + dim * (dim + 1) / 2 * n_q_points * shift; // Integration path if ((evaluation_flag & EvaluationFlags::values) && @@ -4265,14 +4210,12 @@ namespace internal (n_q_points_1d < 1) ? 1 : n_q_points_1d, Number>:: - template evaluate_or_integrate_in_face( - evaluation_flag, - shape_info, - temp, - fe_eval, - scratch_data, - subface_index, - fe_eval.get_face_no()); + template evaluate_or_integrate_in_face(evaluation_flag, + temp, + fe_eval, + scratch_data, + subface_index, + fe_eval.get_face_no()); } else if (fe_degree > -1 && subface_index >= GeometryInfo::max_children_per_cell && @@ -4514,7 +4457,6 @@ namespace internal n_q_points_1d, Number>:: template evaluate_or_integrate_in_face(integration_flag, - shape_info, temp, fe_eval, scratch_data, diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 8102c5dbfb..516b6cb4c5 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -1005,6 +1005,12 @@ public: using BaseClass = FEEvaluationBase; + /** + * @copydoc FEEvaluationBase::get_value() + */ + value_type + get_value(const unsigned int q_point) const; + /** * @copydoc FEEvaluationBase::get_gradient() */ @@ -1046,6 +1052,13 @@ public: gradient_type get_hessian_diagonal(const unsigned int q_point) const; + /** + * @copydoc FEEvaluationBase::submit_value() + */ + void + submit_value(const Tensor<1, dim, VectorizedArrayType> val_in, + const unsigned int q_point); + /** * @copydoc FEEvaluationBase::submit_gradient() */ @@ -5707,13 +5720,152 @@ FEEvaluationAccess::operator=( } +template +inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType> +FEEvaluationAccess::get_value( + const unsigned int q_point) const +{ + // Check if Piola transform is required + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { +# ifdef DEBUG + Assert(this->values_quad_initialized == true, + internal::ExcAccessToUninitializedField()); +# endif + + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_values")); + const std::size_t nqp = this->n_quadrature_points; + Tensor<1, dim, VectorizedArrayType> value_out; + + // Cartesian cell + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, dealii::VectorizedArray> jac = + this->jacobian[1]; + const VectorizedArrayType inv_det = determinant(this->jacobian[0]); + + for (unsigned int comp = 0; comp < n_components; ++comp) + value_out[comp] = this->values_quad[comp * nqp + q_point] * + jac[comp][comp] * + inv_det; // / this->jacobian[0][comp][comp]; + } + + // Affine or general cell + else + { + const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->jacobian[q_point] : + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + transpose(invert(inv_t_jac)) : + this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac); + // J * u * det(J^-1) + for (unsigned int comp = 0; comp < n_components; ++comp) + { + value_out[comp] = + this->values_quad[q_point] * jac[comp][0] * inv_det; + for (unsigned int e = 1; e < dim; ++e) + value_out[comp] += + this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det; + } + } + return value_out; + } + else + { + return BaseClass::get_value(q_point); + } +} template inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> FEEvaluationAccess:: get_gradient(const unsigned int q_point) const { - return BaseClass::get_gradient(q_point); + // Check if Piola transform is required + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { +# ifdef DEBUG + Assert(this->gradients_quad_initialized == true, + internal::ExcAccessToUninitializedField()); +# endif + + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + const std::size_t nqp = this->n_quadrature_points; + Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_out; + + // Cartesian cell + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + const VectorizedArrayType inv_det = determinant(inv_t_jac); + + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int comp = 0; comp < n_components; ++comp) + grad_out[comp][d] = + this->gradients_quad[(comp * dim + d) * nqp + q_point] * + inv_t_jac[d][d] * jac[comp][comp] * inv_det; + } + // Affine cell + else if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac); + + VectorizedArrayType tmp; + // J * grad_quad * J^-1 * det(J^-1) + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + tmp = 0; + for (unsigned int f = 0; f < dim; ++f) + for (unsigned int e = 0; e < dim; ++e) + tmp += jac[comp][f] * inv_t_jac[d][e] * inv_det * + this->gradients_quad[(f * dim + e) * nqp + q_point]; + + grad_out[comp][d] = tmp; + } + } + // General cell TODO + else + { + // Here we need the jacobian gradient and not the inverse which is + // stored in this->jacobian_gradients + AssertThrow(false, ExcNotImplemented()); + } + return grad_out; + } + else + { + return BaseClass::get_gradient(q_point); + } } @@ -5735,28 +5887,55 @@ FEEvaluationAccess:: VectorizedArrayType divergence; const std::size_t nqp = this->n_quadrature_points; - // Cartesian cell - if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) { - divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0]; - for (unsigned int d = 1; d < dim; ++d) - divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] * - this->jacobian[0][d][d]; + // Affine cell + if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + // Derivatives are reordered for faces. Need to take this into account + const VectorizedArrayType inv_det = + (is_face && dim == 2 && this->get_face_no() < 2) ? + -determinant(this->jacobian[0]) : + determinant(this->jacobian[0]); + + divergence = this->gradients_quad[q_point] * inv_det; + for (unsigned int d = 1; d < dim; ++d) + divergence += + this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det; + } + // General cell TODO + else + { + Assert(false, ExcNotImplemented()); + } } - // cell with general/constant Jacobian else { - const Tensor<2, dim, VectorizedArrayType> &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : - this->jacobian[0]; - divergence = jac[0][0] * this->gradients_quad[q_point]; - for (unsigned int e = 1; e < dim; ++e) - divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point]; - for (unsigned int d = 1; d < dim; ++d) - for (unsigned int e = 0; e < dim; ++e) - divergence += - jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point]; + // Cartesian cell + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + divergence = this->gradients_quad[q_point] * this->jacobian[0][0][0]; + for (unsigned int d = 1; d < dim; ++d) + divergence += this->gradients_quad[(dim * d + d) * nqp + q_point] * + this->jacobian[0][d][d]; + } + // cell with general/constant Jacobian + else + { + const Tensor<2, dim, VectorizedArrayType> &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : + this->jacobian[0]; + divergence = jac[0][0] * this->gradients_quad[q_point]; + for (unsigned int e = 1; e < dim; ++e) + divergence += jac[0][e] * this->gradients_quad[e * nqp + q_point]; + for (unsigned int d = 1; d < dim; ++d) + for (unsigned int e = 0; e < dim; ++e) + divergence += + jac[d][e] * this->gradients_quad[(d * dim + e) * nqp + q_point]; + } } return divergence; } @@ -5854,6 +6033,79 @@ FEEvaluationAccess::get_hessian( } +template +inline DEAL_II_ALWAYS_INLINE void +FEEvaluationAccess:: + submit_value(const Tensor<1, dim, VectorizedArrayType> val_in, + const unsigned int q_point) +{ + // Check if Piola transform is required + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { + AssertIndexRange(q_point, this->n_quadrature_points); + + // This is not needed, but might be good to check anyway? + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_value")); +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); + this->values_quad_submitted = true; +# endif + + const std::size_t nqp = this->n_quadrature_points; + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, dealii::VectorizedArray> jac = + this->jacobian[1]; + const VectorizedArrayType weight = this->quadrature_weights[q_point]; + + for (unsigned int comp = 0; comp < n_components; ++comp) + this->values_quad[comp * nqp + q_point] = + val_in[comp] * weight * jac[comp][comp]; + } + // Affine or general cell + else + { + const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->jacobian[q_point] : + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = + (this->cell_type > internal::MatrixFreeFunctions::affine) ? + invert(inv_t_jac) : + this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into account + // and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? + this->quadrature_weights[q_point] : + (((this->cell_type > internal::MatrixFreeFunctions::affine) ? + this->J_value[q_point] : + this->J_value[0] * this->quadrature_weights[q_point]) * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac))); + + // J^T * u * w + for (unsigned int comp = 0; comp < n_components; ++comp) + { + this->values_quad[comp * nqp + q_point] = + val_in[0] * jac[0][comp] * fac; + for (unsigned int e = 1; e < dim; ++e) + this->values_quad[comp * nqp + q_point] += + val_in[e] * jac[e][comp] * fac; + } + } + } + else + { + BaseClass::submit_value(val_in, q_point); + } +} template inline DEAL_II_ALWAYS_INLINE void @@ -5861,7 +6113,76 @@ FEEvaluationAccess:: submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in, const unsigned int q_point) { - BaseClass::submit_gradient(grad_in, q_point); + // Check if Piola transform is required + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); +# endif + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); +# ifdef DEBUG + this->gradients_quad_submitted = true; +# endif + + const std::size_t nqp = this->n_quadrature_points; + // Cartesian cell + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + const VectorizedArrayType weight = this->quadrature_weights[q_point]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int comp = 0; comp < n_components; ++comp) + this->gradients_quad[(comp * dim + d) * nqp + q_point] = + grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight; + } + // Affine cell + else if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into account + // and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? this->quadrature_weights[q_point] : + this->J_value[0] * this->quadrature_weights[q_point] * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac)); + + // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType tmp = 0; + for (unsigned int f = 0; f < dim; ++f) + for (unsigned int e = 0; e < dim; ++e) + tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac; + + this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp; + } + } + // General cell TODO + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + else + { + BaseClass::submit_gradient(grad_in, q_point); + } } @@ -5873,7 +6194,76 @@ FEEvaluationAccess:: const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in, const unsigned int q_point) { - BaseClass::submit_gradient(grad_in, q_point); + // Check if Piola transform is required + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) + { +# ifdef DEBUG + Assert(this->is_reinitialized, ExcNotInitialized()); +# endif + AssertIndexRange(q_point, this->n_quadrature_points); + Assert(this->J_value != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); + Assert(this->jacobian != nullptr, + internal::ExcMatrixFreeAccessToUninitializedMappingField( + "update_gradients")); +# ifdef DEBUG + this->gradients_quad_submitted = true; +# endif + + const std::size_t nqp = this->n_quadrature_points; + // Cartesian cell + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) + { + const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + const VectorizedArrayType weight = this->quadrature_weights[q_point]; + for (unsigned int d = 0; d < dim; ++d) + for (unsigned int comp = 0; comp < n_components; ++comp) + this->gradients_quad[(comp * dim + d) * nqp + q_point] = + grad_in[comp][d] * inv_t_jac[d][d] * jac[comp][comp] * weight; + } + // Affine cell + else if (this->cell_type <= internal::MatrixFreeFunctions::affine) + { + const Tensor<2, dim, dealii::VectorizedArray> &inv_t_jac = + this->jacobian[0]; + const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1]; + + // Derivatives are reordered for faces. Need to take this into account + // and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + (!is_face) ? this->quadrature_weights[q_point] : + this->J_value[0] * this->quadrature_weights[q_point] * + ((dim == 2 && this->get_face_no() < 2) ? + -determinant(inv_t_jac) : + determinant(inv_t_jac)); + + // J_{j,i} * J^{-1}_{k,m} * grad_in_{j,m} * factor + for (unsigned int comp = 0; comp < n_components; ++comp) + for (unsigned int d = 0; d < dim; ++d) + { + VectorizedArrayType tmp = 0; + for (unsigned int f = 0; f < dim; ++f) + for (unsigned int e = 0; e < dim; ++e) + tmp += jac[f][comp] * inv_t_jac[e][d] * grad_in[f][e] * fac; + + this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp; + } + } + // General cell TODO + else + { + AssertThrow(false, ExcNotImplemented()); + } + } + else + { + BaseClass::submit_gradient(grad_in, q_point); + } } @@ -5899,39 +6289,77 @@ FEEvaluationAccess:: # endif const std::size_t nqp = this->n_quadrature_points; - if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian) + if (this->data->element_type == + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) { - const VectorizedArrayType fac = - this->J_value[0] * this->quadrature_weights[q_point] * div_in; - for (unsigned int d = 0; d < dim; ++d) + // Affine cell + if (this->cell_type <= internal::MatrixFreeFunctions::affine) { - this->gradients_quad[(d * dim + d) * nqp + q_point] = - (fac * this->jacobian[0][d][d]); - for (unsigned int e = d + 1; e < dim; ++e) + // Derivatives are reordered for faces. Need to take this into account + // and 1/inv_det != J_value for faces + const VectorizedArrayType fac = + ((!is_face) ? + 1 : + this->J_value[0] * ((dim == 2 && this->get_face_no() < 2) ? + -determinant(this->jacobian[0]) : + determinant(this->jacobian[0]))) * + this->quadrature_weights[q_point] * div_in; + + for (unsigned int d = 0; d < dim; ++d) { - this->gradients_quad[(d * dim + e) * nqp + q_point] = - VectorizedArrayType(); - this->gradients_quad[(e * dim + d) * nqp + q_point] = - VectorizedArrayType(); + this->gradients_quad[(dim * d + d) * nqp + q_point] = fac; + for (unsigned int e = d + 1; e < dim; ++e) + { + this->gradients_quad[(dim * d + e) * nqp + q_point] = + VectorizedArrayType(); + this->gradients_quad[(dim * e + d) * nqp + q_point] = + VectorizedArrayType(); + } } } + // General cell TODO + else + { + AssertThrow(false, ExcNotImplemented()); + } } else { - const Tensor<2, dim, VectorizedArrayType> &jac = - this->cell_type == internal::MatrixFreeFunctions::general ? - this->jacobian[q_point] : - this->jacobian[0]; - const VectorizedArrayType fac = - (this->cell_type == internal::MatrixFreeFunctions::general ? - this->J_value[q_point] : - this->J_value[0] * this->quadrature_weights[q_point]) * - div_in; - for (unsigned int d = 0; d < dim; ++d) + if (!is_face && + this->cell_type == internal::MatrixFreeFunctions::cartesian) { - for (unsigned int e = 0; e < dim; ++e) - this->gradients_quad[(d * dim + e) * nqp + q_point] = - jac[d][e] * fac; + const VectorizedArrayType fac = + this->J_value[0] * this->quadrature_weights[q_point] * div_in; + for (unsigned int d = 0; d < dim; ++d) + { + this->gradients_quad[(d * dim + d) * nqp + q_point] = + (fac * this->jacobian[0][d][d]); + for (unsigned int e = d + 1; e < dim; ++e) + { + this->gradients_quad[(d * dim + e) * nqp + q_point] = + VectorizedArrayType(); + this->gradients_quad[(e * dim + d) * nqp + q_point] = + VectorizedArrayType(); + } + } + } + else + { + const Tensor<2, dim, VectorizedArrayType> &jac = + this->cell_type == internal::MatrixFreeFunctions::general ? + this->jacobian[q_point] : + this->jacobian[0]; + const VectorizedArrayType fac = + (this->cell_type == internal::MatrixFreeFunctions::general ? + this->J_value[q_point] : + this->J_value[0] * this->quadrature_weights[q_point]) * + div_in; + for (unsigned int d = 0; d < dim; ++d) + { + for (unsigned int e = 0; e < dim; ++e) + this->gradients_quad[(d * dim + e) * nqp + q_point] = + jac[d][e] * fac; + } } } } @@ -5945,6 +6373,12 @@ FEEvaluationAccess:: const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad, const unsigned int q_point) { + // TODO + AssertThrow( + this->data->element_type != + internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas, + ExcNotImplemented()); + // could have used base class operator, but that involves some overhead // which is inefficient. it is nice to have the symmetric tensor because // that saves some operations diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index 26df97bf9e..f28639e420 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -735,8 +735,14 @@ protected: const Point *quadrature_points; /** - * A pointer to the Jacobian information of the present cell. Only set to a - * useful value if on a non-Cartesian cell. + * A pointer to the inverse transpose Jacobian information of the present + * cell. Only set to a useful value if on a non-Cartesian cell. If the cell is + * Cartesian/affine then the Jacobian is stored at index 1. For faces on + * hypercube elements, the derivatives are reorder s.t the derivative + * orthogonal to the face is stored last, i.e for dim = 3 and face_no = 0 or + * 1, the derivatives are ordered as [dy, dz, dx], face_no = 2 or 3: [dz, dx, + * dy], and face_no = 5 or 6: [dx, dy, dz]. If the Jacobian also is stored, + * the components are instead reordered in the same way. */ const Tensor<2, dim, Number> *jacobian; diff --git a/include/deal.II/matrix_free/mapping_info.templates.h b/include/deal.II/matrix_free/mapping_info.templates.h index 209d4e65ac..57a1244736 100644 --- a/include/deal.II/matrix_free/mapping_info.templates.h +++ b/include/deal.II/matrix_free/mapping_info.templates.h @@ -2177,6 +2177,18 @@ namespace internal vv, my_data.jacobians[0][offset + q][d][e]); } + if (face_type[face] <= affine) + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = + ExtractFaceHelper::reorder_face_derivative_indices< + dim>(interior_face_no, e); + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + jac[d][ee], + vv, + my_data.jacobians[0][offset + q + 1][d][e]); + } if (update_flags_faces & update_jacobian_grads) { @@ -2281,6 +2293,18 @@ namespace internal vv, my_data.jacobians[1][offset + q][d][e]); } + if (face_type[face] <= affine) + for (unsigned int e = 0; e < dim; ++e) + { + const unsigned int ee = ExtractFaceHelper:: + reorder_face_derivative_indices( + exterior_face_no, e); + for (unsigned int d = 0; d < dim; ++d) + store_vectorized_array( + jac[d][ee], + vv, + my_data.jacobians[1][offset + q + 1][d][e]); + } if (update_flags_faces & update_jacobian_grads) { @@ -2762,7 +2786,7 @@ namespace internal max_size = std::max(max_size, my_data.data_index_offsets[face] + - (face_type[face] <= affine ? 1 : n_q_points)); + (face_type[face] <= affine ? 2 : n_q_points)); } const UpdateFlags update_flags_common = diff --git a/tests/matrix_free/matrix_vector_rt_01.cc b/tests/matrix_free/matrix_vector_rt_01.cc index d358717dc3..94b361e139 100644 --- a/tests/matrix_free/matrix_vector_rt_01.cc +++ b/tests/matrix_free/matrix_vector_rt_01.cc @@ -14,10 +14,10 @@ // --------------------------------------------------------------------- // 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 +// of the FE_RaviartThomasNodal element by evaluating values + gradients +// as well as the divergence 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" @@ -40,5 +40,13 @@ test() AffineConstraints constraints; constraints.close(); - do_test(dof, constraints); + + deallog << "Using " << 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; + do_test(dof, constraints, TestType::values); + do_test(dof, constraints, TestType::gradients); + do_test(dof, constraints, TestType::divergence); } diff --git a/tests/matrix_free/matrix_vector_rt_01.output b/tests/matrix_free/matrix_vector_rt_01.output index 0e53c41bb8..4c1b0e9e9a 100644 --- a/tests/matrix_free/matrix_vector_rt_01.output +++ b/tests/matrix_free/matrix_vector_rt_01.output @@ -1,25 +1,53 @@ -DEAL:2d::Testing FE_RaviartThomasNodal<2>(1) +DEAL:2d::Using 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::Testing Values +DEAL:2d::Norm of difference: 6.05894e-16 DEAL:2d:: -DEAL:2d::Testing FE_RaviartThomasNodal<2>(2) +DEAL:2d::Testing Gradients +DEAL:2d::Norm of difference: 6.33970e-16 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 5.24419e-16 +DEAL:2d:: +DEAL:2d::Using 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::Testing Values +DEAL:2d::Norm of difference: 6.67799e-16 +DEAL:2d:: +DEAL:2d::Testing Gradients +DEAL:2d::Norm of difference: 1.06191e-15 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 5.38383e-16 DEAL:2d:: -DEAL:3d::Testing FE_RaviartThomasNodal<3>(1) +DEAL:3d::Using 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::Testing Values +DEAL:3d::Norm of difference: 6.43004e-16 DEAL:3d:: -DEAL:3d::Testing FE_RaviartThomasNodal<3>(2) +DEAL:3d::Testing Gradients +DEAL:3d::Norm of difference: 8.38689e-16 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 8.84868e-16 +DEAL:3d:: +DEAL:3d::Using 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::Testing Values +DEAL:3d::Norm of difference: 1.10927e-15 +DEAL:3d:: +DEAL:3d::Testing Gradients +DEAL:3d::Norm of difference: 1.50806e-15 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 1.72044e-15 DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_02.cc b/tests/matrix_free/matrix_vector_rt_02.cc new file mode 100644 index 0000000000..e9e03006f3 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_02.cc @@ -0,0 +1,53 @@ +// --------------------------------------------------------------------- +// +// 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 it the same as matrix_vector_rt_01.cc however with +// non-Cartesian (but still affine) cells. + +#include + +#include "../tests.h" + +#include "matrix_vector_rt_common.h" + +template +void +test() +{ + Triangulation tria; + const unsigned int n_subdivisions = 2; + Point corners[dim]; + corners[0] = (dim == 2) ? Point(1, 0) : Point(1, 0, 0); + corners[1] = (dim == 2) ? Point(0.5, 0.5) : Point(0.5, 1, 0); + if (dim == 3) + corners[2] = Point(0.5, 0, 1); + GridGenerator::subdivided_parallelepiped(tria, n_subdivisions, corners); + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + + deallog << "Using " << 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; + + do_test(dof, constraints, TestType::values_gradients); + do_test(dof, constraints, TestType::divergence); +} diff --git a/tests/matrix_free/matrix_vector_rt_02.output b/tests/matrix_free/matrix_vector_rt_02.output new file mode 100644 index 0000000000..ecae4aa430 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_02.output @@ -0,0 +1,41 @@ + +DEAL:2d::Using FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 40 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 3.77395e-16 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 7.00345e-16 +DEAL:2d:: +DEAL:2d::Using FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 84 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 7.63491e-16 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 7.46533e-16 +DEAL:2d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 8 +DEAL:3d::Number of degrees of freedom: 240 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 1.22886e-15 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 1.16405e-15 +DEAL:3d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(2) +DEAL:3d::Number of cells: 8 +DEAL:3d::Number of degrees of freedom: 756 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 2.06049e-15 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 1.23490e-15 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_common.h b/tests/matrix_free/matrix_vector_rt_common.h index e2d0ae2324..9f2ee3d0e2 100644 --- a/tests/matrix_free/matrix_vector_rt_common.h +++ b/tests/matrix_free/matrix_vector_rt_common.h @@ -61,7 +61,38 @@ template void test(); +enum TestType : unsigned char +{ + values = 0, + values_gradients = 1, + gradients = 2, + divergence = 3 +}; +std::string +enum_to_string(TestType const enum_type) +{ + std::string string_type; + switch (enum_type) + { + case TestType::values: + string_type = "Values "; + break; + case TestType::gradients: + string_type = "Gradients "; + break; + case TestType::values_gradients: + string_type = "Values and Gradients "; + break; + case TestType::divergence: + string_type = "Divergence "; + break; + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + return string_type; +} template &data_in) - : data(data_in){}; + MatrixFreeTest(const MatrixFree &data_in, + const TestType test_type) + : data(data_in) + , test_type(test_type) + { + evaluation_flag = + (test_type == TestType::values) ? + EvaluationFlags::values : + ((test_type == TestType::gradients) ? + EvaluationFlags::gradients : + ((test_type == TestType::values_gradients) ? + EvaluationFlags::values | EvaluationFlags::gradients : + EvaluationFlags::gradients)); + }; virtual ~MatrixFreeTest(){}; @@ -83,29 +126,23 @@ public: { 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); + fe_eval.gather_evaluate(src, evaluation_flag); 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); + if (test_type < TestType::gradients) + fe_eval.submit_value(Number(10) * fe_eval.get_value(q), q); + if (test_type == TestType::gradients || + test_type == TestType::values_gradients) + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + else if (test_type == TestType::divergence) + fe_eval.submit_divergence(fe_eval.get_divergence(q), q); } - fe_eval.integrate_scatter(EvaluationFlags::values | - EvaluationFlags::gradients, - dst); + fe_eval.integrate_scatter(evaluation_flag, dst); } }; @@ -117,21 +154,19 @@ public: }; protected: - const MatrixFree &data; + const MatrixFree & data; + EvaluationFlags::EvaluationFlags evaluation_flag; + const TestType test_type; }; template void do_test(const DoFHandler & dof, - const AffineConstraints &constraints) + const AffineConstraints &constraints, + const TestType test_type) { - 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; - + deallog << "Testing " << enum_to_string(test_type) << std::endl; // constraints.distribute(solution); MatrixFree mf_data; @@ -157,7 +192,7 @@ do_test(const DoFHandler & dof, } // MatrixFree solution - MatrixFreeTest mf(mf_data); + MatrixFreeTest mf(mf_data, test_type); mf.test_functions(solution, initial_condition); @@ -191,16 +226,23 @@ do_test(const DoFHandler & dof, { 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); + const Number div_phi_i = fe_val[velocities].divergence(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); + const Number div_phi_j = fe_val[velocities].divergence(j, q); + + if (test_type < TestType::gradients) + local_matrix(i, j) += phi_j * phi_i * fe_val.JxW(q); + if (test_type == TestType::gradients || + test_type == TestType::values_gradients) + local_matrix(i, j) += + scalar_product(grad_phi_i, grad_phi_j) * fe_val.JxW(q); + else if (test_type == TestType::divergence) + local_matrix(i, j) += div_phi_i * div_phi_j * fe_val.JxW(q); } } cell->get_dof_indices(local_dof_indices); diff --git a/tests/matrix_free/matrix_vector_rt_face_01.cc b/tests/matrix_free/matrix_vector_rt_face_01.cc index efb1e82ee6..be2e9a78fa 100644 --- a/tests/matrix_free/matrix_vector_rt_face_01.cc +++ b/tests/matrix_free/matrix_vector_rt_face_01.cc @@ -14,10 +14,10 @@ // --------------------------------------------------------------------- // 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 +// of the FE_RaviartThomasNodal element by evaluating values + gradients +// as well as the divergence on faces and comparing the result with +// FEFaceVaules which is considered the reference. The mesh is a hypercube +// mesh with no hanging nodes and no other constraints. #include "../tests.h" @@ -40,5 +40,11 @@ test() AffineConstraints constraints; constraints.close(); - do_test(dof, constraints); + deallog << "Using " << 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; + do_test(dof, constraints, TestType::values_gradients); + do_test(dof, constraints, TestType::divergence); } diff --git a/tests/matrix_free/matrix_vector_rt_face_01.output b/tests/matrix_free/matrix_vector_rt_face_01.output index 99a9ae91f1..2d1aea3333 100644 --- a/tests/matrix_free/matrix_vector_rt_face_01.output +++ b/tests/matrix_free/matrix_vector_rt_face_01.output @@ -1,19 +1,31 @@ -DEAL:2d::Testing FE_RaviartThomasNodal<2>(1) +DEAL:2d::Using 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::Testing Values and Gradients +DEAL:2d::Norm of difference: 2.36910e-16 DEAL:2d:: -DEAL:2d::Testing FE_RaviartThomasNodal<2>(2) +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 4.52039e-16 +DEAL:2d:: +DEAL:2d::Using 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::Testing Values and Gradients +DEAL:2d::Norm of difference: 7.31531e-15 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 2.97155e-15 DEAL:2d:: -DEAL:3d::Testing FE_RaviartThomasNodal<3>(1) +DEAL:3d::Using 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::Testing Values and Gradients +DEAL:3d::Norm of difference: 3.73456e-15 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 4.51598e-15 DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_rt_face_02.cc b/tests/matrix_free/matrix_vector_rt_face_02.cc new file mode 100644 index 0000000000..060aa8bba4 --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_02.cc @@ -0,0 +1,51 @@ +// --------------------------------------------------------------------- +// +// 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 it the same as matrix_vector_rt_face_01.cc however with +// non-Cartesian (but still affine) cells. + +#include "../tests.h" + +#include "matrix_vector_rt_face_common.h" + + +template +void +test() +{ + Triangulation tria; + const unsigned int n_subdivisions = 2; + Point corners[dim]; + corners[0] = (dim == 2) ? Point(1, 0) : Point(1, 0, 0); + corners[1] = (dim == 2) ? Point(0.5, 0.5) : Point(0.5, 1, 0.25); + if (dim == 3) + corners[2] = Point(0.5, 0, 1); + GridGenerator::subdivided_parallelepiped(tria, n_subdivisions, corners); + + FE_RaviartThomasNodal fe(fe_degree - 1); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + + deallog << "Using " << 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; + do_test(dof, constraints, TestType::values_gradients); + do_test(dof, constraints, TestType::divergence); +} diff --git a/tests/matrix_free/matrix_vector_rt_face_02.output b/tests/matrix_free/matrix_vector_rt_face_02.output new file mode 100644 index 0000000000..d107fc460f --- /dev/null +++ b/tests/matrix_free/matrix_vector_rt_face_02.output @@ -0,0 +1,31 @@ + +DEAL:2d::Using FE_RaviartThomasNodal<2>(1) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 40 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 6.92435e-16 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 2.84791e-16 +DEAL:2d:: +DEAL:2d::Using FE_RaviartThomasNodal<2>(2) +DEAL:2d::Number of cells: 4 +DEAL:2d::Number of degrees of freedom: 84 +DEAL:2d:: +DEAL:2d::Testing Values and Gradients +DEAL:2d::Norm of difference: 4.97871e-15 +DEAL:2d:: +DEAL:2d::Testing Divergence +DEAL:2d::Norm of difference: 4.05232e-15 +DEAL:2d:: +DEAL:3d::Using FE_RaviartThomasNodal<3>(1) +DEAL:3d::Number of cells: 8 +DEAL:3d::Number of degrees of freedom: 240 +DEAL:3d:: +DEAL:3d::Testing Values and Gradients +DEAL:3d::Norm of difference: 2.96549e-15 +DEAL:3d:: +DEAL:3d::Testing Divergence +DEAL:3d::Norm of difference: 4.91752e-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 index 922ab008dc..b9a1cd1cef 100644 --- a/tests/matrix_free/matrix_vector_rt_face_common.h +++ b/tests/matrix_free/matrix_vector_rt_face_common.h @@ -61,7 +61,38 @@ template void test(); +enum TestType : unsigned char +{ + values = 0, + values_gradients = 1, + gradients = 2, + divergence = 3 +}; +std::string +enum_to_string(TestType const enum_type) +{ + std::string string_type; + switch (enum_type) + { + case TestType::values: + string_type = "Values "; + break; + case TestType::gradients: + string_type = "Gradients "; + break; + case TestType::values_gradients: + string_type = "Values and Gradients "; + break; + case TestType::divergence: + string_type = "Divergence "; + break; + default: + AssertThrow(false, ExcNotImplemented()); + break; + } + return string_type; +} template &data_in) - : data(data_in){}; + MatrixFreeTest(const MatrixFree &data_in, + const TestType test_type) + : data(data_in) + , test_type(test_type) + { + evaluation_flag = + (test_type == TestType::values) ? + EvaluationFlags::values : + ((test_type == TestType::gradients) ? + EvaluationFlags::gradients : + ((test_type == TestType::values_gradients) ? + EvaluationFlags::values | EvaluationFlags::gradients : + EvaluationFlags::gradients)); + }; virtual ~MatrixFreeTest(){}; @@ -93,42 +136,34 @@ public: 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); - + fe_eval.gather_evaluate(src, evaluation_flag); + fe_eval_n.gather_evaluate(src, evaluation_flag); 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); + if (test_type < TestType::gradients) + { + fe_eval.submit_value(10. * fe_eval.get_value(q), q); + fe_eval_n.submit_value(10. * fe_eval_n.get_value(q), q); + } + if (test_type == TestType::gradients || + test_type == TestType::values_gradients) + { + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + fe_eval_n.submit_gradient(fe_eval_n.get_gradient(q), q); + } + else if (test_type == TestType::divergence) + { + fe_eval.submit_divergence(fe_eval.get_divergence(q), q); + fe_eval_n.submit_divergence(fe_eval_n.get_divergence(q), q); + } } - - fe_eval.integrate_scatter(EvaluationFlags::values | - EvaluationFlags::gradients, - dst); - fe_eval_n.integrate_scatter(EvaluationFlags::values | - EvaluationFlags::gradients, - dst); + fe_eval.integrate_scatter(evaluation_flag, dst); + fe_eval_n.integrate_scatter(evaluation_flag, dst); } }; @@ -142,29 +177,23 @@ public: 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); + fe_eval.gather_evaluate(src, evaluation_flag); 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); + if (test_type < TestType::gradients) + fe_eval.submit_value(10. * fe_eval.get_value(q), q); + if (test_type == TestType::gradients || + test_type == TestType::values_gradients) + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + else if (test_type == TestType::divergence) + fe_eval.submit_divergence(fe_eval.get_divergence(q), q); } - fe_eval.integrate_scatter(EvaluationFlags::values | - EvaluationFlags::gradients, - dst); + fe_eval.integrate_scatter(evaluation_flag, dst); } }; @@ -181,7 +210,9 @@ public: }; protected: - const MatrixFree &data; + const MatrixFree & data; + EvaluationFlags::EvaluationFlags evaluation_flag; + const TestType test_type; }; @@ -189,16 +220,11 @@ protected: template void do_test(const DoFHandler & dof, - const AffineConstraints &constraints) + const AffineConstraints &constraints, + const TestType test_type) { - 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; - + deallog << "Testing " << enum_to_string(test_type) << std::endl; - // constraints.distribute(solution); MatrixFree mf_data; { const QGaussLobatto<1> quad(fe_degree + 2); @@ -226,7 +252,7 @@ do_test(const DoFHandler & dof, initial_condition[i] = random_value(); } - MatrixFreeTest mf(mf_data); + MatrixFreeTest mf(mf_data, test_type); mf.test_functions(solution, initial_condition); @@ -262,15 +288,23 @@ do_test(const DoFHandler & dof, 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); + const Number div_phi_i = fe_val[velocities].divergence(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); + const Number div_phi_j = fe_val[velocities].divergence(j, q); + + if (test_type < TestType::gradients) + local_matrix(i, j) += phi_j * phi_i * fe_val.JxW(q); + if (test_type == TestType::gradients || + test_type == TestType::values_gradients) + local_matrix(i, j) += + scalar_product(grad_phi_i, grad_phi_j) * fe_val.JxW(q); + else if (test_type == TestType::divergence) + local_matrix(i, j) += div_phi_i * div_phi_j * fe_val.JxW(q); } } cell->get_dof_indices(local_dof_indices);