Actual source code: plexcgns2.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>
  3: #include <petsc/private/viewercgnsimpl.h>

  5: #include <pcgnslib.h>
  6: #include <cgns_io.h>

  8: #if !defined(CGNS_ENUMT)
  9:   #define CGNS_ENUMT(a) a
 10: #endif
 11: #if !defined(CGNS_ENUMV)
 12:   #define CGNS_ENUMV(a) a
 13: #endif

 15: PetscErrorCode DMPlexCreateCGNSFromFile_Internal(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 16: {
 17:   PetscMPIInt rank;
 18:   int         cgid = -1;

 21:   MPI_Comm_rank(comm, &rank);
 22:   if (rank == 0) {
 23:     cg_open(filename, CG_MODE_READ, &cgid);
 25:   }
 26:   DMPlexCreateCGNS(comm, cgid, interpolate, dm);
 27:   if (rank == 0) cg_close(cgid);
 28:   return 0;
 29: }

 31: PetscErrorCode DMPlexCreateCGNS_Internal(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm)
 32: {
 33:   PetscMPIInt  num_proc, rank;
 34:   DM           cdm;
 35:   DMLabel      label;
 36:   PetscSection coordSection;
 37:   Vec          coordinates;
 38:   PetscScalar *coords;
 39:   PetscInt    *cellStart, *vertStart, v;
 40:   PetscInt     labelIdRange[2], labelId;
 41:   /* Read from file */
 42:   char basename[CGIO_MAX_NAME_LENGTH + 1];
 43:   char buffer[CGIO_MAX_NAME_LENGTH + 1];
 44:   int  dim = 0, physDim = 0, coordDim = 0, numVertices = 0, numCells = 0;
 45:   int  nzones = 0;

 47:   MPI_Comm_rank(comm, &rank);
 48:   MPI_Comm_size(comm, &num_proc);
 49:   DMCreate(comm, dm);
 50:   DMSetType(*dm, DMPLEX);

 52:   /* Open CGNS II file and read basic information on rank 0, then broadcast to all processors */
 53:   if (rank == 0) {
 54:     int nbases, z;

 56:     cg_nbases(cgid, &nbases);
 58:     cg_base_read(cgid, 1, basename, &dim, &physDim);
 59:     cg_nzones(cgid, 1, &nzones);
 60:     PetscCalloc2(nzones + 1, &cellStart, nzones + 1, &vertStart);
 61:     for (z = 1; z <= nzones; ++z) {
 62:       cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */

 64:       cg_zone_read(cgid, 1, z, buffer, sizes);
 65:       numVertices += sizes[0];
 66:       numCells += sizes[1];
 67:       cellStart[z] += sizes[1] + cellStart[z - 1];
 68:       vertStart[z] += sizes[0] + vertStart[z - 1];
 69:     }
 70:     for (z = 1; z <= nzones; ++z) vertStart[z] += numCells;
 71:     coordDim = dim;
 72:   }
 73:   MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH + 1, MPI_CHAR, 0, comm);
 74:   MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
 75:   MPI_Bcast(&coordDim, 1, MPI_INT, 0, comm);
 76:   MPI_Bcast(&nzones, 1, MPI_INT, 0, comm);

 78:   PetscObjectSetName((PetscObject)*dm, basename);
 79:   DMSetDimension(*dm, dim);
 80:   DMCreateLabel(*dm, "celltype");
 81:   DMPlexSetChart(*dm, 0, numCells + numVertices);

 83:   /* Read zone information */
 84:   if (rank == 0) {
 85:     int z, c, c_loc;

 87:     /* Read the cell set connectivity table and build mesh topology
 88:        CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */
 89:     /* First set sizes */
 90:     for (z = 1, c = 0; z <= nzones; ++z) {
 91:       CGNS_ENUMT(ZoneType_t) zonetype;
 92:       int nsections;
 93:       CGNS_ENUMT(ElementType_t) cellType;
 94:       cgsize_t       start, end;
 95:       int            nbndry, parentFlag;
 96:       PetscInt       numCorners;
 97:       DMPolytopeType ctype;

 99:       cg_zone_type(cgid, 1, z, &zonetype);
101:       cg_nsections(cgid, 1, z, &nsections);
103:       cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
104:       /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
105:       if (cellType == CGNS_ENUMV(MIXED)) {
106:         cgsize_t elementDataSize, *elements;
107:         PetscInt off;

109:         cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
110:         PetscMalloc1(elementDataSize, &elements);
111:         cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL);
112:         for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) {
113:           switch (elements[off]) {
114:           case CGNS_ENUMV(BAR_2):
115:             numCorners = 2;
116:             ctype      = DM_POLYTOPE_SEGMENT;
117:             break;
118:           case CGNS_ENUMV(TRI_3):
119:             numCorners = 3;
120:             ctype      = DM_POLYTOPE_TRIANGLE;
121:             break;
122:           case CGNS_ENUMV(QUAD_4):
123:             numCorners = 4;
124:             ctype      = DM_POLYTOPE_QUADRILATERAL;
125:             break;
126:           case CGNS_ENUMV(TETRA_4):
127:             numCorners = 4;
128:             ctype      = DM_POLYTOPE_TETRAHEDRON;
129:             break;
130:           case CGNS_ENUMV(PENTA_6):
131:             numCorners = 6;
132:             ctype      = DM_POLYTOPE_TRI_PRISM;
133:             break;
134:           case CGNS_ENUMV(HEXA_8):
135:             numCorners = 8;
136:             ctype      = DM_POLYTOPE_HEXAHEDRON;
137:             break;
138:           default:
139:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int)elements[off]);
140:           }
141:           DMPlexSetConeSize(*dm, c, numCorners);
142:           DMPlexSetCellType(*dm, c, ctype);
143:           off += numCorners + 1;
144:         }
145:         PetscFree(elements);
146:       } else {
147:         switch (cellType) {
148:         case CGNS_ENUMV(BAR_2):
149:           numCorners = 2;
150:           ctype      = DM_POLYTOPE_SEGMENT;
151:           break;
152:         case CGNS_ENUMV(TRI_3):
153:           numCorners = 3;
154:           ctype      = DM_POLYTOPE_TRIANGLE;
155:           break;
156:         case CGNS_ENUMV(QUAD_4):
157:           numCorners = 4;
158:           ctype      = DM_POLYTOPE_QUADRILATERAL;
159:           break;
160:         case CGNS_ENUMV(TETRA_4):
161:           numCorners = 4;
162:           ctype      = DM_POLYTOPE_TETRAHEDRON;
163:           break;
164:         case CGNS_ENUMV(PENTA_6):
165:           numCorners = 6;
166:           ctype      = DM_POLYTOPE_TRI_PRISM;
167:           break;
168:         case CGNS_ENUMV(HEXA_8):
169:           numCorners = 8;
170:           ctype      = DM_POLYTOPE_HEXAHEDRON;
171:           break;
172:         default:
173:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int)cellType);
174:         }
175:         for (c_loc = start; c_loc <= end; ++c_loc, ++c) {
176:           DMPlexSetConeSize(*dm, c, numCorners);
177:           DMPlexSetCellType(*dm, c, ctype);
178:         }
179:       }
180:     }
181:     for (v = numCells; v < numCells + numVertices; ++v) DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT);
182:   }

