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