Use ndarray instead of c-array

Update comments and documentation

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 00c9d4c7e518212c3122720e125fe8c02e15fe2d..fe15e2a34f3c7f3ecde946b78b4efc98ce00f8a0 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -19,6 +19,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/vectorization.h>
 
@@ -3158,9 +3159,8 @@ namespace internal
         (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
                                                     n_dofs_tangent;
 
-      const unsigned int component_table[3][3] = {{1, 2, 0},
-                                                  {2, 0, 1},
-                                                  {0, 1, 2}};
+      static constexpr dealii::ndarray<unsigned int, 3, 3> component_table = {
+        {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
       const unsigned int component =
         (dim == 2 && normal_dir == 0 && face_direction == 1) ?
           0 :
@@ -3335,9 +3335,8 @@ namespace internal
         (std::is_same<Eval0, EvalGeneral>::value) ? n_dofs_normal :
                                                     n_dofs_tangent;
 
-      const unsigned int component_table[3][3] = {{1, 2, 0},
-                                                  {2, 0, 1},
-                                                  {0, 1, 2}};
+      static constexpr dealii::ndarray<unsigned int, 3, 3> component_table = {
+        {{{1, 2, 0}}, {{2, 0, 1}}, {{0, 1, 2}}}};
       const unsigned int component =
         (dim == 2 && normal_dir == 0 && face_direction == 1) ?
           0 :

diff --git a/include/deal.II/matrix_free/fe_evaluation.h b/include/deal.II/matrix_free/fe_evaluation.h
index 516b6cb4c5d1fdf07102d0e667aceeecac61c482..bad7578625cd64f8b282d6edfe472ea1af22a4fc 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation.h
+++ b/include/deal.II/matrix_free/fe_evaluation.h
@@ -5725,10 +5725,10 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
 FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::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)
     {
+      // Piola transform is required
 #  ifdef DEBUG
       Assert(this->values_quad_initialized == true,
              internal::ExcAccessToUninitializedField());
@@ -5741,23 +5741,22 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
       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)
         {
+          // Cartesian cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> jac =
             this->jacobian[1];
           const VectorizedArrayType inv_det = determinant(this->jacobian[0]);
 
+          // J * u * det(J^-1)
           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];
+                              jac[comp][comp] * inv_det;
         }
-
-      // Affine or general cell
       else
         {
+          // Affine or general cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
             (this->cell_type > internal::MatrixFreeFunctions::affine) ?
               this->jacobian[q_point] :
@@ -5786,6 +5785,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::get_value(
     }
   else
     {
+      // No Piola needed
       return BaseClass::get_value(q_point);
     }
 }
@@ -5795,10 +5795,10 @@ inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
 FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   get_gradient(const unsigned int q_point) const
 {
-  // Check if Piola transform is required
   if (this->data->element_type ==
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
     {
+      // Piola transform is required
 #  ifdef DEBUG
       Assert(this->gradients_quad_initialized == true,
              internal::ExcAccessToUninitializedField());
@@ -5811,24 +5811,25 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
       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)
         {
+          // Cartesian cell
           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);
 
+          // J * grad_quad * J^-1 * det(J^-1)
           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)
         {
+          // Affine cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
             this->jacobian[0];
           const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
@@ -5853,9 +5854,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 grad_out[comp][d] = tmp;
               }
         }
-      // General cell TODO
       else
         {
+          // General cell
           // Here we need the jacobian gradient and not the inverse which is
           // stored in this->jacobian_gradients
           AssertThrow(false, ExcNotImplemented());
@@ -5890,40 +5891,41 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   if (this->data->element_type ==
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
     {
-      // Affine cell
       if (this->cell_type <= internal::MatrixFreeFunctions::affine)
         {
+          // Affine cell
           // 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]);
 
+          // div * det(J^-1)
           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
         {
+          // General cell
           Assert(false, ExcNotImplemented());
         }
     }
   else
     {
-      // Cartesian cell
       if (!is_face &&
           this->cell_type == internal::MatrixFreeFunctions::cartesian)
         {
+          // Cartesian cell
           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
         {
+          // cell with general/constant Jacobian
           const Tensor<2, dim, VectorizedArrayType> &jac =
             this->cell_type == internal::MatrixFreeFunctions::general ?
               this->jacobian[q_point] :
@@ -6039,13 +6041,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   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)
     {
+      // Piola transform is required
       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"));
@@ -6066,9 +6066,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
             this->values_quad[comp * nqp + q_point] =
               val_in[comp] * weight * jac[comp][comp];
         }
-      // Affine or general cell
       else
         {
+          // Affine or general cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
             (this->cell_type > internal::MatrixFreeFunctions::affine) ?
               this->jacobian[q_point] :
@@ -6090,7 +6090,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                   -determinant(inv_t_jac) :
                   determinant(inv_t_jac)));
 
-          // J^T * u * w
+          // J^T * u * factor
           for (unsigned int comp = 0; comp < n_components; ++comp)
             {
               this->values_quad[comp * nqp + q_point] =
@@ -6103,6 +6103,7 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
     }
   else
     {
+      // No Piola transform
       BaseClass::submit_value(val_in, q_point);
     }
 }
@@ -6113,10 +6114,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   submit_gradient(const Tensor<2, dim, VectorizedArrayType> grad_in,
                   const unsigned int                        q_point)
 {
-  // Check if Piola transform is required
   if (this->data->element_type ==
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
     {
+      // Piola transform is required
+
 #  ifdef DEBUG
       Assert(this->is_reinitialized, ExcNotInitialized());
 #  endif
@@ -6132,10 +6134,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
 #  endif
 
       const std::size_t nqp = this->n_quadrature_points;
-      // Cartesian cell
       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];
@@ -6145,9 +6147,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
               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)
         {
+          // Affine cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
             this->jacobian[0];
           const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
@@ -6173,9 +6175,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
               }
         }
-      // General cell TODO
       else
         {
+          // General cell
           AssertThrow(false, ExcNotImplemented());
         }
     }
@@ -6194,10 +6196,11 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
     const Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>> grad_in,
     const unsigned int                                        q_point)
 {
-  // Check if Piola transform is required
   if (this->data->element_type ==
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
     {
+      // Piola transform is required
+
 #  ifdef DEBUG
       Assert(this->is_reinitialized, ExcNotInitialized());
 #  endif
@@ -6213,10 +6216,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
 #  endif
 
       const std::size_t nqp = this->n_quadrature_points;
-      // Cartesian cell
       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];
@@ -6226,9 +6229,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
               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)
         {
+          // Affine cell
           const Tensor<2, dim, dealii::VectorizedArray<Number>> &inv_t_jac =
             this->jacobian[0];
           const Tensor<2, dim, VectorizedArrayType> &jac = this->jacobian[1];
@@ -6254,9 +6257,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 this->gradients_quad[(comp * dim + d) * nqp + q_point] = tmp;
               }
         }
-      // General cell TODO
       else
         {
+          // General cell
           AssertThrow(false, ExcNotImplemented());
         }
     }
@@ -6292,9 +6295,10 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
   if (this->data->element_type ==
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas)
     {
-      // Affine cell
       if (this->cell_type <= internal::MatrixFreeFunctions::affine)
         {
+          // Affine cell
+
           // Derivatives are reordered for faces. Need to take this into account
           // and 1/inv_det != J_value for faces
           const VectorizedArrayType fac =
@@ -6317,9 +6321,9 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
                 }
             }
         }