184:   DMSetUp(*dm);

186:   DMCreateLabel(*dm, "zone");
187:   if (rank == 0) {
188:     int z, c, c_loc, v_loc;

190:     DMGetLabel(*dm, "zone", &label);
191:     for (z = 1, c = 0; z <= nzones; ++z) {
192:       CGNS_ENUMT(ElementType_t) cellType;
193:       cgsize_t  elementDataSize, *elements, start, end;
194:       int       nbndry, parentFlag;
195:       PetscInt *cone, numc, numCorners, maxCorners = 27;

197:       cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
198:       numc = end - start;
199:       /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
200:       cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
201:       PetscMalloc2(elementDataSize, &elements, maxCorners, &cone);
202:       cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL);
203:       if (cellType == CGNS_ENUMV(MIXED)) {
204:         /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
205:         for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
206:           switch (elements[v]) {
207:           case CGNS_ENUMV(BAR_2):
208:             numCorners = 2;
209:             break;
210:           case CGNS_ENUMV(TRI_3):
211:             numCorners = 3;
212:             break;
213:           case CGNS_ENUMV(QUAD_4):
214:             numCorners = 4;
215:             break;
216:           case CGNS_ENUMV(TETRA_4):
217:             numCorners = 4;
218:             break;
219:           case CGNS_ENUMV(PENTA_6):
220:             numCorners = 6;
221:             break;
222:           case CGNS_ENUMV(HEXA_8):
223:             numCorners = 8;
224:             break;
225:           default:
226:             SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int)elements[v]);
227:           }
228:           ++v;
229:           for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1;
230:           DMPlexReorderCell(*dm, c, cone);
231:           DMPlexSetCone(*dm, c, cone);
232:           DMLabelSetValue(label, c, z);
233:         }
234:       } else {
235:         switch (cellType) {
236:         case CGNS_ENUMV(BAR_2):
237:           numCorners = 2;
238:           break;
239:         case CGNS_ENUMV(TRI_3):
240:           numCorners = 3;
241:           break;
242:         case CGNS_ENUMV(QUAD_4):
243:           numCorners = 4;
244:           break;
245:         case CGNS_ENUMV(TETRA_4):
246:           numCorners = 4;
247:           break;
248:         case CGNS_ENUMV(PENTA_6):
249:           numCorners = 6;
250:           break;
251:         case CGNS_ENUMV(HEXA_8):
252:           numCorners = 8;
253:           break;
254:         default:
255:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int)cellType);
256:         }
257:         /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */
258:         for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
259:           for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) cone[v_loc] = elements[v] + numCells - 1;
260:           DMPlexReorderCell(*dm, c, cone);
261:           DMPlexSetCone(*dm, c, cone);
262:           DMLabelSetValue(label, c, z);
263:         }
264:       }
265:       PetscFree2(elements, cone);
266:     }
267:   }

