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