]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Vectorize evaluation of hanging node constraints 12994/head
authorMartin Kronbichler <martin.kronbichler@it.uu.se>
Thu, 25 Nov 2021 19:12:56 +0000 (20:12 +0100)
committerMartin Kronbichler <martin.kronbichler@it.uu.se>
Fri, 26 Nov 2021 12:24:25 +0000 (13:24 +0100)
include/deal.II/matrix_free/evaluation_kernels_hanging_nodes.h

index 4e00304ccc62277f58fe96ce585179531f0234bd..0d8d9d35b7394e55eb0d1f4857323a827538419a 100644 (file)
@@ -65,60 +65,22 @@ namespace internal
     }
 
   private:
-    template <int fe_degree, unsigned int side, bool transpose>
-    static void
-    interpolate_2D(const unsigned int given_degree,
-                   const unsigned int v,
-                   const Number *     weight,
-                   Number *           values)
-    {
-      typename Number::value_type temp[40];
-
-      const unsigned int points =
-        (fe_degree != -1 ? fe_degree : given_degree) + 1;
-
-      AssertIndexRange(points, 40);
-
-      const unsigned int d = side / 2; // direction
-      const unsigned int s = side % 2; // left or right surface
-
-      const unsigned int offset = dealii::Utilities::pow(points, d + 1);
-      const unsigned int stride =
-        (s == 0 ? 0 : (points - 1)) * dealii::Utilities::pow(points, d);
-
-      const unsigned int r1 = dealii::Utilities::pow(points, dim - d - 1);
-      const unsigned int r2 = dealii::Utilities::pow(points, d);
-
-      // copy result back
-      for (unsigned int i = 0, k = 0; i < r1; ++i)
-        for (unsigned int j = 0; j < r2; ++j, ++k)
-          temp[k] = values[i * offset + stride + j][v];
-
-      // perform interpolation point by point (note: r1 * r2 == points^(dim-1))
-      for (unsigned int i = 0, k = 0; i < r1; ++i)
-        for (unsigned int j = 0; j < r2; ++j, ++k)
-          {
-            typename Number::value_type sum = 0.0;
-            for (unsigned int h = 0; h < points; ++h)
-              sum += weight[(transpose ? 1 : points) * k +
-                            (transpose ? points : 1) * h][v] *
-                     temp[h];
-            values[i * offset + stride + j][v] = sum;
-          }
-    }
-
-    template <int          fe_degree,
+    template <int          structdim,
+              int          fe_degree,
               unsigned int direction,
-              unsigned int d,
               bool         transpose>
     static void
-    interpolate_3D_face(const unsigned int dof_offset,
-                        const unsigned int given_degree,
-                        const unsigned int v,
-                        const Number *     weight,
-                        Number *           values)
+    interpolate(const unsigned int             offset,
+                const unsigned int             outer_stride,
+                const unsigned int             given_degree,
+                const Number                   mask_weight,
+                const Number                   mask_write,
+                const Number *DEAL_II_RESTRICT weights,
+                Number *DEAL_II_RESTRICT       values)
     {
-      typename Number::value_type temp[40];
+      static_assert(structdim == 1 || structdim == 2,
+                    "Only 1D and 2D interpolation implemented");
+      Number temp[fe_degree != -1 ? fe_degree + 1 : 40];
 
       const unsigned int points =
         (fe_degree != -1 ? fe_degree : given_degree) + 1;
@@ -127,71 +89,40 @@ namespace internal
 
       const unsigned int stride = Utilities::pow(points, direction);
 
-      // direction   side0   side1   side2
-      // 0             -      p^2      p
-      // 1            p^2      -       1
-      // 2             p       -       1
-      const unsigned int stride2 =
-        ((direction == 0 && d == 1) || (direction == 1 && d == 0)) ?
-          (points * points) :
-          (((direction == 0 && d == 2) || (direction == 2 && d == 0)) ? points :
-                                                                        1);
-
-      for (unsigned int g = 1; g < points - 1; ++g)
+      const unsigned int end_of_outer_loop = structdim == 1 ? 2 : points - 1;
+      for (unsigned int g = 1; g < end_of_outer_loop; ++g)
         {
-          // copy result back
+          const unsigned int my_offset =
+            offset + (structdim > 1 ? g * outer_stride : 0);
+
+          // extract values to interpolate
           for (unsigned int k = 0; k < points; ++k)
-            temp[k] = values[dof_offset + k * stride + stride2 * g][v];
+            temp[k] = values[my_offset + k * stride];
 
-          // perform interpolation point by point
+          // perform interpolation point by point and write back
           for (unsigned int k = 0; k < points; ++k)
             {
-              typename Number::value_type sum = 0.0;
+              Number sum0 = Number(), sum1 = Number();
               for (unsigned int h = 0; h < points; ++h)
-                sum += weight[(transpose ? 1 : points) * k +
-                              (transpose ? points : 1) * h][v] *
-                       temp[h];
-              values[dof_offset + k * stride + stride2 * g][v] = sum;
+                {
+                  // sum0 belongs to the other interpolation matrix, but we
+                  // can use the symmetry here to reduce the number of loads
+                  sum0 += temp[h] *
+                          weights[(transpose ? 1 : points) * (points - 1 - k) +
+                                  (transpose ? points : 1) * (points - 1 - h)];
+                  sum1 += temp[h] * weights[(transpose ? 1 : points) * k +
+                                            (transpose ? points : 1) * h];
+                }
+              values[my_offset + k * stride] =
+                temp[k] +
+                mask_write * (sum0 + mask_weight * (sum1 - sum0) - temp[k]);
             }
         }
     }
 
