Actual source code: parmmgadapt.c
1: #include "../mmgcommon.h" /*I "petscdmplex.h" I*/
2: #include <parmmg/libparmmg.h>
4: PETSC_EXTERN PetscErrorCode DMAdaptMetric_ParMmg_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
5: {
6: MPI_Comm comm;
7: const char *bdName = "_boundary_";
8: const char *rgName = "_regions_";
9: DM udm, cdm;
10: DMLabel bdLabelNew, rgLabelNew;
11: const char *bdLabelName, *rgLabelName;
12: IS globalVertexNum;
13: PetscSection coordSection;
14: Vec coordinates;
15: PetscSF sf;
16: const PetscScalar *coords, *met;
17: PetscReal *vertices, *metric, *verticesNew, *verticesNewLoc, gradationFactor, hausdorffNumber;
18: PetscInt *cells, *cellsNew, *cellTags, *cellTagsNew, *verTags, *verTagsNew;
19: PetscInt *bdFaces, *faceTags, *facesNew, *faceTagsNew;
20: PetscInt *corners, *requiredCells, *requiredVer, *ridges, *requiredFaces;
21: PetscInt cStart, cEnd, c, numCells, fStart, fEnd, f, numFaceTags, vStart, vEnd, v, numVertices;
22: PetscInt dim, off, coff, maxConeSize, bdSize, i, j, k, Neq, verbosity, numIter;
23: PetscInt *numVerInterfaces, *ngbRanks, *verNgbRank, *interfaces_lv, *interfaces_gv, *intOffset;
24: PetscInt niranks, nrranks, numNgbRanks, numVerNgbRanksTotal, count, sliceSize, p, r, n, lv, gv;
25: PetscInt *gv_new, *owners, *verticesNewSorted, pStart, pEnd;
26: PetscInt numCellsNew, numVerticesNew, numCornersNew, numFacesNew, numVerticesNewLoc;
27: const PetscInt *gV, *ioffset, *irootloc, *roffset, *rmine, *rremote;
28: PetscBool flg = PETSC_FALSE, noInsert, noSwap, noMove, noSurf, isotropic, uniform;
29: const PetscMPIInt *iranks, *rranks;
30: PetscMPIInt numProcs, rank;
31: PMMG_pParMesh parmesh = NULL;
33: PetscObjectGetComm((PetscObject) dm, &comm);
34: MPI_Comm_size(comm, &numProcs);
35: MPI_Comm_rank(comm, &rank);
36: if (bdLabel) {
37: PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);
38: PetscStrcmp(bdLabelName, bdName, &flg);
40: }
41: if (rgLabel) {
42: PetscObjectGetName((PetscObject) rgLabel, &rgLabelName);
43: PetscStrcmp(rgLabelName, rgName, &flg);
45: }
47: /* Get mesh information */
48: DMGetDimension(dm, &dim);
50: Neq = (dim*(dim+1))/2;
51: DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
52: DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
53: DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
54: DMPlexUninterpolate(dm, &udm);
55: DMPlexGetMaxSizes(udm, &maxConeSize, NULL);
56: numCells = cEnd - cStart;
57: numVertices = vEnd - vStart;
59: /* Get cell offsets */
60: PetscMalloc1(numCells*maxConeSize, &cells);
61: for (c = 0, coff = 0; c < numCells; ++c) {
62: const PetscInt *cone;
63: PetscInt coneSize, cl;
65: DMPlexGetConeSize(udm, c, &coneSize);
66: DMPlexGetCone(udm, c, &cone);
67: for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart+1;
68: }
70: /* Get vertex coordinate array */
71: DMGetCoordinateDM(dm, &cdm);
72: DMGetLocalSection(cdm, &coordSection);
73: DMGetCoordinatesLocal(dm, &coordinates);
74: VecGetArrayRead(coordinates, &coords);
75: PetscMalloc2(numVertices*Neq, &metric, dim*numVertices, &vertices);
76: for (v = 0; v < vEnd-vStart; ++v) {
77: PetscSectionGetOffset(coordSection, v+vStart, &off);
78: for (i = 0; i < dim; ++i) vertices[dim*v+i] = PetscRealPart(coords[off+i]);
79: }
80: VecRestoreArrayRead(coordinates, &coords);
82: /* Get face tags */
83: if (!bdLabel) {
84: flg = PETSC_TRUE;
85: DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabel);
86: DMPlexMarkBoundaryFaces(dm, 1, bdLabel);
87: }
88: DMLabelGetBounds(bdLabel, &pStart, &pEnd);
89: for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
90: PetscBool hasPoint;
91: PetscInt *closure = NULL, closureSize, cl;
93: DMLabelHasPoint(bdLabel, f, &hasPoint);
94: if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
95: numFaceTags++;
97: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
98: for (cl = 0; cl < closureSize*2; cl += 2) {
99: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
100: }
101: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
102: }
103: PetscMalloc2(bdSize, &bdFaces, numFaceTags, &faceTags);
104: for (f = pStart, bdSize = 0, numFaceTags = 0; f < pEnd; ++f) {
105: PetscBool hasPoint;
106: PetscInt *closure = NULL, closureSize, cl;
108: DMLabelHasPoint(bdLabel, f, &hasPoint);
109: if ((!hasPoint) || (f < fStart) || (f >= fEnd)) continue;
111: DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
112: for (cl = 0; cl < closureSize*2; cl += 2) {
113: if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart+1;
114: }
115: DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &closureSize, &closure);
116: DMLabelGetValue(bdLabel, f, &faceTags[numFaceTags++]);
117: }
119: /* Get cell tags */
120: PetscCalloc2(numVertices, &verTags, numCells, &cellTags);
121: if (rgLabel) {
122: for (c = cStart; c < cEnd; ++c) DMLabelGetValue(rgLabel, c, &cellTags[c]);
123: }
125: /* Get metric */
126: VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");
127: VecGetArrayRead(vertexMetric, &met);
128: DMPlexMetricIsIsotropic(dm, &isotropic);
129: DMPlexMetricIsUniform(dm, &uniform);
130: for (v = 0; v < (vEnd-vStart); ++v) {
131: for (i = 0, k = 0; i < dim; ++i) {
132: for (j = i; j < dim; ++j, ++k) {
133: if (isotropic) {
134: if (i == j) {
135: if (uniform) metric[Neq*v+k] = PetscRealPart(met[0]);
136: else metric[Neq*v+k] = PetscRealPart(met[v]);
137: } else metric[Neq*v+k] = 0.0;
138: } else metric[Neq*v+k] = PetscRealPart(met[dim*dim*v+dim*i+j]);
139: }
140: }
141: }
142: VecRestoreArrayRead(vertexMetric, &met);
144: /* Build ParMMG communicators: the list of vertices between two partitions */
145: niranks = nrranks = 0;
146: numNgbRanks = 0;
147: if (numProcs > 1) {
148: DMGetPointSF(dm, &sf);
149: PetscSFSetUp(sf);
150: PetscSFGetLeafRanks(sf, &niranks, &iranks, &ioffset, &irootloc);
151: PetscSFGetRootRanks(sf, &nrranks, &rranks, &roffset, &rmine, &rremote);
152: PetscCalloc1(numProcs, &numVerInterfaces);
154: /* Count number of roots associated with each leaf */
155: for (r = 0; r < niranks; ++r) {
156: for (i=ioffset[r], count=0; i<ioffset[r+1]; ++i) {
157: if (irootloc[i] >= vStart && irootloc[i] < vEnd) count++;
158: }
159: numVerInterfaces[iranks[r]] += count;
160: }
162: /* Count number of leaves associated with each root */
163: for (r = 0; r < nrranks; ++r) {
164: for (i=roffset[r], count=0; i<roffset[r+1]; ++i) {
165: if (rmine[i] >= vStart && rmine[i] < vEnd) count++;
166: }
167: numVerInterfaces[rranks[r]] += count;
168: }
170: /* Count global number of ranks */
171: for (p = 0; p < numProcs; ++p) {
172: if (numVerInterfaces[p]) numNgbRanks++;
173: }
175: /* Provide numbers of vertex interfaces */
176: PetscMalloc2(numNgbRanks, &ngbRanks, numNgbRanks, &verNgbRank);
177: for (p = 0, n = 0; p < numProcs; ++p) {
178: if (numVerInterfaces[p]) {
179: ngbRanks[n] = p;
180: verNgbRank[n] = numVerInterfaces[p];
181: n++;
182: }
183: }
184: numVerNgbRanksTotal = 0;
185: for (p = 0; p < numNgbRanks; ++p) numVerNgbRanksTotal += verNgbRank[p];
187: /* For each neighbor, fill in interface arrays */
188: PetscMalloc3(numVerNgbRanksTotal, &interfaces_lv, numVerNgbRanksTotal, &interfaces_gv, numNgbRanks+1, &intOffset);
189: intOffset[0] = 0;
190: for (p = 0, r = 0, i = 0; p < numNgbRanks; ++p) {
191: intOffset[p+1] = intOffset[p];
193: /* Leaf case */
194: if (iranks && iranks[i] == ngbRanks[p]) {
196: /* Add the right slice of irootloc at the right place */
197: sliceSize = ioffset[i+1]-ioffset[i];
198: for (j = 0, count = 0; j < sliceSize; ++j) {
200: v = irootloc[ioffset[i]+j];
201: if (v >= vStart && v < vEnd) {
203: interfaces_lv[intOffset[p+1]+count] = v-vStart;
204: count++;
205: }
206: }
207: intOffset[p+1] += count;
208: i++;
209: }
211: /* Root case */
212: if (rranks && rranks[r] == ngbRanks[p]) {
214: /* Add the right slice of rmine at the right place */
215: sliceSize = roffset[r+1]-roffset[r];
216: for (j = 0, count = 0; j < sliceSize; ++j) {
218: v = rmine[roffset[r]+j];
219: if (v >= vStart && v < vEnd) {
221: interfaces_lv[intOffset[p+1]+count] = v-vStart;
222: count++;
223: }
224: }
225: intOffset[p+1] += count;
226: r++;
227: }
229: /* Check validity of offsets */
231: }
232: DMPlexGetVertexNumbering(udm, &globalVertexNum);
233: ISGetIndices(globalVertexNum, &gV);
234: for (i = 0; i < numVerNgbRanksTotal; ++i) {
235: v = gV[interfaces_lv[i]];
236: interfaces_gv[i] = v < 0 ? -v-1 : v;
237: interfaces_lv[i] += 1;
238: interfaces_gv[i] += 1;
239: }
240: ISRestoreIndices(globalVertexNum, &gV);
241: PetscFree(numVerInterfaces);
242: }
243: DMDestroy(&udm);
245: /* Send the data to ParMmg and remesh */
246: DMPlexMetricNoInsertion(dm, &noInsert);
247: DMPlexMetricNoSwapping(dm, &noSwap);
248: DMPlexMetricNoMovement(dm, &noMove);
249: DMPlexMetricNoSurf(dm, &noSurf);
250: DMPlexMetricGetVerbosity(dm, &verbosity);
251: DMPlexMetricGetNumIterations(dm, &numIter);
252: DMPlexMetricGetGradationFactor(dm, &gradationFactor);
253: DMPlexMetricGetHausdorffNumber(dm, &hausdorffNumber);
254: 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);
255: PMMG_Set_meshSize(parmesh, numVertices, numCells, 0, numFaceTags, 0, 0);
256: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_APImode, PMMG_APIDISTRIB_nodes);
257: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noinsert, noInsert);
258: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_noswap, noSwap);
259: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nomove, noMove);
260: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_nosurf, noSurf);
261: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_verbose, verbosity);
262: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
263: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_niter, numIter);
264: PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hgrad, gradationFactor);
265: PMMG_Set_dparameter(parmesh, PMMG_DPARAM_hausd, hausdorffNumber);
266: PMMG_Set_vertices(parmesh, vertices, verTags);
267: PMMG_Set_tetrahedra(parmesh, cells, cellTags);
268: PMMG_Set_triangles(parmesh, bdFaces, faceTags);
269: PMMG_Set_metSize(parmesh, MMG5_Vertex, numVertices, MMG5_Tensor);
270: PMMG_Set_tensorMets(parmesh, metric);
271: PMMG_Set_numberOfNodeCommunicators(parmesh, numNgbRanks);
272: for (c = 0; c < numNgbRanks; ++c) {
273: PMMG_Set_ithNodeCommunicatorSize(parmesh, c, ngbRanks[c], intOffset[c+1]-intOffset[c]);
274: PMMG_Set_ithNodeCommunicator_nodes(parmesh, c, &interfaces_lv[intOffset[c]], &interfaces_gv[intOffset[c]], 1);
275: }
276: PMMG_parmmglib_distributed(parmesh);
277: PetscFree(cells);
278: PetscFree2(metric, vertices);
279: PetscFree2(bdFaces, faceTags);
280: PetscFree2(verTags, cellTags);
281: if (numProcs > 1) {
282: PetscFree2(ngbRanks, verNgbRank);
283: PetscFree3(interfaces_lv, interfaces_gv, intOffset);
284: }
286: /* Retrieve mesh from Mmg */
287: numCornersNew = 4;
288: PMMG_Get_meshSize(parmesh, &numVerticesNew, &numCellsNew, 0, &numFacesNew, 0, 0);
289: PetscMalloc4(dim*numVerticesNew, &verticesNew, numVerticesNew, &verTagsNew, numVerticesNew, &corners, numVerticesNew, &requiredVer);
290: PetscMalloc3((dim+1)*numCellsNew, &cellsNew, numCellsNew, &cellTagsNew, numCellsNew, &requiredCells);
291: PetscMalloc4(dim*numFacesNew, &facesNew, numFacesNew, &faceTagsNew, numFacesNew, &ridges, numFacesNew, &requiredFaces);
292: PMMG_Get_vertices(parmesh, verticesNew, verTagsNew, corners, requiredVer);
293: PMMG_Get_tetrahedra(parmesh, cellsNew, cellTagsNew, requiredCells);
294: PMMG_Get_triangles(parmesh, facesNew, faceTagsNew, requiredFaces);
295: PetscMalloc2(numVerticesNew, &owners, numVerticesNew, &gv_new);
296: PMMG_Set_iparameter(parmesh, PMMG_IPARAM_globalNum, 1);
297: PMMG_Get_verticesGloNum(parmesh, gv_new, owners);
298: for (i = 0; i < dim*numFacesNew; ++i) facesNew[i] -= 1;
299: for (i = 0; i < (dim+1)*numCellsNew; ++i) cellsNew[i] = gv_new[cellsNew[i]-1]-1;
300: for (i = 0, numVerticesNewLoc = 0; i < numVerticesNew; ++i) {
301: if (owners[i] == rank) numVerticesNewLoc++;
302: }
303: PetscMalloc2(numVerticesNewLoc*dim, &verticesNewLoc, numVerticesNew, &verticesNewSorted);
304: for (i = 0, c = 0; i < numVerticesNew; i++) {
305: if (owners[i] == rank) {
306: for (j=0; j<dim; ++j) verticesNewLoc[dim*c+j] = verticesNew[dim*i+j];
307: c++;
308: }
309: }
311: /* Reorder for consistency with DMPlex */
312: for (i = 0; i < numCellsNew; ++i) DMPlexInvertCell(DM_POLYTOPE_TETRAHEDRON, &cellsNew[4*i]);
314: /* Create new plex */
315: DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNewLoc, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, verticesNewLoc, NULL, &verticesNewSorted, dmNew);
316: PMMG_Free_all(PMMG_ARG_start, PMMG_ARG_ppParMesh, &parmesh, PMMG_ARG_end);
317: PetscFree4(verticesNew, verTagsNew, corners, requiredVer);
319: /* Get adapted mesh information */
320: DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);
321: DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);
322: DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);
324: /* Rebuild boundary label */
325: DMCreateLabel(*dmNew, flg ? bdName : bdLabelName);
326: DMGetLabel(*dmNew, flg ? bdName : bdLabelName, &bdLabelNew);
327: for (i = 0; i < numFacesNew; i++) {
328: PetscBool hasTag = PETSC_FALSE;
329: PetscInt numCoveredPoints, numFaces = 0, facePoints[3];
330: const PetscInt *coveredPoints = NULL;
332: for (j = 0; j < dim; ++j) {
333: lv = facesNew[i*dim+j];
334: gv = gv_new[lv]-1;
335: PetscFindInt(gv, numVerticesNew, verticesNewSorted, &lv);
336: facePoints[j] = lv+vStart;
337: }
338: DMPlexGetFullJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);
339: for (j = 0; j < numCoveredPoints; ++j) {
340: if (coveredPoints[j] >= fStart && coveredPoints[j] < fEnd) {
341: numFaces++;
342: f = j;
343: }
344: }
346: DMLabelHasStratum(bdLabel, faceTagsNew[i], &hasTag);
347: if (hasTag) DMLabelSetValue(bdLabelNew, coveredPoints[f], faceTagsNew[i]);
348: DMPlexRestoreJoin(*dmNew, dim, facePoints, &numCoveredPoints, &coveredPoints);
349: }
350: PetscFree4(facesNew, faceTagsNew, ridges, requiredFaces);
351: PetscFree2(owners, gv_new);
352: PetscFree2(verticesNewLoc, verticesNewSorted);
353: if (flg) DMLabelDestroy(&bdLabel);
355: /* Rebuild cell labels */
356: DMCreateLabel(*dmNew, rgLabel ? rgLabelName : rgName);
357: DMGetLabel(*dmNew, rgLabel ? rgLabelName : rgName, &rgLabelNew);
358: for (c = cStart; c < cEnd; ++c) DMLabelSetValue(rgLabelNew, c, cellTagsNew[c-cStart]);
359: PetscFree3(cellsNew, cellTagsNew, requiredCells);
361: return 0;
362: }