Actual source code: plexmed.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  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], &ltogVertexNumbering);
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(&ltogVertexNumbering);
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: }