Actual source code: ex1.c

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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:   PetscBool     testShape;                    /* Test the cell shape quality */
 18: } AppCtx;

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

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

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

 49:   PetscLogEventRegister("CreateMesh",          DM_CLASSID,   &options->createMeshEvent);
 50:   return(0);
 51: };

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

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

 85:     if (user->testPartition) {
 86:       const PetscInt  *sizes = NULL;
 87:       const PetscInt  *points = NULL;
 88:       PetscPartitioner part;

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

136: typedef struct ex1_stats
137: {
138:   PetscReal min, max, sum, squaresum;
139:   PetscInt  count;
140: }
141: ex1_stats_t;

143: static void ex1_stats_reduce(void *a, void *b, int * len, MPI_Datatype *datatype)
144: {
145:   PetscInt i, N = *len;

147:   for (i = 0; i < N; i++) {
148:     ex1_stats_t *A = (ex1_stats_t *) a;
149:     ex1_stats_t *B = (ex1_stats_t *) b;

151:     B->min = PetscMin(A->min,B->min);
152:     B->max = PetscMax(A->max,B->max);
153:     B->sum += A->sum;
154:     B->squaresum += A->squaresum;
155:     B->count += A->count;
156:   }
157: }

161: static PetscErrorCode TestCellShape(DM dm)
162: {
163:   PetscMPIInt    rank;
164:   PetscInt       dim, c, cStart, cEnd, count = 0;
165:   ex1_stats_t    stats, globalStats;
166:   PetscReal      *J, *invJ, min = 0, max = 0, mean = 0, stdev = 0;
167:   MPI_Comm       comm = PetscObjectComm((PetscObject)dm);
168:   DM             dmCoarse;

172:   stats.min = PETSC_MAX_REAL;
173:   stats.max = PETSC_MIN_REAL;
174:   stats.sum = stats.squaresum = 0.;
175:   stats.count = 0;

177:   DMGetDimension(dm,&dim);

179:   PetscMalloc2(dim * dim, &J, dim * dim, &invJ);

181:   DMPlexGetHeightStratum(dm,0,&cStart,&cEnd);
182:   for (c = cStart; c < cEnd; c++) {
183:     PetscInt  i;
184:     PetscReal frobJ = 0., frobInvJ = 0., cond2, cond, detJ;

186:     DMPlexComputeCellGeometryAffineFEM(dm,c,NULL,J,invJ,&detJ);

188:     for (i = 0; i < dim * dim; i++) {
189:       frobJ += J[i] * J[i];
190:       frobInvJ += invJ[i] * invJ[i];
191:     }
192:     cond2 = frobJ * frobInvJ;
193:     cond  = PetscSqrtReal(cond2);

195:     stats.min = PetscMin(stats.min,cond);
196:     stats.max = PetscMax(stats.max,cond);
197:     stats.sum += cond;
198:     stats.squaresum += cond2;
199:     stats.count++;
200:   }

202:   {
203:     PetscMPIInt    blockLengths[2] = {4,1};
204:     MPI_Aint       blockOffsets[2] = {offsetof(ex1_stats_t,min),offsetof(ex1_stats_t,count)};
205:     MPI_Datatype   blockTypes[2]   = {MPIU_REAL,MPIU_INT}, statType;
206:     MPI_Op         statReduce;

208:     MPI_Type_create_struct(2,blockLengths,blockOffsets,blockTypes,&statType);
209:     MPI_Type_commit(&statType);
210:     MPI_Op_create(ex1_stats_reduce, PETSC_TRUE, &statReduce);
211:     MPI_Reduce(&stats,&globalStats,1,statType,statReduce,0,comm);
212:     MPI_Op_free(&statReduce);
213:     MPI_Type_free(&statType);
214:   }

216:   MPI_Comm_rank(comm,&rank);
217:   if (!rank) {
218:     count = globalStats.count;
219:     min = globalStats.min;
220:     max = globalStats.max;
221:     mean = globalStats.sum / globalStats.count;
222:     stdev = PetscSqrtReal(globalStats.squaresum / globalStats.count - mean * mean);
223:   }
224:   PetscPrintf(comm,"Mesh with %d cells, shape condition numbers: min = %g, max = %g, mean = %g, stddev = %g\n", count, (double) min, (double) max, (double) mean, (double) stdev);

226:   PetscFree2(J,invJ);

228:   DMGetCoarseDM(dm,&dmCoarse);
229:   if (dmCoarse) {
230:     TestCellShape(dmCoarse);
231:   }

233:   return(0);
234: }

238: int main(int argc, char **argv)
239: {
240:   AppCtx         user;                 /* user-defined work context */

243:   PetscInitialize(&argc, &argv, NULL, help);
244:   ProcessOptions(PETSC_COMM_WORLD, &user);
245:   CreateMesh(PETSC_COMM_WORLD, &user, &user.dm);
246:   if (user.testShape) {
247:     TestCellShape(user.dm);
248:   }
249:   DMDestroy(&user.dm);
250:   PetscFinalize();
251:   return 0;
252: }