Avoid copy paste

Fix documentation

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 5b15ee71e2cf6bc9a15157e89344f1a3345e5b00..4e36995fee9a498466aa7bd38c78d8b7b3e4e755 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -6221,68 +6221,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
       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<Number>> &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<dim, dim, Number, is_face, VectorizedArrayType>::
+        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 7ab59709879481acad00e8b25e69cb0e7e69dac5..1983d8bb0ce00cf94c7edd73005caeb7cb2d1e9d 100644 (file)
--- 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
    */