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:     PetscCallExternal(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, MED_COORDINATE, MED_NO_CMODE, &coordinatechangement, &geotransformation);
 80:   PetscLayoutCreate(comm, &vLayout);
 81:   PetscLayoutSetSize(vLayout, numVertices);
 82:   PetscLayoutSetBlockSize(vLayout, 1);
 83:   PetscLayoutSetUp(vLayout);
 84:   PetscLayoutGetRanges(vLayout, &vrange);
 85:   numVerticesLocal = vrange[rank + 1] - vrange[rank];
 86:   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);
 87:   /* Read mesh coordinates */
 89:   PetscMalloc1(numVerticesLocal * spaceDim, &coordinates);
 90:   PetscCallExternal(MEDmeshNodeCoordinateAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, &vfilter, coordinates);
 91:   /* Read the types of entity sets in the mesh */
 92:   ngeo = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, MED_GEO_ALL, MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
 95:   PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 1, geotypename, &(geotype[0]));
 96:   if (ngeo > 1) PetscCallExternal(MEDmeshEntityInfo, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, 2, geotypename, &(geotype[1]));
 97:   else geotype[1] = 0;
 98:   /* Determine topological dim and set ID for cells */
 99:   cellID     = geotype[0] / 100 > geotype[1] / 100 ? 0 : 1;
100:   facetID    = geotype[0] / 100 > geotype[1] / 100 ? 1 : 0;
101:   meshDim    = geotype[cellID] / 100;
102:   numCorners = geotype[cellID] % 100;
103:   /* Partition cells */
104:   numCells = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
105:   PetscLayoutCreate(comm, &cLayout);
106:   PetscLayoutSetSize(cLayout, numCells);
107:   PetscLayoutSetBlockSize(cLayout, 1);
108:   PetscLayoutSetUp(cLayout);
109:   PetscLayoutGetRanges(cLayout, &crange);
110:   numCellsLocal = crange[rank + 1] - crange[rank];
111:   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);
112:   /* Read cell connectivity */
114:   PetscMalloc1(numCellsLocal * numCorners, &medCellList);
115:   PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[cellID], MED_NODAL, &cfilter, medCellList);
117:   PetscMalloc1(numCellsLocal * numCorners, &cellList);
118:   for (i = 0; i < numCellsLocal * numCorners; i++) { cellList[i] = ((PetscInt)medCellList[i]) - 1; /* Correct entity counting */ }
119:   PetscFree(medCellList);
120:   /* Generate the DM */
121:   if (sizeof(med_float) == sizeof(PetscReal)) {
122:     vertcoords = (PetscReal *)coordinates;
123:   } else {
124:     PetscMalloc1(numVerticesLocal * spaceDim, &vertcoords);
125:     for (i = 0; i < numVerticesLocal * spaceDim; i++) vertcoords[i] = (PetscReal)coordinates[i];
126:   }
127:   /* Account for cell inversion */
128:   for (c = 0; c < numCellsLocal; ++c) {
129:     PetscInt *pcone = &cellList[c * numCorners];

131:     if (meshDim == 3) {
132:       /* Hexahedra are inverted */
133:       if (numCorners == 8) {
134:         PetscInt tmp = pcone[4 + 1];
135:         pcone[4 + 1] = pcone[4 + 3];
136:         pcone[4 + 3] = tmp;
137:       }
138:     }
139:   }
140:   DMPlexCreateFromCellListParallelPetsc(comm, meshDim, numCellsLocal, numVerticesLocal, numVertices, numCorners, interpolate, cellList, spaceDim, vertcoords, &sfVertices, NULL, dm);
141:   if (sizeof(med_float) == sizeof(PetscReal)) {
142:     vertcoords = NULL;
143:   } else {
144:     PetscFree(vertcoords);
145:   }
146:   if (ngeo > 1) {
147:     PetscInt        numFacets = 0, numFacetsLocal, numFacetCorners, numFacetsRendezvous;
148:     PetscInt        c, f, v, vStart, joinSize, vertices[8];
149:     PetscInt       *facetList, *facetListRendezvous, *facetIDs, *facetIDsRendezvous, *facetListRemote, *facetIDsRemote;
150:     const PetscInt *frange, *join;
151:     PetscLayout     fLayout;
152:     med_filter      ffilter = MED_FILTER_INIT, fidfilter = MED_FILTER_INIT;
153:     PetscSection    facetSectionRemote, facetSectionIDsRemote;
154:     /* Partition facets */
155:     numFacets       = MEDmeshnEntity(fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_CONNECTIVITY, MED_NODAL, &coordinatechangement, &geotransformation);
156:     numFacetCorners = geotype[facetID] % 100;
157:     PetscLayoutCreate(comm, &fLayout);
158:     PetscLayoutSetSize(fLayout, numFacets);
159:     PetscLayoutSetBlockSize(fLayout, 1);
160:     PetscLayoutSetUp(fLayout);
161:     PetscLayoutGetRanges(fLayout, &frange);
162:     numFacetsLocal = frange[rank + 1] - frange[rank];
163:     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);
164:     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);
165:     DMPlexGetDepthStratum(*dm, 0, &vStart, NULL);
166:     PetscMalloc1(numFacetsLocal, &facetIDs);
167:     PetscMalloc1(numFacetsLocal * numFacetCorners, &facetList);