269:   DMPlexSymmetrize(*dm);
270:   DMPlexStratify(*dm);
271:   if (interpolate) {
272:     DM idm;

274:     DMPlexInterpolate(*dm, &idm);
275:     DMDestroy(dm);
276:     *dm = idm;
277:   }

279:   /* Read coordinates */
280:   DMSetCoordinateDim(*dm, coordDim);
281:   DMGetCoordinateDM(*dm, &cdm);
282:   DMGetLocalSection(cdm, &coordSection);
283:   PetscSectionSetNumFields(coordSection, 1);
284:   PetscSectionSetFieldComponents(coordSection, 0, coordDim);
285:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
286:   for (v = numCells; v < numCells + numVertices; ++v) {
287:     PetscSectionSetDof(coordSection, v, dim);
288:     PetscSectionSetFieldDof(coordSection, v, 0, coordDim);
289:   }
290:   PetscSectionSetUp(coordSection);

292:   DMCreateLocalVector(cdm, &coordinates);
293:   VecGetArray(coordinates, &coords);
294:   if (rank == 0) {
295:     PetscInt off = 0;
296:     float   *x[3];
297:     int      z, d;

299:     PetscMalloc3(numVertices, &x[0], numVertices, &x[1], numVertices, &x[2]);
300:     for (z = 1; z <= nzones; ++z) {
301:       CGNS_ENUMT(DataType_t) datatype;
302:       cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */
303:       cgsize_t range_min[3] = {1, 1, 1};
304:       cgsize_t range_max[3] = {1, 1, 1};
305:       int      ngrids, ncoords;

307:       cg_zone_read(cgid, 1, z, buffer, sizes);
308:       range_max[0] = sizes[0];
309:       cg_ngrids(cgid, 1, z, &ngrids);
311:       cg_ncoords(cgid, 1, z, &ncoords);
313:       for (d = 0; d < coordDim; ++d) {
314:         cg_coord_info(cgid, 1, z, 1 + d, &datatype, buffer);
315:         cg_coord_read(cgid, 1, z, buffer, CGNS_ENUMV(RealSingle), range_min, range_max, x[d]);
316:       }
317:       if (coordDim >= 1) {
318:         for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 0] = x[0][v];
319:       }
320:       if (coordDim >= 2) {
321:         for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 1] = x[1][v];
322:       }
323:       if (coordDim >= 3) {
324:         for (v = 0; v < sizes[0]; ++v) coords[(v + off) * coordDim + 2] = x[2][v];
325:       }
326:       off += sizes[0];
327:     }
328:     PetscFree3(x[0], x[1], x[2]);
329:   }
330:   VecRestoreArray(coordinates, &coords);

332:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
333:   VecSetBlockSize(coordinates, coordDim);
334:   DMSetCoordinatesLocal(*dm, coordinates);
335:   VecDestroy(&coordinates);

337:   /* Read boundary conditions */
338:   DMGetNumLabels(*dm, &labelIdRange[0]);
339:   if (rank == 0) {
340:     CGNS_ENUMT(BCType_t) bctype;
341:     CGNS_ENUMT(DataType_t) datatype;
342:     CGNS_ENUMT(PointSetType_t) pointtype;
343:     cgsize_t  *points;
344:     PetscReal *normals;
345:     int        normal[3];
346:     char      *bcname = buffer;
347:     cgsize_t   npoints, nnormals;
348:     int        z, nbc, bc, c, ndatasets;

350:     for (z = 1; z <= nzones; ++z) {
351:       cg_nbocos(cgid, 1, z, &nbc);
352:       for (bc = 1; bc <= nbc; ++bc) {
353:         cg_boco_info(cgid, 1, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets);
354:         DMCreateLabel(*dm, bcname);
355:         DMGetLabel(*dm, bcname, &label);
356:         PetscMalloc2(npoints, &points, nnormals, &normals);
357:         cg_boco_read(cgid, 1, z, bc, points, (void *)normals);
358:         if (pointtype == CGNS_ENUMV(ElementRange)) {
359:           /* Range of cells: assuming half-open interval since the documentation sucks */
360:           for (c = points[0]; c < points[1]; ++c) DMLabelSetValue(label, c - cellStart[z - 1], 1);
361:         } else if (pointtype == CGNS_ENUMV(ElementList)) {
362:           /* List of cells */
363:           for (c = 0; c < npoints; ++c) DMLabelSetValue(label, points[c] - cellStart[z - 1], 1);
364:         } else if (pointtype == CGNS_ENUMV(PointRange)) {
365:           CGNS_ENUMT(GridLocation_t) gridloc;

367:           /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */
368:           cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end");
369:           cg_gridlocation_read(&gridloc);
370:           /* Range of points: assuming half-open interval since the documentation sucks */
371:           for (c = points[0]; c < points[1]; ++c) {
372:             if (gridloc == CGNS_ENUMV(Vertex)) DMLabelSetValue(label, c - vertStart[z - 1], 1);
373:             else DMLabelSetValue(label, c - cellStart[z - 1], 1);
374:           }
375:         } else if (pointtype == CGNS_ENUMV(PointList)) {
376:           CGNS_ENUMT(GridLocation_t) gridloc;

378:           /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */
379:           cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end");
380:           cg_gridlocation_read(&gridloc);
381:           for (c = 0; c < npoints; ++c) {
382:             if (gridloc == CGNS_ENUMV(Vertex)) DMLabelSetValue(label, points[c] - vertStart[z - 1], 1);
383:             else DMLabelSetValue(label, points[c] - cellStart[z - 1], 1);
384:           }
385:         } else SETERRQ(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int)pointtype);
386:         PetscFree2(points, normals);
387:       }
388:     }
389:     PetscFree2(cellStart, vertStart);
390:   }
391:   DMGetNumLabels(*dm, &labelIdRange[1]);
392:   MPI_Bcast(labelIdRange, 2, MPIU_INT, 0, comm);

