Refactor libCEED operator construction to happen once rather than at …

…every error estimate calculation
awslabs · Mar 26, 2024 · 2eb4ebe · 2eb4ebe
1 parent 21cd51a
commit 2eb4ebe
Show file tree

Hide file tree

Showing 5 changed files with 190 additions and 102 deletions.
diff --git a/palace/fem/libceed/basis.cpp b/palace/fem/libceed/basis.cpp
@@ -120,11 +120,10 @@ void InitMfemInterpolatorBasis(const mfem::FiniteElement &trial_fe,
   MFEM_VERIFY(trial_num_comp == test_num_comp && trial_num_comp == 1,
               "libCEED discrete linear operator requires same vdim = 1 for trial and test "
               "FE spaces!");
-  const int dim = trial_fe.GetDim();
   const int trial_P = trial_fe.GetDof();
   const int test_P = test_fe.GetDof();
-  mfem::DenseMatrix qX(dim, test_P), Gt(trial_P, test_P * dim), Bt;
-  mfem::Vector qW(test_P);
+  mfem::DenseMatrix Bt, Gt(trial_P, test_P);
+  mfem::Vector qX(test_P), qW(test_P);
   mfem::IsoparametricTransformation dummy;
   dummy.SetIdentityTransformation(trial_fe.GetGeomType());
   if (trial_fe.GetMapType() == test_fe.GetMapType())
@@ -159,9 +158,10 @@ void InitMfemInterpolatorBasis(const mfem::FiniteElement &trial_fe,
   qX = 0.0;
   qW = 0.0;
 
-  PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, GetCeedTopology(trial_fe.GetGeomType()),
-                                         trial_num_comp, trial_P, test_P, Bt.GetData(),
-                                         Gt.GetData(), qX.GetData(), qW.GetData(), basis));
+  // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
+  PalaceCeedCall(ceed, CeedBasisCreateH1(ceed, CEED_TOPOLOGY_LINE, trial_num_comp, trial_P,
+                                         test_P, Bt.GetData(), Gt.GetData(), qX.GetData(),
+                                         qW.GetData(), basis));
 }
 
 }  // namespace

diff --git a/palace/fem/libceed/integrator.cpp b/palace/fem/libceed/integrator.cpp
@@ -589,6 +589,7 @@ void AssembleCeedElementErrorIntegrator(
   PalaceCeedCall(ceed, CeedBasisGetNumQuadraturePoints(input1_basis, &num_qpts));
   CeedBasis mesh_elem_basis;
   {
+    // Note: ceed::GetCeedTopology(CEED_TOPOLOGY_LINE) == 1.
     mfem::Vector Bt(num_qpts), Gt(num_qpts), qX(num_qpts), qW(num_qpts);
     Bt = 1.0;
     Gt = 0.0;

diff --git a/palace/fem/libceed/operator.cpp b/palace/fem/libceed/operator.cpp
@@ -67,8 +67,8 @@ void Operator::AddOper(CeedOperator sub_op, CeedOperator sub_op_t)
   PalaceCeedCallBackend(CeedOperatorGetCeed(sub_op, &ceed));
   CeedSize l_in, l_out;
   PalaceCeedCall(ceed, CeedOperatorGetActiveVectorLengths(sub_op, &l_in, &l_out));
-  MFEM_VERIFY(mfem::internal::to_int(l_in) == width &&
-                  mfem::internal::to_int(l_out) == height,
+  MFEM_VERIFY((l_in < 0 || mfem::internal::to_int(l_in) == width) &&
+                  (l_out < 0 || mfem::internal::to_int(l_out) == height),
               "Dimensions mismatch for CeedOperator!");
   PalaceCeedCall(ceed, CeedCompositeOperatorAddSub(op[id], sub_op));
   PalaceCeedCall(ceed, CeedOperatorDestroy(&sub_op));