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], &ltogVertexNumbering);
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(&ltogVertexNumbering);
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: }