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, ¢erOfGravity, &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, ¢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: }
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*/