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