PetscFE determines basis for examples/fluids/navierstokes (#1264)

* breaking lots of DL stuff to test PetscFE

* fixed missing header

* fluids: Remove unused `Involute` function

* fluids: Set quadrature points for PetscFE

* fluids: Fix GetRestrictionForDomain, Use CreateBasisFromPlex

* fluids: Remove Q_sur arg from SetupStrongBC_Ceed

* fluids: Remove extra SetClosurePermutationTensor

* fluids: Correct quadrature order to auxillary DMs

- Should factor in q_extra

* fluids: Fix parent Restriction creation for spanstats

* fluids: Correct qdata creation for strong bcs

* fluids: Set tensor closer only for tensor-product elements

Co-authored-by: Jeremy L Thompson <jeremy@jeremylt.org>

* fluids: Make shocktube test less brittle

- Small, insignificant changes to the coordinate locations make a
  dramatic difference to the initial condition
- For the test case, there is a node right at the midpoint of the
  domain. Adding the epsilon fudge factor makes the if statement more
  stable

* fluids: Fix SGS test

 - Change the test so that it passes on all hardware. It's a very
   sentitive test and the changes to the basis creation changed the
   initial conditions *just* enough to make the test not happy. Should
   really try and make this more robust at some point...

* fluids: Address misc PR comments

* fluids: Remove quad_mode arg from CreateBasisFromPlex

---------

Co-authored-by: James Wright <james@jameswright.xyz>
Co-authored-by: Jeremy L Thompson <jeremy@jeremylt.org>

Author: 3 people

examples/fluids/navierstokes.c

@@ -18,9 +18,9 @@
 //     ./navierstokes -ceed /cpu/self -problem density_current -degree 1
 //     ./navierstokes -ceed /gpu/cuda -problem advection -degree 1
 //
+//TESTARGS(name="blasius_SGS_DataDriven") -ceed {ceed_resource} -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 2e-9 -state_var primitive -ksp_rtol 1e-12 -snes_rtol 1e-12 -stg_mean_only -stg_fluctuating_IC -test_type solver -compare_final_state_atol 1e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin
 //TESTARGS(name="blasius_aniso_diff_filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 5e-10 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_aniso_vandriest.bin -diff_filter_monitor -ts_max_steps 0 -state_var primitive -diff_filter_friction_length 1e-5 -diff_filter_wall_damping_function van_driest -diff_filter_ksp_rtol 1e-8 -diff_filter_grid_based_width -diff_filter_width_scaling 1,0.7,1
 //TESTARGS(name="blasius_iso_diff_filter") -ceed {ceed_resource} -test_type diff_filter -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius_diff_filter_iso.bin -diff_filter_monitor -ts_max_steps 0 -diff_filter_width_scaling 4.2e-5,4.2e-5,4.2e-5 -diff_filter_ksp_atol 1e-14 -diff_filter_ksp_rtol 1e-16
-//TESTARGS(name="blasius_SGS_DataDriven") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_stgtest.yaml -sgs_model_type data_driven -sgs_model_dd_leakyrelu_alpha 0.3 -sgs_model_dd_parameter_dir examples/fluids/dd_sgs_data -ts_dt 1e-9 -compare_final_state_atol 2e-12 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius-sgs-data-driven.bin -state_var primitive
 //TESTARGS(name="gaussianwave_idl") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/gaussianwave.yaml -compare_final_state_atol 2e-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-gaussianwave-IDL.bin -dm_plex_box_faces 5,5,1 -ts_max_steps 5 -idl_decay_time 2e-3 -idl_length 0.25 -idl_start 0
 //TESTARGS(name="turb_spanstats") -ceed {ceed_resource} -test_type turb_spanstats -options_file examples/fluids/tests-output/stats_test.yaml -compare_final_state_atol 1E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-turb-spanstats-stats.bin
 //TESTARGS(name="blasius") -ceed {ceed_resource} -test_type solver -options_file examples/fluids/tests-output/blasius_test.yaml -compare_final_state_atol 2E-11 -compare_final_state_filename examples/fluids/tests-output/fluids-navierstokes-blasius.bin
@@ -161,7 +161,7 @@ int main(int argc, char **argv) {
   }
 
   // -- Set up DM
-  PetscCall(SetUpDM(dm, problem, app_ctx->degree, bc, phys_ctx));
+  PetscCall(SetUpDM(dm, problem, app_ctx->degree, app_ctx->q_extra, bc, phys_ctx));
 
   // -- Refine DM for high-order viz
   if (app_ctx->viz_refine) PetscCall(VizRefineDM(dm, user, problem, bc, phys_ctx));

examples/fluids/navierstokes.h

@@ -314,9 +314,6 @@ PetscErrorCode PrintRunInfo(User user, Physics phys_ctx, ProblemData *problem, M
 // -----------------------------------------------------------------------------
 // libCEED functions
 // -----------------------------------------------------------------------------
-// Utility function - essential BC dofs are encoded in closure indices as -(i+1).
-PetscInt Involute(PetscInt i);
-
 // Utility function to create local CEED restriction
 PetscErrorCode CreateRestrictionFromPlex(Ceed ceed, DM dm, CeedInt height, DMLabel domain_label, CeedInt label_value, PetscInt dm_field,
                                          CeedElemRestriction *elem_restr);
@@ -326,6 +323,8 @@ PetscErrorCode GetRestrictionForDomain(Ceed ceed, DM dm, CeedInt height, DMLabel
                                        CeedInt q_data_size, CeedElemRestriction *elem_restr_q, CeedElemRestriction *elem_restr_x,
                                        CeedElemRestriction *elem_restr_qd_i);
 
+PetscErrorCode CreateBasisFromPlex(Ceed ceed, DM dm, DMLabel domain_label, CeedInt label_value, CeedInt height, CeedInt dm_field, CeedBasis *basis);
+
 // Utility function to create CEED Composite Operator for the entire domain
 PetscErrorCode CreateOperatorForDomain(Ceed ceed, DM dm, SimpleBC bc, CeedData ceed_data, Physics phys, CeedOperator op_apply_vol,
                                        CeedOperator op_apply_ijacobian_vol, CeedInt height, CeedInt P_sur, CeedInt Q_sur, CeedInt q_data_size_sur,
@@ -361,7 +360,7 @@ PetscErrorCode UpdateBoundaryValues(User user, Vec Q_loc, PetscReal t);
 PetscErrorCode CreateDM(MPI_Comm comm, ProblemData *problem, MatType, VecType, DM *dm);
 
 // Set up DM
-PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc, Physics phys);
+PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, PetscInt q_extra, SimpleBC bc, Physics phys);
 
 // Refine DM for high-order viz
 PetscErrorCode VizRefineDM(DM dm, User user, ProblemData *problem, SimpleBC bc, Physics phys);
@@ -445,8 +444,7 @@ PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User use
 // -----------------------------------------------------------------------------
 
 // Setup StrongBCs that use QFunctions
-PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt Q_sur,
-                                  CeedInt q_data_size_sur);
+PetscErrorCode SetupStrongBC_Ceed(Ceed ceed, CeedData ceed_data, DM dm, User user, ProblemData *problem, SimpleBC bc, CeedInt q_data_size_sur);
 
 PetscErrorCode FreestreamBCSetup(ProblemData *problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference);
 PetscErrorCode OutflowBCSetup(ProblemData *problem, DM dm, void *ctx, NewtonianIdealGasContext newtonian_ig_ctx, const StatePrimitive *reference);

