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: /*
14: Coupling code for the ParMmg metric-based mesh adaptation package.
15: */
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;
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.
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: }
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;
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: }
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;
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: }
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));
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;
144: PetscCall(DMLabelHasPoint(bdLabel, f, &hasPoint));
145: if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
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;
155: numFaceTags++;
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;
167: if (!fIsIncluded[f - pStart]) continue;
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));
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));
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));
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;
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];
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: }
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));
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));
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));
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));
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));
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];
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: }
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));
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)); }
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(PetscMalloc2(numVerticesNewLoc * dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted));
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: }
404: /* Reorder for consistency with DMPlex */
405: for (i = 0; i < numCellsNew; ++i) PetscCall(DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4 * i]));
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));
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));
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;
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(PetscFree2(verticesNewLoc, verticesNewSorted));
446: if (flg) PetscCall(DMLabelDestroy(&bdLabel));
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: }