169:     /* Read facet connectivity */
170:     {
171:       med_int *medFacetList;

173:       PetscMalloc1(numFacetsLocal * numFacetCorners, &medFacetList);
174:       PetscCallExternal(MEDmeshElementConnectivityAdvancedRd, fileID, meshname, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], MED_NODAL, &ffilter, medFacetList);
175:       for (i = 0; i < numFacetsLocal * numFacetCorners; i++) { facetList[i] = ((PetscInt)medFacetList[i]) - 1; /* Correct entity counting */ }
176:       PetscFree(medFacetList);
177:     }

179:     /* Read facet IDs */
180:     {
181:       med_int *medFacetIDs;

183:       PetscMalloc1(numFacetsLocal, &medFacetIDs);
184:       PetscCallExternal(MEDmeshEntityAttributeAdvancedRd, fileID, meshname, MED_FAMILY_NUMBER, MED_NO_DT, MED_NO_IT, MED_CELL, geotype[facetID], &fidfilter, medFacetIDs);
185:       for (i = 0; i < numFacetsLocal; i++) facetIDs[i] = (PetscInt)medFacetIDs[i];
186:       PetscFree(medFacetIDs);
187:     }

189:     /* Send facets and IDs to a rendezvous partition that is based on the initial vertex partitioning. */
190:     {
191:       PetscInt           r;
192:       DMLabel            lblFacetRendezvous, lblFacetMigration;
193:       PetscSection       facetSection, facetSectionRendezvous;
194:       PetscSF            sfProcess, sfFacetMigration;
195:       const PetscSFNode *remoteVertices;
196:       DMLabelCreate(PETSC_COMM_SELF, "Facet Rendezvous", &lblFacetRendezvous);
197:       DMLabelCreate(PETSC_COMM_SELF, "Facet Migration", &lblFacetMigration);
198:       PetscSFGetGraph(sfVertices, NULL, NULL, NULL, &remoteVertices);
199:       for (f = 0; f < numFacetsLocal; f++) {
200:         for (v = 0; v < numFacetCorners; v++) {
201:           /* Find vertex owner on rendezvous partition and mark in label */
202:           const PetscInt vertex = facetList[f * numFacetCorners + v];
203:           r                     = rank;
204:           while (vrange[r] > vertex) r--;
205:           while (vrange[r + 1] < vertex) r++;
206:           DMLabelSetValue(lblFacetRendezvous, f, r);
207:         }
208:       }
209:       /* Build a global process SF */
210:       PetscSFCreate(comm, &sfProcess);
211:       PetscSFSetGraphWithPattern(sfProcess, NULL, PETSCSF_PATTERN_ALLTOALL);
212:       PetscObjectSetName((PetscObject)sfProcess, "Process SF");
213:       /* Convert facet rendezvous label into SF for migration */
214:       DMPlexPartitionLabelInvert(*dm, lblFacetRendezvous, sfProcess, lblFacetMigration);
215:       DMPlexPartitionLabelCreateSF(*dm, lblFacetMigration, &sfFacetMigration);
216:       /* Migrate facet connectivity data */
217:       PetscSectionCreate(comm, &facetSection);
218:       PetscSectionSetChart(facetSection, 0, numFacetsLocal);
219:       for (f = 0; f < numFacetsLocal; f++) PetscSectionSetDof(facetSection, f, numFacetCorners);
220:       PetscSectionSetUp(facetSection);
221:       PetscSectionCreate(comm, &facetSectionRendezvous);
222:       DMPlexDistributeData(*dm, sfFacetMigration, facetSection, MPIU_INT, facetList, facetSectionRendezvous, (void **)&facetListRendezvous);
223:       /* Migrate facet IDs */
224:       PetscSFGetGraph(sfFacetMigration, NULL, &numFacetsRendezvous, NULL, NULL);
225:       PetscMalloc1(numFacetsRendezvous, &facetIDsRendezvous);
226:       PetscSFBcastBegin(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE);
227:       PetscSFBcastEnd(sfFacetMigration, MPIU_INT, facetIDs, facetIDsRendezvous, MPI_REPLACE);
228:       /* Clean up */
229:       DMLabelDestroy(&lblFacetRendezvous);
230:       DMLabelDestroy(&lblFacetMigration);
231:       PetscSFDestroy(&sfProcess);
232:       PetscSFDestroy(&sfFacetMigration);
233:       PetscSectionDestroy(&facetSection);
234:       PetscSectionDestroy(&facetSectionRendezvous);
235:     }