examples/fluids/problems/sgs_dd_model.c

@@ -29,19 +29,20 @@ PetscErrorCode SGS_DD_ModelSetupDataDestroy(SGS_DD_ModelSetupData sgs_dd_setup_d
 }
 
 // @brief Create DM for storing subgrid stress at nodes
-PetscErrorCode SGS_DD_ModelCreateDM(DM dm_source, DM *dm_sgs, PetscInt degree, PetscInt *num_components) {
+PetscErrorCode SGS_DD_ModelCreateDM(DM dm_source, DM *dm_sgs, PetscInt degree, PetscInt q_extra, PetscInt *num_components) {
   PetscFE      fe;
   PetscSection section;
   PetscInt     dim;
 
   PetscFunctionBeginUser;
-  *num_components = 6;
+  *num_components  = 6;
+  PetscInt q_order = degree + q_extra;
 
   PetscCall(DMClone(dm_source, dm_sgs));
   PetscCall(DMGetDimension(*dm_sgs, &dim));
   PetscCall(PetscObjectSetName((PetscObject)*dm_sgs, "Subgrid Stress Projection"));
 
-  PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, *num_components, PETSC_FALSE, degree, PETSC_DECIDE, &fe));
+  PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, *num_components, PETSC_FALSE, degree, q_order, &fe));
   PetscCall(PetscObjectSetName((PetscObject)fe, "Subgrid Stress Projection"));
   PetscCall(DMAddField(*dm_sgs, NULL, (PetscObject)fe));
   PetscCall(DMCreateDS(*dm_sgs));
