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 MPI_Comm comm; 6 const char *bdName = "_boundary_"; 7 #if 0 8 DM odm = dm; 9 #endif 10 DM udm, cdm; 11 DMLabel bdLabelFull; 12 const char *bdLabelName; 13 IS bdIS, globalVertexNum; 14 PetscSection coordSection; 15 Vec coordinates; 16 const PetscScalar *coords, *met; 17 const PetscInt *bdFacesFull, *gV; 18 PetscInt *bdFaces, *bdFaceIds, *l2gv; 19 PetscReal *x, *y, *z, *metric; 20 PetscInt *cells; 21 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v; 22 PetscInt off, maxConeSize, numBdFaces, f, bdSize, i, j, Nd; 23 PetscBool flg, isotropic, uniform; 24 DMLabel bdLabelNew; 25 PetscReal *coordsNew; 26 PetscInt *bdTags; 27 PetscReal *xNew[3] = {NULL, NULL, NULL}; 28 PetscInt *cellsNew; 29 PetscInt d, numCellsNew, numVerticesNew; 30 PetscInt numCornersNew, fStart, fEnd; 31 PetscMPIInt numProcs; 32 33 PetscFunctionBegin; 34 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 switch (dim) { 160 case 2: pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm); break; 161 case 3: pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm); break; 162 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim); 163 } 164 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 165 pragmatic_set_metric(metric); 166 pragmatic_adapt(((DM_Plex *)dm->data)->remeshBd ? 1 : 0); 167 PetscCall(PetscFree(l2gv)); 168 169 /* Retrieve mesh from Pragmatic and create new plex */ 170 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 171 PetscCall(PetscMalloc1(numVerticesNew * dim, &coordsNew)); 172 switch (dim) { 173 case 2: 174 numCornersNew = 3; 175 PetscCall(PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1])); 176 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]); 177 break; 178 case 3: 179 numCornersNew = 4; 180 PetscCall(PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2])); 181 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]); 182 break; 183 default: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %" PetscInt_FMT, dim); 184 } 185 for (v = 0; v < numVerticesNew; ++v) { 186 for (d = 0; d < dim; ++d) coordsNew[v * dim + d] = xNew[d][v]; 187 } 188 PetscCall(PetscMalloc1(numCellsNew * (dim + 1), &cellsNew)); 189 pragmatic_get_elements(cellsNew); 190 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew)); 191 192 /* Rebuild boundary label */ 193 pragmatic_get_boundaryTags(&bdTags); 194 PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName)); 195 PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew)); 196 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 197 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 198 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 199 for (c = cStart; c < cEnd; ++c) { 200 /* Only for simplicial meshes */ 201 coff = (c - cStart) * (dim + 1); 202 203 /* d is the local cell number of the vertex opposite to the face we are marking */ 204 for (d = 0; d < dim + 1; ++d) { 205 if (bdTags[coff + d]) { 206 const PetscInt perm[4][4] = { 207 {-1, -1, -1, -1}, 208 {-1, -1, -1, -1}, 209 {1, 2, 0, -1}, 210 {3, 2, 1, 0 } 211 }; /* perm[d] = face opposite */ 212 const PetscInt *cone; 213 214 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 215 PetscCall(DMPlexGetCone(*dmNew, c, &cone)); 216 PetscCall(DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff + d])); 217 } 218 } 219 } 220 221 /* Clean up */ 222 switch (dim) { 223 case 2: PetscCall(PetscFree2(xNew[0], xNew[1])); break; 224 case 3: PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2])); break; 225 } 226 PetscCall(PetscFree(cellsNew)); 227 PetscCall(PetscFree5(x, y, z, metric, cells)); 228 PetscCall(PetscFree2(bdFaces, bdFaceIds)); 229 PetscCall(PetscFree(coordsNew)); 230 pragmatic_finalize(); 231 PetscFunctionReturn(0); 232 } 233