1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #if !defined(ANSI_DECLARATORS) 4 #define ANSI_DECLARATORS 5 #endif 6 #include <triangle.h> 7 8 static PetscErrorCode InitInput_Triangle(struct triangulateio *inputCtx) 9 { 10 PetscFunctionBegin; 11 inputCtx->numberofpoints = 0; 12 inputCtx->numberofpointattributes = 0; 13 inputCtx->pointlist = NULL; 14 inputCtx->pointattributelist = NULL; 15 inputCtx->pointmarkerlist = NULL; 16 inputCtx->numberofsegments = 0; 17 inputCtx->segmentlist = NULL; 18 inputCtx->segmentmarkerlist = NULL; 19 inputCtx->numberoftriangleattributes = 0; 20 inputCtx->trianglelist = NULL; 21 inputCtx->numberofholes = 0; 22 inputCtx->holelist = NULL; 23 inputCtx->numberofregions = 0; 24 inputCtx->regionlist = NULL; 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode InitOutput_Triangle(struct triangulateio *outputCtx) 29 { 30 PetscFunctionBegin; 31 outputCtx->numberofpoints = 0; 32 outputCtx->pointlist = NULL; 33 outputCtx->pointattributelist = NULL; 34 outputCtx->pointmarkerlist = NULL; 35 outputCtx->numberoftriangles = 0; 36 outputCtx->trianglelist = NULL; 37 outputCtx->triangleattributelist = NULL; 38 outputCtx->neighborlist = NULL; 39 outputCtx->segmentlist = NULL; 40 outputCtx->segmentmarkerlist = NULL; 41 outputCtx->numberofedges = 0; 42 outputCtx->edgelist = NULL; 43 outputCtx->edgemarkerlist = NULL; 44 PetscFunctionReturn(0); 45 } 46 47 static PetscErrorCode FiniOutput_Triangle(struct triangulateio *outputCtx) 48 { 49 PetscFunctionBegin; 50 free(outputCtx->pointlist); 51 free(outputCtx->pointmarkerlist); 52 free(outputCtx->segmentlist); 53 free(outputCtx->segmentmarkerlist); 54 free(outputCtx->edgelist); 55 free(outputCtx->edgemarkerlist); 56 free(outputCtx->trianglelist); 57 free(outputCtx->neighborlist); 58 PetscFunctionReturn(0); 59 } 60 61 PETSC_EXTERN PetscErrorCode DMPlexGenerate_Triangle(DM boundary, PetscBool interpolate, DM *dm) 62 { 63 MPI_Comm comm; 64 DM_Plex *mesh = (DM_Plex *) boundary->data; 65 PetscInt dim = 2; 66 const PetscBool createConvexHull = PETSC_FALSE; 67 const PetscBool constrained = PETSC_FALSE; 68 const char *labelName = "marker"; 69 const char *labelName2 = "Face Sets"; 70 struct triangulateio in; 71 struct triangulateio out; 72 DMLabel label, label2; 73 PetscInt vStart, vEnd, v, eStart, eEnd, e; 74 PetscMPIInt rank; 75 PetscErrorCode ierr; 76 77 PetscFunctionBegin; 78 ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr); 79 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 80 ierr = InitInput_Triangle(&in);CHKERRQ(ierr); 81 ierr = InitOutput_Triangle(&out);CHKERRQ(ierr); 82 ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr); 83 ierr = DMGetLabel(boundary, labelName, &label);CHKERRQ(ierr); 84 ierr = DMGetLabel(boundary, labelName2, &label2);CHKERRQ(ierr); 85 86 in.numberofpoints = vEnd - vStart; 87 if (in.numberofpoints > 0) { 88 PetscSection coordSection; 89 Vec coordinates; 90 PetscScalar *array; 91 92 ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr); 93 ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr); 94 ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr); 95 ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr); 96 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 97 for (v = vStart; v < vEnd; ++v) { 98 const PetscInt idx = v - vStart; 99 PetscInt val, off, d; 100 101 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 102 for (d = 0; d < dim; ++d) { 103 in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 104 } 105 if (label) { 106 ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr); 107 in.pointmarkerlist[idx] = val; 108 } 109 } 110 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 111 } 112 ierr = DMPlexGetHeightStratum(boundary, 0, &eStart, &eEnd);CHKERRQ(ierr); 113 in.numberofsegments = eEnd - eStart; 114 if (in.numberofsegments > 0) { 115 ierr = PetscMalloc1(in.numberofsegments*2, &in.segmentlist);CHKERRQ(ierr); 116 ierr = PetscMalloc1(in.numberofsegments, &in.segmentmarkerlist);CHKERRQ(ierr); 117 for (e = eStart; e < eEnd; ++e) { 118 const PetscInt idx = e - eStart; 119 const PetscInt *cone; 120 PetscInt val; 121 122 ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr); 123 124 in.segmentlist[idx*2+0] = cone[0] - vStart; 125 in.segmentlist[idx*2+1] = cone[1] - vStart; 126 127 if (label) { 128 ierr = DMLabelGetValue(label, e, &val);CHKERRQ(ierr); 129 in.segmentmarkerlist[idx] = val; 130 } 131 } 132 } 133 #if 0 /* Do not currently support holes */ 134 PetscReal *holeCoords; 135 PetscInt h, d; 136 137 ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr); 138 if (in.numberofholes > 0) { 139 ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr); 140 for (h = 0; h < in.numberofholes; ++h) { 141 for (d = 0; d < dim; ++d) { 142 in.holelist[h*dim+d] = holeCoords[h*dim+d]; 143 } 144 } 145 } 146 #endif 147 if (!rank) { 148 char args[32]; 149 150 /* Take away 'Q' for verbose output */ 151 ierr = PetscStrcpy(args, "pqezQ");CHKERRQ(ierr); 152 if (createConvexHull) {ierr = PetscStrcat(args, "c");CHKERRQ(ierr);} 153 if (constrained) {ierr = PetscStrcpy(args, "zepDQ");CHKERRQ(ierr);} 154 if (mesh->triangleOpts) {triangulate(mesh->triangleOpts, &in, &out, NULL);} 155 else {triangulate(args, &in, &out, NULL);} 156 } 157 ierr = PetscFree(in.pointlist);CHKERRQ(ierr); 158 ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr); 159 ierr = PetscFree(in.segmentlist);CHKERRQ(ierr); 160 ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr); 161 ierr = PetscFree(in.holelist);CHKERRQ(ierr); 162 163 { 164 DMLabel glabel = NULL; 165 DMLabel glabel2 = NULL; 166 const PetscInt numCorners = 3; 167 const PetscInt numCells = out.numberoftriangles; 168 const PetscInt numVertices = out.numberofpoints; 169 PetscInt *cells; 170 PetscReal *meshCoords; 171 172 if (sizeof (PetscReal) == sizeof (out.pointlist[0])) { 173 meshCoords = (PetscReal *) out.pointlist; 174 } else { 175 PetscInt i; 176 177 ierr = PetscMalloc1(dim * numVertices,&meshCoords);CHKERRQ(ierr); 178 for (i = 0; i < dim * numVertices; i++) { 179 meshCoords[i] = (PetscReal) out.pointlist[i]; 180 } 181 } 182 if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) { 183 cells = (PetscInt *) out.trianglelist; 184 } else { 185 PetscInt i; 186 187 ierr = PetscMalloc1(numCells * numCorners, &cells);CHKERRQ(ierr); 188 for (i = 0; i < numCells * numCorners; i++) { 189 cells[i] = (PetscInt) out.trianglelist[i]; 190 } 191 } 192 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 193 if (sizeof (PetscReal) != sizeof (out.pointlist[0])) { 194 ierr = PetscFree(meshCoords);CHKERRQ(ierr); 195 } 196 if (sizeof (PetscInt) != sizeof (out.trianglelist[0])) { 197 ierr = PetscFree(cells);CHKERRQ(ierr); 198 } 199 if (label) { 200 ierr = DMCreateLabel(*dm, labelName);CHKERRQ(ierr); 201 ierr = DMGetLabel(*dm, labelName, &glabel);CHKERRQ(ierr); 202 } 203 if (label2) { 204 ierr = DMCreateLabel(*dm, labelName2);CHKERRQ(ierr); 205 ierr = DMGetLabel(*dm, labelName2, &glabel2);CHKERRQ(ierr); 206 } 207 /* Set labels */ 208 for (v = 0; v < numVertices; ++v) { 209 if (out.pointmarkerlist[v]) { 210 if (glabel) {ierr = DMLabelSetValue(glabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 211 } 212 } 213 if (interpolate) { 214 for (e = 0; e < out.numberofedges; e++) { 215 if (out.edgemarkerlist[e]) { 216 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 217 const PetscInt *edges; 218 PetscInt numEdges; 219 220 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 221 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 222 if (glabel) {ierr = DMLabelSetValue(glabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 223 if (glabel2) {ierr = DMLabelSetValue(glabel2, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 224 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 225 } 226 } 227 } 228 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 229 } 230 #if 0 /* Do not currently support holes */ 231 ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr); 232 #endif 233 ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr); 234 PetscFunctionReturn(0); 235 } 236 237 PETSC_EXTERN PetscErrorCode DMPlexRefine_Triangle(DM dm, PetscReal *inmaxVolumes, DM *dmRefined) 238 { 239 MPI_Comm comm; 240 PetscInt dim = 2; 241 const char *labelName = "marker"; 242 struct triangulateio in; 243 struct triangulateio out; 244 DMLabel label; 245 PetscInt vStart, vEnd, v, gcStart, cStart, cEnd, c, depth, depthGlobal; 246 PetscMPIInt rank; 247 PetscErrorCode ierr; 248 double *maxVolumes; 249 250 PetscFunctionBegin; 251 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 252 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 253 ierr = InitInput_Triangle(&in);CHKERRQ(ierr); 254 ierr = InitOutput_Triangle(&out);CHKERRQ(ierr); 255 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 256 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRMPI(ierr); 257 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 258 ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr); 259 260 in.numberofpoints = vEnd - vStart; 261 if (in.numberofpoints > 0) { 262 PetscSection coordSection; 263 Vec coordinates; 264 PetscScalar *array; 265 266 ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr); 267 ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr); 268 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 269 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 270 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 271 for (v = vStart; v < vEnd; ++v) { 272 const PetscInt idx = v - vStart; 273 PetscInt off, d, val; 274 275 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 276 for (d = 0; d < dim; ++d) { 277 in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 278 } 279 if (label) { 280 ierr = DMLabelGetValue(label, v, &val);CHKERRQ(ierr); 281 in.pointmarkerlist[idx] = val; 282 } 283 } 284 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 285 } 286 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 287 ierr = DMPlexGetGhostCellStratum(dm, &gcStart, NULL);CHKERRQ(ierr); 288 if (gcStart >= 0) cEnd = gcStart; 289 290 in.numberofcorners = 3; 291 in.numberoftriangles = cEnd - cStart; 292 293 #if !defined(PETSC_USE_REAL_DOUBLE) 294 ierr = PetscMalloc1(cEnd - cStart,&maxVolumes);CHKERRQ(ierr); 295 for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c]; 296 #else 297 maxVolumes = inmaxVolumes; 298 #endif 299 300 in.trianglearealist = (double*) maxVolumes; 301 if (in.numberoftriangles > 0) { 302 ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr); 303 for (c = cStart; c < cEnd; ++c) { 304 const PetscInt idx = c - cStart; 305 PetscInt *closure = NULL; 306 PetscInt closureSize; 307 308 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 309 if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize); 310 for (v = 0; v < 3; ++v) { 311 in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart; 312 } 313 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 314 } 315 } 316 /* TODO: Segment markers are missing on input */ 317 #if 0 /* Do not currently support holes */ 318 PetscReal *holeCoords; 319 PetscInt h, d; 320 321 ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr); 322 if (in.numberofholes > 0) { 323 ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr); 324 for (h = 0; h < in.numberofholes; ++h) { 325 for (d = 0; d < dim; ++d) { 326 in.holelist[h*dim+d] = holeCoords[h*dim+d]; 327 } 328 } 329 } 330 #endif 331 if (!rank) { 332 char args[32]; 333 334 /* Take away 'Q' for verbose output */ 335 ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr); 336 triangulate(args, &in, &out, NULL); 337 } 338 ierr = PetscFree(in.pointlist);CHKERRQ(ierr); 339 ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr); 340 ierr = PetscFree(in.segmentlist);CHKERRQ(ierr); 341 ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr); 342 ierr = PetscFree(in.trianglelist);CHKERRQ(ierr); 343 344 { 345 DMLabel rlabel = NULL; 346 const PetscInt numCorners = 3; 347 const PetscInt numCells = out.numberoftriangles; 348 const PetscInt numVertices = out.numberofpoints; 349 PetscInt *cells; 350 PetscReal *meshCoords; 351 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 352 353 if (sizeof (PetscReal) == sizeof (out.pointlist[0])) { 354 meshCoords = (PetscReal *) out.pointlist; 355 } else { 356 PetscInt i; 357 358 ierr = PetscMalloc1(dim * numVertices,&meshCoords);CHKERRQ(ierr); 359 for (i = 0; i < dim * numVertices; i++) { 360 meshCoords[i] = (PetscReal) out.pointlist[i]; 361 } 362 } 363 if (sizeof (PetscInt) == sizeof (out.trianglelist[0])) { 364 cells = (PetscInt *) out.trianglelist; 365 } else { 366 PetscInt i; 367 368 ierr = PetscMalloc1(numCells * numCorners, &cells);CHKERRQ(ierr); 369 for (i = 0; i < numCells * numCorners; i++) { 370 cells[i] = (PetscInt) out.trianglelist[i]; 371 } 372 } 373 374 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 375 if (label) { 376 ierr = DMCreateLabel(*dmRefined, labelName);CHKERRQ(ierr); 377 ierr = DMGetLabel(*dmRefined, labelName, &rlabel);CHKERRQ(ierr); 378 } 379 if (sizeof (PetscReal) != sizeof (out.pointlist[0])) { 380 ierr = PetscFree(meshCoords);CHKERRQ(ierr); 381 } 382 if (sizeof (PetscInt) != sizeof (out.trianglelist[0])) { 383 ierr = PetscFree(cells);CHKERRQ(ierr); 384 } 385 /* Set labels */ 386 for (v = 0; v < numVertices; ++v) { 387 if (out.pointmarkerlist[v]) { 388 if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 389 } 390 } 391 if (interpolate) { 392 PetscInt e; 393 394 for (e = 0; e < out.numberofedges; e++) { 395 if (out.edgemarkerlist[e]) { 396 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 397 const PetscInt *edges; 398 PetscInt numEdges; 399 400 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 401 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 402 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 403 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 404 } 405 } 406 } 407 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 408 } 409 #if 0 /* Do not currently support holes */ 410 ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr); 411 #endif 412 ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr); 413 #if !defined(PETSC_USE_REAL_DOUBLE) 414 ierr = PetscFree(maxVolumes);CHKERRQ(ierr); 415 #endif 416 PetscFunctionReturn(0); 417 } 418