]> https://gitweb.dealii.org/ - dealii.git/commitdiff
refactoring ShapeInfo's 1D data into UnivariateShapeData
authorJulius Witte <julius.witte@iwr.uni-heidelberg.de>
Thu, 12 Mar 2020 15:56:30 +0000 (16:56 +0100)
committerJulius Witte <julius.witte@iwr.uni-heidelberg.de>
Thu, 12 Mar 2020 20:36:34 +0000 (21:36 +0100)
include/deal.II/matrix_free/cuda_matrix_free.templates.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_selector.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/matrix_free.templates.h
include/deal.II/matrix_free/operators.h
include/deal.II/matrix_free/shape_info.h
include/deal.II/matrix_free/shape_info.templates.h
tests/matrix_free/shape_info_inverse.cc

index dc91f7d444a325aee4360249f404148fc3cf0afd..df97cbe8e15e50a81c2b250f15fbb49e99baf317 100644 (file)
@@ -875,7 +875,7 @@ namespace CUDAWrappers
 
     cudaError_t cuda_error =
       cudaMemcpyToSymbol(internal::get_global_shape_values<Number>(),
-                         shape_info.shape_values.data(),
+                         shape_info.data.front().shape_values.data(),
                          size_shape_values,
                          0,
                          cudaMemcpyHostToDevice);
@@ -885,7 +885,7 @@ namespace CUDAWrappers
       {
         cuda_error =
           cudaMemcpyToSymbol(internal::get_global_shape_gradients<Number>(),
-                             shape_info.shape_gradients.data(),
+                             shape_info.data.front().shape_gradients.data(),
                              size_shape_values,
                              0,
                              cudaMemcpyHostToDevice);
index 79623b3e8ad805d799d1b50f7c884644b3de3ef6..d559935ff21783343825a868c6739d3d024292c6 100644 (file)
@@ -163,14 +163,17 @@ namespace internal
                                         fe_degree + 1,
                                         n_q_points_1d,
                                         Number>;
-    Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
-                                            shape_info.shape_values,
-              variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
-                                            shape_info.shape_gradients,
-              variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
-                                            shape_info.shape_hessians,
-              shape_info.fe_degree + 1,
-              shape_info.n_q_points_1d);
+    Eval eval(variant == evaluate_evenodd ?
+                shape_info.data.front().shape_values_eo :
+                shape_info.data.front().shape_values,
+              variant == evaluate_evenodd ?
+                shape_info.data.front().shape_gradients_eo :
+                shape_info.data.front().shape_gradients,
+              variant == evaluate_evenodd ?
+                shape_info.data.front().shape_hessians_eo :
+                shape_info.data.front().shape_hessians,
+              shape_info.data.front().fe_degree + 1,
+              shape_info.data.front().n_q_points_1d);
 
     const unsigned int temp_size =
       Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
@@ -183,9 +186,10 @@ namespace internal
     if (temp_size == 0)
       {
         temp1 = scratch_data;
-        temp2 = temp1 +
-                std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
-                         Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+        temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+                                   shape_info.data.front().fe_degree + 1),
+                                 Utilities::fixed_power<dim>(
+                                   shape_info.data.front().n_q_points_1d));
       }
     else
       {
@@ -197,7 +201,7 @@ namespace internal
       temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
     const unsigned int dofs_per_comp =
       (type == MatrixFreeFunctions::truncated_tensor) ?
-        Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
+        Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
         shape_info.dofs_per_component_on_cell;
     const Number *values_dofs = values_dofs_actual;
     if (type == MatrixFreeFunctions::truncated_tensor)
@@ -205,7 +209,8 @@ namespace internal
         Number *values_dofs_tmp =
           scratch_data + 2 * (std::max(shape_info.dofs_per_component_on_cell,
                                        shape_info.n_q_points));
-        const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+        const int degree =
+          fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
         unsigned int count_p = 0, count_q = 0;
         for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
           {
@@ -434,14 +439,17 @@ namespace internal
                                         fe_degree + 1,
                                         n_q_points_1d,
                                         Number>;
-    Eval eval(variant == evaluate_evenodd ? shape_info.shape_values_eo :
-                                            shape_info.shape_values,
-              variant == evaluate_evenodd ? shape_info.shape_gradients_eo :
-                                            shape_info.shape_gradients,
-              variant == evaluate_evenodd ? shape_info.shape_hessians_eo :
-                                            shape_info.shape_hessians,
-              shape_info.fe_degree + 1,
-              shape_info.n_q_points_1d);
+    Eval eval(variant == evaluate_evenodd ?
+                shape_info.data.front().shape_values_eo :
+                shape_info.data.front().shape_values,
+              variant == evaluate_evenodd ?
+                shape_info.data.front().shape_gradients_eo :
+                shape_info.data.front().shape_gradients,
+              variant == evaluate_evenodd ?
+                shape_info.data.front().shape_hessians_eo :
+                shape_info.data.front().shape_hessians,
+              shape_info.data.front().fe_degree + 1,
+              shape_info.data.front().n_q_points_1d);
 
     const unsigned int temp_size =
       Eval::n_rows_of_product == numbers::invalid_unsigned_int ?
@@ -454,9 +462,10 @@ namespace internal
     if (temp_size == 0)
       {
         temp1 = scratch_data;
-        temp2 = temp1 +
-                std::max(Utilities::fixed_power<dim>(shape_info.fe_degree + 1),
-                         Utilities::fixed_power<dim>(shape_info.n_q_points_1d));
+        temp2 = temp1 + std::max(Utilities::fixed_power<dim>(
+                                   shape_info.data.front().fe_degree + 1),
+                                 Utilities::fixed_power<dim>(
+                                   shape_info.data.front().n_q_points_1d));
       }
     else
       {
@@ -468,7 +477,7 @@ namespace internal
       temp_size == 0 ? shape_info.n_q_points : Eval::n_columns_of_product;
     const unsigned int dofs_per_comp =
       (type == MatrixFreeFunctions::truncated_tensor) ?
-        Utilities::fixed_power<dim>(shape_info.fe_degree + 1) :
+        Utilities::fixed_power<dim>(shape_info.data.front().fe_degree + 1) :
         shape_info.dofs_per_component_on_cell;
     // expand dof_values to tensor product for truncated tensor products
     Number *values_dofs =
@@ -612,7 +621,8 @@ namespace internal
       {
         values_dofs -= dofs_per_comp * n_components;
         unsigned int count_p = 0, count_q = 0;
-        const int degree = fe_degree != -1 ? fe_degree : shape_info.fe_degree;
+        const int    degree =
+          fe_degree != -1 ? fe_degree : shape_info.data.front().fe_degree;
         for (int i = 0; i < (dim > 2 ? degree + 1 : 1); ++i)
           {
             for (int j = 0; j < (dim > 1 ? degree + 1 - i : 1); ++j)
@@ -630,7 +640,8 @@ namespace internal
             count_q += i * (degree + 1);
           }
         AssertDimension(count_q,
-                        Utilities::fixed_power<dim>(shape_info.fe_degree + 1));
+                        Utilities::fixed_power<dim>(
+                          shape_info.data.front().fe_degree + 1));
       }
   }
 
@@ -1039,8 +1050,9 @@ namespace internal
     const bool evaluate_gradients,
     const bool evaluate_hessians)
   {
-    AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
-                    (fe_degree + 2) / 2 * (fe_degree + 1));
+    AssertDimension(
+      shape_info.data.front().shape_gradients_collocation_eo.size(),
+      (fe_degree + 2) / 2 * (fe_degree + 1));
 
     EvaluatorTensorProduct<evaluate_evenodd,
                            dim,
