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