Actual source code: pragmaticadapt.c

  1: #include <petsc/private/dmpleximpl.h>
  2: #include <pragmatic/cpragmatic.h>

  4: PETSC_EXTERN PetscErrorCode DMAdaptMetric_Pragmatic_Plex(DM dm, Vec vertexMetric, DMLabel bdLabel, DMLabel rgLabel, DM *dmNew)
  5: {
  6:   MPI_Comm    comm;
  7:   const char *bdName = "_boundary_";
  8: #if 0
  9:   DM                 odm = dm;
 10: #endif
 11:   DM                 udm, cdm;
 12:   DMLabel            bdLabelFull;
 13:   const char        *bdLabelName;
 14:   IS                 bdIS, globalVertexNum;
 15:   PetscSection       coordSection;
 16:   Vec                coordinates;
 17:   const PetscScalar *coords, *met;
 18:   const PetscInt    *bdFacesFull, *gV;
 19:   PetscInt          *bdFaces, *bdFaceIds, *l2gv;
 20:   PetscReal         *x, *y, *z, *metric;
 21:   PetscInt          *cells;
 22:   PetscInt           dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v;
 23:   PetscInt           off, maxConeSize, numBdFaces, f, bdSize, i, j, Nd;
 24:   PetscBool          flg, isotropic, uniform;
 25:   DMLabel            bdLabelNew;
 26:   PetscReal         *coordsNew;
 27:   PetscInt          *bdTags;
 28:   PetscReal         *xNew[3] = {NULL, NULL, NULL};
 29:   PetscInt          *cellsNew;
 30:   PetscInt           d, numCellsNew, numVerticesNew;
 31:   PetscInt           numCornersNew, fStart, fEnd;
 32:   PetscMPIInt        numProcs;

 34:   PetscFunctionBegin;

 36:   /* Check for FEM adjacency flags */
 37:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
 38:   PetscCallMPI(MPI_Comm_size(comm, &numProcs));
 39:   if (bdLabel) {
 40:     PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName));
 41:     PetscCall(PetscStrcmp(bdLabelName, bdName, &flg));
 42:     PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName);
 43:   }
 44:   PetscCheck(!rgLabel, comm, PETSC_ERR_ARG_WRONG, "Cannot currently preserve cell tags with Pragmatic");
 45: #if 0
 46:   /* Check for overlap by looking for cell in the SF */
 47:   if (!overlapped) {
 48:     PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm));
 49:     if (!dm) {dm = odm; PetscCall(PetscObjectReference((PetscObject) dm));}
 50:   }
 51: #endif

 53:   /* Get mesh information */
 54:   PetscCall(DMGetDimension(dm, &dim));
 55:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
 56:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
 57:   PetscCall(DMPlexUninterpolate(dm, &udm));
 58:   PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL));
 59:   numCells = cEnd - cStart;
 60:   if (numCells == 0) {
 61:     PetscMPIInt rank;

 63:     PetscCallMPI(MPI_Comm_rank(comm, &rank));
 64:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank);
 65:   }
 66:   numVertices = vEnd - vStart;
 67:   PetscCall(PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices * PetscSqr(dim), &metric, numCells * maxConeSize, &cells));

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

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

 79:   /* Get local-to-global vertex map */
 80:   PetscCall(PetscCalloc1(numVertices, &l2gv));
 81:   PetscCall(DMPlexGetVertexNumbering(udm, &globalVertexNum));
 82:   PetscCall(ISGetIndices(globalVertexNum, &gV));
 83:   for (v = 0, numLocVertices = 0; v < numVertices; ++v) {
 84:     if (gV[v] >= 0) ++numLocVertices;
 85:     l2gv[v] = gV[v] < 0 ? -(gV[v] + 1) : gV[v];
 86:   }
 87:   PetscCall(ISRestoreIndices(globalVertexNum, &gV));
 88:   PetscCall(DMDestroy(&udm));

 90:   /* Get vertex coordinate arrays */
 91:   PetscCall(DMGetCoordinateDM(dm, &cdm));
 92:   PetscCall(DMGetLocalSection(cdm, &coordSection));
 93:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
 94:   PetscCall(VecGetArrayRead(coordinates, &coords));
 95:   for (v = vStart; v < vEnd; ++v) {
 96:     PetscCall(PetscSectionGetOffset(coordSection, v, &off));
 97:     x[v - vStart] = PetscRealPart(coords[off + 0]);
 98:     if (dim > 1) y[v - vStart] = PetscRealPart(coords[off + 1]);
 99:     if (dim > 2) z[v - vStart] = PetscRealPart(coords[off + 2]);
100:   }
101:   PetscCall(VecRestoreArrayRead(coordinates, &coords));

