Actual source code: ex18.c
1: #include "petscsys.h"
2: static const char help[] = "Test of PETSc/CAD Shape Modification Technology";
4: #include <petscdmplexegads.h>
5: #include <petsc/private/hashmapi.h>
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;
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));
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: }
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: }
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: }
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: }
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;
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));
89: PetscCall(PetscObjectSetName((PetscObject)dm, "Original Surface"));
90: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
91: PetscCall(surfArea(dm));
93: // Expose Geometry Definition Data and Calculate Surface Gradients
94: PetscCall(DMPlexGeomDataAndGrads(dm, PETSC_FALSE));
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));
100: // Get attached EGADS model (pointer)
101: PetscCall(PetscContainerGetPointer(modelObj, (void **)&model));
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));
111: // Get attached EGADS model Control Point and Weights Hash Tables and Data Arrays (pointer)
112: PetscCall(PetscContainerGetPointer(cpHashTableObj, (void **)&cpHashTable));
113: PetscCall(VecGetArrayWrite(cntrlPtCoordsVec, &cpCoordData));
114: PetscCall(PetscContainerGetPointer(cpCoordDataLengthObj, (void **)&cpCoordDataLengthPtr));
115: PetscCall(PetscContainerGetPointer(wHashTableObj, (void **)&wHashTable));
116: PetscCall(VecGetArrayWrite(cntrlPtWeightsVec, &wData));
117: PetscCall(PetscContainerGetPointer(wDataLengthObj, (void **)&wDataLengthPtr));
119: cpCoordDataLength = *cpCoordDataLengthPtr;
120: wDataLength = *wDataLengthPtr;
121: PetscCheck(cpCoordDataLength && wDataLength, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Data sizes must be positive");
123: // Get the number of bodies and body objects in the model
124: PetscCall(DMPlexGetGeomModelBodies(dm, &bodies, &Nb));
126: // Get all FACES of the Body
127: PetscGeom body = bodies[0];
128: PetscCall(DMPlexGetGeomModelBodyFaces(dm, body, &fobjs, &Nf));
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;
135: PetscCall(DMPlexGetGeomID(dm, body, face, &id));
136: PetscCall(DMPlexGetGeomFaceNumOfControlPoints(dm, face, &numCntrlPnts));
138: // Update Control Points
139: PetscHashIter CPiter, Witer;
140: PetscBool CPfound, Wfound;
141: PetscInt faceCPStartRow, faceWStartRow;
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));
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));
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));
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));
173: // Inflate Mesh to Geometry
174: PetscCall(DMPlexInflateToGeomModel(dm, PETSC_FALSE));
175: PetscCall(surfArea(dm));
177: // Output .hdf5 file
178: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
180: for (PetscInt r = 0; r < Nr; ++r) {
181: char name[PETSC_MAX_PATH_LEN];
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: }
191: PetscCall(DMDestroy(&dm));
192: PetscCall(PetscFinalize());
193: return 0;
194: }
196: /*TEST
197: build:
198: requires: egads
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
208: TEST*/