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, 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_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 const double *meshCoords = out.pointlist; 114 int *cells = out.tetrahedronlist; 115 116 ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr); 117 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 118 if (label) {ierr = DMCreateLabel(*dm, labelName);CHKERRQ(ierr); ierr = DMGetLabel(*dm, labelName, &glabel);CHKERRQ(ierr);} 119 /* Set labels */ 120 for (v = 0; v < numVertices; ++v) { 121 if (out.pointmarkerlist[v]) { 122 if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 123 } 124 } 125 if (interpolate) { 126 #if 0 127 PetscInt e; 128 129 /* This check is never actually executed for ctetgen (which never returns edgemarkers) and seems to be broken for 130 * tetgen */ 131 for (e = 0; e < out.numberofedges; e++) { 132 if (out.edgemarkerlist[e]) { 133 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 134 const PetscInt *edges; 135 PetscInt numEdges; 136 137 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 138 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 139 if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 140 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 141 } 142 } 143 #endif 144 for (f = 0; f < out.numberoftrifaces; f++) { 145 if (out.trifacemarkerlist[f]) { 146 const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells}; 147 const PetscInt *faces; 148 PetscInt numFaces; 149 150 ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 151 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 152 if (glabel) {ierr = DMLabelSetValue(glabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);} 153 ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 154 } 155 } 156 } 157 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 PETSC_EXTERN PetscErrorCode DMPlexRefine_Tetgen(DM dm, double *maxVolumes, DM *dmRefined) 163 { 164 MPI_Comm comm; 165 const PetscInt dim = 3; 166 const char *labelName = "marker"; 167 ::tetgenio in; 168 ::tetgenio out; 169 DMLabel label; 170 PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal; 171 PetscMPIInt rank; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 176 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 177 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 178 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 179 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 180 ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr); 181 182 in.numberofpoints = vEnd - vStart; 183 if (in.numberofpoints > 0) { 184 PetscSection coordSection; 185 Vec coordinates; 186 PetscScalar *array; 187 188 in.pointlist = new double[in.numberofpoints*dim]; 189 in.pointmarkerlist = new int[in.numberofpoints]; 190 191 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 192 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 193 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 194 for (v = vStart; v < vEnd; ++v) { 195 const PetscInt idx = v - vStart; 196 PetscInt off, d; 197 198 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 199 for (d = 0; d < dim; ++d) in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 200 if (label) { 201 PetscInt val; 202 203 ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr); 204 in.pointmarkerlist[idx] = (int) val; 205 } 206 } 207 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 208 } 209 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 210 211 in.numberofcorners = 4; 212 in.numberoftetrahedra = cEnd - cStart; 213 in.tetrahedronvolumelist = (double*) maxVolumes; 214 if (in.numberoftetrahedra > 0) { 215 in.tetrahedronlist = new int[in.numberoftetrahedra*in.numberofcorners]; 216 for (c = cStart; c < cEnd; ++c) { 217 const PetscInt idx = c - cStart; 218 PetscInt *closure = NULL; 219 PetscInt closureSize; 220 221 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 222 if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize); 223 for (v = 0; v < 4; ++v) { 224 in.tetrahedronlist[idx*in.numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart; 225 } 226 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 227 } 228 } 229 /* TODO: Put in boundary faces with markers */ 230 if (!rank) { 231 char args[32]; 232 233 #if 1 234 /* Take away 'Q' for verbose output */ 235 ierr = PetscStrcpy(args, "qezQra");CHKERRQ(ierr); 236 #else 237 ierr = PetscStrcpy(args, "qezraVVVV");CHKERRQ(ierr); 238 #endif 239 ::tetrahedralize(args, &in, &out); 240 } 241 in.tetrahedronvolumelist = NULL; 242 243 { 244 DMLabel rlabel = NULL; 245 const PetscInt numCorners = 4; 246 const PetscInt numCells = out.numberoftetrahedra; 247 const PetscInt numVertices = out.numberofpoints; 248 const double *meshCoords = out.pointlist; 249 int *cells = out.tetrahedronlist; 250 251 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 252 253 ierr = DMPlexInvertCells_Tetgen(numCells, numCorners, cells);CHKERRQ(ierr); 254 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 255 if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);} 256 /* Set labels */ 257 for (v = 0; v < numVertices; ++v) { 258 if (out.pointmarkerlist[v]) { 259 if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 260 } 261 } 262 if (interpolate) { 263 PetscInt f; 264 #if 0 265 PetscInt e; 266 267 for (e = 0; e < out.numberofedges; e++) { 268 if (out.edgemarkerlist[e]) { 269 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 270 const PetscInt *edges; 271 PetscInt numEdges; 272 273 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 274 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 275 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 276 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 277 } 278 } 279 #endif 280 for (f = 0; f < out.numberoftrifaces; f++) { 281 if (out.trifacemarkerlist[f]) { 282 const PetscInt vertices[3] = {out.trifacelist[f*3+0]+numCells, out.trifacelist[f*3+1]+numCells, out.trifacelist[f*3+2]+numCells}; 283 const PetscInt *faces; 284 PetscInt numFaces; 285 286 ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 287 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 288 if (rlabel) {ierr = DMLabelSetValue(rlabel, faces[0], out.trifacemarkerlist[f]);CHKERRQ(ierr);} 289 ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 290 } 291 } 292 } 293 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 294 } 295 PetscFunctionReturn(0); 296 } 297