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