]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Unify the different evaluate() functions into one templated function.
authorSebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Wed, 17 Apr 2024 15:17:36 +0000 (17:17 +0200)
committerSebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Wed, 17 Apr 2024 15:21:55 +0000 (17:21 +0200)
source/fe/fe_nedelec_sz.cc

index 5fe73762451a531c97766936a04454e79d204fd3..472b7446f8e7a53faa41ecd30eb849407074d303 100644 (file)
@@ -262,37 +262,25 @@ FE_NedelecSZ<dim, spacedim>::shape_grad_grad_component(
 template <int dim, int spacedim>
 void
 FE_NedelecSZ<dim, spacedim>::evaluate(
-  const std::vector<Point<dim>> /*p_list*/,
-  const UpdateFlags /*update_flags*/,
+  const std::vector<Point<dim>> &p_list,
+  const UpdateFlags              update_flags,
   std::unique_ptr<
-    typename dealii::FiniteElement<dim, spacedim>::InternalDataBase>
-    & /*data_ptr*/) const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-}
-
-
-
-template <>
-void
-FE_NedelecSZ<2, 2>::evaluate(
-  const std::vector<Point<2>> &p_list,
-  const UpdateFlags            update_flags,
-  std::unique_ptr<typename dealii::FiniteElement<2, 2>::InternalDataBase>
-    &data_ptr) const
+    typename dealii::FiniteElement<dim, spacedim>::InternalDataBase> &data_ptr)
+  const
 {
   auto &data       = dynamic_cast<InternalData &>(*data_ptr);
   data.update_each = requires_update_flags(update_flags);
 
   // Useful quantities:
-  const unsigned int dim = 2;
   const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
 
-  const unsigned int vertices_per_cell =
-    GeometryInfo<2 /*dim*/>::vertices_per_cell;
-  const unsigned int lines_per_cell = GeometryInfo<2 /*dim*/>::lines_per_cell;
+  const unsigned int vertices_per_cell = GeometryInfo<dim>::vertices_per_cell;
+  const unsigned int lines_per_cell    = GeometryInfo<dim>::lines_per_cell;
+  const unsigned int faces_per_cell    = GeometryInfo<dim>::faces_per_cell;
 
   const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
+  // we assume that all quads have the same number of dofs
+  const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
 
   const unsigned int n_q_points = p_list.size();
 
@@ -330,1104 +318,1128 @@ FE_NedelecSZ<2, 2>::evaluate(
   std::vector<std::vector<double>> lambda(n_q_points,
                                           std::vector<double>(lines_per_cell));
 
-  for (unsigned int q = 0; q < n_q_points; ++q)
+  switch (dim)
     {
-      sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
-      sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
-      sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
-      sigma[q][3] = p_list[q][0] + p_list[q][1];
-
-      lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
-      lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
-      lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
-      lambda[q][3] = p_list[q][0] * p_list[q][1];
-      for (unsigned int i = 0; i < vertices_per_cell; ++i)
-        for (unsigned int j = 0; j < vertices_per_cell; ++j)
-          data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
-    }
+      case 2:
+        {
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]);
+              sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]);
+              sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1];
+              sigma[q][3] = p_list[q][0] + p_list[q][1];
+
+              lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]);
+              lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]);
+              lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1];
+              lambda[q][3] = p_list[q][0] * p_list[q][1];
+              for (unsigned int i = 0; i < vertices_per_cell; ++i)
+                for (unsigned int j = 0; j < vertices_per_cell; ++j)
+                  data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+            }
 
-  // Calculate the gradient of sigma_imj_values[q][i][j] =
-  // sigma[q][i]-sigma[q][j]
-  //   - this depends on the component and the direction of the
-  //   corresponding edge.
-  //   - the direction of the edge is determined by
-  //   sigma_imj_sign[i][j].
-  // Helper arrays:
-  const int    sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1},
-                                                           {1, -1},
-                                                           {-1, 1},
-                                                           {1, 1}};
-  int          sigma_imj_sign[vertices_per_cell][vertices_per_cell];
-  unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
-
-  for (unsigned int i = 0; i < vertices_per_cell; ++i)
-    for (unsigned int j = 0; j < vertices_per_cell; ++j)
-      {
-        // sigma_imj_sign is the sign (+/-) of the coefficient of
-        // x/y/z in sigma_imj_values Due to the numbering of vertices
-        // on the reference element it is easy to find edges in the
-        // positive direction are from smaller to higher local vertex
-        // numbering.
-        sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
-        sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
-
-        // Now store the component which the sigma_i - sigma_j
-        // corresponds to:
-        sigma_imj_component[i][j] = 0;
-        for (unsigned int d = 0; d < dim; ++d)
-          {
-            int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
-            // Only interested in the first non-zero
-            // as if there is a second, it can not be a valid edge.
-            if (temp_imj != 0)
+          // Calculate the gradient of sigma_imj_values[q][i][j] =
+          // sigma[q][i]-sigma[q][j]
+          //   - this depends on the component and the direction of the
+          //   corresponding edge.
+          //   - the direction of the edge is determined by
+          //   sigma_imj_sign[i][j].
+          // Helper arrays:
+          const int sigma_comp_signs[GeometryInfo<2>::vertices_per_cell][2] = {
+            {-1, -1}, {1, -1}, {-1, 1}, {1, 1}};
+          int          sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+          unsigned int sigma_imj_component[vertices_per_cell]
+                                          [vertices_per_cell];
+
+          for (unsigned int i = 0; i < vertices_per_cell; ++i)
+            for (unsigned int j = 0; j < vertices_per_cell; ++j)
               {
-                sigma_imj_component[i][j] = d;
-                break;
+                // sigma_imj_sign is the sign (+/-) of the coefficient of
+                // x/y/z in sigma_imj_values Due to the numbering of vertices
+                // on the reference element it is easy to find edges in the
+                // positive direction are from smaller to higher local vertex
+                // numbering.
+                sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+                sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+                // Now store the component which the sigma_i - sigma_j
+                // corresponds to:
+                sigma_imj_component[i][j] = 0;
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    int temp_imj =
+                      sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+                    // Only interested in the first non-zero
+                    // as if there is a second, it can not be a valid edge.
+                    if (temp_imj != 0)
+                      {
+                        sigma_imj_component[i][j] = d;
+                        break;
+                      }
+                  }
+                // Can now calculate the gradient, only non-zero in the
+                // component given: Note some i,j combinations will be
+                // incorrect, but only on invalid edges.
+                data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
+                  2.0 * sigma_imj_sign[i][j];
               }
-          }
-        // Can now calculate the gradient, only non-zero in the
-        // component given: Note some i,j combinations will be
-        // incorrect, but only on invalid edges.
-        data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
-          2.0 * sigma_imj_sign[i][j];
-      }
 
-  // Now compute the edge parameterisations for a single element
-  // with global numbering matching that of the reference element:
+          // Now compute the edge parameterisations for a single element
+          // with global numbering matching that of the reference element:
 
-  // Resize the edge parameterisations
-  data.edge_sigma_values.resize(lines_per_cell,
-                                std::vector<double>(n_q_points));
-  data.edge_sigma_grads.resize(lines_per_cell, std::vector<double>(dim));
+          // Resize the edge parameterisations
+          data.edge_sigma_values.resize(lines_per_cell,
+                                        std::vector<double>(n_q_points));
+          data.edge_sigma_grads.resize(lines_per_cell,
+                                       std::vector<double>(dim));
 
-  // Fill the values for edge lambda and edge sigma:
-  const unsigned int edge_sigma_direction[lines_per_cell] = {1, 1, 0, 0};
+          // Fill the values for edge lambda and edge sigma:
+          const unsigned int
+            edge_sigma_direction[GeometryInfo<2>::lines_per_cell] = {1,
+                                                                     1,
+                                                                     0,
+                                                                     0};
 
-  data.edge_lambda_values.resize(lines_per_cell,
-                                 std::vector<double>(n_q_points));
+          data.edge_lambda_values.resize(lines_per_cell,
+                                         std::vector<double>(n_q_points));
 
-  data.edge_lambda_grads_2d.resize(lines_per_cell, std::vector<double>(dim));
+          data.edge_lambda_grads_2d.resize(lines_per_cell,
+                                           std::vector<double>(dim));
 
-  for (unsigned int m = 0; m < lines_per_cell; ++m)
-    {
-      // e1=max(reference vertex numbering on this edge)
-      // e2=min(reference vertex numbering on this edge)
-      // Which is guaranteed to be:
-      const unsigned int e1(GeometryInfo<dim>::line_to_cell_vertices(m, 1));
-      const unsigned int e2(GeometryInfo<dim>::line_to_cell_vertices(m, 0));
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          data.edge_sigma_values[m][q]  = data.sigma_imj_values[q][e2][e1];
-          data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
-        }
-
-      data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
-    }
-
-  data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
-  data.edge_lambda_grads_2d[1] = {1.0, 0.0};
-  data.edge_lambda_grads_2d[2] = {0.0, -1.0};
-  data.edge_lambda_grads_2d[3] = {0.0, 1.0};
+          for (unsigned int m = 0; m < lines_per_cell; ++m)
+            {
+              // e1=max(reference vertex numbering on this edge)
+              // e2=min(reference vertex numbering on this edge)
+              // Which is guaranteed to be:
+              const unsigned int e1(
+                GeometryInfo<dim>::line_to_cell_vertices(m, 1));
+              const unsigned int e2(
+                GeometryInfo<dim>::line_to_cell_vertices(m, 0));
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                {
+                  data.edge_sigma_values[m][q] =
+                    data.sigma_imj_values[q][e2][e1];
+                  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+                }
 
-  // If the polynomial order is 0, then no more work to do:
-  if (degree < 1)
-    return;
+              data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+            }
 
-  // Otherwise, we can compute the non-cell dependent shape functions.
-  //
-  // Note: the local dof numberings follow the usual order of lines ->
-  // faces -> cells
-  //       (we have no vertex-based DoFs in this element).
-  // For a given cell we have:
-  //      n_line_dofs = dofs_per_line*lines_per_cell.
-  //      n_face_dofs = dofs_per_face*faces_per_cell.
-  //      n_cell_dofs = dofs_per_quad (2d)
-  //                  = dofs_per_hex (3d)
-  //
-  // i.e. For the local dof numbering:
-  //      the first line dof is 0,
-  //      the first face dof is n_line_dofs,
-  //      the first cell dof is n_line_dofs + n_face_dofs.
-  //
-  // On a line, DoFs are ordered first by line_dof and then line_index:
-  // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
-  //
-  // and similarly for faces:
-  // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
-  //
-  // HOWEVER, we have different types of DoFs on a line/face/cell.
-  // On a line we have two types, lowest order and higher order
-  // gradients.
-  //    - The numbering is such the lowest order is first, then higher
-  //    order.
-  //      This is simple enough as there is only 1 lowest order and
-  //      degree higher orders DoFs per line.
-  //
-  // On a 2d cell, we have 3 types: Type 1/2/3:
-  //    - The ordering done by type:
-  //      - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
-  //        Numbered: ij1 = i1 + j1*(degree).        i.e. cell_dof_index
-  //        = ij1.
-  //      - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
-  //        Numbered: ij2 = i2 + j2*(degree).        i.e. cell_dof_index
-  //        = degree^2 + ij2
-  //      - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
-  //        Numbered: ij3 = i3.                      i.e. cell_dof_index
-  //        =  2*(degree^2) + ij3.
-  //
-  // These then fit into the local dof numbering described above:
-  // - local dof numberings are:
-  //   line_dofs: local_dof = line_dof_index.    0 <= local_dof <
-  //   dofs_per_line*lines_per_cell face_dofs: local_dof =
-  //   n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
-  //   = n_lines_dof + n_face_dofs + cell_dof_index.
-  //
-  // The cell-based shape functions are:
-  //
-  // Type 1 (gradients):
-  // \phi^{C_{1}}_{ij} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
-  //
-  // 0 <= i,j < degree.
-  //
-  // NOTE: The derivative produced by IntegratedLegendrePolynomials does
-  // not account for the
-  //       (2*x-1) or (2*y-1) so we must take this into account when
-  //       taking derivatives.
-  const unsigned int cell_type1_offset = n_line_dofs;
-
-  // Type 2:
-  // \phi^{C_{2}}_{ij} = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
-  //                     - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
-  //
-  // 0 <= i,j < degree.
-  const unsigned int cell_type2_offset = cell_type1_offset + degree * degree;
+          data.edge_lambda_grads_2d[0] = {-1.0, 0.0};
+          data.edge_lambda_grads_2d[1] = {1.0, 0.0};
+          data.edge_lambda_grads_2d[2] = {0.0, -1.0};
+          data.edge_lambda_grads_2d[3] = {0.0, 1.0};
 
