1 #include "../mmgcommon.h" /*I "petscdmplex.h" I*/
2 #include <parmmg/libparmmg.h>
3
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";
12
13 /*
14 Coupling code for the ParMmg metric-based mesh adaptation package.
15 */
DMAdaptMetric_ParMmg_Plex(DM dm,Vec vertexMetric,DMLabel bdLabel,DMLabel rgLabel,DM * dmNew)16 PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
17 {
18 MPI_Comm comm;
19 const char *bdName = "_boundary_";
20 const char *rgName = "_regions_";
21 DM udm, cdm;
22 DMLabel bdLabelNew, rgLabelNew;
23 const char *bdLabelName, *rgLabelName;
24 IS globalVertexNum;
25 PetscSection coordSection;
26 Vec coordinates;
27 PetscSF sf;
28 const PetscScalar *coords, *met;
29 PetscReal *vertices, *metric, *verticesNew, *verticesNewLoc, gradationFactor, hausdorffNumber;
30 PetscInt *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew;
31 PetscInt *bdFaces, *faceTags, *facesNew, *faceTagsNew;
32 PetscInt *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces;
33 PetscInt cStart, cEnd, c, numCells, fStart, fEnd, f, numFaceTags, vStart, vEnd, v, numVertices;
34 PetscInt numCellsNotShared, *cIsLeaf, numUsedVertices, *vertexNumber, *fIsIncluded;
35 PetscInt dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, numIter;
36 PetscInt *interfaces_lv, *interfaces_gv, *interfacesOffset;
37 PetscMPIInt niranks, nrranks, numNgbRanks;
38 PetscInt r, lv, gv;
39 PetscInt *gv_new, *owners, *verticesNewSorted, pStart, pEnd;
40 PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc;
41 const PetscInt *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote;
42 PetscBool flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform;
43 const PetscMPIInt *iranks, *rranks;
44 PetscMPIInt numProcs, rank;
45 PMMG_pParMesh parmesh = NULL;
46
47 // DEVELOPER NOTE: ParMmg wants to know the rank of every process which is sharing a given point and
48 // for this information to be conveyed to every process that is sharing that point.
49
50 PetscFunctionBegin;
51 PetscCall(PetscCitationsRegister(ParMmgCitation, &ParMmgCite));
52 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
53 PetscCallMPI(MPI_Comm_size(comm, &numProcs));
54 PetscCallMPI(MPI_Comm_rank(comm, &rank));
55 if (bdLabel) {
56 PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName));
57 PetscCall(PetscStrcmp(bdLabelName, bdName, &flg));
58 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
59 }
60 if (rgLabel) {
61 PetscCall(PetscObjectGetName((PetscObject)rgLabel, &rgLabelName));
62 PetscCall(PetscStrcmp(rgLabelName, rgName, &flg));
63 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for element tags", rgLabelName);
64 }
65
66 /* Get mesh information */
67 PetscCall(DMGetDimension(dm, &dim));
68 PetscCheck(dim == 3, comm, PETSC_ERR_ARG_OUTOFRANGE, "ParMmg only works in 3D.");
69 Neq = (dim * (dim + 1)) / 2;
70 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
71 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
72 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
73 PetscCall(DMPlexUninterpolate(dm, &udm));
74 PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL));
75 numCells = cEnd - cStart;
76 numVertices = vEnd - vStart;
77
78 /* Get parallel data; work out which cells are owned and which are leaves */
79 PetscCall(PetscCalloc1(numCells, &cIsLeaf));
80 numCellsNotShared = numCells;
81 niranks = nrranks = 0;
82 if (numProcs > 1) {
83 PetscCall(DMGetPointSF(dm, &sf));
84 PetscCall(PetscSFSetUp(sf));
85 PetscCall(PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc));
86 PetscCall(PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote));
87 for (r = 0; r < nrranks; ++r) {
88 for (i = roffset[r]; i < roffset[r + 1]; ++i) {
89 if (rmine[i] >= cStart && rmine[i] < cEnd) {
90 cIsLeaf[rmine[i] - cStart] = 1;
91 numCellsNotShared--;
92 }
93 }
94 }
95 }
96
97 /*
98 Create a vertex numbering for ParMmg starting at 1. Vertices not included in any
99 owned cell remain 0 and will be removed. Using this numbering, create cells.
100 */
101 numUsedVertices = 0;
102 PetscCall(PetscCalloc1(numVertices, &vertexNumber));
103 PetscCall(PetscMalloc1(numCellsNotShared * maxConeSize, &cells));
104 for (c = 0, coff = 0; c < numCells; ++c) {
105 const PetscInt *cone;
106 PetscInt coneSize, cl;
107
108 if (!cIsLeaf[c]) {
109 PetscCall(DMPlexGetConeSize(udm, cStart + c, &coneSize));
110 PetscCall(DMPlexGetCone(udm, cStart + c, &cone));
111 for (cl = 0; cl < coneSize; ++cl) {
112 if (!vertexNumber[cone[cl] - vStart]) vertexNumber[cone[cl] - vStart] = ++numUsedVertices;
113 cells[coff++] = vertexNumber[cone[cl] - vStart];
114 }
115 }
116 }
117
118 /* Get array of vertex coordinates */
119 PetscCall(DMGetCoordinateDM(dm, &cdm));
120 PetscCall(DMGetLocalSection(cdm, &coordSection));
121 PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
122 PetscCall(VecGetArrayRead(coordinates, &coords));
123 PetscCall(PetscMalloc2(numUsedVertices * Neq, &metric, dim * numUsedVertices, &vertices));
124 for (v = 0; v < vEnd - vStart; ++v) {
125 PetscCall(PetscSectionGetOffset(coordSection, v + vStart, &off));
126 if (vertexNumber[v]) {
127 for (i = 0; i < dim; ++i) vertices[dim * (vertexNumber[v] - 1) + i] = PetscRealPart(coords[off + i]);
128 }
129 }
130 PetscCall(VecRestoreArrayRead(coordinates, &coords));
131
132 /* Get face tags */
133 if (!bdLabel) {
134 flg = PETSC_TRUE;
135 PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel));
136 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabel));
137 }
138 PetscCall(DMLabelGetBounds(bdLabel, &pStart, &pEnd));
139 PetscCall(PetscCalloc1(pEnd - pStart, &fIsIncluded));
140 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
141 PetscBool hasPoint;
142 PetscInt *closure = NULL, closureSize, cl;
143
144 PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
145 if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
146
147 /* Only faces adjacent to an owned (non-leaf) cell are included */
148 PetscInt nnbrs;
149 const PetscInt *nbrs;
150 PetscCall(DMPlexGetSupportSize(dm, f, &nnbrs));
151 PetscCall(DMPlexGetSupport(dm, f, &nbrs));
152 for (c = 0; c < nnbrs; ++c) fIsIncluded[f - pStart] = fIsIncluded[f - pStart] || !cIsLeaf[nbrs[c]];
153 if (!fIsIncluded[f - pStart]) continue;
154
155 numFaceTags++;
156
157 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
158 for (cl = 0; cl < closureSize * 2; cl += 2) {
159 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
160 }
161 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
162 }
163 PetscCall(PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags));
164 for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
165 PetscInt *closure = NULL, closureSize, cl;
166
167 if (!fIsIncluded[f - pStart]) continue;
168
169 PetscCall(DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
170 for (cl = 0; cl < closureSize * 2; cl += 2) {
171 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = vertexNumber[closure[cl] - vStart];
172 }
173 PetscCall(DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure));
174 PetscCall(DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++]));
175 }
176 PetscCall(PetscFree(fIsIncluded));
177
178 /* Get cell tags */
179 PetscCall(PetscCalloc2(numUsedVertices, &verTags, numCellsNotShared, &cellTags));
180 if (rgLabel) {
181 for (c = cStart, coff = 0; c < cEnd; ++c) {
182 if (!cIsLeaf[c - cStart]) PetscCall(DMLabelGetValue(rgLabel, c, &cellTags[coff++]));
183 }
184 }
185 PetscCall(PetscFree(cIsLeaf));
186
187 /* Get metric, using only the upper triangular part */
188 PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view"));
189 PetscCall(VecGetArrayRead(vertexMetric, &met));
190 PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic));
191 PetscCall(DMPlexMetricIsUniform(dm, &uniform));
192 for (v = 0; v < (vEnd - vStart); ++v) {
193 PetscInt vv = vertexNumber[v];
194 if (!vv--) continue;
195 for (i = 0, k = 0; i < dim; ++i) {
196 for (j = i; j < dim; ++j, ++k) {
197 if (isotropic) {
198 if (i == j) {
199 if (uniform) metric[Neq * vv + k] = PetscRealPart(met[0]);
200 else metric[Neq * vv + k] = PetscRealPart(met[v]);
201 } else metric[Neq * vv + k] = 0.0;
202 } else metric[Neq * vv + k] = PetscRealPart(met[dim * dim * v + dim * i + j]);
203 }
204 }
205 }
206 PetscCall(VecRestoreArrayRead(vertexMetric, &met));
207
208 /* Build ParMmg communicators: the list of vertices between two partitions */
209 numNgbRanks = 0;
210 if (numProcs > 1) {
211 DM rankdm;
212 PetscSection rankSection, rankGlobalSection;
213 PetscSF rankSF;
214 const PetscInt *degree;
215 PetscInt *rankOfUsedVertices, *rankOfUsedMultiRootLeaves, *usedCopies;
216 PetscInt *rankArray, *rankGlobalArray, *interfacesPerRank;
217 PetscInt offset, mrl, rootDegreeCnt, s, shareCnt, gv;
218
219 PetscCall(PetscSFComputeDegreeBegin(sf, °ree));
220 PetscCall(PetscSFComputeDegreeEnd(sf, °ree));
221 PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
222 for (i = 0, rootDegreeCnt = 0; i < pEnd - pStart; ++i) rootDegreeCnt += degree[i];
223
224 /* rankOfUsedVertices, point-array: rank+1 if vertex and in use */
225 PetscCall(PetscCalloc1(pEnd - pStart, &rankOfUsedVertices));
226 for (i = 0; i < pEnd - pStart; ++i) rankOfUsedVertices[i] = -1;
227 for (i = vStart; i < vEnd; ++i) {
228 if (vertexNumber[i - vStart]) rankOfUsedVertices[i] = rank;
229 }
230
231 /* rankOfUsedMultiRootLeaves, multiroot-array: rank if vertex and in use, else -1 */
232 PetscCall(PetscMalloc1(rootDegreeCnt, &rankOfUsedMultiRootLeaves));
233 PetscCall(PetscSFGatherBegin(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves));
234 PetscCall(PetscSFGatherEnd(sf, MPIU_INT, rankOfUsedVertices, rankOfUsedMultiRootLeaves));
235 PetscCall(PetscFree(rankOfUsedVertices));
236
237 /* usedCopies, point-array: if vertex, shared by how many processes */
238 PetscCall(PetscCalloc1(pEnd - pStart, &usedCopies));
239 for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i];
240 for (i = vStart - pStart; i < vEnd - pStart; ++i) {
241 for (j = 0; j < degree[i]; j++, mrl++) {
242 if (rankOfUsedMultiRootLeaves[mrl] != -1) usedCopies[i]++;
243 }
244 if (vertexNumber[i - vStart + pStart]) usedCopies[i]++;
245 }
246 PetscCall(PetscSFBcastBegin(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE));
247 PetscCall(PetscSFBcastEnd(sf, MPIU_INT, usedCopies, usedCopies, MPI_REPLACE));
248
249 /* Create a section to store ranks of vertices shared by more than one process */
250 PetscCall(PetscSectionCreate(comm, &rankSection));
251 PetscCall(PetscSectionSetNumFields(rankSection, 1));
252 PetscCall(PetscSectionSetChart(rankSection, pStart, pEnd));
253 for (i = vStart - pStart; i < vEnd - pStart; ++i) {
254 if (usedCopies[i] > 1) PetscCall(PetscSectionSetDof(rankSection, i + pStart, usedCopies[i]));
255 }
256 PetscCall(PetscSectionSetUp(rankSection));
257 PetscCall(PetscSectionCreateGlobalSection(rankSection, sf, PETSC_FALSE, PETSC_FALSE, PETSC_TRUE, &rankGlobalSection));
258
259 PetscCall(PetscSectionGetStorageSize(rankGlobalSection, &s));
260 PetscCall(PetscMalloc1(s, &rankGlobalArray));
261 for (i = 0, mrl = 0; i < vStart - pStart; i++) mrl += degree[i];
262 for (i = vStart - pStart, k = 0; i < vEnd - pStart; ++i) {
263 if (usedCopies[i] > 1 && degree[i]) {
264 PetscCall(PetscSectionGetOffset(rankSection, k, &offset));
265 if (vertexNumber[i - vStart + pStart]) rankGlobalArray[k++] = rank;
266 for (j = 0; j < degree[i]; j++, mrl++) {
267 if (rankOfUsedMultiRootLeaves[mrl] != -1) rankGlobalArray[k++] = rankOfUsedMultiRootLeaves[mrl];
268 }
269 } else mrl += degree[i];
270 }
271 PetscCall(PetscFree(rankOfUsedMultiRootLeaves));
272 PetscCall(PetscFree(usedCopies));
273 PetscCall(PetscSectionDestroy(&rankGlobalSection));
274
275 /*
276 Broadcast the array of ranks.
277 (We want all processes to know all the ranks that are looking at each point.
278 Above, we tell the roots. Here, the roots tell the leaves.)
279 */
280 PetscCall(DMClone(dm, &rankdm));
281 PetscCall(DMSetLocalSection(rankdm, rankSection));
282 PetscCall(DMGetSectionSF(rankdm, &rankSF));
283 PetscCall(PetscSectionGetStorageSize(rankSection, &s));
284 PetscCall(PetscMalloc1(s, &rankArray));
285 PetscCall(PetscSFBcastBegin(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE));
286 PetscCall(PetscSFBcastEnd(rankSF, MPI_INT, rankGlobalArray, rankArray, MPI_REPLACE));
287 PetscCall(PetscFree(rankGlobalArray));
288 PetscCall(DMDestroy(&rankdm));
289
290 /* Count the number of interfaces per rank, not including those on the root */
291 PetscCall(PetscCalloc1(numProcs, &interfacesPerRank));
292 for (v = vStart; v < vEnd; v++) {
293 if (vertexNumber[v - vStart]) {
294 PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt));
295 if (shareCnt) {
296 PetscCall(PetscSectionGetOffset(rankSection, v, &offset));
297 for (j = 0; j < shareCnt; j++) interfacesPerRank[rankArray[offset + j]]++;
298 }
299 }
300 }
301 for (r = 0, k = 0, interfacesPerRank[rank] = 0; r < numProcs; r++) k += interfacesPerRank[r];
302
303 /* Get the degree of the vertex */
304 PetscCall(PetscMalloc3(k, &interfaces_lv, k, &interfaces_gv, numProcs + 1, &interfacesOffset));
305 interfacesOffset[0] = 0;
306 for (r = 0; r < numProcs; r++) {
307 interfacesOffset[r + 1] = interfacesOffset[r] + interfacesPerRank[r];
308 if (interfacesPerRank[r]) numNgbRanks++;
309 interfacesPerRank[r] = 0;
310 }
311
312 /* Get the local and global vertex numbers at interfaces */
313 PetscCall(DMPlexGetVertexNumbering(dm, &globalVertexNum));
314 PetscCall(ISGetIndices(globalVertexNum, &gV));
315 for (v = vStart; v < vEnd; v++) {
316 if (vertexNumber[v - vStart]) {
317 PetscCall(PetscSectionGetDof(rankSection, v, &shareCnt));
318 if (shareCnt) {
319 PetscCall(PetscSectionGetOffset(rankSection, v, &offset));
320 for (j = 0; j < shareCnt; j++) {
321 r = rankArray[offset + j];
322 if (r == rank) continue;
323 k = interfacesOffset[r] + interfacesPerRank[r]++;
324 interfaces_lv[k] = vertexNumber[v - vStart];
325 gv = gV[v - vStart];
326 interfaces_gv[k] = gv < 0 ? -gv : gv + 1;
327 }
328 }
329 }
330 }
331 PetscCall(PetscFree(interfacesPerRank));
332 PetscCall(PetscFree(rankArray));
333 PetscCall(ISRestoreIndices(globalVertexNum, &gV));
334 PetscCall(PetscSectionDestroy(&rankSection));
335 }
336 PetscCall(DMDestroy(&udm));
337 PetscCall(PetscFree(vertexNumber));
338
339 /* Send the data to ParMmg and remesh */
340 PetscCall(DMPlexMetricNoInsertion(dm, &noInsert));
341 PetscCall(DMPlexMetricNoSwapping(dm, &noSwap));
342 PetscCall(DMPlexMetricNoMovement(dm, &noMove));
343 PetscCall(DMPlexMetricNoSurf(dm, &noSurf));
344 PetscCall(DMPlexMetricGetVerbosity(dm, &verbosity));
345 PetscCall(DMPlexMetricGetNumIterations(dm, &numIter));
346 PetscCall(DMPlexMetricGetGradationFactor(dm, &gradationFactor));
347 PetscCall(DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber));
348 PetscCallMMG_NONSTANDARD(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);
349 PetscCallMMG_NONSTANDARD(PMMG_Set_meshSize, parmesh, numUsedVertices, numCellsNotShared, 0, numFaceTags, 0, 0);
350 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
351 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noinsert, noInsert);
352 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_noswap, noSwap);
353 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nomove, noMove);
354 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_nosurf, noSurf);
355 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_verbose, verbosity);
356 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
357 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_niter, numIter);
358 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hgrad, gradationFactor);
359 PetscCallMMG_NONSTANDARD(PMMG_Set_dparameter, parmesh, PMMG_DPARAM_hausd, hausdorffNumber);
360 PetscCallMMG_NONSTANDARD(PMMG_Set_vertices, parmesh, vertices, verTags);
361 PetscCallMMG_NONSTANDARD(PMMG_Set_tetrahedra, parmesh, cells, cellTags);
362 PetscCallMMG_NONSTANDARD(PMMG_Set_triangles, parmesh, bdFaces, faceTags);
363 PetscCallMMG_NONSTANDARD(PMMG_Set_metSize, parmesh, MMG5_Vertex, numUsedVertices, MMG5_Tensor);
364 PetscCallMMG_NONSTANDARD(PMMG_Set_tensorMets, parmesh, metric);
365 PetscCallMMG_NONSTANDARD(PMMG_Set_numberOfNodeCommunicators, parmesh, numNgbRanks);
366 for (r = 0, c = 0; r < numProcs; ++r) {
367 if (interfacesOffset[r + 1] > interfacesOffset[r]) {
368 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicatorSize, parmesh, c, r, interfacesOffset[r + 1] - interfacesOffset[r]);
369 PetscCallMMG_NONSTANDARD(PMMG_Set_ithNodeCommunicator_nodes, parmesh, c++, &interfaces_lv[interfacesOffset[r]], &interfaces_gv[interfacesOffset[r]], 1);
370 }
371 }
372 PetscCallMMG(PMMG_parmmglib_distributed, parmesh);
373 PetscCall(PetscFree(cells));
374 PetscCall(PetscFree2(metric, vertices));
375 PetscCall(PetscFree2(bdFaces, faceTags));
376 PetscCall(PetscFree2(verTags, cellTags));
377 if (numProcs > 1) PetscCall(PetscFree3(interfaces_lv, interfaces_gv, interfacesOffset));
378
379 /* Retrieve mesh from Mmg */
380 numCornersNew = 4;
381 PetscCallMMG_NONSTANDARD(PMMG_Get_meshSize, parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
382 PetscCall(PetscMalloc4(dim * numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer));
383 PetscCall(PetscMalloc3((dim + 1) * numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells));
384 PetscCall(PetscMalloc4(dim * numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces));
385 PetscCallMMG_NONSTANDARD(PMMG_Get_vertices, parmesh, verticesNew, verTagsNew, corners, requiredVer);
386 PetscCallMMG_NONSTANDARD(PMMG_Get_tetrahedra, parmesh, cellsNew, cellTagsNew, requiredCells);
387 PetscCallMMG_NONSTANDARD(PMMG_Get_triangles, parmesh, facesNew, faceTagsNew, requiredFaces);
388 PetscCall(PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new));
389 PetscCallMMG_NONSTANDARD(PMMG_Set_iparameter, parmesh, PMMG_IPARAM_globalNum, 1);
390 PetscCallMMG_NONSTANDARD(PMMG_Get_verticesGloNum, parmesh, gv_new, owners);
391 for (i = 0; i < dim * numFacesNew; ++i) facesNew[i] -= 1;
392 for (i = 0; i < (dim + 1) * numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i] - 1] - 1;
393 for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) {
394 if (owners[i] == rank) numVerticesNewLoc++;
395 }
396 PetscCall(PetscMalloc1(numVerticesNewLoc * dim, &verticesNewLoc));
397 for (i = 0, c = 0; i < numVerticesNew; i++) {
398 if (owners[i] == rank) {
399 for (j = 0; j < dim; ++j) verticesNewLoc[dim * c + j] = verticesNew[dim * i + j];
400 c++;
401 }
402 }
403
404 /* Reorder for consistency with DMPlex */
405 for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i]));
406
407 /* Create new plex */
408 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew));
409 PetscCallMMG_NONSTANDARD(PMMG_Free_all, PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end);
410 PetscCall(PetscFree4(verticesNew, verTagsNew, corners, requiredVer));
411
412 /* Get adapted mesh information */
413 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd));
414 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd));
415 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd));
416
417 /* Rebuild boundary label */
418 PetscCall(DMCreateLabel(*dmNew, flg ? bdName : bdLabelName));
419 PetscCall(DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew));
420 for (i = 0; i < numFacesNew; i++) {
421 PetscBool hasTag = PETSC_FALSE;
422 PetscInt numCoveredPoints, numFaces = 0, facePoints[3];
423 const PetscInt *coveredPoints = NULL;
424
425 for (j = 0; j < dim; ++j) {
426 lv = facesNew[i * dim + j];
427 gv = gv_new[lv] - 1;
428 PetscCall(PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv));
429 facePoints[j] = lv + vStart;
430 }
431 PetscCall(DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
432 for (j = 0; j < numCoveredPoints; ++j) {
433 if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
434 numFaces++;
435 f = j;
436 }
437 }
438 PetscCheck(numFaces == 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "%" PetscInt_FMT " vertices cannot define more than 1 facet (%" PetscInt_FMT ")", dim, numFaces);
439 PetscCall(DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag));
440 if (hasTag) PetscCall(DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]));
441 PetscCall(DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints));
442 }
443 PetscCall(PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces));
444 PetscCall(PetscFree2(owners, gv_new));
445 PetscCall(PetscFree(verticesNewLoc));
446 if (flg) PetscCall(DMLabelDestroy(&bdLabel));
447
448 /* Rebuild cell labels */
449 PetscCall(DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName));
450 PetscCall(DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew));
451 for (c = cStart; c < cEnd; ++c) PetscCall(DMLabelSetValue(rgLabelNew, c, cellTagsNew[c - cStart]));
452 PetscCall(PetscFree3(cellsNew, cellTagsNew, requiredCells));
453 PetscFunctionReturn(PETSC_SUCCESS);
454 }
455