1: static char help[] = "Run C version of TetGen to construct and refine a mesh\n\n";

  3: #include <petscdmplex.h>

  5: typedef struct {
  6:   DM            dm;                /* REQUIRED in order to use SNES evaluation functions */
  7:   PetscInt      debug;             /* The debugging level */
  8:   PetscLogEvent createMeshEvent;
  9:   /* Domain and mesh definition */
 10:   PetscInt  dim;                   /* The topological mesh dimension */
 11:   PetscBool interpolate;           /* Generate intermediate mesh elements */
 12:   PetscBool refinementUniform;     /* Uniformly refine the mesh */
 13:   PetscReal refinementLimit;       /* The largest allowable cell volume */
 14:   PetscBool cellSimplex;           /* Use simplices or hexes */
 15:   char      filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
 16: } AppCtx;

 20: PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 21: {

 25:   options->debug             = 0;
 26:   options->dim               = 2;
 27:   options->interpolate       = PETSC_FALSE;
 28:   options->refinementUniform = PETSC_FALSE;
 29:   options->refinementLimit   = 0.0;
 30:   options->cellSimplex       = PETSC_TRUE;
 31:   options->filename[0]       = '\0';

 33:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
 34:   PetscOptionsInt("-debug", "The debugging level", "ex1.c", options->debug, &options->debug, NULL);
 35:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);
 36:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex1.c", options->interpolate, &options->interpolate, NULL);
 37:   PetscOptionsBool("-refinement_uniform", "Uniformly refine the mesh", "ex1.c", options->refinementUniform, &options->refinementUniform, NULL);
 38:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex1.c", options->refinementLimit, &options->refinementLimit, NULL);
 39:   PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex1.c", options->cellSimplex, &options->cellSimplex, NULL);
 40:   PetscOptionsString("-filename", "The mesh file", "ex7.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
 41:   PetscOptionsEnd();

 43:   PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);
 44:   return(0);
 45: };

 49: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 50: {
 51:   PetscInt       dim               = user->dim;
 52:   PetscBool      interpolate       = user->interpolate;
 53:   PetscBool      refinementUniform = user->refinementUniform;
 54:   PetscReal      refinementLimit   = user->refinementLimit;
 55:   PetscBool      cellSimplex       = user->cellSimplex;
 56:   const char    *filename          = user->filename;
 57:   const char    *partitioner       = "chaco";
 58:   size_t         len;
 59:   PetscMPIInt    rank;

 63:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
 64:   MPI_Comm_rank(comm, &rank);
 65:   PetscStrlen(filename, &len);
 66:   if (len) {
 67:     const char *extGmsh = ".msh";
 68:     size_t      slen;
 69:     PetscBool   isGmsh;

 71:     PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);
 72:     if (isGmsh) {
 73:       PetscViewer viewer;

 75:       PetscViewerCreate(comm, &viewer);
 76:       PetscViewerSetType(viewer, PETSCVIEWERASCII);
 77:       PetscViewerFileSetMode(viewer, FILE_MODE_READ);
 78:       PetscViewerFileSetName(viewer, filename);
 79:       DMPlexCreateGmsh(comm, viewer, interpolate, dm);
 80:       PetscViewerDestroy(&viewer);
 81:     } else {
 82:       DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);
 83:     }
 84:   } else if (cellSimplex) {
 85:     DMPlexCreateBoxMesh(comm, dim, interpolate, dm);
 86:   } else {
 87:     const PetscInt cells[3] = {2, 2, 2};

 89:     DMPlexCreateHexBoxMesh(comm, dim, cells, PETSC_FALSE, PETSC_FALSE, PETSC_FALSE, dm);
 90:   }
 91:   {
 92:     DM refinedMesh     = NULL;
 93:     DM distributedMesh = NULL;

 95:     /* Refine mesh using a volume constraint */
 96:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
 97:     DMPlexSetRefinementLimit(*dm, refinementLimit);
 98:     DMRefine(*dm, comm, &refinedMesh);
 99:     if (refinedMesh) {
100:       DMDestroy(dm);
101:       *dm  = refinedMesh;
102:     }
103:     /* Distribute mesh over processes */
104:     DMPlexDistribute(*dm, partitioner, 0, NULL, &distributedMesh);
105:     if (distributedMesh) {
106:       DMViewFromOptions(distributedMesh, NULL, "-dm_view");
107:       DMDestroy(dm);
108:       *dm  = distributedMesh;
109:     }
110:     if (refinementUniform) {
111:       DMPlexSetRefinementUniform(*dm, refinementUniform);
112:       DMRefine(*dm, comm, &refinedMesh);
113:       if (refinedMesh) {
114:         DMDestroy(dm);
115:         *dm  = refinedMesh;
116:       }
117:     }
118:   }
119:   PetscObjectSetName((PetscObject) *dm, "Simplical Mesh");
120:   DMViewFromOptions(*dm, NULL, "-dm_view");
121:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
122:   user->dm = *dm;
123:   return(0);
124: }

128: int main(int argc, char **argv)
129: {
130:   AppCtx         user;                 /* user-defined work context */

133:   PetscInitialize(&argc, &argv, NULL, help);
134:   ProcessOptions(PETSC_COMM_WORLD, &user);
135:   CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);
136:   DMDestroy(&user.dm);
137:   PetscFinalize();
138:   return 0;
139: }