wip - framing out new approach

CEED · Apr 7, 2023 · c1dae63 · c1dae63
1 parent 8a8689a
commit c1dae63
Show file tree

Hide file tree

Showing 2 changed files with 37 additions and 42 deletions.
diff --git a/include/ceed-impl.h b/include/ceed-impl.h
@@ -187,9 +187,10 @@ struct CeedBasis_private {
   CeedScalar      *interp;       /* row-major matrix of shape [Q_comp*Q, P] expressing the values of nodal basis functions at quadrature points */
   CeedScalar      *interp_1d;    /* row-major matrix of shape [Q1d, P1d] expressing the values of nodal basis functions at quadrature points */
   CeedScalar      *grad; /* row-major matrix of shape [dim*Q_comp*Q, P] matrix expressing derivatives of nodal basis functions at quadrature points */
-  CeedScalar      *grad_1d;    /* row-major matrix of shape [Q1d, P1d] matrix expressing derivatives of nodal basis functions at quadrature points */
-  CeedTensorContract contract; /* tensor contraction object */
-  CeedFESpace        basis_space; /* Initialize in basis constructor with 1,2 for H^1, H(div) FE space */
+  CeedScalar      *grad_1d; /* row-major matrix of shape [Q1d, P1d] matrix expressing derivatives of nodal basis functions at quadrature points */
+  CeedScalar      *legendre_coeffs_1d; /* row-major matrix mapping from quadrature point values to Legendre polynomial coefficients */
+  CeedTensorContract contract;         /* tensor contraction object */
+  CeedFESpace        basis_space;      /* Initialize in basis constructor with 1,2 for H^1, H(div) FE space */
   CeedScalar
        *div; /* row-major matrix of shape [Q, P] expressing the divergence of nodal basis functions at quadrature points for H(div) discretizations */
   void *data; /* place for the backend to store any data */

diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
@@ -1320,9 +1320,7 @@ int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode,
   CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
   CeedCall(CeedBasisGetNumQuadraturePoints(basis, &num_qpts));
   CeedCall(CeedVectorGetLength(v, &v_length));
-  if (u) {
-    CeedCall(CeedVectorGetLength(u, &u_length));
-  }
+  if (u) CeedCall(CeedVectorGetLength(u, &u_length));
 
   if (!basis->Apply) {
     // LCOV_EXCL_START
@@ -1400,9 +1398,7 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
   CeedCall(CeedBasisGetNumNodes(basis, &num_nodes));
   CeedCall(CeedVectorGetLength(x_ref, &x_length));
   CeedCall(CeedVectorGetLength(v, &v_length));
-  if (u) {
-    CeedCall(CeedVectorGetLength(u, &u_length));
-  }
+  if (u) CeedCall(CeedVectorGetLength(u, &u_length));
 
   // Check compatibility of topological and geometrical dimensions
   if ((t_mode == CEED_TRANSPOSE && (v_length % num_nodes != 0 || u_length != (num_points * num_comp))) ||
@@ -1445,39 +1441,36 @@ int CeedBasisApplyAtPoints(CeedBasis basis, CeedInt num_points, CeedTransposeMod
     // LCOV_EXCL_STOP
   }
 
-  if (basis->Apply) CeedCall(basis->ApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v));
-
-  // Create and evaluate basis
-  {
-    CeedScalar c1, c2, c3, c4, dx;
-    CeedScalar interp_1d[num_nodes * num_points], grad_1d[num_nodes * num_points * dim];
-
-    // TODO: need basis node coordinates, not currently provided
-    // TODO: need to de-tensor snippit below
-
-    // Build B, D matrix
-    // Fornberg, 1998
-    for (CeedInt i = 0; i < num_points; i++) {
-      c1                   = 1.0;
-      c3                   = nodes[0] - x_ref[i + d * num_points];
-      interp_1d[i * num_nodes + 0] = 1.0;
-      for (j = 1; j < num_nodes; j++) {
-        c2 = 1.0;
-        c4 = c3;
-        c3 = nodes[j] - x_ref[i + d * num_points];
-        for (CeedInt k = 0; k < j; k++) {
-          dx = nodes[j] - nodes[k];
-          c2 *= dx;
-          if (k == j - 1) {
-            grad_1d[i * num_nodes + j]   = c1 * (interp_1d[i * num_nodes + k] - c4 * grad_1d[i * num_nodes + k]) / c2;
-            interp_1d[i * num_nodes + j] = -c1 * c4 * interp_1d[i * num_nodes + k] / c2;
-          }
-          grad_1d[i * num_nodes + k]   = (c3 * grad_1d[i * num_nodes + k] - interp_1d[i * num_nodes + k]) / dx;
-          interp_1d[i * num_nodes + k] = c3 * interp_1d[i * num_nodes + k] / dx;
-        }
-        c1 = c2;
-      }
-    }
+  // Backend method
+  if (basis->ApplyAtPoints) {
+    CeedCall(basis->ApplyAtPoints(basis, num_points, t_mode, eval_mode, x_ref, u, v));
+    return CEED_ERROR_SUCCESS;
+  }
+
+  // Default implementation
+  if (!basis->legendre_coeffs_1d) {
+    // Build matrix mapping from quadrature point values to Legendre coefficients
+  }
+  if (t_mode == CEED_NOTRANSPOSE) {
+    // Nodes to arbitrary points
+    CeedVector q;
+
+    // -- Interpolate to quadrature points
+    CeedCall(CeedVectorCreate(basis->ceed, basis->Q * basis->Q_comp, &q));
+    CeedCall(CeedBasisApply(basis, 1, t_mode, CEED_EVAL_INTERP, u, q));
+
+    // -- Map to Legendre coefficients
+
+    // -- Evaluate Legendre polynomial at arbitrary points
+  } else {
+    // Arbitrary points to nodes
+    CeedVector q;
+
+    // -- Arbitrary points to quadrature points
+    CeedCall(CeedVectorCreate(basis->ceed, basis->Q * basis->Q_comp, &q));
+
+    // -- Transpose interpolation from quadrature points
+    CeedCall(CeedBasisApply(basis, 1, t_mode, CEED_EVAL_INTERP, q, u));
   }
 
   return CEED_ERROR_SUCCESS;
@@ -1813,6 +1806,7 @@ int CeedBasisDestroy(CeedBasis *basis) {
   CeedCall(CeedFree(&(*basis)->grad));
   CeedCall(CeedFree(&(*basis)->div));
   CeedCall(CeedFree(&(*basis)->grad_1d));
+  CeedCall(CeedFree(&(*basis)->legendre_coeffs_1d));
   CeedCall(CeedFree(&(*basis)->q_ref_1d));
   CeedCall(CeedFree(&(*basis)->q_weight_1d));
   CeedCall(CeedDestroy(&(*basis)->ceed));