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 CHKERRCGNS(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 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 51 #if defined(PETSC_HAVE_CGNS) 52 if (rank == 0) { 53 CHKERRCGNS(cg_open(filename, CG_MODE_READ, &cgid)); 54 PetscCheckFalse(cgid <= 0,PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename); 55 } 56 CHKERRQ(DMPlexCreateCGNS(comm, cgid, interpolate, dm)); 57 if (rank == 0) CHKERRCGNS(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 PetscErrorCode ierr; 95 /* Read from file */ 96 char basename[CGIO_MAX_NAME_LENGTH+1]; 97 char buffer[CGIO_MAX_NAME_LENGTH+1]; 98 int dim = 0, physDim = 0, coordDim =0, numVertices = 0, numCells = 0; 99 int nzones = 0; 100 #endif 101 102 PetscFunctionBegin; 103 #if defined(PETSC_HAVE_CGNS) 104 CHKERRMPI(MPI_Comm_rank(comm, &rank)); 105 CHKERRMPI(MPI_Comm_size(comm, &num_proc)); 106 CHKERRQ(DMCreate(comm, dm)); 107 CHKERRQ(DMSetType(*dm, DMPLEX)); 108 109 /* Open CGNS II file and read basic information on rank 0, then broadcast to all processors */ 110 if (rank == 0) { 111 int nbases, z; 112 113 CHKERRCGNS(cg_nbases(cgid, &nbases)); 114 PetscCheckFalse(nbases > 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single base, not %d",nbases); 115 CHKERRCGNS(cg_base_read(cgid, 1, basename, &dim, &physDim)); 116 CHKERRCGNS(cg_nzones(cgid, 1, &nzones)); 117 CHKERRQ(PetscCalloc2(nzones+1, &cellStart, nzones+1, &vertStart)); 118 for (z = 1; z <= nzones; ++z) { 119 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 120 121 CHKERRCGNS(cg_zone_read(cgid, 1, z, buffer, sizes)); 122 numVertices += sizes[0]; 123 numCells += sizes[1]; 124 cellStart[z] += sizes[1] + cellStart[z-1]; 125 vertStart[z] += sizes[0] + vertStart[z-1]; 126 } 127 for (z = 1; z <= nzones; ++z) { 128 vertStart[z] += numCells; 129 } 130 coordDim = dim; 131 } 132 CHKERRMPI(MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH+1, MPI_CHAR, 0, comm)); 133 CHKERRMPI(MPI_Bcast(&dim, 1, MPI_INT, 0, comm)); 134 CHKERRMPI(MPI_Bcast(&coordDim, 1, MPI_INT, 0, comm)); 135 CHKERRMPI(MPI_Bcast(&nzones, 1, MPI_INT, 0, comm)); 136 137 CHKERRQ(PetscObjectSetName((PetscObject) *dm, basename)); 138 CHKERRQ(DMSetDimension(*dm, dim)); 139 CHKERRQ(DMCreateLabel(*dm, "celltype")); 140 CHKERRQ(DMPlexSetChart(*dm, 0, numCells+numVertices)); 141 142 /* Read zone information */ 143 if (rank == 0) { 144 int z, c, c_loc; 145 146 /* Read the cell set connectivity table and build mesh topology 147 CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */ 148 /* First set sizes */ 149 for (z = 1, c = 0; z <= nzones; ++z) { 150 CGNS_ENUMT(ZoneType_t) zonetype; 151 int nsections; 152 CGNS_ENUMT(ElementType_t) cellType; 153 cgsize_t start, end; 154 int nbndry, parentFlag; 155 PetscInt numCorners; 156 DMPolytopeType ctype; 157 158 CHKERRCGNS(cg_zone_type(cgid, 1, z, &zonetype)); 159 PetscCheckFalse(zonetype == CGNS_ENUMV(Structured),PETSC_COMM_SELF,PETSC_ERR_LIB,"Can only handle Unstructured zones for CGNS"); 160 CHKERRCGNS(cg_nsections(cgid, 1, z, &nsections)); 161 PetscCheckFalse(nsections > 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single section, not %d",nsections); 162 CHKERRCGNS(cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag)); 163 /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */ 164 if (cellType == CGNS_ENUMV(MIXED)) { 165 cgsize_t elementDataSize, *elements; 166 PetscInt off; 167 168 CHKERRCGNS(cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize)); 169 CHKERRQ(PetscMalloc1(elementDataSize, &elements)); 170 CHKERRCGNS(cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL)); 171 for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) { 172 switch (elements[off]) { 173 case CGNS_ENUMV(BAR_2): numCorners = 2; ctype = DM_POLYTOPE_SEGMENT; break; 174 case CGNS_ENUMV(TRI_3): numCorners = 3; ctype = DM_POLYTOPE_TRIANGLE; break; 175 case CGNS_ENUMV(QUAD_4): numCorners = 4; ctype = DM_POLYTOPE_QUADRILATERAL; break; 176 case CGNS_ENUMV(TETRA_4): numCorners = 4; ctype = DM_POLYTOPE_TETRAHEDRON; break; 177 case CGNS_ENUMV(PENTA_6): numCorners = 6; ctype = DM_POLYTOPE_TRI_PRISM; break; 178 case CGNS_ENUMV(HEXA_8): numCorners = 8; ctype = DM_POLYTOPE_HEXAHEDRON; break; 179 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[off]); 180 } 181 CHKERRQ(DMPlexSetConeSize(*dm, c, numCorners)); 182 CHKERRQ(DMPlexSetCellType(*dm, c, ctype)); 183 off += numCorners+1; 184 } 185 CHKERRQ(PetscFree(elements)); 186 } else { 187 switch (cellType) { 188 case CGNS_ENUMV(BAR_2): numCorners = 2; ctype = DM_POLYTOPE_SEGMENT; break; 189 case CGNS_ENUMV(TRI_3): numCorners = 3; ctype = DM_POLYTOPE_TRIANGLE; break; 190 case CGNS_ENUMV(QUAD_4): numCorners = 4; ctype = DM_POLYTOPE_QUADRILATERAL; break; 191 case CGNS_ENUMV(TETRA_4): numCorners = 4; ctype = DM_POLYTOPE_TETRAHEDRON; break; 192 case CGNS_ENUMV(PENTA_6): numCorners = 6; ctype = DM_POLYTOPE_TRI_PRISM; break; 193 case CGNS_ENUMV(HEXA_8): numCorners = 8; ctype = DM_POLYTOPE_HEXAHEDRON; break; 194 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType); 195 } 196 for (c_loc = start; c_loc <= end; ++c_loc, ++c) { 197 CHKERRQ(DMPlexSetConeSize(*dm, c, numCorners)); 198 CHKERRQ(DMPlexSetCellType(*dm, c, ctype)); 199 } 200 } 201 } 202 for (v = numCells; v < numCells+numVertices; ++v) { 203 CHKERRQ(DMPlexSetCellType(*dm, v, DM_POLYTOPE_POINT)); 204 } 205 } 206 207 CHKERRQ(DMSetUp(*dm)); 208 209 CHKERRQ(DMCreateLabel(*dm, "zone")); 210 if (rank == 0) { 211 int z, c, c_loc, v_loc; 212 213 CHKERRQ(DMGetLabel(*dm, "zone", &label)); 214 for (z = 1, c = 0; z <= nzones; ++z) { 215 CGNS_ENUMT(ElementType_t) cellType; 216 cgsize_t elementDataSize, *elements, start, end; 217 int nbndry, parentFlag; 218 PetscInt *cone, numc, numCorners, maxCorners = 27; 219 220 CHKERRCGNS(cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag)); 221 numc = end - start; 222 /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */ 223 CHKERRCGNS(cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize)); 224 CHKERRQ(PetscMalloc2(elementDataSize,&elements,maxCorners,&cone)); 225 CHKERRCGNS(cg_poly_elements_read(cgid, 1, z, 1, elements, NULL, NULL)); 226 if (cellType == CGNS_ENUMV(MIXED)) { 227 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 228 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 229 switch (elements[v]) { 230 case CGNS_ENUMV(BAR_2): numCorners = 2; break; 231 case CGNS_ENUMV(TRI_3): numCorners = 3; break; 232 case CGNS_ENUMV(QUAD_4): numCorners = 4; break; 233 case CGNS_ENUMV(TETRA_4): numCorners = 4; break; 234 case CGNS_ENUMV(PENTA_6): numCorners = 6; break; 235 case CGNS_ENUMV(HEXA_8): numCorners = 8; break; 236 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[v]); 237 } 238 ++v; 239 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) { 240 cone[v_loc] = elements[v]+numCells-1; 241 } 242 CHKERRQ(DMPlexReorderCell(*dm, c, cone)); 243 CHKERRQ(DMPlexSetCone(*dm, c, cone)); 244 CHKERRQ(DMLabelSetValue(label, c, z)); 245 } 246 } else { 247 switch (cellType) { 248 case CGNS_ENUMV(BAR_2): numCorners = 2; break; 249 case CGNS_ENUMV(TRI_3): numCorners = 3; break; 250 case CGNS_ENUMV(QUAD_4): numCorners = 4; break; 251 case CGNS_ENUMV(TETRA_4): numCorners = 4; break; 252 case CGNS_ENUMV(PENTA_6): numCorners = 6; break; 253 case CGNS_ENUMV(HEXA_8): numCorners = 8; break; 254 default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType); 255 } 256 /* CGNS uses Fortran-based indexing, DMPlex uses C-style and numbers cell first then vertices. */ 257 for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) { 258 for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) { 259 cone[v_loc] = elements[v]+numCells-1; 260 } 261 CHKERRQ(DMPlexReorderCell(*dm, c, cone)); 262 CHKERRQ(DMPlexSetCone(*dm, c, cone)); 263 CHKERRQ(DMLabelSetValue(label, c, z)); 264 } 265 } 266 CHKERRQ(PetscFree2(elements,cone)); 267 } 268 } 269 270 CHKERRQ(DMPlexSymmetrize(*dm)); 271 CHKERRQ(DMPlexStratify(*dm)); 272 if (interpolate) { 273 DM idm; 274 275 CHKERRQ(DMPlexInterpolate(*dm, &idm)); 276 CHKERRQ(DMDestroy(dm)); 277 *dm = idm; 278 } 279 280 /* Read coordinates */ 281 CHKERRQ(DMSetCoordinateDim(*dm, coordDim)); 282 CHKERRQ(DMGetCoordinateDM(*dm, &cdm)); 283 CHKERRQ(DMGetLocalSection(cdm, &coordSection)); 284 CHKERRQ(PetscSectionSetNumFields(coordSection, 1)); 285 CHKERRQ(PetscSectionSetFieldComponents(coordSection, 0, coordDim)); 286 CHKERRQ(PetscSectionSetChart(coordSection, numCells, numCells + numVertices)); 287 for (v = numCells; v < numCells+numVertices; ++v) { 288 CHKERRQ(PetscSectionSetDof(coordSection, v, dim)); 289 CHKERRQ(PetscSectionSetFieldDof(coordSection, v, 0, coordDim)); 290 } 291 CHKERRQ(PetscSectionSetUp(coordSection)); 292 293 CHKERRQ(DMCreateLocalVector(cdm, &coordinates)); 294 CHKERRQ(VecGetArray(coordinates, &coords)); 295 if (rank == 0) { 296 PetscInt off = 0; 297 float *x[3]; 298 int z, d; 299 300 CHKERRQ(PetscMalloc3(numVertices,&x[0],numVertices,&x[1],numVertices,&x[2])); 301 for (z = 1; z <= nzones; ++z) { 302 CGNS_ENUMT(DataType_t) datatype; 303 cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */ 304 cgsize_t range_min[3] = {1, 1, 1}; 305 cgsize_t range_max[3] = {1, 1, 1}; 306 int ngrids, ncoords; 307 308 CHKERRCGNS(cg_zone_read(cgid, 1, z, buffer, sizes)); 309 range_max[0] = sizes[0]; 310 CHKERRCGNS(cg_ngrids(cgid, 1, z, &ngrids)); 311 PetscCheckFalse(ngrids > 1,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single grid, not %d",ngrids); 312 CHKERRCGNS(cg_ncoords(cgid, 1, z, &ncoords)); 313 PetscCheckFalse(ncoords != coordDim,PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a coordinate array for each dimension, not %d",ncoords); 314 for (d = 0; d < coordDim; ++d) { 315 CHKERRCGNS(cg_coord_info(cgid, 1, z, 1+d, &datatype, buffer)); 316 CHKERRCGNS(cg_coord_read(cgid, 1, z, buffer, CGNS_ENUMV(RealSingle), range_min, range_max, x[d])); 317 } 318 if (coordDim >= 1) { 319 for (v = 0; v < sizes[0]; ++v) coords[(v+off)*coordDim+0] = x[0][v]; 320 } 321 if (coordDim >= 2) { 322 for (v = 0; v < sizes[0]; ++v) coords[(v+off)*coordDim+1] = x[1][v]; 323 } 324 if (coordDim >= 3) { 325 for (v = 0; v < sizes[0]; ++v) coords[(v+off)*coordDim+2] = x[2][v]; 326 } 327 off += sizes[0]; 328 } 329 CHKERRQ(PetscFree3(x[0],x[1],x[2])); 330 } 331 CHKERRQ(VecRestoreArray(coordinates, &coords)); 332 333 CHKERRQ(PetscObjectSetName((PetscObject) coordinates, "coordinates")); 334 CHKERRQ(VecSetBlockSize(coordinates, coordDim)); 335 CHKERRQ(DMSetCoordinatesLocal(*dm, coordinates)); 336 CHKERRQ(VecDestroy(&coordinates)); 337 338 /* Read boundary conditions */ 339 CHKERRQ(DMGetNumLabels(*dm, &labelIdRange[0])); 340 if (rank == 0) { 341 CGNS_ENUMT(BCType_t) bctype; 342 CGNS_ENUMT(DataType_t) datatype; 343 CGNS_ENUMT(PointSetType_t) pointtype; 344 cgsize_t *points; 345 PetscReal *normals; 346 int normal[3]; 347 char *bcname = buffer; 348 cgsize_t npoints, nnormals; 349 int z, nbc, bc, c, ndatasets; 350 351 for (z = 1; z <= nzones; ++z) { 352 CHKERRCGNS(cg_nbocos(cgid, 1, z, &nbc)); 353 for (bc = 1; bc <= nbc; ++bc) { 354 CHKERRCGNS(cg_boco_info(cgid, 1, z, bc, bcname, &bctype, &pointtype, &npoints, normal, &nnormals, &datatype, &ndatasets)); 355 CHKERRQ(DMCreateLabel(*dm, bcname)); 356 CHKERRQ(DMGetLabel(*dm, bcname, &label)); 357 CHKERRQ(PetscMalloc2(npoints, &points, nnormals, &normals)); 358 CHKERRCGNS(cg_boco_read(cgid, 1, z, bc, points, (void *) normals)); 359 if (pointtype == CGNS_ENUMV(ElementRange)) { 360 /* Range of cells: assuming half-open interval since the documentation sucks */ 361 for (c = points[0]; c < points[1]; ++c) { 362 CHKERRQ(DMLabelSetValue(label, c - cellStart[z-1], 1)); 363 } 364 } else if (pointtype == CGNS_ENUMV(ElementList)) { 365 /* List of cells */ 366 for (c = 0; c < npoints; ++c) { 367 CHKERRQ(DMLabelSetValue(label, points[c] - cellStart[z-1], 1)); 368 } 369 } else if (pointtype == CGNS_ENUMV(PointRange)) { 370 CGNS_ENUMT(GridLocation_t) gridloc; 371 372 /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */ 373 CHKERRCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 374 CHKERRCGNS(cg_gridlocation_read(&gridloc)); 375 /* Range of points: assuming half-open interval since the documentation sucks */ 376 for (c = points[0]; c < points[1]; ++c) { 377 if (gridloc == CGNS_ENUMV(Vertex)) CHKERRQ(DMLabelSetValue(label, c - vertStart[z-1], 1)); 378 else CHKERRQ(DMLabelSetValue(label, c - cellStart[z-1], 1)); 379 } 380 } else if (pointtype == CGNS_ENUMV(PointList)) { 381 CGNS_ENUMT(GridLocation_t) gridloc; 382 383 /* List of points: Oh please, someone get the CGNS developers away from a computer. This is unconscionable. */ 384 CHKERRCGNS(cg_goto(cgid, 1, "Zone_t", z, "BC_t", bc, "end")); 385 CHKERRCGNS(cg_gridlocation_read(&gridloc)); 386 for (c = 0; c < npoints; ++c) { 387 if (gridloc == CGNS_ENUMV(Vertex)) CHKERRQ(DMLabelSetValue(label, points[c] - vertStart[z-1], 1)); 388 else CHKERRQ(DMLabelSetValue(label, points[c] - cellStart[z-1], 1)); 389 } 390 } else SETERRQ(comm, PETSC_ERR_SUP, "Unsupported point set type %d", (int) pointtype); 391 CHKERRQ(PetscFree2(points, normals)); 392 } 393 } 394 CHKERRQ(PetscFree2(cellStart, vertStart)); 395 } 396 CHKERRQ(DMGetNumLabels(*dm, &labelIdRange[1])); 397 CHKERRMPI(MPI_Bcast(labelIdRange, 2, MPIU_INT, 0, comm)); 398 399 /* Create BC labels at all processes */ 400 for (labelId = labelIdRange[0]; labelId < labelIdRange[1]; ++labelId) { 401 char *labelName = buffer; 402 size_t len = sizeof(buffer); 403 const char *locName; 404 405 if (rank == 0) { 406 CHKERRQ(DMGetLabelByNum(*dm, labelId, &label)); 407 CHKERRQ(PetscObjectGetName((PetscObject)label, &locName)); 408 CHKERRQ(PetscStrncpy(labelName, locName, len)); 409 } 410 CHKERRMPI(MPI_Bcast(labelName, (PetscMPIInt)len, MPIU_INT, 0, comm)); 411 CHKERRMPI(DMCreateLabel(*dm, labelName)); 412 } 413 PetscFunctionReturn(0); 414 #else 415 SETERRQ(comm, PETSC_ERR_SUP, "This method requires CGNS support. Reconfigure using --with-cgns-dir"); 416 #endif 417 } 418