1 #define PETSCDM_DLL 2 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 3 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 15 16 #define PetscCallCGNS(ierr) \ 17 do { \ 18 int _cgns_ier = (ierr); \ 19 PetscCheck(!_cgns_ier,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS error %d %s",_cgns_ier,cg_get_error()); \ 20 } while (0) 21 22 /*@C 23 DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file. 24 25 Collective 26 27 Input Parameters: 28 + comm - The MPI communicator 29 . filename - The name of the CGNS file 30 - interpolate - Create faces and edges in the mesh 31 32 Output Parameter: 33 . dm - The DM object representing the mesh 34 35 Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html 36 37 Level: beginner 38 39 .seealso: `DMPlexCreate()`, `DMPlexCreateCGNS()`, `DMPlexCreateExodus()` 40 @*/ 41 PetscErrorCode DMPlexCreateCGNSFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm) 42 { 43 PetscMPIInt rank; 44 #if defined(PETSC_HAVE_CGNS) 45 int cgid = -1; 46 #endif 47 48 PetscFunctionBegin; 49 PetscValidCharPointer(filename, 2); 50 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 51 #if defined(PETSC_HAVE_CGNS) 52 if (rank == 0) { 53 PetscCallCGNS(cg_open(filename, CG_MODE_READ, &cgid)); 54 PetscCheck(cgid > 0,PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename); 55 } 56 PetscCall(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 57 if (rank == 0) PetscCallCGNS(cg_close(cgid)); 58 PetscFunctionReturn(0); 59 #else 60 SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir"); 61 #endif 62 } 63 64 /*@ 65 DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file ID. 66 67 Collective 68 69 Input Parameters: 70 + comm - The MPI communicator 71 . cgid - The CG id associated with a file and obtained using cg_open 72 - interpolate - Create faces and edges in the mesh 73 74 Output Parameter: 75 . dm - The DM object representing the mesh 76 77 Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html 78 79 Level: beginner 80 81 .seealso: `DMPlexCreate()`, `DMPlexCreateExodus()` 82 @*/ 83 PetscErrorCode DMPlexCreateCGNS(MPI_Comm comm, PetscInt cgid, PetscBool interpolate, DM *dm) 84 { 85 #if defined(PETSC_HAVE_CGNS) 86 PetscMPIInt num_proc, rank; 87 DM cdm; 88 DMLabel label; 89 PetscSection coordSection; 90 Vec coordinates; 91 PetscScalar *coords; 92 PetscInt *cellStart, *vertStart, v; 93 PetscInt labelIdRange[2], labelId; 94 /* Read from file */ 95 char basename[CGIO_MAX_NAME_LENGTH+1]; 96 char buffer[CGIO_MAX_NAME_LENGTH+1]; 97 int dim = 0, physDim = 0, coordDim =0, numVertices = 0, numCells = 0; 98 int nzones = 0; 99 #endif 100 101 PetscFunctionBegin; 102 #if defined(PETSC_HAVE_CGNS) 103 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 104 PetscCallMPI(MPI_Comm_size(comm, &num_proc)); 105 PetscCall(DMCreate(comm, dm)); 106 PetscCall(DMSetType(*dm, DMPLEX)); 107 108 /* Open CGNS II file and read basic information on rank 0, then broadcast to all processors */ 109 if (rank == 0) { 110 int nbases, z; 111 112 PetscCallCGNS(cg_nbases(cgid, &nbases)); 113 PetscCheck(nbases <= 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single base, not %d",nbases); 114 PetscCallCGNS(cg_base_read(cgid, 1, basename, &dim, &physDim)); 115 PetscCallCGNS(cg_nzones(cgid, 1, &nzones)); 116 PetscCall(PetscCalloc2(nzones+1, &cellStart, nzones+1, &vertStart)); 117 for (z = 1; z <= nzones; ++z) { 118 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 119 120 PetscCallCGNS(cg_zone_read(cgid, 1, z, buffer, sizes)); 121 numVertices += sizes[0]; 122 numCells += sizes[1]; 123 cellStart[z] += sizes[1] + cellStart[z-1]; 124 vertStart[z] += sizes[0] + vertStart[z-1]; 125 } 126 for (z = 1; z <= nzones; ++z) { 127 vertStart[z] += numCells; 128 } 129 coordDim = dim; 130 } 131 PetscCallMPI(MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH+1, MPI_CHAR, 0, comm)); 132 PetscCallMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 133 PetscCallMPI(MPI_Bcast(&coordDim, 1, MPI_INT, 0, comm)); 134 PetscCallMPI(MPI_Bcast(&nzones, 1, MPI_INT, 0, comm)); 135 136 PetscCall(PetscObjectSetName((PetscObject) *dm, basename)); 137 PetscCall(DMSetDimension(*dm, dim)); 138 PetscCall(DMCreateLabel(*dm, "celltype")); 139 PetscCall(DMPlexSetChart(*dm, 0, numCells+numVertices)); 140 141 /* Read zone information */ 142 if (rank == 0) { 143 int z, c, c_loc; 144 145 /* Read the cell set connectivity table and build mesh topology 146 CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */ 147 /* First set sizes */ 148 for (z = 1, c = 0; z <= nzones; ++z) { 149 CGNS_ENUMT(ZoneType_t) zonetype; 150 int nsections; 151 CGNS_ENUMT(ElementType_t) cellType; 152 cgsize_t start, end; 153 int nbndry, parentFlag; 154 PetscInt numCorners; 155 DMPolytopeType ctype; 156 157 PetscCallCGNS(cg_zone_type(cgid, 1, z, &zonetype)); 158 PetscCheck(zonetype != CGNS_ENUMV(Structured),PETSC_COMM_SELF,PETSC_ERR_LIB,"Can only handle Unstructured zones for CGNS"); 159 PetscCallCGNS(cg_nsections(cgid, 1, z, &nsections)); 160 PetscCheck(nsections <= 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single section, not %d",nsections); 161 PetscCallCGNS(cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag)); 162 /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */ 163 if (cellType == CGNS_ENUMV(MIXED)) { 164 cgsize_t elementDataSize, *elements; 165 PetscInt off; 166 167 PetscCallCGNS(cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize)); 168 PetscCall(PetscMalloc1(elementDataSize, &elements)); 169 PetscCallCGNS(cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL)); 170 for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) { 171 switch (elements[off]) { 172 case CGNS_ENUMV(BAR_2): numCorners = 2; ctype = DM_POLYTOPE_SEGMENT; break; 173 case CGNS_ENUMV(TRI_3): numCorners = 3; ctype = DM_POLYTOPE_TRIANGLE; break; 174 case CGNS_ENUMV(QUAD_4): numCorners = 4; ctype = DM_POLYTOPE_QUADRILATERAL; break; 175 case CGNS_ENUMV(TETRA_4): numCorners = 4; ctype = DM_POLYTOPE_TETRAHEDRON; break; 176 case CGNS_ENUMV(PENTA_6): numCorners = 6; ctype = DM_POLYTOPE_TRI_PRISM; break; 177 case CGNS_ENUMV(HEXA_8): numCorners = 8; ctype = DM_POLYTOPE_HEXAHEDRON; break; 178 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[off]); 179 } 180 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 181 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 182 off += numCorners+1; 183 } 184 PetscCall(PetscFree(elements)); 185 } else { 186 switch (cellType) { 187 case CGNS_ENUMV(BAR_2): numCorners = 2; ctype = DM_POLYTOPE_SEGMENT; break; 188 case CGNS_ENUMV(TRI_3): numCorners = 3; ctype = DM_POLYTOPE_TRIANGLE; break; 189 case CGNS_ENUMV(QUAD_4): numCorners = 4; ctype = DM_POLYTOPE_QUADRILATERAL; break; 190 case CGNS_ENUMV(TETRA_4): numCorners = 4; ctype = DM_POLYTOPE_TETRAHEDRON; break; 191 case CGNS_ENUMV(PENTA_6): numCorners = 6; ctype = DM_POLYTOPE_TRI_PRISM; break; 192 case CGNS_ENUMV(HEXA_8): numCorners = 8; ctype = DM_POLYTOPE_HEXAHEDRON; break; 193 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType); 194 } 195 for (c_loc = start; c_loc <= end; ++c_loc, ++c) { 196 PetscCall(DMPlexSetConeSize(*dm, c, numCorners)); 197 PetscCall(DMPlexSetCellType(*dm, c, ctype)); 198 } 199 } 200 } 201 for (v = numCells; v < numCells+numVertices; ++v) { 202 PetscCall(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 203 } 204 } 205 206 PetscCall(DMSetUp(*dm)); 207 208 PetscCall(DMCreateLabel(*dm, "zone")); 209 if (rank == 0) { 210 int z, c, c_loc, v_loc; 211 212 PetscCall(DMGetLabel(*dm, "zone", &label)); 213 for (z = 1, c = 0; z <= nzones; ++z) { 214 CGNS_ENUMT(ElementType_t) cellType; 215 cgsize_t elementDataSize, *elements, start, end; 216 int nbndry, parentFlag; 217 PetscInt *cone, numc, numCorners, maxCorners = 27; 218 219 PetscCallCGNS(cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag)); 220 numc = end - start; 221 /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */ 222 PetscCallCGNS(cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize)); 223 PetscCall(PetscMalloc2(elementDataSize,&elements,maxCorners,&cone)); 224 PetscCallCGNS(cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL)); 225 if (cellType == CGNS_ENUMV(MIXED)) { 226 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 227 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 228 switch (elements[v]) { 229 case CGNS_ENUMV(BAR_2): numCorners = 2; break; 230 case CGNS_ENUMV(TRI_3): numCorners = 3; break; 231 case CGNS_ENUMV(QUAD_4): numCorners = 4; break; 232 case CGNS_ENUMV(TETRA_4): numCorners = 4; break; 233 case CGNS_ENUMV(PENTA_6): numCorners = 6; break; 234 case CGNS_ENUMV(HEXA_8): numCorners = 8; break; 235 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[v]); 236 } 237 ++v; 238 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) { 239 cone[v_loc] = elements[v]+numCells-1; 240 } 241 PetscCall(DMPlexReorderCell(*dm, c, cone)); 242 PetscCall(DMPlexSetCone(*dm, c, cone)); 243 PetscCall(DMLabelSetValue(label, c, z)); 244 } 245 } else { 246 switch (cellType) { 247 case CGNS_ENUMV(BAR_2): numCorners = 2; break; 248 case CGNS_ENUMV(TRI_3): numCorners = 3; break; 249 case CGNS_ENUMV(QUAD_4): numCorners = 4; break; 250 case CGNS_ENUMV(TETRA_4): numCorners = 4; break; 251 case CGNS_ENUMV(PENTA_6): numCorners = 6; break; 252 case CGNS_ENUMV(HEXA_8): numCorners = 8; break; 253 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType); 254 } 255 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 256 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 257 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) { 258 cone[v_loc] = elements[v]+numCells-1; 259 } 260 PetscCall(DMPlexReorderCell(*dm, c, cone)); 261 PetscCall(DMPlexSetCone(*dm, c, cone)); 262 PetscCall(DMLabelSetValue(label, c, z)); 263 } 264 } 265 PetscCall(PetscFree2(elements,cone)); 266 } 267 } 268 269 PetscCall(DMPlexSymmetrize(*dm)); 270 PetscCall(DMPlexStratify(*dm)); 271 if (interpolate) { 272 DM idm; 273 274 PetscCall(DMPlexInterpolate(*dm, &idm)); 275 PetscCall(DMDestroy(dm)); 276 *dm = idm; 277 } 278 279 /* Read coordinates */ 280 PetscCall(DMSetCoordinateDim(*dm, coordDim)); 281 PetscCall(DMGetCoordinateDM(*dm, &cdm)); 282 PetscCall(DMGetLocalSection(cdm, &coordSection)); 283 PetscCall(PetscSectionSetNumFields(coordSection, 1)); 284 PetscCall(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 285 PetscCall(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 286 for (v = numCells; v < numCells+numVertices; ++v) { 287 PetscCall(PetscSectionSetDof(coordSection, v, dim)); 288 PetscCall(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 289 } 290 PetscCall(PetscSectionSetUp(coordSection)); 291 292 PetscCall(DMCreateLocalVector(cdm, &coordinates)); 293 PetscCall(VecGetArray(coordinates, &coords)); 294 if (rank == 0) { 295 PetscInt off = 0; 296 float *x[3]; 297 int z, d; 298 299 PetscCall(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; 306 307 PetscCallCGNS(cg_zone_read(cgid, 1, z, buffer, sizes)); 308 range_max[0] = sizes[0]; 309 PetscCallCGNS(cg_ngrids(cgid, 1, z, &ngrids)); 310 PetscCheck(ngrids <= 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single grid, not %d",ngrids); 311 PetscCallCGNS(cg_ncoords(cgid, 1, z, &ncoords)); 312 PetscCheck(ncoords == coordDim,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a coordinate array for each dimension, not %d",ncoords); 313 for (d = 0; d < coordDim; ++d) { 314 PetscCallCGNS(cg_coord_info(cgid, 1, z, 1+d, &datatype, buffer)); 315 PetscCallCGNS(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 PetscCall(PetscFree3(x[0],x[1],x[2])); 329 } 330 PetscCall(VecRestoreArray(coordinates, &coords)); 331 332 PetscCall(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 333 PetscCall(VecSetBlockSize(coordinates, coordDim)); 334 PetscCall(DMSetCoordinatesLocal(*dm, coordinates)); 335 PetscCall(VecDestroy(&coordinates)); 336 337 /* Read boundary conditions */ 338 PetscCall(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; 349 350 for (z = 1; z <= nzones; ++z) { 351 PetscCallCGNS(cg_nbocos(cgid, 1, z, &nbc)); 352 for (bc = 1; bc <= nbc; ++bc) { 353 PetscCallCGNS(cg_boco_info(cgid, 1, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets)); 354 PetscCall(DMCreateLabel(*dm, bcname)); 355 PetscCall(DMGetLabel(*dm, bcname, &label)); 356 PetscCall(PetscMalloc2(npoints, &points, nnormals, &normals)); 357 PetscCallCGNS(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) { 361 PetscCall(DMLabelSetValue(label, c - cellStart[z-1], 1)); 362 } 363 } else if (pointtype == CGNS_ENUMV(ElementList)) { 364 /* List of cells */ 365 for (c = 0; c < npoints; ++c) { 366 PetscCall(DMLabelSetValue(label, points[c] - cellStart[z-1], 1)); 367 } 368 } else if (pointtype == CGNS_ENUMV(PointRange)) { 369 CGNS_ENUMT(GridLocation_t) gridloc; 370 371 /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */ 372 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 373 PetscCallCGNS(cg_gridlocation_read(&gridloc)); 374 /* Range of points: assuming half-open interval since the documentation sucks */ 375 for (c = points[0]; c < points[1]; ++c) { 376 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, c - vertStart[z-1], 1)); 377 else PetscCall(DMLabelSetValue(label, c - cellStart[z-1], 1)); 378 } 379 } else if (pointtype == CGNS_ENUMV(PointList)) { 380 CGNS_ENUMT(GridLocation_t) gridloc; 381 382 /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */ 383 PetscCallCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 384 PetscCallCGNS(cg_gridlocation_read(&gridloc)); 385 for (c = 0; c < npoints; ++c) { 386 if (gridloc == CGNS_ENUMV(Vertex)) PetscCall(DMLabelSetValue(label, points[c] - vertStart[z-1], 1)); 387 else PetscCall(DMLabelSetValue(label, points[c] - cellStart[z-1], 1)); 388 } 389 } else SETERRQ(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int) pointtype); 390 PetscCall(PetscFree2(points, normals)); 391 } 392 } 393 PetscCall(PetscFree2(cellStart, vertStart)); 394 } 395 PetscCall(DMGetNumLabels(*dm, &labelIdRange[1])); 396 PetscCallMPI(MPI_Bcast(labelIdRange, 2, MPIU_INT, 0, comm)); 397 398 /* Create BC labels at all processes */ 399 for (labelId = labelIdRange[0]; labelId < labelIdRange[1]; ++labelId) { 400 char *labelName = buffer; 401 size_t len = sizeof(buffer); 402 const char *locName; 403 404 if (rank == 0) { 405 PetscCall(DMGetLabelByNum(*dm, labelId, &label)); 406 PetscCall(PetscObjectGetName((PetscObject)label, &locName)); 407 PetscCall(PetscStrncpy(labelName, locName, len)); 408 } 409 PetscCallMPI(MPI_Bcast(labelName, (PetscMPIInt)len, MPIU_INT, 0, comm)); 410 PetscCallMPI(DMCreateLabel(*dm, labelName)); 411 } 412 PetscFunctionReturn(0); 413 #else 414 SETERRQ(comm, PETSC_ERR_SUP, "This method requires CGNS support. Reconfigure using --with-cgns-dir"); 415 #endif 416 } 417