-  // Type 3 (two subtypes):
-  // \phi^{C_{3}}_{j}        = L_{j+2}(2y-1) \mathbf{e}_{x}
-  //
-  // \phi^{C_{3}}_{j+degree} = L_{j+2}(2x-1) \mathbf{e}_{y}
-  //
-  // 0 <= j < degree
-  const unsigned int cell_type3_offset1 = cell_type2_offset + degree * degree;
-  const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
+          // If the polynomial order is 0, then no more work to do:
+          if (degree < 1)
+            break;
 
-  if (flags & (update_values | update_gradients | update_hessians))
-    {
-      // compute all points we must evaluate the 1d polynomials at:
-      std::vector<Point<dim>> cell_points(n_q_points);
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        for (unsigned int d = 0; d < dim; ++d)
-          cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
-
-      // Loop through quad points:
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          // pre-compute values & required derivatives at this quad
-          // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
+          // Otherwise, we can compute the non-cell dependent shape functions.
+          //
+          // Note: the local dof numberings follow the usual order of lines ->
+          // faces -> cells
+          //       (we have no vertex-based DoFs in this element).
+          // For a given cell we have:
+          //      n_line_dofs = dofs_per_line*lines_per_cell.
+          //      n_face_dofs = dofs_per_face*faces_per_cell.
+          //      n_cell_dofs = dofs_per_quad (2d)
+          //                  = dofs_per_hex (3d)
+          //
+          // i.e. For the local dof numbering:
+          //      the first line dof is 0,
+          //      the first face dof is n_line_dofs,
+          //      the first cell dof is n_line_dofs + n_face_dofs.
+          //
+          // On a line, DoFs are ordered first by line_dof and then line_index:
+          // i.e. line_dof_index = line_dof + line_index*(dofs_per_line)
+          //
+          // and similarly for faces:
+          // i.e. face_dof_index = face_dof + face_index*(dofs_per_face).
+          //
+          // HOWEVER, we have different types of DoFs on a line/face/cell.
+          // On a line we have two types, lowest order and higher order
+          // gradients.
+          //    - The numbering is such the lowest order is first, then higher
+          //    order.
+          //      This is simple enough as there is only 1 lowest order and
+          //      degree higher orders DoFs per line.
+          //
+          // On a 2d cell, we have 3 types: Type 1/2/3:
+          //    - The ordering done by type:
+          //      - Type 1: 0 <= i1,j1 < degree. degree^2 in total.
+          //        Numbered: ij1 = i1 + j1*(degree).        i.e. cell_dof_index
+          //        = ij1.
+          //      - Type 2: 0 <= i2,j2 < degree. degree^2 in total.
+          //        Numbered: ij2 = i2 + j2*(degree).        i.e. cell_dof_index
+          //        = degree^2 + ij2
+          //      - Type 3: 0 <= i3 < 2*degree. 2*degree in total.
+          //        Numbered: ij3 = i3.                      i.e. cell_dof_index
+          //        =  2*(degree^2) + ij3.
+          //
+          // These then fit into the local dof numbering described above:
+          // - local dof numberings are:
+          //   line_dofs: local_dof = line_dof_index.    0 <= local_dof <
+          //   dofs_per_line*lines_per_cell face_dofs: local_dof =
+          //   n_line_dofs*lines_per_cell + face_dof_index. cell dofs: local_dof
+          //   = n_lines_dof + n_face_dofs + cell_dof_index.
           //
-          // for each polyc[d], c=x,y, contains the d-th derivative with
-          // respect to the coordinate c.
+          // The cell-based shape functions are:
+          //
+          // Type 1 (gradients):
+          // \phi^{C_{1}}_{ij} = grad( L_{i+2}(2x-1)L_{j+2}(2y-1) ),
+          //
+          // 0 <= i,j < degree.
+          //
+          // NOTE: The derivative produced by IntegratedLegendrePolynomials does
+          // not account for the
+          //       (2*x-1) or (2*y-1) so we must take this into account when
+          //       taking derivatives.
+          const unsigned int cell_type1_offset = n_line_dofs;
+
+          // Type 2:
+          // \phi^{C_{2}}_{ij} = L'_{i+2}(2x-1) L_{j+2}(2y-1) \mathbf{e}_{x}
+          //                     - L_{i+2}(2x-1) L'_{j+2}(2y-1) \mathbf{e}_{y},
+          //
+          // 0 <= i,j < degree.
+          const unsigned int cell_type2_offset =
+            cell_type1_offset + degree * degree;
 
-          // We only need poly values and 1st derivative for
-          // update_values, but need the 2nd derivative too for
-          // update_gradients. For update_hessians we also need the 3rd
-          // derivatives.
-          const unsigned int poly_length =
-            (flags & update_hessians) ? 4 :
-                                        ((flags & update_gradients) ? 3 : 2);
+          // Type 3 (two subtypes):
+          // \phi^{C_{3}}_{j}        = L_{j+2}(2y-1) \mathbf{e}_{x}
+          //
+          // \phi^{C_{3}}_{j+degree} = L_{j+2}(2x-1) \mathbf{e}_{y}
+          //
+          // 0 <= j < degree
+          const unsigned int cell_type3_offset1 =
+            cell_type2_offset + degree * degree;
+          const unsigned int cell_type3_offset2 = cell_type3_offset1 + degree;
 
-          std::vector<std::vector<double>> polyx(
-            degree, std::vector<double>(poly_length));
-          std::vector<std::vector<double>> polyy(
-            degree, std::vector<double>(poly_length));
-          for (unsigned int i = 0; i < degree; ++i)
-            {
-              // Compute all required 1d polynomials and their
-              // derivatives, starting at degree 2. e.g. to access
-              // L'_{3}(2x-1) use polyx[1][1].
-              IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0],
-                                                         polyx[i]);
-              IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1],
-                                                         polyy[i]);
-            }
-          // Now use these to compute the shape functions:
-          if (flags & update_values)
-            {
-              for (unsigned int j = 0; j < degree; ++j)
-                {
-                  const unsigned int shift_j(j * degree);
-                  for (unsigned int i = 0; i < degree; ++i)
-                    {
-                      const unsigned int shift_ij(i + shift_j);
-
-                      // Type 1:
-                      const unsigned int dof_index1(cell_type1_offset +
-                                                    shift_ij);
-                      data.shape_values[dof_index1][q][0] =
-                        2.0 * polyx[i][1] * polyy[j][0];
-                      data.shape_values[dof_index1][q][1] =
-                        2.0 * polyx[i][0] * polyy[j][1];
-
-                      // Type 2:
-                      const unsigned int dof_index2(cell_type2_offset +
-                                                    shift_ij);
-                      data.shape_values[dof_index2][q][0] =
-                        data.shape_values[dof_index1][q][0];
-                      data.shape_values[dof_index2][q][1] =
-                        -1.0 * data.shape_values[dof_index1][q][1];
-                    }
-                  // Type 3:
-                  const unsigned int dof_index3_1(cell_type3_offset1 + j);
-                  data.shape_values[dof_index3_1][q][0] = polyy[j][0];
-                  data.shape_values[dof_index3_1][q][1] = 0.0;
-
-                  const unsigned int dof_index3_2(cell_type3_offset2 + j);
-                  data.shape_values[dof_index3_2][q][0] = 0.0;
-                  data.shape_values[dof_index3_2][q][1] = polyx[j][0];
-                }
-            }
-          if (flags & update_gradients)
+          if (flags & (update_values | update_gradients | update_hessians))
             {
-              for (unsigned int j = 0; j < degree; ++j)
+              // compute all points we must evaluate the 1d polynomials at:
+              std::vector<Point<dim>> cell_points(n_q_points);
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                for (unsigned int d = 0; d < dim; ++d)
+                  cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
+
+              // Loop through quad points:
+              for (unsigned int q = 0; q < n_q_points; ++q)
                 {
-                  const unsigned int shift_j(j * degree);
+                  // pre-compute values & required derivatives at this quad
+                  // point (x,y): polyx = L_{i+2}(2x-1), polyy = L_{j+2}(2y-1),
+                  //
+                  // for each polyc[d], c=x,y, contains the d-th derivative with
+                  // respect to the coordinate c.
+
+                  // We only need poly values and 1st derivative for
+                  // update_values, but need the 2nd derivative too for
+                  // update_gradients. For update_hessians we also need the 3rd
+                  // derivatives.
+                  const unsigned int poly_length =
+                    (flags & update_hessians) ?
+                      4 :
+                      ((flags & update_gradients) ? 3 : 2);
+
+                  std::vector<std::vector<double>> polyx(
+                    degree, std::vector<double>(poly_length));
+                  std::vector<std::vector<double>> polyy(
+                    degree, std::vector<double>(poly_length));
                   for (unsigned int i = 0; i < degree; ++i)
                     {
-                      const unsigned int shift_ij(i + shift_j);
-
-                      // Type 1:
-                      const unsigned int dof_index1(cell_type1_offset +
-                                                    shift_ij);
-                      data.shape_grads[dof_index1][q][0][0] =
-                        4.0 * polyx[i][2] * polyy[j][0];
-                      data.shape_grads[dof_index1][q][0][1] =
-                        4.0 * polyx[i][1] * polyy[j][1];
-                      data.shape_grads[dof_index1][q][1][0] =
-                        data.shape_grads[dof_index1][q][0][1];
-                      data.shape_grads[dof_index1][q][1][1] =
-                        4.0 * polyx[i][0] * polyy[j][2];
-
-                      // Type 2:
-                      const unsigned int dof_index2(cell_type2_offset +
-                                                    shift_ij);
-                      data.shape_grads[dof_index2][q][0][0] =
-                        data.shape_grads[dof_index1][q][0][0];
-                      data.shape_grads[dof_index2][q][0][1] =
-                        data.shape_grads[dof_index1][q][0][1];
-                      data.shape_grads[dof_index2][q][1][0] =
-                        -1.0 * data.shape_grads[dof_index1][q][1][0];
-                      data.shape_grads[dof_index2][q][1][1] =
-                        -1.0 * data.shape_grads[dof_index1][q][1][1];
+                      // Compute all required 1d polynomials and their
+                      // derivatives, starting at degree 2. e.g. to access
+                      // L'_{3}(2x-1) use polyx[1][1].
+                      IntegratedLegendrePolynomials[i + 2].value(
+                        cell_points[q][0], polyx[i]);
+                      IntegratedLegendrePolynomials[i + 2].value(
+                        cell_points[q][1], polyy[i]);
                     }
-                  // Type 3:
-                  const unsigned int dof_index3_1(cell_type3_offset1 + j);
-                  data.shape_grads[dof_index3_1][q][0][0] = 0.0;
-                  data.shape_grads[dof_index3_1][q][0][1] = 2.0 * polyy[j][1];
-                  data.shape_grads[dof_index3_1][q][1][0] = 0.0;
-                  data.shape_grads[dof_index3_1][q][1][1] = 0.0;
-
-                  const unsigned int dof_index3_2(cell_type3_offset2 + j);
-                  data.shape_grads[dof_index3_2][q][0][0] = 0.0;
-                  data.shape_grads[dof_index3_2][q][0][1] = 0.0;
-                  data.shape_grads[dof_index3_2][q][1][0] = 2.0 * polyx[j][1];
-                  data.shape_grads[dof_index3_2][q][1][1] = 0.0;
-                }
-            }
-          if (flags & update_hessians)
-            {
-              for (unsigned int j = 0; j < degree; ++j)
-                {
-                  const unsigned int shift_j(j * degree);
-                  for (unsigned int i = 0; i < degree; ++i)
+                  // Now use these to compute the shape functions:
+                  if (flags & update_values)
                     {
-                      const unsigned int shift_ij(i + shift_j);
-
-                      // Type 1:
-                      const unsigned int dof_index1(cell_type1_offset +
-                                                    shift_ij);
-                      data.shape_hessians[dof_index1][q][0][0][0] =
-                        8.0 * polyx[i][3] * polyy[j][0];
-                      data.shape_hessians[dof_index1][q][1][0][0] =
-                        8.0 * polyx[i][2] * polyy[j][1];
-
-                      data.shape_hessians[dof_index1][q][0][1][0] =
-                        data.shape_hessians[dof_index1][q][1][0][0];
-                      data.shape_hessians[dof_index1][q][1][1][0] =
-                        8.0 * polyx[i][1] * polyy[j][2];
-
-                      data.shape_hessians[dof_index1][q][0][0][1] =
-                        data.shape_hessians[dof_index1][q][1][0][0];
-                      data.shape_hessians[dof_index1][q][1][0][1] =
-                        data.shape_hessians[dof_index1][q][1][1][0];
-
-                      data.shape_hessians[dof_index1][q][0][1][1] =
-                        data.shape_hessians[dof_index1][q][1][1][0];
-                      data.shape_hessians[dof_index1][q][1][1][1] =
-                        8.0 * polyx[i][0] * polyy[j][3];
+                      for (unsigned int j = 0; j < degree; ++j)
+                        {
+                          const unsigned int shift_j(j * degree);
+                          for (unsigned int i = 0; i < degree; ++i)
+                            {
+                              const unsigned int shift_ij(i + shift_j);
 
+                              // Type 1:
+                              const unsigned int dof_index1(cell_type1_offset +
+                                                            shift_ij);
+                              data.shape_values[dof_index1][q][0] =
+                                2.0 * polyx[i][1] * polyy[j][0];
+                              data.shape_values[dof_index1][q][1] =
+                                2.0 * polyx[i][0] * polyy[j][1];
 
+                              // Type 2:
+                              const unsigned int dof_index2(cell_type2_offset +
+                                                            shift_ij);
+                              data.shape_values[dof_index2][q][0] =
+                                data.shape_values[dof_index1][q][0];
+                              data.shape_values[dof_index2][q][1] =
+                                -1.0 * data.shape_values[dof_index1][q][1];
+                            }
+                          // Type 3:
+                          const unsigned int dof_index3_1(cell_type3_offset1 +
+                                                          j);
+                          data.shape_values[dof_index3_1][q][0] = polyy[j][0];
+                          data.shape_values[dof_index3_1][q][1] = 0.0;
 
-                      // Type 2:
-                      const unsigned int dof_index2(cell_type2_offset +
-                                                    shift_ij);
-                      for (unsigned int d = 0; d < dim; ++d)
-                        {
-                          data.shape_hessians[dof_index2][q][0][0][d] =
-                            data.shape_hessians[dof_index1][q][0][0][d];
-                          data.shape_hessians[dof_index2][q][0][1][d] =
-                            data.shape_hessians[dof_index1][q][0][1][d];
-                          data.shape_hessians[dof_index2][q][1][0][d] =
-                            -1.0 * data.shape_hessians[dof_index1][q][1][0][d];
-                          data.shape_hessians[dof_index2][q][1][1][d] =
-                            -1.0 * data.shape_hessians[dof_index1][q][1][1][d];
+                          const unsigned int dof_index3_2(cell_type3_offset2 +
+                                                          j);
+                          data.shape_values[dof_index3_2][q][0] = 0.0;
+                          data.shape_values[dof_index3_2][q][1] = polyx[j][0];
                         }
                     }
-                  // Type 3:
-                  const unsigned int dof_index3_1(cell_type3_offset1 + j);
-                  data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][0][1][1] =
-                    4.0 * polyy[j][2];
-                  data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
-                  data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
-
-                  const unsigned int dof_index3_2(cell_type3_offset2 + j);
-                  data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][1][0][0] =
-                    4.0 * polyx[j][2];
-                  data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
-                  data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
-                }
-            }
-        }
-    }
-}
-
-
-
-template <>
-void
-FE_NedelecSZ<3, 3>::evaluate(
-  const std::vector<Point<3>> p_list,
-  const UpdateFlags           update_flags,
-  std::unique_ptr<typename dealii::FiniteElement<3, 3>::InternalDataBase>
-    &data_ptr) const
-{
-  auto &data       = dynamic_cast<InternalData &>(*data_ptr);
-  data.update_each = requires_update_flags(update_flags);
+                  if (flags & update_gradients)
+                    {
+                      for (unsigned int j = 0; j < degree; ++j)
+                        {
+                          const unsigned int shift_j(j * degree);
+                          for (unsigned int i = 0; i < degree; ++i)
+                            {
+                              const unsigned int shift_ij(i + shift_j);
 
-  // Useful quantities:
-  const unsigned int dim = 3;
-  const unsigned int degree(this->degree - 1); // Note: FE holds input degree+1
+                              // Type 1:
+                              const unsigned int dof_index1(cell_type1_offset +
+                                                            shift_ij);
+                              data.shape_grads[dof_index1][q][0][0] =
+                                4.0 * polyx[i][2] * polyy[j][0];
+                              data.shape_grads[dof_index1][q][0][1] =
+                                4.0 * polyx[i][1] * polyy[j][1];
+                              data.shape_grads[dof_index1][q][1][0] =
+                                data.shape_grads[dof_index1][q][0][1];
+                              data.shape_grads[dof_index1][q][1][1] =
+                                4.0 * polyx[i][0] * polyy[j][2];
 
-  const unsigned int vertices_per_cell =
-    GeometryInfo<3 /*dim*/>::vertices_per_cell;
-  const unsigned int lines_per_cell = GeometryInfo<3 /*dim*/>::lines_per_cell;
-  const unsigned int faces_per_cell = GeometryInfo<3 /*dim*/>::faces_per_cell;
+                              // Type 2:
+                              const unsigned int dof_index2(cell_type2_offset +
+                                                            shift_ij);
+                              data.shape_grads[dof_index2][q][0][0] =
+                                data.shape_grads[dof_index1][q][0][0];
+                              data.shape_grads[dof_index2][q][0][1] =
+                                data.shape_grads[dof_index1][q][0][1];
+                              data.shape_grads[dof_index2][q][1][0] =
+                                -1.0 * data.shape_grads[dof_index1][q][1][0];
+                              data.shape_grads[dof_index2][q][1][1] =
+                                -1.0 * data.shape_grads[dof_index1][q][1][1];
+                            }
+                          // Type 3:
+                          const unsigned int dof_index3_1(cell_type3_offset1 +
+                                                          j);
+                          data.shape_grads[dof_index3_1][q][0][0] = 0.0;
+                          data.shape_grads[dof_index3_1][q][0][1] =
+                            2.0 * polyy[j][1];
+                          data.shape_grads[dof_index3_1][q][1][0] = 0.0;
+                          data.shape_grads[dof_index3_1][q][1][1] = 0.0;
 
-  const unsigned int n_line_dofs = this->n_dofs_per_line() * lines_per_cell;
+                          const unsigned int dof_index3_2(cell_type3_offset2 +
+                                                          j);
+                          data.shape_grads[dof_index3_2][q][0][0] = 0.0;
+                          data.shape_grads[dof_index3_2][q][0][1] = 0.0;
+                          data.shape_grads[dof_index3_2][q][1][0] =
+                            2.0 * polyx[j][1];
+                          data.shape_grads[dof_index3_2][q][1][1] = 0.0;
+                        }
+                    }
+                  if (flags & update_hessians)
+                    {
+                      for (unsigned int j = 0; j < degree; ++j)
+                        {
+                          const unsigned int shift_j(j * degree);
+                          for (unsigned int i = 0; i < degree; ++i)
+                            {
+                              const unsigned int shift_ij(i + shift_j);
 
-  // we assume that all quads have the same number of dofs
-  const unsigned int n_face_dofs = this->n_dofs_per_quad(0) * faces_per_cell;
-  const unsigned int n_q_points  = p_list.size();
+                              // Type 1:
+                              const unsigned int dof_index1(cell_type1_offset +
+                                                            shift_ij);
+                              data.shape_hessians[dof_index1][q][0][0][0] =
+                                8.0 * polyx[i][3] * polyy[j][0];
+                              data.shape_hessians[dof_index1][q][1][0][0] =
+                                8.0 * polyx[i][2] * polyy[j][1];
 
-  const UpdateFlags flags(data.update_each);
+                              data.shape_hessians[dof_index1][q][0][1][0] =
+                                data.shape_hessians[dof_index1][q][1][0][0];
+                              data.shape_hessians[dof_index1][q][1][1][0] =
+                                8.0 * polyx[i][1] * polyy[j][2];
 
-  // Resize the internal data storage:
-  data.sigma_imj_values.resize(
-    n_q_points,
-    std::vector<std::vector<double>>(vertices_per_cell,
-                                     std::vector<double>(vertices_per_cell)));
+                              data.shape_hessians[dof_index1][q][0][0][1] =
+                                data.shape_hessians[dof_index1][q][1][0][0];
+                              data.shape_hessians[dof_index1][q][1][0][1] =
+                                data.shape_hessians[dof_index1][q][1][1][0];
 
-  data.sigma_imj_grads.resize(vertices_per_cell,
-                              std::vector<std::vector<double>>(
-                                vertices_per_cell, std::vector<double>(dim)));
+                              data.shape_hessians[dof_index1][q][0][1][1] =
+                                data.shape_hessians[dof_index1][q][1][1][0];
+                              data.shape_hessians[dof_index1][q][1][1][1] =
+                                8.0 * polyx[i][0] * polyy[j][3];
 
-  // Resize shape function arrays according to update flags:
-  if (flags & update_values)
-    data.shape_values.resize(this->n_dofs_per_cell(),
-                             std::vector<Tensor<1, dim>>(n_q_points));
 
-  if (flags & update_gradients)
-    data.shape_grads.resize(this->n_dofs_per_cell(),
-                            std::vector<DerivativeForm<1, dim, dim>>(
-                              n_q_points));
 
-  if (flags & update_hessians)
-    data.shape_hessians.resize(this->n_dofs_per_cell(),
-                               std::vector<DerivativeForm<2, dim, dim>>(
-                                 n_q_points));
-
-  // Compute values of sigma & lambda and the sigma differences and
-  // lambda additions.
-  std::vector<std::vector<double>> sigma(n_q_points,
-                                         std::vector<double>(lines_per_cell));
-  std::vector<std::vector<double>> lambda(n_q_points,
-                                          std::vector<double>(lines_per_cell));
+                              // Type 2:
+                              const unsigned int dof_index2(cell_type2_offset +
+                                                            shift_ij);
+                              for (unsigned int d = 0; d < dim; ++d)
+                                {
+                                  data.shape_hessians[dof_index2][q][0][0][d] =
+                                    data.shape_hessians[dof_index1][q][0][0][d];
+                                  data.shape_hessians[dof_index2][q][0][1][d] =
+                                    data.shape_hessians[dof_index1][q][0][1][d];
+                                  data.shape_hessians[dof_index2][q][1][0][d] =
+                                    -1.0 *
+                                    data.shape_hessians[dof_index1][q][1][0][d];
+                                  data.shape_hessians[dof_index2][q][1][1][d] =
+                                    -1.0 *
+                                    data.shape_hessians[dof_index1][q][1][1][d];
+                                }
+                            }
+                          // Type 3:
+                          const unsigned int dof_index3_1(cell_type3_offset1 +
+                                                          j);
+                          data.shape_hessians[dof_index3_1][q][0][0][0] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][0][0][1] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][0][1][0] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][0][1][1] =
+                            4.0 * polyy[j][2];
+                          data.shape_hessians[dof_index3_1][q][1][0][0] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][1][0][1] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][1][1][0] = 0.0;
+                          data.shape_hessians[dof_index3_1][q][1][1][1] = 0.0;
+
+                          const unsigned int dof_index3_2(cell_type3_offset2 +
+                                                          j);
+                          data.shape_hessians[dof_index3_2][q][0][0][0] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][0][0][1] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][0][1][0] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][0][1][1] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][1][0][0] =
+                            4.0 * polyx[j][2];
+                          data.shape_hessians[dof_index3_2][q][1][0][1] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][1][1][0] = 0.0;
+                          data.shape_hessians[dof_index3_2][q][1][1][1] = 0.0;
+                        }
+                    }
+                }
+            }
+          break;
+        }
 
