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:   PetscReal     refinementLimit;              /* The largest allowable cell volume */
 13:   PetscBool     cellSimplex;                  /* Use simplices or hexes */
 14:   char          filename[PETSC_MAX_PATH_LEN]; /* Import mesh from file */
 15:   PetscBool     testPartition;                /* Use a fixed partitioning for testing */
 16:   PetscInt      overlap;                      /* The cell overlap to use during partitioning */
 17: } AppCtx;

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

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

 35:   PetscOptionsBegin(comm, "", "Meshing Problem Options", "DMPLEX");
 36:   PetscOptionsInt("-debug", "The debugging level", "ex1.c", options->debug, &options->debug, NULL);
 37:   PetscOptionsInt("-dim", "The topological mesh dimension", "ex1.c", options->dim, &options->dim, NULL);
 38:   PetscOptionsBool("-interpolate", "Generate intermediate mesh elements", "ex1.c", options->interpolate, &options->interpolate, NULL);
 39:   PetscOptionsReal("-refinement_limit", "The largest allowable cell volume", "ex1.c", options->refinementLimit, &options->refinementLimit, NULL);
 40:   PetscOptionsBool("-cell_simplex", "Use simplices if true, otherwise hexes", "ex1.c", options->cellSimplex, &options->cellSimplex, NULL);
 41:   PetscOptionsString("-filename", "The mesh file", "ex1.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);
 42:   PetscOptionsBool("-test_partition", "Use a fixed partition for testing", "ex1.c", options->testPartition, &options->testPartition, NULL);
 43:   PetscOptionsInt("-overlap", "The cell overlap for partitioning", "ex1.c", options->overlap, &options->overlap, NULL);
 44:   PetscOptionsEnd();

 46:   PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);
 47:   return(0);
 48: };

 52: PetscErrorCode CreateMesh(MPI_Comm comm, AppCtx *user, DM *dm)
 53: {
 54:   PetscInt       dim             = user->dim;
 55:   PetscBool      interpolate     = user->interpolate;
 56:   PetscReal      refinementLimit = user->refinementLimit;
 57:   PetscBool      cellSimplex     = user->cellSimplex;
 58:   const char    *filename        = user->filename;
 59:   PetscInt       triSizes_n2[2]  = {4, 4};
 60:   PetscInt       triPoints_n2[8] = {3, 5, 6, 7, 0, 1, 2, 4};
 61:   PetscInt       triSizes_n8[8]  = {1, 1, 1, 1, 1, 1, 1, 1};
 62:   PetscInt       triPoints_n8[8] = {0, 1, 2, 3, 4, 5, 6, 7};
 63:   PetscInt       quadSizes[2]    = {2, 2};
 64:   PetscInt       quadPoints[4]   = {2, 3, 0, 1};
 65:   const PetscInt cells[3]        = {2, 2, 2};
 66:   size_t         len;
 67:   PetscMPIInt    rank, numProcs;

 71:   PetscLogEventBegin(user->createMeshEvent,0,0,0,0);
 72:   MPI_Comm_rank(comm, &rank);
 73:   MPI_Comm_size(comm, &numProcs);
 74:   PetscStrlen(filename, &len);
 75:   if (len)              {DMPlexCreateFromFile(comm, filename, interpolate, dm);}
 76:   else if (cellSimplex) {DMPlexCreateBoxMesh(comm, dim, interpolate, dm);}
 77:   else                  {DMPlexCreateHexBoxMesh(comm, dim, cells, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, dm);}
 78:   {
 79:     DM refinedMesh     = NULL;
 80:     DM distributedMesh = NULL;

 82:     if (user->testPartition) {
 83:       const PetscInt  *sizes = NULL;
 84:       const PetscInt  *points = NULL;
 85:       PetscPartitioner part;

 87:       if (!rank) {
 88:         if (dim == 2 && cellSimplex && numProcs == 2) {
 89:            sizes = triSizes_n2; points = triPoints_n2;
 90:         } else if (dim == 2 && cellSimplex && numProcs == 8) {
 91:           sizes = triSizes_n8; points = triPoints_n8;
 92:         } else if (dim == 2 && !cellSimplex && numProcs == 2) {
 93:           sizes = quadSizes; points = quadPoints;
 94:         }
 95:       }
 96:       DMPlexGetPartitioner(*dm, &part);
 97:       PetscPartitionerSetType(part, PETSCPARTITIONERSHELL);
 98:       PetscPartitionerShellSetPartition(part, numProcs, sizes, points);
 99:     }
100:     /* Distribute mesh over processes */
101:     DMPlexDistribute(*dm, 0, NULL, &distributedMesh);
102:     if (distributedMesh) {
103:       DMDestroy(dm);
104:       *dm  = distributedMesh;
105:     }
106:     /* Refine mesh using a volume constraint */
107:     DMPlexSetRefinementUniform(*dm, PETSC_FALSE);
108:     DMPlexSetRefinementLimit(*dm, refinementLimit);
109:     DMRefine(*dm, comm, &refinedMesh);
110:     if (refinedMesh) {
111:       DMDestroy(dm);
112:       *dm  = refinedMesh;
113:     }
114:   }
115:   DMSetFromOptions(*dm);
116:   if (user->overlap) {
117:     DM overlapMesh = NULL;
118:     /* Add the level-1 overlap to refined mesh */
119:     DMPlexDistributeOverlap(*dm, 1, NULL, &overlapMesh);
120:     if (overlapMesh) {
121:       DMView(overlapMesh, PETSC_VIEWER_STDOUT_WORLD);
122:       DMDestroy(dm);
123:       *dm = overlapMesh;
124:     }
125:   }
126:   PetscObjectSetName((PetscObject) *dm, "Simplicial Mesh");
127:   DMViewFromOptions(*dm, NULL, "-dm_view");
128:   PetscLogEventEnd(user->createMeshEvent,0,0,0,0);
129:   user->dm = *dm;
130:   return(0);
131: }

135: int main(int argc, char **argv)
136: {
137:   AppCtx         user;                 /* user-defined work context */

140:   PetscInitialize(&argc, &argv, NULL, help);
141:   ProcessOptions(PETSC_COMM_WORLD, &user);
142:   CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);
143:   DMDestroy(&user.dm);
144:   PetscFinalize();
145:   return 0;
146: }