@@ -66,7 +67,7 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c
   SGS_DD_Data         sgs_dd_data = user->sgs_dd_data;
   CeedQFunction       qf_multiplicity, qf_sgs_dd_nodal;
   CeedOperator        op_multiplicity, op_sgs_dd_nodal;
-  CeedInt             num_elem, elem_size, num_comp_q, num_qpts_1d, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso;
+  CeedInt             num_elem, elem_size, num_comp_q, num_comp_grad_velo, num_comp_x, num_comp_grid_aniso;
   PetscInt            dim;
   CeedVector          multiplicity, inv_multiplicity;
   CeedElemRestriction elem_restr_inv_multiplicity, elem_restr_grad_velo, elem_restr_sgs;
@@ -78,16 +79,14 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c
   CeedElemRestrictionGetNumComponents(sgs_dd_setup_data->elem_restr_grid_aniso, &num_comp_grid_aniso);
   CeedElemRestrictionGetNumElements(ceed_data->elem_restr_q, &num_elem);
   CeedElemRestrictionGetElementSize(ceed_data->elem_restr_q, &elem_size);
-  CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d);
 
   {  // Get velocity gradient information
     CeedOperatorField op_field;
     CeedOperatorGetFieldByName(user->grad_velo_proj->l2_rhs_ctx->op, "velocity gradient", &op_field);
     CeedOperatorFieldGetElemRestriction(op_field, &elem_restr_grad_velo);
     CeedElemRestrictionGetNumComponents(elem_restr_grad_velo, &num_comp_grad_velo);
   }
-
-  PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, num_qpts_1d, 0, &elem_restr_sgs, NULL, NULL));
+  PetscCall(GetRestrictionForDomain(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, -1, 0, &elem_restr_sgs, NULL, NULL));
   CeedElemRestrictionCreateVector(elem_restr_sgs, &sgs_dd_data->sgs_nodal_ceed, NULL);
 
   // -- Create inverse multiplicity for correcting nodal assembly
@@ -160,7 +159,7 @@ PetscErrorCode SGS_DD_ModelSetupNodalEvaluation(Ceed ceed, User user, CeedData c
 // @brief Create CeedOperator to compute SGS contribution to the residual
 PetscErrorCode SGS_ModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_data, SGS_DD_ModelSetupData sgs_dd_setup_data) {
   SGS_DD_Data   sgs_dd_data = user->sgs_dd_data;
-  CeedInt       num_comp_q, num_comp_qd, num_comp_x, num_qpts_1d, num_nodes_1d;
+  CeedInt       num_comp_q, num_comp_qd, num_comp_x;
   PetscInt      dim;
   CeedQFunction qf_sgs_apply;
   CeedOperator  op_sgs_apply;
@@ -171,10 +170,8 @@ PetscErrorCode SGS_ModelSetupNodalIFunction(Ceed ceed, User user, CeedData ceed_
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x);
-  CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d);
-  CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d);
 
