From: David Wells <drwells@email.unc.edu>
Date: Tue, 13 Aug 2024 22:04:07 +0000 (-0600)
Subject: TensorProductPolynomialsBubbles: remove some reinterpret_casts.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d6bcb6880d149af2c89d72dbe4e90f9fefe392ee;p=dealii.git

TensorProductPolynomialsBubbles: remove some reinterpret_casts.
---

diff --git a/include/deal.II/base/tensor_product_polynomials_bubbles.h b/include/deal.II/base/tensor_product_polynomials_bubbles.h
index f2a8396978..96083ac2a1 100644
--- a/include/deal.II/base/tensor_product_polynomials_bubbles.h
+++ b/include/deal.II/base/tensor_product_polynomials_bubbles.h
@@ -326,136 +326,129 @@ TensorProductPolynomialsBubbles<dim>::compute_derivative(
   if (i < max_q_indices)
     return tensor_polys.template compute_derivative<order>(i, p);
 
-  const unsigned int comp = i - tensor_polys.n();
+  [[maybe_unused]] const unsigned int comp = i - tensor_polys.n();
 
-  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_1[d] = 1.;
-              // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
-              for (unsigned j = 0; j < dim; ++j)
-                derivative_1[d] *=
-                  (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j]));
-              // and multiply with (2*x_i-1)^{r-1}
-              for (unsigned int i = 0; i < q_degree - 1; ++i)
-                derivative_1[d] *= 2 * p[comp] - 1;
-            }
-
-          if (q_degree >= 2)
-            {
-              // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
-              double value = 1.;
-              for (unsigned int j = 0; j < dim; ++j)
-                value *= 4 * p[j] * (1 - p[j]);
-              // and multiply with grad(2*x_i-1)^{r-1}
-              double tmp = value * 2 * (q_degree - 1);
-              for (unsigned int i = 0; i < q_degree - 2; ++i)
-                tmp *= 2 * p[comp] - 1;
-              derivative_1[comp] += tmp;
-            }
-
-          return derivative;
+          derivative[d] = 1.;
+          // compute grad(4*\prod_{i=1}^d (x_i(1-x_i)))(p)
+          for (unsigned j = 0; j < dim; ++j)
+            derivative[d] *=
+              (d == j ? 4 * (1 - 2 * p[j]) : 4 * p[j] * (1 - p[j]));
+          // and multiply with (2*x_i-1)^{r-1}
+          for (unsigned int i = 0; i < q_degree - 1; ++i)
+            derivative[d] *= 2 * p[comp] - 1;
         }
-      case 2:
+
+      if (q_degree >= 2)
         {
-          Tensor<2, dim> &derivative_2 =
-            *reinterpret_cast<Tensor<2, dim> *>(&derivative);
+          // add \prod_{i=1}^d 4*(x_i(1-x_i))(p)
+          double value = 1.;
+          for (unsigned int j = 0; j < dim; ++j)
+            value *= 4 * p[j] * (1 - p[j]);
+          // and multiply with grad(2*x_i-1)^{r-1}
+          double tmp = value * 2 * (q_degree - 1);
+          for (unsigned int i = 0; i < q_degree - 2; ++i)
+            tmp *= 2 * p[comp] - 1;
+          derivative[comp] += tmp;
+        }
 