-  for (unsigned int q = 0; q < n_q_points; ++q)
-    {
-      sigma[q][0] =
-        (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
-      sigma[q][1] = p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
-      sigma[q][2] = (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
-      sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
-      sigma[q][4] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
-      sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
-      sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
-      sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
-
-      lambda[q][0] =
-        (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
-      lambda[q][1] = p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
-      lambda[q][2] = (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
-      lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
-      lambda[q][4] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
-      lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
-      lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
-      lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
-
-      // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
-      // and lambda_ipj = \lambda_{i} + \lambda_{j}.
-      for (unsigned int i = 0; i < vertices_per_cell; ++i)
-        for (unsigned int j = 0; j < vertices_per_cell; ++j)
-          data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
-    }
+      case 3:
+        {
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              sigma[q][0] = (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) +
+                            (1 - p_list[q][2]);
+              sigma[q][1] =
+                p_list[q][0] + (1.0 - p_list[q][1]) + (1 - p_list[q][2]);
+              sigma[q][2] =
+                (1.0 - p_list[q][0]) + p_list[q][1] + (1 - p_list[q][2]);
+              sigma[q][3] = p_list[q][0] + p_list[q][1] + (1 - p_list[q][2]);
+              sigma[q][4] =
+                (1.0 - p_list[q][0]) + (1.0 - p_list[q][1]) + p_list[q][2];
+              sigma[q][5] = p_list[q][0] + (1.0 - p_list[q][1]) + p_list[q][2];
+              sigma[q][6] = (1.0 - p_list[q][0]) + p_list[q][1] + p_list[q][2];
+              sigma[q][7] = p_list[q][0] + p_list[q][1] + p_list[q][2];
+
+              lambda[q][0] = (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) *
+                             (1.0 - p_list[q][2]);
+              lambda[q][1] =
+                p_list[q][0] * (1.0 - p_list[q][1]) * (1.0 - p_list[q][2]);
+              lambda[q][2] =
+                (1.0 - p_list[q][0]) * p_list[q][1] * (1.0 - p_list[q][2]);
+              lambda[q][3] = p_list[q][0] * p_list[q][1] * (1.0 - p_list[q][2]);
+              lambda[q][4] =
+                (1.0 - p_list[q][0]) * (1.0 - p_list[q][1]) * p_list[q][2];
+              lambda[q][5] = p_list[q][0] * (1.0 - p_list[q][1]) * p_list[q][2];
+              lambda[q][6] = (1.0 - p_list[q][0]) * p_list[q][1] * p_list[q][2];
+              lambda[q][7] = p_list[q][0] * p_list[q][1] * p_list[q][2];
+
+              // Compute values of sigma_imj = \sigma_{i} - \sigma_{j}
+              // and lambda_ipj = \lambda_{i} + \lambda_{j}.
+              for (unsigned int i = 0; i < vertices_per_cell; ++i)
+                for (unsigned int j = 0; j < vertices_per_cell; ++j)
+                  data.sigma_imj_values[q][i][j] = sigma[q][i] - sigma[q][j];
+            }
 
-  // We now want some additional information about
-  // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
-  // calculate values & derivatives of the shape functions we need to
-  // know:
-  // - The component the sigma_imj value corresponds to - this varies
-  // with i & j.
-  // - The gradient of the sigma_imj value
-  //   - this depends on the component and the direction of the
-  //   corresponding edge.
-  //   - the direction of the edge is determined by
-  //   sigma_imj_sign[i][j].
-  //
-  // Note that not every i,j combination is a valid edge (there are only
-  // 12 valid edges in 3d), but we compute them all as it simplifies
-  // things.
-
-  // store the sign of each component x, y, z in the sigma list.
-  // can use this to fill in the sigma_imj_component data.
-  const int sigma_comp_signs[vertices_per_cell][dim] = {{-1, -1, -1},
-                                                        {1, -1, -1},
-                                                        {-1, 1, -1},
-                                                        {1, 1, -1},
-                                                        {-1, -1, 1},
-                                                        {1, -1, 1},
-                                                        {-1, 1, 1},
-                                                        {1, 1, 1}};
-
-  int          sigma_imj_sign[vertices_per_cell][vertices_per_cell];
-  unsigned int sigma_imj_component[vertices_per_cell][vertices_per_cell];
-
-  for (unsigned int i = 0; i < vertices_per_cell; ++i)
-    for (unsigned int j = 0; j < vertices_per_cell; ++j)
-      {
-        // sigma_imj_sign is the sign (+/-) of the coefficient of
-        // x/y/z in sigma_imj. Due to the numbering of vertices on the
-        // reference element this is easy to work out because edges in
-        // the positive direction go from smaller to higher local
-        // vertex numbering.
-        sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
-        sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
-
-        // Now store the component which the sigma_i - sigma_j
-        // corresponds to:
-        sigma_imj_component[i][j] = 0;
-        for (unsigned int d = 0; d < dim; ++d)
-          {
-            int temp_imj = sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
-            // Only interested in the first non-zero
-            // as if there is a second, it will not be a valid edge.
-            if (temp_imj != 0)
+          // We now want some additional information about
+          // sigma_imj_values[q][i][j] = sigma[q][i]-sigma[q][j] In order to
+          // calculate values & derivatives of the shape functions we need to
+          // know:
+          // - The component the sigma_imj value corresponds to - this varies
+          // with i & j.
+          // - The gradient of the sigma_imj value
+          //   - this depends on the component and the direction of the
+          //   corresponding edge.
+          //   - the direction of the edge is determined by
+          //   sigma_imj_sign[i][j].
+          //
+          // Note that not every i,j combination is a valid edge (there are only
+          // 12 valid edges in 3d), but we compute them all as it simplifies
+          // things.
+
+          // store the sign of each component x, y, z in the sigma list.
+          // can use this to fill in the sigma_imj_component data.
+          const int sigma_comp_signs[GeometryInfo<3>::vertices_per_cell][3] = {
+            {-1, -1, -1},
+            {1, -1, -1},
+            {-1, 1, -1},
+            {1, 1, -1},
+            {-1, -1, 1},
+            {1, -1, 1},
+            {-1, 1, 1},
+            {1, 1, 1}};
+
+          int          sigma_imj_sign[vertices_per_cell][vertices_per_cell];
+          unsigned int sigma_imj_component[vertices_per_cell]
+                                          [vertices_per_cell];
+
+          for (unsigned int i = 0; i < vertices_per_cell; ++i)
+            for (unsigned int j = 0; j < vertices_per_cell; ++j)
               {
-                sigma_imj_component[i][j] = d;
-                break;
+                // sigma_imj_sign is the sign (+/-) of the coefficient of
+                // x/y/z in sigma_imj. Due to the numbering of vertices on the
+                // reference element this is easy to work out because edges in
+                // the positive direction go from smaller to higher local
+                // vertex numbering.
+                sigma_imj_sign[i][j] = (i < j) ? -1 : 1;
+                sigma_imj_sign[i][j] = (i == j) ? 0 : sigma_imj_sign[i][j];
+
+                // Now store the component which the sigma_i - sigma_j
+                // corresponds to:
+                sigma_imj_component[i][j] = 0;
+                for (unsigned int d = 0; d < dim; ++d)
+                  {
+                    int temp_imj =
+                      sigma_comp_signs[i][d] - sigma_comp_signs[j][d];
+                    // Only interested in the first non-zero
+                    // as if there is a second, it will not be a valid edge.
+                    if (temp_imj != 0)
+                      {
+                        sigma_imj_component[i][j] = d;
+                        break;
+                      }
+                  }
+                // Can now calculate the gradient, only non-zero in the
+                // component given: Note some i,j combinations will be
+                // incorrect, but only on invalid edges.
+                data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
+                  2.0 * sigma_imj_sign[i][j];
               }
-          }
-        // Can now calculate the gradient, only non-zero in the
-        // component given: Note some i,j combinations will be
-        // incorrect, but only on invalid edges.
-        data.sigma_imj_grads[i][j][sigma_imj_component[i][j]] =
-          2.0 * sigma_imj_sign[i][j];
-      }
-
-  // Now compute the edge parameterisations for a single element
-  // with global numbering matching that of the reference element:
-
-  // resize the edge parameterisations
-  data.edge_sigma_values.resize(lines_per_cell,
-                                std::vector<double>(n_q_points));
-  data.edge_lambda_values.resize(lines_per_cell,
-                                 std::vector<double>(n_q_points));
-  data.edge_sigma_grads.resize(lines_per_cell, std::vector<double>(dim));
-  data.edge_lambda_grads_3d.resize(
-    lines_per_cell,
-    std::vector<std::vector<double>>(n_q_points, std::vector<double>(dim)));
-  data.edge_lambda_gradgrads_3d.resize(
-    lines_per_cell,
-    std::vector<std::vector<double>>(dim, std::vector<double>(dim)));
-
-  // Fill the values:
-  const unsigned int edge_sigma_direction[lines_per_cell] = {
-    1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
-
-  for (unsigned int m = 0; m < lines_per_cell; ++m)
-    {
-      // e1=max(reference vertex numbering on this edge)
-      // e2=min(reference vertex numbering on this edge)
-      // Which is guaranteed to be:
-      const unsigned int e1(GeometryInfo<dim>::line_to_cell_vertices(m, 1));
-      const unsigned int e2(GeometryInfo<dim>::line_to_cell_vertices(m, 0));
-
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          data.edge_sigma_values[m][q]  = data.sigma_imj_values[q][e2][e1];
-          data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
-        }
 
-      data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
-    }
+          // Now compute the edge parameterisations for a single element
+          // with global numbering matching that of the reference element:
+
+          // resize the edge parameterisations
+          data.edge_sigma_values.resize(lines_per_cell,
+                                        std::vector<double>(n_q_points));
+          data.edge_lambda_values.resize(lines_per_cell,
+                                         std::vector<double>(n_q_points));
+          data.edge_sigma_grads.resize(lines_per_cell,
+                                       std::vector<double>(dim));
+          data.edge_lambda_grads_3d.resize(
+            lines_per_cell,
+            std::vector<std::vector<double>>(n_q_points,
+                                             std::vector<double>(dim)));
+          data.edge_lambda_gradgrads_3d.resize(
+            lines_per_cell,
+            std::vector<std::vector<double>>(dim, std::vector<double>(dim)));
+
+          // Fill the values:
+          const unsigned int
+            edge_sigma_direction[GeometryInfo<3>::lines_per_cell] = {
+              1, 1, 0, 0, 1, 1, 0, 0, 2, 2, 2, 2};
 
-  // edge_lambda_grads
-  for (unsigned int q = 0; q < n_q_points; ++q)
-    {
-      double x(p_list[q][0]);
-      double y(p_list[q][1]);
-      double z(p_list[q][2]);
-      data.edge_lambda_grads_3d[0][q]  = {z - 1.0, 0.0, x - 1.0};
-      data.edge_lambda_grads_3d[1][q]  = {1.0 - z, 0.0, -x};
-      data.edge_lambda_grads_3d[2][q]  = {0.0, z - 1.0, y - 1.0};
-      data.edge_lambda_grads_3d[3][q]  = {0.0, 1.0 - z, -y};
-      data.edge_lambda_grads_3d[4][q]  = {-z, 0.0, 1.0 - x};
-      data.edge_lambda_grads_3d[5][q]  = {z, 0.0, x};
-      data.edge_lambda_grads_3d[6][q]  = {0.0, -z, 1.0 - y};
-      data.edge_lambda_grads_3d[7][q]  = {0.0, z, y};
-      data.edge_lambda_grads_3d[8][q]  = {y - 1.0, x - 1.0, 0.0};
-      data.edge_lambda_grads_3d[9][q]  = {1.0 - y, -x, 0.0};
-      data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
-      data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
-    }
+          for (unsigned int m = 0; m < lines_per_cell; ++m)
+            {
+              // e1=max(reference vertex numbering on this edge)
+              // e2=min(reference vertex numbering on this edge)
+              // Which is guaranteed to be:
+              const unsigned int e1(
+                GeometryInfo<dim>::line_to_cell_vertices(m, 1));
+              const unsigned int e2(
+                GeometryInfo<dim>::line_to_cell_vertices(m, 0));
 
-  // edge_lambda gradgrads:
-  const int edge_lambda_sign[lines_per_cell] = {
-    1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1}; // sign of the 2nd derivative for
-                                               // each edge.
-
-  const unsigned int edge_lambda_directions[lines_per_cell][2] = {
-    {0, 2},
-    {0, 2},
-    {1, 2},
-    {1, 2},
-    {0, 2},
-    {0, 2},
-    {1, 2},
-    {1, 2},
-    {0, 1},
-    {0, 1},
-    {0, 1},
-    {0, 1}}; // component which edge_lambda[m] depends on.
-
-  for (unsigned int m = 0; m < lines_per_cell; ++m)
-    {
-      data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
-                                   [edge_lambda_directions[m][1]] =
-        edge_lambda_sign[m];
-      data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
-                                   [edge_lambda_directions[m][0]] =
-        edge_lambda_sign[m];
-    }
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                {
+                  data.edge_sigma_values[m][q] =
+                    data.sigma_imj_values[q][e2][e1];
+                  data.edge_lambda_values[m][q] = lambda[q][e1] + lambda[q][e2];
+                }
 
-  // If the polynomial order is 0, then no more work to do:
-  if (degree < 1)
-    return;
+              data.edge_sigma_grads[m][edge_sigma_direction[m]] = -2.0;
+            }
 
-  // resize required data:
-  data.face_lambda_values.resize(faces_per_cell,
-                                 std::vector<double>(n_q_points));
-  data.face_lambda_grads.resize(faces_per_cell, std::vector<double>(dim));
+          // edge_lambda_grads
+          for (unsigned int q = 0; q < n_q_points; ++q)
+            {
+              double x(p_list[q][0]);
+              double y(p_list[q][1]);
+              double z(p_list[q][2]);
+              data.edge_lambda_grads_3d[0][q]  = {z - 1.0, 0.0, x - 1.0};
+              data.edge_lambda_grads_3d[1][q]  = {1.0 - z, 0.0, -x};
+              data.edge_lambda_grads_3d[2][q]  = {0.0, z - 1.0, y - 1.0};
+              data.edge_lambda_grads_3d[3][q]  = {0.0, 1.0 - z, -y};
+              data.edge_lambda_grads_3d[4][q]  = {-z, 0.0, 1.0 - x};
+              data.edge_lambda_grads_3d[5][q]  = {z, 0.0, x};
+              data.edge_lambda_grads_3d[6][q]  = {0.0, -z, 1.0 - y};
+              data.edge_lambda_grads_3d[7][q]  = {0.0, z, y};
+              data.edge_lambda_grads_3d[8][q]  = {y - 1.0, x - 1.0, 0.0};
+              data.edge_lambda_grads_3d[9][q]  = {1.0 - y, -x, 0.0};
+              data.edge_lambda_grads_3d[10][q] = {-y, 1.0 - x, 0.0};
+              data.edge_lambda_grads_3d[11][q] = {y, x, 0.0};
+            }
 
-  // Fill in the values (these don't change between cells).
-  for (unsigned int q = 0; q < n_q_points; ++q)
-    {
-      double x(p_list[q][0]);
-      double y(p_list[q][1]);
-      double z(p_list[q][2]);
-      data.face_lambda_values[0][q] = 1.0 - x;
-      data.face_lambda_values[1][q] = x;
-      data.face_lambda_values[2][q] = 1.0 - y;
-      data.face_lambda_values[3][q] = y;
-      data.face_lambda_values[4][q] = 1.0 - z;
-      data.face_lambda_values[5][q] = z;
-    }
+          // edge_lambda gradgrads:
+          const int edge_lambda_sign[GeometryInfo<3>::lines_per_cell] = {
+            1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1}; // sign of the 2nd
+                                                       // derivative for each
+                                                       // edge.
 
-  // gradients are constant:
-  data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
-  data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
-  data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
-  data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
-  data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
-  data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
-
-  // for cell-based shape functions:
-  // these don't depend on the cell, so can precompute all here:
-  if (flags & (update_values | update_gradients | update_hessians))
-    {
-      // Cell-based shape functions:
-      //
-      // Type-1 (gradients):
-      // \phi^{C_{1}}_{ijk} = grad(
-      // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
-      //
-      // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
-      const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs);
-      // Type-2:
-      //
-      // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
-      // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
-      // -1)\phi^{C_{1}}_{ijk}
-      //
-      // 0 <= i,j,k < degree. (subtypes in groups of
-      // degree*degree*degree)
-      //
-      // here we order so that all of subtype 1 comes first, then
-      // subtype 2.
-      const unsigned int cell_type2_offset1(cell_type1_offset +
-                                            degree * degree * degree);
-      const unsigned int cell_type2_offset2(cell_type2_offset1 +
-                                            degree * degree * degree);
-      // Type-3
-      // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
-      // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
-      // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
-      //
-      // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
-      //
-      // again we order so we compute all of subtype 1 first, then
-      // subtype 2, etc.
-      const unsigned int cell_type3_offset1(cell_type2_offset2 +
-                                            degree * degree * degree);
-      const unsigned int cell_type3_offset2(cell_type3_offset1 +
-                                            degree * degree);
-      const unsigned int cell_type3_offset3(cell_type3_offset2 +
-                                            degree * degree);
-
-      // compute all points we must evaluate the 1d polynomials at:
-      std::vector<Point<dim>> cell_points(n_q_points);
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          for (unsigned int d = 0; d < dim; ++d)
+          const unsigned int
+            edge_lambda_directions[GeometryInfo<3>::lines_per_cell][2] = {
+              {0, 2},
+              {0, 2},
+              {1, 2},
+              {1, 2},
+              {0, 2},
+              {0, 2},
+              {1, 2},
+              {1, 2},
+              {0, 1},
+              {0, 1},
+              {0, 1},
+              {0, 1}}; // component which edge_lambda[m] depends on.
+
+          for (unsigned int m = 0; m < lines_per_cell; ++m)
             {
-              cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
+              data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][0]]
+                                           [edge_lambda_directions[m][1]] =
+                edge_lambda_sign[m];
+              data.edge_lambda_gradgrads_3d[m][edge_lambda_directions[m][1]]
+                                           [edge_lambda_directions[m][0]] =
+                edge_lambda_sign[m];
             }
-        }
 
-      // We only need poly values and 1st derivative for
-      // update_values, but need the 2nd derivative too for
-      // update_gradients. For update_hessians we also need 3rd
-      // derivative.
-      const unsigned int poly_length =
-        (flags & update_hessians) ? 4 : ((flags & update_gradients) ? 3 : 2);
+          // If the polynomial order is 0, then no more work to do:
+          if (degree < 1)
+            break;
 
-      // Loop through quad points:
-      for (unsigned int q = 0; q < n_q_points; ++q)
-        {
-          // pre-compute values & required derivatives at this quad
-          // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
-          // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
-          //
-          // for each polyc[d], c=x,y,z, contains the d-th
-          // derivative with respect to the coordinate c.
-          std::vector<std::vector<double>> polyx(
-            degree, std::vector<double>(poly_length));
-          std::vector<std::vector<double>> polyy(
-            degree, std::vector<double>(poly_length));
-          std::vector<std::vector<double>> polyz(
-            degree, std::vector<double>(poly_length));
-          for (unsigned int i = 0; i < degree; ++i)
+          // resize required data:
+          data.face_lambda_values.resize(faces_per_cell,
+                                         std::vector<double>(n_q_points));
+          data.face_lambda_grads.resize(faces_per_cell,
+                                        std::vector<double>(dim));
+
+          // Fill in the values (these don't change between cells).
+          for (unsigned int q = 0; q < n_q_points; ++q)
             {
-              // compute all required 1d polynomials for i
-              IntegratedLegendrePolynomials[i + 2].value(cell_points[q][0],
-                                                         polyx[i]);
-              IntegratedLegendrePolynomials[i + 2].value(cell_points[q][1],
-                                                         polyy[i]);
-              IntegratedLegendrePolynomials[i + 2].value(cell_points[q][2],
-                                                         polyz[i]);
+              double x(p_list[q][0]);
+              double y(p_list[q][1]);
+              double z(p_list[q][2]);
+              data.face_lambda_values[0][q] = 1.0 - x;
+              data.face_lambda_values[1][q] = x;
+              data.face_lambda_values[2][q] = 1.0 - y;
+              data.face_lambda_values[3][q] = y;
+              data.face_lambda_values[4][q] = 1.0 - z;
+              data.face_lambda_values[5][q] = z;
             }
-          // Now use these to compute the shape functions:
-          if (flags & update_values)
+
+          // gradients are constant:
+          data.face_lambda_grads[0] = {-1.0, 0.0, 0.0};
+          data.face_lambda_grads[1] = {1.0, 0.0, 0.0};
+          data.face_lambda_grads[2] = {0.0, -1.0, 0.0};
+          data.face_lambda_grads[3] = {0.0, 1.0, 0.0};
+          data.face_lambda_grads[4] = {0.0, 0.0, -1.0};
+          data.face_lambda_grads[5] = {0.0, 0.0, 1.0};
+
+          // for cell-based shape functions:
+          // these don't depend on the cell, so can precompute all here:
+          if (flags & (update_values | update_gradients | update_hessians))
             {
-              for (unsigned int k = 0; k < degree; ++k)
+              // Cell-based shape functions:
+              //
+              // Type-1 (gradients):
+              // \phi^{C_{1}}_{ijk} = grad(
+              // L_{i+2}(2x-1)L_{j+2}(2y-1)L_{k+2}(2z-1) ),
+              //
+              // 0 <= i,j,k < degree. (in a group of degree*degree*degree)
+              const unsigned int cell_type1_offset(n_line_dofs + n_face_dofs);
+              // Type-2:
+              //
+              // \phi^{C_{2}}_{ijk} = diag(1, -1, 1)\phi^{C_{1}}_{ijk}
+              // \phi^{C_{2}}_{ijk + p^3} = diag(1, -1,
+              // -1)\phi^{C_{1}}_{ijk}
+              //
+              // 0 <= i,j,k < degree. (subtypes in groups of
+              // degree*degree*degree)
+              //
+              // here we order so that all of subtype 1 comes first, then
+              // subtype 2.
+              const unsigned int cell_type2_offset1(cell_type1_offset +
+                                                    degree * degree * degree);
+              const unsigned int cell_type2_offset2(cell_type2_offset1 +
+                                                    degree * degree * degree);
+              // Type-3
+              // \phi^{C_{3}}_{jk} = L_{j+2}(2y-1)L_{k+2}(2z-1)e_{x}
+              // \phi^{C_{3}}_{ik} = L_{i+2}(2x-1)L_{k+2}(2z-1)e_{y}
+              // \phi^{C_{3}}_{ij} = L_{i+2}(2x-1)L_{j+2}(2y-1)e_{z}
+              //
+              // 0 <= i,j,k < degree. (subtypes in groups of degree*degree)
+              //
+              // again we order so we compute all of subtype 1 first, then
+              // subtype 2, etc.
+              const unsigned int cell_type3_offset1(cell_type2_offset2 +
+                                                    degree * degree * degree);
+              const unsigned int cell_type3_offset2(cell_type3_offset1 +
+                                                    degree * degree);
+              const unsigned int cell_type3_offset3(cell_type3_offset2 +
+                                                    degree * degree);
+
+              // compute all points we must evaluate the 1d polynomials at:
+              std::vector<Point<dim>> cell_points(n_q_points);
+              for (unsigned int q = 0; q < n_q_points; ++q)
                 {
-                  const unsigned int shift_k(k * degree * degree);
-                  const unsigned int shift_j(k *
-                                             degree); // Used below when subbing
-                                                      // k for j (type 3)
-                  for (unsigned int j = 0; j < degree; ++j)
+                  for (unsigned int d = 0; d < dim; ++d)
                     {
-                      const unsigned int shift_jk(j * degree + shift_k);
-                      for (unsigned int i = 0; i < degree; ++i)
-                        {
-                          const unsigned int shift_ijk(shift_jk + i);
-
-                          // Type 1:
-                          const unsigned int dof_index1(cell_type1_offset +
-                                                        shift_ijk);
-
-                          data.shape_values[dof_index1][q][0] =
-                            2.0 * polyx[i][1] * polyy[j][0] * polyz[k][0];
-                          data.shape_values[dof_index1][q][1] =
-                            2.0 * polyx[i][0] * polyy[j][1] * polyz[k][0];
-                          data.shape_values[dof_index1][q][2] =
-                            2.0 * polyx[i][0] * polyy[j][0] * polyz[k][1];
-
-                          // Type 2:
-                          const unsigned int dof_index2_1(cell_type2_offset1 +
-                                                          shift_ijk);
-                          const unsigned int dof_index2_2(cell_type2_offset2 +
-                                                          shift_ijk);
-
-                          data.shape_values[dof_index2_1][q][0] =
-                            data.shape_values[dof_index1][q][0];
-                          data.shape_values[dof_index2_1][q][1] =
-                            -1.0 * data.shape_values[dof_index1][q][1];
-                          data.shape_values[dof_index2_1][q][2] =
-                            data.shape_values[dof_index1][q][2];
-
-                          data.shape_values[dof_index2_2][q][0] =
-                            data.shape_values[dof_index1][q][0];
-                          data.shape_values[dof_index2_2][q][1] =
-                            -1.0 * data.shape_values[dof_index1][q][1];
-                          data.shape_values[dof_index2_2][q][2] =
-                            -1.0 * data.shape_values[dof_index1][q][2];
-                        }
-                      // Type 3: (note we re-use k and j for
-                      // convenience):
-                      const unsigned int shift_ij(j +
-                                                  shift_j); // here we've subbed
-                                                            // j for i, k for j.
-                      const unsigned int dof_index3_1(cell_type3_offset1 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_2(cell_type3_offset2 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_3(cell_type3_offset3 +
-                                                      shift_ij);
-
-                      data.shape_values[dof_index3_1][q][0] =
-                        polyy[j][0] * polyz[k][0];
-                      data.shape_values[dof_index3_1][q][1] = 0.0;
-                      data.shape_values[dof_index3_1][q][2] = 0.0;
-
-                      data.shape_values[dof_index3_2][q][0] = 0.0;
-                      data.shape_values[dof_index3_2][q][1] =
-                        polyx[j][0] * polyz[k][0];
-                      data.shape_values[dof_index3_2][q][2] = 0.0;
-
-                      data.shape_values[dof_index3_3][q][0] = 0.0;
-                      data.shape_values[dof_index3_3][q][1] = 0.0;
-                      data.shape_values[dof_index3_3][q][2] =
-                        polyx[j][0] * polyy[k][0];
+                      cell_points[q][d] = 2.0 * p_list[q][d] - 1.0;
                     }
                 }
-            }
-          if (flags & update_gradients)
-            {
-              for (unsigned int k = 0; k < degree; ++k)
+
+              // We only need poly values and 1st derivative for
+              // update_values, but need the 2nd derivative too for
+              // update_gradients. For update_hessians we also need 3rd
+              // derivative.
+              const unsigned int poly_length =
+                (flags & update_hessians) ?
+                  4 :
+                  ((flags & update_gradients) ? 3 : 2);
+
+              // Loop through quad points:
+              for (unsigned int q = 0; q < n_q_points; ++q)
                 {
-                  const unsigned int shift_k(k * degree * degree);
-                  const unsigned int shift_j(k *
-                                             degree); // Used below when subbing
-                                                      // k for j (type 3)
-                  for (unsigned int j = 0; j < degree; ++j)
+                  // pre-compute values & required derivatives at this quad
+                  // point, (x,y,z): polyx = L_{i+2}(2x-1), polyy =
+                  // L_{j+2}(2y-1), polyz = L_{k+2}(2z-1).
+                  //
+                  // for each polyc[d], c=x,y,z, contains the d-th
+                  // derivative with respect to the coordinate c.
+                  std::vector<std::vector<double>> polyx(
+                    degree, std::vector<double>(poly_length));
+                  std::vector<std::vector<double>> polyy(
+                    degree, std::vector<double>(poly_length));
+                  std::vector<std::vector<double>> polyz(
+                    degree, std::vector<double>(poly_length));
+                  for (unsigned int i = 0; i < degree; ++i)
                     {
-                      const unsigned int shift_jk(j * degree + shift_k);
-                      for (unsigned int i = 0; i < degree; ++i)
+                      // compute all required 1d polynomials for i
+                      IntegratedLegendrePolynomials[i + 2].value(
+                        cell_points[q][0], polyx[i]);
+                      IntegratedLegendrePolynomials[i + 2].value(
+                        cell_points[q][1], polyy[i]);
+                      IntegratedLegendrePolynomials[i + 2].value(
+                        cell_points[q][2], polyz[i]);
+                    }
+                  // Now use these to compute the shape functions:
+                  if (flags & update_values)
+                    {
+                      for (unsigned int k = 0; k < degree; ++k)
                         {
-                          const unsigned int shift_ijk(shift_jk + i);
-
-                          // Type 1:
-                          const unsigned int dof_index1(cell_type1_offset +
-                                                        shift_ijk);
-
-                          data.shape_grads[dof_index1][q][0][0] =
-                            4.0 * polyx[i][2] * polyy[j][0] * polyz[k][0];
-                          data.shape_grads[dof_index1][q][0][1] =
-                            4.0 * polyx[i][1] * polyy[j][1] * polyz[k][0];
-                          data.shape_grads[dof_index1][q][0][2] =
-                            4.0 * polyx[i][1] * polyy[j][0] * polyz[k][1];
-
-                          data.shape_grads[dof_index1][q][1][0] =
-                            data.shape_grads[dof_index1][q][0][1];
-                          data.shape_grads[dof_index1][q][1][1] =
-                            4.0 * polyx[i][0] * polyy[j][2] * polyz[k][0];
-                          data.shape_grads[dof_index1][q][1][2] =
-                            4.0 * polyx[i][0] * polyy[j][1] * polyz[k][1];
-
-                          data.shape_grads[dof_index1][q][2][0] =
-                            data.shape_grads[dof_index1][q][0][2];
-                          data.shape_grads[dof_index1][q][2][1] =
-                            data.shape_grads[dof_index1][q][1][2];
-                          data.shape_grads[dof_index1][q][2][2] =
-                            4.0 * polyx[i][0] * polyy[j][0] * polyz[k][2];
-
-                          // Type 2:
-                          const unsigned int dof_index2_1(cell_type2_offset1 +
-                                                          shift_ijk);
-                          const unsigned int dof_index2_2(cell_type2_offset2 +
-                                                          shift_ijk);
-
-                          for (unsigned int d = 0; d < dim; ++d)
+                          const unsigned int shift_k(k * degree * degree);
+                          const unsigned int shift_j(
+                            k * degree); // Used below when subbing
+                                         // k for j (type 3)
+                          for (unsigned int j = 0; j < degree; ++j)
                             {
-                              data.shape_grads[dof_index2_1][q][0][d] =
-                                data.shape_grads[dof_index1][q][0][d];
-                              data.shape_grads[dof_index2_1][q][1][d] =
-                                -1.0 * data.shape_grads[dof_index1][q][1][d];
-                              data.shape_grads[dof_index2_1][q][2][d] =
-                                data.shape_grads[dof_index1][q][2][d];
-
-                              data.shape_grads[dof_index2_2][q][0][d] =
-                                data.shape_grads[dof_index1][q][0][d];
-                              data.shape_grads[dof_index2_2][q][1][d] =
-                                -1.0 * data.shape_grads[dof_index1][q][1][d];
-                              data.shape_grads[dof_index2_2][q][2][d] =
-                                -1.0 * data.shape_grads[dof_index1][q][2][d];
+                              const unsigned int shift_jk(j * degree + shift_k);
+                              for (unsigned int i = 0; i < degree; ++i)
+                                {
+                                  const unsigned int shift_ijk(shift_jk + i);
+
+                                  // Type 1:
+                                  const unsigned int dof_index1(
+                                    cell_type1_offset + shift_ijk);
+
+                                  data.shape_values[dof_index1][q][0] =
+                                    2.0 * polyx[i][1] * polyy[j][0] *
+                                    polyz[k][0];
+                                  data.shape_values[dof_index1][q][1] =
+                                    2.0 * polyx[i][0] * polyy[j][1] *
+                                    polyz[k][0];
+                                  data.shape_values[dof_index1][q][2] =
+                                    2.0 * polyx[i][0] * polyy[j][0] *
+                                    polyz[k][1];
+
+                                  // Type 2:
+                                  const unsigned int dof_index2_1(
+                                    cell_type2_offset1 + shift_ijk);
+                                  const unsigned int dof_index2_2(
+                                    cell_type2_offset2 + shift_ijk);
+
+                                  data.shape_values[dof_index2_1][q][0] =
+                                    data.shape_values[dof_index1][q][0];
+                                  data.shape_values[dof_index2_1][q][1] =
+                                    -1.0 * data.shape_values[dof_index1][q][1];
+                                  data.shape_values[dof_index2_1][q][2] =
+                                    data.shape_values[dof_index1][q][2];
+
+                                  data.shape_values[dof_index2_2][q][0] =
+                                    data.shape_values[dof_index1][q][0];
+                                  data.shape_values[dof_index2_2][q][1] =
+                                    -1.0 * data.shape_values[dof_index1][q][1];
+                                  data.shape_values[dof_index2_2][q][2] =
+                                    -1.0 * data.shape_values[dof_index1][q][2];
+                                }
+                              // Type 3: (note we re-use k and j for
+                              // convenience):
+                              const unsigned int shift_ij(
+                                j + shift_j); // here we've subbed
+                                              // j for i, k for j.
+                              const unsigned int dof_index3_1(
+                                cell_type3_offset1 + shift_ij);
+                              const unsigned int dof_index3_2(
+                                cell_type3_offset2 + shift_ij);
+                              const unsigned int dof_index3_3(
+                                cell_type3_offset3 + shift_ij);
+
+                              data.shape_values[dof_index3_1][q][0] =
+                                polyy[j][0] * polyz[k][0];
+                              data.shape_values[dof_index3_1][q][1] = 0.0;
+                              data.shape_values[dof_index3_1][q][2] = 0.0;
+
+                              data.shape_values[dof_index3_2][q][0] = 0.0;
+                              data.shape_values[dof_index3_2][q][1] =
+                                polyx[j][0] * polyz[k][0];
+                              data.shape_values[dof_index3_2][q][2] = 0.0;
+
+                              data.shape_values[dof_index3_3][q][0] = 0.0;
+                              data.shape_values[dof_index3_3][q][1] = 0.0;
+                              data.shape_values[dof_index3_3][q][2] =
+                                polyx[j][0] * polyy[k][0];
                             }
                         }
-                      // Type 3: (note we re-use k and j for
-                      // convenience):
-                      const unsigned int shift_ij(j +
-                                                  shift_j); // here we've subbed
-                                                            // j for i, k for j.
-                      const unsigned int dof_index3_1(cell_type3_offset1 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_2(cell_type3_offset2 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_3(cell_type3_offset3 +
-                                                      shift_ij);
-                      for (unsigned int d1 = 0; d1 < dim; ++d1)
+                    }
+                  if (flags & update_gradients)
+                    {
+                      for (unsigned int k = 0; k < degree; ++k)
                         {
-                          for (unsigned int d2 = 0; d2 < dim; ++d2)
+                          const unsigned int shift_k(k * degree * degree);
+                          const unsigned int shift_j(
+                            k * degree); // Used below when subbing
+                                         // k for j (type 3)
+                          for (unsigned int j = 0; j < degree; ++j)
                             {
-                              data.shape_grads[dof_index3_1][q][d1][d2] = 0.0;
-                              data.shape_grads[dof_index3_2][q][d1][d2] = 0.0;
-                              data.shape_grads[dof_index3_3][q][d1][d2] = 0.0;
+                              const unsigned int shift_jk(j * degree + shift_k);
+                              for (unsigned int i = 0; i < degree; ++i)
+                                {
+                                  const unsigned int shift_ijk(shift_jk + i);
+
+                                  // Type 1:
+                                  const unsigned int dof_index1(
+                                    cell_type1_offset + shift_ijk);
+
+                                  data.shape_grads[dof_index1][q][0][0] =
+                                    4.0 * polyx[i][2] * polyy[j][0] *
+                                    polyz[k][0];
+                                  data.shape_grads[dof_index1][q][0][1] =
+                                    4.0 * polyx[i][1] * polyy[j][1] *
+                                    polyz[k][0];
+                                  data.shape_grads[dof_index1][q][0][2] =
+                                    4.0 * polyx[i][1] * polyy[j][0] *
+                                    polyz[k][1];
+
+                                  data.shape_grads[dof_index1][q][1][0] =
+                                    data.shape_grads[dof_index1][q][0][1];
+                                  data.shape_grads[dof_index1][q][1][1] =
+                                    4.0 * polyx[i][0] * polyy[j][2] *
+                                    polyz[k][0];
+                                  data.shape_grads[dof_index1][q][1][2] =
+                                    4.0 * polyx[i][0] * polyy[j][1] *
+                                    polyz[k][1];
+
+                                  data.shape_grads[dof_index1][q][2][0] =
+                                    data.shape_grads[dof_index1][q][0][2];
+                                  data.shape_grads[dof_index1][q][2][1] =
+                                    data.shape_grads[dof_index1][q][1][2];
+                                  data.shape_grads[dof_index1][q][2][2] =
+                                    4.0 * polyx[i][0] * polyy[j][0] *
+                                    polyz[k][2];
+
+                                  // Type 2:
+                                  const unsigned int dof_index2_1(
+                                    cell_type2_offset1 + shift_ijk);
+                                  const unsigned int dof_index2_2(
+                                    cell_type2_offset2 + shift_ijk);
+
+                                  for (unsigned int d = 0; d < dim; ++d)
+                                    {
+                                      data.shape_grads[dof_index2_1][q][0][d] =
+                                        data.shape_grads[dof_index1][q][0][d];
+                                      data.shape_grads[dof_index2_1][q][1][d] =
+                                        -1.0 *
+                                        data.shape_grads[dof_index1][q][1][d];
+                                      data.shape_grads[dof_index2_1][q][2][d] =
+                                        data.shape_grads[dof_index1][q][2][d];
+
+                                      data.shape_grads[dof_index2_2][q][0][d] =
+                                        data.shape_grads[dof_index1][q][0][d];
+                                      data.shape_grads[dof_index2_2][q][1][d] =
+                                        -1.0 *
+                                        data.shape_grads[dof_index1][q][1][d];
+                                      data.shape_grads[dof_index2_2][q][2][d] =
+                                        -1.0 *
+                                        data.shape_grads[dof_index1][q][2][d];
+                                    }
+                                }
+                              // Type 3: (note we re-use k and j for
+                              // convenience):
+                              const unsigned int shift_ij(
+                                j + shift_j); // here we've subbed
+                                              // j for i, k for j.
+                              const unsigned int dof_index3_1(
+                                cell_type3_offset1 + shift_ij);
+                              const unsigned int dof_index3_2(
+                                cell_type3_offset2 + shift_ij);
+                              const unsigned int dof_index3_3(
+                                cell_type3_offset3 + shift_ij);
+                              for (unsigned int d1 = 0; d1 < dim; ++d1)
+                                {
+                                  for (unsigned int d2 = 0; d2 < dim; ++d2)
+                                    {
+                                      data
+                                        .shape_grads[dof_index3_1][q][d1][d2] =
+                                        0.0;
+                                      data
+                                        .shape_grads[dof_index3_2][q][d1][d2] =
+                                        0.0;
+                                      data
+                                        .shape_grads[dof_index3_3][q][d1][d2] =
+                                        0.0;
+                                    }
+                                }
+                              data.shape_grads[dof_index3_1][q][0][1] =
+                                2.0 * polyy[j][1] * polyz[k][0];
+                              data.shape_grads[dof_index3_1][q][0][2] =
+                                2.0 * polyy[j][0] * polyz[k][1];
+
+                              data.shape_grads[dof_index3_2][q][1][0] =
+                                2.0 * polyx[j][1] * polyz[k][0];
+                              data.shape_grads[dof_index3_2][q][1][2] =
+                                2.0 * polyx[j][0] * polyz[k][1];
+
+                              data.shape_grads[dof_index3_3][q][2][0] =
+                                2.0 * polyx[j][1] * polyy[k][0];
+                              data.shape_grads[dof_index3_3][q][2][1] =
+                                2.0 * polyx[j][0] * polyy[k][1];
                             }
                         }
-                      data.shape_grads[dof_index3_1][q][0][1] =
-                        2.0 * polyy[j][1] * polyz[k][0];
-                      data.shape_grads[dof_index3_1][q][0][2] =
-                        2.0 * polyy[j][0] * polyz[k][1];
-
-                      data.shape_grads[dof_index3_2][q][1][0] =
-                        2.0 * polyx[j][1] * polyz[k][0];
-                      data.shape_grads[dof_index3_2][q][1][2] =
-                        2.0 * polyx[j][0] * polyz[k][1];
-
-                      data.shape_grads[dof_index3_3][q][2][0] =
-                        2.0 * polyx[j][1] * polyy[k][0];
-                      data.shape_grads[dof_index3_3][q][2][1] =
-                        2.0 * polyx[j][0] * polyy[k][1];
                     }
-                }
-            }
-          if (flags & update_hessians)
-            {
-              for (unsigned int k = 0; k < degree; ++k)
-                {
-                  const unsigned int shift_k(k * degree * degree);
-                  const unsigned int shift_j(k *
-                                             degree); // Used below when subbing
-                                                      // k for j type 3
-
-                  for (unsigned int j = 0; j < degree; ++j)
+                  if (flags & update_hessians)
                     {
-                      const unsigned int shift_jk(j * degree + shift_k);
-                      for (unsigned int i = 0; i < degree; ++i)
+                      for (unsigned int k = 0; k < degree; ++k)
                         {
-                          const unsigned int shift_ijk(shift_jk + i);
-
-                          // Type 1:
-                          const unsigned int dof_index1(cell_type1_offset +
-                                                        shift_ijk);
-
-                          data.shape_hessians[dof_index1][q][0][0][0] =
-                            8.0 * polyx[i][3] * polyy[j][0] * polyz[k][0];
-                          data.shape_hessians[dof_index1][q][1][0][0] =
-                            8.0 * polyx[i][2] * polyy[j][1] * polyz[k][0];
-                          data.shape_hessians[dof_index1][q][2][0][0] =
-                            8.0 * polyx[i][2] * polyy[j][0] * polyz[k][1];
-
-                          data.shape_hessians[dof_index1][q][0][1][0] =
-                            data.shape_hessians[dof_index1][q][1][0][0];
-                          data.shape_hessians[dof_index1][q][1][1][0] =
-                            8.0 * polyx[i][1] * polyy[j][2] * polyz[k][0];
-                          data.shape_hessians[dof_index1][q][2][1][0] =
-                            8.0 * polyx[i][1] * polyy[j][1] * polyz[k][1];
-
-                          data.shape_hessians[dof_index1][q][0][2][0] =
-                            data.shape_hessians[dof_index1][q][2][0][0];
-                          data.shape_hessians[dof_index1][q][1][2][0] =
-                            data.shape_hessians[dof_index1][q][2][1][0];
-                          data.shape_hessians[dof_index1][q][2][2][0] =
-                            8.0 * polyx[i][1] * polyy[j][0] * polyz[k][2];
-
-
-                          data.shape_hessians[dof_index1][q][0][0][1] =
-                            data.shape_hessians[dof_index1][q][1][0][0];
-                          data.shape_hessians[dof_index1][q][1][0][1] =
-                            data.shape_hessians[dof_index1][q][1][1][0];
-                          data.shape_hessians[dof_index1][q][2][0][1] =
-                            data.shape_hessians[dof_index1][q][2][1][0];
-
-                          data.shape_hessians[dof_index1][q][0][1][1] =
-                            data.shape_hessians[dof_index1][q][1][1][0];
-                          data.shape_hessians[dof_index1][q][1][1][1] =
-                            8.0 * polyx[i][0] * polyy[j][3] * polyz[k][0];
-                          data.shape_hessians[dof_index1][q][2][1][1] =
-                            8.0 * polyx[i][0] * polyy[j][2] * polyz[k][1];
-
-                          data.shape_hessians[dof_index1][q][0][2][1] =
-                            data.shape_hessians[dof_index1][q][2][1][0];
-                          data.shape_hessians[dof_index1][q][1][2][1] =
-                            data.shape_hessians[dof_index1][q][2][1][1];
-                          data.shape_hessians[dof_index1][q][2][2][1] =
-                            8.0 * polyx[i][0] * polyy[j][1] * polyz[k][2];
-
-
-                          data.shape_hessians[dof_index1][q][0][0][2] =
-                            data.shape_hessians[dof_index1][q][2][0][0];
-                          data.shape_hessians[dof_index1][q][1][0][2] =
-                            data.shape_hessians[dof_index1][q][2][1][0];
-                          data.shape_hessians[dof_index1][q][2][0][2] =
-                            data.shape_hessians[dof_index1][q][2][2][0];
-
-                          data.shape_hessians[dof_index1][q][0][1][2] =
-                            data.shape_hessians[dof_index1][q][2][1][0];
-                          data.shape_hessians[dof_index1][q][1][1][2] =
-                            data.shape_hessians[dof_index1][q][2][1][1];
-                          data.shape_hessians[dof_index1][q][2][1][2] =
-                            data.shape_hessians[dof_index1][q][2][2][1];
-
-                          data.shape_hessians[dof_index1][q][0][2][2] =
-                            data.shape_hessians[dof_index1][q][2][2][0];
-                          data.shape_hessians[dof_index1][q][1][2][2] =
-                            data.shape_hessians[dof_index1][q][2][2][1];
-                          data.shape_hessians[dof_index1][q][2][2][2] =
-                            8.0 * polyx[i][0] * polyy[j][0] * polyz[k][3];
-
-
-                          // Type 2:
-                          const unsigned int dof_index2_1(cell_type2_offset1 +
-                                                          shift_ijk);
-                          const unsigned int dof_index2_2(cell_type2_offset2 +
-                                                          shift_ijk);
+                          const unsigned int shift_k(k * degree * degree);
+                          const unsigned int shift_j(
+                            k * degree); // Used below when subbing
+                                         // k for j type 3
 
-                          for (unsigned int d1 = 0; d1 < dim; ++d1)
+                          for (unsigned int j = 0; j < degree; ++j)
                             {
-                              for (unsigned int d2 = 0; d2 < dim; ++d2)
+                              const unsigned int shift_jk(j * degree + shift_k);
+                              for (unsigned int i = 0; i < degree; ++i)
                                 {
-                                  data.shape_hessians[dof_index2_1][q][0][d1]
-                                                     [d2] =
-                                    data
-                                      .shape_hessians[dof_index1][q][0][d1][d2];
-                                  data.shape_hessians[dof_index2_1][q][1][d1]
-                                                     [d2] =
-                                    -1.0 *
-                                    data
-                                      .shape_hessians[dof_index1][q][1][d1][d2];
-                                  data.shape_hessians[dof_index2_1][q][2][d1]
-                                                     [d2] =
-                                    data
-                                      .shape_hessians[dof_index1][q][2][d1][d2];
-
-                                  data.shape_hessians[dof_index2_2][q][0][d1]
-                                                     [d2] =
-                                    data
-                                      .shape_hessians[dof_index1][q][0][d1][d2];
-                                  data.shape_hessians[dof_index2_2][q][1][d1]
-                                                     [d2] =
-                                    -1.0 *
-                                    data
-                                      .shape_hessians[dof_index1][q][1][d1][d2];
-                                  data.shape_hessians[dof_index2_2][q][2][d1]
-                                                     [d2] =
-                                    -1.0 *
-                                    data
-                                      .shape_hessians[dof_index1][q][2][d1][d2];
+                                  const unsigned int shift_ijk(shift_jk + i);
+
+                                  // Type 1:
+                                  const unsigned int dof_index1(
+                                    cell_type1_offset + shift_ijk);
+
+                                  data.shape_hessians[dof_index1][q][0][0][0] =
+                                    8.0 * polyx[i][3] * polyy[j][0] *
+                                    polyz[k][0];
+                                  data.shape_hessians[dof_index1][q][1][0][0] =
+                                    8.0 * polyx[i][2] * polyy[j][1] *
+                                    polyz[k][0];
+                                  data.shape_hessians[dof_index1][q][2][0][0] =
+                                    8.0 * polyx[i][2] * polyy[j][0] *
+                                    polyz[k][1];
+
+                                  data.shape_hessians[dof_index1][q][0][1][0] =
+                                    data.shape_hessians[dof_index1][q][1][0][0];
+                                  data.shape_hessians[dof_index1][q][1][1][0] =
+                                    8.0 * polyx[i][1] * polyy[j][2] *
+                                    polyz[k][0];
+                                  data.shape_hessians[dof_index1][q][2][1][0] =
+                                    8.0 * polyx[i][1] * polyy[j][1] *
+                                    polyz[k][1];
+
+                                  data.shape_hessians[dof_index1][q][0][2][0] =
+                                    data.shape_hessians[dof_index1][q][2][0][0];
+                                  data.shape_hessians[dof_index1][q][1][2][0] =
+                                    data.shape_hessians[dof_index1][q][2][1][0];
+                                  data.shape_hessians[dof_index1][q][2][2][0] =
+                                    8.0 * polyx[i][1] * polyy[j][0] *
+                                    polyz[k][2];
+
+
+                                  data.shape_hessians[dof_index1][q][0][0][1] =
+                                    data.shape_hessians[dof_index1][q][1][0][0];
+                                  data.shape_hessians[dof_index1][q][1][0][1] =
+                                    data.shape_hessians[dof_index1][q][1][1][0];
+                                  data.shape_hessians[dof_index1][q][2][0][1] =
+                                    data.shape_hessians[dof_index1][q][2][1][0];
+
+                                  data.shape_hessians[dof_index1][q][0][1][1] =
+                                    data.shape_hessians[dof_index1][q][1][1][0];
+                                  data.shape_hessians[dof_index1][q][1][1][1] =
+                                    8.0 * polyx[i][0] * polyy[j][3] *
+                                    polyz[k][0];
+                                  data.shape_hessians[dof_index1][q][2][1][1] =
+                                    8.0 * polyx[i][0] * polyy[j][2] *
+                                    polyz[k][1];
+
+                                  data.shape_hessians[dof_index1][q][0][2][1] =
+                                    data.shape_hessians[dof_index1][q][2][1][0];
+                                  data.shape_hessians[dof_index1][q][1][2][1] =
+                                    data.shape_hessians[dof_index1][q][2][1][1];
+                                  data.shape_hessians[dof_index1][q][2][2][1] =
+                                    8.0 * polyx[i][0] * polyy[j][1] *
+                                    polyz[k][2];
+
+
+                                  data.shape_hessians[dof_index1][q][0][0][2] =
+                                    data.shape_hessians[dof_index1][q][2][0][0];
+                                  data.shape_hessians[dof_index1][q][1][0][2] =
+                                    data.shape_hessians[dof_index1][q][2][1][0];
+                                  data.shape_hessians[dof_index1][q][2][0][2] =
+                                    data.shape_hessians[dof_index1][q][2][2][0];
+
+                                  data.shape_hessians[dof_index1][q][0][1][2] =
+                                    data.shape_hessians[dof_index1][q][2][1][0];
+                                  data.shape_hessians[dof_index1][q][1][1][2] =
+                                    data.shape_hessians[dof_index1][q][2][1][1];
+                                  data.shape_hessians[dof_index1][q][2][1][2] =
+                                    data.shape_hessians[dof_index1][q][2][2][1];
+
+                                  data.shape_hessians[dof_index1][q][0][2][2] =
+                                    data.shape_hessians[dof_index1][q][2][2][0];
+                                  data.shape_hessians[dof_index1][q][1][2][2] =
+                                    data.shape_hessians[dof_index1][q][2][2][1];
+                                  data.shape_hessians[dof_index1][q][2][2][2] =
+                                    8.0 * polyx[i][0] * polyy[j][0] *
+                                    polyz[k][3];
+
+
+                                  // Type 2:
+                                  const unsigned int dof_index2_1(
+                                    cell_type2_offset1 + shift_ijk);
+                                  const unsigned int dof_index2_2(
+                                    cell_type2_offset2 + shift_ijk);
+
+                                  for (unsigned int d1 = 0; d1 < dim; ++d1)
+                                    {
+                                      for (unsigned int d2 = 0; d2 < dim; ++d2)
+                                        {
+                                          data.shape_hessians[dof_index2_1][q]
+                                                             [0][d1][d2] =
+                                            data.shape_hessians[dof_index1][q]
+                                                               [0][d1][d2];
+                                          data.shape_hessians[dof_index2_1][q]
+                                                             [1][d1][d2] =
+                                            -1.0 *
+                                            data.shape_hessians[dof_index1][q]
+                                                               [1][d1][d2];
+                                          data.shape_hessians[dof_index2_1][q]
+                                                             [2][d1][d2] =
+                                            data.shape_hessians[dof_index1][q]
+                                                               [2][d1][d2];
+
+                                          data.shape_hessians[dof_index2_2][q]
+                                                             [0][d1][d2] =
+                                            data.shape_hessians[dof_index1][q]
+                                                               [0][d1][d2];
+                                          data.shape_hessians[dof_index2_2][q]
+                                                             [1][d1][d2] =
+                                            -1.0 *
+                                            data.shape_hessians[dof_index1][q]
+                                                               [1][d1][d2];
+                                          data.shape_hessians[dof_index2_2][q]
+                                                             [2][d1][d2] =
+                                            -1.0 *
+                                            data.shape_hessians[dof_index1][q]
+                                                               [2][d1][d2];
+                                        }
+                                    }
                                 }
-                            }
-                        }
-                      // Type 3: (note we re-use k and j for
-                      // convenience):
-                      const unsigned int shift_ij(j +
-                                                  shift_j); // here we've subbed
-                                                            // j for i, k for j.
-                      const unsigned int dof_index3_1(cell_type3_offset1 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_2(cell_type3_offset2 +
-                                                      shift_ij);
-                      const unsigned int dof_index3_3(cell_type3_offset3 +
-                                                      shift_ij);
-                      for (unsigned int d1 = 0; d1 < dim; ++d1)
-                        {
-                          for (unsigned int d2 = 0; d2 < dim; ++d2)
-                            {
-                              for (unsigned int d3 = 0; d3 < dim; ++d3)
+                              // Type 3: (note we re-use k and j for
+                              // convenience):
+                              const unsigned int shift_ij(
+                                j + shift_j); // here we've subbed
+                                              // j for i, k for j.
+                              const unsigned int dof_index3_1(
+                                cell_type3_offset1 + shift_ij);
+                              const unsigned int dof_index3_2(
+                                cell_type3_offset2 + shift_ij);
+                              const unsigned int dof_index3_3(
+                                cell_type3_offset3 + shift_ij);
+                              for (unsigned int d1 = 0; d1 < dim; ++d1)
                                 {
-                                  data.shape_hessians[dof_index3_1][q][d1][d2]
-                                                     [d3] = 0.0;
-                                  data.shape_hessians[dof_index3_2][q][d1][d2]
-                                                     [d3] = 0.0;
-                                  data.shape_hessians[dof_index3_3][q][d1][d2]
-                                                     [d3] = 0.0;
+                                  for (unsigned int d2 = 0; d2 < dim; ++d2)
+                                    {
+                                      for (unsigned int d3 = 0; d3 < dim; ++d3)
+                                        {
+                                          data.shape_hessians[dof_index3_1][q]
+                                                             [d1][d2][d3] = 0.0;
+                                          data.shape_hessians[dof_index3_2][q]
+                                                             [d1][d2][d3] = 0.0;
+                                          data.shape_hessians[dof_index3_3][q]
+                                                             [d1][d2][d3] = 0.0;
+                                        }
+                                    }
                                 }
-                            }
-                        }
-                      data.shape_hessians[dof_index3_1][q][0][1][1] =
-                        4.0 * polyy[j][2] * polyz[k][0];
-                      data.shape_hessians[dof_index3_1][q][0][1][2] =
-                        4.0 * polyy[j][1] * polyz[k][1];
+                              data.shape_hessians[dof_index3_1][q][0][1][1] =
+                                4.0 * polyy[j][2] * polyz[k][0];
+                              data.shape_hessians[dof_index3_1][q][0][1][2] =
+                                4.0 * polyy[j][1] * polyz[k][1];
 
-                      data.shape_hessians[dof_index3_1][q][0][2][1] =
-                        data.shape_hessians[dof_index3_1][q][0][1][2];
-                      data.shape_hessians[dof_index3_1][q][0][2][2] =
-                        4.0 * polyy[j][0] * polyz[k][2];
+                              data.shape_hessians[dof_index3_1][q][0][2][1] =
+                                data.shape_hessians[dof_index3_1][q][0][1][2];
+                              data.shape_hessians[dof_index3_1][q][0][2][2] =
+                                4.0 * polyy[j][0] * polyz[k][2];
 
 
-                      data.shape_hessians[dof_index3_2][q][1][0][0] =
-                        4.0 * polyx[j][2] * polyz[k][0];
-                      data.shape_hessians[dof_index3_2][q][1][0][2] =
-                        4.0 * polyx[j][1] * polyz[k][1];
+                              data.shape_hessians[dof_index3_2][q][1][0][0] =
+                                4.0 * polyx[j][2] * polyz[k][0];
+                              data.shape_hessians[dof_index3_2][q][1][0][2] =
+                                4.0 * polyx[j][1] * polyz[k][1];
 
-                      data.shape_hessians[dof_index3_2][q][1][2][0] =
-                        data.shape_hessians[dof_index3_2][q][1][0][2];
-                      data.shape_hessians[dof_index3_2][q][1][2][2] =
-                        4.0 * polyx[j][0] * polyz[k][2];
+                              data.shape_hessians[dof_index3_2][q][1][2][0] =
+                                data.shape_hessians[dof_index3_2][q][1][0][2];
+                              data.shape_hessians[dof_index3_2][q][1][2][2] =
+                                4.0 * polyx[j][0] * polyz[k][2];
 
 
-                      data.shape_hessians[dof_index3_3][q][2][0][0] =
-                        4.0 * polyx[j][2] * polyy[k][0];
-                      data.shape_hessians[dof_index3_3][q][2][0][1] =
-                        4.0 * polyx[j][1] * polyy[k][1];
+                              data.shape_hessians[dof_index3_3][q][2][0][0] =
+                                4.0 * polyx[j][2] * polyy[k][0];
+                              data.shape_hessians[dof_index3_3][q][2][0][1] =
+                                4.0 * polyx[j][1] * polyy[k][1];
 
-                      data.shape_hessians[dof_index3_3][q][2][1][0] =
-                        data.shape_hessians[dof_index3_3][q][2][0][1];
-                      data.shape_hessians[dof_index3_3][q][2][1][1] =
-                        4.0 * polyx[j][0] * polyy[k][2];
+                              data.shape_hessians[dof_index3_3][q][2][1][0] =
+                                data.shape_hessians[dof_index3_3][q][2][0][1];
+                              data.shape_hessians[dof_index3_3][q][2][1][1] =
+                                4.0 * polyx[j][0] * polyy[k][2];
+                            }
+                        }
                     }
                 }
             }
+          break;
+        }
+
+      default:
+        {
+          DEAL_II_NOT_IMPLEMENTED();
+          break;
         }
     }
 }
@@ -1448,8 +1460,7 @@ FE_NedelecSZ<dim, spacedim>::get_data(
     typename dealii::FiniteElement<dim, spacedim>::InternalDataBase>
     data_ptr = std::make_unique<InternalData>();
 
-  const unsigned int n_q_points = quadrature.size();
-
+  const unsigned int      n_q_points = quadrature.size();
   std::vector<Point<dim>> p_list(n_q_points);
   p_list = quadrature.get_points();
 
@@ -3310,6 +3321,8 @@ FE_NedelecSZ<dim, spacedim>::create_polynomials(const unsigned int degree)
     IntegratedLegendreSZ::generate_complete_basis(degree + 1);
 }
 
+
+
 // explicit instantiations
 #include "fe_nedelec_sz.inst"
 

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.