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