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, int cells[]) 7 { 8 PetscInt bound = numCells*numCorners, coff; 9 10 PetscFunctionBegin; 11 #define SWAP(a,b) do { int 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 double *meshCoords; 118 int *cells = out->tetrahedronlist; 119 120 if (sizeof (PetscReal) == sizeof (double)) { 121 meshCoords = (double *) out->pointlist; 122 } 123 else { 124 PetscInt i; 125 126 ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr); 127 for (i = 0; i < 3 * numVertices; i++) { 128 meshCoords[i] = (PetscReal) out->pointlist[i]; 129 } 130 } 131 132 ierr = DMPlexInvertCells_CTetgen(numCells, numCorners, cells);CHKERRQ(ierr); 133 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 134 if (sizeof (PetscReal) != sizeof (double)) { 135 ierr = PetscFree(meshCoords);CHKERRQ(ierr); 136 } 137 if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);} 138 /* Set labels */ 139 for (v = 0; v < numVertices; ++v) { 140 if (out->pointmarkerlist[v]) { 141 if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);} 142 } 143 } 144 if (interpolate) { 145 PetscInt e; 146 147 for (e = 0; e < out->numberofedges; e++) { 148 if (out->edgemarkerlist[e]) { 149 const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells}; 150 const PetscInt *edges; 151 PetscInt numEdges; 152 153 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 154 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 155 if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);} 156 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 157 } 158 } 159 for (f = 0; f < out->numberoftrifaces; f++) { 160 if (out->trifacemarkerlist[f]) { 161 const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells}; 162 const PetscInt *faces; 163 PetscInt numFaces; 164 165 ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 166 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 167 if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);} 168 ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 169 } 170 } 171 } 172 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 173 } 174 175 ierr = PLCDestroy(&in);CHKERRQ(ierr); 176 ierr = PLCDestroy(&out);CHKERRQ(ierr); 177 PetscFunctionReturn(0); 178 } 179 180 PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined) 181 { 182 MPI_Comm comm; 183 const PetscInt dim = 3; 184 const char *labelName = "marker"; 185 PLC *in, *out; 186 DMLabel label; 187 PetscInt verbose = 0, vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal; 188 PetscMPIInt rank; 189 PetscErrorCode ierr; 190 191 PetscFunctionBegin; 192 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 193 ierr = PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr); 194 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 195 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 196 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 197 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 198 ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr); 199 ierr = PLCCreate(&in);CHKERRQ(ierr); 200 ierr = PLCCreate(&out);CHKERRQ(ierr); 201 202 in->numberofpoints = vEnd - vStart; 203 if (in->numberofpoints > 0) { 204 PetscSection coordSection; 205 Vec coordinates; 206 PetscScalar *array; 207 208 ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr); 209 ierr = PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);CHKERRQ(ierr); 210 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 211 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 212 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 213 for (v = vStart; v < vEnd; ++v) { 214 const PetscInt idx = v - vStart; 215 PetscInt off, d, m = -1; 216 217 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 218 for (d = 0; d < dim; ++d) { 219 in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 220 } 221 if (label) {ierr = DMLabelGetValue(label, v, &m);CHKERRQ(ierr);} 222 223 in->pointmarkerlist[idx] = (int) m; 224 } 225 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 226 } 227 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 228 229 in->numberofcorners = 4; 230 in->numberoftetrahedra = cEnd - cStart; 231 in->tetrahedronvolumelist = maxVolumes; 232 if (in->numberoftetrahedra > 0) { 233 ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr); 234 for (c = cStart; c < cEnd; ++c) { 235 const PetscInt idx = c - cStart; 236 PetscInt *closure = NULL; 237 PetscInt closureSize; 238 239 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 240 if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize); 241 for (v = 0; v < 4; ++v) { 242 in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart; 243 } 244 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 245 } 246 } 247 if (!rank) { 248 TetGenOpts t; 249 250 ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr); 251 252 t.in = dm; /* Should go away */ 253 t.refine = 1; 254 t.varvolume = 1; 255 t.quality = 1; 256 t.edgesout = 1; 257 t.zeroindex = 1; 258 t.quiet = 1; 259 t.verbose = verbose; /* Change this */ 260 261 ierr = TetGenCheckOpts(&t);CHKERRQ(ierr); 262 ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr); 263 } 264 in->tetrahedronvolumelist = NULL; 265 { 266 DMLabel rlabel = NULL; 267 const PetscInt numCorners = 4; 268 const PetscInt numCells = out->numberoftetrahedra; 269 const PetscInt numVertices = out->numberofpoints; 270 double *meshCoords; 271 int *cells = out->tetrahedronlist; 272 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 273 274 if (sizeof (PetscReal) == sizeof (double)) { 275 meshCoords = (double *) out->pointlist; 276 } 277 else { 278 PetscInt i; 279 280 ierr = PetscMalloc1(3 * numVertices,&meshCoords);CHKERRQ(ierr); 281 for (i = 0; i < 3 * numVertices; i++) { 282 meshCoords[i] = (PetscReal) out->pointlist[i]; 283 } 284 } 285 286 ierr = DMPlexInvertCells_CTetgen(numCells, numCorners, cells);CHKERRQ(ierr); 287 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 288 if (sizeof (PetscReal) != sizeof (double)) { 289 ierr = PetscFree(meshCoords);CHKERRQ(ierr); 290 } 291 if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);} 292 /* Set labels */ 293 for (v = 0; v < numVertices; ++v) { 294 if (out->pointmarkerlist[v]) { 295 if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr);} 296 } 297 } 298 if (interpolate) { 299 PetscInt e, f; 300 301 for (e = 0; e < out->numberofedges; e++) { 302 if (out->edgemarkerlist[e]) { 303 const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells}; 304 const PetscInt *edges; 305 PetscInt numEdges; 306 307 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 308 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 309 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr);} 310 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 311 } 312 } 313 for (f = 0; f < out->numberoftrifaces; f++) { 314 if (out->trifacemarkerlist[f]) { 315 const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells}; 316 const PetscInt *faces; 317 PetscInt numFaces; 318 319 ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 320 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 321 if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr);} 322 ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 323 } 324 } 325 } 326 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 327 } 328 ierr = PLCDestroy(&in);CHKERRQ(ierr); 329 ierr = PLCDestroy(&out);CHKERRQ(ierr); 330 PetscFunctionReturn(0); 331 } 332