@@ -1048,8 +1060,8 @@ namespace internal
                            fe_degree + 1,
                            Number>
                            eval(AlignedVector<Number>(),
-           shape_info.shape_gradients_collocation_eo,
-           shape_info.shape_hessians_collocation_eo);
+           shape_info.data.front().shape_gradients_collocation_eo,
+           shape_info.data.front().shape_hessians_collocation_eo);
     constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
 
     for (unsigned int c = 0; c < n_components; c++)
@@ -1115,8 +1127,9 @@ namespace internal
     const bool integrate_gradients,
     const bool add_into_values_array)
   {
-    AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
-                    (fe_degree + 2) / 2 * (fe_degree + 1));
+    AssertDimension(
+      shape_info.data.front().shape_gradients_collocation_eo.size(),
+      (fe_degree + 2) / 2 * (fe_degree + 1));
 
     EvaluatorTensorProduct<evaluate_evenodd,
                            dim,
@@ -1124,8 +1137,8 @@ namespace internal
                            fe_degree + 1,
                            Number>
                            eval(AlignedVector<Number>(),
-           shape_info.shape_gradients_collocation_eo,
-           shape_info.shape_hessians_collocation_eo);
+           shape_info.data.front().shape_gradients_collocation_eo,
+           shape_info.data.front().shape_hessians_collocation_eo);
     constexpr unsigned int n_q_points = Utilities::pow(fe_degree + 1, dim);
 
     for (unsigned int c = 0; c < n_components; c++)
@@ -1241,7 +1254,7 @@ namespace internal
           n_q_points_1d,
           1,
           Number,
-          Number>::do_forward(shape_info.shape_values_eo,
+          Number>::do_forward(shape_info.data.front().shape_values_eo,
                               values_dofs,
                               values_quad);
 
@@ -1292,8 +1305,9 @@ namespace internal
                       "of lower degree, so the evaluation results would be "
                       "wrong. Thus, this class does not permit the desired "
                       "operation."));
-    AssertDimension(shape_info.shape_gradients_collocation_eo.size(),
-                    (n_q_points_1d + 1) / 2 * n_q_points_1d);
+    AssertDimension(
+      shape_info.data.front().shape_gradients_collocation_eo.size(),
+      (n_q_points_1d + 1) / 2 * n_q_points_1d);
     constexpr unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim);
 
     for (unsigned int c = 0; c < n_components; c++)
@@ -1318,7 +1332,7 @@ namespace internal
           n_q_points_1d,
           1,
           Number,
-          Number>::do_backward(shape_info.shape_values_eo,
+          Number>::do_backward(shape_info.data.front().shape_values_eo,
                                add_into_values_array,
                                values_quad,
                                values_dofs);
@@ -1360,29 +1374,29 @@ namespace internal
     {
       const AlignedVector<Number> &val1 =
         symmetric_evaluate ?
-          data.shape_values_eo :
+          data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_values :
-             data.values_within_subface[subface_index % 2]);
+             data.data.front().shape_values :
+             data.data.front().values_within_subface[subface_index % 2]);
       const AlignedVector<Number> &val2 =
         symmetric_evaluate ?
-          data.shape_values_eo :
+          data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_values :
-             data.values_within_subface[subface_index / 2]);
+             data.data.front().shape_values :
+             data.data.front().values_within_subface[subface_index / 2]);
 
       const AlignedVector<Number> &grad1 =
         symmetric_evaluate ?
-          data.shape_gradients_eo :
+          data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_gradients :
-             data.gradients_within_subface[subface_index % 2]);
+             data.data.front().shape_gradients :
+             data.data.front().gradients_within_subface[subface_index % 2]);
       const AlignedVector<Number> &grad2 =
         symmetric_evaluate ?
-          data.shape_gradients_eo :
+          data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_gradients :
-             data.gradients_within_subface[subface_index / 2]);
+             data.data.front().shape_gradients :
+             data.data.front().gradients_within_subface[subface_index / 2]);
 
       using Eval =
         internal::EvaluatorTensorProduct<symmetric_evaluate ?
@@ -1395,18 +1409,20 @@ namespace internal
       Eval eval1(val1,
                  grad1,
                  AlignedVector<Number>(),
-                 data.fe_degree + 1,
-                 data.n_q_points_1d);
+                 data.data.front().fe_degree + 1,
+                 data.data.front().n_q_points_1d);
       Eval eval2(val2,
                  grad2,
                  AlignedVector<Number>(),
-                 data.fe_degree + 1,
-                 data.n_q_points_1d);
+                 data.data.front().fe_degree + 1,
+                 data.data.front().n_q_points_1d);
 
       const unsigned int size_deg =
         fe_degree > -1 ?
           Utilities::pow(fe_degree + 1, dim - 1) :
