1 static const char help[] = "Test of EGADSLite CAD functionality"; 2 3 #include <petscdmplex.h> 4 5 #include <egads.h> 6 #include <petsc.h> 7 8 typedef struct { 9 char filename[PETSC_MAX_PATH_LEN]; 10 } AppCtx; 11 12 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 13 { 14 PetscErrorCode ierr; 15 16 PetscFunctionBeginUser; 17 options->filename[0] = '\0'; 18 19 ierr = PetscOptionsBegin(comm, "", "EGADSPlex Problem Options", "EGADSLite");CHKERRQ(ierr); 20 ierr = PetscOptionsString("-filename", "The EGADSLite file", "ex9.c", options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL);CHKERRQ(ierr); 21 ierr = PetscOptionsEnd(); 22 PetscFunctionReturn(0); 23 } 24 25 int main(int argc, char *argv[]) 26 { 27 DMLabel bodyLabel, faceLabel, edgeLabel; 28 PetscInt cStart, cEnd, c; 29 /* EGADSLite variables */ 30 ego context, model, geom, *bodies, *objs, *nobjs, *mobjs, *lobjs; 31 int oclass, mtype, nbodies, *senses; 32 int b; 33 /* PETSc variables */ 34 DM dm; 35 PetscInt dim = -1, cdim = -1, numCorners = 0, numVertices = 0, numCells = 0; 36 PetscInt *cells = NULL; 37 PetscReal *coords = NULL; 38 MPI_Comm comm; 39 PetscMPIInt rank; 40 AppCtx ctx; 41 PetscErrorCode ierr; 42 43 ierr = PetscInitialize(&argc, &argv, NULL, help); if (ierr) return ierr; 44 comm = PETSC_COMM_WORLD; 45 ierr = ProcessOptions(comm, &ctx);CHKERRQ(ierr); 46 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 47 if (!rank) { 48 /* Open EGADs file and load EGADs model data */ 49 ierr = EG_open(&context);CHKERRQ(ierr); 50 ierr = EG_loadModel(context, 0, ctx.filename, &model);CHKERRQ(ierr); 51 52 /* test bodyTopo functions */ 53 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);CHKERRQ(ierr); 54 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of BODIES (nbodies): %d \n", nbodies);CHKERRQ(ierr); 55 56 for (b = 0; b < nbodies; ++b) { 57 ego body = bodies[b]; 58 int id, Nsh, Nf, Nl, l, Ne, e, Nv, v; 59 60 /* Output Basic Model Topology */ 61 ierr = EG_getBodyTopos(body, NULL, SHELL, &Nsh, &objs);CHKERRQ(ierr); 62 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of SHELLS: %d \n", Nsh);CHKERRQ(ierr); 63 64 ierr = EG_getBodyTopos(body, NULL, FACE, &Nf, &objs);CHKERRQ(ierr); 65 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of FACES: %d \n", Nf);CHKERRQ(ierr); 66 67 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 68 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of LOOPS: %d \n", Nl);CHKERRQ(ierr); 69 70 ierr = EG_getBodyTopos(body, NULL, EDGE, &Ne, &objs);CHKERRQ(ierr); 71 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of EDGES: %d \n", Ne);CHKERRQ(ierr); 72 73 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &objs);CHKERRQ(ierr); 74 ierr = PetscPrintf(PETSC_COMM_SELF, " Number of NODES: %d \n", Nv);CHKERRQ(ierr); 75 76 for (l = 0; l < Nl; ++l) { 77 ego loop = lobjs[l]; 78 79 id = EG_indexBodyTopo(body, loop); 80 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d\n", id);CHKERRQ(ierr); 81 82 /* Get EDGE info which associated with the current LOOP */ 83 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 84 85 for (e = 0; e < Ne; ++e) { 86 ego edge = objs[e]; 87 88 id = EG_indexBodyTopo(body, edge);CHKERRQ(ierr); 89 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d\n", id);CHKERRQ(ierr); 90 91 double range[4] = {0., 0., 0., 0.}; 92 double point[3] = {0., 0., 0.}; 93 int peri; 94 95 ierr = EG_getRange(objs[e], range, &peri); 96 ierr = PetscPrintf(PETSC_COMM_SELF, " Range = %lf, %lf, %lf, %lf \n", range[0], range[1], range[2], range[3]); 97 98 /* Get NODE info which associated with the current EDGE */ 99 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses);CHKERRQ(ierr); 100 101 for (v = 0; v < Nv; ++v) { 102 ego vertex = nobjs[v]; 103 double limits[4]; 104 int dummy; 105 106 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 107 id = EG_indexBodyTopo(body, vertex); 108 ierr = PetscPrintf(PETSC_COMM_SELF, " NODE ID: %d \n", id);CHKERRQ(ierr); 109 ierr = PetscPrintf(PETSC_COMM_SELF, " (x, y, z) = (%lf, %lf, %lf) \n", limits[0], limits[1], limits[2]); 110 111 point[0] = point[0] + limits[0]; 112 point[1] = point[1] + limits[1]; 113 point[2] = point[2] + limits[2]; 114 } 115 } 116 } 117 } 118 119 /* --------------------------------------------------------------------------------------------------- 120 Generate Petsc Plex 121 Get all Nodes in model, record coordinates in a correctly formatted array 122 Cycle through bodies, cycle through loops, recorde NODE IDs in a correctly formatted array */ 123 124 /* Calculate cell and vertex sizes */ 125 ierr = EG_getTopology(model, &geom, &oclass, &mtype, NULL, &nbodies, &bodies, &senses);CHKERRQ(ierr); 126 numCells = 0; 127 numVertices = 0; 128 for (b = 0; b < nbodies; ++b) { 129 ego body = bodies[b]; 130 int id, Nl, l, Nv, v; 131 132 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 133 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 134 for (l = 0; l < Nl; ++l) { 135 ego loop = lobjs[l]; 136 137 id = EG_indexBodyTopo(body, loop); 138 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 139 numCells = PetscMax(id, numCells); 140 } 141 for (v = 0; v < Nv; ++v) { 142 ego vertex = nobjs[v]; 143 144 id = EG_indexBodyTopo(body, vertex); 145 /* TODO: Instead of assuming contiguous ids, we could use a hash table */ 146 numVertices = PetscMax(id, numVertices); 147 } 148 } 149 ierr = PetscPrintf(PETSC_COMM_SELF, "\nPLEX Input Array Checkouts\n");CHKERRQ(ierr); 150 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Cells = %d \n", numCells);CHKERRQ(ierr); 151 ierr = PetscPrintf(PETSC_COMM_SELF, " Total Number of Unique Vertices = %d \n", numVertices);CHKERRQ(ierr); 152 153 dim = 2; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 154 cdim = 3; /* Assume 3D Models :: Need to update to handle 2D Models in the future */ 155 numCorners = 3; /* TODO Check number of cell corners from EGADSLite */ 156 ierr = PetscMalloc2(numVertices*cdim, &coords, numCells*numCorners, &cells);CHKERRQ(ierr); 157 158 /* Get vertex coordinates */ 159 for (b = 0; b < nbodies; ++b) { 160 ego body = bodies[b]; 161 int id, Nv, v; 162 163 ierr = EG_getBodyTopos(body, NULL, NODE, &Nv, &nobjs);CHKERRQ(ierr); 164 for (v = 0; v < Nv; ++v) { 165 ego vertex = nobjs[v]; 166 double limits[4]; 167 int dummy; 168 169 ierr = EG_getTopology(vertex, &geom, &oclass, &mtype, limits, &dummy, &mobjs, &senses);CHKERRQ(ierr); 170 id = EG_indexBodyTopo(body, vertex);CHKERRQ(ierr); 171 coords[(id-1)*cdim+0] = limits[0]; 172 coords[(id-1)*cdim+1] = limits[1]; 173 coords[(id-1)*cdim+2] = limits[2]; 174 ierr = PetscPrintf(PETSC_COMM_SELF, " Node ID = %d \n", id); 175 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]); 176 } 177 } 178 179 /* Get cell vertices by traversing loops */ 180 for (b = 0; b < nbodies; ++b) { 181 ego body = bodies[b]; 182 int id, Nl, l; 183 184 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 185 for (l = 0; l < Nl; ++l) { 186 ego loop = lobjs[l]; 187 int lid, Ne, e, nc = 0, c; 188 189 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 190 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP ID: %d \n", lid);CHKERRQ(ierr); 191 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 192 193 for (e = 0; e < Ne; ++e) { 194 ego edge = objs[e]; 195 int Nv, v; 196 197 id = EG_indexBodyTopo(body, edge); 198 ierr = PetscPrintf(PETSC_COMM_SELF, " EDGE ID: %d \n", id);CHKERRQ(ierr); 199 if (mtype == DEGENERATE) {ierr = PetscPrintf(PETSC_COMM_SELF, " EGDE %d is DEGENERATE \n", id);CHKERRQ(ierr);} 200 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 201 202 /* Add unique vertices to cells, this handles mtype == DEGENERATE fine */ 203 for (v = 0; v < Nv; ++v) { 204 ego vertex = nobjs[v]; 205 206 id = EG_indexBodyTopo(body, vertex); 207 for (c = 0; c < nc; ++c) if (cells[(lid-1)*numCorners+c] == id-1) break; 208 if (c == nc) cells[(lid-1)*numCorners+nc++] = id-1; 209 } 210 } 211 if (nc != numCorners) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid number of cell corners %D, should be %D", nc, numCorners); 212 ierr = PetscPrintf(PETSC_COMM_SELF, " LOOP Corner NODEs ("); 213 for (c = 0; c < numCorners; ++c) { 214 if (c > 0) {ierr = PetscPrintf(PETSC_COMM_SELF, ", ");} 215 ierr = PetscPrintf(PETSC_COMM_SELF, "%D", cells[(lid-1)*numCorners+c]); 216 } 217 ierr = PetscPrintf(PETSC_COMM_SELF, ")\n"); 218 } 219 } 220 } 221 ierr = DMPlexCreateFromCellList(PETSC_COMM_WORLD, dim, numCells, numVertices, numCorners, PETSC_TRUE, cells, cdim, coords, &dm);CHKERRQ(ierr); 222 ierr = PetscFree2(coords, cells);CHKERRQ(ierr); 223 { 224 PetscContainer modelObj; 225 226 ierr = PetscContainerCreate(PETSC_COMM_SELF, &modelObj);CHKERRQ(ierr); 227 ierr = PetscContainerSetPointer(modelObj, model);CHKERRQ(ierr); 228 ierr = PetscObjectCompose((PetscObject) dm, "EGADS Model", (PetscObject) modelObj);CHKERRQ(ierr); 229 ierr = PetscContainerDestroy(&modelObj);CHKERRQ(ierr); 230 } 231 ierr = DMCreateLabel(dm, "EGADS Body ID");CHKERRQ(ierr); 232 ierr = DMGetLabel(dm, "EGADS Body ID", &bodyLabel);CHKERRQ(ierr); 233 ierr = DMCreateLabel(dm, "EGADS Face ID");CHKERRQ(ierr); 234 ierr = DMGetLabel(dm, "EGADS Face ID", &faceLabel);CHKERRQ(ierr); 235 ierr = DMCreateLabel(dm, "EGADS Edge ID");CHKERRQ(ierr); 236 ierr = DMGetLabel(dm, "EGADS Edge ID", &edgeLabel);CHKERRQ(ierr); 237 for (b = 0; b < nbodies; ++b) { 238 ego body = bodies[b]; 239 int id, Nl, l; 240 241 ierr = EG_getBodyTopos(body, NULL, LOOP, &Nl, &lobjs);CHKERRQ(ierr); 242 for (l = 0; l < Nl; ++l) { 243 ego loop = lobjs[l]; 244 int lid, cell, Ne, e; 245 246 lid = EG_indexBodyTopo(body, loop);CHKERRQ(ierr); 247 cell = lid-1; 248 ierr = DMLabelSetValue(bodyLabel, cell, b);CHKERRQ(ierr); 249 { 250 ego *fobjs; 251 int Nf, fid; 252 253 ierr = EG_getBodyTopos(body, loop, FACE, &Nf, &fobjs);CHKERRQ(ierr); 254 fid = EG_indexBodyTopo(body, fobjs[0]);CHKERRQ(ierr); 255 ierr = DMLabelSetValue(faceLabel, cell, fid);CHKERRQ(ierr); 256 } 257 258 ierr = EG_getTopology(loop, &geom, &oclass, &mtype, NULL, &Ne, &objs, &senses);CHKERRQ(ierr); 259 for (e = 0; e < Ne; ++e) { 260 ego edge = objs[e]; 261 int eid, Nv, v; 262 PetscInt support[2], numEdges; 263 const PetscInt *edges; 264 265 eid = EG_indexBodyTopo(body, edge); 266 ierr = EG_getTopology(edge, &geom, &oclass, &mtype, NULL, &Nv, &nobjs, &senses); 267 if (Nv > 2) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Edge %d has %d vertices > 2", eid, Nv); 268 for (v = 0; v < Nv; ++v) { 269 ego vertex = nobjs[v]; 270 271 id = EG_indexBodyTopo(body, vertex); 272 ierr = DMLabelSetValue(edgeLabel, numCells + id-1, eid);CHKERRQ(ierr); 273 support[v] = numCells + id-1; 274 } 275 if (Nv == 2) { 276 ierr = DMPlexGetJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 277 if (numEdges != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "2 vertices should only bound 1 edge, not %D", numEdges); 278 ierr = DMLabelSetValue(edgeLabel, edges[0], eid);CHKERRQ(ierr); 279 ierr = DMPlexRestoreJoin(dm, 2, support, &numEdges, &edges);CHKERRQ(ierr); 280 } 281 } 282 } 283 } 284 ierr = DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);CHKERRQ(ierr); 285 for (c = cStart; c < cEnd; ++c) { 286 PetscInt *closure = NULL; 287 PetscInt clSize, cl, bval, fval; 288 289 ierr = DMPlexGetTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 290 ierr = DMLabelGetValue(bodyLabel, c, &bval);CHKERRQ(ierr); 291 ierr = DMLabelGetValue(faceLabel, c, &fval);CHKERRQ(ierr); 292 for (cl = 0; cl < clSize*2; cl += 2) { 293 ierr = DMLabelSetValue(bodyLabel, closure[cl], bval);CHKERRQ(ierr); 294 ierr = DMLabelSetValue(faceLabel, closure[cl], fval);CHKERRQ(ierr); 295 } 296 ierr = DMPlexRestoreTransitiveClosure(dm, c, PETSC_TRUE, &clSize, &closure);CHKERRQ(ierr); 297 } 298 ierr = DMLabelView(bodyLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 299 ierr = DMLabelView(faceLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 300 ierr = DMLabelView(edgeLabel, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 301 ierr = DMSetFromOptions(dm);CHKERRQ(ierr); 302 303 ierr = DMViewFromOptions(dm, NULL, "-dm_view");CHKERRQ(ierr); 304 ierr = DMDestroy(&dm);CHKERRQ(ierr); 305 306 /* Close EGADSlite file */ 307 ierr = EG_close(context);CHKERRQ(ierr); 308 ierr = PetscFinalize(); 309 return ierr; 310 } 311 312 /*TEST 313 314 build: 315 requires: egads 316 317 test: 318 suffix: sphere_0 319 filter: sed "s/DM_[0-9a-zA-Z]*_0/DM__0/g" 320 args: -filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/unit_sphere.egadslite -dm_view ::ascii_info_detail 321 322 TEST*/ 323