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