basis - add CeedBasisApplyAtPoints transpose mode

CEED · Jul 12, 2023 · 2a94f45 · 2a94f45
1 parent c8c3fa7
commit 2a94f45
Show file tree

Hide file tree

Showing 6 changed files with 253 additions and 5 deletions.
diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
@@ -1568,8 +1568,47 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
       CeedCall(CeedVectorRestoreArray(v, &v_array));
       break;
     }
-    case CEED_TRANSPOSE:
-      return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, "CEED_TRANSPOSE unsupported for arbitrary basis point evaluation");
+    case CEED_TRANSPOSE: {
+      // Arbitrary points to nodes
+      CeedScalar       *chebyshev_coeffs;
+      const CeedScalar *u_array, *x_array_read;
+
+      // -- Transpose of evaluaton of Chebyshev polynomials at arbitrary points
+      CeedCall(CeedVectorGetArrayWrite(basis->vec_chebyshev, CEED_MEM_HOST, &chebyshev_coeffs));
+      CeedCall(CeedVectorGetArrayRead(x_ref, CEED_MEM_HOST, &x_array_read));
+      CeedCall(CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array));
+      {
+        CeedScalar tmp[2][num_comp * CeedIntPow(Q_1d, dim)], chebyshev_x[Q_1d];
+
+        // ---- Values at point
+        for (CeedInt p = 0; p < num_points; p++) {
+          CeedInt pre = num_comp * 1, post = 1;
+
+          for (CeedInt d = dim - 1; d >= 0; d--) {
+            // ------ Compute Chebyshev polynomial values
+            {
+              const CeedScalar x = x_array_read[p * dim + d];
+
+              chebyshev_x[0] = 1.0;
+              chebyshev_x[1] = 2 * x;
+              for (CeedInt j = 2; j < Q_1d; j++) chebyshev_x[j] = 2 * x * chebyshev_x[j - 1] - chebyshev_x[j - 2];
+            }
+            // ------ Tensor contract
+            CeedCall(CeedTensorContractApply(basis->contract, pre, 1, post, Q_1d, chebyshev_x, t_mode, p > 0 && d == 0,
+                                             d == (dim - 1) ? &u_array[p * num_comp] : tmp[d % 2], d == 0 ? chebyshev_coeffs : tmp[(d + 1) % 2]));
+            pre /= 1;
+            post *= Q_1d;
+          }
+        }
+      }
+      CeedCall(CeedVectorRestoreArray(basis->vec_chebyshev, &chebyshev_coeffs));
+      CeedCall(CeedVectorRestoreArrayRead(x_ref, &x_array_read));
+      CeedCall(CeedVectorRestoreArrayRead(u, &u_array));
+
+      // -- Interpolate transpose from Chebyshev coefficients
+      CeedCall(CeedBasisApply(basis->basis_chebyshev, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, basis->vec_chebyshev, v));
+      break;
+    }
   }
 
   return CEED_ERROR_SUCCESS;

diff --git a/tests/t350-basis.c b/tests/t350-basis.c
@@ -52,7 +52,7 @@ int main(int argc, char **argv) {
   // Interpolate to arbitrary points
   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
   {
-    CeedScalar x_array[4] = {-0.33, -0.65, .16, 0.99};
+    CeedScalar x_array[4] = {-0.33, -0.65, 0.16, 0.99};
 
     CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
   }

diff --git a/tests/t351-basis.c b/tests/t351-basis.c
@@ -61,7 +61,7 @@ int main(int argc, char **argv) {
     // Interpolate to arbitrary points
     CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
     {
-      CeedScalar x_array[12] = {-0.33, -0.65, .16, 0.99, -0.65, .16, 0.99, -0.33, .16, 0.99, -0.33, -0.65};
+      CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65};
 
       CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }

diff --git a/tests/t352-basis.c b/tests/t352-basis.c
@@ -61,7 +61,7 @@ int main(int argc, char **argv) {
     // Interpolate to arbitrary points
     CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp, p, q, CEED_GAUSS, &basis_u);
     {
-      CeedScalar x_array[12] = {-0.33, -0.65, .16, 0.99, -0.65, .16, 0.99, -0.33, .16, 0.99, -0.33, -0.65};
+      CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65};
 
       CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
     }

