From d97d9c31fd7b775e57f6e5b714f436ee9a33acff Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 22 Feb 2021 19:32:19 -0700
Subject: [PATCH] Use the new ndarray class instead of nested std::array.

---
 examples/step-33/step-33.cc                   |  21 +-
 include/deal.II/base/geometry_info.h          |  22 +-
 include/deal.II/base/polynomial_space.h       |   3 +-
 .../deal.II/base/tensor_product_polynomials.h |   3 +-
 include/deal.II/fe/fe_p1nc.h                  |   2 +-
 include/deal.II/grid/connectivity.h           |  81 +++----
 include/deal.II/grid/reference_cell.h         | 207 +++++++++---------
 include/deal.II/meshworker/copy_data.h        |   9 +-
 source/base/geometry_info.cc                  |   7 +-
 source/base/tensor_product_polynomials.cc     |  12 +-
 source/fe/fe_p1nc.cc                          |  13 +-
 source/grid/grid_generator.cc                 | 126 ++++++-----
 source/simplex/polynomials.cc                 |  45 ++--
 13 files changed, 269 insertions(+), 282 deletions(-)

diff --git a/examples/step-33/step-33.cc b/examples/step-33/step-33.cc
index 4c5fcd3422..e9bf0c15e4 100644
--- a/examples/step-33/step-33.cc
+++ b/examples/step-33/step-33.cc
@@ -207,10 +207,10 @@ namespace Step33
     // the function with different input vector data types, so we templatize
     // on it as well:
     template <typename InputVector>
-    static void compute_flux_matrix(
-      const InputVector &                            W,
-      std::array<std::array<typename InputVector::value_type, dim>,
-                 EulerEquations<dim>::n_components> &flux)
+    static void compute_flux_matrix(const InputVector &W,
+                                    ndarray<typename InputVector::value_type,
+                                            EulerEquations<dim>::n_components,
+                                            dim> &     flux)
     {
       // First compute the pressure that appears in the flux matrix, and then
       // compute the first <code>dim</code> columns of the matrix that
@@ -253,8 +253,9 @@ namespace Step33
       const double                                                alpha,
       std::array<typename InputVector::value_type, n_components> &normal_flux)
     {
-      std::array<std::array<typename InputVector::value_type, dim>,
-                 EulerEquations<dim>::n_components>
+      ndarray<typename InputVector::value_type,
+              EulerEquations<dim>::n_components,
+              dim>
         iflux, oflux;
 
       compute_flux_matrix(Wplus, iflux);
@@ -1776,12 +1777,12 @@ namespace Step33
     // terms of autodifferentiation variables, so that the Jacobian
     // contributions can later easily be computed from it:
 
-    std::vector<std::array<std::array<Sacado::Fad::DFad<double>, dim>,
-                           EulerEquations<dim>::n_components>>
+    std::vector<ndarray<Sacado::Fad::DFad<double>,
+                        EulerEquations<dim>::n_components,
+                        dim>>
       flux(n_q_points);
 
-    std::vector<
-      std::array<std::array<double, dim>, EulerEquations<dim>::n_components>>
+    std::vector<ndarray<double, EulerEquations<dim>::n_components, dim>>
       flux_old(n_q_points);
 
     std::vector<
diff --git a/include/deal.II/base/geometry_info.h b/include/deal.II/base/geometry_info.h
index 4bb2e017b2..6f7c4fd39e 100644
--- a/include/deal.II/base/geometry_info.h
+++ b/include/deal.II/base/geometry_info.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/std_cxx20/iota_view.h>
 
@@ -67,7 +68,7 @@ namespace internal
         return {{Tensor<1, 1>{{-1}}, Tensor<1, 1>{{1}}}};
       }
 
-      static constexpr std::array<std::array<Tensor<1, 1>, 0>, 2>
+      static constexpr dealii::ndarray<Tensor<1, 1>, 2, 0>
       unit_tangential_vectors()
       {
         return {{{{}}, {{}}}};
@@ -85,7 +86,7 @@ namespace internal
         return {{0, 1}};
       }
 
-      static constexpr std::array<std::array<unsigned int, 1>, 2>
+      static constexpr dealii::ndarray<unsigned int, 2, 1>
       vertex_to_face()
       {
         return {{{{0}}, {{1}}}};
@@ -122,7 +123,7 @@ namespace internal
                  Tensor<1, 2>{{0., 1.}}}};
       }
 
-      static constexpr std::array<std::array<Tensor<1, 2>, 1>, 4>
+      static constexpr dealii::ndarray<Tensor<1, 2>, 4, 1>
       unit_tangential_vectors()
       {
         return {{{{Tensor<1, 2>{{0, -1}}}},
@@ -143,7 +144,7 @@ namespace internal
         return {{0, 2, 1, 3}};
       }
 
-      static constexpr std::array<std::array<unsigned int, 2>, 4>
+      static constexpr dealii::ndarray<unsigned int, 4, 2>
       vertex_to_face()
       {
         return {{{{0, 2}}, {{1, 2}}, {{0, 3}}, {{1, 3}}}};
@@ -182,7 +183,7 @@ namespace internal
                  Tensor<1, 3>{{0, 0, 1}}}};
       }
 
