Actual source code: plexcgns.c

petsc-3.8.2 2017-11-09
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmpleximpl.h>

  4: #undef I /* Very old CGNS stupidly uses I as a variable, which fails when using complex. Curse you idiot package managers */
  5: #if defined(PETSC_HAVE_CGNS)
  6: #include <cgnslib.h>
  7: #include <cgns_io.h>
  8: #endif
  9: #if !defined(CGNS_ENUMT)
 10: #define CGNS_ENUMT(a) a
 11: #endif
 12: #if !defined(CGNS_ENUMV)
 13: #define CGNS_ENUMV(a) a
 14: #endif

 16: /*@C
 17:   DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file.

 19:   Collective on comm

 21:   Input Parameters:
 22: + comm  - The MPI communicator
 23: . filename - The name of the CGNS file
 24: - interpolate - Create faces and edges in the mesh

 26:   Output Parameter:
 27: . dm  - The DM object representing the mesh

 29:   Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html

 31:   Level: beginner

 33: .keywords: mesh,CGNS
 34: .seealso: DMPlexCreate(), DMPlexCreateCGNS(), DMPlexCreateExodus()
 35: @*/
 36: PetscErrorCode DMPlexCreateCGNSFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 37: {
 38:   PetscMPIInt    rank;
 40: #if defined(PETSC_HAVE_CGNS)
 41:   int cgid = -1;
 42: #endif

 46:   MPI_Comm_rank(comm, &rank);
 47: #if defined(PETSC_HAVE_CGNS)
 48:   if (!rank) {
 49:     cg_open(filename, CG_MODE_READ, &cgid);
 50:     if (cgid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename);
 51:   }
 52:   DMPlexCreateCGNS(comm, cgid, interpolate, dm);
 53:   if (!rank) {cg_close(cgid);}
 54: #else
 55:   SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir");
 56: #endif
 57:   return(0);
 58: }

 60: /*@
 61:   DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file ID.

 63:   Collective on comm

 65:   Input Parameters:
 66: + comm  - The MPI communicator
 67: . cgid - The CG id associated with a file and obtained using cg_open
 68: - interpolate - Create faces and edges in the mesh

 70:   Output Parameter:
 71: . dm  - The DM object representing the mesh

 73:   Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html

 75:   Level: beginner

 77: .keywords: mesh,CGNS
 78: .seealso: DMPlexCreate(), DMPlexCreateExodus()
 79: @*/
 80: PetscErrorCode DMPlexCreateCGNS(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm)
 81: {
 82: #if defined(PETSC_HAVE_CGNS)
 83:   PetscMPIInt    num_proc, rank;
 84:   PetscSection   coordSection;
 85:   Vec            coordinates;
 86:   PetscScalar   *coords;
 87:   PetscInt      *cellStart, *vertStart;
 88:   PetscInt       coordSize, v;
 90:   /* Read from file */
 91:   char basename[CGIO_MAX_NAME_LENGTH+1];
 92:   char buffer[CGIO_MAX_NAME_LENGTH+1];
 93:   int  dim    = 0, physDim = 0, numVertices = 0, numCells = 0;
 94:   int  nzones = 0;
 95: #endif

 98: #if defined(PETSC_HAVE_CGNS)
 99:   MPI_Comm_rank(comm, &rank);
100:   MPI_Comm_size(comm, &num_proc);
101:   DMCreate(comm, dm);
102:   DMSetType(*dm, DMPLEX);
103:   /* Open CGNS II file and read basic informations on rank 0, then broadcast to all processors */
104:   if (!rank) {
105:     int nbases, z;

107:     cg_nbases(cgid, &nbases);
108:     if (nbases > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single base, not %d\n",nbases);
109:     cg_base_read(cgid, 1, basename, &dim, &physDim);
110:     cg_nzones(cgid, 1, &nzones);
111:     PetscCalloc2(nzones+1, &cellStart, nzones+1, &vertStart);
112:     for (z = 1; z <= nzones; ++z) {
113:       cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */

115:       cg_zone_read(cgid, 1, z, buffer, sizes);
116:       numVertices += sizes[0];
117:       numCells    += sizes[1];
118:       cellStart[z] += sizes[1] + cellStart[z-1];
119:       vertStart[z] += sizes[0] + vertStart[z-1];
120:     }
121:     for (z = 1; z <= nzones; ++z) {
122:       vertStart[z] += numCells;
123:     }
124:   }
125:   MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH+1, MPI_CHAR, 0, comm);
126:   MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
127:   MPI_Bcast(&nzones, 1, MPI_INT, 0, comm);
128:   PetscObjectSetName((PetscObject) *dm, basename);
129:   DMSetDimension(*dm, dim);
130:   DMPlexSetChart(*dm, 0, numCells+numVertices);

132:   /* Read zone information */
133:   if (!rank) {
134:     int z, c, c_loc, v, v_loc;

136:     /* Read the cell set connectivity table and build mesh topology
137:        CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */
138:     /* First set sizes */
139:     for (z = 1, c = 0; z <= nzones; ++z) {
140:       CGNS_ENUMT( ZoneType_t )    zonetype;
141:       int                         nsections;
142:       CGNS_ENUMT( ElementType_t ) cellType;
143:       cgsize_t                    start, end;
144:       int                         nbndry, parentFlag;
145:       PetscInt                    numCorners;

147:       cg_zone_type(cgid, 1, z, &zonetype);
148:       if (zonetype == CGNS_ENUMV( Structured )) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can only handle Unstructured zones for CGNS");
149:       cg_nsections(cgid, 1, z, &nsections);
150:       if (nsections > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single section, not %d\n",nsections);
151:       cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
152:       /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
153:       if (cellType == CGNS_ENUMV( MIXED )) {
154:         cgsize_t elementDataSize, *elements;
155:         PetscInt off;

157:         cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
158:         PetscMalloc1(elementDataSize, &elements);
159:         cg_elements_read(cgid, 1, z, 1, elements, NULL);
160:         for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) {
161:           switch (elements[off]) {
162:           case CGNS_ENUMV( TRI_3 ):   numCorners = 3;break;
163:           case CGNS_ENUMV( QUAD_4 ):  numCorners = 4;break;
164:           case CGNS_ENUMV( TETRA_4 ): numCorners = 4;break;
165:           case CGNS_ENUMV( HEXA_8 ):  numCorners = 8;break;
166:           default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[off]);
167:           }
168:           DMPlexSetConeSize(*dm, c, numCorners);
169:           off += numCorners+1;
170:         }
171:         PetscFree(elements);
172:       } else {
173:         switch (cellType) {
174:         case CGNS_ENUMV( TRI_3 ):   numCorners = 3;break;
175:         case CGNS_ENUMV( QUAD_4 ):  numCorners = 4;break;
176:         case CGNS_ENUMV( TETRA_4 ): numCorners = 4;break;
177:         case CGNS_ENUMV( HEXA_8 ):  numCorners = 8;break;
178:         default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType);
179:         }
180:         for (c_loc = start; c_loc <= end; ++c_loc, ++c) {
181:           DMPlexSetConeSize(*dm, c, numCorners);
182:         }
183:       }
184:     }
185:     DMSetUp(*dm);
186:     for (z = 1, c = 0; z <= nzones; ++z) {
187:       CGNS_ENUMT( ElementType_t ) cellType;
188:       cgsize_t                    *elements, elementDataSize, start, end;
189:       int                          nbndry, parentFlag;
190:       PetscInt                    *cone, numc, numCorners, maxCorners = 27;

192:       cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
193:       numc = end - start;
194:       /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
195:       cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
196:       PetscMalloc2(elementDataSize,&elements,maxCorners,&cone);
197:       cg_elements_read(cgid, 1, z, 1, elements, NULL);
198:       if (cellType == CGNS_ENUMV( MIXED )) {
199:         /* CGNS uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
200:         for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
201:           switch (elements[v]) {
202:           case CGNS_ENUMV( TRI_3 ):   numCorners = 3;break;
203:           case CGNS_ENUMV( QUAD_4 ):  numCorners = 4;break;
204:           case CGNS_ENUMV( TETRA_4 ): numCorners = 4;break;
205:           case CGNS_ENUMV( HEXA_8 ):  numCorners = 8;break;
206:           default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[v]);
207:           }
208:           ++v;
209:           for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) {
210:             cone[v_loc] = elements[v]+numCells-1;
211:           }
212:           /* Tetrahedra are inverted */
213:           if (elements[v] == CGNS_ENUMV( TETRA_4 )) {
214:             PetscInt tmp = cone[0];
215:             cone[0] = cone[1];
216:             cone[1] = tmp;
217:           }
218:           /* Hexahedra are inverted */
219:           if (elements[v] == CGNS_ENUMV( HEXA_8 )) {
220:             PetscInt tmp = cone[5];
221:             cone[5] = cone[7];
222:             cone[7] = tmp;
223:           }
224:           DMPlexSetCone(*dm, c, cone);
225:           DMSetLabelValue(*dm, "zone", c, z);
226:         }
227:       } else {
228:         switch (cellType) {
229:         case CGNS_ENUMV( TRI_3 ):   numCorners = 3;break;
230:         case CGNS_ENUMV( QUAD_4 ):  numCorners = 4;break;
231:         case CGNS_ENUMV( TETRA_4 ): numCorners = 4;break;
232:         case CGNS_ENUMV( HEXA_8 ):  numCorners = 8;break;
233:         default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType);
234:         }

