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 /* Check for FEM adjacency flags */ 53 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 54 PetscValidHeaderSpecific(vertexMetric, VEC_CLASSID, 2); 55 PetscValidCharPointer(bdLabelName, 3); 56 PetscValidPointer(dmNew, 5); 57 ierr = PetscObjectGetComm((PetscObject) dm, &comm);CHKERRQ(ierr); 58 ierr = MPI_Comm_size(comm, &numProcs);CHKERRQ(ierr); 59 if (bdLabelName) { 60 size_t len; 61 62 ierr = PetscStrcmp(bdLabelName, bdName, &flg);CHKERRQ(ierr); 63 if (flg) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "\"%s\" cannot be used as label for boundary facets", bdLabelName); 64 ierr = PetscStrlen(bdLabelName, &len);CHKERRQ(ierr); 65 if (len) { 66 ierr = DMGetLabel(dm, bdLabelName, &bdLabel);CHKERRQ(ierr); 67 if (!bdLabel) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Label \"%s\" does not exist in DM", bdLabelName); 68 } 69 } 70 /* Add overlap for Pragmatic */ 71 #if 0 72 /* Check for overlap by looking for cell in the SF */ 73 if (!overlapped) { 74 ierr = DMPlexDistributeOverlap(odm, 1, NULL, &dm);CHKERRQ(ierr); 75 if (!dm) {dm = odm; ierr = PetscObjectReference((PetscObject) dm);CHKERRQ(ierr);} 76 } 77 #endif 78 /* Get mesh information */ 79 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 80 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 81 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 82 ierr = DMPlexUninterpolate(dm, &udm);CHKERRQ(ierr); 83 ierr = DMPlexGetMaxSizes(udm, &maxConeSize, NULL);CHKERRQ(ierr); 84 numCells = cEnd - cStart; 85 numVertices = vEnd - vStart; 86 ierr = PetscCalloc5(numVertices, &x, numVertices, &y, numVertices, &z, numVertices*PetscSqr(dim), &metric, numCells*maxConeSize, &cells);CHKERRQ(ierr); 87 for (c = 0, coff = 0; c < numCells; ++c) { 88 const PetscInt *cone; 89 PetscInt coneSize, cl; 90 91 ierr = DMPlexGetConeSize(udm, c, &coneSize);CHKERRQ(ierr); 92 ierr = DMPlexGetCone(udm, c, &cone);CHKERRQ(ierr); 93 for (cl = 0; cl < coneSize; ++cl) cells[coff++] = cone[cl] - vStart; 94 } 95 ierr = PetscCalloc1(numVertices, &l2gv);CHKERRQ(ierr); 96 ierr = DMPlexGetVertexNumbering(udm, &globalVertexNum);CHKERRQ(ierr); 97 ierr = ISGetIndices(globalVertexNum, &gV);CHKERRQ(ierr); 98 for (v = 0, numLocVertices = 0; v < numVertices; ++v) { 99 if (gV[v] >= 0) ++numLocVertices; 100 l2gv[v] = gV[v] < 0 ? -(gV[v]+1) : gV[v]; 101 } 102 ierr = ISRestoreIndices(globalVertexNum, &gV);CHKERRQ(ierr); 103 ierr = DMDestroy(&udm);CHKERRQ(ierr); 104 ierr = DMGetCoordinateDM(dm, &cdm);CHKERRQ(ierr); 105 ierr = DMGetDefaultSection(cdm, &coordSection);CHKERRQ(ierr); 106 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 107 ierr = VecGetArrayRead(coordinates, &coords);CHKERRQ(ierr); 108 for (v = vStart; v < vEnd; ++v) { 109 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 110 x[v-vStart] = PetscRealPart(coords[off+0]); 111 if (dim > 1) y[v-vStart] = PetscRealPart(coords[off+1]); 112 if (dim > 2) z[v-vStart] = PetscRealPart(coords[off+2]); 113 } 114 ierr = VecRestoreArrayRead(coordinates, &coords);CHKERRQ(ierr); 115 /* Get boundary mesh */ 116 ierr = DMLabelCreate(bdName, &bdLabelFull);CHKERRQ(ierr); 117 ierr = DMPlexMarkBoundaryFaces(dm, bdLabelFull);CHKERRQ(ierr); 118 ierr = DMLabelGetStratumIS(bdLabelFull, 1, &bdIS);CHKERRQ(ierr); 119 ierr = DMLabelGetStratumSize(bdLabelFull, 1, &numBdFaces);CHKERRQ(ierr); 120 ierr = ISGetIndices(bdIS, &bdFacesFull);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)) ++bdSize; 128 } 129 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 130 } 131 ierr = PetscMalloc2(bdSize, &bdFaces, numBdFaces, &bdFaceIds);CHKERRQ(ierr); 132 for (f = 0, bdSize = 0; f < numBdFaces; ++f) { 133 PetscInt *closure = NULL; 134 PetscInt closureSize, cl; 135 136 ierr = DMPlexGetTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 137 for (cl = 0; cl < closureSize*2; cl += 2) { 138 if ((closure[cl] >= vStart) && (closure[cl] < vEnd)) bdFaces[bdSize++] = closure[cl] - vStart; 139 } 140 ierr = DMPlexRestoreTransitiveClosure(dm, bdFacesFull[f], PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 141 if (bdLabel) {ierr = DMLabelGetValue(bdLabel, bdFacesFull[f], &bdFaceIds[f]);CHKERRQ(ierr);} 142 else {bdFaceIds[f] = 1;} 143 } 144 ierr = ISDestroy(&bdIS);CHKERRQ(ierr); 145 ierr = DMLabelDestroy(&bdLabelFull);CHKERRQ(ierr); 146 /* Get metric */ 147 ierr = VecGetArrayRead(vertexMetric, &met);CHKERRQ(ierr); 148 for (v = 0; v < (vEnd-vStart)*PetscSqr(dim); ++v) metric[v] = PetscRealPart(met[v]); 149 ierr = VecRestoreArrayRead(vertexMetric, &met);CHKERRQ(ierr); 150 #if 0 151 /* Destroy overlap mesh */ 152 ierr = DMDestroy(&dm);CHKERRQ(ierr); 153 #endif 154 /* Create new mesh */ 155 #ifdef PETSC_HAVE_PRAGMATIC 156 switch (dim) { 157 case 2: 158 pragmatic_2d_mpi_init(&numVertices, &numCells, cells, x, y, l2gv, numLocVertices, comm);break; 159 case 3: 160 pragmatic_3d_mpi_init(&numVertices, &numCells, cells, x, y, z, l2gv, numLocVertices, comm);break; 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 ierr = DMPlexCreateFromCellListParallel(comm, dim, numCellsNew, numVerticesNew, numCornersNew, PETSC_TRUE, cellsNew, dim, coordsNew, NULL, dmNew);CHKERRQ(ierr); 188 /* Read out boundary label */ 189 pragmatic_get_boundaryTags(&bdTags); 190 ierr = DMCreateLabel(*dmNew, bdLabel ? bdLabelName : bdName);CHKERRQ(ierr); 191 ierr = DMGetLabel(*dmNew, bdLabel ? bdLabelName : bdName, &bdLabelNew);CHKERRQ(ierr); 192 ierr = DMPlexGetHeightStratum(*dmNew, 0, &cStart, &cEnd);CHKERRQ(ierr); 193 ierr = DMPlexGetHeightStratum(*dmNew, 1, &fStart, &fEnd);CHKERRQ(ierr); 194 ierr = DMPlexGetDepthStratum(*dmNew, 0, &vStart, &vEnd);CHKERRQ(ierr); 195 for (c = cStart; c < cEnd; ++c) { 196 /* Only for simplicial meshes */ 197 coff = (c-cStart)*(dim+1); 198 for (d = 0; d < dim+1; ++d) { 199 if (bdTags[coff+d]) { 200 const PetscInt *faces; 201 PetscInt numFaces, vertices[3], p; 202 203 /* Mark face opposite to this vertex */ 204 for (p = 0, v = 0; p < dim+1; ++p) if (p != d) {vertices[v] = cellsNew[coff+p] + vStart; ++v;} 205 ierr = DMPlexGetFullJoin(*dmNew, dim, vertices, &numFaces, &faces);CHKERRQ(ierr); 206 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]); 207 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]); 208 ierr = DMLabelSetValue(bdLabelNew, faces[0], bdTags[coff+d]);CHKERRQ(ierr); 209 ierr = DMPlexRestoreJoin(*dmNew, dim, vertices, &numFaces, &faces);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 /*@ 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