Merge pull request #1632 from CEED/jrwrigh/fix_basis_project

fix(basis): Use basis_from dim for projection matrices
CEED · Jul 8, 2024 · 1088199 · 1088199
2 parents 0de0756 + 706bf52
commit 1088199
Show file tree

Hide file tree

Showing 3 changed files with 109 additions and 1 deletion.
diff --git a/interface/ceed-basis.c b/interface/ceed-basis.c
@@ -234,7 +234,7 @@ static int CeedBasisCreateProjectionMatrices(CeedBasis basis_from, CeedBasis bas
   CeedScalar       *interp_to_inv, *interp_from;
   const CeedScalar *interp_to_source = NULL, *interp_from_source = NULL, *grad_from_source = NULL;
 
-  CeedCall(CeedBasisGetDimension(basis_to, &dim));
+  CeedCall(CeedBasisGetDimension(basis_from, &dim));
   if (are_both_tensor) {
     CeedCall(CeedBasisGetInterp1D(basis_to, &interp_to_source));
     CeedCall(CeedBasisGetInterp1D(basis_from, &interp_from_source));

diff --git a/tests/t319-basis.c b/tests/t319-basis.c
@@ -1,6 +1,7 @@
 /// @file
 /// Test projection interp and grad in multiple dimensions
 /// \test Test projection interp and grad in multiple dimensions
+#include "t319-basis.h"
 #include <ceed.h>
 #include <math.h>
 #include <stdio.h>
@@ -201,6 +202,41 @@ int main(int argc, char **argv) {
     CeedBasisDestroy(&basis_to_nontensor);
     CeedBasisDestroy(&basis_project);
   }
+
+  // Test projection between basis of different topological dimension
+  {
+    CeedInt   face_dim = 2, P_1D = 2;
+    CeedBasis basis_face, basis_cell_to_face, basis_proj;
+
+    CeedScalar       *q_ref = NULL, *q_weights = NULL;
+    const CeedScalar *grad, *interp;
+    CeedInt           P, Q;
+    GetCellToFaceTabulation(CEED_GAUSS, &P, &Q, &interp, &grad);
+
+    CeedBasisCreateTensorH1Lagrange(ceed, face_dim, 1, 2, P_1D, CEED_GAUSS, &basis_face);
+    CeedBasisCreateH1(ceed, CEED_TOPOLOGY_HEX, 1, P, Q, (CeedScalar *)interp, (CeedScalar *)grad, q_ref, q_weights, &basis_cell_to_face);
+    CeedBasisCreateProjection(basis_cell_to_face, basis_face, &basis_proj);
+    const CeedScalar *interp_proj, *grad_proj, *interp_proj_ref, *grad_proj_ref;
+
+    GetCellToFaceTabulation(CEED_GAUSS_LOBATTO, NULL, NULL, &interp_proj_ref, &grad_proj_ref);
+    CeedBasisGetInterp(basis_proj, &interp_proj);
+    CeedBasisGetGrad(basis_proj, &grad_proj);
+    CeedScalar tol = 100 * CEED_EPSILON;
+
+    for (CeedInt i = 0; i < 4 * 8; i++) {
+      if (fabs(interp_proj[i] - ((CeedScalar *)interp_proj_ref)[i]) > tol)
+        printf("Mixed Topology Projection: interp[%" CeedInt_FMT "] expected %f, got %f\n", i, interp_proj[i], ((CeedScalar *)interp_proj_ref)[i]);
+    }
+
+    for (CeedInt i = 0; i < 3 * 4 * 8; i++) {
+      if (fabs(grad_proj[i] - ((CeedScalar *)grad_proj_ref)[i]) > tol)
+        printf("Mixed Topology Projection: grad[%" CeedInt_FMT "] expected %f, got %f\n", i, grad_proj[i], ((CeedScalar *)grad_proj_ref)[i]);
+    }
+
+    CeedBasisDestroy(&basis_face);
+    CeedBasisDestroy(&basis_cell_to_face);
+    CeedBasisDestroy(&basis_proj);
+  }
   CeedDestroy(&ceed);
   return 0;
 }
diff --git a/tests/t319-basis.h b/tests/t319-basis.h
@@ -0,0 +1,72 @@
+// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
+// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
+//
+// SPDX-License-Identifier: BSD-2-Clause
+//
+// This file is part of CEED:  http://github.com/ceed
+
+#include <ceed.h>
+
+// Interpolation matrices for cell-to-face of Q1 hexahedral element onto it's "5" face (in PETSc)
+// Nodes are at Gauss-Lobatto points and quadrature points are Gauss, all over [-1,1] domain range
+const CeedScalar Q1_interp_gauss[4][8] = {
+    {0.62200846792814612,  0, 0.16666666666666669,  0, 0.16666666666666669,  0, 0.044658198738520463, 0},
+    {0.16666666666666669,  0, 0.62200846792814612,  0, 0.044658198738520463, 0, 0.16666666666666669,  0},
+    {0.16666666666666669,  0, 0.044658198738520463, 0, 0.62200846792814612,  0, 0.16666666666666669,  0},
+    {0.044658198738520463, 0, 0.16666666666666669,  0, 0.16666666666666669,  0, 0.62200846792814612,  0}
+};
+const CeedScalar Q1_grad_gauss[3][4][8] = {
+    {{-0.31100423396407312, 0.31100423396407312, -0.083333333333333343, 0.083333333333333343, -0.083333333333333343, 0.083333333333333343,
+      -0.022329099369260232, 0.022329099369260232},
+     {-0.083333333333333343, 0.083333333333333343, -0.31100423396407312, 0.31100423396407312, -0.022329099369260232, 0.022329099369260232,
+      -0.083333333333333343, 0.083333333333333343},
+     {-0.083333333333333343, 0.083333333333333343, -0.022329099369260232, 0.022329099369260232, -0.31100423396407312, 0.31100423396407312,
+      -0.083333333333333343, 0.083333333333333343},
+     {-0.022329099369260232, 0.022329099369260232, -0.083333333333333343, 0.083333333333333343, -0.083333333333333343, 0.083333333333333343,
+      -0.31100423396407312, 0.31100423396407312}                                                       },
+    {{-0.39433756729740643, 0, 0.39433756729740643, 0, -0.10566243270259357, 0, 0.10566243270259357, 0},
+     {-0.39433756729740643, 0, 0.39433756729740643, 0, -0.10566243270259357, 0, 0.10566243270259357, 0},
+     {-0.10566243270259357, 0, 0.10566243270259357, 0, -0.39433756729740643, 0, 0.39433756729740643, 0},
+     {-0.10566243270259357, 0, 0.10566243270259357, 0, -0.39433756729740643, 0, 0.39433756729740643, 0}},
+    {{-0.39433756729740643, 0, -0.10566243270259357, 0, 0.39433756729740643, 0, 0.10566243270259357, 0},
+     {-0.10566243270259357, 0, -0.39433756729740643, 0, 0.10566243270259357, 0, 0.39433756729740643, 0},
+     {-0.39433756729740643, 0, -0.10566243270259357, 0, 0.39433756729740643, 0, 0.10566243270259357, 0},
+     {-0.10566243270259357, 0, -0.39433756729740643, 0, 0.10566243270259357, 0, 0.39433756729740643, 0}}
+};
+
+const CeedScalar Q1_interp_gauss_lobatto[4][8] = {
+    {1, 0, 0, 0, 0, 0, 0, 0},
+    {0, 0, 1, 0, 0, 0, 0, 0},
+    {0, 0, 0, 0, 1, 0, 0, 0},
+    {0, 0, 0, 0, 0, 0, 1, 0}
+};
+/* clang-format off */
+const CeedScalar Q1_grad_gauss_lobatto[3][4][8] = {
+    {{-0.5,  0.5, 0,    0,   0,    0,   0,    0},
+      {0,    0,   -0.5, 0.5, 0,    0,   0,    0},
+      {0,    0,   0,    0,   -0.5, 0.5, 0,    0},
+      {0,    0,   0,    0,   0,    0,   -0.5, 0.5}},
+    {{-0.5,  0,   0.5,  0,   0,    0,   0,    0},
+      {-0.5, 0,   0.5,  0,   0,    0,   0,    0},
+      {0,    0,   0,    0,   -0.5, 0,   0.5,  0},
+      {0,    0,   0,    0,   -0.5, 0,   0.5,  0}},
+    {{-0.5,  0,   0,    0,   0.5,  0,   0,    0},
+      {0,    0,   -0.5, 0,   0,    0,   0.5,  0},
+      {-0.5, 0,   0,    0,   0.5,  0,   0,    0},
+      {0,    0,   -0.5, 0,   0,    0,   0.5,  0}}
+};
+/* clang-format on */
+
+static void GetCellToFaceTabulation(CeedQuadMode quad_mode, CeedInt *P, CeedInt *Q, const CeedScalar **interp, const CeedScalar **grad) {
+  if (P) *P = 8;
+  if (Q) *Q = 4;
+
+  if (quad_mode == CEED_GAUSS) {
+    *interp = (const CeedScalar *)Q1_interp_gauss;
+    *grad   = (const CeedScalar *)Q1_grad_gauss;
+  }
+  if (quad_mode == CEED_GAUSS_LOBATTO) {
+    *interp = (const CeedScalar *)Q1_interp_gauss_lobatto;
+    *grad   = (const CeedScalar *)Q1_grad_gauss_lobatto;
+  }
+}