diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index f4d57aa1..a4e97fed 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -699,6 +699,49 @@ void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
   Kokkos::deep_copy(arrayHost, array_d);
 }
 
+void polympo_setMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, const double** forceArray){
+  //forceArray has 6 entries, each entry must have nVertices values
+  //check mpMesh is valid
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+
+  auto strainRate = p_mesh->getMeshField<polyMPO::MeshF_VtxStrainRate>();
+  auto h_strainRate = Kokkos::create_mirror_view(strainRate);
+  for(int i = 0; i < nVertices; i++){
+    h_strainRate(i,0) = forceArray[0][i];
+    h_strainRate(i,1) = forceArray[1][i];
+    h_strainRate(i,2) = forceArray[2][i];
+    h_strainRate(i,3) = forceArray[3][i];
+    h_strainRate(i,4) = forceArray[4][i];
+    h_strainRate(i,5) = forceArray[5][i];
+  }
+  Kokkos::deep_copy(strainRate, h_strainRate);
+}
+
+void polympo_getMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, double** forceArray){
+  //check mpMesh is valid
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+  
+  auto strainRate = p_mesh->getMeshField<polyMPO::MeshF_VtxStrainRate>();
+  auto h_strainRate = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), strainRate);
+  for(int i = 0; i < nVertices; i++){
+    forceArray[0][i] = h_strainRate(i,0);
+    forceArray[1][i] = h_strainRate(i,1);
+    forceArray[2][i] = h_strainRate(i,2);
+    forceArray[3][i] = h_strainRate(i,3);
+    forceArray[4][i] = h_strainRate(i,4);
+    forceArray[5][i] = h_strainRate(i,5);
+  }
+}
+
+
 void polympo_push_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   ((polyMPO::MPMesh*)p_mpmesh) ->push();
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index c310c05c..b5314034 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -62,6 +62,8 @@ void polympo_setMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, cons
 void polympo_getMeshOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 void polympo_setMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 void polympo_getMeshOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
+void polympo_setMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, const double** forceArray);
+void polympo_getMeshVtxStrainRate_f(MPMesh_ptr p_mpmesh, const int nVertices, double** forceArray);
 
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);
diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90
index d1d14a35..b49c0202 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -537,6 +537,30 @@ subroutine polympo_getMeshOnSurfDispIncr(mpMesh, nComps, nVertices, array) &
     type(c_ptr), value :: array
   end subroutine
   !---------------------------------------------------------------------------
