Pressure boundary condition added. Need to rotate normals.

Author: rezgarshakeri

examples/Hdiv-mixed/basis/Hdiv-hex.h

@@ -1,3 +1,5 @@
+#ifndef Hdiv_hex_h
+#define Hdiv_hex_h
 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
 // All Rights reserved. See files LICENSE and NOTICE for details.
@@ -159,4 +161,4 @@ static void HdivBasisHex(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights,
   }
 }
 
-
+#endif // Hdiv_hex_h

examples/Hdiv-mixed/basis/Hdiv-quad.h

@@ -1,3 +1,5 @@
+#ifndef Hdiv_quad_h
+#define Hdiv_quad_h
 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC.
 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707.
 // All Rights reserved. See files LICENSE and NOTICE for details.
@@ -84,4 +86,4 @@ static void HdivBasisQuad(CeedInt Q, CeedScalar *q_ref, CeedScalar *q_weights,
   }
 }
 
-
+#endif // Hdiv_quad_h

examples/Hdiv-mixed/include/cl-options.h

@@ -1,9 +1,9 @@
 #ifndef cloptions_h
 #define cloptions_h
 
-#include "../include/structs.h"
+#include "structs.h"
 
 // Process general command line options
-PetscErrorCode ProcessCommandLineOptions(MPI_Comm comm, AppCtx app_ctx);
+PetscErrorCode ProcessCommandLineOptions(AppCtx app_ctx);
 
 #endif // cloptions_h

examples/Hdiv-mixed/include/register-problem.h

@@ -1,7 +1,7 @@
 #ifndef register_problems_h
 #define register_problems_h
 
-#include "../include/structs.h"
+#include "structs.h"
 
 // Register problems to be available on the command line
 PetscErrorCode RegisterProblems_Hdiv(AppCtx app_ctx);

examples/Hdiv-mixed/include/setup-boundary.h

@@ -1,11 +1,11 @@
-#ifndef register_boundary_h
-#define register_boundary_h
+#ifndef setup_boundary_h
+#define setup_boundary_h
 
 #include <petsc.h>
 #include <petscdmplex.h>
 #include <petscsys.h>
 #include <ceed.h>
-#include "../include/structs.h"
+#include "structs.h"
 
 // ---------------------------------------------------------------------------
 // Create boundary label
@@ -19,4 +19,7 @@ PetscErrorCode DMAddBoundariesDirichlet(DM dm);
 PetscErrorCode BoundaryDirichletMMS(PetscInt dim, PetscReal t,
                                     const PetscReal coords[],
                                     PetscInt num_comp_u, PetscScalar *u, void *ctx);
-#endif // register_boundary_h
+PetscErrorCode DMAddBoundariesPressure(Ceed ceed, CeedData ceed_data,
+                                       AppCtx app_ctx, ProblemData *problem_data, DM dm,
+                                       CeedVector bc_pressure);
+#endif // setup_boundary_h

examples/Hdiv-mixed/include/setup-dm.h

@@ -5,7 +5,7 @@
 #include <petscdmplex.h>
 #include <petscsys.h>
 #include <ceed.h>
-#include "../include/structs.h"
+#include "structs.h"
 
 // ---------------------------------------------------------------------------
 // Set-up DM

examples/Hdiv-mixed/include/setup-libceed.h

@@ -1,7 +1,7 @@
 #ifndef setuplibceed_h
 #define setuplibceed_h
 
-#include "../include/structs.h"
+#include "structs.h"
 
 // Convert PETSc MemType to libCEED MemType
 CeedMemType MemTypeP2C(PetscMemType mtype);

examples/Hdiv-mixed/include/setup-solvers.h

@@ -6,7 +6,7 @@
 
 #include "structs.h"
 
