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