1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #include <pragmatic/cpragmatic.h> 3 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; 33 34 PetscFunctionBegin; 35 /* Check for FEM adjacency flags */ 36 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 37 PetscCallMPI(MPI_Comm_size(comm, &numProcs)); 38 if (bdLabel) { 39 PetscCall(PetscObjectGetName((PetscObject)bdLabel, &bdLabelName)); 40 PetscCall(PetscStrcmp(bdLabelName, bdName, &flg)); 41 PetscCheck(!flg, comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 42 } 43 PetscCheck(!rgLabel, comm, PETSC_ERR_ARG_WRONG, "Cannot currently preserve cell tags with Pragmatic"); 44 #if 0 45 /* Check for overlap by looking for cell in the SF */ 46 if (!overlapped) { 47 PetscCall(DMPlexDistributeOverlap(odm, 1, NULL, &dm)); 48 if (!dm) {dm = odm; PetscCall(PetscObjectReference((PetscObject) dm));} 49 } 50 #endif 51 52 /* Get mesh information */ 53 PetscCall(DMGetDimension(dm, &dim)); 54 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 55 PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd)); 56 PetscCall(DMPlexUninterpolate(dm, &udm)); 57 PetscCall(DMPlexGetMaxSizes(udm, &maxConeSize, NULL)); 58 numCells = cEnd - cStart; 59 if (numCells == 0) { 60 PetscMPIInt rank; 61 62 PetscCallMPI(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 PetscCall(PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices * PetscSqr(dim), &metric, numCells * maxConeSize, &cells)); 67 68 /* Get cell offsets */ 69 for (c = 0, coff = 0; c < numCells; ++c) { 70 const PetscInt *cone; 71 PetscInt coneSize, cl; 72 73 PetscCall(DMPlexGetConeSize(udm, c, &coneSize)); 74 PetscCall(DMPlexGetCone(udm, c, &cone)); 75 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 76 } 77 78 /* Get local-to-global vertex map */ 79 PetscCall(PetscCalloc1(numVertices, &l2gv)); 80 PetscCall(DMPlexGetVertexNumbering(udm, &globalVertexNum)); 81 PetscCall(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 PetscCall(ISRestoreIndices(globalVertexNum, &gV)); 87 PetscCall(DMDestroy(&udm)); 88 89 /* Get vertex coordinate arrays */ 90 PetscCall(DMGetCoordinateDM(dm, &cdm)); 91 PetscCall(DMGetLocalSection(cdm, &coordSection)); 92 PetscCall(DMGetCoordinatesLocal(dm, &coordinates)); 93 PetscCall(VecGetArrayRead(coordinates, &coords)); 94 for (v = vStart; v < vEnd; ++v) { 95 PetscCall(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 PetscCall(VecRestoreArrayRead(coordinates, &coords)); 101 102 /* Get boundary mesh */ 103 PetscCall(DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull)); 104 PetscCall(DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull)); 105 PetscCall(DMLabelGetStratumIS(bdLabelFull, 1, &bdIS)); 106 PetscCall(DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces)); 107 PetscCall(ISGetIndices(bdIS, &bdFacesFull)); 108 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 109 PetscInt *closure = NULL; 110 PetscInt closureSize, cl; 111 112 PetscCall(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 PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure)); 117 } 118 PetscCall(PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds)); 119 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 120 PetscInt *closure = NULL; 121 PetscInt closureSize, cl; 122 123 PetscCall(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 PetscCall(DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure)); 128 if (bdLabel) PetscCall(DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f])); 129 else bdFaceIds[f] = 1; 130 } 131 PetscCall(ISDestroy(&bdIS)); 132 PetscCall(DMLabelDestroy(&bdLabelFull)); 133 134 /* Get metric */ 135 PetscCall(VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view")); 136 PetscCall(VecGetArrayRead(vertexMetric, &met)); 137 PetscCall(DMPlexMetricIsIsotropic(dm, &isotropic)); 138 PetscCall(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 PetscCall(VecRestoreArrayRead(vertexMetric, &met)); 153 154 #if 0 155 /* Destroy overlap mesh */ 156 PetscCall(DMDestroy(&dm)); 157 #endif 158 /* Send to Pragmatic and remesh */ 159 PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); 160 switch (dim) { 161 case 2: 162 PetscStackCallExternalVoid("pragmatic_2d_mpi_init", pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm)); 163 break; 164 case 3: 165 PetscStackCallExternalVoid("pragmatic_3d_mpi_init", pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm)); 166 break; 167 default: 168 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim); 169 } 170 PetscStackCallExternalVoid("pragmatic_set_boundary", pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds)); 171 PetscStackCallExternalVoid("pragmatic_set_metric", pragmatic_set_metric(metric)); 172 PetscStackCallExternalVoid("pragmatic_adapt", pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0)); 173 PetscCall(PetscFree(l2gv)); 174 PetscCall(PetscFPTrapPop()); 175 176 /* Retrieve mesh from Pragmatic and create new plex */ 177 PetscStackCallExternalVoid("pragmatic_get_info_mpi", pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew)); 178 PetscCall(PetscMalloc1(numVerticesNew * dim, &coordsNew)); 179 switch (dim) { 180 case 2: 181 numCornersNew = 3; 182 PetscCall(PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1])); 183 PetscStackCallExternalVoid("pragmatic_get_coords_2d_mpi", pragmatic_get_coords_2d_mpi(xNew[0], xNew[1])); 184 break; 185 case 3: 186 numCornersNew = 4; 187 PetscCall(PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2])); 188 PetscStackCallExternalVoid("pragmatic_get_coords_3d_mpi", pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2])); 189 break; 190 default: 191 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim); 192 } 193 for (v = 0; v < numVerticesNew; ++v) { 194 for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v]; 195 } 196 PetscCall(PetscMalloc1(numCellsNew * (dim + 1), &cellsNew)); 197 PetscStackCallExternalVoid("pragmatic_get_elements", pragmatic_get_elements(cellsNew)); 198 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew)); 199 200 /* Rebuild boundary label */ 201 PetscStackCallExternalVoid("pragmatic_get_boundaryTags", pragmatic_get_boundaryTags(&bdTags)); 202 PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName)); 203 PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew)); 204 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 205 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 206 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 207 for (c = cStart; c < cEnd; ++c) { 208 /* Only for simplicial meshes */ 209 coff = (c - cStart) * (dim + 1); 210 211 /* d is the local cell number of the vertex opposite to the face we are marking */ 212 for (d = 0; d < dim + 1; ++d) { 213 if (bdTags[coff + d]) { 214 const PetscInt perm[4][4] = { 215 {-1, -1, -1, -1}, 216 {-1, -1, -1, -1}, 217 {1, 2, 0, -1}, 218 {3, 2, 1, 0 } 219 }; /* perm[d] = face opposite */ 220 const PetscInt *cone; 221 222 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 223 PetscCall(DMPlexGetCone(*dmNew, c, &cone)); 224 PetscCall(DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d])); 225 } 226 } 227 } 228 229 /* Clean up */ 230 switch (dim) { 231 case 2: 232 PetscCall(PetscFree2(xNew[0], xNew[1])); 233 break; 234 case 3: 235 PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2])); 236 break; 237 } 238 PetscCall(PetscFree(cellsNew)); 239 PetscCall(PetscFree5(x, y, z, metric, cells)); 240 PetscCall(PetscFree2(bdFaces, bdFaceIds)); 241 PetscCall(PetscFree(coordsNew)); 242 PetscStackCallExternalVoid("pragmatic_finalize", pragmatic_finalize()); 243 PetscFunctionReturn(PETSC_SUCCESS); 244 } 245