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 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 52 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; 62 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)); 68 69 /* Get cell offsets */ 70 for (c = 0, coff = 0; c < numCells; ++c) { 71 const PetscInt *cone; 72 PetscInt coneSize, cl; 73 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 } 78 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)); 89 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)); 102 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; 112 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; 123 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)); 134 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)); 154 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()); 176 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)); 200 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); 211 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; 222 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 } 229 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 } 246