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