Actual source code: plexmed.c

  1: #define PETSCDM_DLL
  2: #include <petsc/private/dmpleximpl.h>

  4: #if defined(PETSC_HAVE_MED)
  5:   #include <med.h>
  6: #endif

  8: /*@C
  9:   DMPlexCreateMedFromFile - Create a `DMPLEX` mesh from a (Salome-)Med file.

 11:   Collective

 13: + comm        - The MPI communicator
 14: . filename    - Name of the .med file
 15: - interpolate - Create faces and edges in the mesh

 17:   Output Parameter:
 18: . dm  - The `DM` object representing the mesh

 20:   Level: beginner

 22:   Reference:
 23: . * -  https://www.salome-platform.org/user-section/about/med, http://docs.salome-platform.org/latest/MED_index.html

 25: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
 26: @*/
 27: PetscErrorCode DMPlexCreateMedFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
 28: {
 29:   PetscMPIInt rank, size;
 30: #if defined(PETSC_HAVE_MED)
 31:   med_idt           fileID;
 32:   PetscInt          i, ngeo, spaceDim, meshDim;
 33:   PetscInt          numVertices = 0, numCells = 0, c, numCorners, numCellsLocal, numVerticesLocal;
 34:   med_int          *medCellList;
 35:   PetscInt         *cellList;
 36:   med_float        *coordinates = NULL;
 37:   PetscReal        *vertcoords  = NULL;
 38:   PetscLayout       vLayout, cLayout;
 39:   const PetscInt   *vrange, *crange;
 40:   PetscSF           sfVertices;
 41:   char             *axisname, *unitname, meshname[MED_NAME_SIZE + 1], geotypename[MED_NAME_SIZE + 1];
 42:   char              meshdescription[MED_COMMENT_SIZE + 1], dtunit[MED_SNAME_SIZE + 1];
 43:   med_sorting_type  sortingtype;
 44:   med_mesh_type     meshtype;
 45:   med_axis_type     axistype;
 46:   med_bool          coordinatechangement, geotransformation, hdfok, medok;
 47:   med_geometry_type geotype[2], cellID, facetID;
 48:   med_filter        vfilter = MED_FILTER_INIT;
 49:   med_filter        cfilter = MED_FILTER_INIT;
 50:   med_err           mederr;
 51:   med_int           major, minor, release;
 52: #endif

 54:   PetscFunctionBegin;
 55:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 56:   PetscCallMPI(MPI_Comm_size(comm, &size));
 57: #if defined(PETSC_HAVE_MED)
 58:   mederr = MEDfileCompatibility(filename, &hdfok, &medok);
 59:   PetscCheck(!mederr, comm, PETSC_ERR_ARG_WRONG, "Cannot determine MED file compatibility: %s", filename);
 60:   PetscCheck(hdfok, comm, PETSC_ERR_ARG_WRONG, "Not a compatible HDF format: %s", filename);
 61:   PetscCheck(medok, comm, PETSC_ERR_ARG_WRONG, "Not a compatible MED format: %s", filename);

 63:   fileID = MEDfileOpen(filename, MED_ACC_RDONLY);
 64:   PetscCheck(fileID >= 0, comm, PETSC_ERR_ARG_WRONG, "Unable to open .med mesh file: %s", filename);
 65:   mederr = MEDfileNumVersionRd(fileID, &major, &minor, &release);
 66:   PetscCheck(MEDnMesh(fileID) >= 1, comm, PETSC_ERR_ARG_WRONG, "No meshes found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
 67:   spaceDim = MEDmeshnAxis(fileID, 1);
 68:   PetscCheck(spaceDim >= 1, comm, PETSC_ERR_ARG_WRONG, "Mesh of unknown space dimension found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
 69:   /* Read general mesh information */
 70:   PetscCall(PetscMalloc1(MED_SNAME_SIZE * spaceDim + 1, &axisname));
 71:   PetscCall(PetscMalloc1(MED_SNAME_SIZE * spaceDim + 1, &unitname));
 72:   {
 73:     med_int medMeshDim, medNstep;
 74:     med_int medSpaceDim = spaceDim;

 76:     PetscCallExternal(MEDmeshInfo, fileID, 1, meshname, &medSpaceDim, &medMeshDim, &meshtype, meshdescription, dtunit, &sortingtype, &medNstep, &axistype, axisname, unitname);
 77:     spaceDim = medSpaceDim;
 78:     meshDim  = medMeshDim;
 79:   }
 80:   PetscCall(PetscFree(axisname));
 81:   PetscCall(PetscFree(unitname));
 82:   /* Partition mesh coordinates */
 83:   numVertices = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE, MED_COORDINATE, MED_NO_CMODE, &coordinatechangement, &geotransformation);
 84:   PetscCall(PetscLayoutCreate(comm, &vLayout));
 85:   PetscCall(PetscLayoutSetSize(vLayout, numVertices));
 86:   PetscCall(PetscLayoutSetBlockSize(vLayout, 1));
 87:   PetscCall(PetscLayoutSetUp(vLayout));
 88:   PetscCall(PetscLayoutGetRanges(vLayout, &vrange));
 89:   numVerticesLocal = vrange[rank + 1] - vrange[rank];
 90:   PetscCallExternal(MEDfilterBlockOfEntityCr, fileID, numVertices, 1, spaceDim, MED_ALL_CONSTITUENT, MED_FULL_INTERLACE, MED_COMPACT_STMODE, MED_NO_PROFILE, vrange[rank] + 1, 1, numVerticesLocal, 1, 1, &vfilter);
 91:   /* Read mesh coordinates */
 92:   PetscCheck(numVertices >= 0, comm, PETSC_ERR_ARG_WRONG, "No nodes found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
 93:   PetscCall(PetscMalloc1(numVerticesLocal * spaceDim, &coordinates));
 94:   PetscCallExternal(MEDmeshNodeCoordinateAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, &vfilter, coordinates);
 95:   /* Read the types of entity sets in the mesh */
 96:   ngeo = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, MED_GEO_ALL, MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
 97:   PetscCheck(ngeo >= 1, comm, PETSC_ERR_ARG_WRONG, "No cells found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
 98:   PetscCheck(ngeo <= 2, comm, PETSC_ERR_ARG_WRONG, "Currently no support for hybrid meshes in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
 99:   PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 1, geotypename, &(geotype[0]));
100:   if (ngeo > 1) PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 2, geotypename, &(geotype[1]));
101:   else geotype[1] = 0;
102:   /* Determine topological dim and set ID for cells */
103:   cellID     = geotype[0] / 100 > geotype[1] / 100 ? 0 : 1;
104:   facetID    = geotype[0] / 100 > geotype[1] / 100 ? 1 : 0;
105:   meshDim    = geotype[cellID] / 100;
106:   numCorners = geotype[cellID] % 100;
107:   /* Partition cells */
108:   numCells = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
109:   PetscCall(PetscLayoutCreate(comm, &cLayout));
110:   PetscCall(PetscLayoutSetSize(cLayout, numCells));
111:   PetscCall(PetscLayoutSetBlockSize(cLayout, 1));
112:   PetscCall(PetscLayoutSetUp(cLayout));
113:   PetscCall(PetscLayoutGetRanges(cLayout, &crange));
114:   numCellsLocal = crange[rank + 1] - crange[rank];
115:   PetscCallExternal(MEDfilterBlockOfEntityCr, fileID, numCells, 1, numCorners, MED_ALL_CONSTITUENT, MED_FULL_INTERLACE, MED_COMPACT_STMODE, MED_NO_PROFILE, crange[rank] + 1, 1, numCellsLocal, 1, 1, &cfilter);
116:   /* Read cell connectivity */
117:   PetscCheck(numCells >= 0, comm, PETSC_ERR_ARG_WRONG, "No cells found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
118:   PetscCall(PetscMalloc1(numCellsLocal * numCorners, &medCellList));
119:   PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_NODAL, &cfilter, medCellList);
120:   PetscCheck(sizeof(med_int) <= sizeof(PetscInt), comm, PETSC_ERR_ARG_SIZ, "Size of PetscInt %zd less than  size of med_int %zd. Reconfigure PETSc --with-64-bit-indices=1", sizeof(PetscInt), sizeof(med_int));
121:   PetscCall(PetscMalloc1(numCellsLocal * numCorners, &cellList));
122:   for (i = 0; i < numCellsLocal * numCorners; i++) { cellList[i] = ((PetscInt)medCellList[i]) - 1; /* Correct entity counting */ }
123:   PetscCall(PetscFree(medCellList));
124:   /* Generate the DM */
125:   if (sizeof(med_float) == sizeof(PetscReal)) {
126:     vertcoords = (PetscReal *)coordinates;
127:   } else {
128:     PetscCall(PetscMalloc1(numVerticesLocal * spaceDim, &vertcoords));
129:     for (i = 0; i < numVerticesLocal * spaceDim; i++) vertcoords[i] = (PetscReal)coordinates[i];
130:   }
131:   /* Account for cell inversion */
132:   for (c = 0; c < numCellsLocal; ++c) {
133:     PetscInt *pcone = &cellList[c * numCorners];

135:     if (meshDim == 3) {
136:       /* Hexahedra are inverted */
137:       if (numCorners == 8) {
138:         PetscInt tmp = pcone[4 + 1];
139:         pcone[4 + 1] = pcone[4 + 3];
140:         pcone[4 + 3] = tmp;
141:       }
142:     }
143:   }
144:   PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, meshDim, numCellsLocal, numVerticesLocal, numVertices, numCorners, interpolate, cellList, spaceDim, vertcoords, &sfVertices, NULL, dm));
145:   if (sizeof(med_float) == sizeof(PetscReal)) {
146:     vertcoords = NULL;
147:   } else {
148:     PetscCall(PetscFree(vertcoords));
149:   }
150:   if (ngeo > 1) {
151:     PetscInt        numFacets = 0, numFacetsLocal, numFacetCorners, numFacetsRendezvous;
152:     PetscInt        c, f, v, vStart, joinSize, vertices[8];
153:     PetscInt       *facetList, *facetListRendezvous, *facetIDs, *facetIDsRendezvous, *facetListRemote, *facetIDsRemote;
154:     const PetscInt *frange, *join;
155:     PetscLayout     fLayout;
156:     med_filter      ffilter = MED_FILTER_INIT, fidfilter = MED_FILTER_INIT;
157:     PetscSection    facetSectionRemote, facetSectionIDsRemote;
158:     /* Partition facets */
159:     numFacets       = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
160:     numFacetCorners = geotype[facetID] % 100;
161:     PetscCall(PetscLayoutCreate(comm, &fLayout));
162:     PetscCall(PetscLayoutSetSize(fLayout, numFacets));
163:     PetscCall(PetscLayoutSetBlockSize(fLayout, 1));
164:     PetscCall(PetscLayoutSetUp(fLayout));
165:     PetscCall(PetscLayoutGetRanges(fLayout, &frange));
166:     numFacetsLocal = frange[rank + 1] - frange[rank];
167:     PetscCallExternal(MEDfilterBlockOfEntityCr, fileID, numFacets, 1, numFacetCorners, MED_ALL_CONSTITUENT, MED_FULL_INTERLACE, MED_COMPACT_STMODE, MED_NO_PROFILE, frange[rank] + 1, 1, numFacetsLocal, 1, 1, &ffilter);
168:     PetscCallExternal(MEDfilterBlockOfEntityCr, fileID, numFacets, 1, 1, MED_ALL_CONSTITUENT, MED_FULL_INTERLACE, MED_COMPACT_STMODE, MED_NO_PROFILE, frange[rank] + 1, 1, numFacetsLocal, 1, 1, &fidfilter);
169:     PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, NULL));
170:     PetscCall(PetscMalloc1(numFacetsLocal, &facetIDs));
171:     PetscCall(PetscMalloc1(numFacetsLocal * numFacetCorners, &facetList));

