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: SETERRQ(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", 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 PetscCall(PetscFree(l2gv)); 173 174 /* Retrieve mesh from Pragmatic and create new plex */ 175 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 176 PetscCall(PetscMalloc1(numVerticesNew*dim, &coordsNew)); 177 switch (dim) { 178 case 2: 179 numCornersNew = 3; 180 PetscCall(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 PetscCall(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 %D", dim); 190 } 191 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];} 192 PetscCall(PetscMalloc1(numCellsNew*(dim+1), &cellsNew)); 193 pragmatic_get_elements(cellsNew); 194 PetscCall(DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew)); 195 196 /* Rebuild boundary label */ 197 pragmatic_get_boundaryTags(&bdTags); 198 PetscCall(DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName)); 199 PetscCall(DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew)); 200 PetscCall(DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd)); 201 PetscCall(DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd)); 202 PetscCall(DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd)); 203 for (c = cStart; c < cEnd; ++c) { 204 205 /* Only for simplicial meshes */ 206 coff = (c-cStart)*(dim+1); 207 208 /* d is the local cell number of the vertex opposite to the face we are marking */ 209 for (d = 0; d < dim+1; ++d) { 210 if (bdTags[coff+d]) { 211 const PetscInt perm[4][4] = {{-1, -1, -1, -1}, {-1, -1, -1, -1}, {1, 2, 0, -1}, {3, 2, 1, 0}}; /* 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])); 224 break; 225 case 3: PetscCall(PetscFree3(xNew[0], xNew[1], xNew[2])); 226 break; 227 } 228 PetscCall(PetscFree(cellsNew)); 229 PetscCall(PetscFree5(x, y, z, metric, cells)); 230 PetscCall(PetscFree2(bdFaces, bdFaceIds)); 231 PetscCall(PetscFree(coordsNew)); 232 pragmatic_finalize(); 233 PetscFunctionReturn(0); 234 } 235