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; 24 PetscBool flg; 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 PetscErrorCode ierr; 34 35 PetscFunctionBegin; 36 37 /* Check for FEM adjacency flags */ 38 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 39 ierr = MPI_Comm_size(comm, &numProcs);CHKERRMPI(ierr); 40 if (bdLabel) { 41 ierr = PetscObjectGetName((PetscObject) bdLabel, &bdLabelName);CHKERRQ(ierr); 42 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 43 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 44 } 45 if (rgLabel) SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Cannot currently preserve cell tags with Pragmatic"); 46 #if 0 47 /* Check for overlap by looking for cell in the SF */ 48 if (!overlapped) { 49 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr); 50 if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);} 51 } 52 #endif 53 54 /* Get mesh information */ 55 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 56 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 57 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 58 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 59 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 60 numCells = cEnd - cStart; 61 if (numCells == 0) { 62 PetscMPIInt rank; 63 64 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 65 SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot perform mesh adaptation because process %d does not own any cells.", rank); 66 } 67 numVertices = vEnd - vStart; 68 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 69 70 /* Get cell offsets */ 71 for (c = 0, coff = 0; c < numCells; ++c) { 72 const PetscInt *cone; 73 PetscInt coneSize, cl; 74 75 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 76 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 77 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 78 } 79 80 /* Get local-to-global vertex map */ 81 ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr); 82 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 83 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 84 for (v = 0, numLocVertices = 0; v < numVertices; ++v) { 85 if (gV[v] >= 0) ++numLocVertices; 86 l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v]; 87 } 88 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 89 ierr = DMDestroy(&udm);CHKERRQ(ierr); 90 91 /* Get vertex coordinate arrays */ 92 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 93 ierr = DMGetLocalSection(cdm, &coordSection);CHKERRQ(ierr); 94 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 95 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 96 for (v = vStart; v < vEnd; ++v) { 97 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 98 x[v-vStart] = PetscRealPart(coords[off+0]); 99 if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]); 100 if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]); 101 } 102 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 103 104 /* Get boundary mesh */ 105 ierr = DMLabelCreate(PETSC_COMM_SELF, bdName, &bdLabelFull);CHKERRQ(ierr); 106 ierr = DMPlexMarkBoundaryFaces(dm, 1, bdLabelFull);CHKERRQ(ierr); 107 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 108 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 109 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 110 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 111 PetscInt *closure = NULL; 112 PetscInt closureSize, cl; 113 114 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 115 for (cl = 0; cl < closureSize*2; cl += 2) { 116 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 117 } 118 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 119 } 120 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 121 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 122 PetscInt *closure = NULL; 123 PetscInt closureSize, cl; 124 125 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 126 for (cl = 0; cl < closureSize*2; cl += 2) { 127 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 128 } 129 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 130 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 131 else {bdFaceIds[f] = 1;} 132 } 133 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 134 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 135 136 /* Get metric */ 137 ierr = VecViewFromOptions(vertexMetric, NULL, "-adapt_metric_view");CHKERRQ(ierr); 138 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 139 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]); 140 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 141 142 #if 0 143 /* Destroy overlap mesh */ 144 ierr = DMDestroy(&dm);CHKERRQ(ierr); 145 #endif 146 /* Send to Pragmatic and remesh */ 147 switch (dim) { 148 case 2: 149 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm); 150 break; 151 case 3: 152 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm); 153 break; 154 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim); 155 } 156 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 157 pragmatic_set_metric(metric); 158 pragmatic_adapt(((DM_Plex *) dm->data)->remeshBd ? 1 : 0); 159 ierr = PetscFree(l2gv);CHKERRQ(ierr); 160 161 /* Retrieve mesh from Pragmatic and create new plex */ 162 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 163 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr); 164 switch (dim) { 165 case 2: 166 numCornersNew = 3; 167 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr); 168 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]); 169 break; 170 case 3: 171 numCornersNew = 4; 172 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr); 173 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]); 174 break; 175 default: 176 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %D", dim); 177 } 178 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];} 179 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr); 180 pragmatic_get_elements(cellsNew); 181 ierr = DMPlexCreateFromCellListParallelPetsc(comm, dim, numCellsNew, numVerticesNew, PETSC_DECIDE, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, NULL, dmNew);CHKERRQ(ierr); 182 183 /* Rebuild boundary label */ 184 pragmatic_get_boundaryTags(&bdTags); 185 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 186 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 187 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 188 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 189 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 190 for (c = cStart; c < cEnd; ++c) { 191 192 /* Only for simplicial meshes */ 193 coff = (c-cStart)*(dim+1); 194 195 /* d is the local cell number of the vertex opposite to the face we are marking */ 196 for (d = 0; d < dim+1; ++d) { 197 if (bdTags[coff+d]) { 198 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 */ 199 const PetscInt *cone; 200 201 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 202 ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr); 203 ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr); 204 } 205 } 206 } 207 208 /* Clean up */ 209 switch (dim) { 210 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr); 211 break; 212 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr); 213 break; 214 } 215 ierr = PetscFree(cellsNew);CHKERRQ(ierr); 216 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr); 217 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 218 ierr = PetscFree(coordsNew);CHKERRQ(ierr); 219 pragmatic_finalize(); 220 PetscFunctionReturn(0); 221 } 222