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: }