+  !> @brief set the verte strain rate from a host array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVertices(in) numVertices
+  !> @param forceArray(in) input mesh strain rate 2D array (6,numVtx)
+  subroutine polympo_setMeshVtxStrainRate(mpMesh, nVertices, forceArray) &
+             bind(C, NAME='polympo_setMeshVtxStrainRate_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: forceArray
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the verte strain rate from polyMPO
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVertices(in) numVertices
+  !> @param forceArray(in/out) output mesh strain rate 2D array (6,numVtx)
+  subroutine polympo_getMeshVtxStrainRate(mpMesh, nVertices, forceArray) &
+             bind(C, NAME='polympo_getMeshVtxStrainRate_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: forceArray
+  end subroutine
+  !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices
   !>        MPs MUST have rotated flag set to True(>0)
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index 78b3a7be..2f7ae6e4 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -37,6 +37,10 @@ namespace polyMPO{
         auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr);
         PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased);
         vtxRotLatLonIncr_ = DoubleVec2dView(vtxRotLatLonIncrMapEntry.second,numVtxs_);
+
+        auto vtxStrainRateMapEntry = meshFields2TypeAndString.at(MeshF_VtxStrainRate);
+        PMT_ALWAYS_ASSERT(vtxStrainRateMapEntry.first == MeshFType_VtxBased);
+        vtxStrainRate_ = DoubleSymMat3dView(vtxStrainRateMapEntry.second,numVtxs_);
     }
     
     void Mesh::computeRotLatLonIncr(){
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index 3f94256a..c0b4df8c 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -25,7 +25,8 @@ enum MeshFieldIndex{
     MeshF_Vel,
     MeshF_OnSurfVeloIncr,
     MeshF_OnSurfDispIncr,
-    MeshF_RotLatLonIncr
+    MeshF_RotLatLonIncr,
+    MeshF_VtxStrainRate
 };
 enum MeshFieldType{
     MeshFType_Invalid = -2,
@@ -42,7 +43,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
                {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
                {MeshF_OnSurfVeloIncr,   {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
                {MeshF_OnSurfDispIncr,   {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
-               {MeshF_RotLatLonIncr,    {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}};
+               {MeshF_RotLatLonIncr,    {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
+               {MeshF_VtxStrainRate,    {MeshFType_VtxBased,"MeshField_StrainRate"}}};
 
 enum mesh_type {mesh_unrecognized_lower = -1,
                 mesh_general_polygonal, //other meshes
@@ -74,6 +76,7 @@ class Mesh {
     DoubleVec2dView vtxOnSurfVeloIncr_;
     DoubleVec2dView vtxOnSurfDispIncr_;
     DoubleVec2dView vtxRotLatLonIncr_;
+    DoubleSymMat3dView vtxStrainRate_;
     //DoubleMat2DView vtxStress_;
 
   public:
@@ -161,6 +164,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_RotLatLonIncr){
         return vtxRotLatLonIncr_;
     }
+    else if constexpr (index==MeshF_VtxStrainRate){
+        return vtxStrainRate_;
+    }
     fprintf(stderr,"Mesh Field Index error!\n");
     exit(1);
 }
diff --git a/test/testFortran.f90 b/test/testFortran.f90
index 8bf0f367..7420cfe0 100644
--- a/test/testFortran.f90
+++ b/test/testFortran.f90
@@ -29,7 +29,8 @@ program main
   integer, dimension(:), pointer :: MPElmID
   real(kind=APP_RKIND), dimension(:,:), pointer :: MParray
   real(kind=APP_RKIND), dimension(:,:), pointer :: MPPositions
-  real(kind=APP_RKIND), dimension(:,:), pointer :: Mesharray
+  real(kind=APP_RKIND), dimension(:,:), pointer :: mesharray_vel
+  real(kind=APP_RKIND), dimension(:,:), pointer :: mesharray_strainrate
   real(kind=APP_RKIND), dimension(:), pointer :: xArray, yArray, zArray
   integer :: ierr, self
   type(c_ptr) :: mpMesh
@@ -52,10 +53,11 @@ program main
   numMPs = 49 !todo use getNumMPs from the MaterialPoints object
   numElms = 10
 
-  allocate(Mesharray(numCompsVel,nverts))
+  allocate(mesharray_vel(numCompsVel,nverts))
   allocate(MParray(numCompsVel,numMPs))
   allocate(MPElmID(numMPs))
   allocate(MPPositions(numCompsCoords,numMPs))
+  allocate(mesharray_strainrate(6,nverts))
   allocate(xArray(nverts))
   allocate(yArray(nverts))
   allocate(zArray(nverts))
@@ -90,25 +92,25 @@ program main
   ! set mesh Fields
   do i = 1,numCompsVel
     do j = 1,nverts 
-        Mesharray(i,j) = (i-1)*nverts + j
+        mesharray_vel(i,j) = (i-1)*nverts + j
     end do
   end do
-  call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
-  call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
+  call polympo_setMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(mesharray_vel))
+  call polympo_setMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(mesharray_vel))
 
   ! check mesh Fields
-  Mesharray = -1
-  call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
+  mesharray_vel = -1
+  call polympo_getMeshOnSurfVeloIncr(mpMesh, numCompsVel, nverts, c_loc(mesharray_vel))
   do i = 1,numCompsVel
     do j = 1,nverts 
-        call assert((Mesharray(i,j) .eq. (i-1)*nverts+j), "Assert MeshOnSurfVeloIncr Fail")
+        call assert((mesharray_vel(i,j) .eq. (i-1)*nverts+j), "Assert MeshOnSurfVeloIncr Fail")
     end do
   end do
-  Mesharray = -1
-  call polympo_getMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(Mesharray))
+  mesharray_vel = -1
+  call polympo_getMeshOnSurfDispIncr(mpMesh, numCompsVel, nverts, c_loc(mesharray_vel))
   do i = 1,numCompsVel
     do j = 1,nverts 
-        call assert((Mesharray(i,j) .eq. (i-1)*nverts+j), "Assert MeshOnSurfDispIncr Fail")
+        call assert((mesharray_vel(i,j) .eq. (i-1)*nverts+j), "Assert MeshOnSurfDispIncr Fail")
     end do
   end do
 
@@ -139,12 +141,28 @@ program main
     call assert((xArray(i) .eq. i+value1), "Assert MeshVel u-component Velocity Fail")
     call assert((yArray(i) .eq. value2-i), "Assert MeshVel v-component Velocity Fail")
   end do 
+  
+  !VtxStrainRate needs 6 arrays, using mesharray_strainrate to hold them
+  do i = 1, 6
+    do j = 1, nverts
+      mesharray_strainrate(i,j) = (i + j) * ((2 * i) + (3 * j))
+    end do
+  end do
+  call polympo_setMeshVtxStrainRate(mpMesh, nverts, c_loc(mesharray_strainrate))
+  mesharray_strainrate = -1
+  call polympo_getMeshVtxStrainRate(mpMesh, nverts, c_loc(mesharray_strainrate))
+  do i = 1, 6
+    do j = 1, nverts
+      call assert((mesharray_strainrate(i,j) .eq. (i + j) * ((2 * i) + (3 * j))), "Assert MeshVtxStrainRate Fail")
+    end do
+  end do
 
   deallocate(MParray)
-  deallocate(Mesharray)
+  deallocate(mesharray_vel)
   deallocate(xArray)
   deallocate(yArray)
   deallocate(zArray)
+  deallocate(mesharray_strainrate)
 
   call polympo_deleteMPMesh(mpMesh)
   call polympo_finalize()