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