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