From: Martin Kronbichler Date: Thu, 31 Aug 2023 07:12:47 +0000 (+0200) Subject: FEEvaluation for Raviart-Thomas: use fewer multiplications X-Git-Tag: relicensing~533^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15950%2Fhead;p=dealii.git FEEvaluation for Raviart-Thomas: use fewer multiplications --- diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h index 3c291b67ef..d1acf28d89 100644 --- a/include/deal.II/matrix_free/fe_evaluation.h +++ b/include/deal.II/matrix_free/fe_evaluation.h @@ -5847,11 +5847,11 @@ FEEvaluationAccess::get_value( // 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; + value_out[comp] = this->values_quad[q_point] * jac[comp][0]; for (unsigned int e = 1; e < dim; ++e) value_out[comp] += - this->values_quad[e * nqp + q_point] * jac[comp][e] * inv_det; + this->values_quad[e * nqp + q_point] * jac[comp][e]; + value_out[comp] *= inv_det; } } return value_out; @@ -5863,6 +5863,8 @@ FEEvaluationAccess::get_value( } } + + template inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType> FEEvaluationAccess:: @@ -5901,7 +5903,7 @@ FEEvaluationAccess:: 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; + inv_t_jac[d][d] * (jac[comp][comp] * inv_det); } else if (this->cell_type <= internal::MatrixFreeFunctions::affine) { @@ -5924,10 +5926,10 @@ FEEvaluationAccess:: 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 * + tmp += jac[comp][f] * inv_t_jac[d][e] * this->gradients_quad[(f * dim + e) * nqp + q_point]; - grad_out[comp][d] = tmp; + grad_out[comp][d] = tmp * inv_det; } } else @@ -6067,10 +6069,10 @@ FEEvaluationAccess:: this->jacobian[0][2][2]; // div * det(J^-1) - divergence = this->gradients_quad[q_point] * inv_det; + divergence = this->gradients_quad[q_point]; for (unsigned int d = 1; d < dim; ++d) - divergence += - this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det; + divergence += this->gradients_quad[(dim * d + d) * nqp + q_point]; + divergence *= inv_det; } else { @@ -6085,10 +6087,10 @@ FEEvaluationAccess:: Number((is_face && dim == 2 && this->get_face_no() < 2) ? -1 : 1); // div * det(J^-1) - divergence = this->gradients_quad[q_point] * inv_det; + divergence = this->gradients_quad[q_point]; for (unsigned int d = 1; d < dim; ++d) - divergence += - this->gradients_quad[(dim * d + d) * nqp + q_point] * inv_det; + divergence += this->gradients_quad[(dim * d + d) * nqp + q_point]; + divergence *= inv_det; } } else @@ -6247,11 +6249,11 @@ FEEvaluationAccess:: else { // Affine or general cell - const Tensor<2, dim, VectorizedArrayType> &inv_t_jac = + const Tensor<2, dim, VectorizedArrayType> inv_t_jac = (this->cell_type > internal::MatrixFreeFunctions::affine) ? this->jacobian[q_point] : this->jacobian[0]; - const Tensor<2, dim, VectorizedArrayType> &jac = + const Tensor<2, dim, VectorizedArrayType> jac = (this->cell_type > internal::MatrixFreeFunctions::affine) ? transpose(invert(inv_t_jac)) : this->jacobian[1]; @@ -6272,10 +6274,11 @@ FEEvaluationAccess:: for (unsigned int comp = 0; comp < n_components; ++comp) { this->values_quad[comp * nqp + q_point] = - val_in[0] * jac[0][comp] * fac; + val_in[0] * jac[0][comp]; for (unsigned int e = 1; e < dim; ++e) this->values_quad[comp * nqp + q_point] += - val_in[e] * jac[e][comp] * fac; + val_in[e] * jac[e][comp]; + this->values_quad[comp * nqp + q_point] *= fac; } } }