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: Input Parameters:
14: + comm - The MPI communicator
15: . filename - Name of the .med file
16: - interpolate - Create faces and edges in the mesh
18: Output Parameter:
19: . dm - The `DM` object representing the mesh
21: Level: beginner
23: References:
24: . * - https://www.salome-platform.org/user-section/about/med, http://docs.salome-platform.org/latest/MED_index.html
26: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexCreateFromFile()`, `DMPlexCreateGmsh()`, `DMPlexCreate()`
27: @*/
28: PetscErrorCode DMPlexCreateMedFromFile(MPI_Comm comm, const char filename[], PetscBool interpolate, DM *dm)
29: {
30: PetscMPIInt rank, size;
31: #if defined(PETSC_HAVE_MED)
32: med_idt fileID;
33: PetscInt i, ngeo, spaceDim, meshDim;
34: PetscInt numVertices = 0, numCells = 0, c, numCorners, numCellsLocal, numVerticesLocal;
35: med_int *medCellList;
36: PetscInt *cellList;
37: med_float *coordinates = NULL;
38: PetscReal *vertcoords = NULL;
39: PetscLayout vLayout, cLayout;
40: const PetscInt *vrange, *crange;
41: PetscSF sfVertices;
42: char *axisname, *unitname, meshname[MED_NAME_SIZE + 1], geotypename[MED_NAME_SIZE + 1];
43: char meshdescription[MED_COMMENT_SIZE + 1], dtunit[MED_SNAME_SIZE + 1];
44: med_sorting_type sortingtype;
45: med_mesh_type meshtype;
46: med_axis_type axistype;
47: med_bool coordinatechangement, geotransformation, hdfok, medok;
48: med_geometry_type geotype[2], cellID, facetID;
49: med_filter vfilter = MED_FILTER_INIT;
50: med_filter cfilter = MED_FILTER_INIT;
51: med_err mederr;
52: med_int major, minor, release;
53: #endif
55: PetscFunctionBegin;
56: PetscCallMPI(MPI_Comm_rank(comm, &rank));
57: PetscCallMPI(MPI_Comm_size(comm, &size));
58: #if defined(PETSC_HAVE_MED)
59: mederr = MEDfileCompatibility(filename, &hdfok, &medok);
60: PetscCheck(!mederr, comm, PETSC_ERR_ARG_WRONG, "Cannot determine MED file compatibility: %s", filename);
61: PetscCheck(hdfok, comm, PETSC_ERR_ARG_WRONG, "Not a compatible HDF format: %s", filename);
62: PetscCheck(medok, comm, PETSC_ERR_ARG_WRONG, "Not a compatible MED format: %s", filename);
64: fileID = MEDfileOpen(filename, MED_ACC_RDONLY);
65: PetscCheck(fileID >= 0, comm, PETSC_ERR_ARG_WRONG, "Unable to open .med mesh file: %s", filename);
66: mederr = MEDfileNumVersionRd(fileID, &major, &minor, &release);
67: 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);
68: spaceDim = MEDmeshnAxis(fileID, 1);
69: 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);
70: /* Read general mesh information */
71: PetscCall(PetscMalloc1(MED_SNAME_SIZE * spaceDim + 1, &axisname));
72: PetscCall(PetscMalloc1(MED_SNAME_SIZE * spaceDim + 1, &unitname));
73: {
74: med_int medMeshDim, medNstep;
75: med_int medSpaceDim = spaceDim;
77: PetscCallExternal(MEDmeshInfo, fileID, 1, meshname, &medSpaceDim, &medMeshDim, &meshtype, meshdescription, dtunit, &sortingtype, &medNstep, &axistype, axisname, unitname);
78: spaceDim = medSpaceDim;
79: meshDim = medMeshDim;
80: }
81: PetscCall(PetscFree(axisname));
82: PetscCall(PetscFree(unitname));
83: /* Partition mesh coordinates */
84: numVertices = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_NODE, MED_NO_GEOTYPE, MED_COORDINATE, MED_NO_CMODE, &coordinatechangement, &geotransformation);
85: PetscCall(PetscLayoutCreate(comm, &vLayout));
86: PetscCall(PetscLayoutSetSize(vLayout, numVertices));
87: PetscCall(PetscLayoutSetBlockSize(vLayout, 1));
88: PetscCall(PetscLayoutSetUp(vLayout));
89: PetscCall(PetscLayoutGetRanges(vLayout, &vrange));
90: numVerticesLocal = vrange[rank + 1] - vrange[rank];
91: 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);
92: /* Read mesh coordinates */
93: PetscCheck(numVertices >= 0, comm, PETSC_ERR_ARG_WRONG, "No nodes found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
94: PetscCall(PetscMalloc1(numVerticesLocal * spaceDim, &coordinates));
95: PetscCallExternal(MEDmeshNodeCoordinateAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, &vfilter, coordinates);
96: /* Read the types of entity sets in the mesh */
97: ngeo = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, MED_GEO_ALL, MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
98: PetscCheck(ngeo >= 1, comm, PETSC_ERR_ARG_WRONG, "No cells found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
99: 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);
100: PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 1, geotypename, &(geotype[0]));
101: if (ngeo > 1) PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 2, geotypename, &(geotype[1]));
102: else geotype[1] = 0;
103: /* Determine topological dim and set ID for cells */
104: cellID = geotype[0] / 100 > geotype[1] / 100 ? 0 : 1;
105: facetID = geotype[0] / 100 > geotype[1] / 100 ? 1 : 0;
106: meshDim = geotype[cellID] / 100;
107: numCorners = geotype[cellID] % 100;
108: /* Partition cells */
109: numCells = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
110: PetscCall(PetscLayoutCreate(comm, &cLayout));
111: PetscCall(PetscLayoutSetSize(cLayout, numCells));
112: PetscCall(PetscLayoutSetBlockSize(cLayout, 1));
113: PetscCall(PetscLayoutSetUp(cLayout));
114: PetscCall(PetscLayoutGetRanges(cLayout, &crange));
115: numCellsLocal = crange[rank + 1] - crange[rank];
116: 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);
117: /* Read cell connectivity */
118: PetscCheck(numCells >= 0, comm, PETSC_ERR_ARG_WRONG, "No cells found in .med v%d.%d.%d mesh file: %s", major, minor, release, filename);
119: PetscCall(PetscMalloc1(numCellsLocal * numCorners, &medCellList));
120: PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_NODAL, &cfilter, medCellList);
121: 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));
122: PetscCall(PetscMalloc1(numCellsLocal * numCorners, &cellList));
123: for (i = 0; i < numCellsLocal * numCorners; i++) { cellList[i] = ((PetscInt)medCellList[i]) - 1; /* Correct entity counting */ }
124: PetscCall(PetscFree(medCellList));
125: /* Generate the DM */
126: if (sizeof(med_float) == sizeof(PetscReal)) {
127: vertcoords = (PetscReal *)coordinates;
128: } else {
129: PetscCall(PetscMalloc1(numVerticesLocal * spaceDim, &vertcoords));
130: for (i = 0; i < numVerticesLocal * spaceDim; i++) vertcoords[i] = (PetscReal)coordinates[i];
131: }
132: /* Account for cell inversion */
133: for (c = 0; c < numCellsLocal; ++c) {
134: PetscInt *pcone = &cellList[c * numCorners];
136: if (meshDim == 3) {
137: /* Hexahedra are inverted */
138: if (numCorners == 8) {
139: PetscInt tmp = pcone[4 + 1];
140: pcone[4 + 1] = pcone[4 + 3];
141: pcone[4 + 3] = tmp;
142: }
143: }
144: }
145: PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, meshDim, numCellsLocal, numVerticesLocal, numVertices, numCorners, interpolate, cellList, spaceDim, vertcoords, &sfVertices, NULL, dm));
146: if (sizeof(med_float) == sizeof(PetscReal)) {
147: vertcoords = NULL;
148: } else {
149: PetscCall(PetscFree(vertcoords));
150: }
151: if (ngeo > 1) {
152: PetscInt numFacets = 0, numFacetsLocal, numFacetCorners, numFacetsRendezvous;
153: PetscInt c, f, v, vStart, joinSize, vertices[8];
154: PetscInt *facetList, *facetListRendezvous, *facetIDs, *facetIDsRendezvous, *facetListRemote, *facetIDsRemote;
155: const PetscInt *frange, *join;
156: PetscLayout fLayout;
157: med_filter ffilter = MED_FILTER_INIT, fidfilter = MED_FILTER_INIT;
158: PetscSection facetSectionRemote, facetSectionIDsRemote;
159: /* Partition facets */
160: numFacets = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
161: numFacetCorners = geotype[facetID] % 100;
162: PetscCall(PetscLayoutCreate(comm, &fLayout));
163: PetscCall(PetscLayoutSetSize(fLayout, numFacets));
164: PetscCall(PetscLayoutSetBlockSize(fLayout, 1));
165: PetscCall(PetscLayoutSetUp(fLayout));
166: PetscCall(PetscLayoutGetRanges(fLayout, &frange));
167: numFacetsLocal = frange[rank + 1] - frange[rank];
168: 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);
169: 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);
170: PetscCall(DMPlexGetDepthStratum(*dm, 0, &vStart, NULL));
171: PetscCall(PetscMalloc1(numFacetsLocal, &facetIDs));
172: PetscCall(PetscMalloc1(numFacetsLocal * numFacetCorners, &facetList));
174: /* Read facet connectivity */
175: {
176: med_int *medFacetList;
178: PetscCall(PetscMalloc1(numFacetsLocal * numFacetCorners, &medFacetList));
179: PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_NODAL, &ffilter, medFacetList);
180: for (i = 0; i < numFacetsLocal * numFacetCorners; i++) { facetList[i] = ((PetscInt)medFacetList[i]) - 1; /* Correct entity counting */ }
181: PetscCall(PetscFree(medFacetList));
182: }
184: /* Read facet IDs */
185: {
186: med_int *medFacetIDs;
188: PetscCall(PetscMalloc1(numFacetsLocal, &medFacetIDs));
189: PetscCallExternal(MEDmeshEntityAttributeAdvancedRd, fileID, meshname, MED_FAMILY_NUMBER, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], &fidfilter, medFacetIDs);
190: for (i = 0; i < numFacetsLocal; i++) facetIDs[i] = (PetscInt)medFacetIDs[i];
191: PetscCall(PetscFree(medFacetIDs));
192: }
194: /* Send facets and IDs to a rendezvous partition that is based on the initial vertex partitioning. */
195: {
196: PetscInt r;
197: DMLabel lblFacetRendezvous, lblFacetMigration;
198: PetscSection facetSection, facetSectionRendezvous;
199: PetscSF sfProcess, sfFacetMigration;
200: const PetscSFNode *remoteVertices;
201: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Facet Rendezvous", &lblFacetRendezvous));
202: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "Facet Migration", &lblFacetMigration));
203: PetscCall(PetscSFGetGraph(sfVertices, NULL, NULL, NULL, &remoteVertices));
204: for (f = 0; f < numFacetsLocal; f++) {
205: for (v = 0; v < numFacetCorners; v++) {
206: /* Find vertex owner on rendezvous partition and mark in label */
207: const PetscInt vertex = facetList[f * numFacetCorners + v];
208: r = rank;
209: while (vrange[r] > vertex) r--;
210: while (vrange[r + 1] < vertex) r++;
211: PetscCall(DMLabelSetValue(lblFacetRendezvous, f, r));
212: }
213: }
214: /* Build a global process SF */
215: PetscCall(PetscSFCreate(comm, &sfProcess));
216: PetscCall(PetscSFSetGraphWithPattern(sfProcess, NULL, PETSCSF_PATTERN_ALLTOALL));
217: PetscCall(PetscObjectSetName((PetscObject)sfProcess, "Process SF"));
218: /* Convert facet rendezvous label into SF for migration */
219: PetscCall(DMPlexPartitionLabelInvert(*dm, lblFacetRendezvous, sfProcess, lblFacetMigration));
220: PetscCall(DMPlexPartitionLabelCreateSF(*dm, lblFacetMigration, &sfFacetMigration));
221: /* Migrate facet connectivity data */
222: PetscCall(PetscSectionCreate(comm, &facetSection));
223: PetscCall(PetscSectionSetChart(facetSection, 0, numFacetsLocal));
224: for (f = 0; f < numFacetsLocal; f++) PetscCall(PetscSectionSetDof(facetSection, f, numFacetCorners));
225: PetscCall(PetscSectionSetUp(facetSection));
226: PetscCall(PetscSectionCreate(comm, &facetSectionRendezvous));
227: PetscCall(DMPlexDistributeData(*dm, sfFacetMigration, facetSection, MPIU_INT, facetList, facetSectionRendezvous, (void **)&facetListRendezvous));
228: /* Migrate facet IDs */
229: PetscCall(PetscSFGetGraph(sfFacetMigration, NULL, &numFacetsRendezvous, NULL, NULL));
230: PetscCall(PetscMalloc1(numFacetsRendezvous, &facetIDsRendezvous));
231: PetscCall(PetscSFBcastBegin(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE));
232: PetscCall(PetscSFBcastEnd(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE));
233: /* Clean up */
234: PetscCall(DMLabelDestroy(&lblFacetRendezvous));
235: PetscCall(DMLabelDestroy(&lblFacetMigration));
236: PetscCall(PetscSFDestroy(&sfProcess));
237: PetscCall(PetscSFDestroy(&sfFacetMigration));
238: PetscCall(PetscSectionDestroy(&facetSection));
239: PetscCall(PetscSectionDestroy(&facetSectionRendezvous));
240: }
242: /* On the rendevouz partition we build a vertex-wise section/array of facets and IDs. */
243: {
244: PetscInt sizeVertexFacets, offset, sizeFacetIDsRemote;
245: PetscInt *vertexFacets, *vertexIdx, *vertexFacetIDs;
246: PetscSection facetSectionVertices, facetSectionIDs;
247: ISLocalToGlobalMapping ltogVertexNumbering;
248: PetscCall(PetscSectionCreate(comm, &facetSectionVertices));
249: PetscCall(PetscSectionSetChart(facetSectionVertices, 0, numVerticesLocal));
250: PetscCall(PetscSectionCreate(comm, &facetSectionIDs));
251: PetscCall(PetscSectionSetChart(facetSectionIDs, 0, numVerticesLocal));
252: for (f = 0; f < numFacetsRendezvous * numFacetCorners; f++) {
253: const PetscInt vertex = facetListRendezvous[f];
254: if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
255: PetscCall(PetscSectionAddDof(facetSectionIDs, vertex - vrange[rank], 1));
256: PetscCall(PetscSectionAddDof(facetSectionVertices, vertex - vrange[rank], numFacetCorners));
257: }
258: }
259: PetscCall(PetscSectionSetUp(facetSectionVertices));
260: PetscCall(PetscSectionSetUp(facetSectionIDs));
261: PetscCall(PetscSectionGetStorageSize(facetSectionVertices, &sizeVertexFacets));
262: PetscCall(PetscSectionGetStorageSize(facetSectionVertices, &sizeFacetIDsRemote));
263: PetscCall(PetscMalloc1(sizeVertexFacets, &vertexFacets));
264: PetscCall(PetscMalloc1(sizeFacetIDsRemote, &vertexFacetIDs));
265: PetscCall(PetscCalloc1(numVerticesLocal, &vertexIdx));
266: for (f = 0; f < numFacetsRendezvous; f++) {
267: for (c = 0; c < numFacetCorners; c++) {
268: const PetscInt vertex = facetListRendezvous[f * numFacetCorners + c];
269: if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
270: /* Flip facet connectivities and IDs to a vertex-wise layout */
271: PetscCall(PetscSectionGetOffset(facetSectionVertices, vertex - vrange[rank], &offset));
272: offset += vertexIdx[vertex - vrange[rank]] * numFacetCorners;
273: PetscCall(PetscArraycpy(&(vertexFacets[offset]), &(facetListRendezvous[f * numFacetCorners]), numFacetCorners));
274: PetscCall(PetscSectionGetOffset(facetSectionIDs, vertex - vrange[rank], &offset));
275: offset += vertexIdx[vertex - vrange[rank]];
276: vertexFacetIDs[offset] = facetIDsRendezvous[f];
277: vertexIdx[vertex - vrange[rank]]++;
278: }
279: }
280: }
281: /* Distribute the vertex-wise facet connectivities over the vertexSF */
282: PetscCall(PetscSectionCreate(comm, &facetSectionRemote));
283: PetscCall(DMPlexDistributeData(*dm, sfVertices, facetSectionVertices, MPIU_INT, vertexFacets, facetSectionRemote, (void **)&facetListRemote));
284: PetscCall(PetscSectionCreate(comm, &facetSectionIDsRemote));
285: PetscCall(DMPlexDistributeData(*dm, sfVertices, facetSectionIDs, MPIU_INT, vertexFacetIDs, facetSectionIDsRemote, (void **)&facetIDsRemote));
286: /* Convert facet connectivities to local vertex numbering */
287: PetscCall(PetscSectionGetStorageSize(facetSectionRemote, &sizeVertexFacets));
288: PetscCall(ISLocalToGlobalMappingCreateSF(sfVertices, vrange[rank], <ogVertexNumbering));
289: PetscCall(ISGlobalToLocalMappingApplyBlock(ltogVertexNumbering, IS_GTOLM_MASK, sizeVertexFacets, facetListRemote, NULL, facetListRemote));
290: /* Clean up */
291: PetscCall(PetscFree(vertexFacets));
292: PetscCall(PetscFree(vertexIdx));
293: PetscCall(PetscFree(vertexFacetIDs));
294: PetscCall(PetscSectionDestroy(&facetSectionVertices));
295: PetscCall(PetscSectionDestroy(&facetSectionIDs));
296: PetscCall(ISLocalToGlobalMappingDestroy(<ogVertexNumbering));
297: }
298: {
299: PetscInt offset, dof;
300: DMLabel lblFaceSets;
301: PetscBool verticesLocal;
302: /* Identify and mark facets locally with facet joins */
303: PetscCall(DMCreateLabel(*dm, "Face Sets"));
304: PetscCall(DMGetLabel(*dm, "Face Sets", &lblFaceSets));
305: /* We need to set a new default value here, since -1 is a legitimate facet ID */
306: PetscCall(DMLabelSetDefaultValue(lblFaceSets, -666666666));
307: for (v = 0; v < numVerticesLocal; v++) {
308: PetscCall(PetscSectionGetOffset(facetSectionRemote, v, &offset));
309: PetscCall(PetscSectionGetDof(facetSectionRemote, v, &dof));
310: for (f = 0; f < dof; f += numFacetCorners) {
311: for (verticesLocal = PETSC_TRUE, c = 0; c < numFacetCorners; ++c) {
312: if (facetListRemote[offset + f + c] < 0) {
313: verticesLocal = PETSC_FALSE;
314: break;
315: }
316: vertices[c] = vStart + facetListRemote[offset + f + c];
317: }
318: if (verticesLocal) {
319: PetscCall(DMPlexGetFullJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join));
320: if (joinSize == 1) PetscCall(DMLabelSetValue(lblFaceSets, join[0], facetIDsRemote[(offset + f) / numFacetCorners]));
321: PetscCall(DMPlexRestoreJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join));
322: }
323: }
324: }
325: }
326: PetscCall(PetscFree(facetList));
327: PetscCall(PetscFree(facetListRendezvous));
328: PetscCall(PetscFree(facetListRemote));
329: PetscCall(PetscFree(facetIDs));
330: PetscCall(PetscFree(facetIDsRendezvous));
331: PetscCall(PetscFree(facetIDsRemote));
332: PetscCall(PetscLayoutDestroy(&fLayout));
333: PetscCall(PetscSectionDestroy(&facetSectionRemote));
334: PetscCall(PetscSectionDestroy(&facetSectionIDsRemote));
335: }
336: PetscCallExternal(MEDfileClose, fileID);
337: PetscCall(PetscFree(coordinates));
338: PetscCall(PetscFree(cellList));
339: PetscCall(PetscLayoutDestroy(&vLayout));
340: PetscCall(PetscLayoutDestroy(&cLayout));
341: PetscCall(PetscSFDestroy(&sfVertices));
342: PetscFunctionReturn(PETSC_SUCCESS);
343: #else
344: SETERRQ(comm, PETSC_ERR_SUP, "This method requires Med mesh reader support. Reconfigure using --download-med");
345: #endif
346: }