237:     /* On the rendevouz partition we build a vertex-wise section/array of facets and IDs. */
238:     {
239:       PetscInt               sizeVertexFacets, offset, sizeFacetIDsRemote;
240:       PetscInt              *vertexFacets, *vertexIdx, *vertexFacetIDs;
241:       PetscSection           facetSectionVertices, facetSectionIDs;
242:       ISLocalToGlobalMapping ltogVertexNumbering;
243:       PetscSectionCreate(comm, &facetSectionVertices);
244:       PetscSectionSetChart(facetSectionVertices, 0, numVerticesLocal);
245:       PetscSectionCreate(comm, &facetSectionIDs);
246:       PetscSectionSetChart(facetSectionIDs, 0, numVerticesLocal);
247:       for (f = 0; f < numFacetsRendezvous * numFacetCorners; f++) {
248:         const PetscInt vertex = facetListRendezvous[f];
249:         if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
250:           PetscSectionAddDof(facetSectionIDs, vertex - vrange[rank], 1);
251:           PetscSectionAddDof(facetSectionVertices, vertex - vrange[rank], numFacetCorners);
252:         }
253:       }
254:       PetscSectionSetUp(facetSectionVertices);
255:       PetscSectionSetUp(facetSectionIDs);
256:       PetscSectionGetStorageSize(facetSectionVertices, &sizeVertexFacets);
257:       PetscSectionGetStorageSize(facetSectionVertices, &sizeFacetIDsRemote);
258:       PetscMalloc1(sizeVertexFacets, &vertexFacets);
259:       PetscMalloc1(sizeFacetIDsRemote, &vertexFacetIDs);
260:       PetscCalloc1(numVerticesLocal, &vertexIdx);
261:       for (f = 0; f < numFacetsRendezvous; f++) {
262:         for (c = 0; c < numFacetCorners; c++) {
263:           const PetscInt vertex = facetListRendezvous[f * numFacetCorners + c];
264:           if (vrange[rank] <= vertex && vertex < vrange[rank + 1]) {
265:             /* Flip facet connectivities and IDs to a vertex-wise layout */
266:             PetscSectionGetOffset(facetSectionVertices, vertex - vrange[rank], &offset);
267:             offset += vertexIdx[vertex - vrange[rank]] * numFacetCorners;
268:             PetscArraycpy(&(vertexFacets[offset]), &(facetListRendezvous[f * numFacetCorners]), numFacetCorners);
269:             PetscSectionGetOffset(facetSectionIDs, vertex - vrange[rank], &offset);
270:             offset += vertexIdx[vertex - vrange[rank]];
271:             vertexFacetIDs[offset] = facetIDsRendezvous[f];
272:             vertexIdx[vertex - vrange[rank]]++;
273:           }
274:         }
275:       }
276:       /* Distribute the vertex-wise facet connectivities over the vertexSF */
277:       PetscSectionCreate(comm, &facetSectionRemote);
278:       DMPlexDistributeData(*dm, sfVertices, facetSectionVertices, MPIU_INT, vertexFacets, facetSectionRemote, (void **)&facetListRemote);
279:       PetscSectionCreate(comm, &facetSectionIDsRemote);
280:       DMPlexDistributeData(*dm, sfVertices, facetSectionIDs, MPIU_INT, vertexFacetIDs, facetSectionIDsRemote, (void **)&facetIDsRemote);
281:       /* Convert facet connectivities to local vertex numbering */
282:       PetscSectionGetStorageSize(facetSectionRemote, &sizeVertexFacets);
283:       ISLocalToGlobalMappingCreateSF(sfVertices, vrange[rank], &ltogVertexNumbering);
284:       ISGlobalToLocalMappingApplyBlock(ltogVertexNumbering, IS_GTOLM_MASK, sizeVertexFacets, facetListRemote, NULL, facetListRemote);
285:       /* Clean up */
286:       PetscFree(vertexFacets);
287:       PetscFree(vertexIdx);
288:       PetscFree(vertexFacetIDs);
289:       PetscSectionDestroy(&facetSectionVertices);
290:       PetscSectionDestroy(&facetSectionIDs);
291:       ISLocalToGlobalMappingDestroy(&ltogVertexNumbering);
292:     }
293:     {
294:       PetscInt  offset, dof;
295:       DMLabel   lblFaceSets;
296:       PetscBool verticesLocal;
297:       /* Identify and mark facets locally with facet joins */
298:       DMCreateLabel(*dm, "Face Sets");
299:       DMGetLabel(*dm, "Face Sets", &lblFaceSets);
300:       /* We need to set a new default value here, since -1 is a legitimate facet ID */
301:       DMLabelSetDefaultValue(lblFaceSets, -666666666);
302:       for (v = 0; v < numVerticesLocal; v++) {
303:         PetscSectionGetOffset(facetSectionRemote, v, &offset);
304:         PetscSectionGetDof(facetSectionRemote, v, &dof);
305:         for (f = 0; f < dof; f += numFacetCorners) {
306:           for (verticesLocal = PETSC_TRUE, c = 0; c < numFacetCorners; ++c) {
307:             if (facetListRemote[offset + f + c] < 0) {
308:               verticesLocal = PETSC_FALSE;
309:               break;
310:             }
311:             vertices[c] = vStart + facetListRemote[offset + f + c];
312:           }
313:           if (verticesLocal) {
314:             DMPlexGetFullJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join);
315:             if (joinSize == 1) DMLabelSetValue(lblFaceSets, join[0], facetIDsRemote[(offset + f) / numFacetCorners]);
316:             DMPlexRestoreJoin(*dm, numFacetCorners, (const PetscInt *)vertices, &joinSize, &join);
317:           }
318:         }
319:       }
320:     }
321:     PetscFree(facetList);
322:     PetscFree(facetListRendezvous);
323:     PetscFree(facetListRemote);
324:     PetscFree(facetIDs);
325:     PetscFree(facetIDsRendezvous);
326:     PetscFree(facetIDsRemote);
327:     PetscLayoutDestroy(&fLayout);
328:     PetscSectionDestroy(&facetSectionRemote);
329:     PetscSectionDestroy(&facetSectionIDsRemote);
330:   }
331:   MEDfileClose(fileID);
332:   PetscFree(coordinates);
333:   PetscFree(cellList);
334:   PetscLayoutDestroy(&vLayout);
335:   PetscLayoutDestroy(&cLayout);
336:   PetscSFDestroy(&sfVertices);
337:   return 0;
338: #else
339:   SETERRQ(comm, PETSC_ERR_SUP, "This method requires Med mesh reader support. Reconfigure using --download-med");
340: #endif
341: }