Start adding pressure bc

Author: rezgarshakeri

examples/Hdiv-mixed/main.c

@@ -158,9 +158,34 @@ int main(int argc, char **argv) {
   // -- Set up libCEED objects
   PetscCall( SetupLibceed(dm, dm_u0, dm_p0, dm_H1, ceed, app_ctx,
                           problem_data, ceed_data) );
-  //CeedVectorView(force_ceed, "%12.8f", stdout);
-  //PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
-  //                                   bc_pressure) );
+
+  // ---------------------------------------------------------------------------
+  // Setup pressure boundary conditions
+  // ---------------------------------------------------------------------------
+  // --Create empty local vector for libCEED
+  Vec P_loc;
+  PetscInt P_loc_size;
+  CeedScalar *p0;
+  CeedVector P_ceed;
+  PetscMemType pressure_mem_type;
+  PetscCall( DMCreateLocalVector(dm, &P_loc) );
+  PetscCall( VecGetSize(P_loc, &P_loc_size) );
+  PetscCall( VecZeroEntries(P_loc) );
+  PetscCall( VecGetArrayAndMemType(P_loc, &p0, &pressure_mem_type) );
+  CeedVectorCreate(ceed, P_loc_size, &P_ceed);
+  CeedVectorSetArray(P_ceed, MemTypeP2C(pressure_mem_type), CEED_USE_POINTER,
+                     p0);
+  // -- Apply operator to create local pressure vector on boundary
+  PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
+                                     P_ceed) );
+  CeedVectorView(P_ceed, "%12.8f", stdout);
+  // -- Map local to global
+  Vec P;
+  CeedVectorTakeArray(P_ceed, MemTypeP2C(pressure_mem_type), NULL);
+  PetscCall( VecRestoreArrayAndMemType(P_loc, &p0) );
+  PetscCall( DMCreateGlobalVector(dm, &P) );
+  PetscCall( VecZeroEntries(P) );
+  PetscCall( DMLocalToGlobal(dm, P_loc, ADD_VALUES, P) );
 
   // ---------------------------------------------------------------------------
   // Setup context for projection problem; post-processing.c

examples/Hdiv-mixed/src/setup-libceed.c

@@ -285,20 +285,20 @@ PetscErrorCode SetupLibceed(DM dm, DM dm_u0, DM dm_p0, DM dm_H1,
     L2BasisP0(dim, Q, q_ref, q_weights, interp_p, problem_data->quadrature_mode);
     CeedBasisCreateH1(ceed, CEED_TOPOLOGY_QUAD, num_comp_p, 1, num_qpts, interp_p,
                       grad, q_ref,q_weights, &ceed_data->basis_p);
-    //HdivBasisQuad(Q, q_ref, q_weights, interp_u, div,
-    //              CEED_GAUSS_LOBATTO);
-    //CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp_u, P_u, num_qpts,
-    //                    interp_u, div, q_ref, q_weights, &ceed_data->basis_u_face);
+    HdivBasisQuad(Q, q_ref, q_weights, interp_u, div,
+                  CEED_GAUSS_LOBATTO);
+    CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_QUAD, num_comp_u, P_u, num_qpts,
+                        interp_u, div, q_ref, q_weights, &ceed_data->basis_u_face);
   } else {
     HdivBasisHex(Q, q_ref, q_weights, interp_u, div, problem_data->quadrature_mode);
     CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_HEX, num_comp_u, P_u, num_qpts,
                         interp_u, div, q_ref, q_weights, &ceed_data->basis_u);
     L2BasisP0(dim, Q, q_ref, q_weights, interp_p, problem_data->quadrature_mode);
     CeedBasisCreateH1(ceed, CEED_TOPOLOGY_HEX, num_comp_p, 1, num_qpts, interp_p,
                       grad, q_ref,q_weights, &ceed_data->basis_p);
-    //HdivBasisHex(Q, q_ref, q_weights, interp_u, div, CEED_GAUSS_LOBATTO);
-    //CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_HEX, num_comp_u, P_u, num_qpts,
-    //                    interp_u, div, q_ref, q_weights, &ceed_data->basis_u_face);
+    HdivBasisHex(Q, q_ref, q_weights, interp_u, div, CEED_GAUSS_LOBATTO);
+    CeedBasisCreateHdiv(ceed, CEED_TOPOLOGY_HEX, num_comp_u, P_u, num_qpts,
+                        interp_u, div, q_ref, q_weights, &ceed_data->basis_u_face);
   }
 
   CeedBasisCreateTensorH1Lagrange(ceed, dim, num_comp_x, 2, Q,