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 switch (dim) { 161 case 2: 162 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm); 163 break; 164 case 3: 165 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 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 171 pragmatic_set_metric(metric); 172 pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0); 173 PetscCall(PetscFree(l2gv)); 174 175 /* Retrieve mesh from Pragmatic and create new plex */ 176 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 177 PetscCall(PetscMalloc1(numVerticesNew * dim, &coordsNew)); 178 switch (dim) { 179 case 2: 180 numCornersNew = 3; 181 PetscCall(PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1])); 182 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]); 183 break; 184 case 3: 185 numCornersNew = 4; 186 PetscCall(PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2])); 187 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]); 188 break; 189 default: 190 SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim); 191 } 192 for (v = 0; v < numVerticesNew; ++v) { 193 for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v]; 194 } 195 PetscCall(PetscMalloc1(numCellsNew * (dim + 1), &cellsNew)); 196 pragmatic_get_elements(cellsNew); 197 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew)); 198 199 /* Rebuild boundary label */ 200 pragmatic_get_boundaryTags(&bdTags); 201 PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName)); 202 PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew)); 203 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 204 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 205 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 206 for (c = cStart; c < cEnd; ++c) { 207 /* Only for simplicial meshes */ 208 coff = (c - cStart) * (dim + 1); 209 210 /* d is the local cell number of the vertex opposite to the face we are marking */ 211 for (d = 0; d < dim + 1; ++d) { 212 if (bdTags[coff + d]) { 213 const PetscInt perm[4][4] = { 214 {-1, -1, -1, -1}, 215 {-1, -1, -1, -1}, 216 {1, 2, 0, -1}, 217 {3, 2, 1, 0 } 218 }; /* perm[d] = face opposite */ 219 const PetscInt *cone; 220 221 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 222 PetscCall(DMPlexGetCone(*dmNew, c, &cone)); 223 PetscCall(DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d])); 224 } 225 } 226 } 227 228 /* Clean up */ 229 switch (dim) { 230 case 2: 231 PetscCall(PetscFree2(xNew[0], xNew[1])); 232 break; 233 case 3: 234 PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2])); 235 break; 236 } 237 PetscCall(PetscFree(cellsNew)); 238 PetscCall(PetscFree5(x, y, z, metric, cells)); 239 PetscCall(PetscFree2(bdFaces, bdFaceIds)); 240 PetscCall(PetscFree(coordsNew)); 241 pragmatic_finalize(); 242 PetscFunctionReturn(0); 243 } 244