-  CeedBasisCreateTensorH1Lagrange(ceed, dim, sgs_dd_data->num_comp_sgs, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_sgs);
+  PetscCall(CreateBasisFromPlex(ceed, sgs_dd_data->dm_sgs, 0, 0, 0, 0, &basis_sgs));
 
   switch (user->phys->state_var) {
     case STATEVAR_PRIMITIVE:
@@ -330,7 +327,8 @@ PetscErrorCode SGS_DD_ModelSetup(Ceed ceed, User user, CeedData ceed_data, Probl
 
   // -- Create DM for storing SGS tensor at nodes
   PetscCall(PetscNew(&user->sgs_dd_data));
-  PetscCall(SGS_DD_ModelCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, &user->sgs_dd_data->num_comp_sgs));
+  PetscCall(
+      SGS_DD_ModelCreateDM(user->dm, &user->sgs_dd_data->dm_sgs, user->app_ctx->degree, user->app_ctx->q_extra, &user->sgs_dd_data->num_comp_sgs));
 
   PetscCall(PetscNew(&sgs_dd_setup_data));
 

examples/fluids/qfunctions/shocktube.h

@@ -94,7 +94,7 @@ CEED_QFUNCTION_HELPER CeedInt Exact_ShockTube(CeedInt dim, CeedScalar time, cons
   CeedScalar rho, P, u[3] = {0.};
 
   // Initial Conditions
-  if (x <= mid_point) {
+  if (x <= mid_point + 200 * CEED_EPSILON) {
     rho = rho_high;
     P   = P_high;
   } else {

examples/fluids/src/differential_filter.c

@@ -17,16 +17,14 @@
 PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData ceed_data, CeedQFunctionContext diff_filter_qfctx) {
   DiffFilterData diff_filter = user->diff_filter;
   DM             dm_filter   = diff_filter->dm_filter;
-  CeedInt        num_comp_q, num_comp_qd, num_qpts_1d, num_nodes_1d, num_comp_x;
+  CeedInt        num_comp_q, num_comp_qd, num_comp_x;
   PetscInt       dim;
 
   PetscFunctionBeginUser;
   PetscCall(DMGetDimension(user->dm, &dim));
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_x, &num_comp_x);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_q, &num_comp_q);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &num_comp_qd);
-  CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d);
-  CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d);
 
   {  // -- Create RHS MatopApplyContext
     CeedQFunction qf_rhs;
@@ -64,9 +62,8 @@ PetscErrorCode DifferentialFilterCreateOperators(Ceed ceed, User user, CeedData
       char                field_name[PETSC_MAX_PATH_LEN];
       CeedElemRestriction elem_restr_filter;
       CeedBasis           basis_filter;
-
-      PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, num_qpts_1d, 0, &elem_restr_filter, NULL, NULL));
-      CeedBasisCreateTensorH1Lagrange(ceed, dim, diff_filter->num_field_components[i], num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_filter);
+      PetscCall(GetRestrictionForDomain(ceed, dm_filter, 0, 0, 0, i, -1, 0, &elem_restr_filter, NULL, NULL));
+      PetscCall(CreateBasisFromPlex(ceed, dm_filter, 0, 0, 0, i, &basis_filter));
 
       PetscCall(PetscSNPrintf(field_name, PETSC_MAX_PATH_LEN, "v%" PetscInt_FMT, i));
       CeedOperatorSetField(op_rhs, field_name, elem_restr_filter, basis_filter, CEED_VECTOR_ACTIVE);
@@ -187,7 +184,7 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data,
   {  // Create DM for filtered quantities
     PetscFE      fe;
     PetscSection section;
-    PetscInt     dim;
+    PetscInt     dim, q_order = user->app_ctx->degree + user->app_ctx->q_extra;
 
     PetscCall(DMClone(user->dm, &diff_filter->dm_filter));
     PetscCall(DMGetDimension(diff_filter->dm_filter, &dim));
@@ -198,8 +195,7 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data,
 
     if (diff_filter->do_mms_test) {
       diff_filter->num_field_components[0] = 1;
-      PetscCall(
-          PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe));
+      PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, q_order, &fe));
       PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - MMS"));
       PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe));
       PetscCall(PetscFEDestroy(&fe));
@@ -209,15 +205,13 @@ PetscErrorCode DifferentialFilterSetup(Ceed ceed, User user, CeedData ceed_data,
       PetscCall(PetscSectionSetComponentName(section, 0, 0, "FilteredPhi"));
     } else {
       diff_filter->num_field_components[0] = DIFF_FILTER_STATE_NUM;
-      PetscCall(
-          PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe));
+      PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[0], PETSC_FALSE, user->app_ctx->degree, q_order, &fe));
       PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - Primitive State Variables"));
       PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe));
       PetscCall(PetscFEDestroy(&fe));
 
       diff_filter->num_field_components[1] = DIFF_FILTER_VELOCITY_SQUARED_NUM;
-      PetscCall(
-          PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[1], PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe));
+      PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, diff_filter->num_field_components[1], PETSC_FALSE, user->app_ctx->degree, q_order, &fe));
       PetscCall(PetscObjectSetName((PetscObject)fe, "Differential Filtering - Velocity Products"));
       PetscCall(DMAddField(diff_filter->dm_filter, NULL, (PetscObject)fe));
       PetscCall(PetscFEDestroy(&fe));

