1 static const char help[] = "Test of PETSc CAD Shape Optimization & Mesh Modification Technology"; 2 3 #include <petscdmplexegads.h> 4 #include <petsc/private/hashmapi.h> 5 6 typedef struct { 7 char filename[PETSC_MAX_PATH_LEN]; 8 PetscInt saMaxIter; // Maximum number of iterations of shape optimization loop 9 } AppCtx; 10 11 static PetscErrorCode surfArea(DM dm) 12 { 13 DMLabel bodyLabel, faceLabel; 14 double surfaceArea = 0., volume = 0.; 15 PetscReal vol, centroid[3], normal[3]; 16 PetscInt dim, cStart, cEnd, fStart, fEnd; 17 PetscInt bodyID, faceID; 18 MPI_Comm comm; 19 const char *name; 20 21 PetscFunctionBeginUser; 22 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm)); 23 PetscCall(PetscObjectGetName((PetscObject)dm, &name)); 24 PetscCall(DMGetDimension(dm, &dim)); 25 PetscCall(PetscPrintf(comm, " dim = %" PetscInt_FMT "\n", dim)); 26 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel)); 27 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel)); 28 29 if (dim == 2) { 30 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 31 for (PetscInt ii = cStart; ii < cEnd; ++ii) { 32 PetscCall(DMLabelGetValue(faceLabel, ii, &faceID)); 33 if (faceID >= 0) { 34 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal)); 35 surfaceArea += vol; 36 } 37 } 38 } 39 40 if (dim == 3) { 41 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd)); 42 for (PetscInt ii = fStart; ii < fEnd; ++ii) { 43 PetscCall(DMLabelGetValue(faceLabel, ii, &faceID)); 44 if (faceID >= 0) { 45 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal)); 46 surfaceArea += vol; 47 } 48 } 49 50 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd)); 51 for (PetscInt ii = cStart; ii < cEnd; ++ii) { 52 PetscCall(DMLabelGetValue(bodyLabel, ii, &bodyID)); 53 if (bodyID >= 0) { 54 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal)); 55 volume += vol; 56 } 57 } 58 } 59 60 if (dim == 2) { 61 PetscCall(PetscPrintf(comm, "%s Surface Area = %.6e \n\n", name, (double)surfaceArea)); 62 } else if (dim == 3) { 63 PetscCall(PetscPrintf(comm, "%s Volume = %.6e \n", name, (double)volume)); 64 PetscCall(PetscPrintf(comm, "%s Surface Area = %.6e \n\n", name, (double)surfaceArea)); 65 } 66 PetscFunctionReturn(PETSC_SUCCESS); 67 } 68 69 static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options) 70 { 71 PetscFunctionBeginUser; 72 options->filename[0] = '\0'; 73 options->saMaxIter = 200; 74 75 PetscOptionsBegin(comm, "", "DMPlex w/CAD Options", "DMPlex w/CAD"); 76 PetscCall(PetscOptionsString("-filename", "The CAD/Geometry file", __FILE__, options->filename, options->filename, PETSC_MAX_PATH_LEN, NULL)); 77 PetscCall(PetscOptionsBoundedInt("-sa_max_iter", "The maximum number of iterates for the shape optimization loop", __FILE__, options->saMaxIter, &options->saMaxIter, NULL, 0)); 78 PetscOptionsEnd(); 79 PetscFunctionReturn(0); 80 } 81 82 int main(int argc, char *argv[]) 83 { 84 /* EGADS/EGADSlite variables */ 85 PetscGeom model, *bodies, *fobjs; 86 int Nb, Nf, id; 87 /* PETSc variables */ 88 DM dmNozzle = NULL; 89 MPI_Comm comm; 90 AppCtx ctx; 91 PetscScalar equivR = 0.0; 92 const char baseName[] = "Nozzle_Mesh"; 93 94 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 95 comm = PETSC_COMM_WORLD; 96 PetscCall(ProcessOptions(comm, &ctx)); 97 98 PetscCall(DMPlexCreateFromFile(comm, ctx.filename, "EGADS", PETSC_TRUE, &dmNozzle)); 99 PetscCall(PetscObjectSetName((PetscObject)dmNozzle, baseName)); 100 //PetscCall(DMCreate(PETSC_COMM_WORLD, &dmNozzle)); 101 //PetscCall(DMSetType(dmNozzle, DMPLEX)); 102 //PetscCall(DMPlexDistributeSetDefault(dmNozzle, PETSC_FALSE)); 103 //PetscCall(DMSetFromOptions(dmNozzle)); 104 105 // Get Common Viewer to store all Mesh Results 106 PetscViewer viewer; 107 PetscViewerFormat format; 108 PetscBool flg; 109 PetscInt num = 0; 110 PetscViewerType viewType; 111 PetscViewerFormat viewFormat; 112 PetscCall(PetscOptionsCreateViewer(PETSC_COMM_WORLD, NULL, NULL, "-dm_view_test", &viewer, &format, &flg)); 113 //PetscCall(PetscViewerPushFormat(viewer, format)); 114 if (flg) { 115 PetscCall(PetscPrintf(PETSC_COMM_SELF, " flg = TRUE \n")); 116 PetscCall(PetscViewerGetType(viewer, &viewType)); 117 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_HDF5_PETSC)); // PetscOptionsGetViewer returns &format as PETSC_VIEWER_DEFAULT need PETSC_VIEWER_HDF5_PETSC to save multiple DMPlexes in a single .h5 file. 118 PetscCall(PetscViewerGetFormat(viewer, &viewFormat)); 119 PetscCall(PetscPrintf(PETSC_COMM_SELF, " viewer type = %s \n", viewType)); 120 PetscCall(PetscPrintf(PETSC_COMM_SELF, " viewer format = %d \n", viewFormat)); 121 } 122 123 // Refines Surface Mesh per option -dm_refine 124 PetscCall(DMSetFromOptions(dmNozzle)); 125 PetscCall(DMViewFromOptions(dmNozzle, NULL, "-dm_view")); 126 //PetscCall(DMGetOutputSequenceNumber(dmNozzle, &num, NULL)); 127 //PetscCall(DMSetOutputSequenceNumber(dmNozzle, num, -1)); 128 num += 1; 129 PetscCall(DMView(dmNozzle, viewer)); 130 PetscCall(surfArea(dmNozzle)); 131 132 for (PetscInt saloop = 0, seqNum = 0; saloop < ctx.saMaxIter; ++saloop) { 133 PetscContainer modelObj; 134 //PetscContainer gradSACPObj, gradSAWObj; 135 PetscScalar *cpCoordData, *wData, *gradSACP, *gradSAW; 136 PetscInt cpCoordDataLength, wDataLength, maxNumEquiv; 137 PetscHMapI cpHashTable, wHashTable; 138 Mat cpEquiv; 139 char stpName[PETSC_MAX_PATH_LEN]; 140 char meshName[PETSC_MAX_PATH_LEN]; 141 142 PetscCall(PetscSNPrintf(meshName, PETSC_MAX_PATH_LEN, "%s_Step_%" PetscInt_FMT, baseName, saloop)); 143 PetscCall(PetscObjectSetName((PetscObject)dmNozzle, meshName)); 144 // Save Step File of Updated Geometry at designated iterations 145 if (saloop == 1 || saloop == 2 || saloop == 5 || saloop == 20 || saloop == 50 || saloop == 100 || saloop == 150 || saloop == 200 || saloop == 300 || saloop == 400 || saloop == 500) { 146 PetscCall(PetscSNPrintf(stpName, sizeof(stpName), "newGeom_clean_%d.stp", saloop)); 147 } 148 149 if (saloop == 0) PetscCall(DMSetFromOptions(dmNozzle)); 150 151 // Expose Geometry Definition Data and Calculate Surface Gradients 152 PetscCall(DMPlexGeomDataAndGrads(dmNozzle, PETSC_FALSE)); 153 154 PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "EGADS Model", (PetscObject *)&modelObj)); 155 if (!modelObj) PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "EGADSlite Model", (PetscObject *)&modelObj)); 156 157 // Get attached EGADS model (pointer) 158 PetscCall(PetscContainerGetPointer(modelObj, &model)); 159 160 // Look to see if DM has Container for Geometry Control Point Data 161 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Hash Table", (PetscObject *)&cpHashTableObj)); 162 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Coordinates", (PetscObject *)&cpCoordDataObj)); 163 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Coordinate Data Length", (PetscObject *)&cpCoordDataLengthObj)); 164 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Weights Hash Table", (PetscObject *)&wHashTableObj)); 165 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Weight Data", (PetscObject *)&wDataObj)); 166 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Weight Data Length", (PetscObject *)&wDataLengthObj)); 167 ////PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Surface Area Control Point Gradient", (PetscObject *)&gradSACPObj)); 168 ////PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Surface Area Weights Gradient", (PetscObject *)&gradSAWObj)); 169 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Control Point Equivalency Matrix", (PetscObject *)&cpEquivObj)); 170 //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Maximum Number Control Point Equivalency", (PetscObject *)&maxNumRelateObj)); 171 172 // Get attached EGADS model Control Point and Weights Hash Tables and Data Arrays (pointer) 173 //PetscCall(PetscContainerGetPointer(cpHashTableObj, &cpHashTable)); 174 //PetscCall(PetscContainerGetPointer(cpCoordDataObj, &cpCoordData)); 175 //PetscCall(PetscContainerGetPointer(cpCoordDataLengthObj, &cpCoordDataLengthPtr)); 176 //PetscCall(PetscContainerGetPointer(wHashTableObj, &wHashTable)); 177 //PetscCall(PetscContainerGetPointer(wDataObj, &wData)); 178 //PetscCall(PetscContainerGetPointer(wDataLengthObj, &wDataLengthPtr)); 179 ////PetscCall(PetscContainerGetPointer(gradSACPObj, &gradSACP)); 180 ////PetscCall(PetscContainerGetPointer(gradSAWObj, &gradSAW)); 181 //PetscCall(PetscContainerGetPointer(cpEquivObj, &cpEquiv)); 182 //PetscCall(PetscContainerGetPointer(maxNumRelateObj, &maxNumRelatePtr)); 183 184 // Trying out new Function 185 PetscInt cpArraySize, wArraySize; 186 PetscHMapI cpSurfGradHT; 187 Mat cpSurfGrad; 188 PetscScalar *gradVolCP, *gradVolW; 189 190 PetscCall(DMPlexGetGeomCntrlPntAndWeightData(dmNozzle, &cpHashTable, &cpCoordDataLength, &cpCoordData, &maxNumEquiv, &cpEquiv, &wHashTable, &wDataLength, &wData)); 191 PetscCall(DMPlexGetGeomGradData(dmNozzle, &cpSurfGradHT, &cpSurfGrad, &cpArraySize, &gradSACP, &gradVolCP, &wArraySize, &gradSAW, &gradVolW)); 192 193 // Get the number of bodies and body objects in the model 194 PetscCall(DMPlexGetGeomModelBodies(dmNozzle, &bodies, &Nb)); 195 196 // Get all FACES of the Body 197 PetscGeom body = bodies[0]; 198 PetscCall(DMPlexGetGeomModelBodyFaces(dmNozzle, body, &fobjs, &Nf)); 199 200 // Print out Surface Area and Volume using EGADS' Function 201 PetscScalar volume, surfaceArea; 202 PetscInt COGsize, IMCOGsize; 203 PetscScalar *centerOfGravity, *inertiaMatrixCOG; 204 PetscCall(DMPlexGetGeomBodyMassProperties(dmNozzle, body, &volume, &surfaceArea, ¢erOfGravity, &COGsize, &inertiaMatrixCOG, &IMCOGsize)); 205 206 // Calculate Radius of Sphere for Equivalent Volume 207 if (saloop == 0) equivR = PetscPowReal(volume * (3. / 4.) / PETSC_PI, 1. / 3.); 208 209 // Get Size of Control Point Equivalency Matrix 210 PetscInt numRows, numCols; 211 PetscCall(MatGetSize(cpEquiv, &numRows, &numCols)); 212 213 // Get max number of relations 214 PetscInt maxRelateNew = 0; 215 for (PetscInt ii = 0; ii < numRows; ++ii) { 216 PetscInt maxRelateNewTemp = 0; 217 for (PetscInt jj = 0; jj < numCols; ++jj) { 218 PetscScalar matValue; 219 PetscCall(MatGetValue(cpEquiv, ii, jj, &matValue)); 220 221 if (matValue > 0.0) maxRelateNewTemp += 1; 222 } 223 if (maxRelateNewTemp > maxRelateNew) maxRelateNew = maxRelateNewTemp; 224 } 225 226 // Update Control Point and Weight definitions for each surface 227 for (PetscInt jj = 0; jj < Nf; ++jj) { 228 PetscGeom face = fobjs[jj]; 229 PetscInt numCntrlPnts = 0; 230 PetscCall(DMPlexGetGeomID(dmNozzle, body, face, &id)); 231 PetscCall(DMPlexGetGeomFaceNumOfControlPoints(dmNozzle, face, &numCntrlPnts)); 232 233 // Update Control Points 234 PetscHashIter CPiter, Witer; 235 PetscBool CPfound, Wfound; 236 PetscInt faceCPStartRow, faceWStartRow; 237 238 PetscCall(PetscHMapIFind(cpHashTable, id, &CPiter, &CPfound)); 239 PetscCheck(CPfound, comm, PETSC_ERR_SUP, "FACE ID not found in Control Point Hash Table"); 240 PetscCall(PetscHMapIGet(cpHashTable, id, &faceCPStartRow)); 241 242 PetscCall(PetscHMapIFind(wHashTable, id, &Witer, &Wfound)); 243 PetscCheck(Wfound, comm, PETSC_ERR_SUP, "FACE ID not found in Control Point Weights Hash Table"); 244 PetscCall(PetscHMapIGet(wHashTable, id, &faceWStartRow)); 245 246 for (PetscInt ii = 0; ii < numCntrlPnts; ++ii) { 247 PetscReal xold, yold, zold, rold, phi, theta; 248 249 // Update Control Points - Original Code 250 xold = cpCoordData[faceCPStartRow + (3 * ii) + 0]; 251 yold = cpCoordData[faceCPStartRow + (3 * ii) + 1]; 252 zold = cpCoordData[faceCPStartRow + (3 * ii) + 2]; 253 rold = PetscSqrtReal(xold * xold + yold * yold + zold * zold); 254 phi = PetscAtan2Real(yold, xold); 255 theta = PetscAtan2Real(PetscSqrtReal(xold * xold + yold * yold), zold); 256 257 // Account for Different Weights for Control Points on Degenerate Edges (multiple control points have same location and different weights 258 // only use the largest weight 259 PetscScalar maxWeight = 0.0; 260 //PetscInt wCntr = 0; 261 for (PetscInt kk = faceWStartRow; kk < faceWStartRow + numCntrlPnts; ++kk) { 262 PetscScalar matValue; 263 PetscCall(MatGetValue(cpEquiv, faceWStartRow + ii, kk, &matValue)); 264 265 PetscScalar weight = 0.0; 266 if (matValue > 0.0) { 267 weight = wData[kk]; 268 269 if (weight > maxWeight) maxWeight = weight; 270 //wCntr += 1; 271 } 272 } 273 274 //Reduce to Constant R = 0.0254 275 PetscScalar deltaR; 276 PetscScalar localTargetR; 277 localTargetR = equivR / maxWeight; 278 deltaR = rold - localTargetR; 279 280 cpCoordData[faceCPStartRow + (3 * ii) + 0] += -0.05 * deltaR * PetscCosReal(phi) * PetscSinReal(theta); 281 cpCoordData[faceCPStartRow + (3 * ii) + 1] += -0.05 * deltaR * PetscSinReal(phi) * PetscSinReal(theta); 282 cpCoordData[faceCPStartRow + (3 * ii) + 2] += -0.05 * deltaR * PetscCosReal(theta); 283 } 284 } 285 PetscCall(DMPlexFreeGeomObject(dmNozzle, fobjs)); 286 PetscBool writeFile = PETSC_FALSE; 287 288 // Modify Control Points of Geometry 289 PetscCall(DMPlexModifyGeomModel(dmNozzle, comm, cpCoordData, wData, PETSC_FALSE, writeFile, stpName)); 290 writeFile = PETSC_FALSE; 291 292 // Get attached EGADS model (pointer) 293 PetscGeom newmodel; 294 PetscCall(PetscContainerGetPointer(modelObj, &newmodel)); 295 296 // Get the number of bodies and body objects in the model 297 PetscCall(DMPlexGetGeomModelBodies(dmNozzle, &bodies, &Nb)); 298 299 // Get all FACES of the Body 300 PetscGeom newbody = bodies[0]; 301 PetscCall(DMPlexGetGeomModelBodyFaces(dmNozzle, newbody, &fobjs, &Nf)); 302 303 // Save Step File of Updated Geometry at designated iterations 304 if (saloop == 1 || saloop == 2 || saloop == 5 || saloop == 20 || saloop == 50 || saloop == 100 || saloop == 150 || saloop == 200 || saloop == 300 || saloop == 400 || saloop == 500) writeFile = PETSC_TRUE; 305 306 // Modify Geometry and Inflate Mesh to New Geoemetry 307 PetscCall(DMPlexModifyGeomModel(dmNozzle, comm, cpCoordData, wData, PETSC_FALSE, writeFile, stpName)); 308 PetscCall(DMPlexInflateToGeomModel(dmNozzle, PETSC_TRUE)); 309 310 // Periodically Refine and Write Mesh to hdf5 file 311 if (saloop == 0) { 312 PetscCall(DMViewFromOptions(dmNozzle, NULL, "-dm_view7")); 313 //PetscCall(DMGetOutputSequenceNumber(dmNozzle, &num, NULL)); 314 PetscCall(PetscPrintf(PETSC_COMM_SELF, " num = %d \n", num)); 315 //PetscCall(DMGetOutputSequenceNumber(dmNozzle, NULL, &num)); 316 //PetscCall(DMSetOutputSequenceNumber(dmNozzle, num, saloop)); 317 num += 1; 318 PetscCall(PetscObjectSetName((PetscObject)dmNozzle, "nozzle_meshes_1")); 319 //PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_APPEND)); 320 PetscCall(DMView(dmNozzle, viewer)); 321 } 322 if (saloop == 1 || saloop == 5 || saloop == 20 || saloop == 50 || saloop == 100 || saloop == 150 || saloop == 200 || saloop == 300 || saloop == 400 || saloop == 500) { 323 PetscCall(DMViewFromOptions(dmNozzle, NULL, "-dm_view")); 324 PetscCall(DMSetOutputSequenceNumber(dmNozzle, seqNum++, saloop)); 325 PetscCall(DMView(dmNozzle, viewer)); 326 if (saloop == 200 || saloop == 500) { 327 PetscCall(DMSetFromOptions(dmNozzle)); 328 PetscCall(DMViewFromOptions(dmNozzle, NULL, "-dm_view")); 329 PetscCall(DMSetOutputSequenceNumber(dmNozzle, seqNum++, saloop)); 330 PetscCall(DMView(dmNozzle, viewer)); 331 332 PetscCall(DMSetFromOptions(dmNozzle)); 333 PetscCall(DMViewFromOptions(dmNozzle, NULL, "-dm_view")); 334 PetscCall(DMSetOutputSequenceNumber(dmNozzle, seqNum++, saloop)); 335 PetscCall(DMView(dmNozzle, viewer)); 336 } 337 } 338 PetscCall(DMPlexRestoreGeomBodyMassProperties(dmNozzle, body, &volume, &surfaceArea, ¢erOfGravity, &COGsize, &inertiaMatrixCOG, &IMCOGsize)); 339 PetscCall(DMPlexRestoreGeomCntrlPntAndWeightData(dmNozzle, &cpHashTable, &cpCoordDataLength, &cpCoordData, &maxNumEquiv, &cpEquiv, &wHashTable, &wDataLength, &wData)); 340 PetscCall(DMPlexRestoreGeomGradData(dmNozzle, &cpSurfGradHT, &cpSurfGrad, &cpArraySize, &gradSACP, &gradVolCP, &wArraySize, &gradSAW, &gradVolW)); 341 PetscCall(DMPlexFreeGeomObject(dmNozzle, fobjs)); 342 } 343 PetscCall(DMDestroy(&dmNozzle)); 344 PetscCall(PetscFinalize()); 345 } 346 347 /*TEST 348 build: 349 requires: egads 350 351 # To view mesh use -dm_plex_view_hdf5_storage_version 3.1.0 -dm_view hdf5:mesh_minSA_abstract.h5:hdf5_petsc:append 352 test: 353 suffix: minSA 354 requires: datafilespath 355 temporaries: newGeom_clean_1.stp newGeom_clean_2.stp newGeom_clean_5.stp 356 args: -filename ${DATAFILESPATH}/meshes/cad/sphere_example.stp \ 357 -dm_refine 1 -dm_plex_geom_print_model 1 -dm_plex_geom_shape_opt 1 \ 358 -sa_max_iter 2 359 360 TEST*/ 361