394:   /* Create BC labels at all processes */
395:   for (labelId = labelIdRange[0]; labelId < labelIdRange[1]; ++labelId) {
396:     char       *labelName = buffer;
397:     size_t      len       = sizeof(buffer);
398:     const char *locName;

400:     if (rank == 0) {
401:       DMGetLabelByNum(*dm, labelId, &label);
402:       PetscObjectGetName((PetscObject)label, &locName);
403:       PetscStrncpy(labelName, locName, len);
404:     }
405:     MPI_Bcast(labelName, (PetscMPIInt)len, MPIU_INT, 0, comm);
406:     DMCreateLabel(*dm, labelName);
407:   }
408:   return 0;
409: }

411: // Permute plex closure ordering to CGNS
412: static PetscErrorCode DMPlexCGNSGetPermutation_Internal(DMPolytopeType cell_type, PetscInt closure_size, CGNS_ENUMT(ElementType_t) * element_type, const int **perm)
413: {
414:   // https://cgns.github.io/CGNS_docs_current/sids/conv.html#unst_example
415:   static const int bar_2[2]   = {0, 1};
416:   static const int bar_3[3]   = {1, 2, 0};
417:   static const int bar_4[4]   = {2, 3, 0, 1};
418:   static const int bar_5[5]   = {3, 4, 0, 1, 2};
419:   static const int tri_3[3]   = {0, 1, 2};
420:   static const int tri_6[6]   = {3, 4, 5, 0, 1, 2};
421:   static const int tri_10[10] = {7, 8, 9, 1, 2, 3, 4, 5, 6, 0};
422:   static const int quad_4[4]  = {0, 1, 2, 3};
423:   static const int quad_9[9]  = {
424:     5, 6, 7, 8, // vertices
425:     1, 2, 3, 4, // edges
426:     0,          // center
427:   };
428:   static const int quad_16[] = {
429:     12, 13, 14, 15,               // vertices
430:     4,  5,  6,  7,  8, 9, 10, 11, // edges
431:     0,  1,  3,  2,                // centers
432:   };
433:   static const int quad_25[] = {
434:     21, 22, 23, 24,                                 // vertices
435:     9,  10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, // edges
436:     0,  1,  2,  5,  8,  7,  6,  3,  4,              // centers
437:   };
438:   static const int tetra_4[4]   = {0, 2, 1, 3};
439:   static const int tetra_10[10] = {6, 8, 7, 9, 2, 1, 0, 3, 5, 4};
440:   static const int tetra_20[20] = {
441:     16, 18, 17, 19,         // vertices
442:     9,  8,  7,  6,  5,  4,  // bottom edges
443:     10, 11, 14, 15, 13, 12, // side edges
444:     0,  2,  3,  1,          // faces
445:   };
446:   static const int hexa_8[8]   = {0, 3, 2, 1, 4, 5, 6, 7};
447:   static const int hexa_27[27] = {
448:     19, 22, 21, 20, 23, 24, 25, 26, // vertices
449:     10, 9,  8,  7,                  // bottom edges
450:     16, 15, 18, 17,                 // mid edges
451:     11, 12, 13, 14,                 // top edges
452:     1,  3,  5,  4,  6,  2,          // faces
453:     0,                              // center
454:   };
455:   static const int hexa_64[64] = {
456:     // debug with $PETSC_ARCH/tests/dm/impls/plex/tests/ex49 -dm_plex_simplex 0 -dm_plex_dim 3 -dm_plex_box_faces 1,1,1 -dm_coord_petscspace_degree 3
457:     56, 59, 58, 57, 60, 61, 62, 63, // vertices
458:     39, 38, 37, 36, 35, 34, 33, 32, // bottom edges
459:     51, 50, 48, 49, 52, 53, 55, 54, // mid edges; Paraview needs edge 21-22 swapped with 23-24
460:     40, 41, 42, 43, 44, 45, 46, 47, // top edges
461:     8,  10, 11, 9,                  // z-minus face
462:     16, 17, 19, 18,                 // y-minus face
463:     24, 25, 27, 26,                 // x-plus face
464:     20, 21, 23, 22,                 // y-plus face
465:     30, 28, 29, 31,                 // x-minus face
466:     12, 13, 15, 14,                 // z-plus face
467:     0,  1,  3,  2,  4,  5,  7,  6,  // center
468:   };

470:   *element_type = CGNS_ENUMV(ElementTypeNull);
471:   *perm         = NULL;
472:   switch (cell_type) {
473:   case DM_POLYTOPE_SEGMENT:
474:     switch (closure_size) {
475:     case 2:
476:       *element_type = CGNS_ENUMV(BAR_2);
477:       *perm         = bar_2;
478:     case 3:
479:       *element_type = CGNS_ENUMV(BAR_3);
480:       *perm         = bar_3;
481:     case 4:
482:       *element_type = CGNS_ENUMV(BAR_4);
483:       *perm         = bar_4;
484:       break;
485:     case 5:
486:       *element_type = CGNS_ENUMV(BAR_5);
487:       *perm         = bar_5;
488:       break;
489:     default:
490:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
491:     }
492:     break;
493:   case DM_POLYTOPE_TRIANGLE:
494:     switch (closure_size) {
495:     case 3:
496:       *element_type = CGNS_ENUMV(TRI_3);
497:       *perm         = tri_3;
498:       break;
499:     case 6:
500:       *element_type = CGNS_ENUMV(TRI_6);
501:       *perm         = tri_6;
502:       break;
503:     case 10:
504:       *element_type = CGNS_ENUMV(TRI_10);
505:       *perm         = tri_10;
506:       break;
507:     default:
508:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
509:     }
510:     break;
511:   case DM_POLYTOPE_QUADRILATERAL:
512:     switch (closure_size) {
513:     case 4:
514:       *element_type = CGNS_ENUMV(QUAD_4);
515:       *perm         = quad_4;
516:       break;
517:     case 9:
518:       *element_type = CGNS_ENUMV(QUAD_9);
519:       *perm         = quad_9;
520:       break;
521:     case 16:
522:       *element_type = CGNS_ENUMV(QUAD_16);
523:       *perm         = quad_16;
524:       break;
525:     case 25:
526:       *element_type = CGNS_ENUMV(QUAD_25);
527:       *perm         = quad_25;
528:       break;
529:     default:
530:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
531:     }
532:     break;
533:   case DM_POLYTOPE_TETRAHEDRON:
534:     switch (closure_size) {
535:     case 4:
536:       *element_type = CGNS_ENUMV(TETRA_4);
537:       *perm         = tetra_4;
538:       break;
539:     case 10:
540:       *element_type = CGNS_ENUMV(TETRA_10);
541:       *perm         = tetra_10;
542:       break;
543:     case 20:
544:       *element_type = CGNS_ENUMV(TETRA_20);
545:       *perm         = tetra_20;
546:       break;
547:     default:
548:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
549:     }
550:     break;
551:   case DM_POLYTOPE_HEXAHEDRON:
552:     switch (closure_size) {
553:     case 8:
554:       *element_type = CGNS_ENUMV(HEXA_8);
555:       *perm         = hexa_8;
556:       break;
557:     case 27:
558:       *element_type = CGNS_ENUMV(HEXA_27);
559:       *perm         = hexa_27;
560:       break;
561:     case 64:
562:       *element_type = CGNS_ENUMV(HEXA_64);
563:       *perm         = hexa_64;
564:       break;
565:     default:
566:       SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
567:     }
568:     break;
569:   default:
570:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cell type %s with closure size %" PetscInt_FMT, DMPolytopeTypes[cell_type], closure_size);
571:   }
572:   return 0;
573: }

