From 41c9c46e52e40de4c8acbf355bc076bb733cb7ab Mon Sep 17 00:00:00 2001 From: Niklas Wik Date: Mon, 23 May 2022 13:04:08 +0200 Subject: [PATCH] Avoid copy paste Fix documentation --- include/deal.II/matrix_free/fe_evaluation.h | 65 +------------------ .../deal.II/matrix_free/fe_evaluation_data.h | 14 ++-- 2 files changed, 10 insertions(+), 69 deletions(-) diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 5b15ee71e2..4e36995fee 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -6221,68 +6221,9 @@ FEEvaluationAccess:: internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas) { // Piola transform is required - -# 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; - if (!is_face && - this->cell_type == internal::MatrixFreeFunctions::cartesian) - { - // Cartesian cell - 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; - } - else if (this->cell_type <= internal::MatrixFreeFunctions::affine) - { - // Affine cell - 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; - } - } - else - { - // General cell - AssertThrow(false, ExcNotImplemented()); - } + const Tensor<2, dim, VectorizedArrayType> &grad = grad_in; + FEEvaluationAccess:: + submit_gradient(grad, q_point); } else { diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h index 7ab5970987..1983d8bb0c 100644 --- a/include/deal.II/matrix_free/fe_evaluation_data.h +++ b/include/deal.II/matrix_free/fe_evaluation_data.h @@ -736,13 +736,13 @@ protected: /** * 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. + * cell. Only the first inverse transpose Jacobian (q_point = 0) is set for + * Cartesian/affine cells, and the actual Jacobian is stored at index 1 + * instead. 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. * Filled from MappingInfoStorage.jacobians in * include/deal.II/matrix_free/mapping_info.templates.h */ -- 2.39.5