-PetscErrorCode SetupCommonCtx(MPI_Comm comm, DM dm, Ceed ceed,
+PetscErrorCode SetupCommonCtx(DM dm, Ceed ceed,
                               CeedData ceed_data,
                               OperatorApplyContext op_apply_ctx);
 PetscErrorCode SetupJacobianOperatorCtx(CeedData ceed_data,
@@ -24,7 +24,7 @@ PetscErrorCode PDESolver(CeedData ceed_data, VecType vec_type, SNES snes,
 PetscErrorCode ComputeL2Error(CeedData ceed_data, Vec U, CeedVector target,
                               CeedScalar *l2_error_u, CeedScalar *l2_error_p,
                               OperatorApplyContext op_apply_ctx);
-PetscErrorCode PrintOutput(MPI_Comm comm, Ceed ceed,
+PetscErrorCode PrintOutput(Ceed ceed,
                            CeedMemType mem_type_backend,
                            SNES snes, KSP ksp,
                            Vec U, CeedScalar l2_error_u,

examples/Hdiv-mixed/include/structs.h

@@ -7,13 +7,14 @@
 // Application context from user command line options
 typedef struct AppCtx_ *AppCtx;
 struct AppCtx_ {
+  MPI_Comm          comm;
   // Degree of polynomial (1 only), extra quadrature pts
   PetscInt          degree;
   PetscInt          q_extra;
   // For applying traction BCs
-  PetscInt          bc_traction_count;
-  PetscInt          bc_traction_faces[16];
-  PetscScalar       bc_traction_vector[16][3];
+  PetscInt          bc_pressure_count;
+  PetscInt          bc_pressure_faces[16];
+  PetscScalar       bc_pressure_vector[16][3];
   // Problem type arguments
   PetscFunctionList problems;
   char              problem_name[PETSC_MAX_PATH_LEN];
@@ -23,12 +24,12 @@ struct AppCtx_ {
 // libCEED data struct
 typedef struct CeedData_ *CeedData;
 struct CeedData_ {
-  CeedBasis            basis_x, basis_u, basis_p;
+  CeedBasis            basis_x, basis_u, basis_p, basis_u_face;
   CeedElemRestriction  elem_restr_x, elem_restr_u, elem_restr_U_i,
                        elem_restr_p;
   CeedQFunction        qf_residual, qf_jacobian, qf_error;
   CeedOperator         op_residual, op_jacobian, op_error;
-  CeedVector           x_ceed, y_ceed;
+  CeedVector           x_ceed, y_ceed, x_coord;
   CeedQFunctionContext pq2d_context;
 };
 
@@ -72,13 +73,12 @@ struct OperatorApplyContext_ {
 
 // Problem specific data
 typedef struct {
-  CeedQFunctionUser setup_rhs, residual, jacobian, setup_error,
-                    setup_true, setup_face_geo;
-  const char        *setup_rhs_loc, *residual_loc, *jacobian_loc,
-        *setup_error_loc, *setup_true_loc, *setup_face_geo_loc;
+  CeedQFunctionUser force, residual, jacobian, error,
+                    setup_true, bc_pressure;
+  const char        *force_loc, *residual_loc, *jacobian_loc,
+        *error_loc, *setup_true_loc, *bc_pressure_loc;
   CeedQuadMode      quadrature_mode;
   CeedInt           elem_node, dim, q_data_size_face;
-  PetscErrorCode    (*setup_ctx)(Ceed, CeedData, Physics);
 
 } ProblemData;
 

examples/Hdiv-mixed/main.c

@@ -23,9 +23,11 @@
 //
 // Build with: make
 // Run with:
-//     ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
-//     ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4
-//     ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
+//   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4
+//   ./main -pc_type svd -problem darcy3d -dm_plex_dim 3 -dm_plex_box_faces 4,4,4
+//   ./main -pc_type svd -problem darcy3d -dm_plex_filename /path to the mesh file
+//   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1
+//   ./main -pc_type svd -problem darcy2d -dm_plex_dim 2 -dm_plex_box_faces 4,4 -bc_pressure 1,2,3,4
 const char help[] = "Solve H(div)-mixed problem using PETSc and libCEED\n";
 
 #include "main.h"
@@ -62,7 +64,8 @@ int main(int argc, char **argv) {
 
   // -- Process general command line options
   MPI_Comm comm = PETSC_COMM_WORLD;
-  PetscCall( ProcessCommandLineOptions(comm, app_ctx) );
+  app_ctx->comm = comm;
+  PetscCall( ProcessCommandLineOptions(app_ctx) );
 
   // ---------------------------------------------------------------------------
   // Choose the problem from the list of registered problems
@@ -107,12 +110,12 @@ int main(int argc, char **argv) {
 
 
   // ---------------------------------------------------------------------------
-  // Create local RHS vector
+  // Create local Force vector
   // ---------------------------------------------------------------------------
   Vec F_loc;
   PetscInt F_loc_size;
   CeedScalar *f;
-  CeedVector force_ceed, target;
+  CeedVector force_ceed, target, bc_pressure;
   PetscMemType force_mem_type;
   PetscCall( DMCreateLocalVector(dm, &F_loc) );
   // Local size for libCEED
@@ -121,15 +124,18 @@ int main(int argc, char **argv) {
   PetscCall( VecGetArrayAndMemType(F_loc, &f, &force_mem_type) );
   CeedVectorCreate(ceed, F_loc_size, &force_ceed);
   CeedVectorSetArray(force_ceed, MemTypeP2C(force_mem_type), CEED_USE_POINTER, f);
-
+  CeedVectorCreate(ceed, F_loc_size, &bc_pressure);
+  CeedVectorSetArray(bc_pressure, MemTypeP2C(force_mem_type), CEED_USE_POINTER,
+                     f);
   // ---------------------------------------------------------------------------
   // Setup libCEED - Compute local F and true solution (target)
   // ---------------------------------------------------------------------------
   // -- Set up libCEED objects
   PetscCall( SetupLibceed(dm, ceed, app_ctx, problem_data,
                           F_loc_size, ceed_data, force_ceed, &target) );
   //CeedVectorView(force_ceed, "%12.8f", stdout);
-
+  PetscCall( DMAddBoundariesPressure(ceed, ceed_data, app_ctx, problem_data, dm,
+                                     bc_pressure) );
   // ---------------------------------------------------------------------------
   // Create global F
   // ---------------------------------------------------------------------------
@@ -147,9 +153,10 @@ int main(int argc, char **argv) {
   SNES snes;
   KSP ksp;
   Vec U;
+  op_apply_ctx->comm = comm;
   PetscCall( SNESCreate(comm, &snes) );
   PetscCall( SNESGetKSP(snes, &ksp) );
-  PetscCall( SetupCommonCtx(comm, dm, ceed, ceed_data, op_apply_ctx) );
+  PetscCall( SetupCommonCtx(dm, ceed, ceed_data, op_apply_ctx) );
   PetscCall( PDESolver(ceed_data, vec_type, snes, ksp, F, &U, op_apply_ctx) );
   //VecView(U, PETSC_VIEWER_STDOUT_WORLD);
 
@@ -163,7 +170,7 @@ int main(int argc, char **argv) {
   // ---------------------------------------------------------------------------
   // Print output results
   // ---------------------------------------------------------------------------
-  PetscCall( PrintOutput(comm, ceed, mem_type_backend,
+  PetscCall( PrintOutput(ceed, mem_type_backend,
                          snes, ksp, U, l2_error_u, l2_error_p, app_ctx) );
 
   // ---------------------------------------------------------------------------
@@ -201,6 +208,7 @@ int main(int argc, char **argv) {
 
   // Free libCEED objects
   CeedVectorDestroy(&force_ceed);
+  CeedVectorDestroy(&bc_pressure);
   CeedVectorDestroy(&target);
   PetscCall( CeedDataDestroy(ceed_data) );
   CeedDestroy(&ceed);