1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #include <triangle.h> 4 5 static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx) 6 { 7 PetscFunctionBegin; 8 inputCtx->numberofpoints = 0; 9 inputCtx->numberofpointattributes = 0; 10 inputCtx->pointlist = NULL; 11 inputCtx->pointattributelist = NULL; 12 inputCtx->pointmarkerlist = NULL; 13 inputCtx->numberofsegments = 0; 14 inputCtx->segmentlist = NULL; 15 inputCtx->segmentmarkerlist = NULL; 16 inputCtx->numberoftriangleattributes = 0; 17 inputCtx->trianglelist = NULL; 18 inputCtx->numberofholes = 0; 19 inputCtx->holelist = NULL; 20 inputCtx->numberofregions = 0; 21 inputCtx->regionlist = NULL; 22 PetscFunctionReturn(0); 23 } 24 25 static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx) 26 { 27 PetscFunctionBegin; 28 outputCtx->numberofpoints = 0; 29 outputCtx->pointlist = NULL; 30 outputCtx->pointattributelist = NULL; 31 outputCtx->pointmarkerlist = NULL; 32 outputCtx->numberoftriangles = 0; 33 outputCtx->trianglelist = NULL; 34 outputCtx->triangleattributelist = NULL; 35 outputCtx->neighborlist = NULL; 36 outputCtx->segmentlist = NULL; 37 outputCtx->segmentmarkerlist = NULL; 38 outputCtx->numberofedges = 0; 39 outputCtx->edgelist = NULL; 40 outputCtx->edgemarkerlist = NULL; 41 PetscFunctionReturn(0); 42 } 43 44 static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx) 45 { 46 PetscFunctionBegin; 47 free(outputCtx->pointlist); 48 free(outputCtx->pointmarkerlist); 49 free(outputCtx->segmentlist); 50 free(outputCtx->segmentmarkerlist); 51 free(outputCtx->edgelist); 52 free(outputCtx->edgemarkerlist); 53 free(outputCtx->trianglelist); 54 free(outputCtx->neighborlist); 55 PetscFunctionReturn(0); 56 } 57 58 PETSC_EXTERN PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm) 59 { 60 MPI_Comm comm; 61 DM_Plex *mesh = (DM_Plex *) boundary->data; 62 PetscInt dim = 2; 63 const PetscBool createConvexHull = PETSC_FALSE; 64 const PetscBool constrained = PETSC_FALSE; 65 const char *labelName = "marker"; 66 const char *labelName2 = "Face Sets"; 67 struct triangulateio in; 68 struct triangulateio out; 69 DMLabel label, label2; 70 PetscInt vStart, vEnd, v, eStart, eEnd, e; 71 PetscMPIInt rank; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr); 76 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 77 ierr = InitInput_Triangle(&in);CHKERRQ(ierr); 78 ierr = InitOutput_Triangle(&out);CHKERRQ(ierr); 79 ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr); 80 ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr); 81 ierr = DMGetLabel(boundary, labelName2, &label2);CHKERRQ(ierr); 82 83 in.numberofpoints = vEnd - vStart; 84 if (in.numberofpoints > 0) { 85 PetscSection coordSection; 86 Vec coordinates; 87 PetscScalar *array; 88 89 ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr); 90 ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr); 91 ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr); 92 ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr); 93 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 94 for (v = vStart; v < vEnd; ++v) { 95 const PetscInt idx = v - vStart; 96 PetscInt off, d; 97 98 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 99 for (d = 0; d < dim; ++d) { 100 in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 101 } 102 if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);} 103 } 104 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 105 } 106 ierr = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr); 107 in.numberofsegments = eEnd - eStart; 108 if (in.numberofsegments > 0) { 109 ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr); 110 ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr); 111 for (e = eStart; e < eEnd; ++e) { 112 const PetscInt idx = e - eStart; 113 const PetscInt *cone; 114 115 ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr); 116 117 in.segmentlist[idx*2+0] = cone[0] - vStart; 118 in.segmentlist[idx*2+1] = cone[1] - vStart; 119 120 if (label) {ierr = DMLabelGetValue(label, e, &in.segmentmarkerlist[idx]);CHKERRQ(ierr);} 121 } 122 } 123 #if 0 /* Do not currently support holes */ 124 PetscReal *holeCoords; 125 PetscInt h, d; 126 127 ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr); 128 if (in.numberofholes > 0) { 129 ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr); 130 for (h = 0; h < in.numberofholes; ++h) { 131 for (d = 0; d < dim; ++d) { 132 in.holelist[h*dim+d] = holeCoords[h*dim+d]; 133 } 134 } 135 } 136 #endif 137 if (!rank) { 138 char args[32]; 139 140 /* Take away 'Q' for verbose output */ 141 ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr); 142 if (createConvexHull) {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);} 143 if (constrained) {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);} 144 if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);} 145 else {triangulate(args, &in, &out, NULL);} 146 } 147 ierr = PetscFree(in.pointlist);CHKERRQ(ierr); 148 ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr); 149 ierr = PetscFree(in.segmentlist);CHKERRQ(ierr); 150 ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr); 151 ierr = PetscFree(in.holelist);CHKERRQ(ierr); 152 153 { 154 DMLabel glabel = NULL; 155 DMLabel glabel2 = NULL; 156 const PetscInt numCorners = 3; 157 const PetscInt numCells = out.numberoftriangles; 158 const PetscInt numVertices = out.numberofpoints; 159 const int *cells = out.trianglelist; 160 const double *meshCoords = out.pointlist; 161 162 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 163 if (label) {ierr = DMCreateLabel(*dm, labelName); ierr = DMGetLabel(*dm, labelName, &glabel);} 164 if (label2) {ierr = DMCreateLabel(*dm, labelName2); ierr = DMGetLabel(*dm, labelName2, &glabel2);} 165 /* Set labels */ 166 for (v = 0; v < numVertices; ++v) { 167 if (out.pointmarkerlist[v]) { 168 if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 169 } 170 } 171 if (interpolate) { 172 for (e = 0; e < out.numberofedges; e++) { 173 if (out.edgemarkerlist[e]) { 174 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 175 const PetscInt *edges; 176 PetscInt numEdges; 177 178 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 179 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 180 if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 181 if (glabel2) {ierr = DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 182 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 183 } 184 } 185 } 186 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 187 } 188 #if 0 /* Do not currently support holes */ 189 ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr); 190 #endif 191 ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, double *maxVolumes, DM *dmRefined) 196 { 197 MPI_Comm comm; 198 PetscInt dim = 2; 199 const char *labelName = "marker"; 200 struct triangulateio in; 201 struct triangulateio out; 202 DMLabel label; 203 PetscInt vStart, vEnd, v, cStart, cEnd, c, depth, depthGlobal; 204 PetscMPIInt rank; 205 PetscErrorCode ierr; 206 207 PetscFunctionBegin; 208 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 209 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 210 ierr = InitInput_Triangle(&in);CHKERRQ(ierr); 211 ierr = InitOutput_Triangle(&out);CHKERRQ(ierr); 212 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 213 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 214 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 215 ierr = DMGetLabel(dm, labelName, &label);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; 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, &in.pointmarkerlist[idx]);CHKERRQ(ierr);} 237 } 238 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 239 } 240 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 241 242 in.numberofcorners = 3; 243 in.numberoftriangles = cEnd - cStart; 244 245 in.trianglearealist = (double*) maxVolumes; 246 if (in.numberoftriangles > 0) { 247 ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr); 248 for (c = cStart; c < cEnd; ++c) { 249 const PetscInt idx = c - cStart; 250 PetscInt *closure = NULL; 251 PetscInt closureSize; 252 253 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 254 if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize); 255 for (v = 0; v < 3; ++v) { 256 in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart; 257 } 258 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 259 } 260 } 261 /* TODO: Segment markers are missing on input */ 262 #if 0 /* Do not currently support holes */ 263 PetscReal *holeCoords; 264 PetscInt h, d; 265 266 ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr); 267 if (in.numberofholes > 0) { 268 ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr); 269 for (h = 0; h < in.numberofholes; ++h) { 270 for (d = 0; d < dim; ++d) { 271 in.holelist[h*dim+d] = holeCoords[h*dim+d]; 272 } 273 } 274 } 275 #endif 276 if (!rank) { 277 char args[32]; 278 279 /* Take away 'Q' for verbose output */ 280 ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr); 281 triangulate(args, &in, &out, NULL); 282 } 283 ierr = PetscFree(in.pointlist);CHKERRQ(ierr); 284 ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr); 285 ierr = PetscFree(in.segmentlist);CHKERRQ(ierr); 286 ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr); 287 ierr = PetscFree(in.trianglelist);CHKERRQ(ierr); 288 289 { 290 DMLabel rlabel = NULL; 291 const PetscInt numCorners = 3; 292 const PetscInt numCells = out.numberoftriangles; 293 const PetscInt numVertices = out.numberofpoints; 294 const int *cells = out.trianglelist; 295 const double *meshCoords = out.pointlist; 296 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 297 298 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 299 if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);} 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 e; 308 309 for (e = 0; e < out.numberofedges; e++) { 310 if (out.edgemarkerlist[e]) { 311 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 312 const PetscInt *edges; 313 PetscInt numEdges; 314 315 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 316 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 317 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 318 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 319 } 320 } 321 } 322 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 323 } 324 #if 0 /* Do not currently support holes */ 325 ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr); 326 #endif 327 ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330