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