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 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break; 151 case 3: 152 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);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(remeshBd ? 1 : 0); 158 ierr = PetscFree(l2gv);CHKERRQ(ierr); 159 /* Read out mesh */ 160 pragmatic_get_info(&numVerticesNew, &numCellsNew); 161 ierr = PetscMalloc1(numVerticesNew*dim, &coordsNew);CHKERRQ(ierr); 162 switch (dim) { 163 case 2: 164 numCornersNew = 3; 165 ierr = PetscMalloc2(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1]);CHKERRQ(ierr); 166 pragmatic_get_coords_2d(xNew[0], xNew[1]); 167 break; 168 case 3: 169 numCornersNew = 4; 170 ierr = PetscMalloc3(numVerticesNew, &xNew[0], numVerticesNew, &xNew[1], numVerticesNew, &xNew[2]);CHKERRQ(ierr); 171 pragmatic_get_coords_3d(xNew[0], xNew[1], xNew[2]); 172 break; 173 default: 174 SETERRQ1(comm, PETSC_ERR_ARG_OUTOFRANGE, "No Pragmatic adaptation defined for dimension %d", dim); 175 } 176 for (v = 0; v < numVerticesNew; ++v) {for (d = 0; d < dim; ++d) coordsNew[v*dim+d] = xNew[d][v];} 177 ierr = PetscMalloc1(numCellsNew*(dim+1), &cellsNew);CHKERRQ(ierr); 178 pragmatic_get_elements(cellsNew); 179 ierr = DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr); 180 /* Read out boundary label */ 181 pragmatic_get_boundaryTags(&bdTags); 182 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 183 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 184 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 185 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 186 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 187 for (c = cStart; c < cEnd; ++c) { 188 /* Only for simplicial meshes */ 189 coff = (c-cStart)*(dim+1); 190 for (d = 0; d < dim+1; ++d) { 191 if (bdTags[coff+d]) { 192 const PetscInt *faces; 193 PetscInt numFaces, vertices[3], p; 194 195 /* Mark face opposite to this vertex */ 196 for (p = 0, v = 0; p < dim+1; ++p) if (p != d) {vertices[v] = cellsNew[coff+p] + vStart; ++v;} 197 ierr = DMPlexGetFullJoin(*dmNew, dim, vertices, &numFaces, &faces);CHKERRQ(ierr); 198 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]); 199 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]); 200 ierr = DMLabelSetValue(bdLabelNew, faces[0], bdTags[coff+d]);CHKERRQ(ierr); 201 ierr = DMPlexRestoreJoin(*dmNew, dim, vertices, &numFaces, &faces);CHKERRQ(ierr); 202 } 203 } 204 } 205 /* Cleanup */ 206 switch (dim) { 207 case 2: ierr = PetscFree2(xNew[0], xNew[1]);CHKERRQ(ierr);break; 208 case 3: ierr = PetscFree3(xNew[0], xNew[1], xNew[2]);CHKERRQ(ierr);break; 209 } 210 ierr = PetscFree(cellsNew);CHKERRQ(ierr); 211 ierr = PetscFree5(x, y, z, metric, cells);CHKERRQ(ierr); 212 ierr = PetscFree2(bdFaces, bdFaceIds);CHKERRQ(ierr); 213 ierr = PetscFree(coordsNew);CHKERRQ(ierr); 214 pragmatic_finalize(); 215 #else 216 SETERRQ(comm, PETSC_ERR_SUP, "Remeshing needs external package support.\nPlease reconfigure with --download-pragmatic."); 217 #endif 218 PetscFunctionReturn(0); 219 } 220 221 /*@ 222 DMPlexAdapt - Generates a mesh adapted to the specified metric field using the pragmatic library. 223 224 Input Parameters: 225 + dm - The DM object 226 . metric - The metric to which the mesh is adapted, defined vertex-wise. 227 - 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". 228 229 Output Parameter: 230 . dmAdapt - Pointer to the DM object containing the adapted mesh 231 232 Level: advanced 233 234 .seealso: DMCoarsen(), DMRefine() 235 @*/ 236 PetscErrorCode DMPlexAdapt(DM dm, Vec metric, const char bdLabelName[], DM *dmAdapt) 237 { 238 DM_Plex *mesh = (DM_Plex *) dm->data; 239 PetscErrorCode ierr; 240 241 PetscFunctionBegin; 242 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 243 PetscValidHeaderSpecific(metric, VEC_CLASSID, 2); 244 PetscValidCharPointer(bdLabelName, 3); 245 PetscValidPointer(dmAdapt, 4); 246 ierr = DMPlexRemesh_Internal(dm, metric, bdLabelName, mesh->remeshBd, dmAdapt);CHKERRQ(ierr); 247 PetscFunctionReturn(0); 248 } 249