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