Constitutive model number passed

Author: dhyan1272

src/pmpo_MPMesh.cpp

@@ -61,7 +61,7 @@ void MPMesh::calculateStrain(){
   p_MPs->parallel_for(setMPStrainRate, "setMPStrainRate");
 }
 
-void MPMesh::calculateStress(){
+void MPMesh::calculateStress(const int constitutive_relation){
   //MeshFields  
   auto solveStress = p_mesh->getMeshField<polyMPO::MeshF_SolveStress>();
   auto elasticTimeStep = p_mesh->getElasticTimeStep();
@@ -75,16 +75,14 @@ void MPMesh::calculateStress(){
   auto MPsIcePressure = p_MPs->getData<polyMPO::MPF_IcePressure>();
   auto MPsRepPressure = p_MPs->getData<polyMPO::MPF_ReplacementPressure>();
 
-  int model_no = 2; //TODO get from MPAS
-  
   auto setMPStress = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
     if(mask){
       Vec3d strain_rate (MPsStrainRate(mp, 0), MPsStrainRate(mp, 1), MPsStrainRate(mp, 2));
       Vec3d stress(MPsStress(mp, 0), MPsStress(mp, 1), MPsStress(mp, 2));
-      if (model_no == 1)
-        constitutive_linear(strain_rate, stress);
-      else if(model_no == 2)
+      if (constitutive_relation == 1)
         constitutive_evp(strain_rate, stress, MPsIcePressure(mp, 0), MPsRepPressure(mp, 0), MPsArea(mp, 0), elasticTimeStep, dampingTimescale);
+      else if(constitutive_relation == 3)
+	constitutive_linear(strain_rate, stress);
       for (int m=0 ; m<3; m++)
         MPsStress(mp, m) = stress[m];
     }

src/pmpo_MPMesh.hpp

@@ -89,7 +89,7 @@ class MPMesh{
 
     void printVTP_mesh(int printVTPIndex);
     void calculateStrain();
-    void calculateStress();
+    void calculateStress(const int constitutive_relation);
     void calculateStressDivergence();
 };
 

src/pmpo_c.cpp

@@ -614,10 +614,10 @@ void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh){
   mpMesh->calculateStrain(); 
 }
 
-void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh){
+void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_model){
   checkMPMeshValid(p_mpmesh);
   auto mpMesh = ((polyMPO::MPMesh*)p_mpmesh);
-  mpMesh->calculateStress(); 
+  mpMesh->calculateStress(constitutive_model); 
 }
 
 void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn) {

src/pmpo_c.h

@@ -47,7 +47,7 @@ void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs
 void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn);
 void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpVelHost);
 void polympo_calculateMPStrainRate_f(MPMesh_ptr p_mpmesh);
-void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh);
+void polympo_calculateMPStress_f(MPMesh_ptr p_mpmesh, const int constitutive_model);
 void polympo_setMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpStressIn);
 void polympo_getMPStress_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpStressHost);
 void polympo_setAreaMP_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* areaMPHost);

src/pmpo_fortran.f90

@@ -371,10 +371,11 @@ subroutine polympo_calculateMPStrainRate(mpMesh) &
   end subroutine
 
   !MP Stress
-  subroutine polympo_calculateMPStress(mpMesh) &
+  subroutine polympo_calculateMPStress(mpMesh, constitutive_model) &
              bind(C, NAME='polympo_calculateMPStress_f')
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
+    integer(c_int), value :: constitutive_model
   end subroutine
 
   subroutine polympo_setMPStress(mpMesh, nComps, numMPs, array) &

src/pmpo_mesh.cpp

@@ -154,7 +154,6 @@ namespace polyMPO{
   }
 
   void Mesh::gridSolveGPU(){
-    std::cout<<__FUNCTION__<<std::endl;
     //Mesh Fields
     int numVerticesOwned = getNumVerticesOwned(); 
     auto totalMassVtx = getMeshField<MeshF_TotalMassVtx>();