1 static const char help[] = "Test of EGADSLite CAD functionality"; 2 3 #include <petscdmplex.h> 4 #include <petsc/private/hashmapi.h> 5 6 #include <egads.h> 7 #include <petsc.h> 8 9 typedef struct { 10 char filename[PETSC_MAX_PATH_LEN]; 11 } AppCtx; 12 13 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 14 { 15 PetscErrorCode ierr; 16 17 PetscFunctionBeginUser; 18 options->filename[0] = '\0'; 19 20 ierr = PetscOptionsBegin(comm, "", "EGADSPlex Problem Options", "EGADSLite");CHKERRQ(ierr); 21 ierr = PetscOptionsString("-filename", "The EGADSLite file", "ex9.c", options->filename, options->filename, sizeof(options->filename), NULL);CHKERRQ(ierr); 22 ierr = PetscOptionsEnd(); 23 PetscFunctionReturn(0); 24 } 25 26 static PetscErrorCode PrintEGADS(ego model) 27 { 28 ego geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 29 int oclass, mtype, *senses; 30 int Nb, b; 31 PetscErrorCode ierr; 32 33 PetscFunctionBeginUser; 34 /* test bodyTopo functions */ 35 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &Nb, &bodies, &senses);CHKERRQ(ierr); 36 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", Nb);CHKERRQ(ierr); 37 38 for (b = 0; b < Nb; ++b) { 39 ego body = bodies[b]; 40 int id, Nsh, Nf, Nl, l, Ne, e, Nv, v; 41 42 /* Output Basic Model Topology */ 43 ierr = EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);CHKERRQ(ierr); 44 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh);CHKERRQ(ierr); 45 EG_free(objs); 46 47 ierr = EG_getBodyTopos(body, NULL, FACE, &Nf, &objs);CHKERRQ(ierr); 48 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf);CHKERRQ(ierr); 49 EG_free(objs); 50 51 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 52 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl);CHKERRQ(ierr); 53 54 ierr = EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs);CHKERRQ(ierr); 55 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne);CHKERRQ(ierr); 56 EG_free(objs); 57 58 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &objs);CHKERRQ(ierr); 59 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv);CHKERRQ(ierr); 60 EG_free(objs); 61 62 for (l = 0; l < Nl; ++l) { 63 ego loop = lobjs[l]; 64 65 id = EG_indexBodyTopo(body, loop); 66 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id);CHKERRQ(ierr); 67 68 /* Get EDGE info which associated with the current LOOP */ 69 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 70 71 for (e = 0; e < Ne; ++e) { 72 ego edge = objs[e]; 73 74 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 75 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d\n", id);CHKERRQ(ierr); 76 77 double range[4] = {0., 0., 0., 0.}; 78 double point[3] = {0., 0., 0.}; 79 int peri; 80 81 ierr = EG_getRange(objs[e], range, &peri); 82 ierr = PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]); 83 84 /* Get NODE info which associated with the current EDGE */ 85 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);CHKERRQ(ierr); 86 87 for (v = 0; v < Nv; ++v) { 88 ego vertex = nobjs[v]; 89 double limits[4]; 90 int dummy; 91 92 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 93 id = EG_indexBodyTopo(body, vertex); 94 ierr = PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);CHKERRQ(ierr); 95 ierr = PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]); 96 97 point[0] = point[0] + limits[0]; 98 point[1] = point[1] + limits[1]; 99 point[2] = point[2] + limits[2]; 100 } 101 } 102 } 103 EG_free(lobjs); 104 } 105 PetscFunctionReturn(0); 106 } 107 108 int main(int argc, char *argv[]) 109 { 110 DMLabel bodyLabel, faceLabel, edgeLabel; 111 PetscInt cStart, cEnd, c; 112 /* EGADSLite variables */ 113 ego context, model, geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 114 int oclass, mtype, nbodies, *senses; 115 int b; 116 /* PETSc variables */ 117 DM dm; 118 PetscHMapI edgeMap = NULL; 119 PetscInt dim = -1, cdim = -1, numCorners = 0, maxCorners = 0, numVertices = 0, newVertices = 0, numEdges = 0, numCells = 0, newCells = 0, numQuads = 0, cOff = 0, fOff = 0; 120 PetscInt *cells = NULL, *cone = NULL; 121 PetscReal *coords = NULL; 122 MPI_Comm comm; 123 PetscMPIInt rank; 124 AppCtx ctx; 125 PetscErrorCode ierr; 126 127 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 128 comm = PETSC_COMM_WORLD; 129 ierr = ProcessOptions(comm, &ctx);CHKERRQ(ierr); 130 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 131 if (!rank) { 132 /* Open EGADs file and load EGADs model data */ 133 ierr = EG_open(&context);CHKERRQ(ierr); 134 ierr = EG_loadModel(context, 0, ctx.filename, &model);CHKERRQ(ierr); 135 136 ierr = PrintEGADS(model);CHKERRQ(ierr); 137 138 /* --------------------------------------------------------------------------------------------------- 139 Generate Petsc Plex 140 Get all Nodes in model, record coordinates in a correctly formatted array 141 Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array 142 We need to uniformly refine the initial geometry to guarantee a valid mesh 143 */ 144 145 /* Calculate cell and vertex sizes */ 146 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);CHKERRQ(ierr); 147 ierr = PetscHMapICreate(&edgeMap);CHKERRQ(ierr); 148 numEdges = 0; 149 for (b = 0; b < nbodies; ++b) { 150 ego body = bodies[b]; 151 int id, Nl, l, Nv, v; 152 153 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 154 for (l = 0; l < Nl; ++l) { 155 ego loop = lobjs[l]; 156 int Ner = 0, Ne, e, Nc; 157 158 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 159 for (e = 0; e < Ne; ++e) { 160 ego edge = objs[e]; 161 int Nv, id; 162 PetscHashIter iter; 163 PetscBool found; 164 165 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 166 if (mtype == DEGENERATE) continue; 167 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 168 ierr = PetscHMapIFind(edgeMap, id-1, &iter, &found);CHKERRQ(ierr); 169 if (!found) {ierr = PetscHMapISet(edgeMap, id-1, numEdges++);CHKERRQ(ierr);} 170 ++Ner; 171 } 172 if (Ner == 2) {Nc = 2;} 173 else if (Ner == 3) {Nc = 4;} 174 else if (Ner == 4) {Nc = 8; ++numQuads;} 175 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot support loop with %d edges", Ner); 176 numCells += Nc; 177 newCells += Nc-1; 178 maxCorners = PetscMax(Ner*2+1, maxCorners); 179 } 180 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 181 for (v = 0; v < Nv; ++v) { 182 ego vertex = nobjs[v]; 183 184 id = EG_indexBodyTopo(body, vertex); 185 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 186 numVertices = PetscMax(id, numVertices); 187 } 188 EG_free(lobjs); 189 EG_free(nobjs); 190 } 191 ierr = PetscHMapIGetSize(edgeMap, &numEdges);CHKERRQ(ierr); 192 newVertices = numEdges + numQuads; 193 numVertices += newVertices; 194 ierr = PetscPrintf(PETSC_COMM_SELF, "\nPLEX Input Array Checkouts\n");CHKERRQ(ierr); 195 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Cells = %D (%D)\n", numCells, newCells);CHKERRQ(ierr); 196 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Vertices = %D (%D) \n", numVertices, newVertices);CHKERRQ(ierr); 197 198 dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 199 cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 200 numCorners = 3; /* Split cells into triangles */ 201 ierr = PetscMalloc3(numVertices*cdim, &coords, numCells*numCorners, &cells, maxCorners, &cone);CHKERRQ(ierr); 202 203 /* Get vertex coordinates */ 204 for (b = 0; b < nbodies; ++b) { 205 ego body = bodies[b]; 206 int id, Nv, v; 207 208 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 209 for (v = 0; v < Nv; ++v) { 210 ego vertex = nobjs[v]; 211 double limits[4]; 212 int dummy; 213 214 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 215 id = EG_indexBodyTopo(body, vertex);CHKERRQ(ierr); 216 coords[(id-1)*cdim+0] = limits[0]; 217 coords[(id-1)*cdim+1] = limits[1]; 218 coords[(id-1)*cdim+2] = limits[2]; 219 ierr = PetscPrintf(PETSC_COMM_SELF, " Node ID = %d (%d)\n", id, id-1); 220 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf) \n \n", coords[(id-1)*cdim+0], coords[(id-1)*cdim+1],coords[(id-1)*cdim+2]); 221 } 222 EG_free(nobjs); 223 } 224 ierr = PetscHMapIClear(edgeMap);CHKERRQ(ierr); 225 fOff = numVertices - newVertices + numEdges; 226 numEdges = 0; 227 numQuads = 0; 228 for (b = 0; b < nbodies; ++b) { 229 ego body = bodies[b]; 230 int Nl, l; 231 232 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 233 for (l = 0; l < Nl; ++l) { 234 ego loop = lobjs[l]; 235 int lid, Ner = 0, Ne, e; 236 237 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 238 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 239 for (e = 0; e < Ne; ++e) { 240 ego edge = objs[e]; 241 int eid, Nv; 242 PetscHashIter iter; 243 PetscBool found; 244 245 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 246 if (mtype == DEGENERATE) continue; 247 ++Ner; 248 eid = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 249 ierr = PetscHMapIFind(edgeMap, eid-1, &iter, &found);CHKERRQ(ierr); 250 if (!found) { 251 PetscInt v = numVertices - newVertices + numEdges; 252 double range[4], params[3] = {0., 0., 0.}, result[18]; 253 int periodic[2]; 254 255 ierr = PetscHMapISet(edgeMap, eid-1, numEdges++);CHKERRQ(ierr); 256 ierr = EG_getRange(edge, range, periodic); CHKERRQ(ierr); 257 params[0] = 0.5*(range[0] + range[1]); 258 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 259 coords[v*cdim+0] = result[0]; 260 coords[v*cdim+1] = result[1]; 261 coords[v*cdim+2] = result[2]; 262 ierr = PetscPrintf(PETSC_COMM_SELF, " Edge ID = %d (%D) \n", eid, v); 263 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf)\n", coords[v*cdim+0], coords[v*cdim+1],coords[v*cdim+2]); 264 params[0] = range[0]; 265 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 266 ierr = PetscPrintf(PETSC_COMM_SELF, " between (%lf, %lf, %lf)", result[0], result[1], result[2]); 267 params[0] = range[1]; 268 ierr = EG_evaluate(edge, params, result);CHKERRQ(ierr); 269 ierr = PetscPrintf(PETSC_COMM_SELF, " and (%lf, %lf, %lf)\n\n", result[0], result[1], result[2]); 270 } 271 } 272 if (Ner == 4) { 273 PetscInt v = fOff + numQuads++; 274 ego *fobjs; 275 double range[4], params[3] = {0., 0., 0.}, result[18]; 276 int Nf, face, periodic[2]; 277 278 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 279 if (Nf != 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Loop %d has %d faces, instead of 1", lid-1, Nf); 280 face = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 281 ierr = EG_getRange(fobjs[0], range, periodic); CHKERRQ(ierr); 282 params[0] = 0.5*(range[0] + range[1]); 283 params[1] = 0.5*(range[2] + range[3]); 284 ierr = EG_evaluate(fobjs[0], params, result);CHKERRQ(ierr); 285 coords[v*cdim+0] = result[0]; 286 coords[v*cdim+1] = result[1]; 287 coords[v*cdim+2] = result[2]; 288 ierr = PetscPrintf(PETSC_COMM_SELF, " Face ID = %d (%D) \n", face-1, v); 289 ierr = PetscPrintf(PETSC_COMM_SELF, " (x,y,z) = (%lf, %lf, %lf) \n \n", coords[v*cdim+0], coords[v*cdim+1],coords[v*cdim+2]); 290 } 291 } 292 } 293 if (numEdges + numQuads != newVertices) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of new vertices %D != %D previous count", numEdges + numQuads, newVertices); 294 CHKMEMQ; 295 296 /* Get cell vertices by traversing loops */ 297 numQuads = 0; 298 cOff = 0; 299 for (b = 0; b < nbodies; ++b) { 300 ego body = bodies[b]; 301 int id, Nl, l; 302 303 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 304 for (l = 0; l < Nl; ++l) { 305 ego loop = lobjs[l]; 306 int lid, Ner = 0, Ne, e, nc = 0, c, Nt, t; 307 308 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 309 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d \n", lid);CHKERRQ(ierr); 310 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 311 312 for (e = 0; e < Ne; ++e) { 313 ego edge = objs[e]; 314 int points[3]; 315 int eid, Nv, v, tmp; 316 317 eid = EG_indexBodyTopo(body, edge); 318 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d \n", eid);CHKERRQ(ierr); 319 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 320 if (mtype == DEGENERATE) {ierr = PetscPrintf(PETSC_COMM_SELF, " EGDE %d is DEGENERATE \n", eid);CHKERRQ(ierr); continue;} 321 else {++Ner;} 322 if (Nv != 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices != 2", eid, Nv); 323 324 for (v = 0; v < Nv; ++v) { 325 ego vertex = nobjs[v]; 326 327 id = EG_indexBodyTopo(body, vertex); 328 points[v*2] = id-1; 329 } 330 { 331 PetscInt edgeNum; 332 333 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 334 points[1] = numVertices - newVertices + edgeNum; 335 } 336 /* EGADS loops are not oriented, but seem to be in order, so we must piece them together */ 337 if (!nc) { 338 for (v = 0; v < Nv+1; ++v) cone[nc++] = points[v]; 339 } else { 340 if (cone[nc-1] == points[0]) {cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];} 341 else if (cone[nc-1] == points[2]) {cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];} 342 else if (cone[nc-3] == points[0]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[2]) cone[nc++] = points[2];} 343 else if (cone[nc-3] == points[2]) {tmp = cone[nc-3]; cone[nc-3] = cone[nc-1]; cone[nc-1] = tmp; cone[nc++] = points[1]; if (cone[0] != points[0]) cone[nc++] = points[0];} 344 else SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Edge %d does not match its predecessor", eid); 345 } 346 } 347 if (nc != 2*Ner) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D != %D", nc, 2*Ner); 348 if (Ner == 4) {cone[nc++] = numVertices - newVertices + numEdges + numQuads++;} 349 if (nc > maxCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of corners %D > %D max", nc, maxCorners); 350 /* Triangulate the loop */ 351 switch (Ner) { 352 case 2: /* Bi-Segment -> 2 triangles */ 353 Nt = 2; 354 cells[cOff*numCorners+0] = cone[0]; 355 cells[cOff*numCorners+1] = cone[1]; 356 cells[cOff*numCorners+2] = cone[2]; 357 ++cOff; 358 cells[cOff*numCorners+0] = cone[0]; 359 cells[cOff*numCorners+1] = cone[2]; 360 cells[cOff*numCorners+2] = cone[3]; 361 ++cOff; 362 break; 363 case 3: /* Triangle -> 4 triangles */ 364 Nt = 4; 365 cells[cOff*numCorners+0] = cone[0]; 366 cells[cOff*numCorners+1] = cone[1]; 367 cells[cOff*numCorners+2] = cone[5]; 368 ++cOff; 369 cells[cOff*numCorners+0] = cone[1]; 370 cells[cOff*numCorners+1] = cone[2]; 371 cells[cOff*numCorners+2] = cone[3]; 372 ++cOff; 373 cells[cOff*numCorners+0] = cone[5]; 374 cells[cOff*numCorners+1] = cone[3]; 375 cells[cOff*numCorners+2] = cone[4]; 376 ++cOff; 377 cells[cOff*numCorners+0] = cone[1]; 378 cells[cOff*numCorners+1] = cone[3]; 379 cells[cOff*numCorners+2] = cone[5]; 380 ++cOff; 381 break; 382 case 4: /* Quad -> 8 triangles */ 383 Nt = 8; 384 cells[cOff*numCorners+0] = cone[0]; 385 cells[cOff*numCorners+1] = cone[1]; 386 cells[cOff*numCorners+2] = cone[7]; 387 ++cOff; 388 cells[cOff*numCorners+0] = cone[1]; 389 cells[cOff*numCorners+1] = cone[2]; 390 cells[cOff*numCorners+2] = cone[3]; 391 ++cOff; 392 cells[cOff*numCorners+0] = cone[3]; 393 cells[cOff*numCorners+1] = cone[4]; 394 cells[cOff*numCorners+2] = cone[5]; 395 ++cOff; 396 cells[cOff*numCorners+0] = cone[5]; 397 cells[cOff*numCorners+1] = cone[6]; 398 cells[cOff*numCorners+2] = cone[7]; 399 ++cOff; 400 cells[cOff*numCorners+0] = cone[8]; 401 cells[cOff*numCorners+1] = cone[1]; 402 cells[cOff*numCorners+2] = cone[3]; 403 ++cOff; 404 cells[cOff*numCorners+0] = cone[8]; 405 cells[cOff*numCorners+1] = cone[3]; 406 cells[cOff*numCorners+2] = cone[5]; 407 ++cOff; 408 cells[cOff*numCorners+0] = cone[8]; 409 cells[cOff*numCorners+1] = cone[5]; 410 cells[cOff*numCorners+2] = cone[7]; 411 ++cOff; 412 cells[cOff*numCorners+0] = cone[8]; 413 cells[cOff*numCorners+1] = cone[7]; 414 cells[cOff*numCorners+2] = cone[1]; 415 ++cOff; 416 break; 417 default: SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d edges, which we do not support", lid, Ner); 418 } 419 for (t = 0; t < Nt; ++t) { 420 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs Triangle %D (", t); 421 for (c = 0; c < numCorners; ++c) { 422 if (c > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");} 423 ierr = PetscPrintf(PETSC_COMM_SELF, "%D", cells[(cOff-Nt+t)*numCorners+c]); 424 } 425 ierr = PetscPrintf(PETSC_COMM_SELF, ")\n"); 426 } 427 } 428 EG_free(lobjs); 429 } 430 } 431 if (cOff != numCells) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Count of total cells %D != %D previous count", cOff, numCells); 432 ierr = DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);CHKERRQ(ierr); 433 ierr = PetscFree3(coords, cells, cone);CHKERRQ(ierr); 434 /* Embed EGADS model in DM */ 435 { 436 PetscContainer modelObj; 437 438 ierr = PetscContainerCreate(PETSC_COMM_SELF, &modelObj);CHKERRQ(ierr); 439 ierr = PetscContainerSetPointer(modelObj, model);CHKERRQ(ierr); 440 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 441 ierr = PetscContainerDestroy(&modelObj);CHKERRQ(ierr); 442 } 443 /* Label points */ 444 ierr = DMCreateLabel(dm, "EGADS Body ID");CHKERRQ(ierr); 445 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 446 ierr = DMCreateLabel(dm, "EGADS Face ID");CHKERRQ(ierr); 447 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 448 ierr = DMCreateLabel(dm, "EGADS Edge ID");CHKERRQ(ierr); 449 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 450 cOff = 0; 451 for (b = 0; b < nbodies; ++b) { 452 ego body = bodies[b]; 453 int id, Nl, l; 454 455 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 456 for (l = 0; l < Nl; ++l) { 457 ego loop = lobjs[l]; 458 ego *fobjs; 459 int lid, Nf, fid, Ner = 0, Ne, e, Nt = 0, t; 460 461 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 462 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 463 if (Nf > 1) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Loop %d has %d > 1 faces, which is not supported", lid, Nf);CHKERRQ(ierr); 464 fid = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 465 EG_free(fobjs); 466 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 467 for (e = 0; e < Ne; ++e) { 468 ego edge = objs[e]; 469 int eid, Nv, v; 470 PetscInt points[3], support[2], numEdges, edgeNum; 471 const PetscInt *edges; 472 473 eid = EG_indexBodyTopo(body, edge); 474 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 475 if (mtype == DEGENERATE) {continue;} 476 else {++Ner;} 477 for (v = 0; v < Nv; ++v) { 478 ego vertex = nobjs[v]; 479 480 id = EG_indexBodyTopo(body, vertex); 481 ierr = DMLabelSetValue(edgeLabel, numCells + id-1, eid);CHKERRQ(ierr); 482 points[v*2] = numCells + id-1; 483 } 484 ierr = PetscHMapIGet(edgeMap, eid-1, &edgeNum);CHKERRQ(ierr); 485 points[1] = numCells + numVertices - newVertices + edgeNum; 486 487 ierr = DMLabelSetValue(edgeLabel, points[1], eid);CHKERRQ(ierr); 488 support[0] = points[0]; 489 support[1] = points[1]; 490 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 491 if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges); 492 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 493 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 494 support[0] = points[1]; 495 support[1] = points[2]; 496 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 497 if (numEdges != 1) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Vertices (%D, %D) should only bound 1 edge, not %D", support[0], support[1], numEdges); 498 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 499 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 500 } 501 switch (Ner) { 502 case 2: Nt = 2;break; 503 case 3: Nt = 4;break; 504 case 4: Nt = 8;break; 505 default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Loop with %d edges is unsupported", Ner); 506 } 507 for (t = 0; t < Nt; ++t) { 508 ierr = DMLabelSetValue(bodyLabel, cOff+t, b);CHKERRQ(ierr); 509 ierr = DMLabelSetValue(faceLabel, cOff+t, fid);CHKERRQ(ierr); 510 } 511 cOff += Nt; 512 } 513 EG_free(lobjs); 514 } 515 ierr = PetscHMapIDestroy(&edgeMap);CHKERRQ(ierr); 516 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 517 for (c = cStart; c < cEnd; ++c) { 518 PetscInt *closure = NULL; 519 PetscInt clSize, cl, bval, fval; 520 521 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 522 ierr = DMLabelGetValue(bodyLabel, c, &bval);CHKERRQ(ierr); 523 ierr = DMLabelGetValue(faceLabel, c, &fval);CHKERRQ(ierr); 524 for (cl = 0; cl < clSize*2; cl += 2) { 525 ierr = DMLabelSetValue(bodyLabel, closure[cl], bval);CHKERRQ(ierr); 526 ierr = DMLabelSetValue(faceLabel, closure[cl], fval);CHKERRQ(ierr); 527 } 528 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 529 } 530 ierr = DMLabelView(bodyLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 531 ierr = DMLabelView(faceLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 532 ierr = DMLabelView(edgeLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 533 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 534 535 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 536 ierr = DMDestroy(&dm);CHKERRQ(ierr); 537 538 /* Close EGADSlite file */ 539 ierr = EG_close(context);CHKERRQ(ierr); 540 ierr = PetscFinalize(); 541 return ierr; 542 } 543 544 /*TEST 545 546 build: 547 requires: egads 548 549 test: 550 suffix: sphere_0 551 filter: sed "s/DM_[0-9a-zA-Z]*_0/DM__0/g" 552 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_refine 1 -dm_view -dm_plex_check_all 553 554 test: 555 suffix: nozzle_0 556 filter: sed "s/DM_[0-9a-zA-Z]*_0/DM__0/g" 557 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/nozzle.egadslite -dm_refine 1 -dm_view -dm_plex_check_all 558 559 TEST*/ 560