-    template <int fe_degree, unsigned int direction, bool transpose>
-    static void
-    interpolate_3D_edge(const unsigned int p,
-                        const unsigned int given_degree,
-                        const unsigned int v,
-                        const Number *     weight,
-                        Number *           values)
-    {
-      typename Number::value_type temp[40];
-
-      const unsigned int points =
-        (fe_degree != -1 ? fe_degree : given_degree) + 1;
-
-      AssertIndexRange(points, 40);
-
-      const unsigned int stride = Utilities::pow(points, direction);
-
-      // copy result back
-      for (unsigned int k = 0; k < points; ++k)
-        temp[k] = values[p + k * stride][v];
-
-      // perform interpolation point by point
-      for (unsigned int k = 0; k < points; ++k)
-        {
-          typename Number::value_type sum = 0.0;
-          for (unsigned int h = 0; h < points; ++h)
-            sum += weight[(transpose ? 1 : points) * k +
-                          (transpose ? points : 1) * h][v] *
-                   temp[h];
-          values[p + k * stride][v] = sum;
-        }
-    }
-
     template <int fe_degree, bool transpose>
     static void
-    run_internal(const unsigned int                  n_desired_components,
+    run_internal(const unsigned int                  n_components,
                  const FEEvaluationBaseData<dim,
                                             typename Number::value_type,
                                             is_face,
@@ -205,8 +136,11 @@ namespace internal
         fe_degree != -1 ? fe_degree :
                           fe_eval.get_shape_info().data.front().fe_degree;
 
-      const auto &interpolation_matrices =
-        fe_eval.get_shape_info().data.front().subface_interpolation_matrices;
+      const Number *DEAL_II_RESTRICT weights =
+        fe_eval.get_shape_info()
+          .data.front()
+          .subface_interpolation_matrices[0]
+          .data();
 
       const auto is_set = [](const auto a, const auto b) -> bool {
         return (a & b) == b;
@@ -217,210 +151,252 @@ namespace internal
       };
 
       const unsigned int points = given_degree + 1;
+      const unsigned int n_dofs =
+        fe_eval.get_shape_info().dofs_per_component_on_cell;
 
-      for (unsigned int c = 0; c < n_desired_components; ++c)
+      if (dim == 2)
         {
+          // x direction
+          {
+            Number                mask_weights = {};
+            std::array<Number, 2> mask_write   = {};
+            std::array<bool, 2>   do_face      = {};
+            for (unsigned int v = 0; v < Number::size(); ++v)
+              {
+                const Kinds mask = constraint_mask[v];
+                if (is_set(mask, Kinds::face_y))
+                  {
+                    const unsigned int side = not_set(mask, Kinds::subcell_y);
+                    mask_write[side][v]     = 1;
+                    do_face[side]           = true;
+                    mask_weights[v]         = is_set(mask, Kinds::subcell_x);
+                  }
+              }
+            unsigned int offsets[2] = {0, (points - 1) * points};
+            for (unsigned int c = 0; c < n_components; ++c)
+              for (unsigned int face = 0; face < 2; ++face)
+                if (do_face[face])
+                  interpolate<1, fe_degree, 0, transpose>(offsets[face],
+                                                          0,
+                                                          given_degree,
+                                                          mask_weights,
+                                                          mask_write[face],
+                                                          weights,
+                                                          values + c * n_dofs);
+          }
+          // y direction
+          {
+            Number                mask_weights = {};
+            std::array<Number, 2> mask_write   = {};
+            std::array<bool, 2>   do_face      = {};
+            for (unsigned int v = 0; v < Number::size(); ++v)
+              {
+                const Kinds mask = constraint_mask[v];
+                if (is_set(mask, Kinds::face_x))
+                  {
+                    const unsigned int side = not_set(mask, Kinds::subcell_x);
+                    mask_write[side][v]     = 1;
+                    do_face[side]           = true;
+                    mask_weights[v]         = is_set(mask, Kinds::subcell_y);
+                  }
+              }
+            unsigned int offsets[2] = {0, points - 1};
+            for (unsigned int c = 0; c < n_components; ++c)
+              for (unsigned int face = 0; face < 2; ++face)
+                if (do_face[face])
+                  interpolate<1, fe_degree, 1, transpose>(offsets[face],
+                                                          0,
+                                                          given_degree,
+                                                          mask_weights,
+                                                          mask_write[face],
+                                                          weights,
+                                                          values + c * n_dofs);
+          }
+        }
+      else if (dim == 3)
+        {
+          const unsigned int p0 = 0;
+          const unsigned int p1 = points - 1;
+          const unsigned int p2 = points * points - points;
+          const unsigned int p3 = points * points - 1;
+          const unsigned int p4 = points * points * points - points * points;
+          const unsigned int p5 =
+            points * points * points - points * points + points - 1;
+          const unsigned int p6 = points * points * points - points;
+
+          std::array<std::array<bool, 4>, 3>   process_edge = {};
+          std::array<std::array<bool, 4>, 3>   process_face = {};
+          std::array<std::array<Number, 4>, 3> mask_edge    = {};
+          std::array<std::array<Number, 4>, 3> mask_face    = {};
+          std::array<Number, 3>                mask_weights = {};
           for (unsigned int v = 0; v < Number::size(); ++v)
             {
-              const auto mask = constraint_mask[v];
+              const Kinds mask = constraint_mask[v];
 
-              if (mask == Kinds::unconstrained)
-                continue;
+              if (is_set(mask, Kinds::subcell_x))
+                mask_weights[0][v] = 1;
+              if (is_set(mask, Kinds::subcell_y))
+                mask_weights[1][v] = 1;
+              if (is_set(mask, Kinds::subcell_z))
+                mask_weights[2][v] = 1;
 
-              if (dim == 2) // 2D: only faces
+              if (is_set(mask, Kinds::face_x))
                 {
-                  // direction 0:
-                  if ((mask & Kinds::face_y) != Kinds::unconstrained)
-                    {
-                      const bool is_subface_0 =
-                        (mask & Kinds::subcell_x) == Kinds::unconstrained;
-
-                      const Number *weights =
-                        interpolation_matrices[is_subface_0].data();
-
-                      if (is_set(mask, Kinds::subcell_y))
-                        interpolate_2D<fe_degree, 2, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 2
-                      else
-                        interpolate_2D<fe_degree, 3, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 3
-                    }
-
-                  // direction 1:
-                  if ((mask & Kinds::face_x) != Kinds::unconstrained)
-                    {
-                      const bool is_subface_0 =
-                        (mask & Kinds::subcell_y) == Kinds::unconstrained;
-
-                      const Number *weights =
-                        interpolation_matrices[is_subface_0].data();
-
-                      if (is_set(mask, Kinds::subcell_x))
-                        interpolate_2D<fe_degree, 0, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 0
-                      else
-                        interpolate_2D<fe_degree, 1, transpose>(
-                          given_degree,
-                          v,
-                          weights,
-                          values); // face 1
-                    }
+                  const unsigned int side = not_set(mask, Kinds::subcell_x);
+
+                  mask_face[1][side][v] = process_face[1][side] = true;
+                  mask_edge[1][side][v] = process_edge[1][side] = true;
+                  mask_edge[1][2 + side][v] = process_edge[1][2 + side] = true;
+                  mask_face[2][side][v] = process_face[2][side] = true;
+                  mask_edge[2][side][v] = process_edge[2][side] = true;
+                  mask_edge[2][2 + side][v] = process_edge[2][2 + side] = true;
                 }
-              else if (dim == 3) // 3D faces and edges
+              if (is_set(mask, Kinds::face_y))
                 {
-                  const unsigned int p0 = 0;
-                  const unsigned int p1 = points - 1;
-                  const unsigned int p2 = points * points - points;
-                  const unsigned int p3 = points * points - 1;
-                  const unsigned int p4 =
-                    points * points * points - points * points;
-                  const unsigned int p5 =
-                    points * points * points - points * points + points - 1;
-                  const unsigned int p6 = points * points * points - points;
-
-                  std::array<std::array<char, 4>, 3> process_edge = {};
-                  std::array<std::array<char, 4>, 3> process_face = {};
-
-                  if (is_set(mask, Kinds::face_x))
-                    {
-                      const unsigned int side = not_set(mask, Kinds::subcell_x);
-                      process_face[1][side]   = 1;
-                      process_edge[1][side]   = 1;
-                      process_edge[1][2 + side] = 1;
-                      process_face[2][side]     = 1;
-                      process_edge[2][side]     = 1;
-                      process_edge[2][2 + side] = 1;
-                    }
-                  if (is_set(mask, Kinds::face_y))
-                    {
-                      const unsigned int side = not_set(mask, Kinds::subcell_y);
-                      process_face[0][side]   = 1;
-                      process_edge[0][side]   = 1;
-                      process_edge[0][2 + side]     = 1;
-                      process_face[2][2 + side]     = 1;
-                      process_edge[2][2 * side]     = 1;
-                      process_edge[2][2 * side + 1] = 1;
-                    }
-                  if (is_set(mask, Kinds::face_z))
-                    {
-                      const unsigned int side = not_set(mask, Kinds::subcell_z);
-                      process_face[0][2 + side]     = 1;
-                      process_edge[0][2 * side]     = 1;
-                      process_edge[0][2 * side + 1] = 1;
-                      process_face[1][2 + side]     = 1;
-                      process_edge[1][2 * side]     = 1;
-                      process_edge[1][2 * side + 1] = 1;
-                    }
-                  if (is_set(mask, Kinds::edge_x))
-                    process_edge[0][not_set(mask, Kinds::subcell_z) * 2 +
-                                    not_set(mask, Kinds::subcell_y)] = 1;
-                  if (is_set(mask, Kinds::edge_y))
-                    process_edge[1][not_set(mask, Kinds::subcell_z) * 2 +
-                                    not_set(mask, Kinds::subcell_x)] = 1;
-                  if (is_set(mask, Kinds::edge_z))
-                    process_edge[2][not_set(mask, Kinds::subcell_y) * 2 +
-                                    not_set(mask, Kinds::subcell_x)] = 1;
-
-                  // direction 0:
-                  {
-                    const bool is_subface_0 =
-                      (mask & Kinds::subcell_x) == Kinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    unsigned int face_offsets[4] = {p0, p2, p0, p4};
-                    // face 2, 3
-                    for (unsigned int face = 0; face < 2; ++face)
-                      if (process_face[0][face])
-                        interpolate_3D_face<fe_degree, 0, 1, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-                    // face 4, 5
-                    for (unsigned int face = 2; face < 4; ++face)
-                      if (process_face[0][face])
-                        interpolate_3D_face<fe_degree, 0, 2, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-
-                    // edges
-                    unsigned int edge_offsets[4] = {p0, p2, p4, p6};
-                    for (unsigned int edge = 0; edge < 4; ++edge)
-                      if (process_edge[0][edge])
-                        interpolate_3D_edge<fe_degree, 0, transpose>(
-                          edge_offsets[edge], given_degree, v, weights, values);
-                  }
-
-                  // direction 1:
-                  {
-                    const bool is_subface_0 =
-                      (mask & Kinds::subcell_y) == Kinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    unsigned int face_offsets[4] = {p0, p1, p0, p4};
-                    // face 0, 1
-                    for (unsigned int face = 0; face < 2; ++face)
-                      if (process_face[1][face])
-                        interpolate_3D_face<fe_degree, 1, 0, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-                    // face 4, 5
-                    for (unsigned int face = 2; face < 4; ++face)
-                      if (process_face[1][face])
-                        interpolate_3D_face<fe_degree, 1, 2, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-
-                    // edges
-                    unsigned int edge_offsets[4] = {p0, p1, p4, p5};
-                    for (unsigned int edge = 0; edge < 4; ++edge)
-                      if (process_edge[1][edge])
-                        interpolate_3D_edge<fe_degree, 1, transpose>(
-                          edge_offsets[edge], given_degree, v, weights, values);
-                  }
-
-                  // direction 2:
-                  {
-                    const bool is_subface_0 =
-                      (mask & Kinds::subcell_z) == Kinds::unconstrained;
-
-                    const Number *weights =
-                      interpolation_matrices[is_subface_0].data();
-
-                    unsigned int face_offsets[4] = {p0, p1, p0, p2};
-                    // face 0, 1
-                    for (unsigned int face = 0; face < 2; ++face)
-                      if (process_face[2][face])
-                        interpolate_3D_face<fe_degree, 2, 0, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-                    // face 2, 3
-                    for (unsigned int face = 2; face < 4; ++face)
-                      if (process_face[2][face])
-                        interpolate_3D_face<fe_degree, 2, 1, transpose>(
-                          face_offsets[face], given_degree, v, weights, values);
-
-                    // edges
-                    unsigned int edge_offsets[4] = {p0, p1, p2, p3};
-                    for (unsigned int edge = 0; edge < 4; ++edge)
-                      if (process_edge[2][edge])
-                        interpolate_3D_edge<fe_degree, 2, transpose>(
-                          edge_offsets[edge], given_degree, v, weights, values);
-                  }
+                  const unsigned int side = not_set(mask, Kinds::subcell_y);
+
+                  mask_face[0][side][v] = process_face[0][side] = true;
+                  mask_edge[0][side][v] = process_edge[0][side] = true;
+                  mask_edge[0][2 + side][v] = process_edge[0][2 + side] = true;
+                  mask_face[2][2 + side][v] = process_face[2][2 + side] = true;
+                  mask_edge[2][2 * side][v] = process_edge[2][2 * side] = true;
+                  mask_edge[2][2 * side + 1][v] =
+                    process_edge[2][2 * side + 1] = true;
+                }
+              if (is_set(mask, Kinds::face_z))
+                {
+                  const unsigned int side = not_set(mask, Kinds::subcell_z);
+
+                  mask_face[0][2 + side][v] = process_face[0][2 + side] = true;
+                  mask_edge[0][2 * side][v] = process_edge[0][2 * side] = true;
+                  mask_edge[0][2 * side + 1][v] =
+                    process_edge[0][2 * side + 1] = true;
+                  mask_face[1][2 + side][v] = process_face[1][2 + side] = true;
+                  mask_edge[1][2 * side][v] = process_edge[1][2 * side] = true;
+                  mask_edge[1][2 * side + 1][v] =
+                    process_edge[1][2 * side + 1] = true;
+                }
+              if (is_set(mask, Kinds::edge_x))
+                {
+                  const unsigned int index =
+                    not_set(mask, Kinds::subcell_z) * 2 +
+                    not_set(mask, Kinds::subcell_y);
+                  mask_edge[0][index][v] = process_edge[0][index] = true;
+                }
+              if (is_set(mask, Kinds::edge_y))
+                {
+                  const unsigned int index =
+                    not_set(mask, Kinds::subcell_z) * 2 +
+                    not_set(mask, Kinds::subcell_x);
+                  mask_edge[1][index][v] = process_edge[1][index] = true;
                 }
-              else
+              if (is_set(mask, Kinds::edge_z))
                 {
-                  Assert(false, ExcNotImplemented());
+                  const unsigned int index =
+                    not_set(mask, Kinds::subcell_y) * 2 +
+                    not_set(mask, Kinds::subcell_x);
+                  mask_edge[2][index][v] = process_edge[2][index] = true;
                 }
             }
 
-          values += fe_eval.get_shape_info().dofs_per_component_on_cell;
+          // direction 0:
+          if (given_degree > 1)
+            {
+              unsigned int face_offsets[4]  = {p0, p2, p0, p4};
+              unsigned int outer_strides[2] = {points * points, points};
+              for (unsigned int c = 0; c < n_components; ++c)
+                for (unsigned int face = 0; face < 4; ++face)
+                  if (process_face[0][face])
+                    interpolate<2, fe_degree, 0, transpose>(
+                      face_offsets[face],
+                      outer_strides[face / 2],
+                      given_degree,
+                      mask_weights[0],
+                      mask_face[0][face],
+                      weights,
+                      values + c * n_dofs);
+            }
+          {
+            unsigned int edge_offsets[4] = {p0, p2, p4, p6};
+            for (unsigned int c = 0; c < n_components; ++c)
+              for (unsigned int edge = 0; edge < 4; ++edge)
+                if (process_edge[0][edge])
+                  interpolate<1, fe_degree, 0, transpose>(edge_offsets[edge],
+                                                          0,
+                                                          given_degree,
+                                                          mask_weights[0],
+                                                          mask_edge[0][edge],
+                                                          weights,
+                                                          values + c * n_dofs);
+          }
+
+          // direction 1:
+          if (given_degree > 1)
+            {
+              unsigned int face_offsets[4]  = {p0, p1, p0, p4};
+              unsigned int outer_strides[2] = {points * points, 1};
+              for (unsigned int c = 0; c < n_components; ++c)
+                for (unsigned int face = 0; face < 4; ++face)
+                  if (process_face[1][face])
+                    interpolate<2, fe_degree, 1, transpose>(
+                      face_offsets[face],
+                      outer_strides[face / 2],
+                      given_degree,
+                      mask_weights[1],
+                      mask_face[1][face],
+                      weights,
+                      values + c * n_dofs);
+            }
+
+          {
+            unsigned int edge_offsets[4] = {p0, p1, p4, p5};
+            for (unsigned int c = 0; c < n_components; ++c)
+              for (unsigned int edge = 0; edge < 4; ++edge)
+                if (process_edge[1][edge])
+                  interpolate<1, fe_degree, 1, transpose>(edge_offsets[edge],
+                                                          0,
+                                                          given_degree,
+                                                          mask_weights[1],
+                                                          mask_edge[1][edge],
+                                                          weights,
+                                                          values + c * n_dofs);
+          }
+
+          // direction 2:
+          if (given_degree > 1)
+            {
+              unsigned int face_offsets[4]  = {p0, p1, p0, p2};
+              unsigned int outer_strides[2] = {points, 1};
+              for (unsigned int c = 0; c < n_components; ++c)
+                for (unsigned int face = 0; face < 4; ++face)
+                  if (process_face[2][face])
+                    interpolate<2, fe_degree, 2, transpose>(
+                      face_offsets[face],
+                      outer_strides[face / 2],
+                      given_degree,
+                      mask_weights[2],
+                      mask_face[2][face],
+                      weights,
+                      values + c * n_dofs);
+            }
+
+          {
+            unsigned int edge_offsets[4] = {p0, p1, p2, p3};
+            for (unsigned int c = 0; c < n_components; ++c)
+              for (unsigned int edge = 0; edge < 4; ++edge)
+                if (process_edge[2][edge])
+                  interpolate<1, fe_degree, 2, transpose>(edge_offsets[edge],
+                                                          0,
+                                                          given_degree,
+                                                          mask_weights[2],
+                                                          mask_edge[2][edge],
+                                                          weights,
+                                                          values + c * n_dofs);
+          }
+        }
+      else
+        {
+          Assert(false, ExcNotImplemented());
         }
     }
   };

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.