diff --git a/tests/t353-basis.c b/tests/t353-basis.c
@@ -0,0 +1,95 @@
+/// @file
+/// Test polynomial interpolation transpose from arbirtary points in 1D
+/// \test Test polynomial interpolation transpose from arbitrary points in 1D
+#include <ceed.h>
+#include <math.h>
+#include <stdio.h>
+
+#define ALEN(a) (sizeof(a) / sizeof((a)[0]))
+
+static CeedScalar Eval(CeedScalar x, CeedInt n, const CeedScalar *c) {
+  CeedScalar y = c[n - 1];
+  for (CeedInt i = n - 2; i >= 0; i--) y = y * x + c[i];
+  return y;
+}
+
+int main(int argc, char **argv) {
+  Ceed             ceed;
+  CeedVector       x, x_nodes, x_points, x_point, u, v, u_point, v_point;
+  CeedBasis        basis_x, basis_u;
+  const CeedInt    p = 5, q = 5, num_points = 4;
+  const CeedScalar c[4] = {1, 2, 3, 4};  // 1 + 2x + 3x^2 + ...
+
+  CeedInit(argv[1], &ceed);
+
+  CeedVectorCreate(ceed, 2, &x);
+  CeedVectorCreate(ceed, p, &x_nodes);
+  CeedVectorCreate(ceed, num_points, &x_points);
+  CeedVectorCreate(ceed, 1, &x_point);
+  CeedVectorCreate(ceed, p, &u);
+  CeedVectorCreate(ceed, num_points, &v);
+  CeedVectorCreate(ceed, p, &u_point);
+  CeedVectorCreate(ceed, 1, &v_point);
+  CeedVectorSetValue(v_point, 1.0);
+
+  // Get nodal coordinates
+  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
+  {
+    CeedScalar x_array[2];
+
+    for (CeedInt i = 0; i < 2; i++) x_array[i] = CeedIntPow(-1, i + 1);
+    CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+  }
+  CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
+
+  // Set values of u at nodes
+  {
+    const CeedScalar *x_array;
+    CeedScalar        u_array[p];
+
+    CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
+    for (CeedInt i = 0; i < p; i++) u_array[i] = Eval(x_array[i], ALEN(c), c);
+    CeedVectorRestoreArrayRead(x_nodes, &x_array);
+    CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
+  }
+
+  // Interpolate to arbitrary points
+  CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis_u);
+  {
+    CeedScalar x_array[4] = {-0.33, -0.65, 0.16, 0.99};
+
+    CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+  }
+  CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v);
+
+  for (CeedInt i = 0; i < num_points; i++) {
+    CeedScalar        fx = 0.0;
+    const CeedScalar *x_array, *u_array, *v_array, *u_point_array;
+
+    CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array);
+    CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
+    CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
+    CeedVectorSetValue(x_point, x_array[i]);
+    CeedBasisApplyAtPoints(basis_u, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point);
+    CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array);
+    for (CeedInt j = 0; j < p; j++) fx += u_array[j] * u_point_array[j];
+    if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) printf("%f != %f = f(%f)\n", v_array[i], fx, x_array[i]);
+    CeedVectorRestoreArrayRead(u_point, &u_point_array);
+    CeedVectorRestoreArrayRead(x_points, &x_array);
+    CeedVectorRestoreArrayRead(u, &u_array);
+    CeedVectorRestoreArrayRead(v, &v_array);
+  }
+
+  CeedVectorDestroy(&x);
+  CeedVectorDestroy(&x_nodes);
+  CeedVectorDestroy(&x_points);
+  CeedVectorDestroy(&x_point);
+  CeedVectorDestroy(&u);
+  CeedVectorDestroy(&v);
+  CeedVectorDestroy(&u_point);
+  CeedVectorDestroy(&v_point);
+  CeedBasisDestroy(&basis_x);
+  CeedBasisDestroy(&basis_u);
+  CeedDestroy(&ceed);
+  return 0;
+}
diff --git a/tests/t354-basis.c b/tests/t354-basis.c
@@ -0,0 +1,114 @@
+/// @file
+/// Test polynomial interpolation to arbirtary points in multiple dimensions
+/// \test Test polynomial interpolation to arbitrary points in multiple dimensions
+#include <ceed.h>
+#include <math.h>
+#include <stdio.h>
+
+static CeedScalar Eval(CeedInt dim, const CeedScalar x[]) {
+  CeedScalar result = 1, center = 0.1;
+  for (CeedInt d = 0; d < dim; d++) {
+    result *= tanh(x[d] - center);
+    center += 0.1;
+  }
+  return result;
+}
+
+int main(int argc, char **argv) {
+  Ceed ceed;
+
+  CeedInit(argv[1], &ceed);
+
+  for (CeedInt dim = 1; dim <= 3; dim++) {
+    CeedVector    x, x_nodes, x_points, x_point, u, v, u_point, v_point;
+    CeedBasis     basis_x, basis_u;
+    const CeedInt p = 9, q = 9, num_points = 4, x_dim = CeedIntPow(2, dim), p_dim = CeedIntPow(p, dim);
+
+    CeedVectorCreate(ceed, x_dim * dim, &x);
+    CeedVectorCreate(ceed, p_dim * dim, &x_nodes);
+    CeedVectorCreate(ceed, num_points * dim, &x_points);
+    CeedVectorCreate(ceed, dim, &x_point);
+    CeedVectorCreate(ceed, p_dim, &u);
+    CeedVectorCreate(ceed, num_points, &v);
+    CeedVectorCreate(ceed, p_dim, &u_point);
+    CeedVectorCreate(ceed, 1, &v_point);
+    CeedVectorSetValue(v_point, 1.0);
+
+    // Get nodal coordinates
+    CeedBasisCreateTensorH1Lagrange(ceed, dim, dim, 2, p, CEED_GAUSS_LOBATTO, &basis_x);
+    {
+      CeedScalar x_array[x_dim * dim];
+
+      for (CeedInt d = 0; d < dim; d++) {
+        for (CeedInt i = 0; i < x_dim; i++) x_array[d * x_dim + i] = (i % CeedIntPow(2, dim - d)) / CeedIntPow(2, dim - d - 1) ? 1 : -1;
+      }
+      CeedVectorSetArray(x, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+    }
+    CeedBasisApply(basis_x, 1, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x, x_nodes);
+
+    // Set values of u at nodes
+    {
+      const CeedScalar *x_array;
+      CeedScalar        u_array[p_dim];
+
+      CeedVectorGetArrayRead(x_nodes, CEED_MEM_HOST, &x_array);
+      for (CeedInt i = 0; i < p_dim; i++) {
+        CeedScalar coord[dim];
+
+        for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d * p_dim + i];
+        u_array[i] = Eval(dim, coord);
+      }
+      CeedVectorRestoreArrayRead(x_nodes, &x_array);
+      CeedVectorSetArray(u, CEED_MEM_HOST, CEED_COPY_VALUES, (CeedScalar *)&u_array);
+    }
+
+    // Interpolate to arbitrary points
+    CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, p, q, CEED_GAUSS, &basis_u);
+    {
+      CeedScalar x_array[12] = {-0.33, -0.65, 0.16, 0.99, -0.65, 0.16, 0.99, -0.33, 0.16, 0.99, -0.33, -0.65};
+
+      CeedVectorSetArray(x_points, CEED_MEM_HOST, CEED_COPY_VALUES, x_array);
+    }
+    CeedBasisApplyAtPoints(basis_u, num_points, CEED_NOTRANSPOSE, CEED_EVAL_INTERP, x_points, u, v);
+
+    for (CeedInt i = 0; i < num_points; i++) {
+      CeedScalar        fx = 0.0;
+      CeedScalar        coord[dim];
+      const CeedScalar *x_array, *u_array, *v_array, *u_point_array;
+
+      CeedVectorGetArrayRead(x_points, CEED_MEM_HOST, &x_array);
+      CeedVectorGetArrayRead(u, CEED_MEM_HOST, &u_array);
+      CeedVectorGetArrayRead(v, CEED_MEM_HOST, &v_array);
+      for (CeedInt d = 0; d < dim; d++) coord[d] = x_array[d + i * dim];
+      CeedVectorSetArray(x_point, CEED_MEM_HOST, CEED_COPY_VALUES, coord);
+      CeedBasisApplyAtPoints(basis_u, 1, CEED_TRANSPOSE, CEED_EVAL_INTERP, x_point, v_point, u_point);
+      CeedVectorGetArrayRead(u_point, CEED_MEM_HOST, &u_point_array);
+      for (CeedInt j = 0; j < p_dim; j++) fx += u_array[j] * u_point_array[j];
+      if (fabs(v_array[i] - fx) > 100. * CEED_EPSILON) {
+        // LCOV_EXCL_START
+        printf("[%" CeedInt_FMT "] %f != %f = f(%f", dim, v_array[i], fx, coord[0]);
+        for (CeedInt d = 1; d < dim; d++) printf(", %f", coord[d]);
+        puts(")");
+        // LCOV_EXCL_STOP
+      }
+      CeedVectorRestoreArrayRead(u_point, &u_point_array);
+      CeedVectorRestoreArrayRead(x_points, &x_array);
+      CeedVectorRestoreArrayRead(u, &u_array);
+      CeedVectorRestoreArrayRead(v, &v_array);
+    }
+
+    CeedVectorDestroy(&x);
+    CeedVectorDestroy(&x_nodes);
+    CeedVectorDestroy(&x_points);
+    CeedVectorDestroy(&x_point);
+    CeedVectorDestroy(&u);
+    CeedVectorDestroy(&v);
+    CeedVectorDestroy(&u_point);
+    CeedVectorDestroy(&v_point);
+    CeedBasisDestroy(&basis_x);
+    CeedBasisDestroy(&basis_u);
+  }
+
+  CeedDestroy(&ceed);
+  return 0;
+}