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;


 35:   /* Check for FEM adjacency flags */
 36:   PetscObjectGetComm((PetscObject)dm, &comm);
 37:   MPI_Comm_size(comm, &numProcs);
 38:   if (bdLabel) {
 39:     PetscObjectGetName((PetscObject)bdLabel, &bdLabelName);
 40:     PetscStrcmp(bdLabelName, bdName, &flg);
 42:   }
 44: #if 0
 45:   /* Check for overlap by looking for cell in the SF */
 46:   if (!overlapped) {
 47:     DMPlexDistributeOverlap(odm, 1, NULL, &dm);
 48:     if (!dm) {dm = odm; PetscObjectReference((PetscObject) dm);}
 49:   }
 50: #endif

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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