Actual source code: parmmgadapt.c

  1: #include "../mmgcommon.h" /*I      "petscdmplex.h"   I*/
  2: #include <parmmg/libparmmg.h>

  4: PetscBool  ParMmgCite       = PETSC_FALSE;
  5: const char ParMmgCitation[] = "@techreport{cirrottola:hal-02386837,\n"
  6:                               "  title       = {Parallel unstructured mesh adaptation using iterative remeshing and repartitioning},\n"
  7:                               "  institution = {Inria Bordeaux},\n"
  8:                               "  author      = {L. Cirrottola and A. Froehly},\n"
  9:                               "  number      = {9307},\n"
 10:                               "  note        = {\\url{https://hal.inria.fr/hal-02386837}},\n"
 11:                               "  year        = {2019}\n}\n";

 13: PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
 14: {
 15:   MPI_Comm           comm;
 16:   const char        *bdName = "_boundary_";
 17:   const char        *rgName = "_regions_";
 18:   DM                 udm, cdm;
 19:   DMLabel            bdLabelNew, rgLabelNew;
 20:   const char        *bdLabelName, *rgLabelName;
 21:   IS                 globalVertexNum;
 22:   PetscSection       coordSection;
 23:   Vec                coordinates;
 24:   PetscSF            sf;
 25:   const PetscScalar *coords, *met;
 26:   PetscReal         *vertices, *metric, *verticesNew, *verticesNewLoc, gradationFactor, hausdorffNumber;
 27:   PetscInt          *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew;
 28:   PetscInt          *bdFaces, *faceTags, *facesNew, *faceTagsNew;
 29:   PetscInt          *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces;
 30:   PetscInt           cStart, cEnd, c, numCells, fStart, fEnd, f, numFaceTags, vStart, vEnd, v, numVertices;
 31:   PetscInt           dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, numIter;
 32:   PetscInt          *numVerInterfaces, *ngbRanks, *verNgbRank, *interfaces_lv, *interfaces_gv, *intOffset;
 33:   PetscInt           niranks, nrranks, numNgbRanks, numVerNgbRanksTotal, count, sliceSize, p, r, n, lv, gv;
 34:   PetscInt          *gv_new, *owners, *verticesNewSorted, pStart, pEnd;
 35:   PetscInt           numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc;
 36:   const PetscInt    *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote;
 37:   PetscBool          flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform;
 38:   const PetscMPIInt *iranks, *rranks;
 39:   PetscMPIInt        numProcs, rank;
 40:   PMMG_pParMesh      parmesh = NULL;

 42:   PetscCitationsRegister(ParMmgCitation, &ParMmgCite);
 43:   PetscObjectGetComm((PetscObject)dm, &comm);
 44:   MPI_Comm_size(comm, &numProcs);
 45:   MPI_Comm_rank(comm, &rank);
 46:   if (bdLabel) {
 47:     PetscObjectGetName((PetscObject)bdLabel, &bdLabelName);
 48:     PetscStrcmp(bdLabelName, bdName, &flg);
 50:   }
 51:   if (rgLabel) {
 52:     PetscObjectGetName((PetscObject)rgLabel, &rgLabelName);
 53:     PetscStrcmp(rgLabelName, rgName, &flg);
 55:   }

 57:   /* Get mesh information */
 58:   DMGetDimension(dm, &dim);
 60:   Neq = (dim * (dim + 1)) / 2;
 61:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
 62:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
 63:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
 64:   DMPlexUninterpolate(dm, &udm);
 65:   DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
 66:   numCells    = cEnd - cStart;
 67:   numVertices = vEnd - vStart;

 69:   /* Get cell offsets */
 70:   PetscMalloc1(numCells * maxConeSize, &cells);
 71:   for (c = 0, coff = 0; c < numCells; ++c) {
 72:     const PetscInt *cone;
 73:     PetscInt        coneSize, cl;

 75:     DMPlexGetConeSize(udm, c, &coneSize);
 76:     DMPlexGetCone(udm, c, &cone);
 77:     for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart + 1;
 78:   }

 80:   /* Get vertex coordinate array */
 81:   DMGetCoordinateDM(dm, &cdm);
 82:   DMGetLocalSection(cdm, &coordSection);
 83:   DMGetCoordinatesLocal(dm, &coordinates);
 84:   VecGetArrayRead(coordinates, &coords);
 85:   PetscMalloc2(numVertices * Neq, &metric, dim * numVertices, &vertices);
 86:   for (v = 0; v < vEnd - vStart; ++v) {
 87:     PetscSectionGetOffset(coordSection, v + vStart, &off);
 88:     for (i = 0; i < dim; ++i) vertices[dim * v + i] = PetscRealPart(coords[off + i]);
 89:   }
 90:   VecRestoreArrayRead(coordinates, &coords);

 92:   /* Get face tags */
 93:   if (!bdLabel) {
 94:     flg = PETSC_TRUE;
 95:     DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel);
 96:     DMPlexMarkBoundaryFaces(dm, 1, bdLabel);
 97:   }
 98:   DMLabelGetBounds(bdLabel, &pStart, &pEnd);
 99:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
