1 #include <petsc/private/dmpleximpl.h> /*I "petscdmplex.h" I*/ 2 3 #ifdef PETSC_HAVE_EGADS 4 #include <egads.h> 5 #endif 6 7 #include <ctetgen.h> 8 9 /* This is to fix the tetrahedron orientation from TetGen */ 10 static PetscErrorCode DMPlexInvertCells_CTetgen(PetscInt numCells, PetscInt numCorners, PetscInt cells[]) 11 { 12 PetscInt bound = numCells*numCorners, coff; 13 14 PetscFunctionBegin; 15 #define SWAP(a,b) do { PetscInt tmp = (a); (a) = (b); (b) = tmp; } while (0) 16 for (coff = 0; coff < bound; coff += numCorners) SWAP(cells[coff],cells[coff+1]); 17 #undef SWAP 18 PetscFunctionReturn(0); 19 } 20 21 PETSC_EXTERN PetscErrorCode DMPlexGenerate_CTetgen(DM boundary, PetscBool interpolate, DM *dm) 22 { 23 MPI_Comm comm; 24 const PetscInt dim = 3; 25 PLC *in, *out; 26 DMUniversalLabel universal; 27 PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, verbose = 0; 28 DMPlexInterpolatedFlag isInterpolated; 29 PetscMPIInt rank; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 ierr = PetscOptionsGetInt(NULL,((PetscObject) boundary)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr); 34 ierr = PetscObjectGetComm((PetscObject)boundary,&comm);CHKERRQ(ierr); 35 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 36 ierr = DMPlexIsInterpolatedCollective(boundary, &isInterpolated);CHKERRQ(ierr); 37 ierr = DMUniversalLabelCreate(boundary, &universal);CHKERRQ(ierr); 38 39 ierr = PLCCreate(&in);CHKERRQ(ierr); 40 ierr = PLCCreate(&out);CHKERRQ(ierr); 41 42 ierr = DMPlexGetDepthStratum(boundary, 0, &vStart, &vEnd);CHKERRQ(ierr); 43 in->numberofpoints = vEnd - vStart; 44 if (in->numberofpoints > 0) { 45 PetscSection coordSection; 46 Vec coordinates; 47 const PetscScalar *array; 48 49 ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr); 50 ierr = PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);CHKERRQ(ierr); 51 ierr = DMGetCoordinatesLocal(boundary, &coordinates);CHKERRQ(ierr); 52 ierr = DMGetCoordinateSection(boundary, &coordSection);CHKERRQ(ierr); 53 ierr = VecGetArrayRead(coordinates, &array);CHKERRQ(ierr); 54 for (v = vStart; v < vEnd; ++v) { 55 const PetscInt idx = v - vStart; 56 PetscInt off, d, m; 57 58 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 59 for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 60 ierr = DMLabelGetValue(universal->label, v, &m);CHKERRQ(ierr); 61 in->pointmarkerlist[idx] = (int) m; 62 } 63 ierr = VecRestoreArrayRead(coordinates, &array);CHKERRQ(ierr); 64 } 65 66 ierr = DMPlexGetHeightStratum(boundary, 1, &eStart, &eEnd);CHKERRQ(ierr); 67 in->numberofedges = eEnd - eStart; 68 if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) { 69 ierr = PetscMalloc1(in->numberofedges*2, &in->edgelist);CHKERRQ(ierr); 70 ierr = PetscMalloc1(in->numberofedges, &in->edgemarkerlist);CHKERRQ(ierr); 71 for (e = eStart; e < eEnd; ++e) { 72 const PetscInt idx = e - eStart; 73 const PetscInt *cone; 74 PetscInt coneSize, val; 75 76 ierr = DMPlexGetConeSize(boundary, e, &coneSize);CHKERRQ(ierr); 77 ierr = DMPlexGetCone(boundary, e, &cone);CHKERRQ(ierr); 78 in->edgelist[idx*2] = cone[0] - vStart; 79 in->edgelist[idx*2 + 1] = cone[1] - vStart; 80 81 ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr); 82 in->edgemarkerlist[idx] = (int) val; 83 } 84 } 85 86 ierr = DMPlexGetHeightStratum(boundary, 0, &fStart, &fEnd);CHKERRQ(ierr); 87 in->numberoffacets = fEnd - fStart; 88 if (in->numberoffacets > 0) { 89 ierr = PetscMalloc1(in->numberoffacets, &in->facetlist);CHKERRQ(ierr); 90 ierr = PetscMalloc1(in->numberoffacets, &in->facetmarkerlist);CHKERRQ(ierr); 91 for (f = fStart; f < fEnd; ++f) { 92 const PetscInt idx = f - fStart; 93 PetscInt *points = NULL, numPoints, p, numVertices = 0, v, m = -1; 94 polygon *poly; 95 96 in->facetlist[idx].numberofpolygons = 1; 97 ierr = PetscMalloc1(in->facetlist[idx].numberofpolygons, &in->facetlist[idx].polygonlist);CHKERRQ(ierr); 98 in->facetlist[idx].numberofholes = 0; 99 in->facetlist[idx].holelist = NULL; 100 101 ierr = DMPlexGetTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 102 for (p = 0; p < numPoints*2; p += 2) { 103 const PetscInt point = points[p]; 104 if ((point >= vStart) && (point < vEnd)) points[numVertices++] = point; 105 } 106 107 poly = in->facetlist[idx].polygonlist; 108 poly->numberofvertices = numVertices; 109 ierr = PetscMalloc1(poly->numberofvertices, &poly->vertexlist);CHKERRQ(ierr); 110 for (v = 0; v < numVertices; ++v) { 111 const PetscInt vIdx = points[v] - vStart; 112 poly->vertexlist[v] = vIdx; 113 } 114 ierr = DMLabelGetValue(universal->label, f, &m);CHKERRQ(ierr); 115 in->facetmarkerlist[idx] = (int) m; 116 ierr = DMPlexRestoreTransitiveClosure(boundary, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 117 } 118 } 119 if (rank == 0) { 120 TetGenOpts t; 121 122 ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr); 123 t.in = boundary; /* Should go away */ 124 t.plc = 1; 125 t.quality = 1; 126 t.edgesout = 1; 127 t.zeroindex = 1; 128 t.quiet = 1; 129 t.verbose = verbose; 130 #if 0 131 #ifdef PETSC_HAVE_EGADS 132 /* Need to add in more TetGen code */ 133 t.nobisect = 1; /* Add Y to preserve Surface Mesh for EGADS */ 134 #endif 135 #endif 136 137 ierr = TetGenCheckOpts(&t);CHKERRQ(ierr); 138 ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr); 139 } 140 { 141 const PetscInt numCorners = 4; 142 const PetscInt numCells = out->numberoftetrahedra; 143 const PetscInt numVertices = out->numberofpoints; 144 PetscReal *meshCoords = NULL; 145 PetscInt *cells = NULL; 146 147 if (sizeof (PetscReal) == sizeof (out->pointlist[0])) { 148 meshCoords = (PetscReal *) out->pointlist; 149 } else { 150 PetscInt i; 151 152 ierr = PetscMalloc1(dim * numVertices, &meshCoords);CHKERRQ(ierr); 153 for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i]; 154 } 155 if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) { 156 cells = (PetscInt *) out->tetrahedronlist; 157 } else { 158 PetscInt i; 159 160 ierr = PetscMalloc1(numCells * numCorners, &cells);CHKERRQ(ierr); 161 for (i = 0; i < numCells * numCorners; i++) cells[i] = (PetscInt) out->tetrahedronlist[i]; 162 } 163 164 ierr = DMPlexInvertCells_CTetgen(numCells, numCorners, cells);CHKERRQ(ierr); 165 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dm);CHKERRQ(ierr); 166 if (sizeof (PetscReal) != sizeof (out->pointlist[0])) { 167 ierr = PetscFree(meshCoords);CHKERRQ(ierr); 168 } 169 if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) { 170 ierr = PetscFree(cells);CHKERRQ(ierr); 171 } 172 173 /* Set labels */ 174 ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dm);CHKERRQ(ierr); 175 for (v = 0; v < numVertices; ++v) { 176 if (out->pointmarkerlist[v]) { 177 ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr); 178 } 179 } 180 if (interpolate) { 181 PetscInt e; 182 183 for (e = 0; e < out->numberofedges; e++) { 184 if (out->edgemarkerlist[e]) { 185 const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells}; 186 const PetscInt *edges; 187 PetscInt numEdges; 188 189 ierr = DMPlexGetJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 190 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 191 ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr); 192 ierr = DMPlexRestoreJoin(*dm, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 193 } 194 } 195 for (f = 0; f < out->numberoftrifaces; f++) { 196 if (out->trifacemarkerlist[f]) { 197 const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells}; 198 const PetscInt *faces; 199 PetscInt numFaces; 200 201 ierr = DMPlexGetFullJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 202 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 203 ierr = DMUniversalLabelSetLabelValue(universal, *dm, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr); 204 ierr = DMPlexRestoreJoin(*dm, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 205 } 206 } 207 } 208 209 #ifdef PETSC_HAVE_EGADS 210 { 211 DMLabel bodyLabel; 212 PetscContainer modelObj; 213 PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd; 214 ego *bodies; 215 ego model, geom; 216 int Nb, oclass, mtype, *senses; 217 218 /* Get Attached EGADS Model from Original DMPlex */ 219 ierr = PetscObjectQuery((PetscObject) boundary, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr); 220 if (modelObj) { 221 ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr); 222 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 223 /* Transfer EGADS Model to Volumetric Mesh */ 224 ierr = PetscObjectCompose((PetscObject) *dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 225 226 /* Set Cell Labels */ 227 ierr = DMGetLabel(*dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 228 ierr = DMPlexGetHeightStratum(*dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 229 ierr = DMPlexGetHeightStratum(*dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 230 ierr = DMPlexGetDepthStratum(*dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 231 232 for (c = cStart; c < cEnd; ++c) { 233 PetscReal centroid[3] = {0., 0., 0.}; 234 PetscInt b; 235 236 /* Deterimine what body the cell's centroid is located in */ 237 if (!interpolate) { 238 PetscSection coordSection; 239 Vec coordinates; 240 PetscScalar *coords = NULL; 241 PetscInt coordSize, s, d; 242 243 ierr = DMGetCoordinatesLocal(*dm, &coordinates);CHKERRQ(ierr); 244 ierr = DMGetCoordinateSection(*dm, &coordSection);CHKERRQ(ierr); 245 ierr = DMPlexVecGetClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr); 246 for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d]; 247 ierr = DMPlexVecRestoreClosure(*dm, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr); 248 } else { 249 ierr = DMPlexComputeCellGeometryFVM(*dm, c, NULL, centroid, NULL);CHKERRQ(ierr); 250 } 251 for (b = 0; b < Nb; ++b) { 252 if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break; 253 } 254 if (b < Nb) { 255 PetscInt cval = b, eVal, fVal; 256 PetscInt *closure = NULL, Ncl, cl; 257 258 ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr); 259 ierr = DMPlexGetTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr); 260 for (cl = 0; cl < Ncl; ++cl) { 261 const PetscInt p = closure[cl*2]; 262 263 if (p >= eStart && p < eEnd) { 264 ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr); 265 if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);} 266 } 267 if (p >= fStart && p < fEnd) { 268 ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr); 269 if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);} 270 } 271 } 272 ierr = DMPlexRestoreTransitiveClosure(*dm, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr); 273 } 274 } 275 } 276 } 277 #endif 278 ierr = DMPlexSetRefinementUniform(*dm, PETSC_FALSE);CHKERRQ(ierr); 279 } 280 281 ierr = DMUniversalLabelDestroy(&universal);CHKERRQ(ierr); 282 ierr = PLCDestroy(&in);CHKERRQ(ierr); 283 ierr = PLCDestroy(&out);CHKERRQ(ierr); 284 PetscFunctionReturn(0); 285 } 286 287 PETSC_EXTERN PetscErrorCode DMPlexRefine_CTetgen(DM dm, PetscReal *maxVolumes, DM *dmRefined) 288 { 289 MPI_Comm comm; 290 const PetscInt dim = 3; 291 PLC *in, *out; 292 DMUniversalLabel universal; 293 PetscInt vStart, vEnd, v, eStart, eEnd, e, fStart, fEnd, f, cStart, cEnd, c, verbose = 0; 294 DMPlexInterpolatedFlag isInterpolated; 295 PetscMPIInt rank; 296 PetscErrorCode ierr; 297 298 PetscFunctionBegin; 299 ierr = PetscOptionsGetInt(NULL,((PetscObject) dm)->prefix, "-ctetgen_verbose", &verbose, NULL);CHKERRQ(ierr); 300 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr); 301 ierr = MPI_Comm_rank(comm, &rank);CHKERRMPI(ierr); 302 ierr = DMPlexIsInterpolatedCollective(dm, &isInterpolated);CHKERRQ(ierr); 303 ierr = DMUniversalLabelCreate(dm, &universal);CHKERRQ(ierr); 304 305 ierr = PLCCreate(&in);CHKERRQ(ierr); 306 ierr = PLCCreate(&out);CHKERRQ(ierr); 307 308 ierr = DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);CHKERRQ(ierr); 309 in->numberofpoints = vEnd - vStart; 310 if (in->numberofpoints > 0) { 311 PetscSection coordSection; 312 Vec coordinates; 313 PetscScalar *array; 314 315 ierr = PetscMalloc1(in->numberofpoints*dim, &in->pointlist);CHKERRQ(ierr); 316 ierr = PetscMalloc1(in->numberofpoints, &in->pointmarkerlist);CHKERRQ(ierr); 317 ierr = DMGetCoordinatesLocal(dm, &coordinates);CHKERRQ(ierr); 318 ierr = DMGetCoordinateSection(dm, &coordSection);CHKERRQ(ierr); 319 ierr = VecGetArray(coordinates, &array);CHKERRQ(ierr); 320 for (v = vStart; v < vEnd; ++v) { 321 const PetscInt idx = v - vStart; 322 PetscInt off, d, m; 323 324 ierr = PetscSectionGetOffset(coordSection, v, &off);CHKERRQ(ierr); 325 for (d = 0; d < dim; ++d) in->pointlist[idx*dim + d] = PetscRealPart(array[off+d]); 326 ierr = DMLabelGetValue(universal->label, v, &m);CHKERRQ(ierr); 327 in->pointmarkerlist[idx] = (int) m; 328 } 329 ierr = VecRestoreArray(coordinates, &array);CHKERRQ(ierr); 330 } 331 332 ierr = DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd);CHKERRQ(ierr); 333 in->numberofedges = eEnd - eStart; 334 if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberofedges > 0) { 335 ierr = PetscMalloc1(in->numberofedges * 2, &in->edgelist);CHKERRQ(ierr); 336 ierr = PetscMalloc1(in->numberofedges, &in->edgemarkerlist);CHKERRQ(ierr); 337 for (e = eStart; e < eEnd; ++e) { 338 const PetscInt idx = e - eStart; 339 const PetscInt *cone; 340 PetscInt coneSize, val; 341 342 ierr = DMPlexGetConeSize(dm, e, &coneSize);CHKERRQ(ierr); 343 ierr = DMPlexGetCone(dm, e, &cone);CHKERRQ(ierr); 344 in->edgelist[idx*2] = cone[0] - vStart; 345 in->edgelist[idx*2 + 1] = cone[1] - vStart; 346 347 ierr = DMLabelGetValue(universal->label, e, &val);CHKERRQ(ierr); 348 in->edgemarkerlist[idx] = (int) val; 349 } 350 } 351 352 ierr = DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);CHKERRQ(ierr); 353 in->numberoftrifaces = 0; 354 for (f = fStart; f < fEnd; ++f) { 355 PetscInt supportSize; 356 357 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 358 if (supportSize == 1) ++in->numberoftrifaces; 359 } 360 if (isInterpolated == DMPLEX_INTERPOLATED_FULL && in->numberoftrifaces > 0) { 361 PetscInt tf = 0; 362 363 ierr = PetscMalloc1(in->numberoftrifaces*3, &in->trifacelist);CHKERRQ(ierr); 364 ierr = PetscMalloc1(in->numberoftrifaces, &in->trifacemarkerlist);CHKERRQ(ierr); 365 for (f = fStart; f < fEnd; ++f) { 366 PetscInt *points = NULL; 367 PetscInt supportSize, numPoints, p, Nv = 0, val; 368 369 ierr = DMPlexGetSupportSize(dm, f, &supportSize);CHKERRQ(ierr); 370 if (supportSize != 1) continue; 371 ierr = DMPlexGetTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 372 for (p = 0; p < numPoints*2; p += 2) { 373 const PetscInt point = points[p]; 374 if ((point >= vStart) && (point < vEnd)) in->trifacelist[tf*3 + Nv++] = point - vStart; 375 } 376 ierr = DMPlexRestoreTransitiveClosure(dm, f, PETSC_TRUE, &numPoints, &points);CHKERRQ(ierr); 377 if (Nv != 3) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Face %D has %D vertices, not 3", f, Nv); 378 ierr = DMLabelGetValue(universal->label, f, &val);CHKERRQ(ierr); 379 in->trifacemarkerlist[tf] = (int) val; 380 ++tf; 381 } 382 } 383 384 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 385 in->numberofcorners = 4; 386 in->numberoftetrahedra = cEnd - cStart; 387 in->tetrahedronvolumelist = maxVolumes; 388 if (in->numberoftetrahedra > 0) { 389 ierr = PetscMalloc1(in->numberoftetrahedra*in->numberofcorners, &in->tetrahedronlist);CHKERRQ(ierr); 390 for (c = cStart; c < cEnd; ++c) { 391 const PetscInt idx = c - cStart; 392 PetscInt *closure = NULL; 393 PetscInt closureSize; 394 395 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 396 if ((closureSize != 5) && (closureSize != 15)) SETERRQ1(comm, PETSC_ERR_ARG_WRONG, "Mesh has cell which is not a tetrahedron, %D vertices in closure", closureSize); 397 for (v = 0; v < 4; ++v) in->tetrahedronlist[idx*in->numberofcorners + v] = closure[(v+closureSize-4)*2] - vStart; 398 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &closureSize, &closure);CHKERRQ(ierr); 399 } 400 } 401 402 if (rank == 0) { 403 TetGenOpts t; 404 405 ierr = TetGenOptsInitialize(&t);CHKERRQ(ierr); 406 407 t.in = dm; /* Should go away */ 408 t.refine = 1; 409 t.varvolume = 1; 410 t.quality = 1; 411 t.edgesout = 1; 412 t.zeroindex = 1; 413 t.quiet = 1; 414 t.verbose = verbose; /* Change this */ 415 416 ierr = TetGenCheckOpts(&t);CHKERRQ(ierr); 417 ierr = TetGenTetrahedralize(&t, in, out);CHKERRQ(ierr); 418 } 419 420 in->tetrahedronvolumelist = NULL; 421 { 422 const PetscInt numCorners = 4; 423 const PetscInt numCells = out->numberoftetrahedra; 424 const PetscInt numVertices = out->numberofpoints; 425 PetscReal *meshCoords = NULL; 426 PetscInt *cells = NULL; 427 PetscBool interpolate = isInterpolated == DMPLEX_INTERPOLATED_FULL ? PETSC_TRUE : PETSC_FALSE; 428 429 if (sizeof (PetscReal) == sizeof (out->pointlist[0])) { 430 meshCoords = (PetscReal *) out->pointlist; 431 } else { 432 PetscInt i; 433 434 ierr = PetscMalloc1(dim * numVertices, &meshCoords);CHKERRQ(ierr); 435 for (i = 0; i < dim * numVertices; ++i) meshCoords[i] = (PetscReal) out->pointlist[i]; 436 } 437 if (sizeof (PetscInt) == sizeof (out->tetrahedronlist[0])) { 438 cells = (PetscInt *) out->tetrahedronlist; 439 } else { 440 PetscInt i; 441 442 ierr = PetscMalloc1(numCells * numCorners, &cells);CHKERRQ(ierr); 443 for (i = 0; i < numCells * numCorners; ++i) cells[i] = (PetscInt) out->tetrahedronlist[i]; 444 } 445 446 ierr = DMPlexInvertCells_CTetgen(numCells, numCorners, cells);CHKERRQ(ierr); 447 ierr = DMPlexCreateFromCellListPetsc(comm, dim, numCells, numVertices, numCorners, interpolate, cells, dim, meshCoords, dmRefined);CHKERRQ(ierr); 448 if (sizeof (PetscReal) != sizeof (out->pointlist[0])) {ierr = PetscFree(meshCoords);CHKERRQ(ierr);} 449 if (sizeof (PetscInt) != sizeof (out->tetrahedronlist[0])) {ierr = PetscFree(cells);CHKERRQ(ierr);} 450 451 /* Set labels */ 452 ierr = DMUniversalLabelCreateLabels(universal, PETSC_TRUE, *dmRefined);CHKERRQ(ierr); 453 for (v = 0; v < numVertices; ++v) { 454 if (out->pointmarkerlist[v]) { 455 ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, v+numCells, out->pointmarkerlist[v]);CHKERRQ(ierr); 456 } 457 } 458 if (interpolate) { 459 PetscInt e, f; 460 461 for (e = 0; e < out->numberofedges; e++) { 462 if (out->edgemarkerlist[e]) { 463 const PetscInt vertices[2] = {out->edgelist[e*2+0]+numCells, out->edgelist[e*2+1]+numCells}; 464 const PetscInt *edges; 465 PetscInt numEdges; 466 467 ierr = DMPlexGetJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 468 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Two vertices must cover only one edge, not %D", numEdges); 469 ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, edges[0], out->edgemarkerlist[e]);CHKERRQ(ierr); 470 ierr = DMPlexRestoreJoin(*dmRefined, 2, vertices, &numEdges, &edges);CHKERRQ(ierr); 471 } 472 } 473 for (f = 0; f < out->numberoftrifaces; f++) { 474 if (out->trifacemarkerlist[f]) { 475 const PetscInt vertices[3] = {out->trifacelist[f*3+0]+numCells, out->trifacelist[f*3+1]+numCells, out->trifacelist[f*3+2]+numCells}; 476 const PetscInt *faces; 477 PetscInt numFaces; 478 479 ierr = DMPlexGetFullJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 480 if (numFaces != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Three vertices must cover only one face, not %D", numFaces); 481 ierr = DMUniversalLabelSetLabelValue(universal, *dmRefined, PETSC_TRUE, faces[0], out->trifacemarkerlist[f]);CHKERRQ(ierr); 482 ierr = DMPlexRestoreJoin(*dmRefined, 3, vertices, &numFaces, &faces);CHKERRQ(ierr); 483 } 484 } 485 } 486 487 #ifdef PETSC_HAVE_EGADS 488 { 489 DMLabel bodyLabel; 490 PetscContainer modelObj; 491 PetscInt cStart, cEnd, c, eStart, eEnd, fStart, fEnd; 492 ego *bodies; 493 ego model, geom; 494 int Nb, oclass, mtype, *senses; 495 496 /* Get Attached EGADS Model from Original DMPlex */ 497 ierr = PetscObjectQuery((PetscObject) dm, "EGADS Model", (PetscObject *) &modelObj);CHKERRQ(ierr); 498 if (modelObj) { 499 ierr = PetscContainerGetPointer(modelObj, (void **) &model);CHKERRQ(ierr); 500 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 501 /* Transfer EGADS Model to Volumetric Mesh */ 502 ierr = PetscObjectCompose((PetscObject) *dmRefined, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 503 504 /* Set Cell Labels */ 505 ierr = DMGetLabel(*dmRefined, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 506 ierr = DMPlexGetHeightStratum(*dmRefined, 0, &cStart, &cEnd);CHKERRQ(ierr); 507 ierr = DMPlexGetHeightStratum(*dmRefined, 1, &fStart, &fEnd);CHKERRQ(ierr); 508 ierr = DMPlexGetDepthStratum(*dmRefined, 1, &eStart, &eEnd);CHKERRQ(ierr); 509 510 for (c = cStart; c < cEnd; ++c) { 511 PetscReal centroid[3] = {0., 0., 0.}; 512 PetscInt b; 513 514 /* Deterimine what body the cell's centroid is located in */ 515 if (!interpolate) { 516 PetscSection coordSection; 517 Vec coordinates; 518 PetscScalar *coords = NULL; 519 PetscInt coordSize, s, d; 520 521 ierr = DMGetCoordinatesLocal(*dmRefined, &coordinates);CHKERRQ(ierr); 522 ierr = DMGetCoordinateSection(*dmRefined, &coordSection);CHKERRQ(ierr); 523 ierr = DMPlexVecGetClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr); 524 for (s = 0; s < coordSize; ++s) for (d = 0; d < dim; ++d) centroid[d] += coords[s*dim+d]; 525 ierr = DMPlexVecRestoreClosure(*dmRefined, coordSection, coordinates, c, &coordSize, &coords);CHKERRQ(ierr); 526 } else { 527 ierr = DMPlexComputeCellGeometryFVM(*dmRefined, c, NULL, centroid, NULL);CHKERRQ(ierr); 528 } 529 for (b = 0; b < Nb; ++b) { 530 if (EG_inTopology(bodies[b], centroid) == EGADS_SUCCESS) break; 531 } 532 if (b < Nb) { 533 PetscInt cval = b, eVal, fVal; 534 PetscInt *closure = NULL, Ncl, cl; 535 536 ierr = DMLabelSetValue(bodyLabel, c, cval);CHKERRQ(ierr); 537 ierr = DMPlexGetTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr); 538 for (cl = 0; cl < Ncl; cl += 2) { 539 const PetscInt p = closure[cl]; 540 541 if (p >= eStart && p < eEnd) { 542 ierr = DMLabelGetValue(bodyLabel, p, &eVal);CHKERRQ(ierr); 543 if (eVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);} 544 } 545 if (p >= fStart && p < fEnd) { 546 ierr = DMLabelGetValue(bodyLabel, p, &fVal);CHKERRQ(ierr); 547 if (fVal < 0) {ierr = DMLabelSetValue(bodyLabel, p, cval);CHKERRQ(ierr);} 548 } 549 } 550 ierr = DMPlexRestoreTransitiveClosure(*dmRefined, c, PETSC_TRUE, &Ncl, &closure);CHKERRQ(ierr); 551 } 552 } 553 } 554 } 555 #endif 556 ierr = DMPlexSetRefinementUniform(*dmRefined, PETSC_FALSE);CHKERRQ(ierr); 557 } 558 ierr = DMUniversalLabelDestroy(&universal);CHKERRQ(ierr); 559 ierr = PLCDestroy(&in);CHKERRQ(ierr); 560 ierr = PLCDestroy(&out);CHKERRQ(ierr); 561 PetscFunctionReturn(0); 562 } 563