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