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