examples/fluids/src/grid_anisotropy_tensor.c

@@ -19,7 +19,7 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
   CeedQFunction        qf_rhs_assemble, qf_mass;
   CeedBasis            basis_grid_aniso;
   PetscInt             dim;
-  CeedInt              num_qpts_1d, num_nodes_1d, q_data_size;
+  CeedInt              q_data_size;
   MPI_Comm             comm = PetscObjectComm((PetscObject)user->dm);
   KSP                  ksp;
 
@@ -35,7 +35,8 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
   {  // -- Setup DM
     PetscFE      fe;
     PetscSection section;
-    PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grid_aniso_proj->num_comp, PETSC_FALSE, user->app_ctx->degree, PETSC_DECIDE, &fe));
+    PetscInt     q_order = user->app_ctx->degree + user->app_ctx->q_extra;
+    PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, grid_aniso_proj->num_comp, PETSC_FALSE, user->app_ctx->degree, q_order, &fe));
     PetscCall(PetscObjectSetName((PetscObject)fe, "Grid Anisotropy Tensor Projection"));
     PetscCall(DMAddField(grid_aniso_proj->dm, NULL, (PetscObject)fe));
     PetscCall(DMCreateDS(grid_aniso_proj->dm));
@@ -55,13 +56,10 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
   }
 
   // -- Get Pre-requisite things
-  CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d);
-  CeedBasisGetNumNodes1D(ceed_data->basis_q, &num_nodes_1d);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size);
 
-  PetscCall(
-      GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, num_qpts_1d, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL));
-  CeedBasisCreateTensorH1Lagrange(ceed, dim, grid_aniso_proj->num_comp, num_nodes_1d, num_qpts_1d, CEED_GAUSS, &basis_grid_aniso);
+  PetscCall(GetRestrictionForDomain(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, -1, grid_aniso_proj->num_comp, elem_restr_grid_aniso, NULL, NULL));
+  PetscCall(CreateBasisFromPlex(ceed, grid_aniso_proj->dm, 0, 0, 0, 0, &basis_grid_aniso));
 
   // -- Build RHS operator
   CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorProjection, AnisotropyTensorProjection_loc, &qf_rhs_assemble);
@@ -136,17 +134,15 @@ PetscErrorCode GridAnisotropyTensorProjectionSetupApply(Ceed ceed, User user, Ce
 
 PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User user, CeedData ceed_data, CeedElemRestriction *elem_restr_grid_aniso,
                                                              CeedVector *aniso_colloc_ceed, PetscInt *num_comp_aniso) {
-  CeedInt       q_data_size, num_qpts_1d, num_qpts;
+  CeedInt       q_data_size, num_nodes;
   CeedQFunction qf_colloc;
   CeedOperator  op_colloc;
 
   PetscFunctionBeginUser;
-  // -- Get Pre-requisite things
   *num_comp_aniso = 7;
-  CeedBasisGetNumQuadraturePoints1D(ceed_data->basis_q, &num_qpts_1d);
+  CeedBasisGetNumNodes(ceed_data->basis_q, &num_nodes);
   CeedElemRestrictionGetNumComponents(ceed_data->elem_restr_qd_i, &q_data_size);
-
-  PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_qpts_1d, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso));
+  PetscCall(GetRestrictionForDomain(ceed, user->dm, 0, 0, 0, 0, num_nodes, *num_comp_aniso, NULL, NULL, elem_restr_grid_aniso));
 
   // -- Build collocation operator
   CeedQFunctionCreateInterior(ceed, 1, AnisotropyTensorCollocate, AnisotropyTensorCollocate_loc, &qf_colloc);
@@ -156,7 +152,6 @@ PetscErrorCode GridAnisotropyTensorCalculateCollocatedVector(Ceed ceed, User use
   CeedOperatorCreate(ceed, qf_colloc, NULL, NULL, &op_colloc);
   CeedOperatorSetField(op_colloc, "qdata", ceed_data->elem_restr_qd_i, CEED_BASIS_COLLOCATED, ceed_data->q_data);
   CeedOperatorSetField(op_colloc, "v", *elem_restr_grid_aniso, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
-  CeedBasisGetNumQuadraturePoints(ceed_data->basis_q, &num_qpts);
 
   CeedElemRestrictionCreateVector(*elem_restr_grid_aniso, aniso_colloc_ceed, NULL);