From c217eb4fffaa16fc549588a814d069be1130b3c9 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Wed, 14 Aug 2024 16:43:57 -0600
Subject: [PATCH] TensorProductPolynomials: remove some reinterpret_casts.

---
 .../deal.II/base/tensor_product_polynomials.h | 150 +++++++++---------
 1 file changed, 71 insertions(+), 79 deletions(-)

diff --git a/include/deal.II/base/tensor_product_polynomials.h b/include/deal.II/base/tensor_product_polynomials.h
index c535fefbfd..86e4b962f3 100644
--- a/include/deal.II/base/tensor_product_polynomials.h
+++ b/include/deal.II/base/tensor_product_polynomials.h
@@ -761,36 +761,78 @@ AnisotropicPolynomials<dim>::compute_derivative(const unsigned int i,
   for (unsigned int d = 0; d < dim; ++d)
     polynomials[d][indices[d]].value(p[d], v[d]);
 
-  Tensor<order, dim> derivative;
-  switch (order)
+  if constexpr (order == 1)
     {
-      case 1:
+      Tensor<1, dim> derivative;
+      for (unsigned int d = 0; d < dim; ++d)
         {
-          Tensor<1, dim> &derivative_1 =
-            *reinterpret_cast<Tensor<1, dim> *>(&derivative);
-          for (unsigned int d = 0; d < dim; ++d)
+          derivative[d] = 1.;
+          for (unsigned int x = 0; x < dim; ++x)
+            {
+              unsigned int x_order = 0;
+              if (d == x)
+                ++x_order;
+
+              derivative[d] *= v[x][x_order];
+            }
+        }
+
+      return derivative;
+    }
+  else if constexpr (order == 2)
+    {
+      Tensor<2, dim> derivative;
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          {
+            derivative[d1][d2] = 1.;
+            for (unsigned int x = 0; x < dim; ++x)
+              {
+                unsigned int x_order = 0;
+                if (d1 == x)
+                  ++x_order;
+                if (d2 == x)
+                  ++x_order;
+
+                derivative[d1][d2] *= v[x][x_order];
+              }
+          }
+
+      return derivative;
+    }
+  else if constexpr (order == 3)
+    {
+      Tensor<3, dim> derivative;
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          for (unsigned int d3 = 0; d3 < dim; ++d3)
             {
-              derivative_1[d] = 1.;
+              derivative[d1][d2][d3] = 1.;
               for (unsigned int x = 0; x < dim; ++x)
                 {
                   unsigned int x_order = 0;
-                  if (d == x)
+                  if (d1 == x)
+                    ++x_order;
+                  if (d2 == x)
+                    ++x_order;
+                  if (d3 == x)
                     ++x_order;
 
-                  derivative_1[d] *= v[x][x_order];
+                  derivative[d1][d2][d3] *= v[x][x_order];
                 }
             }
 
-          return derivative;
-        }
-      case 2:
-        {
-          Tensor<2, dim> &derivative_2 =
-            *reinterpret_cast<Tensor<2, dim> *>(&derivative);
-          for (unsigned int d1 = 0; d1 < dim; ++d1)
-            for (unsigned int d2 = 0; d2 < dim; ++d2)
+      return derivative;
+    }
+  else if constexpr (order == 4)
+    {
+      Tensor<4, dim> derivative;
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          for (unsigned int d3 = 0; d3 < dim; ++d3)
+            for (unsigned int d4 = 0; d4 < dim; ++d4)
               {
-                derivative_2[d1][d2] = 1.;
+                derivative[d1][d2][d3][d4] = 1.;
                 for (unsigned int x = 0; x < dim; ++x)
                   {
                     unsigned int x_order = 0;
@@ -798,71 +840,21 @@ AnisotropicPolynomials<dim>::compute_derivative(const unsigned int i,
                       ++x_order;
                     if (d2 == x)
                       ++x_order;
+                    if (d3 == x)
+                      ++x_order;
+                    if (d4 == x)
+                      ++x_order;
 
-                    derivative_2[d1][d2] *= v[x][x_order];
+                    derivative[d1][d2][d3][d4] *= v[x][x_order];
                   }
               }
 
-          return derivative;
-        }
-      case 3:
-        {
-          Tensor<3, dim> &derivative_3 =
-            *reinterpret_cast<Tensor<3, dim> *>(&derivative);
-          for (unsigned int d1 = 0; d1 < dim; ++d1)
-            for (unsigned int d2 = 0; d2 < dim; ++d2)
-              for (unsigned int d3 = 0; d3 < dim; ++d3)
-                {
-                  derivative_3[d1][d2][d3] = 1.;
-                  for (unsigned int x = 0; x < dim; ++x)
-                    {
-                      unsigned int x_order = 0;
-                      if (d1 == x)
-                        ++x_order;
-                      if (d2 == x)
-                        ++x_order;
-                      if (d3 == x)
-                        ++x_order;
-
-                      derivative_3[d1][d2][d3] *= v[x][x_order];
-                    }
-                }
-
-          return derivative;
-        }
-      case 4:
-        {
-          Tensor<4, dim> &derivative_4 =
-            *reinterpret_cast<Tensor<4, dim> *>(&derivative);
-          for (unsigned int d1 = 0; d1 < dim; ++d1)
-            for (unsigned int d2 = 0; d2 < dim; ++d2)
-              for (unsigned int d3 = 0; d3 < dim; ++d3)
-                for (unsigned int d4 = 0; d4 < dim; ++d4)
-                  {
-                    derivative_4[d1][d2][d3][d4] = 1.;
-                    for (unsigned int x = 0; x < dim; ++x)
-                      {
-                        unsigned int x_order = 0;
-                        if (d1 == x)
-                          ++x_order;
-                        if (d2 == x)
-                          ++x_order;
-                        if (d3 == x)
-                          ++x_order;
-                        if (d4 == x)
-                          ++x_order;
-
-                        derivative_4[d1][d2][d3][d4] *= v[x][x_order];
-                      }
-                  }
-
-          return derivative;
-        }
-      default:
-        {
-          DEAL_II_NOT_IMPLEMENTED();
-          return derivative;
-        }
+      return derivative;
+    }
+  else
+    {
+      DEAL_II_NOT_IMPLEMENTED();
+      return {};
     }
 }
 
-- 
2.39.5