575: // node_l2g must be freed
576: static PetscErrorCode DMPlexCreateNodeNumbering(DM dm, PetscInt *num_local_nodes, PetscInt *num_global_nodes, PetscInt *nStart, PetscInt *nEnd, const PetscInt **node_l2g)
577: {
578:   PetscSection    local_section;
579:   PetscSF         point_sf;
580:   PetscInt        pStart, pEnd, spStart, spEnd, *points, nleaves, ncomp, *nodes;
581:   PetscMPIInt     comm_size;
582:   const PetscInt *ilocal, field = 0;

584:   *num_local_nodes  = -1;
585:   *num_global_nodes = -1;
586:   *nStart           = -1;
587:   *nEnd             = -1;
588:   *node_l2g         = NULL;

590:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &comm_size);
591:   DMGetLocalSection(dm, &local_section);
592:   DMPlexGetChart(dm, &pStart, &pEnd);
593:   PetscSectionGetChart(local_section, &spStart, &spEnd);
594:   DMGetPointSF(dm, &point_sf);
595:   if (comm_size == 1) nleaves = 0;
596:   else PetscSFGetGraph(point_sf, NULL, &nleaves, &ilocal, NULL);
597:   PetscSectionGetFieldComponents(local_section, field, &ncomp);

599:   PetscInt local_node = 0, owned_node = 0, owned_start = 0;
600:   for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) {
601:     PetscInt dof;
602:     PetscSectionGetFieldDof(local_section, p, field, &dof);
603:     PetscAssert(dof % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field dof %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, dof, ncomp);
604:     local_node += dof / ncomp;
605:     if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process
606:       leaf++;
607:     } else {
608:       owned_node += dof / ncomp;
609:     }
610:   }
611:   MPI_Exscan(&owned_node, &owned_start, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
612:   PetscMalloc1(pEnd - pStart, &points);
613:   owned_node = 0;
614:   for (PetscInt p = spStart, leaf = 0; p < spEnd; p++) {
615:     if (leaf < nleaves && p == ilocal[leaf]) { // skip points owned by a different process
616:       points[p - pStart] = -1;
617:       leaf++;
618:       continue;
619:     }
620:     PetscInt dof, offset;
621:     PetscSectionGetFieldDof(local_section, p, field, &dof);
622:     PetscSectionGetFieldOffset(local_section, p, field, &offset);
623:     PetscAssert(offset % ncomp == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Field offset %" PetscInt_FMT " must be divisible by components %" PetscInt_FMT, offset, ncomp);
624:     points[p - pStart] = owned_start + owned_node;
625:     owned_node += dof / ncomp;
626:   }
627:   if (comm_size > 1) {
628:     PetscSFBcastBegin(point_sf, MPIU_INT, points, points, MPI_REPLACE);
629:     PetscSFBcastEnd(point_sf, MPIU_INT, points, points, MPI_REPLACE);
630:   }

632:   // Set up global indices for each local node
633:   PetscMalloc1(local_node, &nodes);
634:   for (PetscInt p = spStart; p < spEnd; p++) {
635:     PetscInt dof, offset;
636:     PetscSectionGetFieldDof(local_section, p, field, &dof);
637:     PetscSectionGetFieldOffset(local_section, p, field, &offset);
638:     for (PetscInt n = 0; n < dof / ncomp; n++) nodes[offset / ncomp + n] = points[p - pStart] + n;
639:   }
640:   PetscFree(points);
641:   *num_local_nodes = local_node;
642:   *nStart          = owned_start;
643:   *nEnd            = owned_start + owned_node;
644:   MPI_Allreduce(&owned_node, num_global_nodes, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
645:   *node_l2g = nodes;
646:   return 0;
647: }

649: PetscErrorCode DMView_PlexCGNS(DM dm, PetscViewer viewer)
650: {
651:   PetscViewer_CGNS *cgv = (PetscViewer_CGNS *)viewer->data;
652:   PetscInt          topo_dim, coord_dim, num_global_elems;
653:   PetscInt          cStart, cEnd, num_local_nodes, num_global_nodes, nStart, nEnd;
654:   const PetscInt   *node_l2g;
655:   Vec               coord;
656:   DM                colloc_dm, cdm;
657:   PetscMPIInt       size;
658:   const char       *dm_name;
659:   int               base, zone;
660:   cgsize_t          isize[3];

662:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
663:   DMGetDimension(dm, &topo_dim);
664:   DMGetCoordinateDim(dm, &coord_dim);
665:   PetscObjectGetName((PetscObject)dm, &dm_name);
666:   cg_base_write(cgv->file_num, dm_name, topo_dim, coord_dim, &base);
667:   cg_goto(cgv->file_num, base, NULL);
668:   cg_dataclass_write(CGNS_ENUMV(NormalizedByDimensional));

670:   {
671:     PetscFE        fe, fe_coord;
672:     PetscDualSpace dual_space, dual_space_coord;
673:     PetscInt       field_order, field_order_coord;
674:     PetscBool      is_simplex;
675:     DMGetField(dm, 0, NULL, (PetscObject *)&fe);
676:     if (fe) {
677:       PetscFEGetDualSpace(fe, &dual_space);
678:       PetscDualSpaceGetOrder(dual_space, &field_order);
679:     } else field_order = 1;
680:     DMGetCoordinateDM(dm, &cdm);
681:     DMGetField(cdm, 0, NULL, (PetscObject *)&fe_coord);
682:     if (fe_coord) {
683:       PetscFEGetDualSpace(fe_coord, &dual_space_coord);
684:       PetscDualSpaceGetOrder(dual_space_coord, &field_order_coord);
685:     } else field_order_coord = 1;
686:     if (field_order != field_order_coord) {
687:       PetscInt quadrature_order = field_order;
688:       DMClone(dm, &colloc_dm);
689:       DMPlexIsSimplex(dm, &is_simplex);
690:       PetscFECreateLagrange(PetscObjectComm((PetscObject)dm), topo_dim, coord_dim, is_simplex, field_order, quadrature_order, &fe);
691:       DMProjectCoordinates(colloc_dm, fe);
692:       PetscFEDestroy(&fe);
693:     } else {
694:       PetscObjectReference((PetscObject)dm);
695:       colloc_dm = dm;
696:     }
697:   }
698:   DMGetCoordinateDM(colloc_dm, &cdm);
699:   DMPlexCreateNodeNumbering(cdm, &num_local_nodes, &num_global_nodes, &nStart, &nEnd, &node_l2g);
700:   DMGetCoordinatesLocal(colloc_dm, &coord);
701:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
702:   num_global_elems = cEnd - cStart;
703:   MPI_Allreduce(MPI_IN_PLACE, &num_global_elems, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
704:   isize[0] = num_global_nodes;
705:   isize[1] = num_global_elems;
706:   isize[2] = 0;
707:   cg_zone_write(cgv->file_num, base, "Zone", isize, CGNS_ENUMV(Unstructured), &zone);

709:   {
710:     const PetscScalar *X;
711:     PetscScalar       *x;
712:     int                coord_ids[3];

714:     VecGetArrayRead(coord, &X);
715:     for (PetscInt d = 0; d < coord_dim; d++) {
716:       const double exponents[] = {0, 1, 0, 0, 0};
717:       char         coord_name[64];
718:       PetscSNPrintf(coord_name, sizeof coord_name, "Coordinate%c", 'X' + (int)d);
719:       cgp_coord_write(cgv->file_num, base, zone, CGNS_ENUMV(RealDouble), coord_name, &coord_ids[d]);
720:       cg_goto(cgv->file_num, base, "Zone_t", zone, "GridCoordinates", 0, coord_name, 0, NULL);
721:       cg_exponents_write(CGNS_ENUMV(RealDouble), exponents);
722:     }

724:     DMPolytopeType cell_type;
725:     int            section;
726:     cgsize_t       e_owned, e_global, e_start, *conn = NULL;
727:     const int     *perm;
728:     CGNS_ENUMT(ElementType_t) element_type = CGNS_ENUMV(ElementTypeNull);
729:     {
730:       PetscMalloc1(nEnd - nStart, &x);
731:       for (PetscInt d = 0; d < coord_dim; d++) {
732:         for (PetscInt n = 0; n < num_local_nodes; n++) {
733:           PetscInt gn = node_l2g[n];
734:           if (gn < nStart || nEnd <= gn) continue;
735:           x[gn - nStart] = X[n * coord_dim + d];
736:         }
737:         // CGNS nodes use 1-based indexing
738:         cgsize_t start = nStart + 1, end = nEnd;
739:         cgp_coord_write_data(cgv->file_num, base, zone, coord_ids[d], &start, &end, x);
740:       }
741:       PetscFree(x);
742:       VecRestoreArrayRead(coord, &X);
743:     }

745:     DMPlexGetCellType(dm, cStart, &cell_type);
746:     for (PetscInt i = cStart, c = 0; i < cEnd; i++) {
747:       PetscInt closure_dof, *closure_indices, elem_size;
748:       DMPlexGetClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL);
749:       elem_size = closure_dof / coord_dim;
750:       if (!conn) PetscMalloc1((cEnd - cStart) * elem_size, &conn);
751:       DMPlexCGNSGetPermutation_Internal(cell_type, closure_dof / coord_dim, &element_type, &perm);
752:       for (PetscInt j = 0; j < elem_size; j++) conn[c++] = node_l2g[closure_indices[perm[j] * coord_dim] / coord_dim] + 1;
753:       DMPlexRestoreClosureIndices(cdm, cdm->localSection, cdm->localSection, i, PETSC_FALSE, &closure_dof, &closure_indices, NULL, NULL);
754:     }
755:     e_owned = cEnd - cStart;
756:     MPI_Allreduce(&e_owned, &e_global, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)dm));
758:     e_start = 0;
759:     MPI_Exscan(&e_owned, &e_start, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)dm));
760:     cgp_section_write(cgv->file_num, base, zone, "Elem", element_type, 1, e_global, 0, &section);
761:     cgp_elements_write_data(cgv->file_num, base, zone, section, e_start + 1, e_start + e_owned, conn);
762:     PetscFree(conn);

