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