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
surfArea(DM dm)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
ProcessOptions(MPI_Comm comm,AppCtx * options)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
main(int argc,char * argv[])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