Actual source code: ex1.c
petsc-3.7.7 2017-09-25
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: }