764:     cgv->base            = base;
765:     cgv->zone            = zone;
766:     cgv->node_l2g        = node_l2g;
767:     cgv->num_local_nodes = num_local_nodes;
768:     cgv->nStart          = nStart;
769:     cgv->nEnd            = nEnd;
770:     if (1) {
771:       PetscMPIInt rank;
772:       int        *efield;
773:       int         sol, field;
774:       DMLabel     label;
775:       MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
776:       PetscMalloc1(e_owned, &efield);
777:       for (PetscInt i = 0; i < e_owned; i++) efield[i] = rank;
778:       cg_sol_write(cgv->file_num, base, zone, "CellInfo", CGNS_ENUMV(CellCenter), &sol);
779:       cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "Rank", &field);
780:       cgsize_t start = e_start + 1, end = e_start + e_owned;
781:       cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield);
782:       DMGetLabel(dm, "Cell Sets", &label);
783:       if (label) {
784:         for (PetscInt c = cStart; c < cEnd; c++) {
785:           PetscInt value;
786:           DMLabelGetValue(label, c, &value);
787:           efield[c - cStart] = value;
788:         }
789:         cgp_field_write(cgv->file_num, base, zone, sol, CGNS_ENUMV(Integer), "CellSet", &field);
790:         cgp_field_write_data(cgv->file_num, base, zone, sol, field, &start, &end, efield);
791:       }
792:       PetscFree(efield);
793:     }
794:   }
795:   DMDestroy(&colloc_dm);
796:   return 0;
797: }

