1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <tetgen.h> 4 5 /* This is to fix the tetrahedron orientation from TetGen */ 6 static PetscErrorCode DMPlexInvertCells_Internal(PetscInt dim, PetscInt numCells, PetscInt numCorners, int cells[]) 7 { 8 PetscInt bound = numCells*numCorners, coff; 9 PetscErrorCode ierr; 10 11 PetscFunctionBegin; 12 for (coff = 0; coff < bound; coff += numCorners) { 13 ierr = DMPlexInvertCell(dim, numCorners, &cells[coff]);CHKERRQ(ierr); 14 } 15 PetscFunctionReturn(0); 16 } 17 18 PETSC_EXTERN PetscErrorCode DMPlexGenerate_Tetgen(DM boundary, PetscBool interpolate, DM *dm) 19 { 20 MPI_Comm comm; 21 DM_Plex *mesh = (DM_Plex *) boundary->data; 22 const PetscInt dim = 3; 23 const char *labelName = "marker"; 24 ::tetgenio in; 25 ::tetgenio out; 26 DMLabel label; 27 PetscInt vStart, vEnd, v, fStart, fEnd, f; 28 PetscMPIInt rank; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr); 33 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 34 ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr); 35 ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr); 36 37 in.numberofpoints = vEnd - vStart; 38 if (in.numberofpoints > 0) { 39 PetscSection coordSection; 40 Vec coordinates; 41 PetscScalar *array; 42 43 in.pointlist = new double[in.numberofpoints*dim]; 44 in.pointmarkerlist = new int[in.numberofpoints]; 45 46 ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr); 47 ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr); 48 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 49 for (v = vStart; v < vEnd; ++v) { 50 const PetscInt idx = v - vStart; 51 PetscInt off, d; 52 53 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 54 for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 55 if (label) { 56 PetscInt val; 57 58 ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr); 59 in.pointmarkerlist[idx] = (int) val; 60 } 61 } 62 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 63 } 64 ierr = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr); 65 66 in.numberoffacets = fEnd - fStart; 67 if (in.numberoffacets > 0) { 68 in.facetlist = new tetgenio::facet[in.numberoffacets]; 69 in.facetmarkerlist = new int[in.numberoffacets]; 70 for (f = fStart; f < fEnd; ++f) { 71 const PetscInt idx = f - fStart; 72 PetscInt *points = NULL, numPoints, p, numVertices = 0, v; 73 74 in.facetlist[idx].numberofpolygons = 1; 75 in.facetlist[idx].polygonlist = new tetgenio::polygon[in.facetlist[idx].numberofpolygons]; 76 in.facetlist[idx].numberofholes = 0; 77 in.facetlist[idx].holelist = NULL; 78 79 ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 80 for (p = 0; p < numPoints*2; p += 2) { 81 const PetscInt point = points[p]; 82 if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point; 83 } 84 85 tetgenio::polygon *poly = in.facetlist[idx].polygonlist; 86 poly->numberofvertices = numVertices; 87 poly->vertexlist = new int[poly->numberofvertices]; 88 for (v = 0; v < numVertices; ++v) { 89 const PetscInt vIdx = points[v] - vStart; 90 poly->vertexlist[v] = vIdx; 91 } 92 if (label) { 93 PetscInt val; 94 95 ierr = DMLabelGetValue(label, f, &val);CHKERRQ(ierr); 96 in.facetmarkerlist[idx] = (int) val; 97 } 98 ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 99 } 100 } 101 if (!rank) { 102 char args[32]; 103 104 /* Take away 'Q' for verbose output */ 105 ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr); 106 if (mesh->tetgenOpts) {::tetrahedralize(mesh->tetgenOpts, &in, &out);} 107 else {::tetrahedralize(args, &in, &out);} 108 } 109 { 110 DMLabel glabel = NULL; 111 const PetscInt numCorners = 4; 112 const PetscInt numCells = out.numberoftetrahedra; 113 const PetscInt numVertices = out.numberofpoints; 114 const double *meshCoords = out.pointlist; 115 int *cells = out.tetrahedronlist; 116 117 ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr); 118 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 119 if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);} 120 /* Set labels */ 121 for (v = 0; v < numVertices; ++v) { 122 if (out.pointmarkerlist[v]) { 123 if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 124 } 125 } 126 if (interpolate) { 127 #if 0 128 PetscInt e; 129 130 /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for 131 * tetgen */ 132 for (e = 0; e < out.numberofedges; e++) { 133 if (out.edgemarkerlist[e]) { 134 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 135 const PetscInt *edges; 136 PetscInt numEdges; 137 138 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 139 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 140 if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 141 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 142 } 143 } 144 #endif 145 for (f = 0; f < out.numberoftrifaces; f++) { 146 if (out.trifacemarkerlist[f]) { 147 const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells}; 148 const PetscInt *faces; 149 PetscInt numFaces; 150 151 ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 152 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 153 if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);} 154 ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 155 } 156 } 157 } 158 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 159 } 160 PetscFunctionReturn(0); 161 } 162 163 PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined) 164 { 165 MPI_Comm comm; 166 const PetscInt dim = 3; 167 const char *labelName = "marker"; 168 ::tetgenio in; 169 ::tetgenio out; 170 DMLabel label; 171 PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal; 172 PetscMPIInt rank; 173 PetscErrorCode ierr; 174 175 PetscFunctionBegin; 176 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 177 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 178 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 179 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 180 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 181 ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr); 182 183 in.numberofpoints = vEnd - vStart; 184 if (in.numberofpoints > 0) { 185 PetscSection coordSection; 186 Vec coordinates; 187 PetscScalar *array; 188 189 in.pointlist = new double[in.numberofpoints*dim]; 190 in.pointmarkerlist = new int[in.numberofpoints]; 191 192 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 193 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 194 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 195 for (v = vStart; v < vEnd; ++v) { 196 const PetscInt idx = v - vStart; 197 PetscInt off, d; 198 199 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 200 for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 201 if (label) { 202 PetscInt val; 203 204 ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr); 205 in.pointmarkerlist[idx] = (int) val; 206 } 207 } 208 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 209 } 210 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 211 212 in.numberofcorners = 4; 213 in.numberoftetrahedra = cEnd - cStart; 214 in.tetrahedronvolumelist = (double*) maxVolumes; 215 if (in.numberoftetrahedra > 0) { 216 in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners]; 217 for (c = cStart; c < cEnd; ++c) { 218 const PetscInt idx = c - cStart; 219 PetscInt *closure = NULL; 220 PetscInt closureSize; 221 222 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 223 if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize); 224 for (v = 0; v < 4; ++v) { 225 in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart; 226 } 227 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 228 } 229 } 230 /* TODO: Put in boundary faces with markers */ 231 if (!rank) { 232 char args[32]; 233 234 #if 1 235 /* Take away 'Q' for verbose output */ 236 ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr); 237 #else 238 ierr = PetscStrcpy(args, "qezraVVVV");CHKERRQ(ierr); 239 #endif 240 ::tetrahedralize(args, &in, &out); 241 } 242 in.tetrahedronvolumelist = NULL; 243 244 { 245 DMLabel rlabel = NULL; 246 const PetscInt numCorners = 4; 247 const PetscInt numCells = out.numberoftetrahedra; 248 const PetscInt numVertices = out.numberofpoints; 249 const double *meshCoords = out.pointlist; 250 int *cells = out.tetrahedronlist; 251 252 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 253 254 ierr = DMPlexInvertCells_Internal(dim, numCells, numCorners, cells);CHKERRQ(ierr); 255 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 256 if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);} 257 /* Set labels */ 258 for (v = 0; v < numVertices; ++v) { 259 if (out.pointmarkerlist[v]) { 260 if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 261 } 262 } 263 if (interpolate) { 264 PetscInt f; 265 #if 0 266 PetscInt e; 267 268 for (e = 0; e < out.numberofedges; e++) { 269 if (out.edgemarkerlist[e]) { 270 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 271 const PetscInt *edges; 272 PetscInt numEdges; 273 274 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 275 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 276 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 277 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 278 } 279 } 280 #endif 281 for (f = 0; f < out.numberoftrifaces; f++) { 282 if (out.trifacemarkerlist[f]) { 283 const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells}; 284 const PetscInt *faces; 285 PetscInt numFaces; 286 287 ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 288 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 289 if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);} 290 ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 291 } 292 } 293 } 294 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 295 } 296 PetscFunctionReturn(0); 297 } 298