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