Actual source code: ex19.c

  1: static const char help[] = "Test of PETSc CAD Shape Optimization & Mesh Modification Technology";

  3: #include <petscdmplexegads.h>
  4: #include <petsc/private/hashmapi.h>

  6: typedef struct {
  7:   char     filename[PETSC_MAX_PATH_LEN];
  8:   PetscInt saMaxIter; // Maximum number of iterations of shape optimization loop
  9: } AppCtx;

 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;

 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));

 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:   }

 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:     }

 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:   }

 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: }

 69: static PetscErrorCode ProcessOptions(MPI_Comm comm, AppCtx *options)
 70: {
 71:   PetscFunctionBeginUser;
 72:   options->filename[0] = '\0';
 73:   options->saMaxIter   = 200;

 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: }

 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";

 94:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 95:   comm = PETSC_COMM_WORLD;
 96:   PetscCall(ProcessOptions(comm, &ctx));

 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));

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:   }

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));

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];

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:     }

149:     if (saloop == 0) PetscCall(DMSetFromOptions(dmNozzle));

151:     // Expose Geometry Definition Data and Calculate Surface Gradients
152:     PetscCall(DMPlexGeomDataAndGrads(dmNozzle, PETSC_FALSE));

154:     PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "EGADS Model", (PetscObject *)&modelObj));
155:     if (!modelObj) PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "EGADSlite Model", (PetscObject *)&modelObj));

157:     // Get attached EGADS model (pointer)
158:     PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));

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 Equivalancy Matrix", (PetscObject *)&cpEquivObj));
170:     //PetscCall(PetscObjectQuery((PetscObject)dmNozzle, "Maximum Number Control Point Equivalency", (PetscObject *)&maxNumRelateObj));

172:     // Get attached EGADS model Control Point and Weights Hash Tables and Data Arrays (pointer)
173:     //PetscCall(PetscContainerGetPointer(cpHashTableObj, (void **)&cpHashTable));
174:     //PetscCall(PetscContainerGetPointer(cpCoordDataObj, (void **)&cpCoordData));
175:     //PetscCall(PetscContainerGetPointer(cpCoordDataLengthObj, (void **)&cpCoordDataLengthPtr));
176:     //PetscCall(PetscContainerGetPointer(wHashTableObj, (void **)&wHashTable));
177:     //PetscCall(PetscContainerGetPointer(wDataObj, (void **)&wData));
178:     //PetscCall(PetscContainerGetPointer(wDataLengthObj, (void **)&wDataLengthPtr));
179:     ////PetscCall(PetscContainerGetPointer(gradSACPObj, (void **)&gradSACP));
180:     ////PetscCall(PetscContainerGetPointer(gradSAWObj, (void **)&gradSAW));
181:     //PetscCall(PetscContainerGetPointer(cpEquivObj, (void **)&cpEquiv));
182:     //PetscCall(PetscContainerGetPointer(maxNumRelateObj, (void **)&maxNumRelatePtr));

184:     // Trying out new Function
185:     PetscInt     cpArraySize, wArraySize;
186:     PetscHMapI   cpSurfGradHT;
187:     Mat          cpSurfGrad;
188:     PetscScalar *gradVolCP, *gradVolW;

190:     PetscCall(DMPlexGetGeomCntrlPntAndWeightData(dmNozzle, &cpHashTable, &cpCoordDataLength, &cpCoordData, &maxNumEquiv, &cpEquiv, &wHashTable, &wDataLength, &wData));
191:     PetscCall(DMPlexGetGeomGradData(dmNozzle, &cpSurfGradHT, &cpSurfGrad, &cpArraySize, &gradSACP, &gradVolCP, &wArraySize, &gradSAW, &gradVolW));

193:     // Get the number of bodies and body objects in the model
194:     PetscCall(DMPlexGetGeomModelBodies(dmNozzle, &bodies, &Nb));

196:     // Get all FACES of the Body
197:     PetscGeom body = bodies[0];
198:     PetscCall(DMPlexGetGeomModelBodyFaces(dmNozzle, body, &fobjs, &Nf));

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, &centerOfGravity, &COGsize, &inertiaMatrixCOG, &IMCOGsize));

206:     // Calculate Radius of Sphere for Equivalent Volume
207:     if (saloop == 0) equivR = PetscPowReal(volume * (3. / 4.) / PETSC_PI, 1. / 3.);

209:     // Get Size of Control Point Equivalancy Matrix
210:     PetscInt numRows, numCols;
211:     PetscCall(MatGetSize(cpEquiv, &numRows, &numCols));

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));

221:         if (matValue > 0.0) maxRelateNewTemp += 1;
222:       }
223:       if (maxRelateNewTemp > maxRelateNew) maxRelateNew = maxRelateNewTemp;
224:     }

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));

233:       // Update Control Points
234:       PetscHashIter CPiter, Witer;
235:       PetscBool     CPfound, Wfound;
236:       PetscInt      faceCPStartRow, faceWStartRow;

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));

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));

246:       for (PetscInt ii = 0; ii < numCntrlPnts; ++ii) {
247:         PetscReal xold, yold, zold, rold, phi, theta;

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);

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));

265:           PetscScalar weight = 0.0;
266:           if (matValue > 0.0) {
267:             weight = wData[kk];

269:             if (weight > maxWeight) { maxWeight = weight; }
270:             //wCntr += 1;
271:           }
272:         }

274:         //Reduce to Constant R = 0.0254
275:         PetscScalar deltaR;
276:         PetscScalar localTargetR;
277:         localTargetR = equivR / maxWeight;
278:         deltaR       = rold - localTargetR;

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;

288:     // Modify Control Points of Geometry
289:     PetscCall(DMPlexModifyGeomModel(dmNozzle, comm, cpCoordData, wData, PETSC_FALSE, writeFile, stpName));
290:     writeFile = PETSC_FALSE;

292:     // Get attached EGADS model (pointer)
293:     PetscGeom newmodel;
294:     PetscCall(PetscContainerGetPointer(modelObj, (void **)&newmodel));

296:     // Get the number of bodies and body objects in the model
297:     PetscCall(DMPlexGetGeomModelBodies(dmNozzle, &bodies, &Nb));

299:     // Get all FACES of the Body
300:     PetscGeom newbody = bodies[0];
301:     PetscCall(DMPlexGetGeomModelBodyFaces(dmNozzle, newbody, &fobjs, &Nf));

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; }

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));

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));

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, &centerOfGravity, &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: }

347: /*TEST
348:   build:
349:     requires: egads

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

360: TEST*/