-      // General cell TODO
       else
         {
+          // General cell
           AssertThrow(false, ExcNotImplemented());
         }
     }
@@ -6373,7 +6377,6 @@ FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
     const SymmetricTensor<2, dim, VectorizedArrayType> sym_grad,
     const unsigned int                                 q_point)
 {
-  // TODO
   AssertThrow(
     this->data->element_type !=
       internal::MatrixFreeFunctions::ElementType::tensor_raviart_thomas,

diff --git a/include/deal.II/matrix_free/fe_evaluation_data.h b/include/deal.II/matrix_free/fe_evaluation_data.h
index f28639e420ac65641840e081d05689fee2300d25..7ab59709879481acad00e8b25e69cb0e7e69dac5 100644 (file)
--- a/include/deal.II/matrix_free/fe_evaluation_data.h
+++ b/include/deal.II/matrix_free/fe_evaluation_data.h
@@ -743,6 +743,8 @@ protected:
    * 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
    */
   const Tensor<2, dim, Number> *jacobian;
 

diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h
index e70c65033b0e4c62fccefdbc47d0010d2f0b6ad0..18bbba346251d5870ca0b16438b55322c5861b1a 100644 (file)
--- a/include/deal.II/matrix_free/mapping_info_storage.h
+++ b/include/deal.II/matrix_free/mapping_info_storage.h
@@ -227,6 +227,14 @@ namespace internal
        * Contains two fields for access from both sides for interior faces,
        * but the default case (cell integrals or boundary integrals) only
        * fills the zeroth component and ignores the first one.
+       *
+       * If the cell is Cartesian/affine then the Jacobian is stored at index 1
+       * of the AlignedVector. 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.
        */
       std::array<AlignedVector<Tensor<2, spacedim, Number>>, 2> jacobians;