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