-          (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
+          (dim > 1 ?
+             Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
+             1);
 
       const unsigned int n_q_points = fe_degree > -1 ?
                                         Utilities::pow(n_q_points_1d, dim - 1) :
@@ -1454,9 +1470,10 @@ namespace internal
                         n_q_points_1d,
                         n_q_points_1d,
                         Number>
-                        eval_grad(AlignedVector<Number>(),
-                                  data.shape_gradients_collocation_eo,
-                                  AlignedVector<Number>());
+                        eval_grad(
+                          AlignedVector<Number>(),
+                          data.data.front().shape_gradients_collocation_eo,
+                          AlignedVector<Number>());
                       eval_grad.template gradients<0, true, false>(
                         values_quad, gradients_quad);
                       eval_grad.template gradients<1, true, false>(
@@ -1521,29 +1538,29 @@ namespace internal
     {
       const AlignedVector<Number> &val1 =
         symmetric_evaluate ?
-          data.shape_values_eo :
+          data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_values :
-             data.values_within_subface[subface_index % 2]);
+             data.data.front().shape_values :
+             data.data.front().values_within_subface[subface_index % 2]);
       const AlignedVector<Number> &val2 =
         symmetric_evaluate ?
-          data.shape_values_eo :
+          data.data.front().shape_values_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_values :
-             data.values_within_subface[subface_index / 2]);
+             data.data.front().shape_values :
+             data.data.front().values_within_subface[subface_index / 2]);
 
       const AlignedVector<Number> &grad1 =
         symmetric_evaluate ?
-          data.shape_gradients_eo :
+          data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_gradients :
-             data.gradients_within_subface[subface_index % 2]);
+             data.data.front().shape_gradients :
+             data.data.front().gradients_within_subface[subface_index % 2]);
       const AlignedVector<Number> &grad2 =
         symmetric_evaluate ?
-          data.shape_gradients_eo :
+          data.data.front().shape_gradients_eo :
           (subface_index >= GeometryInfo<dim>::max_children_per_cell ?
-             data.shape_gradients :
-             data.gradients_within_subface[subface_index / 2]);
+             data.data.front().shape_gradients :
+             data.data.front().gradients_within_subface[subface_index / 2]);
 
       using Eval =
         internal::EvaluatorTensorProduct<symmetric_evaluate ?
@@ -1553,13 +1570,23 @@ namespace internal
                                          fe_degree + 1,
                                          n_q_points_1d,
                                          Number>;
-      Eval eval1(val1, grad1, val1, data.fe_degree + 1, data.n_q_points_1d);
-      Eval eval2(val2, grad2, val1, data.fe_degree + 1, data.n_q_points_1d);
+      Eval eval1(val1,
+                 grad1,
+                 val1,
+                 data.data.front().fe_degree + 1,
+                 data.data.front().n_q_points_1d);
+      Eval eval2(val2,
+                 grad2,
+                 val1,
+                 data.data.front().fe_degree + 1,
+                 data.data.front().n_q_points_1d);
 
       const unsigned int size_deg =
         fe_degree > -1 ?
           Utilities::pow(fe_degree + 1, dim - 1) :
-          (dim > 1 ? Utilities::fixed_power<dim - 1>(data.fe_degree + 1) : 1);
+          (dim > 1 ?
+             Utilities::fixed_power<dim - 1>(data.data.front().fe_degree + 1) :
+             1);
 
       const unsigned int n_q_points = fe_degree > -1 ?
                                         Utilities::pow(n_q_points_1d, dim - 1) :
@@ -1609,9 +1636,10 @@ namespace internal
                         n_q_points_1d,
                         n_q_points_1d,
                         Number>
-                        eval_grad(AlignedVector<Number>(),
-                                  data.shape_gradients_collocation_eo,
-                                  AlignedVector<Number>());
+                        eval_grad(
+                          AlignedVector<Number>(),
+                          data.data.front().shape_gradients_collocation_eo,
+                          AlignedVector<Number>());
                       if (integrate_val)
                         eval_grad.template gradients<1, false, true>(
                           gradients_quad + n_q_points, values_quad);
@@ -1687,10 +1715,10 @@ namespace internal
                                        fe_degree + 1,
                                        0,
                                        Number>
-        evalf(data.shape_data_on_face[face_no % 2],
+        evalf(data.data.front().shape_data_on_face[face_no % 2],
               AlignedVector<Number>(),
               AlignedVector<Number>(),
-              data.fe_degree + 1,
+              data.data.front().fe_degree + 1,
               0);
 
       const unsigned int in_stride = do_evaluate ?
@@ -1865,8 +1893,9 @@ namespace internal
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                          numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ? static_dofs_per_face :
-                         Utilities::pow(data.fe_degree + 1, dim - 1);
+        fe_degree > -1 ?
+          static_dofs_per_face :
+          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
 
       // we allocate small amounts of data on the stack to signal the compiler
       // that this temporary data is only needed for the calculations but the
@@ -1972,8 +2001,9 @@ namespace internal
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                          numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ? static_dofs_per_face :
-                         Utilities::pow(data.fe_degree + 1, dim - 1);
+        fe_degree > -1 ?
+          static_dofs_per_face :
+          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
 
       constexpr unsigned int stack_array_size_threshold = 100;
 
@@ -2062,8 +2092,9 @@ namespace internal
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
                          numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ? static_dofs_per_face :
-                         Utilities::pow(data.fe_degree + 1, dim - 1);
+        fe_degree > -1 ?
+          static_dofs_per_face :
+          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
 
       constexpr unsigned int stack_array_size_threshold = 100;
 
@@ -2079,7 +2110,7 @@ namespace internal
 
       // case 1: contiguous and interleaved indices
       if (((evaluate_gradients == false &&
-            data.nodal_at_cell_boundaries == true) ||
+            data.data.front().nodal_at_cell_boundaries == true) ||
            (data.element_type ==
               MatrixFreeFunctions::tensor_symmetric_hermite &&
             fe_degree > 1)) &&
@@ -2103,7 +2134,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 1 + side];
+                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
               const unsigned int *index_array =