100:     PetscBool hasPoint;
101:     PetscInt *closure = NULL, closureSize, cl;

103:     DMLabelHasPoint(bdLabel, f, &hasPoint);
104:     if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
105:     numFaceTags++;

107:     DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
108:     for (cl = 0; cl < closureSize * 2; cl += 2) {
109:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
110:     }
111:     DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
112:   }
113:   PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags);
114:   for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
115:     PetscBool hasPoint;
116:     PetscInt *closure = NULL, closureSize, cl;

118:     DMLabelHasPoint(bdLabel, f, &hasPoint);
119:     if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;

121:     DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
122:     for (cl = 0; cl < closureSize * 2; cl += 2) {
123:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart + 1;
124:     }
125:     DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
126:     DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++]);
127:   }

129:   /* Get cell tags */
130:   PetscCalloc2(numVertices, &verTags, numCells, &cellTags);
131:   if (rgLabel) {
132:     for (c = cStart; c < cEnd; ++c) DMLabelGetValue(rgLabel, c, &cellTags[c]);
133:   }

135:   /* Get metric */
136:   VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
137:   VecGetArrayRead(vertexMetric, &met);
138:   DMPlexMetricIsIsotropic(dm, &isotropic);
139:   DMPlexMetricIsUniform(dm, &uniform);
140:   for (v = 0; v < (vEnd - vStart); ++v) {
141:     for (i = 0, k = 0; i < dim; ++i) {
142:       for (j = i; j < dim; ++j, ++k) {
143:         if (isotropic) {
144:           if (i == j) {
145:             if (uniform) metric[Neq * v + k] = PetscRealPart(met[0]);
146:             else metric[Neq * v + k] = PetscRealPart(met[v]);
147:           } else metric[Neq * v + k] = 0.0;
148:         } else metric[Neq * v + k] = PetscRealPart(met[dim * dim * v + dim * i + j]);
149:       }
150:     }
151:   }
152:   VecRestoreArrayRead(vertexMetric, &met);

154:   /* Build ParMMG communicators: the list of vertices between two partitions  */
155:   niranks = nrranks = 0;
156:   numNgbRanks       = 0;
157:   if (numProcs > 1) {
158:     DMGetPointSF(dm, &sf);
159:     PetscSFSetUp(sf);
160:     PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc);
161:     PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote);
162:     PetscCalloc1(numProcs, &numVerInterfaces);

164:     /* Count number of roots associated with each leaf */
165:     for (r = 0; r < niranks; ++r) {
166:       for (i = ioffset[r], count = 0; i < ioffset[r + 1]; ++i) {
167:         if (irootloc[i] >= vStart && irootloc[i] < vEnd) count++;
168:       }
169:       numVerInterfaces[iranks[r]] += count;
170:     }

172:     /* Count number of leaves associated with each root */
173:     for (r = 0; r < nrranks; ++r) {
174:       for (i = roffset[r], count = 0; i < roffset[r + 1]; ++i) {
175:         if (rmine[i] >= vStart && rmine[i] < vEnd) count++;
176:       }
177:       numVerInterfaces[rranks[r]] += count;
178:     }

180:     /* Count global number of ranks */
181:     for (p = 0; p < numProcs; ++p) {
182:       if (numVerInterfaces[p]) numNgbRanks++;
183:     }