173:     /* Read facet connectivity */
174:     {
175:       med_int *medFacetList;

177:       PetscCall(PetscMalloc1(numFacetsLocal * numFacetCorners, &medFacetList));
178:       PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_NODAL, &ffilter, medFacetList);
179:       for (i = 0; i < numFacetsLocal * numFacetCorners; i++) { facetList[i] = ((PetscInt)medFacetList[i]) - 1; /* Correct entity counting */ }
180:       PetscCall(PetscFree(medFacetList));
181:     }

183:     /* Read facet IDs */
184:     {
185:       med_int *medFacetIDs;

187:       PetscCall(PetscMalloc1(numFacetsLocal, &medFacetIDs));
188:       PetscCallExternal(MEDmeshEntityAttributeAdvancedRd, fileID, meshname, MED_FAMILY_NUMBER, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], &fidfilter, medFacetIDs);
189:       for (i = 0; i < numFacetsLocal; i++) facetIDs[i] = (PetscInt)medFacetIDs[i];
190:       PetscCall(PetscFree(medFacetIDs));
191:     }

193:     /* Send facets and IDs to a rendezvous partition that is based on the initial vertex partitioning. */
194:     {
195:       PetscInt           r;
196:       DMLabel            lblFacetRendezvous, lblFacetMigration;
197:       PetscSection       facetSection, facetSectionRendezvous;
198:       PetscSF            sfProcess, sfFacetMigration;
199:       const PetscSFNode *remoteVertices;
200:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Facet Rendezvous", &lblFacetRendezvous));
201:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Facet Migration", &lblFacetMigration));
202:       PetscCall(PetscSFGetGraph(sfVertices, NULL, NULL, NULL, &remoteVertices));
203:       for (f = 0; f < numFacetsLocal; f++) {
204:         for (v = 0; v < numFacetCorners; v++) {
205:           /* Find vertex owner on rendezvous partition and mark in label */
206:           const PetscInt vertex = facetList[f * numFacetCorners + v];
207:           r                     = rank;
208:           while (vrange[r] > vertex) r--;
209:           while (vrange[r + 1] < vertex) r++;
210:           PetscCall(DMLabelSetValue(lblFacetRendezvous, f, r));
211:         }
212:       }
213:       /* Build a global process SF */
214:       PetscCall(PetscSFCreate(comm, &sfProcess));
215:       PetscCall(PetscSFSetGraphWithPattern(sfProcess, NULL, PETSCSF_PATTERN_ALLTOALL));
216:       PetscCall(PetscObjectSetName((PetscObject)sfProcess, "Process SF"));
217:       /* Convert facet rendezvous label into SF for migration */
218:       PetscCall(DMPlexPartitionLabelInvert(*dm, lblFacetRendezvous, sfProcess, lblFacetMigration));
219:       PetscCall(DMPlexPartitionLabelCreateSF(*dm, lblFacetMigration, &sfFacetMigration));
220:       /* Migrate facet connectivity data */
221:       PetscCall(PetscSectionCreate(comm, &facetSection));
222:       PetscCall(PetscSectionSetChart(facetSection, 0, numFacetsLocal));
223:       for (f = 0; f < numFacetsLocal; f++) PetscCall(PetscSectionSetDof(facetSection, f, numFacetCorners));
224:       PetscCall(PetscSectionSetUp(facetSection));
225:       PetscCall(PetscSectionCreate(comm, &facetSectionRendezvous));
226:       PetscCall(DMPlexDistributeData(*dm, sfFacetMigration, facetSection, MPIU_INT, facetList, facetSectionRendezvous, (void **)&facetListRendezvous));
227:       /* Migrate facet IDs */
228:       PetscCall(PetscSFGetGraph(sfFacetMigration, NULL, &numFacetsRendezvous, NULL, NULL));
229:       PetscCall(PetscMalloc1(numFacetsRendezvous, &facetIDsRendezvous));
230:       PetscCall(PetscSFBcastBegin(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE));
231:       PetscCall(PetscSFBcastEnd(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE));
232:       /* Clean up */
233:       PetscCall(DMLabelDestroy(&lblFacetRendezvous));
234:       PetscCall(DMLabelDestroy(&lblFacetMigration));
235:       PetscCall(PetscSFDestroy(&sfProcess));
236:       PetscCall(PetscSFDestroy(&sfFacetMigration));
237:       PetscCall(PetscSectionDestroy(&facetSection));
238:       PetscCall(PetscSectionDestroy(&facetSectionRendezvous));
239:     }