236:         /* CGNS uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
237:         for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
238:           for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) {
239:             cone[v_loc] = elements[v]+numCells-1;
240:           }
241:           /* Tetrahedra are inverted */
242:           if (cellType == CGNS_ENUMV( TETRA_4 )) {
243:             PetscInt tmp = cone[0];
244:             cone[0] = cone[1];
245:             cone[1] = tmp;
246:           }
247:           /* Hexahedra are inverted, and they give the top first */
248:           if (cellType == CGNS_ENUMV( HEXA_8 )) {
249:             PetscInt tmp = cone[5];
250:             cone[5] = cone[7];
251:             cone[7] = tmp;
252:           }
253:           DMPlexSetCone(*dm, c, cone);
254:           DMSetLabelValue(*dm, "zone", c, z);
255:         }
256:       }
257:       PetscFree2(elements,cone);
258:     }
259:   }
260:   DMPlexSymmetrize(*dm);
261:   DMPlexStratify(*dm);
262:   if (interpolate) {
263:     DM idm = NULL;

265:     DMPlexInterpolate(*dm, &idm);
266:     /* Maintain zone label */
267:     {
268:       DMLabel label;

270:       DMRemoveLabel(*dm, "zone", &label);
271:       if (label) {DMAddLabel(idm, label);}
272:     }
273:     DMDestroy(dm);
274:     *dm  = idm;
275:   }