185:     /* Provide numbers of vertex interfaces */
186:     PetscMalloc2(numNgbRanks, &ngbRanks, numNgbRanks, &verNgbRank);
187:     for (p = 0, n = 0; p < numProcs; ++p) {
188:       if (numVerInterfaces[p]) {
189:         ngbRanks[n]   = p;
190:         verNgbRank[n] = numVerInterfaces[p];
191:         n++;
192:       }
193:     }
194:     numVerNgbRanksTotal = 0;
195:     for (p = 0; p < numNgbRanks; ++p) numVerNgbRanksTotal += verNgbRank[p];

197:     /* For each neighbor, fill in interface arrays */
198:     PetscMalloc3(numVerNgbRanksTotal, &interfaces_lv, numVerNgbRanksTotal, &interfaces_gv, numNgbRanks + 1, &intOffset);
199:     intOffset[0] = 0;
200:     for (p = 0, r = 0, i = 0; p < numNgbRanks; ++p) {
201:       intOffset[p + 1] = intOffset[p];

203:       /* Leaf case */
204:       if (iranks && iranks[i] == ngbRanks[p]) {
205:         /* Add the right slice of irootloc at the right place */
206:         sliceSize = ioffset[i + 1] - ioffset[i];
207:         for (j = 0, count = 0; j < sliceSize; ++j) {
209:           v = irootloc[ioffset[i] + j];
210:           if (v >= vStart && v < vEnd) {
212:             interfaces_lv[intOffset[p + 1] + count] = v - vStart;
213:             count++;
214:           }
215:         }
216:         intOffset[p + 1] += count;
217:         i++;
218:       }

220:       /* Root case */
221:       if (rranks && rranks[r] == ngbRanks[p]) {
222:         /* Add the right slice of rmine at the right place */
223:         sliceSize = roffset[r + 1] - roffset[r];
224:         for (j = 0, count = 0; j < sliceSize; ++j) {
226:           v = rmine[roffset[r] + j];
227:           if (v >= vStart && v < vEnd) {
229:             interfaces_lv[intOffset[p + 1] + count] = v - vStart;
230:             count++;
231:           }
232:         }
233:         intOffset[p + 1] += count;
234:         r++;
235:       }

237:       /* Check validity of offsets */
239:     }
240:     DMPlexGetVertexNumbering(udm, &globalVertexNum);
241:     ISGetIndices(globalVertexNum, &gV);
242:     for (i = 0; i < numVerNgbRanksTotal; ++i) {
243:       v                = gV[interfaces_lv[i]];
244:       interfaces_gv[i] = v < 0 ? -v - 1 : v;
245:       interfaces_lv[i] += 1;
246:       interfaces_gv[i] += 1;
247:     }
248:     ISRestoreIndices(globalVertexNum, &gV);
249:     PetscFree(numVerInterfaces);
250:   }
251:   DMDestroy(&udm);

253:   /* Send the data to ParMmg and remesh */
254:   DMPlexMetricNoInsertion(dm, &noInsert);
255:   DMPlexMetricNoSwapping(dm, &noSwap);
256:   DMPlexMetricNoMovement(dm, &noMove);
257:   DMPlexMetricNoSurf(dm, &noSurf);
258:   DMPlexMetricGetVerbosity(dm, &verbosity);
259:   DMPlexMetricGetNumIterations(dm, &numIter);
260:   DMPlexMetricGetGradationFactor(dm, &gradationFactor);
261:   DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber);
262:   PMMG_Init_parMesh(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_pMesh, PMMG_ARG_pMet, PMMG_ARG_dim, 3, PMMG_ARG_MPIComm, comm, PMMG_ARG_end);
263:   PMMG_Set_meshSize(parmesh, numVertices, numCells, 0, numFaceTags, 0, 0);
264:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
265:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noinsert, noInsert);
266:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noswap, noSwap);
267:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nomove, noMove);
268:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nosurf, noSurf);
269:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_verbose, verbosity);
270:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
271:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_niter, numIter);
272:   PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hgrad, gradationFactor);
273:   PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hausd, hausdorffNumber);
274:   PMMG_Set_vertices(parmesh, vertices, verTags);
275:   PMMG_Set_tetrahedra(parmesh, cells, cellTags);
276:   PMMG_Set_triangles(parmesh, bdFaces, faceTags);
277:   PMMG_Set_metSize(parmesh, MMG5_Vertex, numVertices, MMG5_Tensor);
278:   PMMG_Set_tensorMets(parmesh, metric);
279:   PMMG_Set_numberOfNodeCommunicators(parmesh, numNgbRanks);
280:   for (c = 0; c < numNgbRanks; ++c) {
281:     PMMG_Set_ithNodeCommunicatorSize(parmesh, c, ngbRanks[c], intOffset[c + 1] - intOffset[c]);
282:     PMMG_Set_ithNodeCommunicator_nodes(parmesh, c, &interfaces_lv[intOffset[c]], &interfaces_gv[intOffset[c]], 1);
283:   }
284:   PMMG_parmmglib_distributed(parmesh);
285:   PetscFree(cells);
286:   PetscFree2(metric, vertices);
287:   PetscFree2(bdFaces, faceTags);
288:   PetscFree2(verTags, cellTags);
289:   if (numProcs > 1) {
290:     PetscFree2(ngbRanks, verNgbRank);
291:     PetscFree3(interfaces_lv, interfaces_gv, intOffset);
292:   }

