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