1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 #ifdef PETSC_HAVE_PRAGMATIC 3 #include <pragmatic/cpragmatic.h> 4 #endif 5 6 /* 7 DMPlexRemesh_Internal - Generates a new mesh conforming to a metric field. 8 9 Input Parameters: 10 + dm - The DM object 11 . vertexMetric - The metric to which the mesh is adapted, defined vertex-wise in a LOCAL vector 12 . remeshBd - Flag to allow boundary changes 13 - bdLabelName - Label name for boundary tags which are preserved in dmNew, or NULL. Should not be "_boundary_". 14 15 Output Parameter: 16 . dmNew - the new DM 17 18 Level: advanced 19 20 .seealso: DMCoarsen(), DMRefine() 21 */ 22 PetscErrorCode DMPlexRemesh_Internal(DM dm, Vec vertexMetric, const char bdLabelName[], PetscBool remeshBd, DM *dmNew) 23 { 24 MPI_Comm comm; 25 const char *bdName = "_boundary_"; 26 #if 0 27 DM odm = dm; 28 #endif 29 DM udm, cdm; 30 DMLabel bdLabel = NULL, bdLabelFull; 31 IS bdIS, globalVertexNum; 32 PetscSection coordSection; 33 Vec coordinates; 34 const PetscScalar *coords, *met; 35 const PetscInt *bdFacesFull, *gV; 36 PetscInt *bdFaces, *bdFaceIds, *l2gv; 37 PetscReal *x, *y, *z, *metric; 38 PetscInt *cells; 39 PetscInt dim, cStart, cEnd, numCells, c, coff, vStart, vEnd, numVertices, numLocVertices, v; 40 PetscInt off, maxConeSize, numBdFaces, f, bdSize; 41 PetscBool flg; 42 #ifdef PETSC_HAVE_PRAGMATIC 43 DMLabel bdLabelNew; 44 double *coordsNew; 45 PetscInt *bdTags; 46 PetscReal *xNew[3] = {NULL, NULL, NULL}; 47 PetscInt *cellsNew; 48 PetscInt d, numCellsNew, numVerticesNew; 49 PetscInt numCornersNew, fStart, fEnd; 50 #endif 51 PetscMPIInt numProcs; 52 PetscErrorCode ierr; 53 54 PetscFunctionBegin; 55 /* Check for FEM adjacency flags */ 56 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 57 PetscValidHeaderSpecific(vertexMetric, VEC_CLASSID, 2); 58 PetscValidCharPointer(bdLabelName, 3); 59 PetscValidPointer(dmNew, 5); 60 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 61 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr); 62 if (bdLabelName) { 63 size_t len; 64 65 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 66 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 67 ierr = PetscStrlen(bdLabelName, &len);CHKERRQ(ierr); 68 if (len) { 69 ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr); 70 if (!bdLabel) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Label \"%s\" does not exist in DM", bdLabelName); 71 } 72 } 73 /* Add overlap for Pragmatic */ 74 #if 0 75 /* Check for overlap by looking for cell in the SF */ 76 if (!overlapped) { 77 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr); 78 if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);} 79 } 80 #endif 81 /* Get mesh information */ 82 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 83 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 84 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 85 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 86 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 87 numCells = cEnd - cStart; 88 numVertices = vEnd - vStart; 89 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 90 for (c = 0, coff = 0; c < numCells; ++c) { 91 const PetscInt *cone; 92 PetscInt coneSize, cl; 93 94 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 95 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 96 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 97 } 98 ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr); 99 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 100 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 101 for (v = 0, numLocVertices = 0; v < numVertices; ++v) { 102 if (gV[v] >= 0) ++numLocVertices; 103 l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v]; 104 } 105 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 106 ierr = DMDestroy(&udm);CHKERRQ(ierr); 107 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 108 ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr); 109 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 110 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 111 for (v = vStart; v < vEnd; ++v) { 112 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 113 x[v-vStart] = PetscRealPart(coords[off+0]); 114 if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]); 115 if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]); 116 } 117 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 118 /* Get boundary mesh */ 119 ierr = DMLabelCreate(bdName, &bdLabelFull);CHKERRQ(ierr); 120 ierr = DMPlexMarkBoundaryFaces(dm, bdLabelFull);CHKERRQ(ierr); 121 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 122 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 123 ierr = ISGetIndices(bdIS, &bdFacesFull);CHKERRQ(ierr); 124 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 125 PetscInt *closure = NULL; 126 PetscInt closureSize, cl; 127 128 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 129 for (cl = 0; cl < closureSize*2; cl += 2) { 130 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) ++bdSize; 131 } 132 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 133 } 134 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 135 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 136 PetscInt *closure = NULL; 137 PetscInt closureSize, cl; 138 139 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 140 for (cl = 0; cl < closureSize*2; cl += 2) { 141 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 142 } 143 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 144 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 145 else {bdFaceIds[f] = 1;} 146 } 147 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 148 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 149 /* Get metric */ 150 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 151 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]); 152 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 153 #if 0 154 /* Destroy overlap mesh */ 155 ierr = DMDestroy(&dm);CHKERRQ(ierr); 156 #endif 157 /* Create new mesh */ 158 #ifdef PETSC_HAVE_PRAGMATIC 159 switch (dim) { 160 case 2: 161 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break; 162 case 3: 163 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break; 164 default: SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 165 } 166 pragmatic_set_boundary(&numBdFaces, bdFaces, bdFaceIds); 167 pragmatic_set_metric(metric); 168 pragmatic_adapt(remeshBd ? 1 : 0); 169 ierr = PetscFree(l2gv);CHKERRQ(ierr); 170 /* Read out mesh */ 171 pragmatic_get_info_mpi(&numVerticesNew, &numCellsNew); 172 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr); 173 switch (dim) { 174 case 2: 175 numCornersNew = 3; 176 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr); 177 pragmatic_get_coords_2d_mpi(xNew[0], xNew[1]); 178 break; 179 case 3: 180 numCornersNew = 4; 181 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr); 182 pragmatic_get_coords_3d_mpi(xNew[0], xNew[1], xNew[2]); 183 break; 184 default: 185 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 186 } 187 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = (double) xNew[d][v];} 188 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr); 189 pragmatic_get_elements(cellsNew); 190 ierr = DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr); 191 /* Read out boundary label */ 192 pragmatic_get_boundaryTags(&bdTags); 193 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 194 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 195 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 196 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 197 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 198 for (c = cStart; c < cEnd; ++c) { 199 /* Only for simplicial meshes */ 200 coff = (c-cStart)*(dim+1); 201 /* d is the local cell number of the vertex opposite to the face we are marking */ 202 for (d = 0; d < dim+1; ++d) { 203 if (bdTags[coff+d]) { 204 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 */ 205 const PetscInt *cone; 206 207 /* Mark face opposite to this vertex: This pattern is specified in DMPlexGetRawFaces_Internal() */ 208 ierr = DMPlexGetCone(*dmNew, c, &cone);CHKERRQ(ierr); 209 ierr = DMLabelSetValue(bdLabelNew, cone[perm[dim][d]], bdTags[coff+d]);CHKERRQ(ierr); 210 } 211 } 212 } 213 /* Cleanup */ 214 switch (dim) { 215 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break; 216 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break; 217 } 218 ierr = PetscFree(cellsNew);CHKERRQ(ierr); 219 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr); 220 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 221 ierr = PetscFree(coordsNew);CHKERRQ(ierr); 222 pragmatic_finalize(); 223 #else 224 SETERRQ(comm, PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic."); 225 #endif 226 PetscFunctionReturn(0); 227 } 228 229 /*@C 230 DMPlexAdapt - Generates a mesh adapted to the specified metric field using the pragmatic library. 231 232 Input Parameters: 233 + dm - The DM object 234 . metric - The metric to which the mesh is adapted, defined vertex-wise. 235 - bdyLabelName - Label name for boundary tags. These will be preserved in the output mesh. bdyLabelName should be "" (empty string) if there is no such label, and should be different from "boundary". 236 237 Output Parameter: 238 . dmAdapt - Pointer to the DM object containing the adapted mesh 239 240 Level: advanced 241 242 .seealso: DMCoarsen(), DMRefine() 243 @*/ 244 PetscErrorCode DMPlexAdapt(DM dm, Vec metric, const char bdLabelName[], DM *dmAdapt) 245 { 246 DM_Plex *mesh = (DM_Plex *) dm->data; 247 PetscErrorCode ierr; 248 249 PetscFunctionBegin; 250 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 251 PetscValidHeaderSpecific(metric, VEC_CLASSID, 2); 252 PetscValidCharPointer(bdLabelName, 3); 253 PetscValidPointer(dmAdapt, 4); 254 ierr = DMPlexRemesh_Internal(dm, metric, bdLabelName, mesh->remeshBd, dmAdapt);CHKERRQ(ierr); 255 PetscFunctionReturn(0); 256 } 257