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