Strain rate, Stress and StressDivergence Calculations in polyMPO

src/CMakeLists.txt

@@ -7,6 +7,7 @@ set(HEADERS
   pmpo_MPMesh_assembly.hpp
   pmpo_c.h
   pmpo_createTestMPMesh.hpp
+  pmpo_const_relation.hpp
 )
 
 set(SOURCES
@@ -53,4 +54,4 @@ add_library(polyMPO INTERFACE)
 target_link_libraries(polyMPO INTERFACE ${polyMPO_EXPORTED_TARGETS})
 bob_export_target(polyMPO)
 
-bob_end_subdir()
+bob_end_subdir()

src/pmpo_MPMesh.hpp

@@ -18,30 +18,36 @@ template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OnSurfVeloI
 #define maxMPsPerElm 8
 
 class MPMesh{
-  private:
-
-    bool isPreComputed;
-
+
   public:
-
-    MPMesh() : isPreComputed(false){};
-    void computeMatricesAndSolve(); 
-    void resetPreComputeFlag();
-    Kokkos::View<double*[vec4d_nEntries]> precomputedVtxCoeffs;
 
     Mesh* p_mesh;
     MaterialPoints* p_MPs;
 
-    std::map<MeshFieldIndex, std::function<void()>> reconstructSlice = std::map<MeshFieldIndex, std::function<void()>>();
-
+    //For MPI Communication
+    int numOwnersTot, numHalosTot;
+    std::vector<int> numOwnersOnOtherProcs;
+    std::vector<int> numHalosOnOtherProcs;
+    std::vector<int>haloOwnerProcs;
+    std::vector<std::vector<int>> haloOwnerLocalIDs;
+    std::vector<std::vector<int>> ownerOwnerLocalIDs;
+    std::vector<std::vector<int>> ownerHaloLocalIDs;
+
+    void startCommunication();
+    void communicateFields(const std::vector<std::vector<double>>& fieldData, const int numEntities, const int numEntries, int mode,
+                           std::vector<std::vector<int>>& recvIDVec,  std::vector<std::vector<double>>& recvDataVec);
+    void communicate_and_take_halo_contributions(const Kokkos::View<double**>& meshField, int nEntities, int numEntries, int mode, int op);
+
     MPMesh(Mesh* inMesh, MaterialPoints* inMPs):
-        p_mesh(inMesh), p_MPs(inMPs) {
+      p_mesh(inMesh), p_MPs(inMPs) {
     };
+
     ~MPMesh() {
       delete p_mesh;
       delete p_MPs;
     }
 
+    //MP advection and tracking
     void CVTTrackingEdgeCenterBased(Vec2dView dx);
     void CVTTrackingElmCenterBased(const int printVTPIndex = -1);
     void T2LTracking(Vec2dView dx);
@@ -50,38 +56,41 @@ class MPMesh{
     void push_swap();
     void push_swap_pos();
     void push();
+
+    //Used before advection to interpolate fields from mesh to MPs
+    //And also before reconstruction
     void calcBasis();
 
+    //Reconstruction
     DoubleView assemblyV0();
-    template <MaterialPointSlice index>
-    DoubleView wtScaAssembly();
-    template <MaterialPointSlice index>
-    Vec2dView wtVec2Assembly();
-    template <MeshFieldIndex meshFieldIndex>
-    void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag);
     template <MeshFieldIndex meshFieldIndex>
     void assemblyVtx0();
     template <MeshFieldIndex meshFieldIndex>
     void assemblyElm0();
     template <MeshFieldIndex meshFieldIndex>
     void assemblyVtx1();
-    template <MeshFieldIndex meshFieldIndex>
-    void subAssemblyVtx1(int vtxPerElm, int nCells, int comp, double* array);
-
-    void subAssemblyCoeffs(int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14, 
-                                                           double* m22, double* m23, double* m24, 
-                                                           double* m33, double* m34, 
-                                                           double* m44);
-    void solveMatrixAndRegularize(int nVertices, double* m11, double* m12, double* m13, double* m14, 
-                                                      double* m22, double* m23, double* m24, 
-                                                      double* m33, double* m34,
-                                                      double* m44);
+    void reconstruct_coeff_full();
+    void invertMatrix(const Kokkos::View<double**>& vtxMatrices, const double& radius);
+    Kokkos::View<double*[vec3d_nEntries][vec4d_nEntries]> precomputedVtxCoeffs_new;
+    Kokkos::View<double*> nearAnEdge;   
+    Kokkos::View<double*> vtxMatrixMass;
 
+    //Not used currently
+    std::map<MeshFieldIndex, std::function<void()>> reconstructSlice = std::map<MeshFieldIndex, std::function<void()>>();
+    template <MaterialPointSlice index>
+    DoubleView wtScaAssembly();
+    template <MaterialPointSlice index>
+    Vec2dView wtVec2Assembly();
+    template <MeshFieldIndex meshFieldIndex>
+    void assembly(int order, MeshFieldType type, bool basisWeightFlag, bool massWeightFlag);
     template<MeshFieldIndex meshFieldIndex>
     void setReconstructSlice(int order, MeshFieldType type);
     void reconstructSlices();
 
     void printVTP_mesh(int printVTPIndex);
+    void calculateStrain();
+    void calculateStress();
+    void calculateStressDivergence();
 };
 
 }//namespace polyMPO end