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