799: static PetscErrorCode PetscCGNSDataType(PetscDataType pd, CGNS_ENUMT(DataType_t) * cd)
800: {
801:   switch (pd) {
802:   case PETSC_FLOAT:
803:     *cd = CGNS_ENUMV(RealSingle);
804:     break;
805:   case PETSC_DOUBLE:
806:     *cd = CGNS_ENUMV(RealDouble);
807:     break;
808:   case PETSC_COMPLEX:
809:     *cd = CGNS_ENUMV(ComplexDouble);
810:     break;
811:   default:
812:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Data type %s", PetscDataTypes[pd]);
813:   }
814:   return 0;
815: }

817: PetscErrorCode VecView_Plex_Local_CGNS(Vec V, PetscViewer viewer)
818: {
819:   PetscViewer_CGNS  *cgv = (PetscViewer_CGNS *)viewer->data;
820:   DM                 dm;
821:   PetscSection       section;
822:   PetscInt           ncomp, time_step;
823:   PetscReal          time, *time_slot;
824:   const PetscInt     field = 0;
825:   const PetscScalar *v;
826:   char               solution_name[PETSC_MAX_PATH_LEN];
827:   int                sol;

829:   VecGetDM(V, &dm);
830:   if (!cgv->node_l2g) DMView(dm, viewer);
831:   if (!cgv->nodal_field) PetscMalloc1(cgv->nEnd - cgv->nStart, &cgv->nodal_field);
832:   if (!cgv->output_times) PetscSegBufferCreate(sizeof(PetscReal), 20, &cgv->output_times);

834:   DMGetOutputSequenceNumber(dm, &time_step, &time);
835:   if (time_step < 0) {
836:     time_step = 0;
837:     time      = 0.;
838:   }
839:   PetscSegBufferGet(cgv->output_times, 1, &time_slot);
840:   *time_slot = time;
841:   PetscSNPrintf(solution_name, sizeof solution_name, "FlowSolution%" PetscInt_FMT, time_step);
842:   cg_sol_write(cgv->file_num, cgv->base, cgv->zone, solution_name, CGNS_ENUMV(Vertex), &sol);
843:   DMGetLocalSection(dm, &section);
844:   PetscSectionGetFieldComponents(section, field, &ncomp);
845:   VecGetArrayRead(V, &v);
846:   for (PetscInt comp = 0; comp < ncomp; comp++) {
847:     int         cgfield;
848:     const char *comp_name;
849:     CGNS_ENUMT(DataType_t) datatype;
850:     PetscSectionGetComponentName(section, field, comp, &comp_name);
851:     PetscCGNSDataType(PETSC_SCALAR, &datatype);
852:     cgp_field_write(cgv->file_num, cgv->base, cgv->zone, sol, datatype, comp_name, &cgfield);
853:     for (PetscInt n = 0; n < cgv->num_local_nodes; n++) {
854:       PetscInt gn = cgv->node_l2g[n];
855:       if (gn < cgv->nStart || cgv->nEnd <= gn) continue;
856:       cgv->nodal_field[gn - cgv->nStart] = v[n * ncomp + comp];
857:     }
858:     // CGNS nodes use 1-based indexing
859:     cgsize_t start = cgv->nStart + 1, end = cgv->nEnd;
860:     cgp_field_write_data(cgv->file_num, cgv->base, cgv->zone, sol, cgfield, &start, &end, cgv->nodal_field);
861:   }
862:   VecRestoreArrayRead(V, &v);
863:   return 0;
864: }