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