Actual source code: plexcgns.c
petsc-3.5.4 2015-05-23
1: #define PETSCDM_DLL
2: #include <petsc-private/dmpleximpl.h> /*I "petscdmplex.h" I*/
4: #if defined(PETSC_HAVE_CGNS)
5: #include <cgnslib.h>
6: #include <cgns_io.h>
7: #endif
11: /*@C
12: DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file.
14: Collective on comm
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
21: Output Parameter:
22: . dm - The DM object representing the mesh
24: Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html
26: Level: beginner
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;
35: #if defined(PETSC_HAVE_CGNS)
36: int cgid = -1;
37: #endif
41: MPI_Comm_rank(comm, &rank);
42: #if defined(PETSC_HAVE_CGNS)
43: if (!rank) {
44: cg_open(filename, CG_MODE_READ, &cgid);
45: if (cgid <= 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_LIB, "cg_open(\"%s\",...) did not return a valid file ID", filename);
46: }
47: DMPlexCreateCGNS(comm, cgid, interpolate, dm);
48: if (!rank) {cg_close(cgid);}
49: #else
50: SETERRQ(comm, PETSC_ERR_SUP, "Loading meshes requires CGNS support. Reconfigure using --with-cgns-dir");
51: #endif
52: return(0);
53: }
57: /*@
58: DMPlexCreateCGNS - Create a DMPlex mesh from a CGNS file ID.
60: Collective on comm
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
67: Output Parameter:
68: . dm - The DM object representing the mesh
70: Note: http://www.grc.nasa.gov/WWW/cgns/CGNS_docs_current/index.html
72: Level: beginner
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 coordSize, v;
86: /* Read from file */
87: char basename[CGIO_MAX_NAME_LENGTH+1];
88: char buffer[CGIO_MAX_NAME_LENGTH+1];
89: int dim = 0, physDim = 0, numVertices = 0, numCells = 0;
90: int nzones = 0;
91: #endif
94: #if defined(PETSC_HAVE_CGNS)
95: MPI_Comm_rank(comm, &rank);
96: MPI_Comm_size(comm, &num_proc);
97: DMCreate(comm, dm);
98: DMSetType(*dm, DMPLEX);
99: /* Open CGNS II file and read basic informations on rank 0, then broadcast to all processors */
100: if (!rank) {
101: int nbases, z;
103: cg_nbases(cgid, &nbases);
104: if (nbases > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single base, not %d\n",nbases);
105: cg_base_read(cgid, 1, basename, &dim, &physDim);
106: cg_nzones(cgid, 1, &nzones);
107: for (z = 1; z <= nzones; ++z) {
108: cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */
110: cg_zone_read(cgid, 1, z, buffer, sizes);
111: numVertices += sizes[0];
112: numCells += sizes[1];
113: }
114: }
115: MPI_Bcast(basename, CGIO_MAX_NAME_LENGTH+1, MPI_CHAR, 0, comm);
116: MPI_Bcast(&dim, 1, MPI_INT, 0, comm);
117: MPI_Bcast(&nzones, 1, MPI_INT, 0, comm);
118: PetscObjectSetName((PetscObject) *dm, basename);
119: DMPlexSetDimension(*dm, dim);
120: DMPlexSetChart(*dm, 0, numCells+numVertices);
122: /* Read zone information */
123: if (!rank) {
124: int z, c, c_loc, v, v_loc;
126: /* Read the cell set connectivity table and build mesh topology
127: CGNS standard requires that cells in a zone be numbered sequentially and be pairwise disjoint. */
128: /* First set sizes */
129: for (z = 1, c = 0; z <= nzones; ++z) {
130: ZoneType_t zonetype;
131: int nsections;
132: ElementType_t cellType;
133: cgsize_t start, end;
134: int nbndry, parentFlag;
135: PetscInt numCorners;
137: cg_zone_type(cgid, 1, z, &zonetype);
138: if (zonetype == Structured) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Can only handle Unstructured zones for CGNS");
139: cg_nsections(cgid, 1, z, &nsections);
140: if (nsections > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single section, not %d\n",nsections);
141: cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
142: /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
143: if (cellType == MIXED) {
144: cgsize_t elementDataSize, *elements;
145: PetscInt off;
147: cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
148: PetscMalloc1(elementDataSize, &elements);
149: cg_elements_read(cgid, 1, z, 1, elements, NULL);
150: for (c_loc = start, off = 0; c_loc <= end; ++c_loc, ++c) {
151: switch (elements[off]) {
152: case TRI_3: numCorners = 3;break;
153: case QUAD_4: numCorners = 4;break;
154: case TETRA_4: numCorners = 4;break;
155: case HEXA_8: numCorners = 8;break;
156: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[off]);
157: }
158: DMPlexSetConeSize(*dm, c, numCorners);
159: off += numCorners+1;
160: }
161: PetscFree(elements);
162: } else {
163: switch (cellType) {
164: case TRI_3: numCorners = 3;break;
165: case QUAD_4: numCorners = 4;break;
166: case TETRA_4: numCorners = 4;break;
167: case HEXA_8: numCorners = 8;break;
168: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType);
169: }
170: for (c_loc = start; c_loc <= end; ++c_loc, ++c) {
171: DMPlexSetConeSize(*dm, c, numCorners);
172: }
173: }
174: }
175: DMSetUp(*dm);
176: for (z = 1, c = 0; z <= nzones; ++z) {
177: ElementType_t cellType;
178: cgsize_t *elements, elementDataSize, start, end;
179: int nbndry, parentFlag;
180: PetscInt *cone, numc, numCorners, maxCorners = 27;
182: cg_section_read(cgid, 1, z, 1, buffer, &cellType, &start, &end, &nbndry, &parentFlag);
183: numc = end - start;
184: /* This alone is reason enough to bludgeon every single CGNDS developer, this must be what they describe as the "idiocy of crowds" */
185: cg_ElementDataSize(cgid, 1, z, 1, &elementDataSize);
186: PetscMalloc2(elementDataSize,&elements,maxCorners,&cone);
187: cg_elements_read(cgid, 1, z, 1, elements, NULL);
188: if (cellType == MIXED) {
189: /* CGNS uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
190: for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
191: switch (elements[v]) {
192: case TRI_3: numCorners = 3;break;
193: case QUAD_4: numCorners = 4;break;
194: case TETRA_4: numCorners = 4;break;
195: case HEXA_8: numCorners = 8;break;
196: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) elements[v]);
197: }
198: ++v;
199: for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) {
200: cone[v_loc] = elements[v]+numCells-1;
201: }
202: /* Tetrahedra are inverted */
203: if (elements[v] == TETRA_4) {
204: PetscInt tmp = cone[0];
205: cone[0] = cone[1];
206: cone[1] = tmp;
207: }
208: /* Hexahedra are inverted */
209: if (elements[v] == HEXA_8) {
210: PetscInt tmp = cone[1];
211: cone[1] = cone[3];
212: cone[3] = tmp;
213: }
214: DMPlexSetCone(*dm, c, cone);
215: DMPlexSetLabelValue(*dm, "zone", c, z);
216: }
217: } else {
218: switch (cellType) {
219: case TRI_3: numCorners = 3;break;
220: case QUAD_4: numCorners = 4;break;
221: case TETRA_4: numCorners = 4;break;
222: case HEXA_8: numCorners = 8;break;
223: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid cell type %d", (int) cellType);
224: }
226: /* CGNS uses Fortran-based indexing, sieve uses C-style and numbers cell first then vertices. */
227: for (c_loc = 0, v = 0; c_loc <= numc; ++c_loc, ++c) {
228: for (v_loc = 0; v_loc < numCorners; ++v_loc, ++v) {
229: cone[v_loc] = elements[v]+numCells-1;
230: }
231: /* Tetrahedra are inverted */
232: if (cellType == TETRA_4) {
233: PetscInt tmp = cone[0];
234: cone[0] = cone[1];
235: cone[1] = tmp;
236: }
237: /* Hexahedra are inverted */
238: if (cellType == HEXA_8) {
239: PetscInt tmp = cone[1];
240: cone[1] = cone[3];
241: cone[3] = tmp;
242: }
243: DMPlexSetCone(*dm, c, cone);
244: DMPlexSetLabelValue(*dm, "zone", c, z);
245: }
246: }
247: PetscFree2(elements,cone);
248: }
249: }
250: DMPlexSymmetrize(*dm);
251: DMPlexStratify(*dm);
252: if (interpolate) {
253: DM idm = NULL;
255: DMPlexInterpolate(*dm, &idm);
256: /* Maintain zone label */
257: {
258: DMLabel label;
260: DMPlexRemoveLabel(*dm, "zone", &label);
261: if (label) {DMPlexAddLabel(idm, label);}
262: }
263: DMDestroy(dm);
264: *dm = idm;
265: }
267: /* Read coordinates */
268: DMGetCoordinateSection(*dm, &coordSection);
269: PetscSectionSetNumFields(coordSection, 1);
270: PetscSectionSetFieldComponents(coordSection, 0, dim);
271: PetscSectionSetChart(coordSection, numCells, numCells + numVertices);
272: for (v = numCells; v < numCells+numVertices; ++v) {
273: PetscSectionSetDof(coordSection, v, dim);
274: PetscSectionSetFieldDof(coordSection, v, 0, dim);
275: }
276: PetscSectionSetUp(coordSection);
277: PetscSectionGetStorageSize(coordSection, &coordSize);
278: VecCreate(comm, &coordinates);
279: PetscObjectSetName((PetscObject) coordinates, "coordinates");
280: VecSetSizes(coordinates, coordSize, PETSC_DETERMINE);
281: VecSetType(coordinates,VECSTANDARD);
282: VecGetArray(coordinates, &coords);
283: if (!rank) {
284: PetscInt off = 0;
285: float *x[3];
286: int z, d;
288: PetscMalloc3(numVertices,&x[0],numVertices,&x[1],numVertices,&x[2]);
289: for (z = 1; z <= nzones; ++z) {
290: DataType_t datatype;
291: cgsize_t sizes[3]; /* Number of vertices, number of cells, number of boundary vertices */
292: cgsize_t range_min[3] = {1, 1, 1};
293: cgsize_t range_max[3] = {1, 1, 1};
294: int ngrids, ncoords;
297: cg_zone_read(cgid, 1, z, buffer, sizes);
298: range_max[0] = sizes[0];
299: cg_ngrids(cgid, 1, z, &ngrids);
300: if (ngrids > 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a single grid, not %d\n",ngrids);
301: cg_ncoords(cgid, 1, z, &ncoords);
302: if (ncoords != dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"CGNS file must have a coordinate array for each dimension, not %d\n",ncoords);
303: for (d = 0; d < dim; ++d) {
304: cg_coord_info(cgid, 1, z, 1+d, &datatype, buffer);
305: cg_coord_read(cgid, 1, z, buffer, RealSingle, range_min, range_max, x[d]);
306: }
307: if (dim > 0) {
308: for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+0] = x[0][v];
309: }
310: if (dim > 1) {
311: for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+1] = x[1][v];
312: }
313: if (dim > 2) {
314: for (v = 0; v < sizes[0]; ++v) coords[(v+off)*dim+2] = x[2][v];
315: }
316: off += sizes[0];
317: }
318: PetscFree3(x[0],x[1],x[2]);
319: }
320: VecRestoreArray(coordinates, &coords);
321: DMSetCoordinatesLocal(*dm, coordinates);
322: VecDestroy(&coordinates);
323: #else
324: SETERRQ(comm, PETSC_ERR_SUP, "This method requires CGNS support. Reconfigure using --with-cgns-dir");
325: #endif
326: return(0);
327: }