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