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, PetscReal *inmaxVolumes, 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 double *maxVolumes; 207 208 PetscFunctionBegin; 209 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 210 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 211 ierr = InitInput_Triangle(&in);CHKERRQ(ierr); 212 ierr = InitOutput_Triangle(&out);CHKERRQ(ierr); 213 ierr = DMPlexGetDepth(dm, &depth);CHKERRQ(ierr); 214 ierr = MPIU_Allreduce(&depth, &depthGlobal, 1, MPIU_INT, MPI_MAX, comm);CHKERRQ(ierr); 215 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 216 ierr = DMGetLabel(dm, labelName, &label);CHKERRQ(ierr); 217 218 in.numberofpoints = vEnd - vStart; 219 if (in.numberofpoints > 0) { 220 PetscSection coordSection; 221 Vec coordinates; 222 PetscScalar *array; 223 224 ierr = PetscMalloc1(in.numberofpoints*dim, &in.pointlist);CHKERRQ(ierr); 225 ierr = PetscMalloc1(in.numberofpoints, &in.pointmarkerlist);CHKERRQ(ierr); 226 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 227 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 228 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 229 for (v = vStart; v < vEnd; ++v) { 230 const PetscInt idx = v - vStart; 231 PetscInt off, d; 232 233 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 234 for (d = 0; d < dim; ++d) { 235 in.pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 236 } 237 if (label) {ierr = DMLabelGetValue(label, v, &in.pointmarkerlist[idx]);CHKERRQ(ierr);} 238 } 239 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 240 } 241 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 242 243 in.numberofcorners = 3; 244 in.numberoftriangles = cEnd - cStart; 245 246 #if !defined(PETSC_USE_REAL_DOUBLE) 247 ierr = PetscMalloc1(cEnd - cStart,&maxVolumes);CHKERRQ(ierr); 248 for (c = 0; c < cEnd-cStart; ++c) maxVolumes[c] = (double)inmaxVolumes[c]; 249 #else 250 maxVolumes = inmaxVolumes; 251 #endif 252 253 in.trianglearealist = (double*) maxVolumes; 254 if (in.numberoftriangles > 0) { 255 ierr = PetscMalloc1(in.numberoftriangles*in.numberofcorners, &in.trianglelist);CHKERRQ(ierr); 256 for (c = cStart; c < cEnd; ++c) { 257 const PetscInt idx = c - cStart; 258 PetscInt *closure = NULL; 259 PetscInt closureSize; 260 261 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 262 if ((closureSize != 4) && (closureSize != 7)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a triangle, %D vertices in closure", closureSize); 263 for (v = 0; v < 3; ++v) { 264 in.trianglelist[idx*in.numberofcorners + v] = closure[(v+closureSize-3)*2] - vStart; 265 } 266 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 267 } 268 } 269 /* TODO: Segment markers are missing on input */ 270 #if 0 /* Do not currently support holes */ 271 PetscReal *holeCoords; 272 PetscInt h, d; 273 274 ierr = DMPlexGetHoles(boundary, &in.numberofholes, &holeCords);CHKERRQ(ierr); 275 if (in.numberofholes > 0) { 276 ierr = PetscMalloc1(in.numberofholes*dim, &in.holelist);CHKERRQ(ierr); 277 for (h = 0; h < in.numberofholes; ++h) { 278 for (d = 0; d < dim; ++d) { 279 in.holelist[h*dim+d] = holeCoords[h*dim+d]; 280 } 281 } 282 } 283 #endif 284 if (!rank) { 285 char args[32]; 286 287 /* Take away 'Q' for verbose output */ 288 ierr = PetscStrcpy(args, "pqezQra");CHKERRQ(ierr); 289 triangulate(args, &in, &out, NULL); 290 } 291 ierr = PetscFree(in.pointlist);CHKERRQ(ierr); 292 ierr = PetscFree(in.pointmarkerlist);CHKERRQ(ierr); 293 ierr = PetscFree(in.segmentlist);CHKERRQ(ierr); 294 ierr = PetscFree(in.segmentmarkerlist);CHKERRQ(ierr); 295 ierr = PetscFree(in.trianglelist);CHKERRQ(ierr); 296 297 { 298 DMLabel rlabel = NULL; 299 const PetscInt numCorners = 3; 300 const PetscInt numCells = out.numberoftriangles; 301 const PetscInt numVertices = out.numberofpoints; 302 const int *cells = out.trianglelist; 303 const double *meshCoords = out.pointlist; 304 PetscBool interpolate = depthGlobal > 1 ? PETSC_TRUE : PETSC_FALSE; 305 306 ierr = DMPlexCreateFromCellList(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 307 if (label) {ierr = DMCreateLabel(*dmRefined, labelName); ierr = DMGetLabel(*dmRefined, labelName, &rlabel);} 308 /* Set labels */ 309 for (v = 0; v < numVertices; ++v) { 310 if (out.pointmarkerlist[v]) { 311 if (rlabel) {ierr = DMLabelSetValue(rlabel, v+numCells, out.pointmarkerlist[v]);CHKERRQ(ierr);} 312 } 313 } 314 if (interpolate) { 315 PetscInt e; 316 317 for (e = 0; e < out.numberofedges; e++) { 318 if (out.edgemarkerlist[e]) { 319 const PetscInt vertices[2] = {out.edgelist[e*2+0]+numCells, out.edgelist[e*2+1]+numCells}; 320 const PetscInt *edges; 321 PetscInt numEdges; 322 323 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 324 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 325 if (rlabel) {ierr = DMLabelSetValue(rlabel, edges[0], out.edgemarkerlist[e]);CHKERRQ(ierr);} 326 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 327 } 328 } 329 } 330 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 331 } 332 #if 0 /* Do not currently support holes */ 333 ierr = DMPlexCopyHoles(*dm, boundary);CHKERRQ(ierr); 334 #endif 335 ierr = FiniOutput_Triangle(&out);CHKERRQ(ierr); 336 #if !defined(PETSC_USE_REAL_DOUBLE) 337 ierr = PetscFree(maxVolumes);CHKERRQ(ierr); 338 #endif 339 PetscFunctionReturn(0); 340 } 341