-      static constexpr std::array<std::array<Tensor<1, 3>, 2>, 6>
+      static constexpr dealii::ndarray<Tensor<1, 3>, 6, 2>
       unit_tangential_vectors()
       {
         return {{{{Tensor<1, 3>{{0, -1, 0}}, Tensor<1, 3>{{0, 0, 1}}}},
@@ -205,7 +206,7 @@ namespace internal
         return {{0, 4, 2, 6, 1, 5, 3, 7}};
       }
 
-      static constexpr std::array<std::array<unsigned int, 3>, 8>
+      static constexpr dealii::ndarray<unsigned int, 8, 3>
       vertex_to_face()
       {
         return {{{{0, 2, 4}},
@@ -268,7 +269,7 @@ namespace internal
                  Tensor<1, 4>{{0, 0, 0, 1}}}};
       }
 
-      static constexpr std::array<std::array<Tensor<1, 4>, 3>, 8>
+      static constexpr dealii::ndarray<Tensor<1, 4>, 8, 3>
       unit_tangential_vectors()
       {
         return {{{{Tensor<1, 4>{{0, -1, 0, 0}},
@@ -324,7 +325,7 @@ namespace internal
                  numbers::invalid_unsigned_int}};
       }
 
-      static constexpr std::array<std::array<unsigned int, 4>, 16>
+      static constexpr dealii::ndarray<unsigned int, 16, 4>
       vertex_to_face()
       {
         return {{{{numbers::invalid_unsigned_int,
@@ -2116,7 +2117,7 @@ struct GeometryInfo
    * bounds the reference cell in <i>x</i> direction, the second in <i>y</i>
    * direction, and so on.
    */
-  static constexpr std::array<std::array<unsigned int, dim>, vertices_per_cell>
+  static constexpr ndarray<unsigned int, vertices_per_cell, dim>
     vertex_to_face =
       internal::GeometryInfoHelper::Initializers<dim>::vertex_to_face();
 
@@ -2618,8 +2619,7 @@ struct GeometryInfo
    * @ref GlossFaceOrientation "glossary"
    * entry on face orientation.
    */
-  static constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
-                              faces_per_cell>
+  static constexpr ndarray<Tensor<1, dim>, faces_per_cell, dim - 1>
     unit_tangential_vectors = internal::GeometryInfoHelper::Initializers<
       dim>::unit_tangential_vectors();
 
diff --git a/include/deal.II/base/polynomial_space.h b/include/deal.II/base/polynomial_space.h
index 336e65c1fe..c7424d8bb7 100644
--- a/include/deal.II/base/polynomial_space.h
+++ b/include/deal.II/base/polynomial_space.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/scalar_polynomials_base.h>
@@ -334,7 +335,7 @@ PolynomialSpace<dim>::compute_derivative(const unsigned int i,
 {
   const std::array<unsigned int, dim> indices = compute_index(i);
 
-  std::array<std::array<double, order + 1>, dim> v;
+  ndarray<double, dim, order + 1> v;
   {
     std::vector<double> tmp(order + 1);
     for (unsigned int d = 0; d < dim; ++d)
diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h
index 9351d5a8d3..7ed3c46170 100644
--- a/include/deal.II/base/tensor_product_polynomials.h
+++ b/include/deal.II/base/tensor_product_polynomials.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/exceptions.h>
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/point.h>
 #include <deal.II/base/polynomial.h>
 #include <deal.II/base/scalar_polynomials_base.h>
@@ -542,7 +543,7 @@ TensorProductPolynomials<dim, PolynomialType>::compute_derivative(
   std::array<unsigned int, dim> indices;
   compute_index(i, indices);
 
-  std::array<std::array<double, 5>, dim> v;
+  ndarray<double, dim, 5> v;
   {
     std::vector<double> tmp(5);
     for (unsigned int d = 0; d < dim; ++d)
diff --git a/include/deal.II/fe/fe_p1nc.h b/include/deal.II/fe/fe_p1nc.h
index 9144516ac0..968943c230 100644
--- a/include/deal.II/fe/fe_p1nc.h
+++ b/include/deal.II/fe/fe_p1nc.h
@@ -292,7 +292,7 @@ private:
    * x + b y + c$ on given cell. For each local shape function, the array
    * consists of three coefficients is in order of a,b and c.
    */
-  static std::array<std::array<double, 3>, 4>
+  static ndarray<double, 4, 3>
   get_linear_shape_coefficients(const Triangulation<2, 2>::cell_iterator &cell);
 
   /**
diff --git a/include/deal.II/grid/connectivity.h b/include/deal.II/grid/connectivity.h
index 04c951c527..03399eb880 100644
--- a/include/deal.II/grid/connectivity.h
+++ b/include/deal.II/grid/connectivity.h
@@ -19,6 +19,7 @@
 #include <deal.II/base/config.h>
 
 #include <deal.II/base/array_view.h>
+#include <deal.II/base/ndarray.h>
 
 #include <deal.II/grid/reference_cell.h>
 #include <deal.II/grid/tria_description.h>
@@ -193,7 +194,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 3> table = {
+            static const dealii::ndarray<unsigned int, 3, 2> table = {
               {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
             return {table[e]};
@@ -250,7 +251,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 4> table = {
+            static const dealii::ndarray<unsigned int, 4, 2> table = {
               {{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}};
 
             return {table[e]};
@@ -307,7 +308,7 @@ namespace internal
 
         if (d == 2)
           {
-            static const std::array<std::array<unsigned int, 3>, 4> table = {
+            static const dealii::ndarray<unsigned int, 4, 3> table = {
               {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
 
             return {table[e]};
@@ -315,7 +316,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 6> table = {
+            static const dealii::ndarray<unsigned int, 6, 2> table = {
               {{{0, 1}}, {{1, 2}}, {{2, 0}}, {{0, 3}}, {{1, 3}}, {{2, 3}}}};
 
             return {table[e]};
@@ -363,7 +364,7 @@ namespace internal
       nth_line_of_surface(const unsigned int line,
                           const unsigned int face) const override
       {
-        const static std::array<std::array<unsigned int, 3>, 4> table = {
+        const static dealii::ndarray<unsigned int, 4, 3> table = {
           {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
 
         return table[face][line];
@@ -373,11 +374,11 @@ namespace internal
       vertices_of_nth_line_of_surface(const unsigned int line,
                                       const unsigned int face) const override
       {
-        const static std::array<std::array<std::array<unsigned int, 2>, 3>, 4>
-          table = {{{{{{0, 1}}, {{1, 2}}, {{2, 0}}}},
-                    {{{{1, 0}}, {{0, 3}}, {{3, 1}}}},
-                    {{{{0, 2}}, {{2, 3}}, {{3, 0}}}},
-                    {{{{2, 1}}, {{1, 3}}, {{3, 2}}}}}};
+        const static dealii::ndarray<unsigned int, 4, 3, 2> table = {
+          {{{{{0, 1}}, {{1, 2}}, {{2, 0}}}},
+           {{{{1, 0}}, {{0, 3}}, {{3, 1}}}},
+           {{{{0, 2}}, {{2, 3}}, {{3, 0}}}},
+           {{{{2, 1}}, {{1, 3}}, {{3, 2}}}}}};
 
         return table[face][line];
       }
@@ -411,7 +412,7 @@ namespace internal
                 return {table};
               }
 
-            static const std::array<std::array<unsigned int, 3>, 4> table = {
+            static const dealii::ndarray<unsigned int, 4, 3> table = {
               {{{0, 2, 4}}, {{3, 1, 4}}, {{1, 0, 4}}, {{2, 3, 4}}}};
 
             return {table[e - 1]};
@@ -419,7 +420,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 8> table = {
+            static const dealii::ndarray<unsigned int, 8, 2> table = {
               {{{0, 2}},
                {{1, 3}},
                {{0, 1}},
@@ -478,7 +479,7 @@ namespace internal
       nth_line_of_surface(const unsigned int line,
                           const unsigned int face) const override
       {
-        const static std::array<std::array<unsigned int, 4>, 5> table = {
+        const static dealii::ndarray<unsigned int, 5, 4> table = {
           {{{0, 1, 2, 3}},
            {{0, 6, 4, numbers::invalid_unsigned_int}},
            {{1, 5, 7, numbers::invalid_unsigned_int}},
@@ -494,12 +495,12 @@ namespace internal
       {
         static const unsigned int X = static_cast<unsigned int>(-1);
 
-        const static std::array<std::array<std::array<unsigned int, 2>, 4>, 5>
-          table = {{{{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
-                    {{{{0, 2}}, {{2, 4}}, {{4, 0}}, {{X, X}}}},
-                    {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
-                    {{{{1, 0}}, {{0, 4}}, {{4, 1}}, {{X, X}}}},
-                    {{{{2, 3}}, {{3, 4}}, {{4, 2}}, {{X, X}}}}}};
+        const static dealii::ndarray<unsigned int, 5, 4, 2> table = {
+          {{{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
+           {{{{0, 2}}, {{2, 4}}, {{4, 0}}, {{X, X}}}},
+           {{{{3, 1}}, {{1, 4}}, {{4, 3}}, {{X, X}}}},
+           {{{{1, 0}}, {{0, 4}}, {{4, 1}}, {{X, X}}}},
+           {{{{2, 3}}, {{3, 4}}, {{4, 2}}, {{X, X}}}}}};
 
         return table[face][line];
       }
@@ -530,13 +531,13 @@ namespace internal
           {
             if (e == 0 || e == 1)
               {
-                static const std::array<std::array<unsigned int, 3>, 2> table =
-                  {{{{1, 0, 2}}, {{3, 4, 5}}}};
+                static const dealii::ndarray<unsigned int, 2, 3> table = {
+                  {{{1, 0, 2}}, {{3, 4, 5}}}};
 
                 return {table[e]};
               }
 
-            static const std::array<std::array<unsigned int, 4>, 3> table = {
+            static const dealii::ndarray<unsigned int, 3, 4> table = {
               {{{0, 1, 3, 4}}, {{1, 2, 4, 5}}, {{2, 0, 5, 3}}}};
 
             return {table[e - 2]};
@@ -544,7 +545,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 9> table = {
+            static const dealii::ndarray<unsigned int, 9, 2> table = {
               {{{0, 1}},
                {{1, 2}},
                {{2, 0}},
@@ -606,7 +607,7 @@ namespace internal
       {
         static const unsigned int X = static_cast<unsigned int>(-1);
 
-        const static std::array<std::array<unsigned int, 4>, 5> table = {
+        const static dealii::ndarray<unsigned int, 5, 4> table = {
           {{{0, 2, 1, X}},
            {{3, 4, 5, X}},
            {{6, 7, 0, 3}},
@@ -622,12 +623,12 @@ namespace internal
       {
         static const unsigned int X = static_cast<unsigned int>(-1);
 
-        const static std::array<std::array<std::array<unsigned int, 2>, 4>, 5>
-          table = {{{{{{1, 0}}, {{0, 2}}, {{2, 1}}, {{X, X}}}},
-                    {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}},
-                    {{{{0, 3}}, {{1, 4}}, {{0, 1}}, {{3, 4}}}},
-                    {{{{1, 4}}, {{2, 5}}, {{1, 2}}, {{4, 5}}}},
-                    {{{{2, 5}}, {{0, 3}}, {{2, 0}}, {{5, 3}}}}}};
+        const static dealii::ndarray<unsigned int, 5, 4, 2> table = {
+          {{{{{1, 0}}, {{0, 2}}, {{2, 1}}, {{X, X}}}},
+           {{{{3, 4}}, {{4, 5}}, {{5, 3}}, {{X, X}}}},
+           {{{{0, 3}}, {{1, 4}}, {{0, 1}}, {{3, 4}}}},
+           {{{{1, 4}}, {{2, 5}}, {{1, 2}}, {{4, 5}}}},
+           {{{{2, 5}}, {{0, 3}}, {{2, 0}}, {{5, 3}}}}}};
 
         return table[face][line];
       }
@@ -656,7 +657,7 @@ namespace internal
 
         if (d == 2)
           {
-            static const std::array<std::array<unsigned int, 4>, 6> table = {
+            static const dealii::ndarray<unsigned int, 6, 4> table = {
               {{{0, 2, 4, 6}},
                {{1, 3, 5, 7}},
                {{0, 4, 1, 5}},
@@ -669,7 +670,7 @@ namespace internal
 
         if (d == 1)
           {
-            static const std::array<std::array<unsigned int, 2>, 12> table = {
+            static const dealii::ndarray<unsigned int, 12, 2> table = {
               {{{0, 2}},
                {{1, 3}},
                {{0, 1}},
@@ -728,7 +729,7 @@ namespace internal
       nth_line_of_surface(const unsigned int line,
                           const unsigned int face) const override
       {
-        const static std::array<std::array<unsigned int, 4>, 6> table = {
+        const static dealii::ndarray<unsigned int, 6, 4> table = {
           {{{8, 10, 0, 4}},
            {{9, 11, 1, 5}},
            {{2, 6, 8, 9}},
@@ -743,13 +744,13 @@ namespace internal
       vertices_of_nth_line_of_surface(const unsigned int line,
                                       const unsigned int face) const override
       {
-        const static std::array<std::array<std::array<unsigned int, 2>, 4>, 6>
-          table = {{{{{{0, 4}}, {{2, 6}}, {{0, 2}}, {{4, 6}}}},
-                    {{{{1, 5}}, {{3, 7}}, {{1, 3}}, {{5, 7}}}},
-                    {{{{0, 1}}, {{4, 5}}, {{0, 4}}, {{1, 5}}}},
-                    {{{{2, 3}}, {{6, 7}}, {{2, 6}}, {{3, 7}}}},
-                    {{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
-                    {{{{4, 6}}, {{5, 7}}, {{4, 5}}, {{6, 7}}}}}};
+        const static dealii::ndarray<unsigned int, 6, 4, 2> table = {
+          {{{{{0, 4}}, {{2, 6}}, {{0, 2}}, {{4, 6}}}},
+           {{{{1, 5}}, {{3, 7}}, {{1, 3}}, {{5, 7}}}},
+           {{{{0, 1}}, {{4, 5}}, {{0, 4}}, {{1, 5}}}},
+           {{{{2, 3}}, {{6, 7}}, {{2, 6}}, {{3, 7}}}},
+           {{{{0, 2}}, {{1, 3}}, {{0, 1}}, {{2, 3}}}},
+           {{{{4, 6}}, {{5, 7}}, {{4, 5}}, {{6, 7}}}}}};
 
         return table[face][line];
       }
diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 9afa30c311..22469696a8 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/array_view.h>
 #include <deal.II/base/geometry_info.h>
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/utilities.h>
 
@@ -661,7 +662,7 @@ ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
   else if (*this == ReferenceCells::Triangle)
     {
       AssertIndexRange(vertex, 3);
-      static const std::array<std::array<unsigned int, 2>, 3> table = {
+      static const ndarray<unsigned int, 3, 2> table = {
         {{{0, 2}}, {{0, 1}}, {{1, 2}}}};
 
       return table[vertex];
@@ -669,7 +670,7 @@ ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
   else if (*this == ReferenceCells::Tetrahedron)
     {
       AssertIndexRange(vertex, 4);
-      static const std::array<std::array<unsigned int, 3>, 4> table = {
+      static const ndarray<unsigned int, 4, 3> table = {
         {{{0, 1, 2}}, {{0, 1, 3}}, {{0, 2, 3}}, {{1, 2, 3}}}};
 
       return table[vertex];
@@ -677,13 +678,12 @@ ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
   else if (*this == ReferenceCells::Wedge)
     {
       AssertIndexRange(vertex, 6);
-      static const std::array<std::array<unsigned int, 3>, 6> table = {
-        {{{0, 2, 4}},
-         {{0, 2, 3}},
-         {{0, 3, 4}},
-         {{1, 2, 4}},
-         {{1, 2, 3}},
-         {{1, 3, 4}}}};
+      static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 4}},
+                                                         {{0, 2, 3}},
+                                                         {{0, 3, 4}},
+                                                         {{1, 2, 4}},
+                                                         {{1, 2, 3}},
+                                                         {{1, 3, 4}}}};
 
       return table[vertex];
     }
@@ -691,12 +691,11 @@ ReferenceCell::faces_for_given_vertex(const unsigned int vertex) const
     {
       AssertIndexRange(vertex, 5);
       static const unsigned int X = numbers::invalid_unsigned_int;
-      static const std::array<std::array<unsigned int, 4>, 5> table = {
-        {{{0, 1, 3, X}},
-         {{0, 2, 3, X}},
-         {{0, 1, 4, X}},
-         {{0, 2, 4, X}},
-         {{1, 2, 3, 4}}}};
+      static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 3, X}},
+                                                         {{0, 2, 3, X}},
+                                                         {{0, 1, 4, X}},
+                                                         {{0, 2, 4, X}},
+                                                         {{1, 2, 3, 4}}}};
 
       return {&table[vertex][0], vertex == 4 ? 4u : 3u};
     }
@@ -906,7 +905,7 @@ ReferenceCell::child_cell_on_face(const unsigned int face,
     }
   else if (*this == ReferenceCells::Triangle)
     {
-      static const std::array<std::array<unsigned int, 2>, 3> subcells = {
+      static const ndarray<unsigned int, 3, 2> subcells = {
         {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
       return subcells[face][subface];
@@ -954,7 +953,7 @@ ReferenceCell::standard_vertex_to_face_and_vertex_index(
     }
   else if (*this == ReferenceCells::Triangle)
     {
-      static const std::array<std::array<unsigned int, 2>, 3> table = {
+      static const ndarray<unsigned int, 3, 2> table = {
         {{{0, 0}}, {{0, 1}}, {{1, 1}}}};
 
       return table[vertex];
@@ -965,21 +964,21 @@ ReferenceCell::standard_vertex_to_face_and_vertex_index(
     }
   else if (*this == ReferenceCells::Tetrahedron)
     {
-      static const std::array<std::array<unsigned int, 2>, 4> table = {
+      static const ndarray<unsigned int, 4, 2> table = {
         {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{1, 2}}}};
 
       return table[vertex];
     }
   else if (*this == ReferenceCells::Pyramid)
     {
-      static const std::array<std::array<unsigned int, 2>, 5> table = {
+      static const ndarray<unsigned int, 5, 2> table = {
         {{{0, 0}}, {{0, 1}}, {{0, 2}}, {{0, 3}}, {{1, 2}}}};
 
       return table[vertex];
     }
   else if (*this == ReferenceCells::Wedge)
     {
-      static const std::array<std::array<unsigned int, 2>, 6> table = {
+      static const ndarray<unsigned int, 6, 2> table = {
         {{{0, 1}}, {{0, 0}}, {{0, 2}}, {{1, 0}}, {{1, 1}}, {{1, 2}}}};
 
       return table[vertex];
@@ -1098,7 +1097,7 @@ ReferenceCell::face_to_cell_lines(const unsigned int  face,
     }
   else if (*this == ReferenceCells::Tetrahedron)
     {
-      const static std::array<std::array<unsigned int, 3>, 4> table = {
+      const static ndarray<unsigned int, 4, 3> table = {
         {{{0, 1, 2}}, {{0, 3, 4}}, {{2, 5, 3}}, {{1, 4, 5}}}};
 
       return table[face]
@@ -1151,7 +1150,7 @@ ReferenceCell::face_to_cell_vertices(const unsigned int  face,
     }
   else if (*this == ReferenceCells::Triangle)
     {
-      static const std::array<std::array<unsigned int, 2>, 3> table = {
+      static const ndarray<unsigned int, 3, 2> table = {
         {{{0, 1}}, {{1, 2}}, {{2, 0}}}};
 
       return table[face][face_orientation ? vertex : (1 - vertex)];
@@ -1167,7 +1166,7 @@ ReferenceCell::face_to_cell_vertices(const unsigned int  face,
     }
   else if (*this == ReferenceCells::Tetrahedron)
     {
-      static const std::array<std::array<unsigned int, 3>, 4> table = {
+      static const ndarray<unsigned int, 4, 3> table = {
         {{{0, 1, 2}}, {{1, 0, 3}}, {{0, 2, 3}}, {{2, 1, 3}}}};
 
       return table[face][standard_to_real_face_vertex(
@@ -1176,12 +1175,11 @@ ReferenceCell::face_to_cell_vertices(const unsigned int  face,
   else if (*this == ReferenceCells::Pyramid)
     {
       constexpr auto X = numbers::invalid_unsigned_int;
-      static const std::array<std::array<unsigned int, 4>, 5> table = {
-        {{{0, 1, 2, 3}},
-         {{0, 2, 4, X}},
-         {{3, 1, 4, X}},
-         {{1, 0, 4, X}},
-         {{2, 3, 4, X}}}};
+      static const ndarray<unsigned int, 5, 4> table = {{{{0, 1, 2, 3}},
+                                                         {{0, 2, 4, X}},
+                                                         {{3, 1, 4, X}},
+                                                         {{1, 0, 4, X}},
+                                                         {{2, 3, 4, X}}}};
 
       return table[face][standard_to_real_face_vertex(
         vertex, face, face_orientation)];
@@ -1189,12 +1187,11 @@ ReferenceCell::face_to_cell_vertices(const unsigned int  face,
   else if (*this == ReferenceCells::Wedge)
     {
       constexpr auto X = numbers::invalid_unsigned_int;
-      static const std::array<std::array<unsigned int, 4>, 6> table = {
-        {{{1, 0, 2, X}},
-         {{3, 4, 5, X}},
-         {{0, 1, 3, 4}},
-         {{1, 2, 4, 5}},
-         {{2, 0, 5, 3}}}};
+      static const ndarray<unsigned int, 6, 4> table = {{{{1, 0, 2, X}},
+                                                         {{3, 4, 5, X}},
+                                                         {{0, 1, 3, 4}},
+                                                         {{1, 2, 4, 5}},
+                                                         {{2, 0, 5, 3}}}};
 
       return table[face][standard_to_real_face_vertex(
         vertex, face, face_orientation)];
@@ -1234,8 +1231,7 @@ ReferenceCell::standard_to_real_face_vertex(
     }
   else if (*this == ReferenceCells::Triangle)
     {
-      static const std::array<std::array<unsigned int, 2>, 2> table = {
-        {{{1, 0}}, {{0, 1}}}};
+      static const ndarray<unsigned int, 2, 2> table = {{{{1, 0}}, {{0, 1}}}};
 
       return table[face_orientation][vertex];
     }
@@ -1246,13 +1242,12 @@ ReferenceCell::standard_to_real_face_vertex(
     }
   else if (*this == ReferenceCells::Tetrahedron)
     {
-      static const std::array<std::array<unsigned int, 3>, 6> table = {
-        {{{0, 2, 1}},
-         {{0, 1, 2}},
-         {{2, 1, 0}},
-         {{1, 2, 0}},
-         {{1, 0, 2}},
-         {{2, 0, 1}}}};
+      static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
+                                                         {{0, 1, 2}},
+                                                         {{2, 1, 0}},
+                                                         {{1, 2, 0}},
+                                                         {{1, 0, 2}},
+                                                         {{2, 0, 1}}}};
 
       return table[face_orientation][vertex];
     }
@@ -1268,13 +1263,12 @@ ReferenceCell::standard_to_real_face_vertex(
         }
       else // One of the triangular faces
         {
-          static const std::array<std::array<unsigned int, 3>, 6> table = {
-            {{{0, 2, 1}},
-             {{0, 1, 2}},
-             {{2, 1, 0}},
-             {{1, 2, 0}},
-             {{1, 0, 2}},
-             {{2, 0, 1}}}};
+          static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
+                                                             {{0, 1, 2}},
+                                                             {{2, 1, 0}},
+                                                             {{1, 2, 0}},
+                                                             {{1, 0, 2}},
+                                                             {{2, 0, 1}}}};
 
           return table[face_orientation][vertex];
         }
@@ -1291,13 +1285,12 @@ ReferenceCell::standard_to_real_face_vertex(
         }
       else // One of the triangular faces
         {
-          static const std::array<std::array<unsigned int, 3>, 6> table = {
-            {{{0, 2, 1}},
-             {{0, 1, 2}},
-             {{2, 1, 0}},
-             {{1, 2, 0}},
-             {{1, 0, 2}},
-             {{2, 0, 1}}}};
+          static const ndarray<unsigned int, 6, 3> table = {{{{0, 2, 1}},
+                                                             {{0, 1, 2}},
+                                                             {{2, 1, 0}},
+                                                             {{1, 2, 0}},
+                                                             {{1, 0, 2}},
+                                                             {{2, 0, 1}}}};
 
           return table[face_orientation][vertex];
         }
@@ -1344,13 +1337,12 @@ ReferenceCell::standard_to_real_face_line(
     }
   else if (*this == ReferenceCells::Tetrahedron)
     {
-      static const std::array<std::array<unsigned int, 3>, 6> table = {
-        {{{2, 1, 0}},
-         {{0, 1, 2}},
-         {{1, 0, 2}},
-         {{1, 2, 0}},
-         {{0, 2, 1}},
-         {{2, 0, 1}}}};
+      static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
+                                                         {{0, 1, 2}},
+                                                         {{1, 0, 2}},
+                                                         {{1, 2, 0}},
+                                                         {{0, 2, 1}},
+                                                         {{2, 0, 1}}}};
 
       return table[face_orientation][line];
     }
@@ -1366,13 +1358,12 @@ ReferenceCell::standard_to_real_face_line(
         }
       else // One of the triangular faces
         {
-          static const std::array<std::array<unsigned int, 3>, 6> table = {
-            {{{2, 1, 0}},
-             {{0, 1, 2}},
-             {{1, 0, 2}},
-             {{1, 2, 0}},
-             {{0, 2, 1}},
-             {{2, 0, 1}}}};
+          static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
+                                                             {{0, 1, 2}},
+                                                             {{1, 0, 2}},
+                                                             {{1, 2, 0}},
+                                                             {{0, 2, 1}},
+                                                             {{2, 0, 1}}}};
 
           return table[face_orientation][line];
         }
@@ -1389,13 +1380,12 @@ ReferenceCell::standard_to_real_face_line(
         }
       else // One of the triangular faces
         {
-          static const std::array<std::array<unsigned int, 3>, 6> table = {
-            {{{2, 1, 0}},
-             {{0, 1, 2}},
-             {{1, 0, 2}},
-             {{1, 2, 0}},
-             {{0, 2, 1}},
-             {{2, 0, 1}}}};
+          static const ndarray<unsigned int, 6, 3> table = {{{{2, 1, 0}},
+                                                             {{0, 1, 2}},
+                                                             {{1, 0, 2}},
+                                                             {{1, 2, 0}},
+                                                             {{0, 2, 1}},
+                                                             {{2, 0, 1}}}};
 
           return table[face_orientation][line];
         }
@@ -1467,33 +1457,32 @@ ReferenceCell::n_vertices_to_type(const int dim, const unsigned int n_vertices)
   AssertIndexRange(dim, 4);
   AssertIndexRange(n_vertices, 9);
 
-  const auto X = ReferenceCells::Invalid;
-  static const std::array<std::array<ReferenceCell, 9>,
-                          4>
-    table = {{// dim 0
-              {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
-              // dim 1
-              {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
-              // dim 2
-              {{X,
-                X,
-                X,
-                ReferenceCells::Triangle,
-                ReferenceCells::Quadrilateral,
-                X,
-                X,
-                X,
-                X}},
-              // dim 3
-              {{X,
-                X,
-                X,
-                X,
-                ReferenceCells::Tetrahedron,
-                ReferenceCells::Pyramid,
-                ReferenceCells::Wedge,
-                X,
-                ReferenceCells::Hexahedron}}}};
+  const auto                                X     = ReferenceCells::Invalid;
+  static const ndarray<ReferenceCell, 4, 9> table = {
+    {// dim 0
+     {{X, ReferenceCells::Vertex, X, X, X, X, X, X, X}},
+     // dim 1
+     {{X, X, ReferenceCells::Line, X, X, X, X, X, X}},
+     // dim 2
+     {{X,
+       X,
+       X,
+       ReferenceCells::Triangle,
+       ReferenceCells::Quadrilateral,
+       X,
+       X,
+       X,
+       X}},
+     // dim 3
+     {{X,
+       X,
+       X,
+       X,
+       ReferenceCells::Tetrahedron,
+       ReferenceCells::Pyramid,
+       ReferenceCells::Wedge,
+       X,
+       ReferenceCells::Hexahedron}}}};
   Assert(table[dim][n_vertices] != ReferenceCells::Invalid,
          ExcMessage("The combination of dim = " + std::to_string(dim) +
                     " and n_vertices = " + std::to_string(n_vertices) +
@@ -1653,7 +1642,7 @@ ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
   else if (*this == ReferenceCells::Tetrahedron)
     {
       AssertIndexRange(face_no, 4);
-      static const std::array<std::array<Tensor<1, dim>, 2>, 4> table = {
+      static const ndarray<Tensor<1, dim>, 4, 2> table = {
         {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
          {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
          {{Point<dim>(0, 0, 1), Point<dim>(0, 1, 0)}},
@@ -1669,7 +1658,7 @@ ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
   else if (*this == ReferenceCells::Wedge)
     {
       AssertIndexRange(face_no, 5);
-      static const std::array<std::array<Tensor<1, dim>, 2>, 5> table = {
+      static const ndarray<Tensor<1, dim>, 5, 2> table = {
         {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
          {{Point<dim>(1, 0, 0), Point<dim>(0, 0, 1)}},
          {{Point<dim>(-1 / std::sqrt(2.0), +1 / std::sqrt(2.0), 0),
@@ -1682,7 +1671,7 @@ ReferenceCell::unit_tangential_vectors(const unsigned int face_no,
   else if (*this == ReferenceCells::Pyramid)
     {
       AssertIndexRange(face_no, 5);
-      static const std::array<std::array<Tensor<1, dim>, 2>, 5> table = {
+      static const ndarray<Tensor<1, dim>, 5, 2> table = {
         {{{Point<dim>(0, 1, 0), Point<dim>(1, 0, 0)}},
          {{Point<dim>(+1.0 / sqrt(2.0), 0, +1.0 / sqrt(2.0)),
            Point<dim>(0, 1, 0)}},
diff --git a/include/deal.II/meshworker/copy_data.h b/include/deal.II/meshworker/copy_data.h
index 8571ec401a..13be3a8ecb 100644
--- a/include/deal.II/meshworker/copy_data.h
+++ b/include/deal.II/meshworker/copy_data.h
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/ndarray.h>
 #include <deal.II/base/types.h>
 
 #include <deal.II/lac/full_matrix.h>
@@ -59,8 +60,8 @@ namespace MeshWorker
      * For every object, specify the size they should have.
      */
     explicit CopyData(
-      const std::array<std::array<unsigned int, 2>, n_matrices> &matrix_sizes,
-      const std::array<unsigned int, n_vectors> &                vector_sizes,
+      const ndarray<unsigned int, n_matrices, 2> &   matrix_sizes,
+      const std::array<unsigned int, n_vectors> &    vector_sizes,
       const std::array<unsigned int, n_dof_indices> &dof_indices_sizes);
 
     /**
@@ -107,8 +108,8 @@ namespace MeshWorker
 
   template <int n_matrices, int n_vectors, int n_dof_indices>
   CopyData<n_matrices, n_vectors, n_dof_indices>::CopyData(
-    const std::array<std::array<unsigned int, 2>, n_matrices> &matrix_sizes,
-    const std::array<unsigned int, n_vectors> &                vector_sizes,
+    const ndarray<unsigned int, n_matrices, 2> &   matrix_sizes,
+    const std::array<unsigned int, n_vectors> &    vector_sizes,
     const std::array<unsigned int, n_dof_indices> &dof_indices_sizes)
   {
     for (unsigned int i = 0; i < n_matrices; ++i)
diff --git a/source/base/geometry_info.cc b/source/base/geometry_info.cc
index 435cbd66f4..72c5d3820f 100644
--- a/source/base/geometry_info.cc
+++ b/source/base/geometry_info.cc
@@ -45,8 +45,7 @@ constexpr std::array<int, GeometryInfo<dim>::faces_per_cell>
   GeometryInfo<dim>::unit_normal_orientation;
 
 template <int dim>
-constexpr std::array<std::array<unsigned int, dim>,
-                     GeometryInfo<dim>::vertices_per_cell>
+constexpr ndarray<unsigned int, GeometryInfo<dim>::vertices_per_cell, dim>
   GeometryInfo<dim>::vertex_to_face;
 
 template <int dim>
@@ -58,8 +57,8 @@ constexpr std::array<Tensor<1, dim>, GeometryInfo<dim>::faces_per_cell>
   GeometryInfo<dim>::unit_normal_vector;
 
 template <int dim>
-constexpr std::array<std::array<Tensor<1, dim>, dim - 1>,
-                     GeometryInfo<dim>::faces_per_cell>
+constexpr ndarray<Tensor<1, dim>, GeometryInfo<dim>::faces_per_cell, dim - 1>
+
   GeometryInfo<dim>::unit_tangential_vectors;
 
 template <int dim>
diff --git a/source/base/tensor_product_polynomials.cc b/source/base/tensor_product_polynomials.cc
index 4917076f64..0dfb1dfb23 100644
--- a/source/base/tensor_product_polynomials.cc
+++ b/source/base/tensor_product_polynomials.cc
@@ -211,7 +211,7 @@ TensorProductPolynomials<dim, PolynomialType>::compute_grad(
   // uni-directional derivatives at
   // the given point in each
   // co-ordinate direction
-  std::array<std::array<double, 2>, dim> v;
+  ndarray<double, dim, 2> v;
   {
     std::vector<double> tmp(2);
     for (unsigned int d = 0; d < dim; ++d)
@@ -258,7 +258,7 @@ TensorProductPolynomials<dim, PolynomialType>::compute_grad_grad(
   std::array<unsigned int, dim> indices;
   compute_index(i, indices);
 
-  std::array<std::array<double, 3>, dim> v;
+  ndarray<double, dim, 3> v;
   {
     std::vector<double> tmp(3);
     for (unsigned int d = 0; d < dim; ++d)
@@ -353,8 +353,8 @@ TensorProductPolynomials<dim, PolynomialType>::evaluate(
   // coordinate direction. Once we have those values, we perform the
   // multiplications for the tensor product in the arbitrary dimension.
   const unsigned int n_polynomials = polynomials.size();
-  boost::container::small_vector<std::array<std::array<double, 5>, dim>, 20>
-    values_1d(n_polynomials);
+  boost::container::small_vector<ndarray<double, dim, 5>, 20> values_1d(
+    n_polynomials);
   if (n_values_and_derivatives == 1)
     for (unsigned int i = 0; i < n_polynomials; ++i)
       for (unsigned int d = 0; d < dim; ++d)
@@ -608,7 +608,7 @@ AnisotropicPolynomials<dim>::compute_grad(const unsigned int i,
   // uni-directional derivatives at
   // the given point in each
   // co-ordinate direction
-  std::array<std::array<double, 2>, dim> v;
+  ndarray<double, dim, 2> v;
   for (unsigned int d = 0; d < dim; ++d)
     polynomials[d][indices[d]].value(p(d), 1, v[d].data());
 
@@ -646,7 +646,7 @@ AnisotropicPolynomials<dim>::compute_grad_grad(const unsigned int i,
   std::array<unsigned int, dim> indices;
   compute_index(i, indices);
 
-  std::array<std::array<double, 3>, dim> v;
+  ndarray<double, dim, 3> v;
   for (unsigned int d = 0; d < dim; ++d)
     polynomials[d][indices[d]].value(p(d), 2, v[d].data());
 
diff --git a/source/fe/fe_p1nc.cc b/source/fe/fe_p1nc.cc
index 499ef20b50..babd3a2aa6 100644
--- a/source/fe/fe_p1nc.cc
+++ b/source/fe/fe_p1nc.cc
@@ -85,7 +85,7 @@ FE_P1NC::get_dpo_vector()
 
 
 
-std::array<std::array<double, 3>, 4>
+ndarray<double, 4, 3>
 FE_P1NC::get_linear_shape_coefficients(
   const Triangulation<2, 2>::cell_iterator &cell)
 {
@@ -102,7 +102,7 @@ FE_P1NC::get_linear_shape_coefficients(
   const double det = (mpt[0](0) - mpt[1](0)) * (mpt[2](1) - mpt[3](1)) -
                      (mpt[2](0) - mpt[3](0)) * (mpt[0](1) - mpt[1](1));
 
-  std::array<std::array<double, 3>, 4> coeffs;
+  ndarray<double, 4, 3> coeffs;
   coeffs[0][0] =
     ((mpt[2](1) - mpt[3](1)) * (0.5) - (mpt[0](1) - mpt[1](1)) * (0.5)) / det;
   coeffs[1][0] =
@@ -224,8 +224,7 @@ FE_P1NC::fill_fe_values(
   const unsigned int n_q_points = mapping_data.quadrature_points.size();
 
   // linear shape functions
-  std::array<std::array<double, 3>, 4> coeffs =
-    get_linear_shape_coefficients(cell);
+  ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
   // compute on the cell
   if (flags & update_values)
@@ -261,8 +260,7 @@ FE_P1NC::fill_fe_face_values(
   const UpdateFlags flags(fe_internal.update_each);
 
   // linear shape functions
-  const std::array<std::array<double, 3>, 4> coeffs =
-    get_linear_shape_coefficients(cell);
+  const ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
   // compute on the face
   const Quadrature<2> quadrature_on_face =
@@ -308,8 +306,7 @@ FE_P1NC::fill_fe_subface_values(
   const UpdateFlags flags(fe_internal.update_each);
 
   // linear shape functions
-  const std::array<std::array<double, 3>, 4> coeffs =
-    get_linear_shape_coefficients(cell);
+  const ndarray<double, 4, 3> coeffs = get_linear_shape_coefficients(cell);
 
   // compute on the subface
   const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(
diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc
index 197d3504ba..ce59fc0d2c 100644
--- a/source/grid/grid_generator.cc
+++ b/source/grid/grid_generator.cc
@@ -13,6 +13,8 @@
 //
 // ---------------------------------------------------------------------
 
+#include <deal.II/base/ndarray.h>
+
 #include <deal.II/distributed/fully_distributed_tria.h>
 #include <deal.II/distributed/shared_tria.h>
 #include <deal.II/distributed/tria.h>
@@ -7536,37 +7538,35 @@ namespace GridGenerator
      * A quadrilateral element is converted to 8 simplices elements. Each
      * triangle is defined by 3 vertices.
      */
-    static const std::array<std::array<unsigned int, 3>, 8> table_2D_cell = {
-      {{{0, 6, 4}},
-       {{8, 4, 6}},
-       {{8, 6, 5}},
-       {{1, 5, 6}},
-       {{2, 4, 7}},
-       {{8, 7, 4}},
-       {{8, 5, 7}},
-       {{3, 7, 5}}}};
+    static const ndarray<unsigned int, 8, 3> table_2D_cell = {{{{0, 6, 4}},
+                                                               {{8, 4, 6}},
+                                                               {{8, 6, 5}},
+                                                               {{1, 5, 6}},
+                                                               {{2, 4, 7}},
+                                                               {{8, 7, 4}},
+                                                               {{8, 5, 7}},
+                                                               {{3, 7, 5}}}};
 
     /* Cell definition 3D:
      * A hexahedron element is converted to 24 tetrahedron elements. Each
      * tetrahedron is defined by 4 vertices.
      */
-    static const std::array<std::array<unsigned int, 4>, 24>
-      vertex_ids_for_cells_3d = {
-        {{{0, 1, 12, 10}},  {{2, 3, 11, 12}},  {{7, 6, 11, 13}},
-         {{5, 4, 13, 10}},  {{0, 2, 8, 12}},   {{4, 6, 13, 8}},
-         {{5, 13, 7, 9}},   {{1, 9, 3, 12}},   {{0, 8, 4, 10}},
-         {{1, 5, 9, 10}},   {{3, 7, 11, 9}},   {{2, 6, 8, 11}},
-         {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
-         {{12, 13, 8, 10}}, {{13, 8, 10, 4}},  {{13, 10, 9, 5}},
-         {{13, 9, 11, 7}},  {{13, 11, 8, 6}},  {{10, 12, 9, 1}},
-         {{9, 12, 11, 3}},  {{11, 12, 8, 2}},  {{8, 12, 10, 0}}}};
+    static const ndarray<unsigned int, 24, 4> vertex_ids_for_cells_3d = {
+      {{{0, 1, 12, 10}},  {{2, 3, 11, 12}},  {{7, 6, 11, 13}},
+       {{5, 4, 13, 10}},  {{0, 2, 8, 12}},   {{4, 6, 13, 8}},
+       {{5, 13, 7, 9}},   {{1, 9, 3, 12}},   {{0, 8, 4, 10}},
+       {{1, 5, 9, 10}},   {{3, 7, 11, 9}},   {{2, 6, 8, 11}},
+       {{12, 13, 10, 9}}, {{12, 13, 9, 11}}, {{12, 13, 11, 8}},
+       {{12, 13, 8, 10}}, {{13, 8, 10, 4}},  {{13, 10, 9, 5}},
+       {{13, 9, 11, 7}},  {{13, 11, 8, 6}},  {{10, 12, 9, 1}},
+       {{9, 12, 11, 3}},  {{11, 12, 8, 2}},  {{8, 12, 10, 0}}}};
 
     /* Boundary-faces 2D:
      * After converting, each of the 4 quadrilateral faces is defined by faces
      * of 2 different triangles, i.e., lines. Note that lines are defined by 2
      * vertices.
      */
-    static const std::array<std::array<std::array<unsigned int, 2>, 2>, 4>
+    static const ndarray<unsigned int, 4, 2, 2>
       vertex_ids_for_boundary_faces_2d = {{{{{{0, 4}}, {{4, 2}}}},
                                            {{{{1, 5}}, {{5, 3}}}},
                                            {{{{0, 6}}, {{6, 1}}}},
@@ -7577,7 +7577,7 @@ namespace GridGenerator
      * 4 different tetrahedron faces, i.e., triangles. Note that a triangle is
      * defined by 3 vertices.
      */
-    static const std::array<std::array<std::array<unsigned int, 3>, 4>, 6>
+    static const ndarray<unsigned int, 6, 4, 3>
       vertex_ids_for_boundary_faces_3d = {
         {{{{{0, 4, 8}}, {{4, 8, 6}}, {{8, 6, 2}}, {{0, 2, 8}}}},
          {{{{1, 3, 9}}, {{3, 9, 7}}, {{9, 7, 5}}, {{1, 9, 5}}}},
@@ -7590,66 +7590,62 @@ namespace GridGenerator
      * The converted triangulation based on simplices has 8 faces that do not
      * form the boundary, i.e. inner-faces, each defined by 2 vertices.
      */
-    static const std::array<std::array<unsigned int, 2>, 8>
-      vertex_ids_for_inner_faces_2d = {{{{6, 4}},
-                                        {{6, 8}},
-                                        {{6, 5}},
-                                        {{4, 8}},
-                                        {{8, 5}},
-                                        {{7, 4}},
-                                        {{7, 8}},
-                                        {{7, 5}}}};
+    static const ndarray<unsigned int, 8, 2> vertex_ids_for_inner_faces_2d = {
+      {{{6, 4}},
+       {{6, 8}},
+       {{6, 5}},
+       {{4, 8}},
+       {{8, 5}},
+       {{7, 4}},
+       {{7, 8}},
+       {{7, 5}}}};
 
     /* Inner-faces 3D:
      * The converted triangulation based on simplices has 72 faces that do not
      * form the boundary, i.e. inner-faces, each defined by 3 vertices.
      */
-    static const std::array<std::array<unsigned int, 3>, 72>
-      vertex_ids_for_inner_faces_3d = {
-        {{{0, 12, 10}},  {{12, 1, 10}},  {{12, 1, 9}},  {{12, 3, 9}},
-         {{12, 2, 11}},  {{12, 3, 11}},  {{12, 0, 8}},  {{12, 2, 8}},
-         {{9, 13, 5}},   {{13, 7, 9}},   {{11, 7, 13}}, {{11, 6, 13}},
-         {{4, 8, 13}},   {{6, 8, 13}},   {{4, 13, 10}}, {{13, 5, 10}},
-         {{10, 9, 5}},   {{10, 9, 1}},   {{11, 9, 7}},  {{11, 9, 3}},
-         {{8, 11, 2}},   {{8, 11, 6}},   {{8, 10, 0}},  {{8, 10, 4}},
-         {{12, 3, 9}},   {{12, 9, 11}},  {{12, 3, 11}}, {{3, 9, 11}},
-         {{2, 12, 8}},   {{2, 12, 11}},  {{2, 11, 8}},  {{8, 12, 11}},
-         {{0, 12, 10}},  {{0, 12, 8}},   {{0, 8, 10}},  {{8, 10, 12}},
-         {{12, 1, 10}},  {{12, 1, 9}},   {{1, 10, 9}},  {{10, 9, 12}},
-         {{10, 8, 4}},   {{10, 8, 13}},  {{4, 13, 8}},  {{4, 13, 10}},
-         {{10, 9, 13}},  {{10, 9, 5}},   {{13, 5, 10}}, {{13, 5, 9}},
-         {{13, 7, 9}},   {{13, 7, 11}},  {{9, 11, 13}}, {{9, 11, 7}},
-         {{8, 11, 13}},  {{8, 11, 6}},   {{6, 13, 8}},  {{6, 13, 11}},
-         {{12, 13, 10}}, {{12, 13, 8}},  {{8, 10, 13}}, {{8, 10, 12}},
-         {{12, 13, 10}}, {{12, 13, 9}},  {{10, 9, 13}}, {{10, 9, 12}},
-         {{12, 13, 9}},  {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
-         {{12, 13, 11}}, {{12, 13, 8}},  {{8, 11, 13}}, {{8, 11, 12}}}};
+    static const ndarray<unsigned int, 72, 3> vertex_ids_for_inner_faces_3d = {
+      {{{0, 12, 10}},  {{12, 1, 10}},  {{12, 1, 9}},  {{12, 3, 9}},
+       {{12, 2, 11}},  {{12, 3, 11}},  {{12, 0, 8}},  {{12, 2, 8}},
+       {{9, 13, 5}},   {{13, 7, 9}},   {{11, 7, 13}}, {{11, 6, 13}},
+       {{4, 8, 13}},   {{6, 8, 13}},   {{4, 13, 10}}, {{13, 5, 10}},
+       {{10, 9, 5}},   {{10, 9, 1}},   {{11, 9, 7}},  {{11, 9, 3}},
+       {{8, 11, 2}},   {{8, 11, 6}},   {{8, 10, 0}},  {{8, 10, 4}},
+       {{12, 3, 9}},   {{12, 9, 11}},  {{12, 3, 11}}, {{3, 9, 11}},
+       {{2, 12, 8}},   {{2, 12, 11}},  {{2, 11, 8}},  {{8, 12, 11}},
+       {{0, 12, 10}},  {{0, 12, 8}},   {{0, 8, 10}},  {{8, 10, 12}},
+       {{12, 1, 10}},  {{12, 1, 9}},   {{1, 10, 9}},  {{10, 9, 12}},
+       {{10, 8, 4}},   {{10, 8, 13}},  {{4, 13, 8}},  {{4, 13, 10}},
+       {{10, 9, 13}},  {{10, 9, 5}},   {{13, 5, 10}}, {{13, 5, 9}},
+       {{13, 7, 9}},   {{13, 7, 11}},  {{9, 11, 13}}, {{9, 11, 7}},
+       {{8, 11, 13}},  {{8, 11, 6}},   {{6, 13, 8}},  {{6, 13, 11}},
+       {{12, 13, 10}}, {{12, 13, 8}},  {{8, 10, 13}}, {{8, 10, 12}},
+       {{12, 13, 10}}, {{12, 13, 9}},  {{10, 9, 13}}, {{10, 9, 12}},
+       {{12, 13, 9}},  {{12, 13, 11}}, {{9, 11, 13}}, {{9, 11, 12}},
+       {{12, 13, 11}}, {{12, 13, 8}},  {{8, 11, 13}}, {{8, 11, 12}}}};
 
     /* Inner-edges 3D:
      * The converted triangulation based on simplices has 60 edges that do not
      * coincide with the boundary, i.e. inner-edges, each defined by 2 vertices.
      */
-    static const std::array<std::array<unsigned int, 2>, 60>
-      vertex_ids_for_inner_edges_3d = {
-        {{{12, 10}}, {{12, 9}},  {{12, 11}}, {{12, 8}},  {{9, 13}},
-         {{11, 13}}, {{8, 13}},  {{10, 13}}, {{10, 9}},  {{9, 11}},
-         {{11, 8}},  {{8, 10}},  {{12, 9}},  {{12, 11}}, {{11, 9}},
-         {{12, 8}},  {{12, 11}}, {{11, 8}},  {{12, 8}},  {{12, 10}},
-         {{10, 8}},  {{12, 10}}, {{12, 9}},  {{9, 10}},  {{13, 10}},
-         {{13, 8}},  {{8, 10}},  {{13, 10}}, {{13, 9}},  {{9, 10}},
-         {{13, 11}}, {{13, 9}},  {{11, 9}},  {{13, 11}}, {{13, 8}},
-         {{11, 8}},  {{12, 13}}, {{8, 10}},  {{8, 13}},  {{10, 13}},
-         {{8, 12}},  {{10, 12}}, {{12, 13}}, {{10, 9}},  {{10, 13}},
-         {{9, 13}},  {{10, 12}}, {{9, 12}},  {{12, 13}}, {{9, 11}},
-         {{9, 13}},  {{11, 13}}, {{9, 12}},  {{11, 12}}, {{12, 13}},
-         {{11, 8}},  {{11, 13}}, {{8, 13}},  {{11, 12}}, {{8, 12}}}};
+    static const ndarray<unsigned int, 60, 2> vertex_ids_for_inner_edges_3d = {
+      {{{12, 10}}, {{12, 9}},  {{12, 11}}, {{12, 8}},  {{9, 13}},  {{11, 13}},
+       {{8, 13}},  {{10, 13}}, {{10, 9}},  {{9, 11}},  {{11, 8}},  {{8, 10}},
+       {{12, 9}},  {{12, 11}}, {{11, 9}},  {{12, 8}},  {{12, 11}}, {{11, 8}},
+       {{12, 8}},  {{12, 10}}, {{10, 8}},  {{12, 10}}, {{12, 9}},  {{9, 10}},
+       {{13, 10}}, {{13, 8}},  {{8, 10}},  {{13, 10}}, {{13, 9}},  {{9, 10}},
+       {{13, 11}}, {{13, 9}},  {{11, 9}},  {{13, 11}}, {{13, 8}},  {{11, 8}},
+       {{12, 13}}, {{8, 10}},  {{8, 13}},  {{10, 13}}, {{8, 12}},  {{10, 12}},
+       {{12, 13}}, {{10, 9}},  {{10, 13}}, {{9, 13}},  {{10, 12}}, {{9, 12}},
+       {{12, 13}}, {{9, 11}},  {{9, 13}},  {{11, 13}}, {{9, 12}},  {{11, 12}},
+       {{12, 13}}, {{11, 8}},  {{11, 13}}, {{8, 13}},  {{11, 12}}, {{8, 12}}}};
 
     /* Boundary-edges 3D:
      * For each of the 6 boundary-faces of the hexahedron, there are 8 edges (of
      * different tetrahedrons) that coincide with the boundary, i.e.
      * boundary-edges. Each boundary-edge is defined by 2 vertices.
      */
-    static const std::array<std::array<std::array<unsigned int, 2>, 8>, 6>
+    static const ndarray<unsigned int, 6, 8, 2>
       vertex_ids_for_boundary_edges_3d = {{{{{{4, 6}},
                                              {{4, 8}},
                                              {{6, 8}},
diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc
index a7fcd3dc7f..4808d45584 100644
--- a/source/simplex/polynomials.cc
+++ b/source/simplex/polynomials.cc
@@ -14,6 +14,8 @@
 // ---------------------------------------------------------------------
 
 
+#include <deal.II/base/ndarray.h>
+
 #include <deal.II/simplex/barycentric_polynomials.h>
 #include <deal.II/simplex/polynomials.h>
 
@@ -74,34 +76,33 @@ namespace Simplex
      * to access the right shape function within the triangle and and within
      * the line.
      */
-    static const constexpr std::array<std::array<unsigned int, 2>, 6>
-      wedge_table_1{
-        {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
+    static const constexpr ndarray<unsigned int, 6, 2> wedge_table_1{
+      {{{0, 0}}, {{1, 0}}, {{2, 0}}, {{0, 1}}, {{1, 1}}, {{2, 1}}}};
 
     /**
      * Decompose the shape-function index of a quadratic wedge into an index
      * to access the right shape function within the triangle and and within
      * the line.
      */
-    static const constexpr std::array<std::array<unsigned int, 2>, 18>
-      wedge_table_2{{{{0, 0}},
-                     {{1, 0}},
-                     {{2, 0}},
-                     {{0, 1}},
-                     {{1, 1}},
-                     {{2, 1}},
-                     {{3, 0}},
-                     {{4, 0}},
-                     {{5, 0}},
-                     {{3, 1}},
-                     {{4, 1}},
-                     {{5, 1}},
-                     {{0, 2}},
-                     {{1, 2}},
-                     {{2, 2}},
-                     {{3, 2}},
-                     {{4, 2}},
-                     {{5, 2}}}};
+    static const constexpr ndarray<unsigned int, 18, 2> wedge_table_2{
+      {{{0, 0}},
+       {{1, 0}},
+       {{2, 0}},
+       {{0, 1}},
+       {{1, 1}},
+       {{2, 1}},
+       {{3, 0}},
+       {{4, 0}},
+       {{5, 0}},
+       {{3, 1}},
+       {{4, 1}},
+       {{5, 1}},
+       {{0, 2}},
+       {{1, 2}},
+       {{2, 2}},
+       {{3, 2}},
+       {{4, 2}},
+       {{5, 2}}}};
   } // namespace
 
 
-- 
2.39.5