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