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