294:   /* Retrieve mesh from Mmg */
295:   numCornersNew = 4;
296:   PMMG_Get_meshSize(parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
297:   PetscMalloc4(dim * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);
298:   PetscMalloc3((dim + 1) * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);
299:   PetscMalloc4(dim * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);
300:   PMMG_Get_vertices(parmesh, verticesNew, verTagsNew, corners, requiredVer);
301:   PMMG_Get_tetrahedra(parmesh, cellsNew, cellTagsNew, requiredCells);
302:   PMMG_Get_triangles(parmesh, facesNew, faceTagsNew, requiredFaces);
303:   PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new);
304:   PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
305:   PMMG_Get_verticesGloNum(parmesh, gv_new, owners);
306:   for (i = 0; i < dim * numFacesNew; ++i) facesNew[i] -= 1;
307:   for (i = 0; i < (dim + 1) * numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i] - 1] - 1;
308:   for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) {
309:     if (owners[i] == rank) numVerticesNewLoc++;
310:   }
311:   PetscMalloc2(numVerticesNewLoc * dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted);
312:   for (i = 0, c = 0; i < numVerticesNew; i++) {
313:     if (owners[i] == rank) {
314:       for (j = 0; j < dim; ++j) verticesNewLoc[dim * c + j] = verticesNew[dim * i + j];
315:       c++;
316:     }
317:   }

319:   /* Reorder for consistency with DMPlex */
320:   for (i = 0; i < numCellsNew; ++i) DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i]);

322:   /* Create new plex */
323:   DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew);
324:   PMMG_Free_all(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end);
325:   PetscFree4(verticesNew, verTagsNew, corners, requiredVer);

327:   /* Get adapted mesh information */
328:   DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
329:   DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
330:   DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);

332:   /* Rebuild boundary label */
333:   DMCreateLabel(*dmNew, flg ? bdName : bdLabelName);
334:   DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew);
335:   for (i = 0; i < numFacesNew; i++) {
336:     PetscBool       hasTag = PETSC_FALSE;
337:     PetscInt        numCoveredPoints, numFaces = 0, facePoints[3];
338:     const PetscInt *coveredPoints = NULL;

340:     for (j = 0; j < dim; ++j) {
341:       lv = facesNew[i * dim + j];
342:       gv = gv_new[lv] - 1;
343:       PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv);
344:       facePoints[j] = lv + vStart;
345:     }
346:     DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);
347:     for (j = 0; j < numCoveredPoints; ++j) {
348:       if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
349:         numFaces++;
350:         f = j;
351:       }
352:     }
354:     DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag);
355:     if (hasTag) DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]);
356:     DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);
357:   }
358:   PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces);
359:   PetscFree2(owners, gv_new);
360:   PetscFree2(verticesNewLoc, verticesNewSorted);
361:   if (flg) DMLabelDestroy(&bdLabel);

363:   /* Rebuild cell labels */
364:   DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName);
365:   DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew);
366:   for (c = cStart; c < cEnd; ++c) DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart]);
367:   PetscFree3(cellsNew, cellTagsNew, requiredCells);

369:   return 0;
370: }