-          double v[dim + 1][3];
-          {
-            for (unsigned int c = 0; c < dim; ++c)
-              {
-                v[c][0] = 4 * p[c] * (1 - p[c]);
-                v[c][1] = 4 * (1 - 2 * p[c]);
-                v[c][2] = -8;
-              }
+      return derivative;
+    }
+  else if constexpr (order == 2)
+    {
+      Tensor<2, dim> derivative;
 
-            double tmp = 1.;
-            for (unsigned int i = 0; i < q_degree - 1; ++i)
-              tmp *= 2 * p[comp] - 1;
-            v[dim][0] = tmp;
+      double v[dim + 1][3];
+      {
+        for (unsigned int c = 0; c < dim; ++c)
+          {
+            v[c][0] = 4 * p[c] * (1 - p[c]);
+            v[c][1] = 4 * (1 - 2 * p[c]);
+            v[c][2] = -8;
+          }
 
-            if (q_degree >= 2)
-              {
-                double tmp = 2 * (q_degree - 1);
-                for (unsigned int i = 0; i < q_degree - 2; ++i)
-                  tmp *= 2 * p[comp] - 1;
-                v[dim][1] = tmp;
-              }
-            else
-              v[dim][1] = 0.;
+        double tmp = 1.;
+        for (unsigned int i = 0; i < q_degree - 1; ++i)
+          tmp *= 2 * p[comp] - 1;
+        v[dim][0] = tmp;
 
-            if (q_degree >= 3)
-              {
-                double tmp = 4 * (q_degree - 2) * (q_degree - 1);
-                for (unsigned int i = 0; i < q_degree - 3; ++i)
-                  tmp *= 2 * p[comp] - 1;
-                v[dim][2] = tmp;
-              }
-            else
-              v[dim][2] = 0.;
+        if (q_degree >= 2)
+          {
+            double tmp = 2 * (q_degree - 1);
+            for (unsigned int i = 0; i < q_degree - 2; ++i)
+              tmp *= 2 * p[comp] - 1;
+            v[dim][1] = tmp;
           }
+        else
+          v[dim][1] = 0.;
 
-          // calculate (\partial_j \partial_k \psi) * monomial
-          Tensor<2, dim> grad_grad_1;
-          for (unsigned int d1 = 0; d1 < dim; ++d1)
-            for (unsigned int d2 = 0; d2 < dim; ++d2)
+        if (q_degree >= 3)
+          {
+            double tmp = 4 * (q_degree - 2) * (q_degree - 1);
+            for (unsigned int i = 0; i < q_degree - 3; ++i)
+              tmp *= 2 * p[comp] - 1;
+            v[dim][2] = tmp;
+          }
+        else
+          v[dim][2] = 0.;
+      }
+
+      // calculate (\partial_j \partial_k \psi) * monomial
+      Tensor<2, dim> grad_grad_1;
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          {
+            grad_grad_1[d1][d2] = v[dim][0];
+            for (unsigned int x = 0; x < dim; ++x)
               {
-                grad_grad_1[d1][d2] = v[dim][0];
-                for (unsigned int x = 0; x < dim; ++x)
+                unsigned int derivative = 0;
+                if (d1 == x || d2 == x)
                   {
-                    unsigned int derivative = 0;
-                    if (d1 == x || d2 == x)
-                      {
-                        if (d1 == d2)
-                          derivative = 2;
-                        else
-                          derivative = 1;
-                      }
-                    grad_grad_1[d1][d2] *= v[x][derivative];
+                    if (d1 == d2)
+                      derivative = 2;
+                    else
+                      derivative = 1;
                   }
+                grad_grad_1[d1][d2] *= v[x][derivative];
               }
+          }
 
-          // calculate (\partial_j  \psi) *(\partial_k monomial)
-          // and (\partial_k  \psi) *(\partial_j monomial)
-          Tensor<2, dim> grad_grad_2;
-          Tensor<2, dim> grad_grad_3;
-          for (unsigned int d = 0; d < dim; ++d)
+      // calculate (\partial_j  \psi) *(\partial_k monomial)
+      // and (\partial_k  \psi) *(\partial_j monomial)
+      Tensor<2, dim> grad_grad_2;
+      Tensor<2, dim> grad_grad_3;
+      for (unsigned int d = 0; d < dim; ++d)
+        {
+          grad_grad_2[d][comp] = v[dim][1];
+          grad_grad_3[comp][d] = v[dim][1];
+          for (unsigned int x = 0; x < dim; ++x)
             {
-              grad_grad_2[d][comp] = v[dim][1];
-              grad_grad_3[comp][d] = v[dim][1];
-              for (unsigned int x = 0; x < dim; ++x)
-                {
-                  grad_grad_2[d][comp] *= v[x][d == x];
-                  grad_grad_3[comp][d] *= v[x][d == x];
-                }
+              grad_grad_2[d][comp] *= v[x][d == x];
+              grad_grad_3[comp][d] *= v[x][d == x];
             }
+        }
 
-          // calculate \psi *(\partial j \partial_k monomial) and sum
-          double psi_value = 1.;
-          for (unsigned int x = 0; x < dim; ++x)
-            psi_value *= v[x][0];
+      // calculate \psi *(\partial j \partial_k monomial) and sum
+      double psi_value = 1.;
+      for (unsigned int x = 0; x < dim; ++x)
+        psi_value *= v[x][0];
 
-          for (unsigned int d1 = 0; d1 < dim; ++d1)
-            for (unsigned int d2 = 0; d2 < dim; ++d2)
-              derivative_2[d1][d2] =
-                grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
-          derivative_2[comp][comp] += psi_value * v[dim][2];
+      for (unsigned int d1 = 0; d1 < dim; ++d1)
+        for (unsigned int d2 = 0; d2 < dim; ++d2)
+          derivative[d1][d2] =
+            grad_grad_1[d1][d2] + grad_grad_2[d1][d2] + grad_grad_3[d1][d2];
+      derivative[comp][comp] += psi_value * v[dim][2];
 
-          return derivative;
-        }
-      default:
-        {
-          DEAL_II_NOT_IMPLEMENTED();
-          return derivative;
-        }
+      return derivative;
+    }
+  else
+    {
+      DEAL_II_NOT_IMPLEMENTED();
+      return {};
     }
 }