1 #include "petscsys.h"
2 static const char help[] = "Test of PETSc/CAD Shape Modification Technology";
3
4 #include <petscdmplexegads.h>
5 #include <petsc/private/hashmapi.h>
6
surfArea(DM dm)7 static PetscErrorCode surfArea(DM dm)
8 {
9 DMLabel bodyLabel, faceLabel;
10 double surfaceArea = 0., volume = 0.;
11 PetscReal vol, centroid[3], normal[3];
12 PetscInt dim, cStart, cEnd, fStart, fEnd;
13 PetscInt bodyID, faceID;
14 MPI_Comm comm;
15 const char *name;
16
17 PetscFunctionBeginUser;
18 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
19 PetscCall(PetscObjectGetName((PetscObject)dm, &name));
20 PetscCall(DMGetDimension(dm, &dim));
21 PetscCall(PetscPrintf(comm, " dim = %" PetscInt_FMT "\n", dim));
22 PetscCall(DMGetLabel(dm, "EGADS Body ID", &bodyLabel));
23 PetscCall(DMGetLabel(dm, "EGADS Face ID", &faceLabel));
24 PetscCall(DMGetCoordinatesLocalSetUp(dm));
25
26 if (dim == 2) {
27 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
28 for (PetscInt ii = cStart; ii < cEnd; ++ii) {
29 PetscCall(DMLabelGetValue(faceLabel, ii, &faceID));
30 if (faceID >= 0) {
31 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal));
32 surfaceArea += vol;
33 }
34 }
35 }
36
37 if (dim == 3) {
38 PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
39 for (PetscInt ii = fStart; ii < fEnd; ++ii) {
40 PetscCall(DMLabelGetValue(faceLabel, ii, &faceID));
41 if (faceID >= 0) {
42 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal));
43 surfaceArea += vol;
44 }
45 }
46
47 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
48 for (PetscInt ii = cStart; ii < cEnd; ++ii) {
49 PetscCall(DMLabelGetValue(bodyLabel, ii, &bodyID));
50 if (bodyID >= 0) {
51 PetscCall(DMPlexComputeCellGeometryFVM(dm, ii, &vol, centroid, normal));
52 volume += vol;
53 }
54 }
55 }
56
57 if (dim == 2) {
58 PetscCall(PetscPrintf(comm, "%s Surface Area = %.6e \n\n", name, (double)surfaceArea));
59 } else if (dim == 3) {
60 PetscCall(PetscPrintf(comm, "%s Volume = %.6e \n", name, (double)volume));
61 PetscCall(PetscPrintf(comm, "%s Surface Area = %.6e \n\n", name, (double)surfaceArea));
62 }
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
main(int argc,char * argv[])66 int main(int argc, char *argv[])
67 {
68 /* EGADSlite variables */
69 PetscGeom model, *bodies, *fobjs;
70 int Nb, Nf, id;
71 /* PETSc variables */
72 DM dm;
73 MPI_Comm comm;
74 PetscContainer modelObj, cpHashTableObj, wHashTableObj, cpCoordDataLengthObj, wDataLengthObj;
75 Vec cntrlPtCoordsVec, cntrlPtWeightsVec;
76 PetscScalar *cpCoordData, *wData;
77 PetscInt cpCoordDataLength, wDataLength;
78 PetscInt *cpCoordDataLengthPtr, *wDataLengthPtr;
79 PetscHMapI cpHashTable, wHashTable;
80 PetscInt Nr = 2;
81
82 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
83 comm = PETSC_COMM_WORLD;
84 PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
85 PetscCall(DMSetType(dm, DMPLEX));
86 PetscCall(DMPlexDistributeSetDefault(dm, PETSC_FALSE));
87 PetscCall(DMSetFromOptions(dm));
88
89 PetscCall(PetscObjectSetName((PetscObject)dm, "Original Surface"));
90 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
91 PetscCall(surfArea(dm));
92
93 // Expose Geometry Definition Data and Calculate Surface Gradients
94 PetscCall(DMPlexGeomDataAndGrads(dm, PETSC_FALSE));
95
96 // Get Attached EGADS model
97 PetscCall(PetscObjectQuery((PetscObject)dm, "EGADS Model", (PetscObject *)&modelObj));
98 if (!modelObj) PetscCall(PetscObjectQuery((PetscObject)dm, "EGADSlite Model", (PetscObject *)&modelObj));
99
100 // Get attached EGADS model (pointer)
101 PetscCall(PetscContainerGetPointer(modelObj, &model));
102
103 // Look to see if DM has Container for Geometry Control Point Data
104 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Hash Table", (PetscObject *)&cpHashTableObj));
105 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Coordinates", (PetscObject *)&cntrlPtCoordsVec));
106 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Coordinate Data Length", (PetscObject *)&cpCoordDataLengthObj));
107 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Weights Hash Table", (PetscObject *)&wHashTableObj));
108 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Weight Data", (PetscObject *)&cntrlPtWeightsVec));
109 PetscCall(PetscObjectQuery((PetscObject)dm, "Control Point Weight Data Length", (PetscObject *)&wDataLengthObj));
110
111 // Get attached EGADS model Control Point and Weights Hash Tables and Data Arrays (pointer)
112 PetscCall(PetscContainerGetPointer(cpHashTableObj, &cpHashTable));
113 PetscCall(VecGetArrayWrite(cntrlPtCoordsVec, &cpCoordData));
114 PetscCall(PetscContainerGetPointer(cpCoordDataLengthObj, &cpCoordDataLengthPtr));
115 PetscCall(PetscContainerGetPointer(wHashTableObj, &wHashTable));
116 PetscCall(VecGetArrayWrite(cntrlPtWeightsVec, &wData));
117 PetscCall(PetscContainerGetPointer(wDataLengthObj, &wDataLengthPtr));
118
119 cpCoordDataLength = *cpCoordDataLengthPtr;
120 wDataLength = *wDataLengthPtr;
121 PetscCheck(cpCoordDataLength && wDataLength, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Data sizes must be positive");
122
123 // Get the number of bodies and body objects in the model
124 PetscCall(DMPlexGetGeomModelBodies(dm, &bodies, &Nb));
125
126 // Get all FACES of the Body
127 PetscGeom body = bodies[0];
128 PetscCall(DMPlexGetGeomModelBodyFaces(dm, body, &fobjs, &Nf));
129
130 // Update Control Point and Weight definitions for each surface
131 for (PetscInt jj = 0; jj < Nf; ++jj) {
132 PetscGeom face = fobjs[jj];
133 PetscInt numCntrlPnts = 0;
134
135 PetscCall(DMPlexGetGeomID(dm, body, face, &id));
136 PetscCall(DMPlexGetGeomFaceNumOfControlPoints(dm, face, &numCntrlPnts));
137
138 // Update Control Points
139 PetscHashIter CPiter, Witer;
140 PetscBool CPfound, Wfound;
141 PetscInt faceCPStartRow, faceWStartRow;
142
143 PetscCall(PetscHMapIFind(cpHashTable, id, &CPiter, &CPfound));
144 PetscCheck(CPfound, comm, PETSC_ERR_SUP, "FACE ID not found in Control Point Hash Table");
145 PetscCall(PetscHMapIGet(cpHashTable, id, &faceCPStartRow));
146
147 PetscCall(PetscHMapIFind(wHashTable, id, &Witer, &Wfound));
148 PetscCheck(Wfound, comm, PETSC_ERR_SUP, "FACE ID not found in Control Point Weights Hash Table");
149 PetscCall(PetscHMapIGet(wHashTable, id, &faceWStartRow));
150
151 for (PetscInt ii = 0; ii < numCntrlPnts; ++ii) {
152 if (ii == 4) {
153 // Update Control Points - Change the location of the center control point of the faces
154 // Note :: Modification geometry requires knowledge of how the geometry is defined.
155 cpCoordData[faceCPStartRow + (3 * ii) + 0] = 2.0 * cpCoordData[faceCPStartRow + (3 * ii) + 0];
156 cpCoordData[faceCPStartRow + (3 * ii) + 1] = 2.0 * cpCoordData[faceCPStartRow + (3 * ii) + 1];
157 cpCoordData[faceCPStartRow + (3 * ii) + 2] = 2.0 * cpCoordData[faceCPStartRow + (3 * ii) + 2];
158 } else {
159 // Do Nothing
160 // Note: Could use section the change location of other face control points.
161 }
162 }
163 }
164 PetscCall(DMPlexFreeGeomObject(dm, fobjs));
165
166 // Modify Control Points of Geometry
167 PetscCall(PetscObjectSetName((PetscObject)dm, "Modified Surface"));
168 // TODO Wrap EG_saveModel() in a Plex viewer to manage file access
169 PetscCall(DMPlexModifyGeomModel(dm, comm, cpCoordData, wData, PETSC_FALSE, PETSC_TRUE, "newModel.stp"));
170 PetscCall(VecRestoreArrayWrite(cntrlPtCoordsVec, &cpCoordData));
171 PetscCall(VecRestoreArrayWrite(cntrlPtWeightsVec, &wData));
172
173 // Inflate Mesh to Geometry
174 PetscCall(DMPlexInflateToGeomModel(dm, PETSC_FALSE));
175 PetscCall(surfArea(dm));
176
177 // Output .hdf5 file
178 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
179
180 for (PetscInt r = 0; r < Nr; ++r) {
181 char name[PETSC_MAX_PATH_LEN];
182
183 // Perform refinement on the Mesh attached to the new geometry
184 PetscCall(DMSetFromOptions(dm));
185 PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN, "Modified Surface refinement %" PetscInt_FMT, r));
186 PetscCall(PetscObjectSetName((PetscObject)dm, name));
187 PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
188 PetscCall(surfArea(dm));
189 }
190
191 PetscCall(DMDestroy(&dm));
192 PetscCall(PetscFinalize());
193 return 0;
194 }
195
196 /*TEST
197 build:
198 requires: egads
199
200 # Use -dm_view hdf5:mesh_shapeMod_sphere.h5 to view the mesh
201 test:
202 suffix: sphere_shapeMod
203 requires: datafilespath
204 temporaries: newModel.stp
205 args: -dm_plex_filename ${DATAFILESPATH}/meshes/cad/sphere_example.stp \
206 -dm_refine -dm_plex_geom_print_model 1 -dm_plex_geom_shape_opt 1
207
208 TEST*/
209