277:   /* Read coordinates */
278:   DMGetCoordinateSection(*dm, &coordSection);
279:   PetscSectionSetNumFields(coordSection, 1);
280:   PetscSectionSetFieldComponents(coordSection, 0, dim);
281:   PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
282:   for (v = numCells; v < numCells+numVertices; ++v) {
283:     PetscSectionSetDof(coordSection, v, dim);
284:     PetscSectionSetFieldDof(coordSection, v, 0, dim);
285:   }
286:   PetscSectionSetUp(coordSection);
287:   PetscSectionGetStorageSize(coordSection, &coordSize);
288:   VecCreate(PETSC_COMM_SELF, &coordinates);
289:   PetscObjectSetName((PetscObject) coordinates, "coordinates");
290:   VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
291:   VecSetType(coordinates,VECSTANDARD);
292:   VecGetArray(coordinates, &coords);
293:   if (!rank) {
294:     PetscInt off = 0;
295:     float   *x[3];
296:     int      z, d;

298:     PetscMalloc3(numVertices,&x[0],numVertices,&x[1],numVertices,&x[2]);
299:     for (z = 1; z <= nzones; ++z) {
300:       CGNS_ENUMT( DataType_t ) datatype;
301:       cgsize_t                 sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */
302:       cgsize_t                 range_min[3] = {1, 1, 1};
303:       cgsize_t                 range_max[3] = {1, 1, 1};
304:       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);
310:       if (ngrids > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single grid, not %d\n",ngrids);
311:       cg_ncoords(cgid, 1, z, &ncoords);
312:       if (ncoords != dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a coordinate array for each dimension, not %d\n",ncoords);
313:       for (d = 0; d < dim; ++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 (dim > 0) {
318:         for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+0] = x[0][v];
319:       }
320:       if (dim > 1) {
321:         for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+1] = x[1][v];
322:       }
323:       if (dim > 2) {
324:         for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+2] = x[2][v];
325:       }
326:       off += sizes[0];
327:     }
328:     PetscFree3(x[0],x[1],x[2]);
329:   }
330:   VecRestoreArray(coordinates, &coords);
331:   DMSetCoordinatesLocal(*dm, coordinates);
332:   VecDestroy(&coordinates);
333:   /* Read boundary conditions */
334:   if (!rank) {
335:     DMLabel                       label;
336:     CGNS_ENUMT( BCType_t )        bctype;
337:     CGNS_ENUMT( DataType_t )      datatype;
338:     CGNS_ENUMT( PointSetType_t )  pointtype;
339:     cgsize_t                     *points;
340:     PetscReal                    *normals;
341:     int                           normal[3];
342:     char                         *bcname = buffer;
343:     cgsize_t                      npoints, nnormals;
344:     int                           z, nbc, bc, c, ndatasets;

346:     for (z = 1; z <= nzones; ++z) {
347:       cg_nbocos(cgid, 1, z, &nbc);
348:       for (bc = 1; bc <= nbc; ++bc) {
349:         cg_boco_info(cgid, 1, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets);
350:         DMCreateLabel(*dm, bcname);
351:         DMGetLabel(*dm, bcname, &label);
352:         PetscMalloc2(npoints, &points, nnormals, &normals);
353:         cg_boco_read(cgid, 1, z, bc, points, (void *) normals);
354:         if (pointtype == CGNS_ENUMV( ElementRange )) {
355:           /* Range of cells: assuming half-open interval since the documentation sucks */
356:           for (c = points[0]; c < points[1]; ++c) {
357:             DMLabelSetValue(label, c - cellStart[z-1], 1);
358:           }
359:         } else if (pointtype == CGNS_ENUMV( ElementList )) {
360:           /* List of cells */
361:           for (c = 0; c < npoints; ++c) {
362:             DMLabelSetValue(label, points[c] - cellStart[z-1], 1);
363:           }
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 SETERRQ1(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int) pointtype);
386:         PetscFree2(points, normals);
387:       }
388:     }
389:     PetscFree2(cellStart, vertStart);
390:   }
391: #else
392:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires CGNS support. Reconfigure using --with-cgns-dir");
393: #endif
394:   return(0);
395: }