103:   /* Get boundary mesh */
104:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull));
105:   PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull));
106:   PetscCall(DMLabelGetStratumIS(bdLabelFull, 1, &bdIS));
107:   PetscCall(DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces));
108:   PetscCall(ISGetIndices(bdIS, &bdFacesFull));
109:   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
110:     PetscInt *closure = NULL;
111:     PetscInt  closureSize, cl;

113:     PetscCall(DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
114:     for (cl = 0; cl < closureSize * 2; cl += 2) {
115:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize;
116:     }
117:     PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
118:   }
119:   PetscCall(PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds));
120:   for (f = 0, bdSize = 0; f < numBdFaces; ++f) {
121:     PetscInt *closure = NULL;
122:     PetscInt  closureSize, cl;

124:     PetscCall(DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
125:     for (cl = 0; cl < closureSize * 2; cl += 2) {
126:       if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart;
127:     }
128:     PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure));
129:     if (bdLabel) PetscCall(DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]));
130:     else bdFaceIds[f] = 1;
131:   }
132:   PetscCall(ISDestroy(&bdIS));
133:   PetscCall(DMLabelDestroy(&bdLabelFull));

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

155: #if 0
156:   /* Destroy overlap mesh */
157:   PetscCall(DMDestroy(&dm));
158: #endif
159:   /* Send to Pragmatic and remesh */
160:   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
161:   switch (dim) {
162:   case 2:
163:     PetscStackCallExternalVoid("pragmatic_2d_mpi_init", pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm));
164:     break;
165:   case 3:
166:     PetscStackCallExternalVoid("pragmatic_3d_mpi_init", pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm));
167:     break;
168:   default:
169:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
170:   }
171:   PetscStackCallExternalVoid("pragmatic_set_boundary", pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds));
172:   PetscStackCallExternalVoid("pragmatic_set_metric", pragmatic_set_metric(metric));
173:   PetscStackCallExternalVoid("pragmatic_adapt", pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0));
174:   PetscCall(PetscFree(l2gv));
175:   PetscCall(PetscFPTrapPop());

177:   /* Retrieve mesh from Pragmatic and create new plex */
178:   PetscStackCallExternalVoid("pragmatic_get_info_mpi", pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew));
179:   PetscCall(PetscMalloc1(numVerticesNew * dim, &coordsNew));
180:   switch (dim) {
181:   case 2:
182:     numCornersNew = 3;
183:     PetscCall(PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]));
184:     PetscStackCallExternalVoid("pragmatic_get_coords_2d_mpi", pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]));
185:     break;
186:   case 3:
187:     numCornersNew = 4;
188:     PetscCall(PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]));
189:     PetscStackCallExternalVoid("pragmatic_get_coords_3d_mpi", pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]));
190:     break;
191:   default:
192:     SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim);
193:   }
194:   for (v = 0; v < numVerticesNew; ++v) {
195:     for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v];
196:   }
197:   PetscCall(PetscMalloc1(numCellsNew * (dim + 1), &cellsNew));
198:   PetscStackCallExternalVoid("pragmatic_get_elements", pragmatic_get_elements(cellsNew));
199:   PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew));

201:   /* Rebuild boundary label */
202:   PetscStackCallExternalVoid("pragmatic_get_boundaryTags", pragmatic_get_boundaryTags(&bdTags));
203:   PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName));
204:   PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew));
205:   PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd));
206:   PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd));
207:   PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd));
208:   for (c = cStart; c < cEnd; ++c) {
209:     /* Only for simplicial meshes */
210:     coff = (c - cStart) * (dim + 1);

212:     /* d is the local cell number of the vertex opposite to the face we are marking */
213:     for (d = 0; d < dim + 1; ++d) {
214:       if (bdTags[coff + d]) {
215:         const PetscInt perm[4][4] = {
216:           {-1, -1, -1, -1},
217:           {-1, -1, -1, -1},
218:           {1,  2,  0,  -1},
219:           {3,  2,  1,  0 }
220:         }; /* perm[d] = face opposite */
221:         const PetscInt *cone;

223:         /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */
224:         PetscCall(DMPlexGetCone(*dmNew, c, &cone));
225:         PetscCall(DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d]));
226:       }
227:     }
228:   }

230:   /* Clean up */
231:   switch (dim) {
232:   case 2:
233:     PetscCall(PetscFree2(xNew[0], xNew[1]));
234:     break;
235:   case 3:
236:     PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2]));
237:     break;
238:   }
239:   PetscCall(PetscFree(cellsNew));
240:   PetscCall(PetscFree5(x, y, z, metric, cells));
241:   PetscCall(PetscFree2(bdFaces, bdFaceIds));
242:   PetscCall(PetscFree(coordsNew));
243:   PetscStackCallExternalVoid("pragmatic_finalize", pragmatic_finalize());
244:   PetscFunctionReturn(PETSC_SUCCESS);
245: }