1 #define PETSCDM_DLL
2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/
3 #include <petsc/private/hashseti.h> /*I "petscdmplex.h" I*/
4 #include <petscsf.h>
5
6 PetscLogEvent DMPLEX_CreateFromFile, DMPLEX_BuildFromCellList, DMPLEX_BuildCoordinatesFromCellList;
7
8 /*@
9 DMPlexCreateDoublet - Creates a mesh of two cells of the specified type, optionally with later refinement.
10
11 Collective
12
13 Input Parameters:
14 + comm - The communicator for the DM object
15 . dim - The spatial dimension
16 . simplex - Flag for simplicial cells, otherwise they are tensor product cells
17 . interpolate - Flag to create intermediate mesh pieces (edges, faces)
18 . refinementUniform - Flag for uniform parallel refinement
19 - refinementLimit - A nonzero number indicates the largest admissible volume for a refined cell
20
21 Output Parameter:
22 . dm - The DM object
23
24 Level: beginner
25
26 .seealso: DMSetType(), DMCreate()
27 @*/
DMPlexCreateDoublet(MPI_Comm comm,PetscInt dim,PetscBool simplex,PetscBool interpolate,PetscBool refinementUniform,PetscReal refinementLimit,DM * newdm)28 PetscErrorCode DMPlexCreateDoublet(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscBool interpolate, PetscBool refinementUniform, PetscReal refinementLimit, DM *newdm)
29 {
30 DM dm;
31 PetscInt p;
32 PetscMPIInt rank;
33 PetscErrorCode ierr;
34
35 PetscFunctionBegin;
36 ierr = DMCreate(comm, &dm);CHKERRQ(ierr);
37 ierr = DMSetType(dm, DMPLEX);CHKERRQ(ierr);
38 ierr = DMSetDimension(dm, dim);CHKERRQ(ierr);
39 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
40 switch (dim) {
41 case 2:
42 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "triangular");CHKERRQ(ierr);}
43 else {ierr = PetscObjectSetName((PetscObject) dm, "quadrilateral");CHKERRQ(ierr);}
44 break;
45 case 3:
46 if (simplex) {ierr = PetscObjectSetName((PetscObject) dm, "tetrahedral");CHKERRQ(ierr);}
47 else {ierr = PetscObjectSetName((PetscObject) dm, "hexahedral");CHKERRQ(ierr);}
48 break;
49 default:
50 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim);
51 }
52 if (rank) {
53 PetscInt numPoints[2] = {0, 0};
54 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
55 } else {
56 switch (dim) {
57 case 2:
58 if (simplex) {
59 PetscInt numPoints[2] = {4, 2};
60 PetscInt coneSize[6] = {3, 3, 0, 0, 0, 0};
61 PetscInt cones[6] = {2, 3, 4, 5, 4, 3};
62 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
63 PetscScalar vertexCoords[8] = {-0.5, 0.5, 0.0, 0.0, 0.0, 1.0, 0.5, 0.5};
64 PetscInt markerPoints[8] = {2, 1, 3, 1, 4, 1, 5, 1};
65
66 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
67 for (p = 0; p < 4; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
68 } else {
69 PetscInt numPoints[2] = {6, 2};
70 PetscInt coneSize[8] = {4, 4, 0, 0, 0, 0, 0, 0};
71 PetscInt cones[8] = {2, 3, 4, 5, 3, 6, 7, 4};
72 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
73 PetscScalar vertexCoords[12] = {-1.0, -0.5, 0.0, -0.5, 0.0, 0.5, -1.0, 0.5, 1.0, -0.5, 1.0, 0.5};
74
75 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
76 }
77 break;
78 case 3:
79 if (simplex) {
80 PetscInt numPoints[2] = {5, 2};
81 PetscInt coneSize[7] = {4, 4, 0, 0, 0, 0, 0};
82 PetscInt cones[8] = {4, 3, 5, 2, 5, 3, 4, 6};
83 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
84 PetscScalar vertexCoords[15] = {-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0};
85 PetscInt markerPoints[10] = {2, 1, 3, 1, 4, 1, 5, 1, 6, 1};
86
87 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
88 for (p = 0; p < 5; ++p) {ierr = DMSetLabelValue(dm, "marker", markerPoints[p*2], markerPoints[p*2+1]);CHKERRQ(ierr);}
89 } else {
90 PetscInt numPoints[2] = {12, 2};
91 PetscInt coneSize[14] = {8, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
92 PetscInt cones[16] = {2, 3, 4, 5, 6, 7, 8, 9, 5, 4, 10, 11, 7, 12, 13, 8};
93 PetscInt coneOrientations[16] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
94 PetscScalar vertexCoords[36] = {-1.0, -0.5, -0.5, -1.0, 0.5, -0.5, 0.0, 0.5, -0.5, 0.0, -0.5, -0.5,
95 -1.0, -0.5, 0.5, 0.0, -0.5, 0.5, 0.0, 0.5, 0.5, -1.0, 0.5, 0.5,
96 1.0, 0.5, -0.5, 1.0, -0.5, -0.5, 1.0, -0.5, 0.5, 1.0, 0.5, 0.5};
97
98 ierr = DMPlexCreateFromDAG(dm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
99 }
100 break;
101 default:
102 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Cannot make meshes for dimension %D", dim);
103 }
104 }
105 *newdm = dm;
106 if (refinementLimit > 0.0) {
107 DM rdm;
108 const char *name;
109
110 ierr = DMPlexSetRefinementUniform(*newdm, PETSC_FALSE);CHKERRQ(ierr);
111 ierr = DMPlexSetRefinementLimit(*newdm, refinementLimit);CHKERRQ(ierr);
112 ierr = DMRefine(*newdm, comm, &rdm);CHKERRQ(ierr);
113 ierr = PetscObjectGetName((PetscObject) *newdm, &name);CHKERRQ(ierr);
114 ierr = PetscObjectSetName((PetscObject) rdm, name);CHKERRQ(ierr);
115 ierr = DMDestroy(newdm);CHKERRQ(ierr);
116 *newdm = rdm;
117 }
118 if (interpolate) {
119 DM idm;
120
121 ierr = DMPlexInterpolate(*newdm, &idm);CHKERRQ(ierr);
122 ierr = DMDestroy(newdm);CHKERRQ(ierr);
123 *newdm = idm;
124 }
125 {
126 DM refinedMesh = NULL;
127 DM distributedMesh = NULL;
128
129 /* Distribute mesh over processes */
130 ierr = DMPlexDistribute(*newdm, 0, NULL, &distributedMesh);CHKERRQ(ierr);
131 if (distributedMesh) {
132 ierr = DMDestroy(newdm);CHKERRQ(ierr);
133 *newdm = distributedMesh;
134 }
135 if (refinementUniform) {
136 ierr = DMPlexSetRefinementUniform(*newdm, refinementUniform);CHKERRQ(ierr);
137 ierr = DMRefine(*newdm, comm, &refinedMesh);CHKERRQ(ierr);
138 if (refinedMesh) {
139 ierr = DMDestroy(newdm);CHKERRQ(ierr);
140 *newdm = refinedMesh;
141 }
142 }
143 }
144 PetscFunctionReturn(0);
145 }
146
147 /*@
148 DMPlexCreateSquareBoundary - Creates a 1D mesh the is the boundary of a square lattice.
149
150 Collective
151
152 Input Parameters:
153 + comm - The communicator for the DM object
154 . lower - The lower left corner coordinates
155 . upper - The upper right corner coordinates
156 - edges - The number of cells in each direction
157
158 Output Parameter:
159 . dm - The DM object
160
161 Note: Here is the numbering returned for 2 cells in each direction:
162 $ 18--5-17--4--16
163 $ | | |
164 $ 6 10 3
165 $ | | |
166 $ 19-11-20--9--15
167 $ | | |
168 $ 7 8 2
169 $ | | |
170 $ 12--0-13--1--14
171
172 Level: beginner
173
174 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateCubeBoundary(), DMSetType(), DMCreate()
175 @*/
DMPlexCreateSquareBoundary(DM dm,const PetscReal lower[],const PetscReal upper[],const PetscInt edges[])176 PetscErrorCode DMPlexCreateSquareBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[])
177 {
178 const PetscInt numVertices = (edges[0]+1)*(edges[1]+1);
179 const PetscInt numEdges = edges[0]*(edges[1]+1) + (edges[0]+1)*edges[1];
180 PetscInt markerTop = 1;
181 PetscInt markerBottom = 1;
182 PetscInt markerRight = 1;
183 PetscInt markerLeft = 1;
184 PetscBool markerSeparate = PETSC_FALSE;
185 Vec coordinates;
186 PetscSection coordSection;
187 PetscScalar *coords;
188 PetscInt coordSize;
189 PetscMPIInt rank;
190 PetscInt v, vx, vy;
191 PetscErrorCode ierr;
192
193 PetscFunctionBegin;
194 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr);
195 if (markerSeparate) {
196 markerTop = 3;
197 markerBottom = 1;
198 markerRight = 2;
199 markerLeft = 4;
200 }
201 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
202 if (!rank) {
203 PetscInt e, ex, ey;
204
205 ierr = DMPlexSetChart(dm, 0, numEdges+numVertices);CHKERRQ(ierr);
206 for (e = 0; e < numEdges; ++e) {
207 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr);
208 }
209 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
210 for (vx = 0; vx <= edges[0]; vx++) {
211 for (ey = 0; ey < edges[1]; ey++) {
212 PetscInt edge = vx*edges[1] + ey + edges[0]*(edges[1]+1);
213 PetscInt vertex = ey*(edges[0]+1) + vx + numEdges;
214 PetscInt cone[2];
215
216 cone[0] = vertex; cone[1] = vertex+edges[0]+1;
217 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
218 if (vx == edges[0]) {
219 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr);
220 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr);
221 if (ey == edges[1]-1) {
222 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr);
223 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerRight);CHKERRQ(ierr);
224 }
225 } else if (vx == 0) {
226 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr);
227 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr);
228 if (ey == edges[1]-1) {
229 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr);
230 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerLeft);CHKERRQ(ierr);
231 }
232 }
233 }
234 }
235 for (vy = 0; vy <= edges[1]; vy++) {
236 for (ex = 0; ex < edges[0]; ex++) {
237 PetscInt edge = vy*edges[0] + ex;
238 PetscInt vertex = vy*(edges[0]+1) + ex + numEdges;
239 PetscInt cone[2];
240
241 cone[0] = vertex; cone[1] = vertex+1;
242 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
243 if (vy == edges[1]) {
244 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr);
245 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr);
246 if (ex == edges[0]-1) {
247 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr);
248 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerTop);CHKERRQ(ierr);
249 }
250 } else if (vy == 0) {
251 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr);
252 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr);
253 if (ex == edges[0]-1) {
254 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr);
255 ierr = DMSetLabelValue(dm, "Face Sets", cone[1], markerBottom);CHKERRQ(ierr);
256 }
257 }
258 }
259 }
260 }
261 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
262 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
263 /* Build coordinates */
264 ierr = DMSetCoordinateDim(dm, 2);CHKERRQ(ierr);
265 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
266 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
267 ierr = PetscSectionSetChart(coordSection, numEdges, numEdges + numVertices);CHKERRQ(ierr);
268 ierr = PetscSectionSetFieldComponents(coordSection, 0, 2);CHKERRQ(ierr);
269 for (v = numEdges; v < numEdges+numVertices; ++v) {
270 ierr = PetscSectionSetDof(coordSection, v, 2);CHKERRQ(ierr);
271 ierr = PetscSectionSetFieldDof(coordSection, v, 0, 2);CHKERRQ(ierr);
272 }
273 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
274 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
275 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
276 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
277 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
278 ierr = VecSetBlockSize(coordinates, 2);CHKERRQ(ierr);
279 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
280 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
281 for (vy = 0; vy <= edges[1]; ++vy) {
282 for (vx = 0; vx <= edges[0]; ++vx) {
283 coords[(vy*(edges[0]+1)+vx)*2+0] = lower[0] + ((upper[0] - lower[0])/edges[0])*vx;
284 coords[(vy*(edges[0]+1)+vx)*2+1] = lower[1] + ((upper[1] - lower[1])/edges[1])*vy;
285 }
286 }
287 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
288 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
289 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
290 PetscFunctionReturn(0);
291 }
292
293 /*@
294 DMPlexCreateCubeBoundary - Creates a 2D mesh that is the boundary of a cubic lattice.
295
296 Collective
297
298 Input Parameters:
299 + comm - The communicator for the DM object
300 . lower - The lower left front corner coordinates
301 . upper - The upper right back corner coordinates
302 - faces - The number of faces in each direction (the same as the number of cells)
303
304 Output Parameter:
305 . dm - The DM object
306
307 Level: beginner
308
309 .seealso: DMPlexCreateBoxMesh(), DMPlexCreateSquareBoundary(), DMSetType(), DMCreate()
310 @*/
DMPlexCreateCubeBoundary(DM dm,const PetscReal lower[],const PetscReal upper[],const PetscInt faces[])311 PetscErrorCode DMPlexCreateCubeBoundary(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt faces[])
312 {
313 PetscInt vertices[3], numVertices;
314 PetscInt numFaces = 2*faces[0]*faces[1] + 2*faces[1]*faces[2] + 2*faces[0]*faces[2];
315 Vec coordinates;
316 PetscSection coordSection;
317 PetscScalar *coords;
318 PetscInt coordSize;
319 PetscMPIInt rank;
320 PetscInt v, vx, vy, vz;
321 PetscInt voffset, iface=0, cone[4];
322 PetscErrorCode ierr;
323
324 PetscFunctionBegin;
325 if ((faces[0] < 1) || (faces[1] < 1) || (faces[2] < 1)) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Must have at least 1 face per side");
326 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
327 vertices[0] = faces[0]+1; vertices[1] = faces[1]+1; vertices[2] = faces[2]+1;
328 numVertices = vertices[0]*vertices[1]*vertices[2];
329 if (!rank) {
330 PetscInt f;
331
332 ierr = DMPlexSetChart(dm, 0, numFaces+numVertices);CHKERRQ(ierr);
333 for (f = 0; f < numFaces; ++f) {
334 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr);
335 }
336 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
337
338 /* Side 0 (Top) */
339 for (vy = 0; vy < faces[1]; vy++) {
340 for (vx = 0; vx < faces[0]; vx++) {
341 voffset = numFaces + vertices[0]*vertices[1]*(vertices[2]-1) + vy*vertices[0] + vx;
342 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]+1; cone[3] = voffset+vertices[0];
343 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
344 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
345 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
346 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr);
347 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr);
348 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr);
349 iface++;
350 }
351 }
352
353 /* Side 1 (Bottom) */
354 for (vy = 0; vy < faces[1]; vy++) {
355 for (vx = 0; vx < faces[0]; vx++) {
356 voffset = numFaces + vy*(faces[0]+1) + vx;
357 cone[0] = voffset+1; cone[1] = voffset; cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]+1;
358 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
359 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
360 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
361 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr);
362 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr);
363 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+1, 1);CHKERRQ(ierr);
364 iface++;
365 }
366 }
367
368 /* Side 2 (Front) */
369 for (vz = 0; vz < faces[2]; vz++) {
370 for (vx = 0; vx < faces[0]; vx++) {
371 voffset = numFaces + vz*vertices[0]*vertices[1] + vx;
372 cone[0] = voffset; cone[1] = voffset+1; cone[2] = voffset+vertices[0]*vertices[1]+1; cone[3] = voffset+vertices[0]*vertices[1];
373 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
374 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
375 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
376 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr);
377 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr);
378 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr);
379 iface++;
380 }
381 }
382
383 /* Side 3 (Back) */
384 for (vz = 0; vz < faces[2]; vz++) {
385 for (vx = 0; vx < faces[0]; vx++) {
386 voffset = numFaces + vz*vertices[0]*vertices[1] + vertices[0]*(vertices[1]-1) + vx;
387 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset+vertices[0]*vertices[1]+1;
388 cone[2] = voffset+1; cone[3] = voffset;
389 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
390 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
391 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
392 ierr = DMSetLabelValue(dm, "marker", voffset+1, 1);CHKERRQ(ierr);
393 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr);
394 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+1, 1);CHKERRQ(ierr);
395 iface++;
396 }
397 }
398
399 /* Side 4 (Left) */
400 for (vz = 0; vz < faces[2]; vz++) {
401 for (vy = 0; vy < faces[1]; vy++) {
402 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0];
403 cone[0] = voffset; cone[1] = voffset+vertices[0]*vertices[1];
404 cone[2] = voffset+vertices[0]*vertices[1]+vertices[0]; cone[3] = voffset+vertices[0];
405 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
406 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
407 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
408 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr);
409 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[1]+0, 1);CHKERRQ(ierr);
410 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr);
411 iface++;
412 }
413 }
414
415 /* Side 5 (Right) */
416 for (vz = 0; vz < faces[2]; vz++) {
417 for (vy = 0; vy < faces[1]; vy++) {
418 voffset = numFaces + vz*vertices[0]*vertices[1] + vy*vertices[0] + faces[0];
419 cone[0] = voffset+vertices[0]*vertices[1]; cone[1] = voffset;
420 cone[2] = voffset+vertices[0]; cone[3] = voffset+vertices[0]*vertices[1]+vertices[0];
421 ierr = DMPlexSetCone(dm, iface, cone);CHKERRQ(ierr);
422 ierr = DMSetLabelValue(dm, "marker", iface, 1);CHKERRQ(ierr);
423 ierr = DMSetLabelValue(dm, "marker", voffset+0, 1);CHKERRQ(ierr);
424 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]+0, 1);CHKERRQ(ierr);
425 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+0, 1);CHKERRQ(ierr);
426 ierr = DMSetLabelValue(dm, "marker", voffset+vertices[0]*vertices[1]+vertices[0], 1);CHKERRQ(ierr);
427 iface++;
428 }
429 }
430 }
431 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
432 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
433 /* Build coordinates */
434 ierr = DMSetCoordinateDim(dm, 3);CHKERRQ(ierr);
435 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
436 ierr = PetscSectionSetChart(coordSection, numFaces, numFaces + numVertices);CHKERRQ(ierr);
437 for (v = numFaces; v < numFaces+numVertices; ++v) {
438 ierr = PetscSectionSetDof(coordSection, v, 3);CHKERRQ(ierr);
439 }
440 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
441 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
442 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
443 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
444 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
445 ierr = VecSetBlockSize(coordinates, 3);CHKERRQ(ierr);
446 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
447 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
448 for (vz = 0; vz <= faces[2]; ++vz) {
449 for (vy = 0; vy <= faces[1]; ++vy) {
450 for (vx = 0; vx <= faces[0]; ++vx) {
451 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+0] = lower[0] + ((upper[0] - lower[0])/faces[0])*vx;
452 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+1] = lower[1] + ((upper[1] - lower[1])/faces[1])*vy;
453 coords[((vz*(faces[1]+1)+vy)*(faces[0]+1)+vx)*3+2] = lower[2] + ((upper[2] - lower[2])/faces[2])*vz;
454 }
455 }
456 }
457 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
458 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
459 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
460 PetscFunctionReturn(0);
461 }
462
DMPlexCreateLineMesh_Internal(MPI_Comm comm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd,DM * dm)463 static PetscErrorCode DMPlexCreateLineMesh_Internal(MPI_Comm comm,PetscInt segments,PetscReal lower,PetscReal upper,DMBoundaryType bd,DM *dm)
464 {
465 PetscInt i,fStart,fEnd,numCells = 0,numVerts = 0;
466 PetscInt numPoints[2],*coneSize,*cones,*coneOrientations;
467 PetscScalar *vertexCoords;
468 PetscReal L,maxCell;
469 PetscBool markerSeparate = PETSC_FALSE;
470 PetscInt markerLeft = 1, faceMarkerLeft = 1;
471 PetscInt markerRight = 1, faceMarkerRight = 2;
472 PetscBool wrap = (bd == DM_BOUNDARY_PERIODIC || bd == DM_BOUNDARY_TWIST) ? PETSC_TRUE : PETSC_FALSE;
473 PetscMPIInt rank;
474 PetscErrorCode ierr;
475
476 PetscFunctionBegin;
477 PetscValidPointer(dm,4);
478
479 ierr = DMCreate(comm,dm);CHKERRQ(ierr);
480 ierr = DMSetType(*dm,DMPLEX);CHKERRQ(ierr);
481 ierr = DMSetDimension(*dm,1);CHKERRQ(ierr);
482 ierr = DMCreateLabel(*dm,"marker");CHKERRQ(ierr);
483 ierr = DMCreateLabel(*dm,"Face Sets");CHKERRQ(ierr);
484
485 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
486 if (!rank) numCells = segments;
487 if (!rank) numVerts = segments + (wrap ? 0 : 1);
488
489 numPoints[0] = numVerts ; numPoints[1] = numCells;
490 ierr = PetscMalloc4(numCells+numVerts,&coneSize,numCells*2,&cones,numCells+numVerts,&coneOrientations,numVerts,&vertexCoords);CHKERRQ(ierr);
491 ierr = PetscArrayzero(coneOrientations,numCells+numVerts);CHKERRQ(ierr);
492 for (i = 0; i < numCells; ++i) { coneSize[i] = 2; }
493 for (i = 0; i < numVerts; ++i) { coneSize[numCells+i] = 0; }
494 for (i = 0; i < numCells; ++i) { cones[2*i] = numCells + i%numVerts; cones[2*i+1] = numCells + (i+1)%numVerts; }
495 for (i = 0; i < numVerts; ++i) { vertexCoords[i] = lower + (upper-lower)*((PetscReal)i/(PetscReal)numCells); }
496 ierr = DMPlexCreateFromDAG(*dm,1,numPoints,coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr);
497 ierr = PetscFree4(coneSize,cones,coneOrientations,vertexCoords);CHKERRQ(ierr);
498
499 ierr = PetscOptionsGetBool(((PetscObject)*dm)->options,((PetscObject)*dm)->prefix,"-dm_plex_separate_marker",&markerSeparate,NULL);CHKERRQ(ierr);
500 if (markerSeparate) { markerLeft = faceMarkerLeft; markerRight = faceMarkerRight;}
501 if (!wrap && !rank) {
502 ierr = DMPlexGetHeightStratum(*dm,1,&fStart,&fEnd);CHKERRQ(ierr);
503 ierr = DMSetLabelValue(*dm,"marker",fStart,markerLeft);CHKERRQ(ierr);
504 ierr = DMSetLabelValue(*dm,"marker",fEnd-1,markerRight);CHKERRQ(ierr);
505 ierr = DMSetLabelValue(*dm,"Face Sets",fStart,faceMarkerLeft);CHKERRQ(ierr);
506 ierr = DMSetLabelValue(*dm,"Face Sets",fEnd-1,faceMarkerRight);CHKERRQ(ierr);
507 }
508 if (wrap) {
509 L = upper - lower;
510 maxCell = (PetscReal)1.1*(L/(PetscReal)PetscMax(1,segments));
511 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,&maxCell,&L,&bd);CHKERRQ(ierr);
512 }
513 PetscFunctionReturn(0);
514 }
515
DMPlexCreateBoxMesh_Simplex_Internal(MPI_Comm comm,PetscInt dim,const PetscInt faces[],const PetscReal lower[],const PetscReal upper[],const DMBoundaryType periodicity[],PetscBool interpolate,DM * dm)516 static PetscErrorCode DMPlexCreateBoxMesh_Simplex_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
517 {
518 DM boundary;
519 PetscInt i;
520 PetscErrorCode ierr;
521
522 PetscFunctionBegin;
523 PetscValidPointer(dm, 4);
524 for (i = 0; i < dim; ++i) if (periodicity[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity is not supported for simplex meshes");
525 ierr = DMCreate(comm, &boundary);CHKERRQ(ierr);
526 PetscValidLogicalCollectiveInt(boundary,dim,2);
527 ierr = DMSetType(boundary, DMPLEX);CHKERRQ(ierr);
528 ierr = DMSetDimension(boundary, dim-1);CHKERRQ(ierr);
529 ierr = DMSetCoordinateDim(boundary, dim);CHKERRQ(ierr);
530 switch (dim) {
531 case 2: ierr = DMPlexCreateSquareBoundary(boundary, lower, upper, faces);CHKERRQ(ierr);break;
532 case 3: ierr = DMPlexCreateCubeBoundary(boundary, lower, upper, faces);CHKERRQ(ierr);break;
533 default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %D", dim);
534 }
535 ierr = DMPlexGenerate(boundary, NULL, interpolate, dm);CHKERRQ(ierr);
536 ierr = DMDestroy(&boundary);CHKERRQ(ierr);
537 PetscFunctionReturn(0);
538 }
539
DMPlexCreateCubeMesh_Internal(DM dm,const PetscReal lower[],const PetscReal upper[],const PetscInt edges[],DMBoundaryType bdX,DMBoundaryType bdY,DMBoundaryType bdZ)540 static PetscErrorCode DMPlexCreateCubeMesh_Internal(DM dm, const PetscReal lower[], const PetscReal upper[], const PetscInt edges[], DMBoundaryType bdX, DMBoundaryType bdY, DMBoundaryType bdZ)
541 {
542 DMLabel cutLabel = NULL;
543 PetscInt markerTop = 1, faceMarkerTop = 1;
544 PetscInt markerBottom = 1, faceMarkerBottom = 1;
545 PetscInt markerFront = 1, faceMarkerFront = 1;
546 PetscInt markerBack = 1, faceMarkerBack = 1;
547 PetscInt markerRight = 1, faceMarkerRight = 1;
548 PetscInt markerLeft = 1, faceMarkerLeft = 1;
549 PetscInt dim;
550 PetscBool markerSeparate = PETSC_FALSE, cutMarker = PETSC_FALSE;
551 PetscMPIInt rank;
552 PetscErrorCode ierr;
553
554 PetscFunctionBegin;
555 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
556 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);CHKERRQ(ierr);
557 ierr = DMCreateLabel(dm,"marker");CHKERRQ(ierr);
558 ierr = DMCreateLabel(dm,"Face Sets");CHKERRQ(ierr);
559 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_periodic_cut", &cutMarker, NULL);CHKERRQ(ierr);
560 if (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ||
561 bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ||
562 bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST) {
563
564 if (cutMarker) {ierr = DMCreateLabel(dm, "periodic_cut");CHKERRQ(ierr); ierr = DMGetLabel(dm, "periodic_cut", &cutLabel);CHKERRQ(ierr);}
565 }
566 switch (dim) {
567 case 2:
568 faceMarkerTop = 3;
569 faceMarkerBottom = 1;
570 faceMarkerRight = 2;
571 faceMarkerLeft = 4;
572 break;
573 case 3:
574 faceMarkerBottom = 1;
575 faceMarkerTop = 2;
576 faceMarkerFront = 3;
577 faceMarkerBack = 4;
578 faceMarkerRight = 5;
579 faceMarkerLeft = 6;
580 break;
581 default:
582 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dimension %D not supported",dim);
583 break;
584 }
585 ierr = PetscOptionsGetBool(((PetscObject) dm)->options,((PetscObject) dm)->prefix, "-dm_plex_separate_marker", &markerSeparate, NULL);CHKERRQ(ierr);
586 if (markerSeparate) {
587 markerBottom = faceMarkerBottom;
588 markerTop = faceMarkerTop;
589 markerFront = faceMarkerFront;
590 markerBack = faceMarkerBack;
591 markerRight = faceMarkerRight;
592 markerLeft = faceMarkerLeft;
593 }
594 {
595 const PetscInt numXEdges = !rank ? edges[0] : 0;
596 const PetscInt numYEdges = !rank ? edges[1] : 0;
597 const PetscInt numZEdges = !rank ? edges[2] : 0;
598 const PetscInt numXVertices = !rank ? (bdX == DM_BOUNDARY_PERIODIC || bdX == DM_BOUNDARY_TWIST ? edges[0] : edges[0]+1) : 0;
599 const PetscInt numYVertices = !rank ? (bdY == DM_BOUNDARY_PERIODIC || bdY == DM_BOUNDARY_TWIST ? edges[1] : edges[1]+1) : 0;
600 const PetscInt numZVertices = !rank ? (bdZ == DM_BOUNDARY_PERIODIC || bdZ == DM_BOUNDARY_TWIST ? edges[2] : edges[2]+1) : 0;
601 const PetscInt numCells = numXEdges*numYEdges*numZEdges;
602 const PetscInt numXFaces = numYEdges*numZEdges;
603 const PetscInt numYFaces = numXEdges*numZEdges;
604 const PetscInt numZFaces = numXEdges*numYEdges;
605 const PetscInt numTotXFaces = numXVertices*numXFaces;
606 const PetscInt numTotYFaces = numYVertices*numYFaces;
607 const PetscInt numTotZFaces = numZVertices*numZFaces;
608 const PetscInt numFaces = numTotXFaces + numTotYFaces + numTotZFaces;
609 const PetscInt numTotXEdges = numXEdges*numYVertices*numZVertices;
610 const PetscInt numTotYEdges = numYEdges*numXVertices*numZVertices;
611 const PetscInt numTotZEdges = numZEdges*numXVertices*numYVertices;
612 const PetscInt numVertices = numXVertices*numYVertices*numZVertices;
613 const PetscInt numEdges = numTotXEdges + numTotYEdges + numTotZEdges;
614 const PetscInt firstVertex = (dim == 2) ? numFaces : numCells;
615 const PetscInt firstXFace = (dim == 2) ? 0 : numCells + numVertices;
616 const PetscInt firstYFace = firstXFace + numTotXFaces;
617 const PetscInt firstZFace = firstYFace + numTotYFaces;
618 const PetscInt firstXEdge = numCells + numFaces + numVertices;
619 const PetscInt firstYEdge = firstXEdge + numTotXEdges;
620 const PetscInt firstZEdge = firstYEdge + numTotYEdges;
621 Vec coordinates;
622 PetscSection coordSection;
623 PetscScalar *coords;
624 PetscInt coordSize;
625 PetscInt v, vx, vy, vz;
626 PetscInt c, f, fx, fy, fz, e, ex, ey, ez;
627
628 ierr = DMPlexSetChart(dm, 0, numCells+numFaces+numEdges+numVertices);CHKERRQ(ierr);
629 for (c = 0; c < numCells; c++) {
630 ierr = DMPlexSetConeSize(dm, c, 6);CHKERRQ(ierr);
631 }
632 for (f = firstXFace; f < firstXFace+numFaces; ++f) {
633 ierr = DMPlexSetConeSize(dm, f, 4);CHKERRQ(ierr);
634 }
635 for (e = firstXEdge; e < firstXEdge+numEdges; ++e) {
636 ierr = DMPlexSetConeSize(dm, e, 2);CHKERRQ(ierr);
637 }
638 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
639 /* Build cells */
640 for (fz = 0; fz < numZEdges; ++fz) {
641 for (fy = 0; fy < numYEdges; ++fy) {
642 for (fx = 0; fx < numXEdges; ++fx) {
643 PetscInt cell = (fz*numYEdges + fy)*numXEdges + fx;
644 PetscInt faceB = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
645 PetscInt faceT = firstZFace + (fy*numXEdges+fx)*numZVertices + ((fz+1)%numZVertices);
646 PetscInt faceF = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
647 PetscInt faceK = firstYFace + (fz*numXEdges+fx)*numYVertices + ((fy+1)%numYVertices);
648 PetscInt faceL = firstXFace + (fz*numYEdges+fy)*numXVertices + fx;
649 PetscInt faceR = firstXFace + (fz*numYEdges+fy)*numXVertices + ((fx+1)%numXVertices);
650 /* B, T, F, K, R, L */
651 PetscInt ornt[6] = {-4, 0, 0, -1, 0, -4}; /* ??? */
652 PetscInt cone[6];
653
654 /* no boundary twisting in 3D */
655 cone[0] = faceB; cone[1] = faceT; cone[2] = faceF; cone[3] = faceK; cone[4] = faceR; cone[5] = faceL;
656 ierr = DMPlexSetCone(dm, cell, cone);CHKERRQ(ierr);
657 ierr = DMPlexSetConeOrientation(dm, cell, ornt);CHKERRQ(ierr);
658 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);}
659 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);}
660 if (bdZ != DM_BOUNDARY_NONE && fz == numZEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, cell, 2);CHKERRQ(ierr);}
661 }
662 }
663 }
664 /* Build x faces */
665 for (fz = 0; fz < numZEdges; ++fz) {
666 for (fy = 0; fy < numYEdges; ++fy) {
667 for (fx = 0; fx < numXVertices; ++fx) {
668 PetscInt face = firstXFace + (fz*numYEdges+fy) *numXVertices+fx;
669 PetscInt edgeL = firstZEdge + (fy *numXVertices+fx)*numZEdges + fz;
670 PetscInt edgeR = firstZEdge + (((fy+1)%numYVertices)*numXVertices+fx)*numZEdges + fz;
671 PetscInt edgeB = firstYEdge + (fz *numXVertices+fx)*numYEdges + fy;
672 PetscInt edgeT = firstYEdge + (((fz+1)%numZVertices)*numXVertices+fx)*numYEdges + fy;
673 PetscInt ornt[4] = {0, 0, -2, -2};
674 PetscInt cone[4];
675
676 if (dim == 3) {
677 /* markers */
678 if (bdX != DM_BOUNDARY_PERIODIC) {
679 if (fx == numXVertices-1) {
680 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerRight);CHKERRQ(ierr);
681 ierr = DMSetLabelValue(dm, "marker", face, markerRight);CHKERRQ(ierr);
682 }
683 else if (fx == 0) {
684 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerLeft);CHKERRQ(ierr);
685 ierr = DMSetLabelValue(dm, "marker", face, markerLeft);CHKERRQ(ierr);
686 }
687 }
688 }
689 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
690 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr);
691 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr);
692 }
693 }
694 }
695 /* Build y faces */
696 for (fz = 0; fz < numZEdges; ++fz) {
697 for (fx = 0; fx < numXEdges; ++fx) {
698 for (fy = 0; fy < numYVertices; ++fy) {
699 PetscInt face = firstYFace + (fz*numXEdges+fx)*numYVertices + fy;
700 PetscInt edgeL = firstZEdge + (fy*numXVertices+ fx)*numZEdges + fz;
701 PetscInt edgeR = firstZEdge + (fy*numXVertices+((fx+1)%numXVertices))*numZEdges + fz;
702 PetscInt edgeB = firstXEdge + (fz *numYVertices+fy)*numXEdges + fx;
703 PetscInt edgeT = firstXEdge + (((fz+1)%numZVertices)*numYVertices+fy)*numXEdges + fx;
704 PetscInt ornt[4] = {0, 0, -2, -2};
705 PetscInt cone[4];
706
707 if (dim == 3) {
708 /* markers */
709 if (bdY != DM_BOUNDARY_PERIODIC) {
710 if (fy == numYVertices-1) {
711 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBack);CHKERRQ(ierr);
712 ierr = DMSetLabelValue(dm, "marker", face, markerBack);CHKERRQ(ierr);
713 }
714 else if (fy == 0) {
715 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerFront);CHKERRQ(ierr);
716 ierr = DMSetLabelValue(dm, "marker", face, markerFront);CHKERRQ(ierr);
717 }
718 }
719 }
720 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
721 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr);
722 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr);
723 }
724 }
725 }
726 /* Build z faces */
727 for (fy = 0; fy < numYEdges; ++fy) {
728 for (fx = 0; fx < numXEdges; ++fx) {
729 for (fz = 0; fz < numZVertices; fz++) {
730 PetscInt face = firstZFace + (fy*numXEdges+fx)*numZVertices + fz;
731 PetscInt edgeL = firstYEdge + (fz*numXVertices+ fx)*numYEdges + fy;
732 PetscInt edgeR = firstYEdge + (fz*numXVertices+((fx+1)%numXVertices))*numYEdges + fy;
733 PetscInt edgeB = firstXEdge + (fz*numYVertices+ fy)*numXEdges + fx;
734 PetscInt edgeT = firstXEdge + (fz*numYVertices+((fy+1)%numYVertices))*numXEdges + fx;
735 PetscInt ornt[4] = {0, 0, -2, -2};
736 PetscInt cone[4];
737
738 if (dim == 2) {
739 if (bdX == DM_BOUNDARY_TWIST && fx == numXEdges-1) {edgeR += numYEdges-1-2*fy; ornt[1] = -2;}
740 if (bdY == DM_BOUNDARY_TWIST && fy == numYEdges-1) {edgeT += numXEdges-1-2*fx; ornt[2] = 0;}
741 if (bdX != DM_BOUNDARY_NONE && fx == numXEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);}
742 if (bdY != DM_BOUNDARY_NONE && fy == numYEdges-1 && cutLabel) {ierr = DMLabelSetValue(cutLabel, face, 2);CHKERRQ(ierr);}
743 } else {
744 /* markers */
745 if (bdZ != DM_BOUNDARY_PERIODIC) {
746 if (fz == numZVertices-1) {
747 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerTop);CHKERRQ(ierr);
748 ierr = DMSetLabelValue(dm, "marker", face, markerTop);CHKERRQ(ierr);
749 }
750 else if (fz == 0) {
751 ierr = DMSetLabelValue(dm, "Face Sets", face, faceMarkerBottom);CHKERRQ(ierr);
752 ierr = DMSetLabelValue(dm, "marker", face, markerBottom);CHKERRQ(ierr);
753 }
754 }
755 }
756 cone[0] = edgeB; cone[1] = edgeR; cone[2] = edgeT; cone[3] = edgeL;
757 ierr = DMPlexSetCone(dm, face, cone);CHKERRQ(ierr);
758 ierr = DMPlexSetConeOrientation(dm, face, ornt);CHKERRQ(ierr);
759 }
760 }
761 }
762 /* Build Z edges*/
763 for (vy = 0; vy < numYVertices; vy++) {
764 for (vx = 0; vx < numXVertices; vx++) {
765 for (ez = 0; ez < numZEdges; ez++) {
766 const PetscInt edge = firstZEdge + (vy*numXVertices+vx)*numZEdges + ez;
767 const PetscInt vertexB = firstVertex + (ez *numYVertices+vy)*numXVertices + vx;
768 const PetscInt vertexT = firstVertex + (((ez+1)%numZVertices)*numYVertices+vy)*numXVertices + vx;
769 PetscInt cone[2];
770
771 if (dim == 3) {
772 if (bdX != DM_BOUNDARY_PERIODIC) {
773 if (vx == numXVertices-1) {
774 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr);
775 }
776 else if (vx == 0) {
777 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr);
778 }
779 }
780 if (bdY != DM_BOUNDARY_PERIODIC) {
781 if (vy == numYVertices-1) {
782 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr);
783 }
784 else if (vy == 0) {
785 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr);
786 }
787 }
788 }
789 cone[0] = vertexB; cone[1] = vertexT;
790 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
791 }
792 }
793 }
794 /* Build Y edges*/
795 for (vz = 0; vz < numZVertices; vz++) {
796 for (vx = 0; vx < numXVertices; vx++) {
797 for (ey = 0; ey < numYEdges; ey++) {
798 const PetscInt nextv = (dim == 2 && bdY == DM_BOUNDARY_TWIST && ey == numYEdges-1) ? (numXVertices-vx-1) : (vz*numYVertices+((ey+1)%numYVertices))*numXVertices + vx;
799 const PetscInt edge = firstYEdge + (vz*numXVertices+vx)*numYEdges + ey;
800 const PetscInt vertexF = firstVertex + (vz*numYVertices+ey)*numXVertices + vx;
801 const PetscInt vertexK = firstVertex + nextv;
802 PetscInt cone[2];
803
804 cone[0] = vertexF; cone[1] = vertexK;
805 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
806 if (dim == 2) {
807 if ((bdX != DM_BOUNDARY_PERIODIC) && (bdX != DM_BOUNDARY_TWIST)) {
808 if (vx == numXVertices-1) {
809 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerRight);CHKERRQ(ierr);
810 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr);
811 ierr = DMSetLabelValue(dm, "marker", cone[0], markerRight);CHKERRQ(ierr);
812 if (ey == numYEdges-1) {
813 ierr = DMSetLabelValue(dm, "marker", cone[1], markerRight);CHKERRQ(ierr);
814 }
815 } else if (vx == 0) {
816 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerLeft);CHKERRQ(ierr);
817 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr);
818 ierr = DMSetLabelValue(dm, "marker", cone[0], markerLeft);CHKERRQ(ierr);
819 if (ey == numYEdges-1) {
820 ierr = DMSetLabelValue(dm, "marker", cone[1], markerLeft);CHKERRQ(ierr);
821 }
822 }
823 } else {
824 if (vx == 0 && cutLabel) {
825 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr);
826 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr);
827 if (ey == numYEdges-1) {
828 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr);
829 }
830 }
831 }
832 } else {
833 if (bdX != DM_BOUNDARY_PERIODIC) {
834 if (vx == numXVertices-1) {
835 ierr = DMSetLabelValue(dm, "marker", edge, markerRight);CHKERRQ(ierr);
836 } else if (vx == 0) {
837 ierr = DMSetLabelValue(dm, "marker", edge, markerLeft);CHKERRQ(ierr);
838 }
839 }
840 if (bdZ != DM_BOUNDARY_PERIODIC) {
841 if (vz == numZVertices-1) {
842 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr);
843 } else if (vz == 0) {
844 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr);
845 }
846 }
847 }
848 }
849 }
850 }
851 /* Build X edges*/
852 for (vz = 0; vz < numZVertices; vz++) {
853 for (vy = 0; vy < numYVertices; vy++) {
854 for (ex = 0; ex < numXEdges; ex++) {
855 const PetscInt nextv = (dim == 2 && bdX == DM_BOUNDARY_TWIST && ex == numXEdges-1) ? (numYVertices-vy-1)*numXVertices : (vz*numYVertices+vy)*numXVertices + (ex+1)%numXVertices;
856 const PetscInt edge = firstXEdge + (vz*numYVertices+vy)*numXEdges + ex;
857 const PetscInt vertexL = firstVertex + (vz*numYVertices+vy)*numXVertices + ex;
858 const PetscInt vertexR = firstVertex + nextv;
859 PetscInt cone[2];
860
861 cone[0] = vertexL; cone[1] = vertexR;
862 ierr = DMPlexSetCone(dm, edge, cone);CHKERRQ(ierr);
863 if (dim == 2) {
864 if ((bdY != DM_BOUNDARY_PERIODIC) && (bdY != DM_BOUNDARY_TWIST)) {
865 if (vy == numYVertices-1) {
866 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerTop);CHKERRQ(ierr);
867 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr);
868 ierr = DMSetLabelValue(dm, "marker", cone[0], markerTop);CHKERRQ(ierr);
869 if (ex == numXEdges-1) {
870 ierr = DMSetLabelValue(dm, "marker", cone[1], markerTop);CHKERRQ(ierr);
871 }
872 } else if (vy == 0) {
873 ierr = DMSetLabelValue(dm, "Face Sets", edge, faceMarkerBottom);CHKERRQ(ierr);
874 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr);
875 ierr = DMSetLabelValue(dm, "marker", cone[0], markerBottom);CHKERRQ(ierr);
876 if (ex == numXEdges-1) {
877 ierr = DMSetLabelValue(dm, "marker", cone[1], markerBottom);CHKERRQ(ierr);
878 }
879 }
880 } else {
881 if (vy == 0 && cutLabel) {
882 ierr = DMLabelSetValue(cutLabel, edge, 1);CHKERRQ(ierr);
883 ierr = DMLabelSetValue(cutLabel, cone[0], 1);CHKERRQ(ierr);
884 if (ex == numXEdges-1) {
885 ierr = DMLabelSetValue(cutLabel, cone[1], 1);CHKERRQ(ierr);
886 }
887 }
888 }
889 } else {
890 if (bdY != DM_BOUNDARY_PERIODIC) {
891 if (vy == numYVertices-1) {
892 ierr = DMSetLabelValue(dm, "marker", edge, markerBack);CHKERRQ(ierr);
893 }
894 else if (vy == 0) {
895 ierr = DMSetLabelValue(dm, "marker", edge, markerFront);CHKERRQ(ierr);
896 }
897 }
898 if (bdZ != DM_BOUNDARY_PERIODIC) {
899 if (vz == numZVertices-1) {
900 ierr = DMSetLabelValue(dm, "marker", edge, markerTop);CHKERRQ(ierr);
901 }
902 else if (vz == 0) {
903 ierr = DMSetLabelValue(dm, "marker", edge, markerBottom);CHKERRQ(ierr);
904 }
905 }
906 }
907 }
908 }
909 }
910 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
911 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
912 /* Build coordinates */
913 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
914 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
915 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
916 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVertices);CHKERRQ(ierr);
917 for (v = firstVertex; v < firstVertex+numVertices; ++v) {
918 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
919 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
920 }
921 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
922 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
923 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
924 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
925 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
926 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
927 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
928 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
929 for (vz = 0; vz < numZVertices; ++vz) {
930 for (vy = 0; vy < numYVertices; ++vy) {
931 for (vx = 0; vx < numXVertices; ++vx) {
932 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+0] = lower[0] + ((upper[0] - lower[0])/numXEdges)*vx;
933 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+1] = lower[1] + ((upper[1] - lower[1])/numYEdges)*vy;
934 if (dim == 3) {
935 coords[((vz*numYVertices+vy)*numXVertices+vx)*dim+2] = lower[2] + ((upper[2] - lower[2])/numZEdges)*vz;
936 }
937 }
938 }
939 }
940 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
941 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
942 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
943 }
944 PetscFunctionReturn(0);
945 }
946
DMPlexCreateBoxMesh_Tensor_Internal(MPI_Comm comm,PetscInt dim,const PetscInt faces[],const PetscReal lower[],const PetscReal upper[],const DMBoundaryType periodicity[],PetscBool interpolate,DM * dm)947 static PetscErrorCode DMPlexCreateBoxMesh_Tensor_Internal(MPI_Comm comm, PetscInt dim, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
948 {
949 PetscInt i;
950 PetscErrorCode ierr;
951
952 PetscFunctionBegin;
953 PetscValidPointer(dm, 7);
954 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
955 PetscValidLogicalCollectiveInt(*dm,dim,2);
956 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
957 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
958 switch (dim) {
959 case 2: {ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], DM_BOUNDARY_NONE);CHKERRQ(ierr);break;}
960 case 3: {ierr = DMPlexCreateCubeMesh_Internal(*dm, lower, upper, faces, periodicity[0], periodicity[1], periodicity[2]);CHKERRQ(ierr);break;}
961 default: SETERRQ1(comm, PETSC_ERR_SUP, "Dimension not supported: %D", dim);
962 }
963 if (periodicity[0] == DM_BOUNDARY_PERIODIC || periodicity[0] == DM_BOUNDARY_TWIST ||
964 periodicity[1] == DM_BOUNDARY_PERIODIC || periodicity[1] == DM_BOUNDARY_TWIST ||
965 (dim > 2 && (periodicity[2] == DM_BOUNDARY_PERIODIC || periodicity[2] == DM_BOUNDARY_TWIST))) {
966 PetscReal L[3];
967 PetscReal maxCell[3];
968
969 for (i = 0; i < dim; i++) {
970 L[i] = upper[i] - lower[i];
971 maxCell[i] = 1.1 * (L[i] / PetscMax(1,faces[i]));
972 }
973 ierr = DMSetPeriodicity(*dm,PETSC_TRUE,maxCell,L,periodicity);CHKERRQ(ierr);
974 }
975 if (!interpolate) {
976 DM udm;
977
978 ierr = DMPlexUninterpolate(*dm, &udm);CHKERRQ(ierr);
979 ierr = DMPlexCopyCoordinates(*dm, udm);CHKERRQ(ierr);
980 ierr = DMDestroy(dm);CHKERRQ(ierr);
981 *dm = udm;
982 }
983 PetscFunctionReturn(0);
984 }
985
986 /*@C
987 DMPlexCreateBoxMesh - Creates a mesh on the tensor product of unit intervals (box) using simplices or tensor cells (hexahedra).
988
989 Collective
990
991 Input Parameters:
992 + comm - The communicator for the DM object
993 . dim - The spatial dimension
994 . simplex - PETSC_TRUE for simplices, PETSC_FALSE for tensor cells
995 . faces - Number of faces per dimension, or NULL for (1,) in 1D and (2, 2) in 2D and (1, 1, 1) in 3D
996 . lower - The lower left corner, or NULL for (0, 0, 0)
997 . upper - The upper right corner, or NULL for (1, 1, 1)
998 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
999 - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1000
1001 Output Parameter:
1002 . dm - The DM object
1003
1004 Options Database Keys:
1005 These options override the hard-wired input values.
1006 + -dm_plex_box_dim <dim> - Set the topological dimension
1007 . -dm_plex_box_simplex <bool> - PETSC_TRUE for simplex elements, PETSC_FALSE for tensor elements
1008 . -dm_plex_box_lower <x,y,z> - Specify lower-left-bottom coordinates for the box
1009 . -dm_plex_box_upper <x,y,z> - Specify upper-right-top coordinates for the box
1010 . -dm_plex_box_faces <m,n,p> - Number of faces in each linear direction
1011 . -dm_plex_box_bd <bx,by,bz> - Specify the DMBoundaryType for each direction
1012 - -dm_plex_box_interpolate <bool> - PETSC_TRUE turns on topological interpolation (creating edges and faces)
1013
1014 Notes:
1015 The options database keys above take lists of length d in d dimensions.
1016
1017 Here is the numbering returned for 2 faces in each direction for tensor cells:
1018 $ 10---17---11---18----12
1019 $ | | |
1020 $ | | |
1021 $ 20 2 22 3 24
1022 $ | | |
1023 $ | | |
1024 $ 7---15----8---16----9
1025 $ | | |
1026 $ | | |
1027 $ 19 0 21 1 23
1028 $ | | |
1029 $ | | |
1030 $ 4---13----5---14----6
1031
1032 and for simplicial cells
1033
1034 $ 14----8---15----9----16
1035 $ |\ 5 |\ 7 |
1036 $ | \ | \ |
1037 $ 13 2 14 3 15
1038 $ | 4 \ | 6 \ |
1039 $ | \ | \ |
1040 $ 11----6---12----7----13
1041 $ |\ |\ |
1042 $ | \ 1 | \ 3 |
1043 $ 10 0 11 1 12
1044 $ | 0 \ | 2 \ |
1045 $ | \ | \ |
1046 $ 8----4----9----5----10
1047
1048 Level: beginner
1049
1050 .seealso: DMPlexCreateFromFile(), DMPlexCreateHexCylinderMesh(), DMSetType(), DMCreate()
1051 @*/
DMPlexCreateBoxMesh(MPI_Comm comm,PetscInt dim,PetscBool simplex,const PetscInt faces[],const PetscReal lower[],const PetscReal upper[],const DMBoundaryType periodicity[],PetscBool interpolate,DM * dm)1052 PetscErrorCode DMPlexCreateBoxMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool interpolate, DM *dm)
1053 {
1054 PetscInt fac[3] = {0, 0, 0};
1055 PetscReal low[3] = {0, 0, 0};
1056 PetscReal upp[3] = {1, 1, 1};
1057 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1058 PetscInt i, n;
1059 PetscBool flg;
1060 PetscErrorCode ierr;
1061
1062 PetscFunctionBegin;
1063 ierr = PetscOptionsGetInt(NULL, NULL, "-dm_plex_box_dim", &dim, &flg);CHKERRQ(ierr);
1064 if ((dim < 0) || (dim > 3)) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Dimension %D should be in [1, 3]", dim);
1065 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_box_simplex", &simplex, &flg);CHKERRQ(ierr);
1066 n = dim;
1067 ierr = PetscOptionsGetIntArray(NULL, NULL, "-dm_plex_box_faces", fac, &n, &flg);CHKERRQ(ierr);
1068 for (i = 0; i < dim; ++i) fac[i] = faces ? faces[i] : (flg && i < n ? fac[i] : (dim == 1 ? 1 : 4-dim));
1069 if (lower) for (i = 0; i < dim; ++i) low[i] = lower[i];
1070 if (upper) for (i = 0; i < dim; ++i) upp[i] = upper[i];
1071 if (periodicity) for (i = 0; i < dim; ++i) bdt[i] = periodicity[i];
1072 /* Allow bounds to be specified from the command line */
1073 n = 3;
1074 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_lower", low, &n, &flg);CHKERRQ(ierr);
1075 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Lower box point had %D values, should have been %D", n, dim);
1076 n = 3;
1077 ierr = PetscOptionsGetRealArray(NULL, NULL, "-dm_plex_box_upper", upp, &n, &flg);CHKERRQ(ierr);
1078 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Upper box point had %D values, should have been %D", n, dim);
1079 n = 3;
1080 ierr = PetscOptionsGetEnumArray(NULL, NULL, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *) bdt, &n, &flg);CHKERRQ(ierr);
1081 if (flg && (n != dim)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Box boundary types had %D values, should have been %D", n, dim);
1082 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_box_interpolate", &interpolate, &flg);CHKERRQ(ierr);
1083
1084 if (dim == 1) {ierr = DMPlexCreateLineMesh_Internal(comm, fac[0], low[0], upp[0], bdt[0], dm);CHKERRQ(ierr);}
1085 else if (simplex) {ierr = DMPlexCreateBoxMesh_Simplex_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);}
1086 else {ierr = DMPlexCreateBoxMesh_Tensor_Internal(comm, dim, fac, low, upp, bdt, interpolate, dm);CHKERRQ(ierr);}
1087 PetscFunctionReturn(0);
1088 }
1089
1090 /*@
1091 DMPlexCreateWedgeBoxMesh - Creates a 3-D mesh tesselating the (x,y) plane and extruding in the third direction using wedge cells.
1092
1093 Collective
1094
1095 Input Parameters:
1096 + comm - The communicator for the DM object
1097 . faces - Number of faces per dimension, or NULL for (1, 1, 1)
1098 . lower - The lower left corner, or NULL for (0, 0, 0)
1099 . upper - The upper right corner, or NULL for (1, 1, 1)
1100 . periodicity - The boundary type for the X,Y,Z direction, or NULL for DM_BOUNDARY_NONE
1101 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1102 - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1103
1104 Output Parameter:
1105 . dm - The DM object
1106
1107 Level: beginner
1108
1109 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateWedgeCylinderMesh(), DMPlexExtrude(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1110 @*/
DMPlexCreateWedgeBoxMesh(MPI_Comm comm,const PetscInt faces[],const PetscReal lower[],const PetscReal upper[],const DMBoundaryType periodicity[],PetscBool orderHeight,PetscBool interpolate,DM * dm)1111 PetscErrorCode DMPlexCreateWedgeBoxMesh(MPI_Comm comm, const PetscInt faces[], const PetscReal lower[], const PetscReal upper[], const DMBoundaryType periodicity[], PetscBool orderHeight, PetscBool interpolate, DM *dm)
1112 {
1113 DM bdm, botdm;
1114 PetscInt i;
1115 PetscInt fac[3] = {0, 0, 0};
1116 PetscReal low[3] = {0, 0, 0};
1117 PetscReal upp[3] = {1, 1, 1};
1118 DMBoundaryType bdt[3] = {DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE};
1119 PetscErrorCode ierr;
1120
1121 PetscFunctionBegin;
1122 for (i = 0; i < 3; ++i) fac[i] = faces ? (faces[i] > 0 ? faces[i] : 1) : 1;
1123 if (lower) for (i = 0; i < 3; ++i) low[i] = lower[i];
1124 if (upper) for (i = 0; i < 3; ++i) upp[i] = upper[i];
1125 if (periodicity) for (i = 0; i < 3; ++i) bdt[i] = periodicity[i];
1126 for (i = 0; i < 3; ++i) if (bdt[i] != DM_BOUNDARY_NONE) SETERRQ(comm, PETSC_ERR_SUP, "Periodicity not yet supported");
1127
1128 ierr = DMCreate(comm, &bdm);CHKERRQ(ierr);
1129 ierr = DMSetType(bdm, DMPLEX);CHKERRQ(ierr);
1130 ierr = DMSetDimension(bdm, 1);CHKERRQ(ierr);
1131 ierr = DMSetCoordinateDim(bdm, 2);CHKERRQ(ierr);
1132 ierr = DMPlexCreateSquareBoundary(bdm, low, upp, fac);CHKERRQ(ierr);
1133 ierr = DMPlexGenerate(bdm, NULL, PETSC_FALSE, &botdm);CHKERRQ(ierr);
1134 ierr = DMDestroy(&bdm);CHKERRQ(ierr);
1135 ierr = DMPlexExtrude(botdm, fac[2], upp[2] - low[2], orderHeight, NULL, interpolate, dm);CHKERRQ(ierr);
1136 if (low[2] != 0.0) {
1137 Vec v;
1138 PetscScalar *x;
1139 PetscInt cDim, n;
1140
1141 ierr = DMGetCoordinatesLocal(*dm, &v);CHKERRQ(ierr);
1142 ierr = VecGetBlockSize(v, &cDim);CHKERRQ(ierr);
1143 ierr = VecGetLocalSize(v, &n);CHKERRQ(ierr);
1144 ierr = VecGetArray(v, &x);CHKERRQ(ierr);
1145 x += cDim;
1146 for (i=0; i<n; i+=cDim) x[i] += low[2];
1147 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr);
1148 ierr = DMSetCoordinatesLocal(*dm, v);CHKERRQ(ierr);
1149 }
1150 ierr = DMDestroy(&botdm);CHKERRQ(ierr);
1151 PetscFunctionReturn(0);
1152 }
1153
1154 /*@C
1155 DMPlexExtrude - Creates a (d+1)-D mesh by extruding a d-D mesh in the normal direction using prismatic cells.
1156
1157 Collective on idm
1158
1159 Input Parameters:
1160 + idm - The mesh to be extruded
1161 . layers - The number of layers, or PETSC_DETERMINE to use the default
1162 . height - The total height of the extrusion, or PETSC_DETERMINE to use the default
1163 . orderHeight - If PETSC_TRUE, orders the extruded cells in the height first. Otherwise, orders the cell on the layers first
1164 . extNormal - The normal direction in which the mesh should be extruded, or NULL to extrude using the surface normal
1165 - interpolate - Flag to create intermediate mesh pieces (edges, faces)
1166
1167 Output Parameter:
1168 . dm - The DM object
1169
1170 Notes:
1171 The mesh created has prismatic cells, and the vertex ordering in the cone of the cell is that of the tensor prismatic cells. Not currently supported in Fortran.
1172
1173 Options Database Keys:
1174 + -dm_plex_extrude_layers <k> - Sets the number of layers k
1175 . -dm_plex_extrude_height <h> - Sets the total height of the extrusion
1176 . -dm_plex_extrude_heights <h0,h1,...> - Sets the height of each layer
1177 . -dm_plex_extrude_order_height - If true, order cells by height first
1178 - -dm_plex_extrude_normal <n0,...,nd> - Sets the normal vector along which to extrude
1179
1180 Level: advanced
1181
1182 .seealso: DMPlexCreateWedgeCylinderMesh(), DMPlexCreateWedgeBoxMesh(), DMSetType(), DMCreate()
1183 @*/
DMPlexExtrude(DM idm,PetscInt layers,PetscReal height,PetscBool orderHeight,const PetscReal extNormal[],PetscBool interpolate,DM * dm)1184 PetscErrorCode DMPlexExtrude(DM idm, PetscInt layers, PetscReal height, PetscBool orderHeight, const PetscReal extNormal[], PetscBool interpolate, DM* dm)
1185 {
1186 PetscScalar *coordsB;
1187 const PetscScalar *coordsA;
1188 PetscReal *normals = NULL, *heights = NULL;
1189 PetscReal clNormal[3];
1190 Vec coordinatesA, coordinatesB;
1191 PetscSection coordSectionA, coordSectionB;
1192 PetscInt dim, cDim, cDimB, c, l, v, coordSize, *newCone, nl;
1193 PetscInt cStart, cEnd, vStart, vEnd, cellV, numCells, numVertices;
1194 const char *prefix;
1195 PetscBool haveCLNormal, flg;
1196 PetscErrorCode ierr;
1197
1198 PetscFunctionBegin;
1199 PetscValidHeaderSpecific(idm, DM_CLASSID, 1);
1200 PetscValidLogicalCollectiveInt(idm, layers, 2);
1201 PetscValidLogicalCollectiveReal(idm, height, 3);
1202 PetscValidLogicalCollectiveBool(idm, interpolate, 4);
1203 ierr = DMGetDimension(idm, &dim);CHKERRQ(ierr);
1204 ierr = DMGetCoordinateDim(idm, &cDim);CHKERRQ(ierr);
1205 cDimB = cDim == dim ? cDim+1 : cDim;
1206 if (dim < 1 || dim > 3) SETERRQ1(PetscObjectComm((PetscObject)idm), PETSC_ERR_SUP, "Support for dimension %D not coded", dim);
1207
1208 ierr = PetscObjectGetOptionsPrefix((PetscObject) idm, &prefix);CHKERRQ(ierr);
1209 if (layers < 0) layers = 1;
1210 ierr = PetscOptionsGetInt(NULL, prefix, "-dm_plex_extrude_layers", &layers, NULL);CHKERRQ(ierr);
1211 if (layers <= 0) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Number of layers %D must be positive", layers);
1212 if (height < 0.) height = 1.;
1213 ierr = PetscOptionsGetReal(NULL, prefix, "-dm_plex_extrude_height", &height, NULL);CHKERRQ(ierr);
1214 if (height <= 0.) SETERRQ1(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height of layers %g must be positive", (double) height);
1215 ierr = PetscMalloc1(layers, &heights);CHKERRQ(ierr);
1216 nl = layers;
1217 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_heights", heights, &nl, &flg);CHKERRQ(ierr);
1218 if (flg) {
1219 if (!nl) SETERRQ(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Must give at least one height for -dm_plex_extrude_heights");
1220 for (l = nl; l < layers; ++l) heights[l] = heights[l-1];
1221 for (l = 0; l < layers; ++l) if (heights[l] <= 0.) SETERRQ2(PetscObjectComm((PetscObject) idm), PETSC_ERR_ARG_OUTOFRANGE, "Height %g of layers %D must be positive", (double) heights[l], l);
1222 } else {
1223 for (l = 0; l < layers; ++l) heights[l] = height/layers;
1224 }
1225 ierr = PetscOptionsGetBool(NULL, prefix, "-dm_plex_extrude_order_height", &orderHeight, NULL);CHKERRQ(ierr);
1226 c = 3;
1227 ierr = PetscOptionsGetRealArray(NULL, prefix, "-dm_plex_extrude_normal", clNormal, &c, &haveCLNormal);CHKERRQ(ierr);
1228 if (haveCLNormal && c != cDimB) SETERRQ2(PetscObjectComm((PetscObject)idm), PETSC_ERR_ARG_SIZ, "Input normal has size %D != %D extruded coordinate dimension", c, cDimB);
1229
1230 ierr = DMPlexGetHeightStratum(idm, 0, &cStart, &cEnd);CHKERRQ(ierr);
1231 ierr = DMPlexGetDepthStratum(idm, 0, &vStart, &vEnd);CHKERRQ(ierr);
1232 numCells = (cEnd - cStart)*layers;
1233 numVertices = (vEnd - vStart)*(layers+1);
1234 ierr = DMCreate(PetscObjectComm((PetscObject)idm), dm);CHKERRQ(ierr);
1235 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
1236 ierr = DMSetDimension(*dm, dim+1);CHKERRQ(ierr);
1237 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr);
1238 /* Must create the celltype label here so that we do not automatically try to compute the types */
1239 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr);
1240 for (c = cStart, cellV = 0; c < cEnd; ++c) {
1241 DMPolytopeType ct, nct;
1242 PetscInt *closure = NULL;
1243 PetscInt closureSize, numCorners = 0;
1244
1245 ierr = DMPlexGetCellType(idm, c, &ct);CHKERRQ(ierr);
1246 switch (ct) {
1247 case DM_POLYTOPE_SEGMENT: nct = DM_POLYTOPE_SEG_PRISM_TENSOR;break;
1248 case DM_POLYTOPE_TRIANGLE: nct = DM_POLYTOPE_TRI_PRISM_TENSOR;break;
1249 case DM_POLYTOPE_QUADRILATERAL: nct = DM_POLYTOPE_QUAD_PRISM_TENSOR;break;
1250 default: nct = DM_POLYTOPE_UNKNOWN;
1251 }
1252 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1253 for (v = 0; v < closureSize*2; v += 2) if ((closure[v] >= vStart) && (closure[v] < vEnd)) numCorners++;
1254 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1255 for (l = 0; l < layers; ++l) {
1256 const PetscInt cell = orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart;
1257
1258 ierr = DMPlexSetConeSize(*dm, cell, 2*numCorners);CHKERRQ(ierr);
1259 ierr = DMPlexSetCellType(*dm, cell, nct);CHKERRQ(ierr);
1260 }
1261 cellV = PetscMax(numCorners,cellV);
1262 }
1263 ierr = DMSetUp(*dm);CHKERRQ(ierr);
1264
1265 if (dim != cDim && !(extNormal || haveCLNormal)) {ierr = PetscCalloc1(cDim*(vEnd - vStart), &normals);CHKERRQ(ierr);}
1266 ierr = PetscMalloc1(3*cellV,&newCone);CHKERRQ(ierr);
1267 for (c = cStart; c < cEnd; ++c) {
1268 PetscInt *closure = NULL;
1269 PetscInt closureSize, numCorners = 0, l;
1270 PetscReal normal[3] = {0, 0, 0};
1271
1272 if (normals) {ierr = DMPlexComputeCellGeometryFVM(idm, c, NULL, NULL, normal);CHKERRQ(ierr);}
1273 ierr = DMPlexGetTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1274 for (v = 0; v < closureSize*2; v += 2) {
1275 if ((closure[v] >= vStart) && (closure[v] < vEnd)) {
1276 PetscInt d;
1277
1278 newCone[numCorners++] = closure[v] - vStart;
1279 if (normals) {for (d = 0; d < cDim; ++d) normals[cDim*(closure[v]-vStart)+d] += normal[d];}
1280 }
1281 }
1282 ierr = DMPlexRestoreTransitiveClosure(idm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr);
1283 for (l = 0; l < layers; ++l) {
1284 PetscInt i;
1285
1286 for (i = 0; i < numCorners; ++i) {
1287 newCone[ numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + numCells : l*(vEnd - vStart) + newCone[i] + numCells;
1288 newCone[2*numCorners + i] = orderHeight ? (layers+1)*newCone[i] + l + 1 + numCells : (l+1)*(vEnd - vStart) + newCone[i] + numCells;
1289 }
1290 ierr = DMPlexSetCone(*dm, orderHeight ? layers*(c - cStart) + l : l*(cEnd - cStart) + c - cStart, newCone + numCorners);CHKERRQ(ierr);
1291 }
1292 }
1293 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
1294 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
1295 ierr = PetscFree(newCone);CHKERRQ(ierr);
1296
1297 ierr = DMGetCoordinateSection(*dm, &coordSectionB);CHKERRQ(ierr);
1298 ierr = PetscSectionSetNumFields(coordSectionB, 1);CHKERRQ(ierr);
1299 ierr = PetscSectionSetFieldComponents(coordSectionB, 0, cDimB);CHKERRQ(ierr);
1300 ierr = PetscSectionSetChart(coordSectionB, numCells, numCells+numVertices);CHKERRQ(ierr);
1301 for (v = numCells; v < numCells+numVertices; ++v) {
1302 ierr = PetscSectionSetDof(coordSectionB, v, cDimB);CHKERRQ(ierr);
1303 ierr = PetscSectionSetFieldDof(coordSectionB, v, 0, cDimB);CHKERRQ(ierr);
1304 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);
1305 }
1306 ierr = PetscSectionSetUp(coordSectionB);CHKERRQ(ierr);
1307 ierr = PetscSectionGetStorageSize(coordSectionB, &coordSize);CHKERRQ(ierr);
1308 ierr = VecCreate(PETSC_COMM_SELF, &coordinatesB);CHKERRQ(ierr);
1309 ierr = PetscObjectSetName((PetscObject) coordinatesB, "coordinates");CHKERRQ(ierr);
1310 ierr = VecSetSizes(coordinatesB, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
1311 ierr = VecSetBlockSize(coordinatesB, cDimB);CHKERRQ(ierr);
1312 ierr = VecSetType(coordinatesB,VECSTANDARD);CHKERRQ(ierr);
1313
1314 ierr = DMGetCoordinateSection(idm, &coordSectionA);CHKERRQ(ierr);
1315 ierr = DMGetCoordinatesLocal(idm, &coordinatesA);CHKERRQ(ierr);
1316 ierr = VecGetArray(coordinatesB, &coordsB);CHKERRQ(ierr);
1317 ierr = VecGetArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr);
1318 for (v = vStart; v < vEnd; ++v) {
1319 const PetscScalar *cptr;
1320 PetscReal ones2[2] = { 0., 1.}, ones3[3] = { 0., 0., 1.};
1321 PetscReal normal[3];
1322 PetscReal norm;
1323 PetscInt offA, d, cDimA = cDim;
1324
1325 if (normals) {for (d = 0; d < cDimB; ++d) normal[d] = normals[cDimB*(v - vStart)+d];}
1326 else if (haveCLNormal) {for (d = 0; d < cDimB; ++d) normal[d] = clNormal[d];}
1327 else if (extNormal) {for (d = 0; d < cDimB; ++d) normal[d] = extNormal[d];}
1328 else if (cDimB == 2) {for (d = 0; d < cDimB; ++d) normal[d] = ones2[d];}
1329 else if (cDimB == 3) {for (d = 0; d < cDimB; ++d) normal[d] = ones3[d];}
1330 else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unable to determine normal for extrusion");
1331 for (d = 0, norm = 0.0; d < cDimB; ++d) norm += normal[d]*normal[d];
1332 for (d = 0; d < cDimB; ++d) normal[d] *= 1./PetscSqrtReal(norm);
1333
1334 ierr = PetscSectionGetOffset(coordSectionA, v, &offA);CHKERRQ(ierr);
1335 cptr = coordsA + offA;
1336 for (l = 0; l <= layers; ++l) {
1337 PetscInt offB, d, newV;
1338
1339 newV = orderHeight ? (layers+1)*(v -vStart) + l + numCells : (vEnd -vStart)*l + (v -vStart) + numCells;
1340 ierr = PetscSectionGetOffset(coordSectionB, newV, &offB);CHKERRQ(ierr);
1341 for (d = 0; d < cDimA; ++d) { coordsB[offB+d] = cptr[d]; }
1342 for (d = 0; d < cDimB; ++d) { coordsB[offB+d] += l ? normal[d]*heights[l-1] : 0.0; }
1343 cptr = coordsB + offB;
1344 cDimA = cDimB;
1345 }
1346 }
1347 ierr = VecRestoreArrayRead(coordinatesA, &coordsA);CHKERRQ(ierr);
1348 ierr = VecRestoreArray(coordinatesB, &coordsB);CHKERRQ(ierr);
1349 ierr = DMSetCoordinatesLocal(*dm, coordinatesB);CHKERRQ(ierr);
1350 ierr = VecDestroy(&coordinatesB);CHKERRQ(ierr);
1351 ierr = PetscFree(normals);CHKERRQ(ierr);
1352 ierr = PetscFree(heights);CHKERRQ(ierr);
1353 if (interpolate) {
1354 DM idm;
1355
1356 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
1357 ierr = DMPlexCopyCoordinates(*dm, idm);CHKERRQ(ierr);
1358 ierr = DMDestroy(dm);CHKERRQ(ierr);
1359 *dm = idm;
1360 }
1361 PetscFunctionReturn(0);
1362 }
1363
1364 /*@C
1365 DMPlexSetOptionsPrefix - Sets the prefix used for searching for all DM options in the database.
1366
1367 Logically Collective on dm
1368
1369 Input Parameters:
1370 + dm - the DM context
1371 - prefix - the prefix to prepend to all option names
1372
1373 Notes:
1374 A hyphen (-) must NOT be given at the beginning of the prefix name.
1375 The first character of all runtime options is AUTOMATICALLY the hyphen.
1376
1377 Level: advanced
1378
1379 .seealso: SNESSetFromOptions()
1380 @*/
DMPlexSetOptionsPrefix(DM dm,const char prefix[])1381 PetscErrorCode DMPlexSetOptionsPrefix(DM dm, const char prefix[])
1382 {
1383 DM_Plex *mesh = (DM_Plex *) dm->data;
1384 PetscErrorCode ierr;
1385
1386 PetscFunctionBegin;
1387 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1388 ierr = PetscObjectSetOptionsPrefix((PetscObject) dm, prefix);CHKERRQ(ierr);
1389 ierr = PetscObjectSetOptionsPrefix((PetscObject) mesh->partitioner, prefix);CHKERRQ(ierr);
1390 PetscFunctionReturn(0);
1391 }
1392
1393 /*@
1394 DMPlexCreateHexCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using hexahedra.
1395
1396 Collective
1397
1398 Input Parameters:
1399 + comm - The communicator for the DM object
1400 . numRefine - The number of regular refinements to the basic 5 cell structure
1401 - periodicZ - The boundary type for the Z direction
1402
1403 Output Parameter:
1404 . dm - The DM object
1405
1406 Options Database Keys:
1407 These options override the hard-wired input values.
1408 + -dm_plex_hex_cyl_refine <r> - Refine the sylinder r times
1409 - -dm_plex_hex_cyl_bd <bz> - Specify the DMBoundaryType in the z-direction
1410
1411 Note:
1412 Here is the output numbering looking from the bottom of the cylinder:
1413 $ 17-----14
1414 $ | |
1415 $ | 2 |
1416 $ | |
1417 $ 17-----8-----7-----14
1418 $ | | | |
1419 $ | 3 | 0 | 1 |
1420 $ | | | |
1421 $ 19-----5-----6-----13
1422 $ | |
1423 $ | 4 |
1424 $ | |
1425 $ 19-----13
1426 $
1427 $ and up through the top
1428 $
1429 $ 18-----16
1430 $ | |
1431 $ | 2 |
1432 $ | |
1433 $ 18----10----11-----16
1434 $ | | | |
1435 $ | 3 | 0 | 1 |
1436 $ | | | |
1437 $ 20-----9----12-----15
1438 $ | |
1439 $ | 4 |
1440 $ | |
1441 $ 20-----15
1442
1443 Level: beginner
1444
1445 .seealso: DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1446 @*/
DMPlexCreateHexCylinderMesh(MPI_Comm comm,PetscInt numRefine,DMBoundaryType periodicZ,DM * dm)1447 PetscErrorCode DMPlexCreateHexCylinderMesh(MPI_Comm comm, PetscInt numRefine, DMBoundaryType periodicZ, DM *dm)
1448 {
1449 const PetscInt dim = 3;
1450 PetscInt numCells, numVertices, r;
1451 PetscMPIInt rank;
1452 PetscErrorCode ierr;
1453
1454 PetscFunctionBegin;
1455 PetscValidPointer(dm, 4);
1456 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1457 ierr = PetscOptionsGetInt(NULL, NULL, "-dm_plex_hex_cyl_refine", &numRefine, NULL);CHKERRQ(ierr);
1458 if (numRefine < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of refinements %D cannot be negative", numRefine);
1459 ierr = PetscOptionsGetEnum(NULL, NULL, "-dm_plex_box_bd", DMBoundaryTypes, (PetscEnum *) &periodicZ, NULL);CHKERRQ(ierr);
1460 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
1461 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
1462 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
1463 /* Create topology */
1464 {
1465 PetscInt cone[8], c;
1466
1467 numCells = !rank ? 5 : 0;
1468 numVertices = !rank ? 16 : 0;
1469 if (periodicZ == DM_BOUNDARY_PERIODIC) {
1470 numCells *= 3;
1471 numVertices = !rank ? 24 : 0;
1472 }
1473 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr);
1474 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 8);CHKERRQ(ierr);}
1475 ierr = DMSetUp(*dm);CHKERRQ(ierr);
1476 if (!rank) {
1477 if (periodicZ == DM_BOUNDARY_PERIODIC) {
1478 cone[0] = 15; cone[1] = 18; cone[2] = 17; cone[3] = 16;
1479 cone[4] = 31; cone[5] = 32; cone[6] = 33; cone[7] = 34;
1480 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr);
1481 cone[0] = 16; cone[1] = 17; cone[2] = 24; cone[3] = 23;
1482 cone[4] = 32; cone[5] = 36; cone[6] = 37; cone[7] = 33; /* 22 25 26 21 */
1483 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr);
1484 cone[0] = 18; cone[1] = 27; cone[2] = 24; cone[3] = 17;
1485 cone[4] = 34; cone[5] = 33; cone[6] = 37; cone[7] = 38;
1486 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr);
1487 cone[0] = 29; cone[1] = 27; cone[2] = 18; cone[3] = 15;
1488 cone[4] = 35; cone[5] = 31; cone[6] = 34; cone[7] = 38;
1489 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr);
1490 cone[0] = 29; cone[1] = 15; cone[2] = 16; cone[3] = 23;
1491 cone[4] = 35; cone[5] = 36; cone[6] = 32; cone[7] = 31;
1492 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr);
1493
1494 cone[0] = 31; cone[1] = 34; cone[2] = 33; cone[3] = 32;
1495 cone[4] = 19; cone[5] = 22; cone[6] = 21; cone[7] = 20;
1496 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr);
1497 cone[0] = 32; cone[1] = 33; cone[2] = 37; cone[3] = 36;
1498 cone[4] = 22; cone[5] = 25; cone[6] = 26; cone[7] = 21;
1499 ierr = DMPlexSetCone(*dm, 6, cone);CHKERRQ(ierr);
1500 cone[0] = 34; cone[1] = 38; cone[2] = 37; cone[3] = 33;
1501 cone[4] = 20; cone[5] = 21; cone[6] = 26; cone[7] = 28;
1502 ierr = DMPlexSetCone(*dm, 7, cone);CHKERRQ(ierr);
1503 cone[0] = 35; cone[1] = 38; cone[2] = 34; cone[3] = 31;
1504 cone[4] = 30; cone[5] = 19; cone[6] = 20; cone[7] = 28;
1505 ierr = DMPlexSetCone(*dm, 8, cone);CHKERRQ(ierr);
1506 cone[0] = 35; cone[1] = 31; cone[2] = 32; cone[3] = 36;
1507 cone[4] = 30; cone[5] = 25; cone[6] = 22; cone[7] = 19;
1508 ierr = DMPlexSetCone(*dm, 9, cone);CHKERRQ(ierr);
1509
1510 cone[0] = 19; cone[1] = 20; cone[2] = 21; cone[3] = 22;
1511 cone[4] = 15; cone[5] = 16; cone[6] = 17; cone[7] = 18;
1512 ierr = DMPlexSetCone(*dm, 10, cone);CHKERRQ(ierr);
1513 cone[0] = 22; cone[1] = 21; cone[2] = 26; cone[3] = 25;
1514 cone[4] = 16; cone[5] = 23; cone[6] = 24; cone[7] = 17;
1515 ierr = DMPlexSetCone(*dm, 11, cone);CHKERRQ(ierr);
1516 cone[0] = 20; cone[1] = 28; cone[2] = 26; cone[3] = 21;
1517 cone[4] = 18; cone[5] = 17; cone[6] = 24; cone[7] = 27;
1518 ierr = DMPlexSetCone(*dm, 12, cone);CHKERRQ(ierr);
1519 cone[0] = 30; cone[1] = 28; cone[2] = 20; cone[3] = 19;
1520 cone[4] = 29; cone[5] = 15; cone[6] = 18; cone[7] = 27;
1521 ierr = DMPlexSetCone(*dm, 13, cone);CHKERRQ(ierr);
1522 cone[0] = 30; cone[1] = 19; cone[2] = 22; cone[3] = 25;
1523 cone[4] = 29; cone[5] = 23; cone[6] = 16; cone[7] = 15;
1524 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr);
1525 } else {
1526 cone[0] = 5; cone[1] = 8; cone[2] = 7; cone[3] = 6;
1527 cone[4] = 9; cone[5] = 12; cone[6] = 11; cone[7] = 10;
1528 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr);
1529 cone[0] = 6; cone[1] = 7; cone[2] = 14; cone[3] = 13;
1530 cone[4] = 12; cone[5] = 15; cone[6] = 16; cone[7] = 11;
1531 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr);
1532 cone[0] = 8; cone[1] = 17; cone[2] = 14; cone[3] = 7;
1533 cone[4] = 10; cone[5] = 11; cone[6] = 16; cone[7] = 18;
1534 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr);
1535 cone[0] = 19; cone[1] = 17; cone[2] = 8; cone[3] = 5;
1536 cone[4] = 20; cone[5] = 9; cone[6] = 10; cone[7] = 18;
1537 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr);
1538 cone[0] = 19; cone[1] = 5; cone[2] = 6; cone[3] = 13;
1539 cone[4] = 20; cone[5] = 15; cone[6] = 12; cone[7] = 9;
1540 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr);
1541 }
1542 }
1543 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
1544 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
1545 }
1546 /* Interpolate */
1547 {
1548 DM idm;
1549
1550 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
1551 ierr = DMDestroy(dm);CHKERRQ(ierr);
1552 *dm = idm;
1553 }
1554 /* Create cube geometry */
1555 {
1556 Vec coordinates;
1557 PetscSection coordSection;
1558 PetscScalar *coords;
1559 PetscInt coordSize, v;
1560 const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1561 const PetscReal ds2 = dis/2.0;
1562
1563 /* Build coordinates */
1564 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
1565 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
1566 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
1567 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr);
1568 for (v = numCells; v < numCells+numVertices; ++v) {
1569 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
1570 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
1571 }
1572 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
1573 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
1574 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
1575 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
1576 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
1577 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
1578 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
1579 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
1580 if (!rank) {
1581 coords[0*dim+0] = -ds2; coords[0*dim+1] = -ds2; coords[0*dim+2] = 0.0;
1582 coords[1*dim+0] = ds2; coords[1*dim+1] = -ds2; coords[1*dim+2] = 0.0;
1583 coords[2*dim+0] = ds2; coords[2*dim+1] = ds2; coords[2*dim+2] = 0.0;
1584 coords[3*dim+0] = -ds2; coords[3*dim+1] = ds2; coords[3*dim+2] = 0.0;
1585 coords[4*dim+0] = -ds2; coords[4*dim+1] = -ds2; coords[4*dim+2] = 1.0;
1586 coords[5*dim+0] = -ds2; coords[5*dim+1] = ds2; coords[5*dim+2] = 1.0;
1587 coords[6*dim+0] = ds2; coords[6*dim+1] = ds2; coords[6*dim+2] = 1.0;
1588 coords[7*dim+0] = ds2; coords[7*dim+1] = -ds2; coords[7*dim+2] = 1.0;
1589 coords[ 8*dim+0] = dis; coords[ 8*dim+1] = -dis; coords[ 8*dim+2] = 0.0;
1590 coords[ 9*dim+0] = dis; coords[ 9*dim+1] = dis; coords[ 9*dim+2] = 0.0;
1591 coords[10*dim+0] = dis; coords[10*dim+1] = -dis; coords[10*dim+2] = 1.0;
1592 coords[11*dim+0] = dis; coords[11*dim+1] = dis; coords[11*dim+2] = 1.0;
1593 coords[12*dim+0] = -dis; coords[12*dim+1] = dis; coords[12*dim+2] = 0.0;
1594 coords[13*dim+0] = -dis; coords[13*dim+1] = dis; coords[13*dim+2] = 1.0;
1595 coords[14*dim+0] = -dis; coords[14*dim+1] = -dis; coords[14*dim+2] = 0.0;
1596 coords[15*dim+0] = -dis; coords[15*dim+1] = -dis; coords[15*dim+2] = 1.0;
1597 if (periodicZ == DM_BOUNDARY_PERIODIC) {
1598 /* 15 31 19 */ coords[16*dim+0] = -ds2; coords[16*dim+1] = -ds2; coords[16*dim+2] = 0.5;
1599 /* 16 32 22 */ coords[17*dim+0] = ds2; coords[17*dim+1] = -ds2; coords[17*dim+2] = 0.5;
1600 /* 17 33 21 */ coords[18*dim+0] = ds2; coords[18*dim+1] = ds2; coords[18*dim+2] = 0.5;
1601 /* 18 34 20 */ coords[19*dim+0] = -ds2; coords[19*dim+1] = ds2; coords[19*dim+2] = 0.5;
1602 /* 29 35 30 */ coords[20*dim+0] = -dis; coords[20*dim+1] = -dis; coords[20*dim+2] = 0.5;
1603 /* 23 36 25 */ coords[21*dim+0] = dis; coords[21*dim+1] = -dis; coords[21*dim+2] = 0.5;
1604 /* 24 37 26 */ coords[22*dim+0] = dis; coords[22*dim+1] = dis; coords[22*dim+2] = 0.5;
1605 /* 27 38 28 */ coords[23*dim+0] = -dis; coords[23*dim+1] = dis; coords[23*dim+2] = 0.5;
1606 }
1607 }
1608 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
1609 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
1610 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
1611 }
1612 /* Create periodicity */
1613 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1614 PetscReal L[3];
1615 PetscReal maxCell[3];
1616 DMBoundaryType bdType[3];
1617 PetscReal lower[3] = {0.0, 0.0, 0.0};
1618 PetscReal upper[3] = {1.0, 1.0, 1.5};
1619 PetscInt i, numZCells = 3;
1620
1621 bdType[0] = DM_BOUNDARY_NONE;
1622 bdType[1] = DM_BOUNDARY_NONE;
1623 bdType[2] = periodicZ;
1624 for (i = 0; i < dim; i++) {
1625 L[i] = upper[i] - lower[i];
1626 maxCell[i] = 1.1 * (L[i] / numZCells);
1627 }
1628 ierr = DMSetPeriodicity(*dm, PETSC_TRUE, maxCell, L, bdType);CHKERRQ(ierr);
1629 }
1630 /* Refine topology */
1631 for (r = 0; r < numRefine; ++r) {
1632 DM rdm = NULL;
1633
1634 ierr = DMRefine(*dm, comm, &rdm);CHKERRQ(ierr);
1635 ierr = DMDestroy(dm);CHKERRQ(ierr);
1636 *dm = rdm;
1637 }
1638 /* Remap geometry to cylinder
1639 Interior square: Linear interpolation is correct
1640 The other cells all have vertices on rays from the origin. We want to uniformly expand the spacing
1641 such that the last vertex is on the unit circle. So the closest and farthest vertices are at distance
1642
1643 phi = arctan(y/x)
1644 d_close = sqrt(1/8 + 1/4 sin^2(phi))
1645 d_far = sqrt(1/2 + sin^2(phi))
1646
1647 so we remap them using
1648
1649 x_new = x_close + (x - x_close) (1 - d_close) / (d_far - d_close)
1650 y_new = y_close + (y - y_close) (1 - d_close) / (d_far - d_close)
1651
1652 If pi/4 < phi < 3pi/4 or -3pi/4 < phi < -pi/4, then we switch x and y.
1653 */
1654 {
1655 Vec coordinates;
1656 PetscSection coordSection;
1657 PetscScalar *coords;
1658 PetscInt vStart, vEnd, v;
1659 const PetscReal dis = 1.0/PetscSqrtReal(2.0);
1660 const PetscReal ds2 = 0.5*dis;
1661
1662 ierr = DMPlexGetDepthStratum(*dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
1663 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
1664 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr);
1665 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
1666 for (v = vStart; v < vEnd; ++v) {
1667 PetscReal phi, sinp, cosp, dc, df, x, y, xc, yc;
1668 PetscInt off;
1669
1670 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr);
1671 if ((PetscAbsScalar(coords[off+0]) <= ds2) && (PetscAbsScalar(coords[off+1]) <= ds2)) continue;
1672 x = PetscRealPart(coords[off]);
1673 y = PetscRealPart(coords[off+1]);
1674 phi = PetscAtan2Real(y, x);
1675 sinp = PetscSinReal(phi);
1676 cosp = PetscCosReal(phi);
1677 if ((PetscAbsReal(phi) > PETSC_PI/4.0) && (PetscAbsReal(phi) < 3.0*PETSC_PI/4.0)) {
1678 dc = PetscAbsReal(ds2/sinp);
1679 df = PetscAbsReal(dis/sinp);
1680 xc = ds2*x/PetscAbsReal(y);
1681 yc = ds2*PetscSignReal(y);
1682 } else {
1683 dc = PetscAbsReal(ds2/cosp);
1684 df = PetscAbsReal(dis/cosp);
1685 xc = ds2*PetscSignReal(x);
1686 yc = ds2*y/PetscAbsReal(x);
1687 }
1688 coords[off+0] = xc + (coords[off+0] - xc)*(1.0 - dc)/(df - dc);
1689 coords[off+1] = yc + (coords[off+1] - yc)*(1.0 - dc)/(df - dc);
1690 }
1691 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
1692 if (periodicZ == DM_BOUNDARY_PERIODIC || periodicZ == DM_BOUNDARY_TWIST) {
1693 ierr = DMLocalizeCoordinates(*dm);CHKERRQ(ierr);
1694 }
1695 }
1696 PetscFunctionReturn(0);
1697 }
1698
1699 /*@
1700 DMPlexCreateWedgeCylinderMesh - Creates a mesh on the tensor product of the unit interval with the circle (cylinder) using wedges.
1701
1702 Collective
1703
1704 Input Parameters:
1705 + comm - The communicator for the DM object
1706 . n - The number of wedges around the origin
1707 - interpolate - Create edges and faces
1708
1709 Output Parameter:
1710 . dm - The DM object
1711
1712 Level: beginner
1713
1714 .seealso: DMPlexCreateHexCylinderMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1715 @*/
DMPlexCreateWedgeCylinderMesh(MPI_Comm comm,PetscInt n,PetscBool interpolate,DM * dm)1716 PetscErrorCode DMPlexCreateWedgeCylinderMesh(MPI_Comm comm, PetscInt n, PetscBool interpolate, DM *dm)
1717 {
1718 const PetscInt dim = 3;
1719 PetscInt numCells, numVertices, v;
1720 PetscMPIInt rank;
1721 PetscErrorCode ierr;
1722
1723 PetscFunctionBegin;
1724 PetscValidPointer(dm, 3);
1725 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
1726 if (n < 0) SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "Number of wedges %D cannot be negative", n);
1727 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
1728 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
1729 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
1730 /* Must create the celltype label here so that we do not automatically try to compute the types */
1731 ierr = DMCreateLabel(*dm, "celltype");CHKERRQ(ierr);
1732 /* Create topology */
1733 {
1734 PetscInt cone[6], c;
1735
1736 numCells = !rank ? n : 0;
1737 numVertices = !rank ? 2*(n+1) : 0;
1738 ierr = DMPlexSetChart(*dm, 0, numCells+numVertices);CHKERRQ(ierr);
1739 for (c = 0; c < numCells; c++) {ierr = DMPlexSetConeSize(*dm, c, 6);CHKERRQ(ierr);}
1740 ierr = DMSetUp(*dm);CHKERRQ(ierr);
1741 for (c = 0; c < numCells; c++) {
1742 cone[0] = c+n*1; cone[1] = (c+1)%n+n*1; cone[2] = 0+3*n;
1743 cone[3] = c+n*2; cone[4] = (c+1)%n+n*2; cone[5] = 1+3*n;
1744 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr);
1745 ierr = DMPlexSetCellType(*dm, c, DM_POLYTOPE_TRI_PRISM_TENSOR);CHKERRQ(ierr);
1746 }
1747 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
1748 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
1749 }
1750 for (v = numCells; v < numCells+numVertices; ++v) {
1751 ierr = DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);
1752 }
1753 /* Interpolate */
1754 if (interpolate) {
1755 DM idm;
1756
1757 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
1758 ierr = DMDestroy(dm);CHKERRQ(ierr);
1759 *dm = idm;
1760 }
1761 /* Create cylinder geometry */
1762 {
1763 Vec coordinates;
1764 PetscSection coordSection;
1765 PetscScalar *coords;
1766 PetscInt coordSize, c;
1767
1768 /* Build coordinates */
1769 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
1770 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
1771 ierr = PetscSectionSetFieldComponents(coordSection, 0, dim);CHKERRQ(ierr);
1772 ierr = PetscSectionSetChart(coordSection, numCells, numCells+numVertices);CHKERRQ(ierr);
1773 for (v = numCells; v < numCells+numVertices; ++v) {
1774 ierr = PetscSectionSetDof(coordSection, v, dim);CHKERRQ(ierr);
1775 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dim);CHKERRQ(ierr);
1776 }
1777 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
1778 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
1779 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
1780 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
1781 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
1782 ierr = VecSetBlockSize(coordinates, dim);CHKERRQ(ierr);
1783 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
1784 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
1785 for (c = 0; c < numCells; c++) {
1786 coords[(c+0*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+0*n)*dim+2] = 1.0;
1787 coords[(c+1*n)*dim+0] = PetscCosReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+1] = PetscSinReal(2.0*c*PETSC_PI/n); coords[(c+1*n)*dim+2] = 0.0;
1788 }
1789 if (!rank) {
1790 coords[(2*n+0)*dim+0] = 0.0; coords[(2*n+0)*dim+1] = 0.0; coords[(2*n+0)*dim+2] = 1.0;
1791 coords[(2*n+1)*dim+0] = 0.0; coords[(2*n+1)*dim+1] = 0.0; coords[(2*n+1)*dim+2] = 0.0;
1792 }
1793 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
1794 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
1795 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
1796 }
1797 PetscFunctionReturn(0);
1798 }
1799
DiffNormReal(PetscInt dim,const PetscReal x[],const PetscReal y[])1800 PETSC_STATIC_INLINE PetscReal DiffNormReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1801 {
1802 PetscReal prod = 0.0;
1803 PetscInt i;
1804 for (i = 0; i < dim; ++i) prod += PetscSqr(x[i] - y[i]);
1805 return PetscSqrtReal(prod);
1806 }
DotReal(PetscInt dim,const PetscReal x[],const PetscReal y[])1807 PETSC_STATIC_INLINE PetscReal DotReal(PetscInt dim, const PetscReal x[], const PetscReal y[])
1808 {
1809 PetscReal prod = 0.0;
1810 PetscInt i;
1811 for (i = 0; i < dim; ++i) prod += x[i]*y[i];
1812 return prod;
1813 }
1814
1815 /* The first constant is the sphere radius */
snapToSphere(PetscInt dim,PetscInt Nf,PetscInt NfAux,const PetscInt uOff[],const PetscInt uOff_x[],const PetscScalar u[],const PetscScalar u_t[],const PetscScalar u_x[],const PetscInt aOff[],const PetscInt aOff_x[],const PetscScalar a[],const PetscScalar a_t[],const PetscScalar a_x[],PetscReal t,const PetscReal x[],PetscInt numConstants,const PetscScalar constants[],PetscScalar f0[])1816 static void snapToSphere(PetscInt dim, PetscInt Nf, PetscInt NfAux,
1817 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
1818 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
1819 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
1820 {
1821 PetscReal r = PetscRealPart(constants[0]);
1822 PetscReal norm2 = 0.0, fac;
1823 PetscInt n = uOff[1] - uOff[0], d;
1824
1825 for (d = 0; d < n; ++d) norm2 += PetscSqr(PetscRealPart(u[d]));
1826 fac = r/PetscSqrtReal(norm2);
1827 for (d = 0; d < n; ++d) f0[d] = u[d]*fac;
1828 }
1829
1830 /*@
1831 DMPlexCreateSphereMesh - Creates a mesh on the d-dimensional sphere, S^d.
1832
1833 Collective
1834
1835 Input Parameters:
1836 + comm - The communicator for the DM object
1837 . dim - The dimension
1838 . simplex - Use simplices, or tensor product cells
1839 - R - The radius
1840
1841 Output Parameter:
1842 . dm - The DM object
1843
1844 Level: beginner
1845
1846 .seealso: DMPlexCreateBallMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
1847 @*/
DMPlexCreateSphereMesh(MPI_Comm comm,PetscInt dim,PetscBool simplex,PetscReal R,DM * dm)1848 PetscErrorCode DMPlexCreateSphereMesh(MPI_Comm comm, PetscInt dim, PetscBool simplex, PetscReal R, DM *dm)
1849 {
1850 const PetscInt embedDim = dim+1;
1851 PetscSection coordSection;
1852 Vec coordinates;
1853 PetscScalar *coords;
1854 PetscReal *coordsIn;
1855 PetscInt numCells, numEdges, numVerts, firstVertex, v, firstEdge, coordSize, d, c, e;
1856 PetscMPIInt rank;
1857 PetscErrorCode ierr;
1858
1859 PetscFunctionBegin;
1860 PetscValidPointer(dm, 4);
1861 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
1862 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
1863 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
1864 ierr = DMSetCoordinateDim(*dm, dim+1);CHKERRQ(ierr);
1865 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) *dm), &rank);CHKERRQ(ierr);
1866 switch (dim) {
1867 case 2:
1868 if (simplex) {
1869 DM idm;
1870 const PetscReal radius = PetscSqrtReal(1 + PETSC_PHI*PETSC_PHI)/(1.0 + PETSC_PHI);
1871 const PetscReal edgeLen = 2.0/(1.0 + PETSC_PHI) * (R/radius);
1872 const PetscInt degree = 5;
1873 PetscReal vertex[3] = {0.0, 1.0/(1.0 + PETSC_PHI), PETSC_PHI/(1.0 + PETSC_PHI)};
1874 PetscInt s[3] = {1, 1, 1};
1875 PetscInt cone[3];
1876 PetscInt *graph, p, i, j, k;
1877
1878 vertex[0] *= R/radius; vertex[1] *= R/radius; vertex[2] *= R/radius;
1879 numCells = !rank ? 20 : 0;
1880 numVerts = !rank ? 12 : 0;
1881 firstVertex = numCells;
1882 /* Use icosahedron, which for a R-sphere has coordinates which are all cyclic permutations of
1883
1884 (0, \pm 1/\phi+1, \pm \phi/\phi+1)
1885
1886 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
1887 length is then given by 2/(1+\phi) = 2 * 0.38197 = 0.76393.
1888 */
1889 /* Construct vertices */
1890 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr);
1891 if (!rank) {
1892 for (p = 0, i = 0; p < embedDim; ++p) {
1893 for (s[1] = -1; s[1] < 2; s[1] += 2) {
1894 for (s[2] = -1; s[2] < 2; s[2] += 2) {
1895 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertex[(d+p)%embedDim];
1896 ++i;
1897 }
1898 }
1899 }
1900 }
1901 /* Construct graph */
1902 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr);
1903 for (i = 0; i < numVerts; ++i) {
1904 for (j = 0, k = 0; j < numVerts; ++j) {
1905 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
1906 }
1907 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid icosahedron, vertex %D degree %D != %D", i, k, degree);
1908 }
1909 /* Build Topology */
1910 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr);
1911 for (c = 0; c < numCells; c++) {
1912 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr);
1913 }
1914 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */
1915 /* Cells */
1916 for (i = 0, c = 0; i < numVerts; ++i) {
1917 for (j = 0; j < i; ++j) {
1918 for (k = 0; k < j; ++k) {
1919 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i]) {
1920 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k;
1921 /* Check orientation */
1922 {
1923 const PetscInt epsilon[3][3][3] = {{{0, 0, 0}, {0, 0, 1}, {0, -1, 0}}, {{0, 0, -1}, {0, 0, 0}, {1, 0, 0}}, {{0, 1, 0}, {-1, 0, 0}, {0, 0, 0}}};
1924 PetscReal normal[3];
1925 PetscInt e, f;
1926
1927 for (d = 0; d < embedDim; ++d) {
1928 normal[d] = 0.0;
1929 for (e = 0; e < embedDim; ++e) {
1930 for (f = 0; f < embedDim; ++f) {
1931 normal[d] += epsilon[d][e][f]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f]);
1932 }
1933 }
1934 }
1935 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
1936 }
1937 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr);
1938 }
1939 }
1940 }
1941 }
1942 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
1943 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
1944 ierr = PetscFree(graph);CHKERRQ(ierr);
1945 /* Interpolate mesh */
1946 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
1947 ierr = DMDestroy(dm);CHKERRQ(ierr);
1948 *dm = idm;
1949 } else {
1950 /*
1951 12-21--13
1952 | |
1953 25 4 24
1954 | |
1955 12-25--9-16--8-24--13
1956 | | | |
1957 23 5 17 0 15 3 22
1958 | | | |
1959 10-20--6-14--7-19--11
1960 | |
1961 20 1 19
1962 | |
1963 10-18--11
1964 | |
1965 23 2 22
1966 | |
1967 12-21--13
1968 */
1969 PetscInt cone[4], ornt[4];
1970
1971 numCells = !rank ? 6 : 0;
1972 numEdges = !rank ? 12 : 0;
1973 numVerts = !rank ? 8 : 0;
1974 firstVertex = numCells;
1975 firstEdge = numCells + numVerts;
1976 /* Build Topology */
1977 ierr = DMPlexSetChart(*dm, 0, numCells+numEdges+numVerts);CHKERRQ(ierr);
1978 for (c = 0; c < numCells; c++) {
1979 ierr = DMPlexSetConeSize(*dm, c, 4);CHKERRQ(ierr);
1980 }
1981 for (e = firstEdge; e < firstEdge+numEdges; ++e) {
1982 ierr = DMPlexSetConeSize(*dm, e, 2);CHKERRQ(ierr);
1983 }
1984 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */
1985 if (!rank) {
1986 /* Cell 0 */
1987 cone[0] = 14; cone[1] = 15; cone[2] = 16; cone[3] = 17;
1988 ierr = DMPlexSetCone(*dm, 0, cone);CHKERRQ(ierr);
1989 ornt[0] = 0; ornt[1] = 0; ornt[2] = 0; ornt[3] = 0;
1990 ierr = DMPlexSetConeOrientation(*dm, 0, ornt);CHKERRQ(ierr);
1991 /* Cell 1 */
1992 cone[0] = 18; cone[1] = 19; cone[2] = 14; cone[3] = 20;
1993 ierr = DMPlexSetCone(*dm, 1, cone);CHKERRQ(ierr);
1994 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0;
1995 ierr = DMPlexSetConeOrientation(*dm, 1, ornt);CHKERRQ(ierr);
1996 /* Cell 2 */
1997 cone[0] = 21; cone[1] = 22; cone[2] = 18; cone[3] = 23;
1998 ierr = DMPlexSetCone(*dm, 2, cone);CHKERRQ(ierr);
1999 ornt[0] = 0; ornt[1] = 0; ornt[2] = -2; ornt[3] = 0;
2000 ierr = DMPlexSetConeOrientation(*dm, 2, ornt);CHKERRQ(ierr);
2001 /* Cell 3 */
2002 cone[0] = 19; cone[1] = 22; cone[2] = 24; cone[3] = 15;
2003 ierr = DMPlexSetCone(*dm, 3, cone);CHKERRQ(ierr);
2004 ornt[0] = -2; ornt[1] = -2; ornt[2] = 0; ornt[3] = -2;
2005 ierr = DMPlexSetConeOrientation(*dm, 3, ornt);CHKERRQ(ierr);
2006 /* Cell 4 */
2007 cone[0] = 16; cone[1] = 24; cone[2] = 21; cone[3] = 25;
2008 ierr = DMPlexSetCone(*dm, 4, cone);CHKERRQ(ierr);
2009 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = 0;
2010 ierr = DMPlexSetConeOrientation(*dm, 4, ornt);CHKERRQ(ierr);
2011 /* Cell 5 */
2012 cone[0] = 20; cone[1] = 17; cone[2] = 25; cone[3] = 23;
2013 ierr = DMPlexSetCone(*dm, 5, cone);CHKERRQ(ierr);
2014 ornt[0] = -2; ornt[1] = -2; ornt[2] = -2; ornt[3] = -2;
2015 ierr = DMPlexSetConeOrientation(*dm, 5, ornt);CHKERRQ(ierr);
2016 /* Edges */
2017 cone[0] = 6; cone[1] = 7;
2018 ierr = DMPlexSetCone(*dm, 14, cone);CHKERRQ(ierr);
2019 cone[0] = 7; cone[1] = 8;
2020 ierr = DMPlexSetCone(*dm, 15, cone);CHKERRQ(ierr);
2021 cone[0] = 8; cone[1] = 9;
2022 ierr = DMPlexSetCone(*dm, 16, cone);CHKERRQ(ierr);
2023 cone[0] = 9; cone[1] = 6;
2024 ierr = DMPlexSetCone(*dm, 17, cone);CHKERRQ(ierr);
2025 cone[0] = 10; cone[1] = 11;
2026 ierr = DMPlexSetCone(*dm, 18, cone);CHKERRQ(ierr);
2027 cone[0] = 11; cone[1] = 7;
2028 ierr = DMPlexSetCone(*dm, 19, cone);CHKERRQ(ierr);
2029 cone[0] = 6; cone[1] = 10;
2030 ierr = DMPlexSetCone(*dm, 20, cone);CHKERRQ(ierr);
2031 cone[0] = 12; cone[1] = 13;
2032 ierr = DMPlexSetCone(*dm, 21, cone);CHKERRQ(ierr);
2033 cone[0] = 13; cone[1] = 11;
2034 ierr = DMPlexSetCone(*dm, 22, cone);CHKERRQ(ierr);
2035 cone[0] = 10; cone[1] = 12;
2036 ierr = DMPlexSetCone(*dm, 23, cone);CHKERRQ(ierr);
2037 cone[0] = 13; cone[1] = 8;
2038 ierr = DMPlexSetCone(*dm, 24, cone);CHKERRQ(ierr);
2039 cone[0] = 12; cone[1] = 9;
2040 ierr = DMPlexSetCone(*dm, 25, cone);CHKERRQ(ierr);
2041 }
2042 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
2043 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
2044 /* Build coordinates */
2045 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr);
2046 if (!rank) {
2047 coordsIn[0*embedDim+0] = -R; coordsIn[0*embedDim+1] = R; coordsIn[0*embedDim+2] = -R;
2048 coordsIn[1*embedDim+0] = R; coordsIn[1*embedDim+1] = R; coordsIn[1*embedDim+2] = -R;
2049 coordsIn[2*embedDim+0] = R; coordsIn[2*embedDim+1] = -R; coordsIn[2*embedDim+2] = -R;
2050 coordsIn[3*embedDim+0] = -R; coordsIn[3*embedDim+1] = -R; coordsIn[3*embedDim+2] = -R;
2051 coordsIn[4*embedDim+0] = -R; coordsIn[4*embedDim+1] = R; coordsIn[4*embedDim+2] = R;
2052 coordsIn[5*embedDim+0] = R; coordsIn[5*embedDim+1] = R; coordsIn[5*embedDim+2] = R;
2053 coordsIn[6*embedDim+0] = -R; coordsIn[6*embedDim+1] = -R; coordsIn[6*embedDim+2] = R;
2054 coordsIn[7*embedDim+0] = R; coordsIn[7*embedDim+1] = -R; coordsIn[7*embedDim+2] = R;
2055 }
2056 }
2057 break;
2058 case 3:
2059 if (simplex) {
2060 DM idm;
2061 const PetscReal edgeLen = 1.0/PETSC_PHI;
2062 PetscReal vertexA[4] = {0.5, 0.5, 0.5, 0.5};
2063 PetscReal vertexB[4] = {1.0, 0.0, 0.0, 0.0};
2064 PetscReal vertexC[4] = {0.5, 0.5*PETSC_PHI, 0.5/PETSC_PHI, 0.0};
2065 const PetscInt degree = 12;
2066 PetscInt s[4] = {1, 1, 1};
2067 PetscInt evenPerm[12][4] = {{0, 1, 2, 3}, {0, 2, 3, 1}, {0, 3, 1, 2}, {1, 0, 3, 2}, {1, 2, 0, 3}, {1, 3, 2, 0},
2068 {2, 0, 1, 3}, {2, 1, 3, 0}, {2, 3, 0, 1}, {3, 0, 2, 1}, {3, 1, 0, 2}, {3, 2, 1, 0}};
2069 PetscInt cone[4];
2070 PetscInt *graph, p, i, j, k, l;
2071
2072 vertexA[0] *= R; vertexA[1] *= R; vertexA[2] *= R; vertexA[3] *= R;
2073 vertexB[0] *= R; vertexB[1] *= R; vertexB[2] *= R; vertexB[3] *= R;
2074 vertexC[0] *= R; vertexC[1] *= R; vertexC[2] *= R; vertexC[3] *= R;
2075 numCells = !rank ? 600 : 0;
2076 numVerts = !rank ? 120 : 0;
2077 firstVertex = numCells;
2078 /* Use the 600-cell, which for a unit sphere has coordinates which are
2079
2080 1/2 (\pm 1, \pm 1, \pm 1, \pm 1) 16
2081 (\pm 1, 0, 0, 0) all cyclic permutations 8
2082 1/2 (\pm 1, \pm phi, \pm 1/phi, 0) all even permutations 96
2083
2084 where \phi^2 - \phi - 1 = 0, meaning \phi is the golden ratio \frac{1 + \sqrt{5}}{2}. The edge
2085 length is then given by 1/\phi = 0.61803.
2086
2087 http://buzzard.pugetsound.edu/sage-practice/ch03s03.html
2088 http://mathworld.wolfram.com/600-Cell.html
2089 */
2090 /* Construct vertices */
2091 ierr = PetscCalloc1(numVerts * embedDim, &coordsIn);CHKERRQ(ierr);
2092 i = 0;
2093 if (!rank) {
2094 for (s[0] = -1; s[0] < 2; s[0] += 2) {
2095 for (s[1] = -1; s[1] < 2; s[1] += 2) {
2096 for (s[2] = -1; s[2] < 2; s[2] += 2) {
2097 for (s[3] = -1; s[3] < 2; s[3] += 2) {
2098 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[d]*vertexA[d];
2099 ++i;
2100 }
2101 }
2102 }
2103 }
2104 for (p = 0; p < embedDim; ++p) {
2105 s[1] = s[2] = s[3] = 1;
2106 for (s[0] = -1; s[0] < 2; s[0] += 2) {
2107 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[(d+p)%embedDim]*vertexB[(d+p)%embedDim];
2108 ++i;
2109 }
2110 }
2111 for (p = 0; p < 12; ++p) {
2112 s[3] = 1;
2113 for (s[0] = -1; s[0] < 2; s[0] += 2) {
2114 for (s[1] = -1; s[1] < 2; s[1] += 2) {
2115 for (s[2] = -1; s[2] < 2; s[2] += 2) {
2116 for (d = 0; d < embedDim; ++d) coordsIn[i*embedDim+d] = s[evenPerm[p][d]]*vertexC[evenPerm[p][d]];
2117 ++i;
2118 }
2119 }
2120 }
2121 }
2122 }
2123 if (i != numVerts) SETERRQ2(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertices %D != %D", i, numVerts);
2124 /* Construct graph */
2125 ierr = PetscCalloc1(numVerts * numVerts, &graph);CHKERRQ(ierr);
2126 for (i = 0; i < numVerts; ++i) {
2127 for (j = 0, k = 0; j < numVerts; ++j) {
2128 if (PetscAbsReal(DiffNormReal(embedDim, &coordsIn[i*embedDim], &coordsIn[j*embedDim]) - edgeLen) < PETSC_SMALL) {graph[i*numVerts+j] = 1; ++k;}
2129 }
2130 if (k != degree) SETERRQ3(comm, PETSC_ERR_PLIB, "Invalid 600-cell, vertex %D degree %D != %D", i, k, degree);
2131 }
2132 /* Build Topology */
2133 ierr = DMPlexSetChart(*dm, 0, numCells+numVerts);CHKERRQ(ierr);
2134 for (c = 0; c < numCells; c++) {
2135 ierr = DMPlexSetConeSize(*dm, c, embedDim);CHKERRQ(ierr);
2136 }
2137 ierr = DMSetUp(*dm);CHKERRQ(ierr); /* Allocate space for cones */
2138 /* Cells */
2139 if (!rank) {
2140 for (i = 0, c = 0; i < numVerts; ++i) {
2141 for (j = 0; j < i; ++j) {
2142 for (k = 0; k < j; ++k) {
2143 for (l = 0; l < k; ++l) {
2144 if (graph[i*numVerts+j] && graph[j*numVerts+k] && graph[k*numVerts+i] &&
2145 graph[l*numVerts+i] && graph[l*numVerts+j] && graph[l*numVerts+k]) {
2146 cone[0] = firstVertex+i; cone[1] = firstVertex+j; cone[2] = firstVertex+k; cone[3] = firstVertex+l;
2147 /* Check orientation: https://ef.gy/linear-algebra:normal-vectors-in-higher-dimensional-spaces */
2148 {
2149 const PetscInt epsilon[4][4][4][4] = {{{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2150 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, -1, 0}},
2151 {{0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 1, 0, 0}},
2152 {{0, 0, 0, 0}, { 0, 0, 1, 0}, { 0, -1, 0, 0}, { 0, 0, 0, 0}}},
2153
2154 {{{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, -1}, { 0, 0, 1, 0}},
2155 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2156 {{0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}},
2157 {{0, 0, -1, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}}},
2158
2159 {{{0, 0, 0, 0}, { 0, 0, 0, 1}, { 0, 0, 0, 0}, { 0, -1, 0, 0}},
2160 {{0, 0, 0, -1}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 1, 0, 0, 0}},
2161 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2162 {{0, 1, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}},
2163
2164 {{{0, 0, 0, 0}, { 0, 0, -1, 0}, { 0, 1, 0, 0}, { 0, 0, 0, 0}},
2165 {{0, 0, 1, 0}, { 0, 0, 0, 0}, {-1, 0, 0, 0}, { 0, 0, 0, 0}},
2166 {{0, -1, 0, 0}, { 1, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}},
2167 {{0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}, { 0, 0, 0, 0}}}};
2168 PetscReal normal[4];
2169 PetscInt e, f, g;
2170
2171 for (d = 0; d < embedDim; ++d) {
2172 normal[d] = 0.0;
2173 for (e = 0; e < embedDim; ++e) {
2174 for (f = 0; f < embedDim; ++f) {
2175 for (g = 0; g < embedDim; ++g) {
2176 normal[d] += epsilon[d][e][f][g]*(coordsIn[j*embedDim+e] - coordsIn[i*embedDim+e])*(coordsIn[k*embedDim+f] - coordsIn[i*embedDim+f])*(coordsIn[l*embedDim+f] - coordsIn[i*embedDim+f]);
2177 }
2178 }
2179 }
2180 }
2181 if (DotReal(embedDim, normal, &coordsIn[i*embedDim]) < 0) {PetscInt tmp = cone[1]; cone[1] = cone[2]; cone[2] = tmp;}
2182 }
2183 ierr = DMPlexSetCone(*dm, c++, cone);CHKERRQ(ierr);
2184 }
2185 }
2186 }
2187 }
2188 }
2189 }
2190 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
2191 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
2192 ierr = PetscFree(graph);CHKERRQ(ierr);
2193 /* Interpolate mesh */
2194 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
2195 ierr = DMDestroy(dm);CHKERRQ(ierr);
2196 *dm = idm;
2197 break;
2198 }
2199 default: SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported dimension for sphere: %D", dim);
2200 }
2201 /* Create coordinates */
2202 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
2203 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
2204 ierr = PetscSectionSetFieldComponents(coordSection, 0, embedDim);CHKERRQ(ierr);
2205 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numVerts);CHKERRQ(ierr);
2206 for (v = firstVertex; v < firstVertex+numVerts; ++v) {
2207 ierr = PetscSectionSetDof(coordSection, v, embedDim);CHKERRQ(ierr);
2208 ierr = PetscSectionSetFieldDof(coordSection, v, 0, embedDim);CHKERRQ(ierr);
2209 }
2210 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
2211 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
2212 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
2213 ierr = VecSetBlockSize(coordinates, embedDim);CHKERRQ(ierr);
2214 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
2215 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
2216 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
2217 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
2218 for (v = 0; v < numVerts; ++v) for (d = 0; d < embedDim; ++d) {coords[v*embedDim+d] = coordsIn[v*embedDim+d];}
2219 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
2220 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
2221 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
2222 ierr = PetscFree(coordsIn);CHKERRQ(ierr);
2223 /* Create coordinate function space */
2224 {
2225 DM cdm;
2226 PetscDS cds;
2227 PetscFE fe;
2228 PetscScalar radius = R;
2229 PetscInt dT, dE;
2230
2231 ierr = DMGetCoordinateDM(*dm, &cdm);CHKERRQ(ierr);
2232 ierr = DMGetDimension(*dm, &dT);CHKERRQ(ierr);
2233 ierr = DMGetCoordinateDim(*dm, &dE);CHKERRQ(ierr);
2234 ierr = PetscFECreateLagrange(PETSC_COMM_SELF, dT, dE, simplex, 1, -1, &fe);CHKERRQ(ierr);
2235 ierr = DMSetField(cdm, 0, NULL, (PetscObject) fe);CHKERRQ(ierr);
2236 ierr = PetscFEDestroy(&fe);CHKERRQ(ierr);
2237 ierr = DMCreateDS(cdm);CHKERRQ(ierr);
2238
2239 ierr = DMGetDS(cdm, &cds);CHKERRQ(ierr);
2240 ierr = PetscDSSetConstants(cds, 1, &radius);CHKERRQ(ierr);
2241 }
2242 ((DM_Plex *) (*dm)->data)->coordFunc = snapToSphere;
2243 PetscFunctionReturn(0);
2244 }
2245
2246 /*@
2247 DMPlexCreateBallMesh - Creates a simplex mesh on the d-dimensional ball, B^d.
2248
2249 Collective
2250
2251 Input Parameters:
2252 + comm - The communicator for the DM object
2253 . dim - The dimension
2254 - R - The radius
2255
2256 Output Parameter:
2257 . dm - The DM object
2258
2259 Options Database Keys:
2260 - bd_dm_refine - This will refine the surface mesh preserving the sphere geometry
2261
2262 Level: beginner
2263
2264 .seealso: DMPlexCreateSphereMesh(), DMPlexCreateBoxMesh(), DMSetType(), DMCreate()
2265 @*/
DMPlexCreateBallMesh(MPI_Comm comm,PetscInt dim,PetscReal R,DM * dm)2266 PetscErrorCode DMPlexCreateBallMesh(MPI_Comm comm, PetscInt dim, PetscReal R, DM *dm)
2267 {
2268 DM sdm;
2269 DMLabel bdlabel;
2270 PetscErrorCode ierr;
2271
2272 PetscFunctionBegin;
2273 ierr = DMPlexCreateSphereMesh(comm, dim-1, PETSC_TRUE, R, &sdm);CHKERRQ(ierr);
2274 ierr = PetscObjectSetOptionsPrefix((PetscObject) sdm, "bd_");CHKERRQ(ierr);
2275 ierr = DMSetFromOptions(sdm);CHKERRQ(ierr);
2276 ierr = DMPlexGenerate(sdm, NULL, PETSC_TRUE, dm);CHKERRQ(ierr);
2277 ierr = DMDestroy(&sdm);CHKERRQ(ierr);
2278 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr);
2279 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr);
2280 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr);
2281 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr);
2282 PetscFunctionReturn(0);
2283 }
2284
2285 /* External function declarations here */
2286 extern PetscErrorCode DMCreateInterpolation_Plex(DM dmCoarse, DM dmFine, Mat *interpolation, Vec *scaling);
2287 extern PetscErrorCode DMCreateInjection_Plex(DM dmCoarse, DM dmFine, Mat *mat);
2288 extern PetscErrorCode DMCreateMassMatrix_Plex(DM dmCoarse, DM dmFine, Mat *mat);
2289 extern PetscErrorCode DMCreateLocalSection_Plex(DM dm);
2290 extern PetscErrorCode DMCreateDefaultConstraints_Plex(DM dm);
2291 extern PetscErrorCode DMCreateMatrix_Plex(DM dm, Mat *J);
2292 extern PetscErrorCode DMCreateCoordinateDM_Plex(DM dm, DM *cdm);
2293 extern PetscErrorCode DMCreateCoordinateField_Plex(DM dm, DMField *field);
2294 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm);
2295 extern PetscErrorCode DMSetUp_Plex(DM dm);
2296 extern PetscErrorCode DMDestroy_Plex(DM dm);
2297 extern PetscErrorCode DMView_Plex(DM dm, PetscViewer viewer);
2298 extern PetscErrorCode DMLoad_Plex(DM dm, PetscViewer viewer);
2299 extern PetscErrorCode DMCreateSubDM_Plex(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm);
2300 extern PetscErrorCode DMCreateSuperDM_Plex(DM dms[], PetscInt len, IS **is, DM *superdm);
2301 static PetscErrorCode DMInitialize_Plex(DM dm);
2302
2303 /* Replace dm with the contents of dmNew
2304 - Share the DM_Plex structure
2305 - Share the coordinates
2306 - Share the SF
2307 */
DMPlexReplace_Static(DM dm,DM dmNew)2308 static PetscErrorCode DMPlexReplace_Static(DM dm, DM dmNew)
2309 {
2310 PetscSF sf;
2311 DM coordDM, coarseDM;
2312 DMField coordField;
2313 Vec coords;
2314 PetscBool isper;
2315 const PetscReal *maxCell, *L;
2316 const DMBoundaryType *bd;
2317 PetscErrorCode ierr;
2318
2319 PetscFunctionBegin;
2320 ierr = DMGetPointSF(dmNew, &sf);CHKERRQ(ierr);
2321 ierr = DMSetPointSF(dm, sf);CHKERRQ(ierr);
2322 ierr = DMGetCoordinateDM(dmNew, &coordDM);CHKERRQ(ierr);
2323 ierr = DMGetCoordinatesLocal(dmNew, &coords);CHKERRQ(ierr);
2324 ierr = DMGetCoordinateField(dmNew, &coordField);CHKERRQ(ierr);
2325 ierr = DMSetCoordinateDM(dm, coordDM);CHKERRQ(ierr);
2326 ierr = DMSetCoordinatesLocal(dm, coords);CHKERRQ(ierr);
2327 ierr = DMSetCoordinateField(dm, coordField);CHKERRQ(ierr);
2328 ierr = DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);CHKERRQ(ierr);
2329 ierr = DMSetPeriodicity(dmNew, isper, maxCell, L, bd);CHKERRQ(ierr);
2330 ierr = DMDestroy_Plex(dm);CHKERRQ(ierr);
2331 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr);
2332 dm->data = dmNew->data;
2333 ((DM_Plex *) dmNew->data)->refct++;
2334 ierr = DMDestroyLabelLinkList_Internal(dm);CHKERRQ(ierr);
2335 ierr = DMCopyLabels(dmNew, dm, PETSC_OWN_POINTER, PETSC_TRUE);CHKERRQ(ierr);
2336 ierr = DMGetCoarseDM(dmNew,&coarseDM);CHKERRQ(ierr);
2337 ierr = DMSetCoarseDM(dm,coarseDM);CHKERRQ(ierr);
2338 PetscFunctionReturn(0);
2339 }
2340
2341 /* Swap dm with the contents of dmNew
2342 - Swap the DM_Plex structure
2343 - Swap the coordinates
2344 - Swap the point PetscSF
2345 */
DMPlexSwap_Static(DM dmA,DM dmB)2346 static PetscErrorCode DMPlexSwap_Static(DM dmA, DM dmB)
2347 {
2348 DM coordDMA, coordDMB;
2349 Vec coordsA, coordsB;
2350 PetscSF sfA, sfB;
2351 void *tmp;
2352 DMLabelLink listTmp;
2353 DMLabel depthTmp;
2354 PetscInt tmpI;
2355 PetscErrorCode ierr;
2356
2357 PetscFunctionBegin;
2358 ierr = DMGetPointSF(dmA, &sfA);CHKERRQ(ierr);
2359 ierr = DMGetPointSF(dmB, &sfB);CHKERRQ(ierr);
2360 ierr = PetscObjectReference((PetscObject) sfA);CHKERRQ(ierr);
2361 ierr = DMSetPointSF(dmA, sfB);CHKERRQ(ierr);
2362 ierr = DMSetPointSF(dmB, sfA);CHKERRQ(ierr);
2363 ierr = PetscObjectDereference((PetscObject) sfA);CHKERRQ(ierr);
2364
2365 ierr = DMGetCoordinateDM(dmA, &coordDMA);CHKERRQ(ierr);
2366 ierr = DMGetCoordinateDM(dmB, &coordDMB);CHKERRQ(ierr);
2367 ierr = PetscObjectReference((PetscObject) coordDMA);CHKERRQ(ierr);
2368 ierr = DMSetCoordinateDM(dmA, coordDMB);CHKERRQ(ierr);
2369 ierr = DMSetCoordinateDM(dmB, coordDMA);CHKERRQ(ierr);
2370 ierr = PetscObjectDereference((PetscObject) coordDMA);CHKERRQ(ierr);
2371
2372 ierr = DMGetCoordinatesLocal(dmA, &coordsA);CHKERRQ(ierr);
2373 ierr = DMGetCoordinatesLocal(dmB, &coordsB);CHKERRQ(ierr);
2374 ierr = PetscObjectReference((PetscObject) coordsA);CHKERRQ(ierr);
2375 ierr = DMSetCoordinatesLocal(dmA, coordsB);CHKERRQ(ierr);
2376 ierr = DMSetCoordinatesLocal(dmB, coordsA);CHKERRQ(ierr);
2377 ierr = PetscObjectDereference((PetscObject) coordsA);CHKERRQ(ierr);
2378
2379 tmp = dmA->data;
2380 dmA->data = dmB->data;
2381 dmB->data = tmp;
2382 listTmp = dmA->labels;
2383 dmA->labels = dmB->labels;
2384 dmB->labels = listTmp;
2385 depthTmp = dmA->depthLabel;
2386 dmA->depthLabel = dmB->depthLabel;
2387 dmB->depthLabel = depthTmp;
2388 depthTmp = dmA->celltypeLabel;
2389 dmA->celltypeLabel = dmB->celltypeLabel;
2390 dmB->celltypeLabel = depthTmp;
2391 tmpI = dmA->levelup;
2392 dmA->levelup = dmB->levelup;
2393 dmB->levelup = tmpI;
2394 PetscFunctionReturn(0);
2395 }
2396
DMSetFromOptions_NonRefinement_Plex(PetscOptionItems * PetscOptionsObject,DM dm)2397 PetscErrorCode DMSetFromOptions_NonRefinement_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
2398 {
2399 DM_Plex *mesh = (DM_Plex*) dm->data;
2400 PetscBool flg;
2401 PetscErrorCode ierr;
2402
2403 PetscFunctionBegin;
2404 /* Handle viewing */
2405 ierr = PetscOptionsBool("-dm_plex_print_set_values", "Output all set values info", "DMPlexMatSetClosure", PETSC_FALSE, &mesh->printSetValues, NULL);CHKERRQ(ierr);
2406 ierr = PetscOptionsBoundedInt("-dm_plex_print_fem", "Debug output level all fem computations", "DMPlexSNESComputeResidualFEM", 0, &mesh->printFEM, NULL,0);CHKERRQ(ierr);
2407 ierr = PetscOptionsReal("-dm_plex_print_tol", "Tolerance for FEM output", "DMPlexSNESComputeResidualFEM", mesh->printTol, &mesh->printTol, NULL);CHKERRQ(ierr);
2408 ierr = PetscOptionsBoundedInt("-dm_plex_print_l2", "Debug output level all L2 diff computations", "DMComputeL2Diff", 0, &mesh->printL2, NULL,0);CHKERRQ(ierr);
2409 ierr = DMMonitorSetFromOptions(dm, "-dm_plex_monitor_throughput", "Monitor the simulation throughput", "DMPlexMonitorThroughput", DMPlexMonitorThroughput, NULL, &flg);CHKERRQ(ierr);
2410 if (flg) {ierr = PetscLogDefaultBegin();CHKERRQ(ierr);}
2411 /* Point Location */
2412 ierr = PetscOptionsBool("-dm_plex_hash_location", "Use grid hashing for point location", "DMInterpolate", PETSC_FALSE, &mesh->useHashLocation, NULL);CHKERRQ(ierr);
2413 /* Partitioning and distribution */
2414 ierr = PetscOptionsBool("-dm_plex_partition_balance", "Attempt to evenly divide points on partition boundary between processes", "DMPlexSetPartitionBalance", PETSC_FALSE, &mesh->partitionBalance, NULL);CHKERRQ(ierr);
2415 /* Generation and remeshing */
2416 ierr = PetscOptionsBool("-dm_plex_remesh_bd", "Allow changes to the boundary on remeshing", "DMAdapt", PETSC_FALSE, &mesh->remeshBd, NULL);CHKERRQ(ierr);
2417 /* Projection behavior */
2418 ierr = PetscOptionsBoundedInt("-dm_plex_max_projection_height", "Maxmimum mesh point height used to project locally", "DMPlexSetMaxProjectionHeight", 0, &mesh->maxProjectionHeight, NULL,0);CHKERRQ(ierr);
2419 ierr = PetscOptionsBool("-dm_plex_regular_refinement", "Use special nested projection algorithm for regular refinement", "DMPlexSetRegularRefinement", mesh->regularRefinement, &mesh->regularRefinement, NULL);CHKERRQ(ierr);
2420 ierr = PetscOptionsEnum("-dm_plex_cell_refiner", "Strategy for cell refinment", "ex40.c", DMPlexCellRefinerTypes, (PetscEnum) mesh->cellRefiner, (PetscEnum *) &mesh->cellRefiner, NULL);CHKERRQ(ierr);
2421 /* Checking structure */
2422 {
2423 PetscBool flg = PETSC_FALSE, flg2 = PETSC_FALSE, all = PETSC_FALSE;
2424
2425 ierr = PetscOptionsBool("-dm_plex_check_all", "Perform all checks", NULL, PETSC_FALSE, &all, &flg2);CHKERRQ(ierr);
2426 ierr = PetscOptionsBool("-dm_plex_check_symmetry", "Check that the adjacency information in the mesh is symmetric", "DMPlexCheckSymmetry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2427 if (all || (flg && flg2)) {ierr = DMPlexCheckSymmetry(dm);CHKERRQ(ierr);}
2428 ierr = PetscOptionsBool("-dm_plex_check_skeleton", "Check that each cell has the correct number of vertices (only for homogeneous simplex or tensor meshes)", "DMPlexCheckSkeleton", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2429 if (all || (flg && flg2)) {ierr = DMPlexCheckSkeleton(dm, 0);CHKERRQ(ierr);}
2430 ierr = PetscOptionsBool("-dm_plex_check_faces", "Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type", "DMPlexCheckFaces", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2431 if (all || (flg && flg2)) {ierr = DMPlexCheckFaces(dm, 0);CHKERRQ(ierr);}
2432 ierr = PetscOptionsBool("-dm_plex_check_geometry", "Check that cells have positive volume", "DMPlexCheckGeometry", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2433 if (all || (flg && flg2)) {ierr = DMPlexCheckGeometry(dm);CHKERRQ(ierr);}
2434 ierr = PetscOptionsBool("-dm_plex_check_pointsf", "Check some necessary conditions for PointSF", "DMPlexCheckPointSF", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2435 if (all || (flg && flg2)) {ierr = DMPlexCheckPointSF(dm);CHKERRQ(ierr);}
2436 ierr = PetscOptionsBool("-dm_plex_check_interface_cones", "Check points on inter-partition interfaces have conforming order of cone points", "DMPlexCheckInterfaceCones", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2437 if (all || (flg && flg2)) {ierr = DMPlexCheckInterfaceCones(dm);CHKERRQ(ierr);}
2438 ierr = PetscOptionsBool("-dm_plex_check_cell_shape", "Check cell shape", "DMPlexCheckCellShape", PETSC_FALSE, &flg, &flg2);CHKERRQ(ierr);
2439 if (flg && flg2) {ierr = DMPlexCheckCellShape(dm, PETSC_TRUE, PETSC_DETERMINE);CHKERRQ(ierr);}
2440 }
2441
2442 ierr = PetscPartitionerSetFromOptions(mesh->partitioner);CHKERRQ(ierr);
2443 PetscFunctionReturn(0);
2444 }
2445
DMSetFromOptions_Plex(PetscOptionItems * PetscOptionsObject,DM dm)2446 static PetscErrorCode DMSetFromOptions_Plex(PetscOptionItems *PetscOptionsObject,DM dm)
2447 {
2448 PetscReal volume = -1.0;
2449 PetscInt prerefine = 0, refine = 0, r, coarsen = 0, overlap = 0;
2450 PetscBool uniformOrig, uniform = PETSC_TRUE, distribute = PETSC_FALSE, isHierarchy, flg;
2451 PetscErrorCode ierr;
2452
2453 PetscFunctionBegin;
2454 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2455 ierr = PetscOptionsHead(PetscOptionsObject,"DMPlex Options");CHKERRQ(ierr);
2456 /* Handle DMPlex refinement before distribution */
2457 ierr = DMPlexGetRefinementUniform(dm, &uniformOrig);CHKERRQ(ierr);
2458 ierr = PetscOptionsBool("-dm_refine_uniform_pre", "Flag for uniform refinement before distribution", "DMCreate", uniform, &uniform, &flg);CHKERRQ(ierr);
2459 if (flg) {ierr = DMPlexSetRefinementUniform(dm, uniform);CHKERRQ(ierr);}
2460 ierr = PetscOptionsReal("-dm_refine_volume_limit_pre", "The maximum cell volume after refinement before distribution", "DMCreate", volume, &volume, &flg);CHKERRQ(ierr);
2461 if (flg) {ierr = DMPlexSetRefinementLimit(dm, volume);CHKERRQ(ierr);}
2462 ierr = PetscOptionsBoundedInt("-dm_refine_pre", "The number of refinements before distribution", "DMCreate", prerefine, &prerefine, NULL,0);CHKERRQ(ierr);
2463 for (r = 0; r < prerefine; ++r) {
2464 DM rdm;
2465 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc;
2466
2467 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2468 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &rdm);CHKERRQ(ierr);
2469 /* Total hack since we do not pass in a pointer */
2470 ierr = DMPlexReplace_Static(dm, rdm);CHKERRQ(ierr);
2471 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2472 if (coordFunc) {
2473 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr);
2474 ((DM_Plex*) dm->data)->coordFunc = coordFunc;
2475 }
2476 ierr = DMDestroy(&rdm);CHKERRQ(ierr);
2477 }
2478 ierr = DMPlexSetRefinementUniform(dm, uniformOrig);CHKERRQ(ierr);
2479 /* Handle DMPlex distribution */
2480 ierr = PetscOptionsBool("-dm_distribute", "Flag to redistribute a mesh among processes", "DMCreate", distribute, &distribute, NULL);CHKERRQ(ierr);
2481 ierr = PetscOptionsBoundedInt("-dm_distribute_overlap", "The size of the overlap halo", "DMCreate", overlap, &overlap, NULL, 0);CHKERRQ(ierr);
2482 if (distribute) {
2483 DM pdm = NULL;
2484 PetscPartitioner part;
2485
2486 ierr = DMPlexGetPartitioner(dm, &part);CHKERRQ(ierr);
2487 ierr = PetscPartitionerSetFromOptions(part);CHKERRQ(ierr);
2488 ierr = DMPlexDistribute(dm, overlap, NULL, &pdm);CHKERRQ(ierr);
2489 if (pdm) {
2490 ierr = DMPlexReplace_Static(dm, pdm);CHKERRQ(ierr);
2491 ierr = DMDestroy(&pdm);CHKERRQ(ierr);
2492 }
2493 }
2494 /* Handle DMPlex refinement */
2495 ierr = PetscOptionsBoundedInt("-dm_refine", "The number of uniform refinements", "DMCreate", refine, &refine, NULL,0);CHKERRQ(ierr);
2496 ierr = PetscOptionsBoundedInt("-dm_refine_hierarchy", "The number of uniform refinements", "DMCreate", refine, &refine, &isHierarchy,0);CHKERRQ(ierr);
2497 if (refine) {ierr = DMPlexSetRefinementUniform(dm, PETSC_TRUE);CHKERRQ(ierr);}
2498 if (refine && isHierarchy) {
2499 DM *dms, coarseDM;
2500
2501 ierr = DMGetCoarseDM(dm, &coarseDM);CHKERRQ(ierr);
2502 ierr = PetscObjectReference((PetscObject)coarseDM);CHKERRQ(ierr);
2503 ierr = PetscMalloc1(refine,&dms);CHKERRQ(ierr);
2504 ierr = DMRefineHierarchy(dm, refine, dms);CHKERRQ(ierr);
2505 /* Total hack since we do not pass in a pointer */
2506 ierr = DMPlexSwap_Static(dm, dms[refine-1]);CHKERRQ(ierr);
2507 if (refine == 1) {
2508 ierr = DMSetCoarseDM(dm, dms[0]);CHKERRQ(ierr);
2509 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr);
2510 } else {
2511 ierr = DMSetCoarseDM(dm, dms[refine-2]);CHKERRQ(ierr);
2512 ierr = DMPlexSetRegularRefinement(dm, PETSC_TRUE);CHKERRQ(ierr);
2513 ierr = DMSetCoarseDM(dms[0], dms[refine-1]);CHKERRQ(ierr);
2514 ierr = DMPlexSetRegularRefinement(dms[0], PETSC_TRUE);CHKERRQ(ierr);
2515 }
2516 ierr = DMSetCoarseDM(dms[refine-1], coarseDM);CHKERRQ(ierr);
2517 ierr = PetscObjectDereference((PetscObject)coarseDM);CHKERRQ(ierr);
2518 /* Free DMs */
2519 for (r = 0; r < refine; ++r) {
2520 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr);
2521 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr);
2522 }
2523 ierr = PetscFree(dms);CHKERRQ(ierr);
2524 } else {
2525 for (r = 0; r < refine; ++r) {
2526 DM refinedMesh;
2527 PetscPointFunc coordFunc = ((DM_Plex*) dm->data)->coordFunc;
2528
2529 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2530 ierr = DMRefine(dm, PetscObjectComm((PetscObject) dm), &refinedMesh);CHKERRQ(ierr);
2531 /* Total hack since we do not pass in a pointer */
2532 ierr = DMPlexReplace_Static(dm, refinedMesh);CHKERRQ(ierr);
2533 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2534 if (coordFunc) {
2535 ierr = DMPlexRemapGeometry(dm, 0.0, coordFunc);CHKERRQ(ierr);
2536 ((DM_Plex*) dm->data)->coordFunc = coordFunc;
2537 }
2538 ierr = DMDestroy(&refinedMesh);CHKERRQ(ierr);
2539 }
2540 }
2541 /* Handle DMPlex coarsening */
2542 ierr = PetscOptionsBoundedInt("-dm_coarsen", "Coarsen the mesh", "DMCreate", coarsen, &coarsen, NULL,0);CHKERRQ(ierr);
2543 ierr = PetscOptionsBoundedInt("-dm_coarsen_hierarchy", "The number of coarsenings", "DMCreate", coarsen, &coarsen, &isHierarchy,0);CHKERRQ(ierr);
2544 if (coarsen && isHierarchy) {
2545 DM *dms;
2546
2547 ierr = PetscMalloc1(coarsen, &dms);CHKERRQ(ierr);
2548 ierr = DMCoarsenHierarchy(dm, coarsen, dms);CHKERRQ(ierr);
2549 /* Free DMs */
2550 for (r = 0; r < coarsen; ++r) {
2551 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dms[r]);CHKERRQ(ierr);
2552 ierr = DMDestroy(&dms[r]);CHKERRQ(ierr);
2553 }
2554 ierr = PetscFree(dms);CHKERRQ(ierr);
2555 } else {
2556 for (r = 0; r < coarsen; ++r) {
2557 DM coarseMesh;
2558
2559 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2560 ierr = DMCoarsen(dm, PetscObjectComm((PetscObject) dm), &coarseMesh);CHKERRQ(ierr);
2561 /* Total hack since we do not pass in a pointer */
2562 ierr = DMPlexReplace_Static(dm, coarseMesh);CHKERRQ(ierr);
2563 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2564 ierr = DMDestroy(&coarseMesh);CHKERRQ(ierr);
2565 }
2566 }
2567 /* Handle */
2568 ierr = DMSetFromOptions_NonRefinement_Plex(PetscOptionsObject, dm);CHKERRQ(ierr);
2569 ierr = PetscOptionsTail();CHKERRQ(ierr);
2570 PetscFunctionReturn(0);
2571 }
2572
DMCreateGlobalVector_Plex(DM dm,Vec * vec)2573 static PetscErrorCode DMCreateGlobalVector_Plex(DM dm,Vec *vec)
2574 {
2575 PetscErrorCode ierr;
2576
2577 PetscFunctionBegin;
2578 ierr = DMCreateGlobalVector_Section_Private(dm,vec);CHKERRQ(ierr);
2579 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */
2580 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex);CHKERRQ(ierr);
2581 ierr = VecSetOperation(*vec, VECOP_VIEWNATIVE, (void (*)(void)) VecView_Plex_Native);CHKERRQ(ierr);
2582 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex);CHKERRQ(ierr);
2583 ierr = VecSetOperation(*vec, VECOP_LOADNATIVE, (void (*)(void)) VecLoad_Plex_Native);CHKERRQ(ierr);
2584 PetscFunctionReturn(0);
2585 }
2586
DMCreateLocalVector_Plex(DM dm,Vec * vec)2587 static PetscErrorCode DMCreateLocalVector_Plex(DM dm,Vec *vec)
2588 {
2589 PetscErrorCode ierr;
2590
2591 PetscFunctionBegin;
2592 ierr = DMCreateLocalVector_Section_Private(dm,vec);CHKERRQ(ierr);
2593 ierr = VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Plex_Local);CHKERRQ(ierr);
2594 ierr = VecSetOperation(*vec, VECOP_LOAD, (void (*)(void)) VecLoad_Plex_Local);CHKERRQ(ierr);
2595 PetscFunctionReturn(0);
2596 }
2597
DMGetDimPoints_Plex(DM dm,PetscInt dim,PetscInt * pStart,PetscInt * pEnd)2598 static PetscErrorCode DMGetDimPoints_Plex(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
2599 {
2600 PetscInt depth, d;
2601 PetscErrorCode ierr;
2602
2603 PetscFunctionBegin;
2604 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr);
2605 if (depth == 1) {
2606 ierr = DMGetDimension(dm, &d);CHKERRQ(ierr);
2607 if (dim == 0) {ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);}
2608 else if (dim == d) {ierr = DMPlexGetDepthStratum(dm, 1, pStart, pEnd);CHKERRQ(ierr);}
2609 else {*pStart = 0; *pEnd = 0;}
2610 } else {
2611 ierr = DMPlexGetDepthStratum(dm, dim, pStart, pEnd);CHKERRQ(ierr);
2612 }
2613 PetscFunctionReturn(0);
2614 }
2615
DMGetNeighbors_Plex(DM dm,PetscInt * nranks,const PetscMPIInt * ranks[])2616 static PetscErrorCode DMGetNeighbors_Plex(DM dm, PetscInt *nranks, const PetscMPIInt *ranks[])
2617 {
2618 PetscSF sf;
2619 PetscInt niranks, njranks, n;
2620 const PetscMPIInt *iranks, *jranks;
2621 DM_Plex *data = (DM_Plex*) dm->data;
2622 PetscErrorCode ierr;
2623
2624 PetscFunctionBegin;
2625 ierr = DMGetPointSF(dm, &sf);CHKERRQ(ierr);
2626 if (!data->neighbors) {
2627 ierr = PetscSFGetRootRanks(sf, &njranks, &jranks, NULL, NULL, NULL);CHKERRQ(ierr);
2628 ierr = PetscSFGetLeafRanks(sf, &niranks, &iranks, NULL, NULL);CHKERRQ(ierr);
2629 ierr = PetscMalloc1(njranks + niranks + 1, &data->neighbors);CHKERRQ(ierr);
2630 ierr = PetscArraycpy(data->neighbors + 1, jranks, njranks);CHKERRQ(ierr);
2631 ierr = PetscArraycpy(data->neighbors + njranks + 1, iranks, niranks);CHKERRQ(ierr);
2632 n = njranks + niranks;
2633 ierr = PetscSortRemoveDupsMPIInt(&n, data->neighbors + 1);CHKERRQ(ierr);
2634 /* The following cast should never fail: can't have more neighbors than PETSC_MPI_INT_MAX */
2635 ierr = PetscMPIIntCast(n, data->neighbors);CHKERRQ(ierr);
2636 }
2637 if (nranks) *nranks = data->neighbors[0];
2638 if (ranks) {
2639 if (data->neighbors[0]) *ranks = data->neighbors + 1;
2640 else *ranks = NULL;
2641 }
2642 PetscFunctionReturn(0);
2643 }
2644
DMInitialize_Plex(DM dm)2645 static PetscErrorCode DMInitialize_Plex(DM dm)
2646 {
2647 PetscErrorCode ierr;
2648
2649 PetscFunctionBegin;
2650 dm->ops->view = DMView_Plex;
2651 dm->ops->load = DMLoad_Plex;
2652 dm->ops->setfromoptions = DMSetFromOptions_Plex;
2653 dm->ops->clone = DMClone_Plex;
2654 dm->ops->setup = DMSetUp_Plex;
2655 dm->ops->createlocalsection = DMCreateLocalSection_Plex;
2656 dm->ops->createdefaultconstraints = DMCreateDefaultConstraints_Plex;
2657 dm->ops->createglobalvector = DMCreateGlobalVector_Plex;
2658 dm->ops->createlocalvector = DMCreateLocalVector_Plex;
2659 dm->ops->getlocaltoglobalmapping = NULL;
2660 dm->ops->createfieldis = NULL;
2661 dm->ops->createcoordinatedm = DMCreateCoordinateDM_Plex;
2662 dm->ops->createcoordinatefield = DMCreateCoordinateField_Plex;
2663 dm->ops->getcoloring = NULL;
2664 dm->ops->creatematrix = DMCreateMatrix_Plex;
2665 dm->ops->createinterpolation = DMCreateInterpolation_Plex;
2666 dm->ops->createmassmatrix = DMCreateMassMatrix_Plex;
2667 dm->ops->createinjection = DMCreateInjection_Plex;
2668 dm->ops->refine = DMRefine_Plex;
2669 dm->ops->coarsen = DMCoarsen_Plex;
2670 dm->ops->refinehierarchy = DMRefineHierarchy_Plex;
2671 dm->ops->coarsenhierarchy = DMCoarsenHierarchy_Plex;
2672 dm->ops->adaptlabel = DMAdaptLabel_Plex;
2673 dm->ops->adaptmetric = DMAdaptMetric_Plex;
2674 dm->ops->globaltolocalbegin = NULL;
2675 dm->ops->globaltolocalend = NULL;
2676 dm->ops->localtoglobalbegin = NULL;
2677 dm->ops->localtoglobalend = NULL;
2678 dm->ops->destroy = DMDestroy_Plex;
2679 dm->ops->createsubdm = DMCreateSubDM_Plex;
2680 dm->ops->createsuperdm = DMCreateSuperDM_Plex;
2681 dm->ops->getdimpoints = DMGetDimPoints_Plex;
2682 dm->ops->locatepoints = DMLocatePoints_Plex;
2683 dm->ops->projectfunctionlocal = DMProjectFunctionLocal_Plex;
2684 dm->ops->projectfunctionlabellocal = DMProjectFunctionLabelLocal_Plex;
2685 dm->ops->projectfieldlocal = DMProjectFieldLocal_Plex;
2686 dm->ops->projectfieldlabellocal = DMProjectFieldLabelLocal_Plex;
2687 dm->ops->projectbdfieldlabellocal = DMProjectBdFieldLabelLocal_Plex;
2688 dm->ops->computel2diff = DMComputeL2Diff_Plex;
2689 dm->ops->computel2gradientdiff = DMComputeL2GradientDiff_Plex;
2690 dm->ops->computel2fielddiff = DMComputeL2FieldDiff_Plex;
2691 dm->ops->getneighbors = DMGetNeighbors_Plex;
2692 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertBoundaryValues_C",DMPlexInsertBoundaryValues_Plex);CHKERRQ(ierr);
2693 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexInsertTimeDerviativeBoundaryValues_C",DMPlexInsertTimeDerivativeBoundaryValues_Plex);CHKERRQ(ierr);
2694 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMSetUpGLVisViewer_C",DMSetUpGLVisViewer_Plex);CHKERRQ(ierr);
2695 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMCreateNeumannOverlap_C",DMCreateNeumannOverlap_Plex);CHKERRQ(ierr);
2696 ierr = PetscObjectComposeFunction((PetscObject)dm,"DMPlexGetOverlap_C",DMPlexGetOverlap_Plex);CHKERRQ(ierr);
2697 PetscFunctionReturn(0);
2698 }
2699
DMClone_Plex(DM dm,DM * newdm)2700 PETSC_INTERN PetscErrorCode DMClone_Plex(DM dm, DM *newdm)
2701 {
2702 DM_Plex *mesh = (DM_Plex *) dm->data;
2703 PetscErrorCode ierr;
2704
2705 PetscFunctionBegin;
2706 mesh->refct++;
2707 (*newdm)->data = mesh;
2708 ierr = PetscObjectChangeTypeName((PetscObject) *newdm, DMPLEX);CHKERRQ(ierr);
2709 ierr = DMInitialize_Plex(*newdm);CHKERRQ(ierr);
2710 PetscFunctionReturn(0);
2711 }
2712
2713 /*MC
2714 DMPLEX = "plex" - A DM object that encapsulates an unstructured mesh, or CW Complex, which can be expressed using a Hasse Diagram.
2715 In the local representation, Vecs contain all unknowns in the interior and shared boundary. This is
2716 specified by a PetscSection object. Ownership in the global representation is determined by
2717 ownership of the underlying DMPlex points. This is specified by another PetscSection object.
2718
2719 Options Database Keys:
2720 + -dm_refine_pre - Refine mesh before distribution
2721 + -dm_refine_uniform_pre - Choose uniform or generator-based refinement
2722 + -dm_refine_volume_limit_pre - Cell volume limit after pre-refinement using generator
2723 . -dm_distribute - Distribute mesh across processes
2724 . -dm_distribute_overlap - Number of cells to overlap for distribution
2725 . -dm_refine - Refine mesh after distribution
2726 . -dm_plex_hash_location - Use grid hashing for point location
2727 . -dm_plex_partition_balance - Attempt to evenly divide points on partition boundary between processes
2728 . -dm_plex_remesh_bd - Allow changes to the boundary on remeshing
2729 . -dm_plex_max_projection_height - Maxmimum mesh point height used to project locally
2730 . -dm_plex_regular_refinement - Use special nested projection algorithm for regular refinement
2731 . -dm_plex_check_all - Perform all shecks below
2732 . -dm_plex_check_symmetry - Check that the adjacency information in the mesh is symmetric
2733 . -dm_plex_check_skeleton <celltype> - Check that each cell has the correct number of vertices
2734 . -dm_plex_check_faces <celltype> - Check that the faces of each cell give a vertex order this is consistent with what we expect from the cell type
2735 . -dm_plex_check_geometry - Check that cells have positive volume
2736 . -dm_view :mesh.tex:ascii_latex - View the mesh in LaTeX/TikZ
2737 . -dm_plex_view_scale <num> - Scale the TikZ
2738 - -dm_plex_print_fem <num> - View FEM assembly information, such as element vectors and matrices
2739
2740
2741 Level: intermediate
2742
2743 .seealso: DMType, DMPlexCreate(), DMCreate(), DMSetType()
2744 M*/
2745
DMCreate_Plex(DM dm)2746 PETSC_EXTERN PetscErrorCode DMCreate_Plex(DM dm)
2747 {
2748 DM_Plex *mesh;
2749 PetscInt unit;
2750 PetscErrorCode ierr;
2751
2752 PetscFunctionBegin;
2753 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
2754 ierr = PetscNewLog(dm,&mesh);CHKERRQ(ierr);
2755 dm->dim = 0;
2756 dm->data = mesh;
2757
2758 mesh->refct = 1;
2759 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->coneSection);CHKERRQ(ierr);
2760 mesh->maxConeSize = 0;
2761 mesh->cones = NULL;
2762 mesh->coneOrientations = NULL;
2763 ierr = PetscSectionCreate(PetscObjectComm((PetscObject)dm), &mesh->supportSection);CHKERRQ(ierr);
2764 mesh->maxSupportSize = 0;
2765 mesh->supports = NULL;
2766 mesh->refinementUniform = PETSC_TRUE;
2767 mesh->refinementLimit = -1.0;
2768 mesh->interpolated = DMPLEX_INTERPOLATED_INVALID;
2769 mesh->interpolatedCollective = DMPLEX_INTERPOLATED_INVALID;
2770
2771 mesh->facesTmp = NULL;
2772
2773 mesh->tetgenOpts = NULL;
2774 mesh->triangleOpts = NULL;
2775 ierr = PetscPartitionerCreate(PetscObjectComm((PetscObject)dm), &mesh->partitioner);CHKERRQ(ierr);
2776 mesh->remeshBd = PETSC_FALSE;
2777
2778 mesh->subpointMap = NULL;
2779
2780 for (unit = 0; unit < NUM_PETSC_UNITS; ++unit) mesh->scale[unit] = 1.0;
2781
2782 mesh->regularRefinement = PETSC_FALSE;
2783 mesh->depthState = -1;
2784 mesh->celltypeState = -1;
2785 mesh->globalVertexNumbers = NULL;
2786 mesh->globalCellNumbers = NULL;
2787 mesh->anchorSection = NULL;
2788 mesh->anchorIS = NULL;
2789 mesh->createanchors = NULL;
2790 mesh->computeanchormatrix = NULL;
2791 mesh->parentSection = NULL;
2792 mesh->parents = NULL;
2793 mesh->childIDs = NULL;
2794 mesh->childSection = NULL;
2795 mesh->children = NULL;
2796 mesh->referenceTree = NULL;
2797 mesh->getchildsymmetry = NULL;
2798 mesh->vtkCellHeight = 0;
2799 mesh->useAnchors = PETSC_FALSE;
2800
2801 mesh->maxProjectionHeight = 0;
2802
2803 mesh->neighbors = NULL;
2804
2805 mesh->printSetValues = PETSC_FALSE;
2806 mesh->printFEM = 0;
2807 mesh->printTol = 1.0e-10;
2808
2809 ierr = DMInitialize_Plex(dm);CHKERRQ(ierr);
2810 PetscFunctionReturn(0);
2811 }
2812
2813 /*@
2814 DMPlexCreate - Creates a DMPlex object, which encapsulates an unstructured mesh, or CW complex, which can be expressed using a Hasse Diagram.
2815
2816 Collective
2817
2818 Input Parameter:
2819 . comm - The communicator for the DMPlex object
2820
2821 Output Parameter:
2822 . mesh - The DMPlex object
2823
2824 Level: beginner
2825
2826 @*/
DMPlexCreate(MPI_Comm comm,DM * mesh)2827 PetscErrorCode DMPlexCreate(MPI_Comm comm, DM *mesh)
2828 {
2829 PetscErrorCode ierr;
2830
2831 PetscFunctionBegin;
2832 PetscValidPointer(mesh,2);
2833 ierr = DMCreate(comm, mesh);CHKERRQ(ierr);
2834 ierr = DMSetType(*mesh, DMPLEX);CHKERRQ(ierr);
2835 PetscFunctionReturn(0);
2836 }
2837
2838 /*@C
2839 DMPlexBuildFromCellListParallel - Build distributed DMPLEX topology from a list of vertices for each cell (common mesh generator output)
2840
2841 Input Parameters:
2842 + dm - The DM
2843 . numCells - The number of cells owned by this process
2844 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
2845 . NVertices - The global number of vertices, or PETSC_DECIDE
2846 . numCorners - The number of vertices for each cell
2847 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
2848
2849 Output Parameter:
2850 . vertexSF - (Optional) SF describing complete vertex ownership
2851
2852 Notes:
2853 Two triangles sharing a face
2854 $
2855 $ 2
2856 $ / | \
2857 $ / | \
2858 $ / | \
2859 $ 0 0 | 1 3
2860 $ \ | /
2861 $ \ | /
2862 $ \ | /
2863 $ 1
2864 would have input
2865 $ numCells = 2, numVertices = 4
2866 $ cells = [0 1 2 1 3 2]
2867 $
2868 which would result in the DMPlex
2869 $
2870 $ 4
2871 $ / | \
2872 $ / | \
2873 $ / | \
2874 $ 2 0 | 1 5
2875 $ \ | /
2876 $ \ | /
2877 $ \ | /
2878 $ 3
2879
2880 Vertices are implicitly numbered consecutively 0,...,NVertices.
2881 Each rank owns a chunk of numVertices consecutive vertices.
2882 If numVertices is PETSC_DECIDE, PETSc will distribute them as evenly as possible using PetscLayout.
2883 If both NVertices and numVertices are PETSC_DECIDE, NVertices is computed by PETSc as the maximum vertex index in cells + 1.
2884 If only NVertices is PETSC_DECIDE, it is computed as the sum of numVertices over all ranks.
2885
2886 The cell distribution is arbitrary non-overlapping, independent of the vertex distribution.
2887
2888 Not currently supported in Fortran.
2889
2890 Level: advanced
2891
2892 .seealso: DMPlexBuildFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildCoordinatesFromCellListParallel()
2893 @*/
DMPlexBuildFromCellListParallel(DM dm,PetscInt numCells,PetscInt numVertices,PetscInt NVertices,PetscInt numCorners,const PetscInt cells[],PetscSF * vertexSF)2894 PetscErrorCode DMPlexBuildFromCellListParallel(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, const PetscInt cells[], PetscSF *vertexSF)
2895 {
2896 PetscSF sfPoint, sfVert;
2897 PetscLayout vLayout;
2898 PetscSFNode *vertexLocal, *vertexOwner, *remoteVertex;
2899 PetscInt numVerticesAdj, *verticesAdj, numVerticesGhost = 0, *localVertex, *cones, c, p, v, g, dim;
2900 PetscMPIInt rank, size;
2901 PetscErrorCode ierr;
2902
2903 PetscFunctionBegin;
2904 PetscValidLogicalCollectiveInt(dm,NVertices,4);
2905 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr);
2906 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);CHKERRQ(ierr);
2907 ierr = MPI_Comm_size(PetscObjectComm((PetscObject) dm), &size);CHKERRQ(ierr);
2908 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
2909 /* Get/check global number of vertices */
2910 {
2911 PetscInt NVerticesInCells, i;
2912 const PetscInt len = numCells * numCorners;
2913
2914 /* NVerticesInCells = max(cells) + 1 */
2915 NVerticesInCells = PETSC_MIN_INT;
2916 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
2917 ++NVerticesInCells;
2918 ierr = MPI_Allreduce(MPI_IN_PLACE, &NVerticesInCells, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject) dm));CHKERRQ(ierr);
2919
2920 if (numVertices == PETSC_DECIDE && NVertices == PETSC_DECIDE) NVertices = NVerticesInCells;
2921 else if (NVertices != PETSC_DECIDE && NVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified global number of vertices %D must be greater than or equal to the number of vertices in cells %D",NVertices,NVerticesInCells);
2922 }
2923 /* Partition vertices */
2924 ierr = PetscLayoutCreate(PetscObjectComm((PetscObject) dm), &vLayout);CHKERRQ(ierr);
2925 ierr = PetscLayoutSetLocalSize(vLayout, numVertices);CHKERRQ(ierr);
2926 ierr = PetscLayoutSetSize(vLayout, NVertices);CHKERRQ(ierr);
2927 ierr = PetscLayoutSetBlockSize(vLayout, 1);CHKERRQ(ierr);
2928 ierr = PetscLayoutSetUp(vLayout);CHKERRQ(ierr);
2929 /* Count locally unique vertices */
2930 {
2931 PetscHSetI vhash;
2932 PetscInt off = 0;
2933
2934 ierr = PetscHSetICreate(&vhash);CHKERRQ(ierr);
2935 for (c = 0; c < numCells; ++c) {
2936 for (p = 0; p < numCorners; ++p) {
2937 ierr = PetscHSetIAdd(vhash, cells[c*numCorners+p]);CHKERRQ(ierr);
2938 }
2939 }
2940 ierr = PetscHSetIGetSize(vhash, &numVerticesAdj);CHKERRQ(ierr);
2941 ierr = PetscMalloc1(numVerticesAdj, &verticesAdj);CHKERRQ(ierr);
2942 ierr = PetscHSetIGetElems(vhash, &off, verticesAdj);CHKERRQ(ierr);
2943 ierr = PetscHSetIDestroy(&vhash);CHKERRQ(ierr);
2944 if (off != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of local vertices %D should be %D", off, numVerticesAdj);
2945 }
2946 ierr = PetscSortInt(numVerticesAdj, verticesAdj);CHKERRQ(ierr);
2947 /* Create cones */
2948 ierr = DMPlexSetChart(dm, 0, numCells+numVerticesAdj);CHKERRQ(ierr);
2949 for (c = 0; c < numCells; ++c) {ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);}
2950 ierr = DMSetUp(dm);CHKERRQ(ierr);
2951 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr);
2952 for (c = 0; c < numCells; ++c) {
2953 for (p = 0; p < numCorners; ++p) {
2954 const PetscInt gv = cells[c*numCorners+p];
2955 PetscInt lv;
2956
2957 /* Positions within verticesAdj form 0-based local vertex numbering;
2958 we need to shift it by numCells to get correct DAG points (cells go first) */
2959 ierr = PetscFindInt(gv, numVerticesAdj, verticesAdj, &lv);CHKERRQ(ierr);
2960 if (lv < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Could not find global vertex %D in local connectivity", gv);
2961 cones[c*numCorners+p] = lv+numCells;
2962 }
2963 }
2964 /* Create SF for vertices */
2965 ierr = PetscSFCreate(PetscObjectComm((PetscObject)dm), &sfVert);CHKERRQ(ierr);
2966 ierr = PetscObjectSetName((PetscObject) sfVert, "Vertex Ownership SF");CHKERRQ(ierr);
2967 ierr = PetscSFSetFromOptions(sfVert);CHKERRQ(ierr);
2968 ierr = PetscSFSetGraphLayout(sfVert, vLayout, numVerticesAdj, NULL, PETSC_OWN_POINTER, verticesAdj);CHKERRQ(ierr);
2969 ierr = PetscFree(verticesAdj);CHKERRQ(ierr);
2970 /* Build pointSF */
2971 ierr = PetscMalloc2(numVerticesAdj, &vertexLocal, numVertices, &vertexOwner);CHKERRQ(ierr);
2972 for (v = 0; v < numVerticesAdj; ++v) {vertexLocal[v].index = v+numCells; vertexLocal[v].rank = rank;}
2973 for (v = 0; v < numVertices; ++v) {vertexOwner[v].index = -1; vertexOwner[v].rank = -1;}
2974 ierr = PetscSFReduceBegin(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr);
2975 ierr = PetscSFReduceEnd(sfVert, MPIU_2INT, vertexLocal, vertexOwner, MPI_MAXLOC);CHKERRQ(ierr);
2976 for (v = 0; v < numVertices; ++v) if (vertexOwner[v].rank < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Global vertex %D was unclaimed", v + vLayout->rstart);
2977 ierr = PetscSFBcastBegin(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr);
2978 ierr = PetscSFBcastEnd(sfVert, MPIU_2INT, vertexOwner, vertexLocal);CHKERRQ(ierr);
2979 for (v = 0; v < numVerticesAdj; ++v) if (vertexLocal[v].rank != rank) ++numVerticesGhost;
2980 ierr = PetscMalloc1(numVerticesGhost, &localVertex);CHKERRQ(ierr);
2981 ierr = PetscMalloc1(numVerticesGhost, &remoteVertex);CHKERRQ(ierr);
2982 for (v = 0, g = 0; v < numVerticesAdj; ++v) {
2983 if (vertexLocal[v].rank != rank) {
2984 localVertex[g] = v+numCells;
2985 remoteVertex[g].index = vertexLocal[v].index;
2986 remoteVertex[g].rank = vertexLocal[v].rank;
2987 ++g;
2988 }
2989 }
2990 ierr = PetscFree2(vertexLocal, vertexOwner);CHKERRQ(ierr);
2991 if (g != numVerticesGhost) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid number of vertex ghosts %D should be %D", g, numVerticesGhost);
2992 ierr = DMGetPointSF(dm, &sfPoint);CHKERRQ(ierr);
2993 ierr = PetscObjectSetName((PetscObject) sfPoint, "point SF");CHKERRQ(ierr);
2994 ierr = PetscSFSetGraph(sfPoint, numCells+numVerticesAdj, numVerticesGhost, localVertex, PETSC_OWN_POINTER, remoteVertex, PETSC_OWN_POINTER);CHKERRQ(ierr);
2995 ierr = PetscLayoutDestroy(&vLayout);CHKERRQ(ierr);
2996 /* Fill in the rest of the topology structure */
2997 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
2998 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
2999 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr);
3000 if (vertexSF) *vertexSF = sfVert;
3001 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);}
3002 PetscFunctionReturn(0);
3003 }
3004
3005 /*@C
3006 DMPlexBuildCoordinatesFromCellListParallel - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)
3007
3008 Input Parameters:
3009 + dm - The DM
3010 . spaceDim - The spatial dimension used for coordinates
3011 . sfVert - SF describing complete vertex ownership
3012 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3013
3014 Level: advanced
3015
3016 Notes:
3017 Not currently supported in Fortran.
3018
3019 .seealso: DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellListParallel()
3020 @*/
DMPlexBuildCoordinatesFromCellListParallel(DM dm,PetscInt spaceDim,PetscSF sfVert,const PetscReal vertexCoords[])3021 PetscErrorCode DMPlexBuildCoordinatesFromCellListParallel(DM dm, PetscInt spaceDim, PetscSF sfVert, const PetscReal vertexCoords[])
3022 {
3023 PetscSection coordSection;
3024 Vec coordinates;
3025 PetscScalar *coords;
3026 PetscInt numVertices, numVerticesAdj, coordSize, v, vStart, vEnd;
3027 PetscErrorCode ierr;
3028
3029 PetscFunctionBegin;
3030 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr);
3031 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
3032 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first.");
3033 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr);
3034 ierr = PetscSFGetGraph(sfVert, &numVertices, &numVerticesAdj, NULL, NULL);CHKERRQ(ierr);
3035 if (vEnd - vStart != numVerticesAdj) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Supplied sfVert has wrong number of leaves = %D != %D = vEnd - vStart",numVerticesAdj,vEnd - vStart);
3036 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
3037 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
3038 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
3039 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
3040 for (v = vStart; v < vEnd; ++v) {
3041 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
3042 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
3043 }
3044 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
3045 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
3046 ierr = VecCreate(PetscObjectComm((PetscObject)dm), &coordinates);CHKERRQ(ierr);
3047 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr);
3048 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
3049 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
3050 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
3051 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
3052 {
3053 MPI_Datatype coordtype;
3054
3055 /* Need a temp buffer for coords if we have complex/single */
3056 ierr = MPI_Type_contiguous(spaceDim, MPIU_SCALAR, &coordtype);CHKERRQ(ierr);
3057 ierr = MPI_Type_commit(&coordtype);CHKERRQ(ierr);
3058 #if defined(PETSC_USE_COMPLEX)
3059 {
3060 PetscScalar *svertexCoords;
3061 PetscInt i;
3062 ierr = PetscMalloc1(numVertices*spaceDim,&svertexCoords);CHKERRQ(ierr);
3063 for (i=0; i<numVertices*spaceDim; i++) svertexCoords[i] = vertexCoords[i];
3064 ierr = PetscSFBcastBegin(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr);
3065 ierr = PetscSFBcastEnd(sfVert, coordtype, svertexCoords, coords);CHKERRQ(ierr);
3066 ierr = PetscFree(svertexCoords);CHKERRQ(ierr);
3067 }
3068 #else
3069 ierr = PetscSFBcastBegin(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr);
3070 ierr = PetscSFBcastEnd(sfVert, coordtype, vertexCoords, coords);CHKERRQ(ierr);
3071 #endif
3072 ierr = MPI_Type_free(&coordtype);CHKERRQ(ierr);
3073 }
3074 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
3075 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
3076 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
3077 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr);
3078 PetscFunctionReturn(0);
3079 }
3080
3081 /*@
3082 DMPlexCreateFromCellListParallelPetsc - Create distributed DMPLEX from a list of vertices for each cell (common mesh generator output)
3083
3084 Input Parameters:
3085 + comm - The communicator
3086 . dim - The topological dimension of the mesh
3087 . numCells - The number of cells owned by this process
3088 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
3089 . NVertices - The global number of vertices, or PETSC_DECIDE
3090 . numCorners - The number of vertices for each cell
3091 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
3092 . cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
3093 . spaceDim - The spatial dimension used for coordinates
3094 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3095
3096 Output Parameter:
3097 + dm - The DM
3098 - vertexSF - (Optional) SF describing complete vertex ownership
3099
3100 Notes:
3101 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(),
3102 DMPlexBuildFromCellListParallel(), DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellListParallel()
3103
3104 See DMPlexBuildFromCellListParallel() for an example and details about the topology-related parameters.
3105 See DMPlexBuildCoordinatesFromCellListParallel() for details about the geometry-related parameters.
3106
3107 Level: intermediate
3108
3109 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromDAG(), DMPlexCreate()
3110 @*/
DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm,PetscInt dim,PetscInt numCells,PetscInt numVertices,PetscInt NVertices,PetscInt numCorners,PetscBool interpolate,const PetscInt cells[],PetscInt spaceDim,const PetscReal vertexCoords[],PetscSF * vertexSF,DM * dm)3111 PetscErrorCode DMPlexCreateFromCellListParallelPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt NVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm)
3112 {
3113 PetscSF sfVert;
3114 PetscErrorCode ierr;
3115
3116 PetscFunctionBegin;
3117 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
3118 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
3119 PetscValidLogicalCollectiveInt(*dm, dim, 2);
3120 PetscValidLogicalCollectiveInt(*dm, spaceDim, 8);
3121 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
3122 ierr = DMPlexBuildFromCellListParallel(*dm, numCells, numVertices, NVertices, numCorners, cells, &sfVert);CHKERRQ(ierr);
3123 if (interpolate) {
3124 DM idm;
3125
3126 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
3127 ierr = DMDestroy(dm);CHKERRQ(ierr);
3128 *dm = idm;
3129 }
3130 ierr = DMPlexBuildCoordinatesFromCellListParallel(*dm, spaceDim, sfVert, vertexCoords);CHKERRQ(ierr);
3131 if (vertexSF) *vertexSF = sfVert;
3132 else {ierr = PetscSFDestroy(&sfVert);CHKERRQ(ierr);}
3133 PetscFunctionReturn(0);
3134 }
3135
3136
3137 /*@
3138 DMPlexCreateFromCellListParallel - Deprecated, use DMPlexCreateFromCellListParallelPetsc()
3139
3140 Level: deprecated
3141
3142 .seealso: DMPlexCreateFromCellListParallelPetsc()
3143 @*/
DMPlexCreateFromCellListParallel(MPI_Comm comm,PetscInt dim,PetscInt numCells,PetscInt numVertices,PetscInt numCorners,PetscBool interpolate,const int cells[],PetscInt spaceDim,const PetscReal vertexCoords[],PetscSF * vertexSF,DM * dm)3144 PetscErrorCode DMPlexCreateFromCellListParallel(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const PetscReal vertexCoords[], PetscSF *vertexSF, DM *dm)
3145 {
3146 PetscErrorCode ierr;
3147 PetscInt i;
3148 PetscInt *pintCells;
3149
3150 PetscFunctionBegin;
3151 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt));
3152 if (sizeof(int) == sizeof(PetscInt)) {
3153 pintCells = (PetscInt *) cells;
3154 } else {
3155 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr);
3156 for (i = 0; i < numCells*numCorners; i++) {
3157 pintCells[i] = (PetscInt) cells[i];
3158 }
3159 }
3160 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCells, numVertices, PETSC_DECIDE, numCorners, interpolate, pintCells, spaceDim, vertexCoords, vertexSF, dm);CHKERRQ(ierr);
3161 if (sizeof(int) != sizeof(PetscInt)) {
3162 ierr = PetscFree(pintCells);CHKERRQ(ierr);
3163 }
3164 PetscFunctionReturn(0);
3165 }
3166
3167 /*@C
3168 DMPlexBuildFromCellList - Build DMPLEX topology from a list of vertices for each cell (common mesh generator output)
3169
3170 Input Parameters:
3171 + dm - The DM
3172 . numCells - The number of cells owned by this process
3173 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
3174 . numCorners - The number of vertices for each cell
3175 - cells - An array of numCells*numCorners numbers, the global vertex numbers for each cell
3176
3177 Level: advanced
3178
3179 Notes:
3180 Two triangles sharing a face
3181 $
3182 $ 2
3183 $ / | \
3184 $ / | \
3185 $ / | \
3186 $ 0 0 | 1 3
3187 $ \ | /
3188 $ \ | /
3189 $ \ | /
3190 $ 1
3191 would have input
3192 $ numCells = 2, numVertices = 4
3193 $ cells = [0 1 2 1 3 2]
3194 $
3195 which would result in the DMPlex
3196 $
3197 $ 4
3198 $ / | \
3199 $ / | \
3200 $ / | \
3201 $ 2 0 | 1 5
3202 $ \ | /
3203 $ \ | /
3204 $ \ | /
3205 $ 3
3206
3207 If numVertices is PETSC_DECIDE, it is computed by PETSc as the maximum vertex index in cells + 1.
3208
3209 Not currently supported in Fortran.
3210
3211 .seealso: DMPlexBuildFromCellListParallel(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromCellListPetsc()
3212 @*/
DMPlexBuildFromCellList(DM dm,PetscInt numCells,PetscInt numVertices,PetscInt numCorners,const PetscInt cells[])3213 PetscErrorCode DMPlexBuildFromCellList(DM dm, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, const PetscInt cells[])
3214 {
3215 PetscInt *cones, c, p, dim;
3216 PetscErrorCode ierr;
3217
3218 PetscFunctionBegin;
3219 ierr = PetscLogEventBegin(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr);
3220 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
3221 /* Get/check global number of vertices */
3222 {
3223 PetscInt NVerticesInCells, i;
3224 const PetscInt len = numCells * numCorners;
3225
3226 /* NVerticesInCells = max(cells) + 1 */
3227 NVerticesInCells = PETSC_MIN_INT;
3228 for (i=0; i<len; i++) if (cells[i] > NVerticesInCells) NVerticesInCells = cells[i];
3229 ++NVerticesInCells;
3230
3231 if (numVertices == PETSC_DECIDE) numVertices = NVerticesInCells;
3232 else if (numVertices < NVerticesInCells) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONG, "Specified number of vertices %D must be greater than or equal to the number of vertices in cells %D",numVertices,NVerticesInCells);
3233 }
3234 ierr = DMPlexSetChart(dm, 0, numCells+numVertices);CHKERRQ(ierr);
3235 for (c = 0; c < numCells; ++c) {
3236 ierr = DMPlexSetConeSize(dm, c, numCorners);CHKERRQ(ierr);
3237 }
3238 ierr = DMSetUp(dm);CHKERRQ(ierr);
3239 ierr = DMPlexGetCones(dm,&cones);CHKERRQ(ierr);
3240 for (c = 0; c < numCells; ++c) {
3241 for (p = 0; p < numCorners; ++p) {
3242 cones[c*numCorners+p] = cells[c*numCorners+p]+numCells;
3243 }
3244 }
3245 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
3246 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
3247 ierr = PetscLogEventEnd(DMPLEX_BuildFromCellList,dm,0,0,0);CHKERRQ(ierr);
3248 PetscFunctionReturn(0);
3249 }
3250
3251 /*@C
3252 DMPlexBuildCoordinatesFromCellList - Build DM coordinates from a list of coordinates for each owned vertex (common mesh generator output)
3253
3254 Input Parameters:
3255 + dm - The DM
3256 . spaceDim - The spatial dimension used for coordinates
3257 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3258
3259 Level: advanced
3260
3261 Notes:
3262 Not currently supported in Fortran.
3263
3264 .seealso: DMPlexBuildCoordinatesFromCellListParallel(), DMPlexCreateFromCellListPetsc(), DMPlexBuildFromCellList()
3265 @*/
DMPlexBuildCoordinatesFromCellList(DM dm,PetscInt spaceDim,const PetscReal vertexCoords[])3266 PetscErrorCode DMPlexBuildCoordinatesFromCellList(DM dm, PetscInt spaceDim, const PetscReal vertexCoords[])
3267 {
3268 PetscSection coordSection;
3269 Vec coordinates;
3270 DM cdm;
3271 PetscScalar *coords;
3272 PetscInt v, vStart, vEnd, d;
3273 PetscErrorCode ierr;
3274
3275 PetscFunctionBegin;
3276 ierr = PetscLogEventBegin(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr);
3277 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr);
3278 if (vStart < 0 || vEnd < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "DM is not set up properly. DMPlexBuildFromCellList() should be called first.");
3279 ierr = DMSetCoordinateDim(dm, spaceDim);CHKERRQ(ierr);
3280 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
3281 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
3282 ierr = PetscSectionSetFieldComponents(coordSection, 0, spaceDim);CHKERRQ(ierr);
3283 ierr = PetscSectionSetChart(coordSection, vStart, vEnd);CHKERRQ(ierr);
3284 for (v = vStart; v < vEnd; ++v) {
3285 ierr = PetscSectionSetDof(coordSection, v, spaceDim);CHKERRQ(ierr);
3286 ierr = PetscSectionSetFieldDof(coordSection, v, 0, spaceDim);CHKERRQ(ierr);
3287 }
3288 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
3289
3290 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr);
3291 ierr = DMCreateLocalVector(cdm, &coordinates);CHKERRQ(ierr);
3292 ierr = VecSetBlockSize(coordinates, spaceDim);CHKERRQ(ierr);
3293 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
3294 ierr = VecGetArrayWrite(coordinates, &coords);CHKERRQ(ierr);
3295 for (v = 0; v < vEnd-vStart; ++v) {
3296 for (d = 0; d < spaceDim; ++d) {
3297 coords[v*spaceDim+d] = vertexCoords[v*spaceDim+d];
3298 }
3299 }
3300 ierr = VecRestoreArrayWrite(coordinates, &coords);CHKERRQ(ierr);
3301 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
3302 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
3303 ierr = PetscLogEventEnd(DMPLEX_BuildCoordinatesFromCellList,dm,0,0,0);CHKERRQ(ierr);
3304 PetscFunctionReturn(0);
3305 }
3306
3307 /*@
3308 DMPlexCreateFromCellListPetsc - Create DMPLEX from a list of vertices for each cell (common mesh generator output)
3309
3310 Input Parameters:
3311 + comm - The communicator
3312 . dim - The topological dimension of the mesh
3313 . numCells - The number of cells
3314 . numVertices - The number of vertices owned by this process, or PETSC_DECIDE
3315 . numCorners - The number of vertices for each cell
3316 . interpolate - Flag indicating that intermediate mesh entities (faces, edges) should be created automatically
3317 . cells - An array of numCells*numCorners numbers, the vertices for each cell
3318 . spaceDim - The spatial dimension used for coordinates
3319 - vertexCoords - An array of numVertices*spaceDim numbers, the coordinates of each vertex
3320
3321 Output Parameter:
3322 . dm - The DM
3323
3324 Notes:
3325 This function is just a convenient sequence of DMCreate(), DMSetType(), DMSetDimension(), DMPlexBuildFromCellList(),
3326 DMPlexInterpolate(), DMPlexBuildCoordinatesFromCellList()
3327
3328 See DMPlexBuildFromCellList() for an example and details about the topology-related parameters.
3329 See DMPlexBuildCoordinatesFromCellList() for details about the geometry-related parameters.
3330
3331 Level: intermediate
3332
3333 .seealso: DMPlexCreateFromCellListParallelPetsc(), DMPlexBuildFromCellList(), DMPlexBuildCoordinatesFromCellList(), DMPlexCreateFromDAG(), DMPlexCreate()
3334 @*/
DMPlexCreateFromCellListPetsc(MPI_Comm comm,PetscInt dim,PetscInt numCells,PetscInt numVertices,PetscInt numCorners,PetscBool interpolate,const PetscInt cells[],PetscInt spaceDim,const PetscReal vertexCoords[],DM * dm)3335 PetscErrorCode DMPlexCreateFromCellListPetsc(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const PetscInt cells[], PetscInt spaceDim, const PetscReal vertexCoords[], DM *dm)
3336 {
3337 PetscErrorCode ierr;
3338
3339 PetscFunctionBegin;
3340 if (!dim) SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "This is not appropriate for 0-dimensional meshes. Consider either creating the DM using DMPlexCreateFromDAG(), by hand, or using DMSwarm.");
3341 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
3342 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
3343 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
3344 ierr = DMPlexBuildFromCellList(*dm, numCells, numVertices, numCorners, cells);CHKERRQ(ierr);
3345 if (interpolate) {
3346 DM idm;
3347
3348 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
3349 ierr = DMDestroy(dm);CHKERRQ(ierr);
3350 *dm = idm;
3351 }
3352 ierr = DMPlexBuildCoordinatesFromCellList(*dm, spaceDim, vertexCoords);CHKERRQ(ierr);
3353 PetscFunctionReturn(0);
3354 }
3355
3356 /*@
3357 DMPlexCreateFromCellList - Deprecated, use DMPlexCreateFromCellListPetsc()
3358
3359 Level: deprecated
3360
3361 .seealso: DMPlexCreateFromCellListPetsc()
3362 @*/
DMPlexCreateFromCellList(MPI_Comm comm,PetscInt dim,PetscInt numCells,PetscInt numVertices,PetscInt numCorners,PetscBool interpolate,const int cells[],PetscInt spaceDim,const double vertexCoords[],DM * dm)3363 PetscErrorCode DMPlexCreateFromCellList(MPI_Comm comm, PetscInt dim, PetscInt numCells, PetscInt numVertices, PetscInt numCorners, PetscBool interpolate, const int cells[], PetscInt spaceDim, const double vertexCoords[], DM *dm)
3364 {
3365 PetscErrorCode ierr;
3366 PetscInt i;
3367 PetscInt *pintCells;
3368 PetscReal *prealVC;
3369
3370 PetscFunctionBegin;
3371 if (sizeof(int) > sizeof(PetscInt)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of int %zd greater than size of PetscInt %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(int), sizeof(PetscInt));
3372 if (sizeof(int) == sizeof(PetscInt)) {
3373 pintCells = (PetscInt *) cells;
3374 } else {
3375 ierr = PetscMalloc1(numCells*numCorners, &pintCells);CHKERRQ(ierr);
3376 for (i = 0; i < numCells*numCorners; i++) {
3377 pintCells[i] = (PetscInt) cells[i];
3378 }
3379 }
3380 if (sizeof(double) > sizeof(PetscReal)) SETERRQ2(comm, PETSC_ERR_ARG_SIZ, "Size of double %zd greater than size of PetscReal %zd. Reconfigure PETSc --with-precision=<higher precision>.", sizeof(double), sizeof(PetscReal));
3381 if (sizeof(double) == sizeof(PetscReal)) {
3382 prealVC = (PetscReal *) vertexCoords;
3383 } else {
3384 ierr = PetscMalloc1(numVertices*spaceDim, &prealVC);CHKERRQ(ierr);
3385 for (i = 0; i < numVertices*spaceDim; i++) {
3386 prealVC[i] = (PetscReal) vertexCoords[i];
3387 }
3388 }
3389 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, pintCells, spaceDim, prealVC, dm);CHKERRQ(ierr);
3390 if (sizeof(int) != sizeof(PetscInt)) {
3391 ierr = PetscFree(pintCells);CHKERRQ(ierr);
3392 }
3393 if (sizeof(double) != sizeof(PetscReal)) {
3394 ierr = PetscFree(prealVC);CHKERRQ(ierr);
3395 }
3396 PetscFunctionReturn(0);
3397 }
3398
3399 /*@
3400 DMPlexCreateFromDAG - This takes as input the adjacency-list representation of the Directed Acyclic Graph (Hasse Diagram) encoding a mesh, and produces a DM
3401
3402 Input Parameters:
3403 + dm - The empty DM object, usually from DMCreate() and DMSetDimension()
3404 . depth - The depth of the DAG
3405 . numPoints - Array of size depth + 1 containing the number of points at each depth
3406 . coneSize - The cone size of each point
3407 . cones - The concatenation of the cone points for each point, the cone list must be oriented correctly for each point
3408 . coneOrientations - The orientation of each cone point
3409 - vertexCoords - An array of numPoints[0]*spacedim numbers representing the coordinates of each vertex, with spacedim the value set via DMSetCoordinateDim()
3410
3411 Output Parameter:
3412 . dm - The DM
3413
3414 Note: Two triangles sharing a face would have input
3415 $ depth = 1, numPoints = [4 2], coneSize = [3 3 0 0 0 0]
3416 $ cones = [2 3 4 3 5 4], coneOrientations = [0 0 0 0 0 0]
3417 $ vertexCoords = [-1.0 0.0 0.0 -1.0 0.0 1.0 1.0 0.0]
3418 $
3419 which would result in the DMPlex
3420 $
3421 $ 4
3422 $ / | \
3423 $ / | \
3424 $ / | \
3425 $ 2 0 | 1 5
3426 $ \ | /
3427 $ \ | /
3428 $ \ | /
3429 $ 3
3430 $
3431 $ Notice that all points are numbered consecutively, unlike DMPlexCreateFromCellListPetsc()
3432
3433 Level: advanced
3434
3435 .seealso: DMPlexCreateFromCellListPetsc(), DMPlexCreate()
3436 @*/
DMPlexCreateFromDAG(DM dm,PetscInt depth,const PetscInt numPoints[],const PetscInt coneSize[],const PetscInt cones[],const PetscInt coneOrientations[],const PetscScalar vertexCoords[])3437 PetscErrorCode DMPlexCreateFromDAG(DM dm, PetscInt depth, const PetscInt numPoints[], const PetscInt coneSize[], const PetscInt cones[], const PetscInt coneOrientations[], const PetscScalar vertexCoords[])
3438 {
3439 Vec coordinates;
3440 PetscSection coordSection;
3441 PetscScalar *coords;
3442 PetscInt coordSize, firstVertex = -1, pStart = 0, pEnd = 0, p, v, dim, dimEmbed, d, off;
3443 PetscErrorCode ierr;
3444
3445 PetscFunctionBegin;
3446 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
3447 ierr = DMGetCoordinateDim(dm, &dimEmbed);CHKERRQ(ierr);
3448 if (dimEmbed < dim) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Embedding dimension %D cannot be less than intrinsic dimension %d",dimEmbed,dim);
3449 for (d = 0; d <= depth; ++d) pEnd += numPoints[d];
3450 ierr = DMPlexSetChart(dm, pStart, pEnd);CHKERRQ(ierr);
3451 for (p = pStart; p < pEnd; ++p) {
3452 ierr = DMPlexSetConeSize(dm, p, coneSize[p-pStart]);CHKERRQ(ierr);
3453 if (firstVertex < 0 && !coneSize[p - pStart]) {
3454 firstVertex = p - pStart;
3455 }
3456 }
3457 if (firstVertex < 0 && numPoints[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Expected %D vertices but could not find any", numPoints[0]);
3458 ierr = DMSetUp(dm);CHKERRQ(ierr); /* Allocate space for cones */
3459 for (p = pStart, off = 0; p < pEnd; off += coneSize[p-pStart], ++p) {
3460 ierr = DMPlexSetCone(dm, p, &cones[off]);CHKERRQ(ierr);
3461 ierr = DMPlexSetConeOrientation(dm, p, &coneOrientations[off]);CHKERRQ(ierr);
3462 }
3463 ierr = DMPlexSymmetrize(dm);CHKERRQ(ierr);
3464 ierr = DMPlexStratify(dm);CHKERRQ(ierr);
3465 /* Build coordinates */
3466 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr);
3467 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
3468 ierr = PetscSectionSetFieldComponents(coordSection, 0, dimEmbed);CHKERRQ(ierr);
3469 ierr = PetscSectionSetChart(coordSection, firstVertex, firstVertex+numPoints[0]);CHKERRQ(ierr);
3470 for (v = firstVertex; v < firstVertex+numPoints[0]; ++v) {
3471 ierr = PetscSectionSetDof(coordSection, v, dimEmbed);CHKERRQ(ierr);
3472 ierr = PetscSectionSetFieldDof(coordSection, v, 0, dimEmbed);CHKERRQ(ierr);
3473 }
3474 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
3475 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
3476 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
3477 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
3478 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
3479 ierr = VecSetBlockSize(coordinates, dimEmbed);CHKERRQ(ierr);
3480 ierr = VecSetType(coordinates,VECSTANDARD);CHKERRQ(ierr);
3481 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
3482 for (v = 0; v < numPoints[0]; ++v) {
3483 PetscInt off;
3484
3485 ierr = PetscSectionGetOffset(coordSection, v+firstVertex, &off);CHKERRQ(ierr);
3486 for (d = 0; d < dimEmbed; ++d) {
3487 coords[off+d] = vertexCoords[v*dimEmbed+d];
3488 }
3489 }
3490 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
3491 ierr = DMSetCoordinatesLocal(dm, coordinates);CHKERRQ(ierr);
3492 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
3493 PetscFunctionReturn(0);
3494 }
3495
3496 /*@C
3497 DMPlexCreateCellVertexFromFile - Create a DMPlex mesh from a simple cell-vertex file.
3498
3499 + comm - The MPI communicator
3500 . filename - Name of the .dat file
3501 - interpolate - Create faces and edges in the mesh
3502
3503 Output Parameter:
3504 . dm - The DM object representing the mesh
3505
3506 Note: The format is the simplest possible:
3507 $ Ne
3508 $ v0 v1 ... vk
3509 $ Nv
3510 $ x y z marker
3511
3512 Level: beginner
3513
3514 .seealso: DMPlexCreateFromFile(), DMPlexCreateMedFromFile(), DMPlexCreateGmsh(), DMPlexCreate()
3515 @*/
DMPlexCreateCellVertexFromFile(MPI_Comm comm,const char filename[],PetscBool interpolate,DM * dm)3516 PetscErrorCode DMPlexCreateCellVertexFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
3517 {
3518 DMLabel marker;
3519 PetscViewer viewer;
3520 Vec coordinates;
3521 PetscSection coordSection;
3522 PetscScalar *coords;
3523 char line[PETSC_MAX_PATH_LEN];
3524 PetscInt dim = 3, cdim = 3, coordSize, v, c, d;
3525 PetscMPIInt rank;
3526 int snum, Nv, Nc;
3527 PetscErrorCode ierr;
3528
3529 PetscFunctionBegin;
3530 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
3531 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
3532 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr);
3533 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
3534 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
3535 if (!rank) {
3536 ierr = PetscViewerRead(viewer, line, 2, NULL, PETSC_STRING);CHKERRQ(ierr);
3537 snum = sscanf(line, "%d %d", &Nc, &Nv);
3538 if (snum != 2) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3539 } else {
3540 Nc = Nv = 0;
3541 }
3542 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
3543 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
3544 ierr = DMPlexSetChart(*dm, 0, Nc+Nv);CHKERRQ(ierr);
3545 ierr = DMSetDimension(*dm, dim);CHKERRQ(ierr);
3546 ierr = DMSetCoordinateDim(*dm, cdim);CHKERRQ(ierr);
3547 /* Read topology */
3548 if (!rank) {
3549 PetscInt cone[8], corners = 8;
3550 int vbuf[8], v;
3551
3552 for (c = 0; c < Nc; ++c) {ierr = DMPlexSetConeSize(*dm, c, corners);CHKERRQ(ierr);}
3553 ierr = DMSetUp(*dm);CHKERRQ(ierr);
3554 for (c = 0; c < Nc; ++c) {
3555 ierr = PetscViewerRead(viewer, line, corners, NULL, PETSC_STRING);CHKERRQ(ierr);
3556 snum = sscanf(line, "%d %d %d %d %d %d %d %d", &vbuf[0], &vbuf[1], &vbuf[2], &vbuf[3], &vbuf[4], &vbuf[5], &vbuf[6], &vbuf[7]);
3557 if (snum != corners) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3558 for (v = 0; v < corners; ++v) cone[v] = vbuf[v] + Nc;
3559 /* Hexahedra are inverted */
3560 {
3561 PetscInt tmp = cone[1];
3562 cone[1] = cone[3];
3563 cone[3] = tmp;
3564 }
3565 ierr = DMPlexSetCone(*dm, c, cone);CHKERRQ(ierr);
3566 }
3567 }
3568 ierr = DMPlexSymmetrize(*dm);CHKERRQ(ierr);
3569 ierr = DMPlexStratify(*dm);CHKERRQ(ierr);
3570 /* Read coordinates */
3571 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr);
3572 ierr = PetscSectionSetNumFields(coordSection, 1);CHKERRQ(ierr);
3573 ierr = PetscSectionSetFieldComponents(coordSection, 0, cdim);CHKERRQ(ierr);
3574 ierr = PetscSectionSetChart(coordSection, Nc, Nc + Nv);CHKERRQ(ierr);
3575 for (v = Nc; v < Nc+Nv; ++v) {
3576 ierr = PetscSectionSetDof(coordSection, v, cdim);CHKERRQ(ierr);
3577 ierr = PetscSectionSetFieldDof(coordSection, v, 0, cdim);CHKERRQ(ierr);
3578 }
3579 ierr = PetscSectionSetUp(coordSection);CHKERRQ(ierr);
3580 ierr = PetscSectionGetStorageSize(coordSection, &coordSize);CHKERRQ(ierr);
3581 ierr = VecCreate(PETSC_COMM_SELF, &coordinates);CHKERRQ(ierr);
3582 ierr = PetscObjectSetName((PetscObject) coordinates, "coordinates");CHKERRQ(ierr);
3583 ierr = VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);CHKERRQ(ierr);
3584 ierr = VecSetBlockSize(coordinates, cdim);CHKERRQ(ierr);
3585 ierr = VecSetType(coordinates, VECSTANDARD);CHKERRQ(ierr);
3586 ierr = VecGetArray(coordinates, &coords);CHKERRQ(ierr);
3587 if (!rank) {
3588 double x[3];
3589 int val;
3590
3591 ierr = DMCreateLabel(*dm, "marker");CHKERRQ(ierr);
3592 ierr = DMGetLabel(*dm, "marker", &marker);CHKERRQ(ierr);
3593 for (v = 0; v < Nv; ++v) {
3594 ierr = PetscViewerRead(viewer, line, 4, NULL, PETSC_STRING);CHKERRQ(ierr);
3595 snum = sscanf(line, "%lg %lg %lg %d", &x[0], &x[1], &x[2], &val);
3596 if (snum != 4) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unable to parse cell-vertex file: %s", line);
3597 for (d = 0; d < cdim; ++d) coords[v*cdim+d] = x[d];
3598 if (val) {ierr = DMLabelSetValue(marker, v+Nc, val);CHKERRQ(ierr);}
3599 }
3600 }
3601 ierr = VecRestoreArray(coordinates, &coords);CHKERRQ(ierr);
3602 ierr = DMSetCoordinatesLocal(*dm, coordinates);CHKERRQ(ierr);
3603 ierr = VecDestroy(&coordinates);CHKERRQ(ierr);
3604 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
3605 if (interpolate) {
3606 DM idm;
3607 DMLabel bdlabel;
3608
3609 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
3610 ierr = DMDestroy(dm);CHKERRQ(ierr);
3611 *dm = idm;
3612
3613 ierr = DMGetLabel(*dm, "marker", &bdlabel);CHKERRQ(ierr);
3614 ierr = DMPlexMarkBoundaryFaces(*dm, PETSC_DETERMINE, bdlabel);CHKERRQ(ierr);
3615 ierr = DMPlexLabelComplete(*dm, bdlabel);CHKERRQ(ierr);
3616 }
3617 PetscFunctionReturn(0);
3618 }
3619
3620 /*@C
3621 DMPlexCreateFromFile - This takes a filename and produces a DM
3622
3623 Input Parameters:
3624 + comm - The communicator
3625 . filename - A file name
3626 - interpolate - Flag to create intermediate mesh pieces (edges, faces)
3627
3628 Output Parameter:
3629 . dm - The DM
3630
3631 Options Database Keys:
3632 . -dm_plex_create_from_hdf5_xdmf - use the PETSC_VIEWER_HDF5_XDMF format for reading HDF5
3633
3634 Use -dm_plex_create_ prefix to pass options to the internal PetscViewer, e.g.
3635 $ -dm_plex_create_viewer_hdf5_collective
3636
3637 Level: beginner
3638
3639 .seealso: DMPlexCreateFromDAG(), DMPlexCreateFromCellListPetsc(), DMPlexCreate()
3640 @*/
DMPlexCreateFromFile(MPI_Comm comm,const char filename[],PetscBool interpolate,DM * dm)3641 PetscErrorCode DMPlexCreateFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
3642 {
3643 const char *extGmsh = ".msh";
3644 const char *extGmsh2 = ".msh2";
3645 const char *extGmsh4 = ".msh4";
3646 const char *extCGNS = ".cgns";
3647 const char *extExodus = ".exo";
3648 const char *extGenesis = ".gen";
3649 const char *extFluent = ".cas";
3650 const char *extHDF5 = ".h5";
3651 const char *extMed = ".med";
3652 const char *extPLY = ".ply";
3653 const char *extCV = ".dat";
3654 size_t len;
3655 PetscBool isGmsh, isGmsh2, isGmsh4, isCGNS, isExodus, isGenesis, isFluent, isHDF5, isMed, isPLY, isCV;
3656 PetscMPIInt rank;
3657 PetscErrorCode ierr;
3658
3659 PetscFunctionBegin;
3660 PetscValidCharPointer(filename, 2);
3661 PetscValidPointer(dm, 4);
3662 ierr = DMInitializePackage();CHKERRQ(ierr);
3663 ierr = PetscLogEventBegin(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr);
3664 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
3665 ierr = PetscStrlen(filename, &len);CHKERRQ(ierr);
3666 if (!len) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Filename must be a valid path");
3667 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGmsh, 4, &isGmsh);CHKERRQ(ierr);
3668 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh2, 5, &isGmsh2);CHKERRQ(ierr);
3669 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extGmsh4, 5, &isGmsh4);CHKERRQ(ierr);
3670 ierr = PetscStrncmp(&filename[PetscMax(0,len-5)], extCGNS, 5, &isCGNS);CHKERRQ(ierr);
3671 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extExodus, 4, &isExodus);CHKERRQ(ierr);
3672 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extGenesis, 4, &isGenesis);CHKERRQ(ierr);
3673 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extFluent, 4, &isFluent);CHKERRQ(ierr);
3674 ierr = PetscStrncmp(&filename[PetscMax(0,len-3)], extHDF5, 3, &isHDF5);CHKERRQ(ierr);
3675 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extMed, 4, &isMed);CHKERRQ(ierr);
3676 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extPLY, 4, &isPLY);CHKERRQ(ierr);
3677 ierr = PetscStrncmp(&filename[PetscMax(0,len-4)], extCV, 4, &isCV);CHKERRQ(ierr);
3678 if (isGmsh || isGmsh2 || isGmsh4) {
3679 ierr = DMPlexCreateGmshFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3680 } else if (isCGNS) {
3681 ierr = DMPlexCreateCGNSFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3682 } else if (isExodus || isGenesis) {
3683 ierr = DMPlexCreateExodusFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3684 } else if (isFluent) {
3685 ierr = DMPlexCreateFluentFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3686 } else if (isHDF5) {
3687 PetscBool load_hdf5_xdmf = PETSC_FALSE;
3688 PetscViewer viewer;
3689
3690 /* PETSC_VIEWER_HDF5_XDMF is used if the filename ends with .xdmf.h5, or if -dm_plex_create_from_hdf5_xdmf option is present */
3691 ierr = PetscStrncmp(&filename[PetscMax(0,len-8)], ".xdmf", 5, &load_hdf5_xdmf);CHKERRQ(ierr);
3692 ierr = PetscOptionsGetBool(NULL, NULL, "-dm_plex_create_from_hdf5_xdmf", &load_hdf5_xdmf, NULL);CHKERRQ(ierr);
3693 ierr = PetscViewerCreate(comm, &viewer);CHKERRQ(ierr);
3694 ierr = PetscViewerSetType(viewer, PETSCVIEWERHDF5);CHKERRQ(ierr);
3695 ierr = PetscViewerSetOptionsPrefix(viewer, "dm_plex_create_");CHKERRQ(ierr);
3696 ierr = PetscViewerSetFromOptions(viewer);CHKERRQ(ierr);
3697 ierr = PetscViewerFileSetMode(viewer, FILE_MODE_READ);CHKERRQ(ierr);
3698 ierr = PetscViewerFileSetName(viewer, filename);CHKERRQ(ierr);
3699 ierr = DMCreate(comm, dm);CHKERRQ(ierr);
3700 ierr = DMSetType(*dm, DMPLEX);CHKERRQ(ierr);
3701 if (load_hdf5_xdmf) {ierr = PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_XDMF);CHKERRQ(ierr);}
3702 ierr = DMLoad(*dm, viewer);CHKERRQ(ierr);
3703 if (load_hdf5_xdmf) {ierr = PetscViewerPopFormat(viewer);CHKERRQ(ierr);}
3704 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
3705
3706 if (interpolate) {
3707 DM idm;
3708
3709 ierr = DMPlexInterpolate(*dm, &idm);CHKERRQ(ierr);
3710 ierr = DMDestroy(dm);CHKERRQ(ierr);
3711 *dm = idm;
3712 }
3713 } else if (isMed) {
3714 ierr = DMPlexCreateMedFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3715 } else if (isPLY) {
3716 ierr = DMPlexCreatePLYFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3717 } else if (isCV) {
3718 ierr = DMPlexCreateCellVertexFromFile(comm, filename, interpolate, dm);CHKERRQ(ierr);
3719 } else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot load file %s: unrecognized extension", filename);
3720 ierr = PetscLogEventEnd(DMPLEX_CreateFromFile,0,0,0,0);CHKERRQ(ierr);
3721 PetscFunctionReturn(0);
3722 }
3723
3724 /*@
3725 DMPlexCreateReferenceCellByType - Create a DMPLEX with the appropriate FEM reference cell
3726
3727 Collective
3728
3729 Input Parameters:
3730 + comm - The communicator
3731 - ct - The cell type of the reference cell
3732
3733 Output Parameter:
3734 . refdm - The reference cell
3735
3736 Level: intermediate
3737
3738 .seealso: DMPlexCreateReferenceCell(), DMPlexCreateBoxMesh()
3739 @*/
DMPlexCreateReferenceCellByType(MPI_Comm comm,DMPolytopeType ct,DM * refdm)3740 PetscErrorCode DMPlexCreateReferenceCellByType(MPI_Comm comm, DMPolytopeType ct, DM *refdm)
3741 {
3742 DM rdm;
3743 Vec coords;
3744 PetscErrorCode ierr;
3745
3746 PetscFunctionBegin;
3747 ierr = DMCreate(comm, &rdm);CHKERRQ(ierr);
3748 ierr = DMSetType(rdm, DMPLEX);CHKERRQ(ierr);
3749 switch (ct) {
3750 case DM_POLYTOPE_POINT:
3751 {
3752 PetscInt numPoints[1] = {1};
3753 PetscInt coneSize[1] = {0};
3754 PetscInt cones[1] = {0};
3755 PetscInt coneOrientations[1] = {0};
3756 PetscScalar vertexCoords[1] = {0.0};
3757
3758 ierr = DMSetDimension(rdm, 0);CHKERRQ(ierr);
3759 ierr = DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3760 }
3761 break;
3762 case DM_POLYTOPE_SEGMENT:
3763 {
3764 PetscInt numPoints[2] = {2, 1};
3765 PetscInt coneSize[3] = {2, 0, 0};
3766 PetscInt cones[2] = {1, 2};
3767 PetscInt coneOrientations[2] = {0, 0};
3768 PetscScalar vertexCoords[2] = {-1.0, 1.0};
3769
3770 ierr = DMSetDimension(rdm, 1);CHKERRQ(ierr);
3771 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3772 }
3773 break;
3774 case DM_POLYTOPE_TRIANGLE:
3775 {
3776 PetscInt numPoints[2] = {3, 1};
3777 PetscInt coneSize[4] = {3, 0, 0, 0};
3778 PetscInt cones[3] = {1, 2, 3};
3779 PetscInt coneOrientations[3] = {0, 0, 0};
3780 PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
3781
3782 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr);
3783 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3784 }
3785 break;
3786 case DM_POLYTOPE_QUADRILATERAL:
3787 {
3788 PetscInt numPoints[2] = {4, 1};
3789 PetscInt coneSize[5] = {4, 0, 0, 0, 0};
3790 PetscInt cones[4] = {1, 2, 3, 4};
3791 PetscInt coneOrientations[4] = {0, 0, 0, 0};
3792 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
3793
3794 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr);
3795 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3796 }
3797 break;
3798 case DM_POLYTOPE_SEG_PRISM_TENSOR:
3799 {
3800 PetscInt numPoints[2] = {4, 1};
3801 PetscInt coneSize[5] = {4, 0, 0, 0, 0};
3802 PetscInt cones[4] = {1, 2, 3, 4};
3803 PetscInt coneOrientations[4] = {0, 0, 0, 0};
3804 PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, 1.0};
3805
3806 ierr = DMSetDimension(rdm, 2);CHKERRQ(ierr);
3807 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3808 }
3809 break;
3810 case DM_POLYTOPE_TETRAHEDRON:
3811 {
3812 PetscInt numPoints[2] = {4, 1};
3813 PetscInt coneSize[5] = {4, 0, 0, 0, 0};
3814 PetscInt cones[4] = {1, 3, 2, 4};
3815 PetscInt coneOrientations[4] = {0, 0, 0, 0};
3816 PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0};
3817
3818 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr);
3819 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3820 }
3821 break;
3822 case DM_POLYTOPE_HEXAHEDRON:
3823 {
3824 PetscInt numPoints[2] = {8, 1};
3825 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
3826 PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8};
3827 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3828 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
3829 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
3830
3831 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr);
3832 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3833 }
3834 break;
3835 case DM_POLYTOPE_TRI_PRISM:
3836 {
3837 PetscInt numPoints[2] = {6, 1};
3838 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0};
3839 PetscInt cones[6] = {1, 3, 2, 4, 5, 6};
3840 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
3841 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0,
3842 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
3843
3844 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr);
3845 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3846 }
3847 break;
3848 case DM_POLYTOPE_TRI_PRISM_TENSOR:
3849 {
3850 PetscInt numPoints[2] = {6, 1};
3851 PetscInt coneSize[7] = {6, 0, 0, 0, 0, 0, 0};
3852 PetscInt cones[6] = {1, 2, 3, 4, 5, 6};
3853 PetscInt coneOrientations[6] = {0, 0, 0, 0, 0, 0};
3854 PetscScalar vertexCoords[18] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0,
3855 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, -1.0, 1.0, 1.0};
3856
3857 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr);
3858 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3859 }
3860 break;
3861 case DM_POLYTOPE_QUAD_PRISM_TENSOR:
3862 {
3863 PetscInt numPoints[2] = {8, 1};
3864 PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
3865 PetscInt cones[8] = {1, 2, 3, 4, 5, 6, 7, 8};
3866 PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
3867 PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
3868 -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
3869
3870 ierr = DMSetDimension(rdm, 3);CHKERRQ(ierr);
3871 ierr = DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);CHKERRQ(ierr);
3872 }
3873 break;
3874 default: SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for cell type %s", DMPolytopeTypes[ct]);
3875 }
3876 {
3877 PetscInt Nv, v;
3878
3879 /* Must create the celltype label here so that we do not automatically try to compute the types */
3880 ierr = DMCreateLabel(rdm, "celltype");CHKERRQ(ierr);
3881 ierr = DMPlexSetCellType(rdm, 0, ct);CHKERRQ(ierr);
3882 ierr = DMPlexGetChart(rdm, NULL, &Nv);CHKERRQ(ierr);
3883 for (v = 1; v < Nv; ++v) {ierr = DMPlexSetCellType(rdm, v, DM_POLYTOPE_POINT);CHKERRQ(ierr);}
3884 }
3885 ierr = DMPlexInterpolate(rdm, refdm);CHKERRQ(ierr);
3886 if (rdm->coordinateDM) {
3887 DM ncdm;
3888 PetscSection cs;
3889 PetscInt pEnd = -1;
3890
3891 ierr = DMGetLocalSection(rdm->coordinateDM, &cs);CHKERRQ(ierr);
3892 if (cs) {ierr = PetscSectionGetChart(cs, NULL, &pEnd);CHKERRQ(ierr);}
3893 if (pEnd >= 0) {
3894 ierr = DMClone(*refdm, &ncdm);CHKERRQ(ierr);
3895 ierr = DMCopyDisc(rdm->coordinateDM, ncdm);CHKERRQ(ierr);
3896 ierr = DMSetLocalSection(ncdm, cs);CHKERRQ(ierr);
3897 ierr = DMSetCoordinateDM(*refdm, ncdm);CHKERRQ(ierr);
3898 ierr = DMDestroy(&ncdm);CHKERRQ(ierr);
3899 }
3900 }
3901 ierr = DMGetCoordinatesLocal(rdm, &coords);CHKERRQ(ierr);
3902 if (coords) {
3903 ierr = DMSetCoordinatesLocal(*refdm, coords);CHKERRQ(ierr);
3904 } else {
3905 ierr = DMGetCoordinates(rdm, &coords);CHKERRQ(ierr);
3906 if (coords) {ierr = DMSetCoordinates(*refdm, coords);CHKERRQ(ierr);}
3907 }
3908 ierr = DMDestroy(&rdm);CHKERRQ(ierr);
3909 PetscFunctionReturn(0);
3910 }
3911
3912 /*@
3913 DMPlexCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
3914
3915 Collective
3916
3917 Input Parameters:
3918 + comm - The communicator
3919 . dim - The spatial dimension
3920 - simplex - Flag for simplex, otherwise use a tensor-product cell
3921
3922 Output Parameter:
3923 . refdm - The reference cell
3924
3925 Level: intermediate
3926
3927 .seealso: DMPlexCreateReferenceCellByType(), DMPlexCreateBoxMesh()
3928 @*/
DMPlexCreateReferenceCell(MPI_Comm comm,PetscInt dim,PetscBool simplex,DM * refdm)3929 PetscErrorCode DMPlexCreateReferenceCell(MPI_Comm comm, PetscInt dim, PetscBool simplex, DM *refdm)
3930 {
3931 PetscErrorCode ierr;
3932
3933 PetscFunctionBeginUser;
3934 switch (dim) {
3935 case 0: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_POINT, refdm);CHKERRQ(ierr);break;
3936 case 1: ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_SEGMENT, refdm);CHKERRQ(ierr);break;
3937 case 2:
3938 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TRIANGLE, refdm);CHKERRQ(ierr);}
3939 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_QUADRILATERAL, refdm);CHKERRQ(ierr);}
3940 break;
3941 case 3:
3942 if (simplex) {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_TETRAHEDRON, refdm);CHKERRQ(ierr);}
3943 else {ierr = DMPlexCreateReferenceCellByType(comm, DM_POLYTOPE_HEXAHEDRON, refdm);CHKERRQ(ierr);}
3944 break;
3945 default:
3946 SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %D", dim);
3947 }
3948 PetscFunctionReturn(0);
3949 }
3950