241:     /* On the rendevouz partition we build a vertex-wise section/array of facets and IDs. */
242:     {
243:       PetscInt               sizeVertexFacets, offset, sizeFacetIDsRemote;
244:       PetscInt              *vertexFacets, *vertexIdx, *vertexFacetIDs;
245:       PetscSection           facetSectionVertices, facetSectionIDs;
246:       ISLocalToGlobalMapping ltogVertexNumbering;
247:       PetscCall(PetscSectionCreate(comm, &facetSectionVertices));
248:       PetscCall(PetscSectionSetChart(facetSectionVertices, 0, numVerticesLocal));
249:       PetscCall(PetscSectionCreate(comm, &facetSectionIDs));
250:       PetscCall(PetscSectionSetChart(facetSectionIDs, 0, numVerticesLocal));
251:       for (f = 0; f < numFacetsRendezvous * numFacetCorners; f++) {
252:         const PetscInt vertex = facetListRendezvous[f];
253:         if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
254:           PetscCall(PetscSectionAddDof(facetSectionIDs, vertex - vrange[rank], 1));
255:           PetscCall(PetscSectionAddDof(facetSectionVertices, vertex - vrange[rank], numFacetCorners));
256:         }
257:       }
258:       PetscCall(PetscSectionSetUp(facetSectionVertices));
259:       PetscCall(PetscSectionSetUp(facetSectionIDs));
260:       PetscCall(PetscSectionGetStorageSize(facetSectionVertices, &sizeVertexFacets));
261:       PetscCall(PetscSectionGetStorageSize(facetSectionVertices, &sizeFacetIDsRemote));
262:       PetscCall(PetscMalloc1(sizeVertexFacets, &vertexFacets));
263:       PetscCall(PetscMalloc1(sizeFacetIDsRemote, &vertexFacetIDs));
264:       PetscCall(PetscCalloc1(numVerticesLocal, &vertexIdx));
265:       for (f = 0; f < numFacetsRendezvous; f++) {
266:         for (c = 0; c < numFacetCorners; c++) {
267:           const PetscInt vertex = facetListRendezvous[f * numFacetCorners + c];
268:           if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
269:             /* Flip facet connectivities and IDs to a vertex-wise layout */
270:             PetscCall(PetscSectionGetOffset(facetSectionVertices, vertex - vrange[rank], &offset));
271:             offset += vertexIdx[vertex - vrange[rank]] * numFacetCorners;
272:             PetscCall(PetscArraycpy(&(vertexFacets[offset]), &(facetListRendezvous[f * numFacetCorners]), numFacetCorners));
273:             PetscCall(PetscSectionGetOffset(facetSectionIDs, vertex - vrange[rank], &offset));
274:             offset += vertexIdx[vertex - vrange[rank]];
275:             vertexFacetIDs[offset] = facetIDsRendezvous[f];
276:             vertexIdx[vertex - vrange[rank]]++;
277:           }
278:         }
279:       }
280:       /* Distribute the vertex-wise facet connectivities over the vertexSF */
281:       PetscCall(PetscSectionCreate(comm, &facetSectionRemote));
282:       PetscCall(DMPlexDistributeData(*dm, sfVertices, facetSectionVertices, MPIU_INT, vertexFacets, facetSectionRemote, (void **)&facetListRemote));
283:       PetscCall(PetscSectionCreate(comm, &facetSectionIDsRemote));
284:       PetscCall(DMPlexDistributeData(*dm, sfVertices, facetSectionIDs, MPIU_INT, vertexFacetIDs, facetSectionIDsRemote, (void **)&facetIDsRemote));
285:       /* Convert facet connectivities to local vertex numbering */
286:       PetscCall(PetscSectionGetStorageSize(facetSectionRemote, &sizeVertexFacets));
287:       PetscCall(ISLocalToGlobalMappingCreateSF(sfVertices, vrange[rank], &ltogVertexNumbering));
288:       PetscCall(ISGlobalToLocalMappingApplyBlock(ltogVertexNumbering, IS_GTOLM_MASK, sizeVertexFacets, facetListRemote, NULL, facetListRemote));
289:       /* Clean up */
290:       PetscCall(PetscFree(vertexFacets));
291:       PetscCall(PetscFree(vertexIdx));
292:       PetscCall(PetscFree(vertexFacetIDs));
293:       PetscCall(PetscSectionDestroy(&facetSectionVertices));
294:       PetscCall(PetscSectionDestroy(&facetSectionIDs));
295:       PetscCall(ISLocalToGlobalMappingDestroy(&ltogVertexNumbering));
296:     }
297:     {
298:       PetscInt  offset, dof;
299:       DMLabel   lblFaceSets;
300:       PetscBool verticesLocal;
301:       /* Identify and mark facets locally with facet joins */
302:       PetscCall(DMCreateLabel(*dm, "Face Sets"));
303:       PetscCall(DMGetLabel(*dm, "Face Sets", &lblFaceSets));
304:       /* We need to set a new default value here, since -1 is a legitimate facet ID */
305:       PetscCall(DMLabelSetDefaultValue(lblFaceSets, -666666666));
306:       for (v = 0; v < numVerticesLocal; v++) {
307:         PetscCall(PetscSectionGetOffset(facetSectionRemote, v, &offset));
308:         PetscCall(PetscSectionGetDof(facetSectionRemote, v, &dof));
309:         for (f = 0; f < dof; f += numFacetCorners) {
310:           for (verticesLocal = PETSC_TRUE, c = 0; c < numFacetCorners; ++c) {
311:             if (facetListRemote[offset + f + c] < 0) {
312:               verticesLocal = PETSC_FALSE;
313:               break;
314:             }
315:             vertices[c] = vStart + facetListRemote[offset + f + c];
316:           }
317:           if (verticesLocal) {
318:             PetscCall(DMPlexGetFullJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join));
319:             if (joinSize == 1) PetscCall(DMLabelSetValue(lblFaceSets, join[0], facetIDsRemote[(offset + f) / numFacetCorners]));
320:             PetscCall(DMPlexRestoreJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join));
321:           }
322:         }
323:       }
324:     }
325:     PetscCall(PetscFree(facetList));
326:     PetscCall(PetscFree(facetListRendezvous));
327:     PetscCall(PetscFree(facetListRemote));
328:     PetscCall(PetscFree(facetIDs));
329:     PetscCall(PetscFree(facetIDsRendezvous));
330:     PetscCall(PetscFree(facetIDsRemote));
331:     PetscCall(PetscLayoutDestroy(&fLayout));
332:     PetscCall(PetscSectionDestroy(&facetSectionRemote));
333:     PetscCall(PetscSectionDestroy(&facetSectionIDsRemote));
334:   }
335:   PetscCallExternal(MEDfileClose, fileID);
336:   PetscCall(PetscFree(coordinates));
337:   PetscCall(PetscFree(cellList));
338:   PetscCall(PetscLayoutDestroy(&vLayout));
339:   PetscCall(PetscLayoutDestroy(&cLayout));
340:   PetscCall(PetscSFDestroy(&sfVertices));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: #else
343:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires Med mesh reader support. Reconfigure using --download-med");
344: #endif
345: }