@@ -2154,7 +2185,7 @@ namespace internal
 
       // case 2: contiguous and interleaved indices with fixed stride
       else if (((evaluate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2174,7 +2205,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 1 + side];
+                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2240,7 +2271,7 @@ namespace internal
 
       // case 3: contiguous and interleaved indices with mixed stride
       else if (((evaluate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2270,7 +2301,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 1 + side];
+                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2379,7 +2410,7 @@ namespace internal
 
       // case 4: contiguous indices without interleaving
       else if (((evaluate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2400,7 +2431,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 1 + side];
+                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2544,8 +2575,9 @@ namespace internal
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
                          numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         Utilities::pow(data.fe_degree + 1, dim - 1);
+        fe_degree > -1 ?
+          Utilities::pow(fe_degree + 1, dim - 1) :
+          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
 
       constexpr unsigned int stack_array_size_threshold = 100;
 
@@ -2597,7 +2629,7 @@ namespace internal
 
       // case 1: contiguous and interleaved indices
       if (((integrate_gradients == false &&
-            data.nodal_at_cell_boundaries == true) ||
+            data.data.front().nodal_at_cell_boundaries == true) ||
            (data.element_type ==
               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
             fe_degree > 1)) &&
@@ -2621,7 +2653,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
               const unsigned int *index_array =
@@ -2675,7 +2707,7 @@ namespace internal
 
       // case 2: contiguous and interleaved indices with fixed stride
       else if (((integrate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2695,7 +2727,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2764,7 +2796,7 @@ namespace internal
 
       // case 3: contiguous and interleaved indices with mixed stride
       else if (((integrate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2794,7 +2826,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2899,7 +2931,7 @@ namespace internal
 
       // case 4: contiguous indices without interleaving
       else if (((integrate_gradients == false &&
-                 data.nodal_at_cell_boundaries == true) ||
+                 data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
                    internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
@@ -2921,7 +2953,7 @@ namespace internal
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
               const unsigned int *index_array =
@@ -3062,9 +3094,10 @@ namespace internal
                                        fe_degree + 1,
                                        fe_degree + 1,
                                        Number>
-        evaluator(AlignedVector<Number>(),
-                  AlignedVector<Number>(),
-                  fe_eval.get_shape_info().inverse_shape_values_eo);
+        evaluator(
+          AlignedVector<Number>(),
+          AlignedVector<Number>(),
+          fe_eval.get_shape_info().data.front().inverse_shape_values_eo);
 
       for (unsigned int d = 0; d < n_components; ++d)
         {
index 3f865caaf8318c1dfbee38004025d8da832a0cf9..8f2438167957184f215815e671f5fc2d9b3ea801 100644 (file)
@@ -174,7 +174,7 @@ namespace internal
         const bool evaluate_gradients,
         const bool evaluate_hessians)
       {
-        const unsigned int     runtime_degree   = shape_info.fe_degree;
+        const unsigned int runtime_degree = shape_info.data.front().fe_degree;
         constexpr unsigned int start_n_q_points = degree + 1;
         if (runtime_degree == degree)
           Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
@@ -211,7 +211,7 @@ namespace internal
         const bool integrate_gradients,
         const bool sum_into_values_array = false)
       {
-        const int              runtime_degree   = shape_info.fe_degree;
+        const int runtime_degree = shape_info.data.front().fe_degree;
         constexpr unsigned int start_n_q_points = degree + 1;
         if (runtime_degree == degree)
           Factory<dim, n_components, Number, 1, degree, start_n_q_points>::
@@ -276,7 +276,7 @@ namespace internal
         const bool evaluate_gradients,
         const bool evaluate_hessians)
       {
-        const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+        const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d;
         if (runtime_n_q_points_1d == n_q_points_1d)
           {
             if (n_q_points_1d == degree + 1 &&
@@ -349,7 +349,7 @@ namespace internal
         const bool integrate_gradients,
         const bool sum_into_values_array)
       {
-        const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
+        const int runtime_n_q_points_1d = shape_info.data.front().n_q_points_1d;
         if (runtime_n_q_points_1d == n_q_points_1d)
           {
             if (n_q_points_1d == degree + 1 &&
index 8d4e14d7393a6d4b9f699b7c49876b280e42b876..6bb9519e87a291d6c8da700ce063ea5840bc0679 100644 (file)
@@ -3371,7 +3371,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   Assert(scratch_data_array != nullptr, ExcInternalError());
 
   const unsigned int tensor_dofs_per_component =
-    Utilities::fixed_power<dim>(this->data->fe_degree + 1);
+    Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
   const unsigned int dofs_per_component =
     this->data->dofs_per_component_on_cell;
   const unsigned int n_quadrature_points =
@@ -6677,7 +6677,8 @@ FEEvaluation<dim,
   // print error message when the dimensions do not match. Propose a possible
   // fix
   if ((static_cast<unsigned int>(fe_degree) != numbers::invalid_unsigned_int &&
-       static_cast<unsigned int>(fe_degree) != this->data->fe_degree) ||
+       static_cast<unsigned int>(fe_degree) !=
+         this->data->data.front().fe_degree) ||
       n_q_points != this->n_quadrature_points)
     {
       std::string message =
@@ -6703,7 +6704,8 @@ FEEvaluation<dim,
                    proposed_quad_comp = numbers::invalid_unsigned_int;
       if (dof_no != numbers::invalid_unsigned_int)
         {
-          if (static_cast<unsigned int>(fe_degree) == this->data->fe_degree)
+          if (static_cast<unsigned int>(fe_degree) ==
+              this->data->data.front().fe_degree)
             {
               proposed_dof_comp = dof_no;
               proposed_fe_comp  = first_selected_component;
@@ -6716,6 +6718,7 @@ FEEvaluation<dim,
                    ++nf)
                 if (this->matrix_info
                       ->get_shape_info(no, 0, nf, this->active_fe_index, 0)
+                      .data.front()
                       .fe_degree == static_cast<unsigned int>(fe_degree))
                   {
                     proposed_dof_comp = no;
@@ -6781,7 +6784,8 @@ FEEvaluation<dim,
         std::pow(1.001 * this->n_quadrature_points, 1. / dim));
       message += "Wrong template arguments:\n";
       message += "    Did you mean FEEvaluation<dim,";
-      message += Utilities::int_to_string(this->data->fe_degree) + ",";
+      message +=
+        Utilities::int_to_string(this->data->data.front().fe_degree) + ",";
       message += Utilities::int_to_string(proposed_n_q_points_1d);
       message += "," + Utilities::int_to_string(n_components);
       message += ",Number>(data";
@@ -6793,7 +6797,8 @@ FEEvaluation<dim,
         }
       message += ")?\n";
       std::string correct_pos;
-      if (this->data->fe_degree != static_cast<unsigned int>(fe_degree))
+      if (this->data->data.front().fe_degree !=
+          static_cast<unsigned int>(fe_degree))
         correct_pos = " ^";
       else
         correct_pos = "  ";
@@ -6803,7 +6808,8 @@ FEEvaluation<dim,
         correct_pos += "  \n";
       message += "                                 " + correct_pos;
 
-      Assert(static_cast<unsigned int>(fe_degree) == this->data->fe_degree &&
+      Assert(static_cast<unsigned int>(fe_degree) ==
+                 this->data->data.front().fe_degree &&
                n_q_points == this->n_quadrature_points,
              ExcMessage(message));
     }
@@ -6946,7 +6952,7 @@ FEEvaluation<dim,
   AssertIndexRange(q, n_q_points);
 
   const unsigned int n_q_points_1d_actual =
-    fe_degree == -1 ? this->data->n_q_points_1d : n_q_points_1d;
+    fe_degree == -1 ? this->data->data.front().n_q_points_1d : n_q_points_1d;
 
   // Cartesian mesh: not all quadrature points are stored, only the
   // diagonal. Hence, need to find the tensor product index and retrieve the
index 7cd7dd18c085a77b777a5481a83dbd53378f31a9..eec6794d263267ab82f32abc6bf64bf69c284e26 100644 (file)
@@ -1704,6 +1704,7 @@ MatrixFree<dim, Number, VectorizedArrayType>::initialize_indices(
               bool all_nodal = true;
               for (unsigned int c = 0; c < di.n_base_elements; ++c)
                 if (!shape_info(di.global_base_element_offset + c, 0, 0, 0)
+                       .data.front()
                        .nodal_at_cell_boundaries)
                   all_nodal = false;
               if (all_nodal == false)
index 5f93f4ed2f5a992dc9d37563412a0f89fde1ab7e..e1d17a200f9bf68a8b58bdefd9b0824b1b9b6cd7 100644 (file)
@@ -1048,7 +1048,7 @@ namespace MatrixFreeOperators
                                             fe_degree,
                                             n_components,
                                             VectorizedArrayType>::
-      apply(fe_eval.get_shape_info().inverse_shape_values_eo,
+      apply(fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
             inverse_coefficients,
             n_actual_components,
             in_array,
@@ -1077,7 +1077,7 @@ namespace MatrixFreeOperators
                                             n_components,
                                             VectorizedArrayType>::
       transform_from_q_points_to_basis(
-        fe_eval.get_shape_info().inverse_shape_values_eo,
+        fe_eval.get_shape_info().data.front().inverse_shape_values_eo,
         n_actual_components,
         in_array,
         out_array);
index afc38492bed80d06fd5ce3766895d16d8e9b3ab0..d143e4044177487cdf1eca36e160e756b0ae2485 100644 (file)
@@ -94,46 +94,9 @@ namespace internal
       tensor_symmetric_plus_dg0 = 5
     };
 
-    /**
-     * The class that stores the shape functions, gradients and Hessians
-     * evaluated for a tensor product finite element and tensor product
-     * quadrature formula on the unit cell. Because of this structure, only
-     * one-dimensional data is stored.
-     *
-     * @ingroup matrixfree
-     *
-     * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
-     */
     template <typename Number>
-    struct ShapeInfo
+    struct UnivariateShapeData
     {
-      /**
-       * Empty constructor. Does nothing.
-       */
-      ShapeInfo();
-
-      /**
-       * Constructor that initializes the data fields using the reinit method.
-       */
-      template <int dim>
-      ShapeInfo(const Quadrature<1> &     quad,
-                const FiniteElement<dim> &fe,
-                const unsigned int        base_element = 0);
-
-      /**
-       * Initializes the data fields. Takes a one-dimensional quadrature
-       * formula and a finite element as arguments and evaluates the shape
-       * functions, gradients and Hessians on the one-dimensional unit cell.
-       * This function assumes that the finite element is derived from a one-
-       * dimensional element by a tensor product and that the zeroth shape
-       * function in zero evaluates to one.
-       */
-      template <int dim>
-      void
-      reinit(const Quadrature<1> &     quad,
-             const FiniteElement<dim> &fe_dim,
-             const unsigned int        base_element = 0);
-
       /**
        * Return the memory consumption of this class in bytes.
        */
@@ -145,7 +108,7 @@ namespace internal
        * will select the most efficient algorithm based on the given element
        * type.
        */
-      ElementType element_type;
+      ElementType element_type = ElementType::tensor_general;
 
       /**
        * Stores the shape values of the 1D finite element evaluated on all 1D
@@ -275,6 +238,88 @@ namespace internal
        */
       Quadrature<1> quadrature;
 
+      /**
+       * Stores the degree of the element.
+       */
+      unsigned int fe_degree = 0;
+
+      /**
+       * Stores the number of quadrature points per dimension.
+       */
+      unsigned int n_q_points_1d = 0;
+
+      /**
+       * Indicates whether the basis functions are nodal in 0 and 1, i.e., the
+       * end points of the unit cell.
+       */
+      bool nodal_at_cell_boundaries = false;
+    };
+
+
+
+    /**
+     * The class that stores the shape functions, gradients and Hessians
+     * evaluated for a tensor product finite element and tensor product
+     * quadrature formula on the unit cell. Because of this structure, only
+     * one-dimensional data is stored.
+     *
+     * @ingroup matrixfree
+     *
+     * @author Katharina Kormann and Martin Kronbichler, 2010, 2011
+     */
+    template <typename Number>
+    struct ShapeInfo
+    {
+      /**
+       * Encodes the type of element detected at construction. FEEvaluation
+       * will select the most efficient algorithm based on the given element
+       * type.
+       */
+      ElementType element_type;
+
+      /**
+       * Empty constructor. Does nothing.
+       */
+      ShapeInfo();
+
+      /**
+       * Constructor that initializes the data fields using the reinit method.
+       */
+      template <int dim>
+      ShapeInfo(const Quadrature<1> &     quad,
+                const FiniteElement<dim> &fe,
+                const unsigned int        base_element = 0);
+
+      /**
+       * Initializes the data fields. Takes a one-dimensional quadrature
+       * formula and a finite element as arguments and evaluates the shape
+       * functions, gradients and Hessians on the one-dimensional unit cell.
+       * This function assumes that the finite element is derived from a one-
+       * dimensional element by a tensor product and that the zeroth shape
+       * function in zero evaluates to one.
+       */
+      template <int dim>
+      void
+      reinit(const Quadrature<1> &     quad,
+             const FiniteElement<dim> &fe_dim,
+             const unsigned int        base_element = 0);
+
+      /**
+       * Return data of univariate shape functions which defines the
+       * dimension @p dimension of tensor product shape functions
+       * regarding vector component @p component of the underlying
+       * finite element.
+       */
+      const UnivariateShapeData<Number> &
+      get_shape_data(const unsigned int dimension = 0,
+                     const unsigned int component = 0) const;
+
+      /**
+       * Return the memory consumption of this class in bytes.
+       */
+      std::size_t
+      memory_consumption() const;
+
       /**
        * Renumbering from deal.II's numbering of cell degrees of freedom to
        * lexicographic numbering used inside the FEEvaluation schemes of the
@@ -285,14 +330,28 @@ namespace internal
       std::vector<unsigned int> lexicographic_numbering;
 
       /**
-       * Stores the degree of the element.
+       * Stores data of univariate shape functions defining the
+       * underlying tensor product finite element.
        */
-      unsigned int fe_degree;
+      std::vector<UnivariateShapeData<Number>> data;
 
       /**
-       * Stores the number of quadrature points per dimension.
+       * Grants access to univariate shape function data of given
+       * dimension and vector component. Rows identify dimensions and
+       * columns identify vector components.
+       */
+      dealii::Table<2, UnivariateShapeData<Number> *> data_access;
+
+      /**
+       * Stores the number of space dimensions.
        */
-      unsigned int n_q_points_1d;
+      unsigned int n_dimensions;
+
+      /**
+       * Stores the number of vector components of the underlying
+       * vector-valued finite element.
+       */
+      unsigned int n_components;
 
       /**
        * Stores the number of quadrature points in @p dim dimensions for a
@@ -316,12 +375,6 @@ namespace internal
        */
       unsigned int dofs_per_component_on_face;
 
-      /**
-       * Indicates whether the basis functions are nodal in 0 and 1, i.e., the
-       * end points of the unit cell.
-       */
-      bool nodal_at_cell_boundaries;
-
       /**
        * For nodal basis functions with nodes located at the boundary of the
        * unit cell, face integrals that involve only the values of the shape
@@ -405,7 +458,8 @@ namespace internal
        * also fill the shape_???_eo fields.
        */
       bool
-      check_1d_shapes_symmetric(const unsigned int n_q_points_1d);
+      check_1d_shapes_symmetric(
+        UnivariateShapeData<Number> &univariate_shape_data);
 
       /**
        * Check whether symmetric 1D basis functions are such that the shape
@@ -414,7 +468,8 @@ namespace internal
        * that save some operations in the evaluation.
        */
       bool
-      check_1d_shapes_collocation() const;
+      check_1d_shapes_collocation(
+        const UnivariateShapeData<Number> &univariate_shape_data) const;
     };
 
 
@@ -427,17 +482,28 @@ namespace internal
                                         const FiniteElement<dim> &fe_in,
                                         const unsigned int base_element_number)
       : element_type(tensor_general)
-      , fe_degree(0)
-      , n_q_points_1d(0)
+      , n_dimensions(0)
+      , n_components(0)
       , n_q_points(0)
       , dofs_per_component_on_cell(0)
       , n_q_points_face(0)
       , dofs_per_component_on_face(0)
-      , nodal_at_cell_boundaries(false)
     {
       reinit(quad, fe_in, base_element_number);
     }
 
+    template <typename Number>
+    inline const UnivariateShapeData<Number> &
+    ShapeInfo<Number>::get_shape_data(const unsigned int dimension,
+                                      const unsigned int component) const
+    {
+      AssertDimension(n_dimensions, data_access.size(0));
+      AssertDimension(n_components, data_access.size(1));
+      AssertIndexRange(dimension, n_dimensions);
+      AssertIndexRange(component, n_components);
+      return *(data_access(dimension, component));
+    }
+
   } // end of namespace MatrixFreeFunctions
 
 } // end of namespace internal
index c89ba34f220bc3c471bcd66c48f190c75e2a6523..da01d690184ab0d06af2249f1398393937c529f1 100644 (file)
@@ -61,13 +61,12 @@ namespace internal
     template <typename Number>
     ShapeInfo<Number>::ShapeInfo()
       : element_type(tensor_general)
-      , fe_degree(numbers::invalid_unsigned_int)
-      , n_q_points_1d(0)
+      , n_dimensions(0)
+      , n_components(0)
       , n_q_points(0)
       , dofs_per_component_on_cell(0)
       , n_q_points_face(0)
       , dofs_per_component_on_face(0)
-      , nodal_at_cell_boundaries(false)
     {}
 
 
@@ -79,15 +78,43 @@ namespace internal
                               const FiniteElement<dim> &fe_in,
                               const unsigned int        base_element_number)
     {
-      this->quadrature             = quad;
       const FiniteElement<dim> *fe = &fe_in.base_element(base_element_number);
+      n_dimensions                 = dim;
+      n_components                 = fe_in.n_components();
 
       Assert(fe->n_components() == 1,
              ExcMessage("FEEvaluation only works for scalar finite elements."));
 
-      fe_degree     = fe->degree;
-      n_q_points_1d = quad.size();
-
+      // assuming isotropy of dimensions and components
+      data.resize(1);
+      UnivariateShapeData<Number> *univariate_shape_data = &data.front();
+      data_access.reinit(n_dimensions, n_components);
+      data_access.fill(univariate_shape_data);
+      univariate_shape_data->quadrature    = quad;
+      univariate_shape_data->fe_degree     = fe->degree;
+      univariate_shape_data->n_q_points_1d = quad.size();
+
+      // grant write access to common univariate shape data
+      auto &shape_values    = univariate_shape_data->shape_values;
+      auto &shape_gradients = univariate_shape_data->shape_gradients;
+      auto &shape_hessians  = univariate_shape_data->shape_hessians;
+      auto &shape_gradients_collocation =
+        univariate_shape_data->shape_gradients_collocation;
+      auto &shape_hessians_collocation =
+        univariate_shape_data->shape_hessians_collocation;
+      auto &inverse_shape_values = univariate_shape_data->inverse_shape_values;
+      auto &shape_data_on_face   = univariate_shape_data->shape_data_on_face;
+      auto &values_within_subface =
+        univariate_shape_data->values_within_subface;
+      auto &gradients_within_subface =
+        univariate_shape_data->gradients_within_subface;
+      auto &hessians_within_subface =
+        univariate_shape_data->hessians_within_subface;
+      auto &nodal_at_cell_boundaries =
+        univariate_shape_data->nodal_at_cell_boundaries;
+
+      const unsigned int fe_degree     = fe->degree;
+      const unsigned int n_q_points_1d = quad.size();
       const unsigned int n_dofs_1d = std::min(fe->dofs_per_cell, fe_degree + 1);
 
       // renumber (this is necessary for FE_Q, for example, since there the
@@ -205,18 +232,18 @@ namespace internal
         dim > 1 ? Utilities::fixed_power<dim - 1>(fe_degree + 1) : 1;
 
       const unsigned int array_size = n_dofs_1d * n_q_points_1d;
-      this->shape_gradients.resize_fast(array_size);
-      this->shape_values.resize_fast(array_size);
-      this->shape_hessians.resize_fast(array_size);
-
-      this->shape_data_on_face[0].resize(3 * n_dofs_1d);
-      this->shape_data_on_face[1].resize(3 * n_dofs_1d);
-      this->values_within_subface[0].resize(array_size);
-      this->values_within_subface[1].resize(array_size);
-      this->gradients_within_subface[0].resize(array_size);
-      this->gradients_within_subface[1].resize(array_size);
-      this->hessians_within_subface[0].resize(array_size);
-      this->hessians_within_subface[1].resize(array_size);
+      shape_gradients.resize_fast(array_size);
+      shape_values.resize_fast(array_size);
+      shape_hessians.resize_fast(array_size);
+
+      shape_data_on_face[0].resize(3 * n_dofs_1d);
+      shape_data_on_face[1].resize(3 * n_dofs_1d);
+      values_within_subface[0].resize(array_size);
+      values_within_subface[1].resize(array_size);
+      gradients_within_subface[0].resize(array_size);
+      gradients_within_subface[1].resize(array_size);
+      hessians_within_subface[0].resize(array_size);
+      hessians_within_subface[1].resize(array_size);
 
       for (unsigned int i = 0; i < n_dofs_1d; ++i)
         {
@@ -255,18 +282,18 @@ namespace internal
             }
 
           // evaluate basis functions on the 1D faces, i.e., in zero and one
-          Point<dim> q_point             = unit_point;
-          q_point[0]                     = 0;
-          this->shape_data_on_face[0][i] = fe->shape_value(my_i, q_point);
-          this->shape_data_on_face[0][i + n_dofs_1d] =
+          Point<dim> q_point       = unit_point;
+          q_point[0]               = 0;
+          shape_data_on_face[0][i] = fe->shape_value(my_i, q_point);
+          shape_data_on_face[0][i + n_dofs_1d] =
             fe->shape_grad(my_i, q_point)[0];
-          this->shape_data_on_face[0][i + 2 * n_dofs_1d] =
+          shape_data_on_face[0][i + 2 * n_dofs_1d] =
             fe->shape_grad_grad(my_i, q_point)[0][0];
-          q_point[0]                     = 1;
-          this->shape_data_on_face[1][i] = fe->shape_value(my_i, q_point);
-          this->shape_data_on_face[1][i + n_dofs_1d] =
+          q_point[0]               = 1;
+          shape_data_on_face[1][i] = fe->shape_value(my_i, q_point);
+          shape_data_on_face[1][i + n_dofs_1d] =
             fe->shape_grad(my_i, q_point)[0];
-          this->shape_data_on_face[1][i + 2 * n_dofs_1d] =
+          shape_data_on_face[1][i + 2 * n_dofs_1d] =
             fe->shape_grad_grad(my_i, q_point)[0][0];
         }
 
@@ -391,9 +418,9 @@ namespace internal
         }
 
       if (element_type == tensor_general &&
-          check_1d_shapes_symmetric(n_q_points_1d))
+          check_1d_shapes_symmetric(*univariate_shape_data))
         {
-          if (check_1d_shapes_collocation())
+          if (check_1d_shapes_collocation(*univariate_shape_data))
             element_type = tensor_symmetric_collocation;
           else
             element_type = tensor_symmetric;
@@ -412,7 +439,7 @@ namespace internal
             }
         }
       else if (element_type == tensor_symmetric_plus_dg0)
-        check_1d_shapes_symmetric(n_q_points_1d);
+        check_1d_shapes_symmetric(*univariate_shape_data);
 
       nodal_at_cell_boundaries = true;
       for (unsigned int i = 1; i < n_dofs_1d; ++i)
@@ -494,6 +521,9 @@ namespace internal
                     }
             }
         }
+
+      // TODO !!!
+      univariate_shape_data->element_type = this->element_type;
     }
 
 
@@ -501,11 +531,31 @@ namespace internal
     template <typename Number>
     bool
     ShapeInfo<Number>::check_1d_shapes_symmetric(
-      const unsigned int n_q_points_1d)
+      UnivariateShapeData<Number> &univariate_shape_data)
     {
       if (dofs_per_component_on_cell == 0)
         return false;
 
+      const auto n_q_points_1d   = univariate_shape_data.n_q_points_1d;
+      const auto fe_degree       = univariate_shape_data.fe_degree;
+      auto &     shape_values    = univariate_shape_data.shape_values;
+      auto &     shape_gradients = univariate_shape_data.shape_gradients;
+      auto &     shape_hessians  = univariate_shape_data.shape_hessians;
+      auto &     shape_gradients_collocation =
+        univariate_shape_data.shape_gradients_collocation;
+      auto &shape_hessians_collocation =
+        univariate_shape_data.shape_hessians_collocation;
+      auto &shape_values_eo    = univariate_shape_data.shape_values_eo;
+      auto &shape_gradients_eo = univariate_shape_data.shape_gradients_eo;
+      auto &shape_hessians_eo  = univariate_shape_data.shape_hessians_eo;
+      auto &shape_gradients_collocation_eo =
+        univariate_shape_data.shape_gradients_collocation_eo;
+      auto &shape_hessians_collocation_eo =
+        univariate_shape_data.shape_hessians_collocation_eo;
+      auto &inverse_shape_values = univariate_shape_data.inverse_shape_values;
+      auto &inverse_shape_values_eo =
+        univariate_shape_data.inverse_shape_values_eo;
+
       const double zero_tol =
         std::is_same<Number, double>::value == true ? 1e-12 : 1e-7;
       // symmetry for values
@@ -616,11 +666,15 @@ namespace internal
 
     template <typename Number>
     bool
-    ShapeInfo<Number>::check_1d_shapes_collocation() const
+    ShapeInfo<Number>::check_1d_shapes_collocation(
+      const UnivariateShapeData<Number> &univariate_shape_data) const
     {
       if (dofs_per_component_on_cell != n_q_points)
         return false;
 
+      const auto fe_degree    = univariate_shape_data.fe_degree;
+      auto &     shape_values = univariate_shape_data.shape_values;
+
       const double zero_tol =
         std::is_same<Number, double>::value == true ? 1e-12 : 1e-7;
       // check: identity operation for shape values
@@ -648,6 +702,16 @@ namespace internal
     template <typename Number>
     std::size_t
     ShapeInfo<Number>::memory_consumption() const
+    {
+      std::size_t memory = sizeof(*this);
+      for (const auto &univariate_shape_data : data)
+        memory += univariate_shape_data.memory_consumption();
+      return memory;
+    }
+
+    template <typename Number>
+    std::size_t
+    UnivariateShapeData<Number>::memory_consumption() const
     {
       std::size_t memory = sizeof(*this);
       memory += MemoryConsumption::memory_consumption(shape_values);
index be98ff8a3a20d083012ae420af3fe7d27e8634b1..34f62acfc60d44d392584bbc5fb5705be282c409 100644 (file)
@@ -43,19 +43,23 @@ test(const FiniteElement<dim> &fe,
   deallog << "Testing " << fe.get_name() << " with " << quadrature_name << "("
           << quad.size() << ")" << std::endl;
   deallog << "shape values: " << std::endl;
-  for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+  const auto &univariate_shape_data = shape_info.get_shape_data(0, 0);
+  for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
     {
       for (unsigned int q = 0; q < quad.size(); ++q)
-        deallog << std::setw(15) << shape_info.shape_values[i * quad.size() + q]
+        deallog << std::setw(15)
+                << univariate_shape_data.shape_values[i * quad.size() + q]
                 << " ";
       deallog << std::endl;
     }
   deallog << "inverse shape values: " << std::endl;
-  for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+  for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
     {
       for (unsigned int q = 0; q < quad.size(); ++q)
-        deallog << std::setw(15)
-                << shape_info.inverse_shape_values[i * quad.size() + q] << " ";
+        deallog
+          << std::setw(15)
+          << univariate_shape_data.inverse_shape_values[i * quad.size() + q]
+          << " ";
       deallog << std::endl;
     }
   deallog << "inverse shapes' * shapes: " << std::endl;
@@ -64,22 +68,24 @@ test(const FiniteElement<dim> &fe,
       for (unsigned int j = 0; j < quad.size(); ++j)
         {
           double sum = 0;
-          for (unsigned int k = 0; k < shape_info.fe_degree + 1; ++k)
-            sum += shape_info.inverse_shape_values[k * quad.size() + i] *
-                   shape_info.shape_values[k * quad.size() + j];
+          for (unsigned int k = 0; k < univariate_shape_data.fe_degree + 1; ++k)
+            sum +=
+              univariate_shape_data.inverse_shape_values[k * quad.size() + i] *
+              univariate_shape_data.shape_values[k * quad.size() + j];
           deallog << std::setw(15) << sum << " ";
         }
       deallog << std::endl;
     }
   deallog << "inverse shapes * shapes': " << std::endl;
-  for (unsigned int i = 0; i < shape_info.fe_degree + 1; ++i)
+  for (unsigned int i = 0; i < univariate_shape_data.fe_degree + 1; ++i)
     {
-      for (unsigned int j = 0; j < shape_info.fe_degree + 1; ++j)
+      for (unsigned int j = 0; j < univariate_shape_data.fe_degree + 1; ++j)
         {
           double sum = 0;
           for (unsigned int k = 0; k < quad.size(); ++k)
-            sum += shape_info.inverse_shape_values[i * quad.size() + k] *
-                   shape_info.shape_values[j * quad.size() + k];
+            sum +=
+              univariate_shape_data.inverse_shape_values[i * quad.size() + k] *
+              univariate_shape_data.shape_values[j * quad.size() + k];
           deallog << std::setw(15) << sum << " ";
         }
       deallog << std::endl;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.