diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index 45319be..ba02a82 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -127,7 +127,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
     auto MPs2Elm = p_MPs->getData<MPF_Tgt_Elm_ID>();
     auto MPs2Proc = p_MPs->getData<MPF_Tgt_Proc_ID>();
     auto elm2Process = p_mesh->getElm2Process();
-    
+    auto elm2global = p_mesh->getElmGlobal();
+
     if(printVTPIndex>=0) {
       printVTP_mesh(printVTPIndex);
     }
@@ -155,8 +156,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                 for(int i=1; i<=numConnElms; i++){
                     int elmID = elm2ElmConn(iElm,i)-1;
                     
-	            //New delta
-	            Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2));
+                    //New delta
+                    Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2));
                     delta = MPnew - center;
 
                     double neighborDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
@@ -175,7 +176,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                     iElm = closestElm;
                 }
             }
-            if(printVTPIndex>=0){ 
+            if(printVTPIndex>=0 && numMPs>0){
                 double d1 = dx[0];
                 double d2 = dx[2];
                 double d3 = dx[3];
@@ -301,6 +302,7 @@ void MPMesh::reconstructSlices() {
     if (reconstructSlice.size() == 0) return;
     Kokkos::Timer timer;
     calcBasis();
+    resetPreComputeFlag();
     for (auto const& [index, reconstruct] : reconstructSlice) {
         if (reconstruct) reconstruct();
     }
@@ -322,11 +324,54 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) {
   return anyIsMigrating;
 }
 
+void MPMesh::push_ahead(){
+  Kokkos::Timer timer;
+  //Latitude Longitude increment at mesh vertices and interpolate to particle position
+  p_mesh->computeRotLatLonIncr();   
+
+  //Interpolates latitude longitude increments and mesh velocity increments to
+  //MP positions
+  sphericalInterpolationDispVelIncr(*this);
+  
+  //Push the MPs
+  p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius());
+  pumipic::RecordTime("PolyMPO_interpolateAndPush", timer.seconds());
+}
+
+bool MPMesh::push1P(){
+  Kokkos::Timer timer;
+  //Given target location find the new element or the last element in a partioned mesh
+  //and the process it belongs to so that migration can be checked
+  CVTTrackingElmCenterBased(); 
+  //From the above two inputs find if any particle needs to be migrated
+  bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate());
+  pumipic::RecordTime("PolyMPO_trackAndCheckMigrate", timer.seconds());
+  return anyIsMigrating;
+}
+
+void MPMesh::push_swap(){
+  //current becomes target, target becomes -1
+  p_MPs->updateMPElmID();
+}
+
+void MPMesh::push_swap_pos(){
+  //current becomes target, target becomes -1
+  //Making read for next push_ahead  
+  p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>();
+  p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>();
+}
+
+
 void MPMesh::push(){
+  
   Kokkos::Timer timer;
+  
   p_mesh->computeRotLatLonIncr();
+  
   sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
+  
   p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius()); // set Tgt_XYZ
+  
   auto elm2Process = p_mesh->getElm2Process();
 
   bool anyIsMigrating = false;
@@ -334,10 +379,14 @@ void MPMesh::push(){
     CVTTrackingElmCenterBased(); // move to Tgt_XYZ
     p_MPs->updateMPSlice<MPF_Cur_Pos_XYZ, MPF_Tgt_Pos_XYZ>(); // Tgt_XYZ becomes Cur_XYZ
     p_MPs->updateMPSlice<MPF_Cur_Pos_Rot_Lat_Lon, MPF_Tgt_Pos_Rot_Lat_Lon>(); // Tgt becomes Cur
-    if (elm2Process.size() > 0)
-        anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->migrate());
+    
+    bool anyIsMigrating = getAnyIsMigrating(p_MPs, p_MPs->check_migrate());
+   
+    if(anyIsMigrating)
+      p_MPs->migrate();
     else
-        p_MPs->rebuild(); //rebuild pumi-pic
+      p_MPs->rebuild();
+    
     p_MPs->updateMPElmID(); //update mpElm IDs slices
     reconstructSlices(); 
   } 
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index 5a63ef9..cb5c893 100644
--- a/src/pmpo_MPMesh.hpp
+++ b/src/pmpo_MPMesh.hpp
@@ -7,16 +7,28 @@
 
 namespace polyMPO{
 
-template <MeshFieldIndex> const MaterialPointSlice meshFieldIndexToMPSlice;
+template <MeshFieldIndex>
+const MaterialPointSlice meshFieldIndexToMPSlice;
 template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_Vel            > = MPF_Vel;
 template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_VtxMass        > = MPF_Mass;
 template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_ElmMass        > = MPF_Mass;
 template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_RotLatLonIncr  > = MPF_Rot_Lat_Lon_Incr;
+template <> const MaterialPointSlice meshFieldIndexToMPSlice < MeshF_OnSurfVeloIncr > = MPF_Vel_Incr;
 
 #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;
 
@@ -33,6 +45,10 @@ class MPMesh{
     void CVTTrackingEdgeCenterBased(Vec2dView dx);
     void CVTTrackingElmCenterBased(const int printVTPIndex = -1);
     void T2LTracking(Vec2dView dx);
+    bool push1P();
+    void push_ahead();
+    void push_swap();
+    void push_swap_pos();
     void push();
     void calcBasis();
 
@@ -49,6 +65,17 @@ class MPMesh{
     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);
 
     template<MeshFieldIndex meshFieldIndex>
     void setReconstructSlice(int order, MeshFieldType type);
diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp
index e5dedb5..033386e 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -98,35 +98,37 @@ void MPMesh::assemblyElm0() {
   pumipic::RecordTime("PolyMPO_Reconstruct_Elm0", timer.seconds());
 }
 
-template <MeshFieldIndex meshFieldIndex>
-void MPMesh::assemblyVtx1() {
+
+void MPMesh::resetPreComputeFlag(){
+  isPreComputed = false;
+}
+
+//Method 1 for coefficients
+void MPMesh::computeMatricesAndSolve(){
+  Kokkos::Timer timer;
   //Mesh Information
   auto elm2VtxConn = p_mesh->getElm2VtxConn();  
   int numVtx = p_mesh->getNumVertices();
   auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
 
-  //Mesh Field
-  constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
-  const int numEntries = mpSliceToNumEntries<mpfIndex>();
-  p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
-  auto meshField = p_mesh->getMeshField<meshFieldIndex>();
+  //Dual Element Area for Regularization
+  auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
 
   //Material Points
-  auto mpData = p_MPs->getData<mpfIndex>();
   auto weight = p_MPs->getData<MPF_Basis_Vals>();
   auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
 
   //Matrix for each vertex
   Kokkos::View<double*[vec4d_nEntries][vec4d_nEntries]> VtxMatrices("VtxMatrices", p_mesh->getNumVertices());
 
-  //Reconstructed values
-  Kokkos::View<double**> reconVals("meshField", p_mesh->getNumVertices(), numEntries);
-  
   //Earth Radius
   double radius = 1.0;
   if(p_mesh->getGeomType() == geom_spherical_surf)
     radius=p_mesh->getSphereRadius();
 
+  bool scaling=true;
+  int reg_method = 2;
+  
   //Assemble matrix for each vertex
   auto assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -134,37 +136,336 @@ void MPMesh::assemblyVtx1() {
       for(int i=0; i<nVtxE; i++){
         int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
         double w_vtx=weight(mp,i);
-        
+        double mScale=1;
+        if(scaling)
+          mScale=sqrt(dual_triangle_area(vID,0))/radius;
         double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius, 
                                                 (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
                                                 (vtxCoords(vID,2) - mpPositions(mp,2))/radius};		 
-        for (int k=0; k<vec4d_nEntries; k++)
-          for (int l=0; l<vec4d_nEntries; l++)
+        //All entries except first row and column
+        for (int k=1; k<vec4d_nEntries; k++)
+          for (int l=1; l<vec4d_nEntries; l++)
             Kokkos::atomic_add(&VtxMatrices(vID,k,l), CoordDiffs[k] * CoordDiffs[l] * w_vtx);
+        //First entry
+        Kokkos::atomic_add(&VtxMatrices(vID,0,0), CoordDiffs[0] * CoordDiffs[0] * w_vtx*mScale*mScale); 
+        //First row and column except the first entry
+        for (int k=1; k<vec4d_nEntries; k++){
+          Kokkos::atomic_add(&VtxMatrices(vID,0,k), CoordDiffs[0] * CoordDiffs[k] * w_vtx*mScale);
+          Kokkos::atomic_add(&VtxMatrices(vID,k,0), CoordDiffs[k] * CoordDiffs[0] * w_vtx*mScale);
+        }
       }
     }
   };
   p_MPs->parallel_for(assemble, "assembly");
   
-  //Solve Ax=b for each vertex  
+  //Assemble matrix for each vertex and apply regularization
   Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxCoeffs", p_mesh->getNumVertices());
+
   Kokkos::parallel_for("solving Ax=b", numVtx, KOKKOS_LAMBDA(const int vtx){
     Vec4d v0 = {VtxMatrices(vtx,0,0), VtxMatrices(vtx,0,1), VtxMatrices(vtx,0,2), VtxMatrices(vtx,0,3)};
     Vec4d v1 = {VtxMatrices(vtx,1,0), VtxMatrices(vtx,1,1), VtxMatrices(vtx,1,2), VtxMatrices(vtx,1,3)};
     Vec4d v2 = {VtxMatrices(vtx,2,0), VtxMatrices(vtx,2,1), VtxMatrices(vtx,2,2), VtxMatrices(vtx,2,3)};
     Vec4d v3 = {VtxMatrices(vtx,3,0), VtxMatrices(vtx,3,1), VtxMatrices(vtx,3,2), VtxMatrices(vtx,3,3)};
-     
-    Matrix4d A = {v0,v1,v2,v3}; 
-    double A_trace = A.trace();
+    //Define the matrices
+    Matrix4d A = {v0,v1,v2,v3};
     Matrix4d A_regularized = {v0, v1, v2, v3};
-    A_regularized.addToDiag(A_trace*1e-8);
-   
+    //Regularization
+    switch(reg_method){
+      case 0:{
+        break;
+      }   
+      case 1:{
+        double A_trace = A.trace();
+        A_regularized.addToDiag(A_trace*1e-8);
+        break;
+      }
+      case 2:{
+        double regParam=sqrt(EPSILON)*(VtxMatrices(vtx,0,0)+VtxMatrices(vtx,1,1)+
+			               VtxMatrices(vtx,2,2)+VtxMatrices(vtx,3,3));
+        A_regularized.addToDiag(regParam);
+        break;
+      }
+      default:{
+        printf("Invalid regularization method \n");
+        break;	
+      }
+    }
+    //Solve Ax=b 
     double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
     CholeskySolve4d_UnitRHS(A_regularized, coeff);
+    // Undo scaling
+    double mScale=1;
+    if(scaling)
+      mScale=sqrt(dual_triangle_area(vtx,0))/radius;
+        
+    coeff[0]=coeff[0]*mScale*mScale;
+    coeff[1]=coeff[1]*mScale;
+    coeff[2]=coeff[2]*mScale;
+    coeff[3]=coeff[3]*mScale;
     for (int i=0; i<vec4d_nEntries; i++) 
       VtxCoeffs(vtx,i)=coeff[i];
   });
+  this->precomputedVtxCoeffs = VtxCoeffs;
+  pumipic::RecordTime("PolyMPO_Calculate_MLS_Coeff", timer.seconds());
+}
+
+//Method 2 for coefficients
+void MPMesh::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){
+  
+  Kokkos::Timer timer;
+  //Material Points Information
+  MPI_Comm comm = p_MPs->getMPIComm(); 
+  int comm_rank;
+  MPI_Comm_rank(comm, &comm_rank);
+
+  //Mesh Information
+  auto elm2VtxConn = p_mesh->getElm2VtxConn();  
+  int numVtx = p_mesh->getNumVertices();
+  auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
+  auto elm2Process = p_mesh->getElm2Process();
+  auto elm2global = p_mesh->getElmGlobal();
+
+  //Dual Element Area for Regularization
+  auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
   
+  //Material Points
+  calcBasis();
+  auto weight = p_MPs->getData<MPF_Basis_Vals>();
+  auto mpPos = p_MPs->getData<MPF_Cur_Pos_XYZ>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+
+  //Radius
+  double radius = 1.0;
+  if(p_mesh->getGeomType() == geom_spherical_surf)
+    radius=p_mesh->getSphereRadius();
+
+  Kokkos::View<double**> m11_d("m11", vtxPerElm, nCells);
+  Kokkos::View<double**> m12_d("m12", vtxPerElm, nCells);
+  Kokkos::View<double**> m13_d("m13", vtxPerElm, nCells);
+  Kokkos::View<double**> m14_d("m14", vtxPerElm, nCells);
+  Kokkos::View<double**> m22_d("m22", vtxPerElm, nCells);
+  Kokkos::View<double**> m23_d("m23", vtxPerElm, nCells);
+  Kokkos::View<double**> m24_d("m23", vtxPerElm, nCells);
+  Kokkos::View<double**> m33_d("m33", vtxPerElm, nCells);
+  Kokkos::View<double**> m34_d("m34", vtxPerElm, nCells);
+  Kokkos::View<double**> m44_d("m34", vtxPerElm, nCells);
+ 
+  auto sub_assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
+    if(mask && (elm2Process(elm)==comm_rank)) { //if material point is 'active'/'enabled'
+      int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
+      for(int i=0; i<nVtxE; i++){
+        int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
+          
+          double w_vtx=weight(mp,i);
+          double mScale=sqrt(dual_triangle_area(vID,0))/radius;
+          
+          Kokkos::atomic_add(&m11_d(i,elm), w_vtx*mScale*mScale);
+          Kokkos::atomic_add(&m12_d(i,elm), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))/radius);
+          Kokkos::atomic_add(&m13_d(i,elm), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))/radius);
+          Kokkos::atomic_add(&m14_d(i,elm), w_vtx*mScale*(vtxCoords(vID,2)-mpPos(mp,2))/radius);
+          Kokkos::atomic_add(&m22_d(i,elm), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,0)-mpPos(mp,0))/(radius*radius));
+          Kokkos::atomic_add(&m23_d(i,elm), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,1)-mpPos(mp,1))/(radius*radius));
+          Kokkos::atomic_add(&m24_d(i,elm), w_vtx*mScale*(vtxCoords(vID,0)-mpPos(mp,0))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
+          Kokkos::atomic_add(&m33_d(i,elm), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))*(vtxCoords(vID,1)-mpPos(mp,1))/(radius*radius));
+          Kokkos::atomic_add(&m34_d(i,elm), w_vtx*mScale*(vtxCoords(vID,1)-mpPos(mp,1))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
+          Kokkos::atomic_add(&m44_d(i,elm), w_vtx*mScale*(vtxCoords(vID,2)-mpPos(mp,2))*(vtxCoords(vID,2)-mpPos(mp,2))/(radius*radius));
+      }
+    }
+  };
+  p_MPs->parallel_for(sub_assemble, "sub_assembly");
+  pumipic::RecordTime("VtxSubAssemblyComputeCoeff", timer.seconds());
+
+  Kokkos::Timer timer2;
+  kkDbl2dViewHostU m11_h(m11, vtxPerElm, nCells);
+  kkDbl2dViewHostU m12_h(m12, vtxPerElm, nCells);
+  kkDbl2dViewHostU m13_h(m13, vtxPerElm, nCells);
+  kkDbl2dViewHostU m14_h(m14, vtxPerElm, nCells);
+  kkDbl2dViewHostU m22_h(m22, vtxPerElm, nCells);
+  kkDbl2dViewHostU m23_h(m23, vtxPerElm, nCells);
+  kkDbl2dViewHostU m24_h(m24, vtxPerElm, nCells);
+  kkDbl2dViewHostU m33_h(m33, vtxPerElm, nCells);
+  kkDbl2dViewHostU m34_h(m34, vtxPerElm, nCells);
+  kkDbl2dViewHostU m44_h(m44, vtxPerElm, nCells);
+  
+  Kokkos::deep_copy(m11_h, m11_d); 
+  Kokkos::deep_copy(m12_h, m12_d); 
+  Kokkos::deep_copy(m13_h, m13_d); 
+  Kokkos::deep_copy(m14_h, m14_d); 
+  Kokkos::deep_copy(m22_h, m22_d); 
+  Kokkos::deep_copy(m23_h, m23_d); 
+  Kokkos::deep_copy(m24_h, m24_d); 
+  Kokkos::deep_copy(m33_h, m33_d); 
+  Kokkos::deep_copy(m34_h, m34_d); 
+  Kokkos::deep_copy(m44_h, m44_d); 
+  pumipic::RecordTime("VtxSubAssemblyGetCoeff", timer2.seconds());
+  
+}
+
+//Method 2 for coefficients Solve matrix
+void MPMesh::solveMatrixAndRegularize(int nVertices, double* m11, double* m12, double* m13, double* m14, 
+                                       double* m22, double* m23, double* m24, 
+                                       double* m33, double* m34,
+                                       double* m44){
+
+  Kokkos::Timer timer;
+  kkViewHostU<const double*> m11_h(m11, nVertices);
+  kkViewHostU<const double*> m12_h(m12, nVertices);
+  kkViewHostU<const double*> m13_h(m13, nVertices);
+  kkViewHostU<const double*> m14_h(m14, nVertices);
+  kkViewHostU<const double*> m22_h(m22, nVertices);
+  kkViewHostU<const double*> m23_h(m23, nVertices);
+  kkViewHostU<const double*> m24_h(m24, nVertices);
+  kkViewHostU<const double*> m33_h(m33, nVertices);
+  kkViewHostU<const double*> m34_h(m34, nVertices);
+  kkViewHostU<const double*> m44_h(m44, nVertices);
+  
+  Kokkos::View<double*> m11_d("m11", nVertices);
+  Kokkos::View<double*> m12_d("m12", nVertices);
+  Kokkos::View<double*> m13_d("m13", nVertices);
+  Kokkos::View<double*> m14_d("m14", nVertices); 
+  Kokkos::View<double*> m22_d("m22", nVertices);
+  Kokkos::View<double*> m23_d("m23", nVertices);
+  Kokkos::View<double*> m24_d("m24", nVertices);
+  Kokkos::View<double*> m33_d("m33", nVertices);
+  Kokkos::View<double*> m34_d("m34", nVertices);
+  Kokkos::View<double*> m44_d("m44", nVertices);
+  
+  Kokkos::deep_copy(m11_d, m11_h);
+  Kokkos::deep_copy(m12_d, m12_h);
+  Kokkos::deep_copy(m13_d, m13_h);
+  Kokkos::deep_copy(m14_d, m14_h);
+  Kokkos::deep_copy(m22_d, m22_h);
+  Kokkos::deep_copy(m23_d, m23_h);
+  Kokkos::deep_copy(m24_d, m24_h);
+  Kokkos::deep_copy(m33_d, m33_h);
+  Kokkos::deep_copy(m34_d, m34_h);
+  Kokkos::deep_copy(m44_d, m44_h);
+  pumipic::RecordTime("polyMPOsolveMatrixCoeffSet", timer.seconds());
+
+  Kokkos::Timer timer2;
+  auto dual_triangle_area=p_mesh->getMeshField<MeshF_DualTriangleArea>();
+  Kokkos::View<double*[vec4d_nEntries]> VtxCoeffs("VtxCoeffs", nVertices);
+  double radius=p_mesh->getSphereRadius();
+  Kokkos::parallel_for("fill", nVertices, KOKKOS_LAMBDA(const int vtx){
+    Vec4d v0 = {m11_d(vtx), m12_d(vtx), m13_d(vtx), m14_d(vtx)};
+    Vec4d v1 = {m12_d(vtx), m22_d(vtx), m23_d(vtx), m24_d(vtx)};
+    Vec4d v2 = {m13_d(vtx), m23_d(vtx), m33_d(vtx), m34_d(vtx)};
+    Vec4d v3 = {m14_d(vtx), m24_d(vtx), m34_d(vtx), m44_d(vtx)}; 
+    //Matrix4d A = {v0,v1,v2,v3};
+    Matrix4d A_regularized = {v0, v1, v2, v3};
+    double regParam = sqrt(EPSILON)*(m11_d(vtx) + m22_d(vtx) + m33_d(vtx) + m44_d(vtx));
+    A_regularized.addToDiag(regParam);
+    
+    double coeff[vec4d_nEntries]={0.0, 0.0, 0.0, 0.0};
+    CholeskySolve4d_UnitRHS(A_regularized, coeff);
+    
+    double mScale=sqrt(dual_triangle_area(vtx,0))/radius;
+    coeff[0]=coeff[0]*mScale*mScale;
+    coeff[1]=coeff[1]*mScale;
+    coeff[2]=coeff[2]*mScale;
+    coeff[3]=coeff[3]*mScale;
+
+    for (int i=0; i<vec4d_nEntries; i++) 
+      VtxCoeffs(vtx,i)=coeff[i];
+  });
+  this->precomputedVtxCoeffs = VtxCoeffs;
+  pumipic::RecordTime("polyMPOsolveMatrixCoeffCompute", timer2.seconds());
+
+}
+
+//Method2
+template <MeshFieldIndex meshFieldIndex>
+void MPMesh::subAssemblyVtx1(int vtxPerElm, int nCells, int comp, double* array) {
+  Kokkos::Timer timer; 
+  
+  auto VtxCoeffs=this->precomputedVtxCoeffs; 
+  
+  // MPI Information
+  MPI_Comm comm = p_MPs->getMPIComm(); 
+  int comm_rank;
+  MPI_Comm_rank(comm, &comm_rank);
+
+  //Mesh Information
+  auto elm2VtxConn = p_mesh->getElm2VtxConn();  
+  int numVtx = p_mesh->getNumVertices();
+  auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
+  auto elm2Process = p_mesh->getElm2Process();
+ 
+  // Material Points Information
+  constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
+  auto mpData = p_MPs->getData<mpfIndex>();
+  auto weight = p_MPs->getData<MPF_Basis_Vals>();
+  auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
+ 
+  double radius=p_mesh->getSphereRadius();
+
+  Kokkos::View<double**> array_d("reconstructedIceArea", vtxPerElm, nCells);
+  auto sub_assemble = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
+    if(mask && (elm2Process(elm)==comm_rank)) { 
+      int nVtxE = elm2VtxConn(elm,0); //number of vertices bounding the element
+      for(int i=0; i<nVtxE; i++){
+        int vID = elm2VtxConn(elm,i+1)-1; //vID = vertex id
+        double w_vtx=weight(mp,i);
+        double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius, 
+                                                (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
+                                                (vtxCoords(vID,2) - mpPositions(mp,2))/radius};
+
+        auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + 
+                                                VtxCoeffs(vID,2)*CoordDiffs[2] + 
+                                                VtxCoeffs(vID,3)*CoordDiffs[3]);
+  
+        auto val = factor*mpData(mp, comp);
+        Kokkos::atomic_add(&array_d(i, elm), val);
+      }
+    }
+  };
+  p_MPs->parallel_for(sub_assemble, "sub_assembly"); 
+  pumipic::RecordTime("polyMPOsubAssemblyFieldCompute", timer.seconds());
+  
+  Kokkos::Timer timer2;
+  kkDbl2dViewHostU arrayHost(array, vtxPerElm, nCells);
+  Kokkos::deep_copy(arrayHost, array_d); 
+  pumipic::RecordTime("PolyMPOsubAssemblyFieldGet", timer2.seconds());
+}
+
+//Method 1
+template <MeshFieldIndex meshFieldIndex>
+void MPMesh::assemblyVtx1() {
+  Kokkos::Timer timer; 
+  //If no reconstruction till now calculate the coeffs
+  if (!isPreComputed) {
+    computeMatricesAndSolve();
+    isPreComputed=true;
+  }
+  
+  auto VtxCoeffs=this->precomputedVtxCoeffs;
+  //Mesh Information
+  auto elm2VtxConn = p_mesh->getElm2VtxConn();  
+  int numVtx = p_mesh->getNumVertices();
+  auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
+
+  //Mesh Field
+  constexpr MaterialPointSlice mpfIndex = meshFieldIndexToMPSlice<meshFieldIndex>;
+  const int numEntries = mpSliceToNumEntries<mpfIndex>();
+  p_mesh->fillMeshField<meshFieldIndex>(numVtx, numEntries, 0.0);
+  auto meshField = p_mesh->getMeshField<meshFieldIndex>();
+
+  //Material Points
+  auto mpData = p_MPs->getData<mpfIndex>();
+  auto weight = p_MPs->getData<MPF_Basis_Vals>();
+  auto mpPositions = p_MPs->getData<MPF_Cur_Pos_XYZ>();
+
+  //Earth Radius
+  double radius = 1.0;
+  if(p_mesh->getGeomType() == geom_spherical_surf)
+    radius=p_mesh->getSphereRadius();
+
+  //Reconstructed values
+  Kokkos::View<double**> reconVals("meshField", p_mesh->getNumVertices(), numEntries);
+
   //Reconstruct
   auto reconstruct = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
     if(mask) { //if material point is 'active'/'enabled'
@@ -175,7 +476,7 @@ void MPMesh::assemblyVtx1() {
         double CoordDiffs[vec4d_nEntries] = {1, (vtxCoords(vID,0) - mpPositions(mp,0))/radius,
                                                 (vtxCoords(vID,1) - mpPositions(mp,1))/radius, 
                                                 (vtxCoords(vID,2) - mpPositions(mp,2))/radius};
-        
+
         auto factor = w_vtx*(VtxCoeffs(vID,0) + VtxCoeffs(vID,1)*CoordDiffs[1] + 
                                                 VtxCoeffs(vID,2)*CoordDiffs[2] + 
                                                 VtxCoeffs(vID,3)*CoordDiffs[3]);
@@ -188,11 +489,13 @@ void MPMesh::assemblyVtx1() {
     }
   };
   p_MPs->parallel_for(reconstruct, "reconstruct");
- 
+
+  //Assign as a field 
   Kokkos::parallel_for("assigning", numVtx, KOKKOS_LAMBDA(const int vtx){
     for(int k=0; k<numEntries; k++)
       meshField(vtx, k) = reconVals(vtx,k);
   });
+  pumipic::RecordTime("PolyMPO_Reconstruct_Vtx1", timer.seconds());
 }
 
 template <MeshFieldIndex meshFieldIndex>
diff --git a/src/pmpo_c.cpp b/src/pmpo_c.cpp
index d5d0650..2369d57 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -61,41 +61,48 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm){
 
 void polympo_createMPs_f(MPMesh_ptr p_mpmesh,
                        const int numElms,
-                       const int numMPs, // total number of MPs which is GREATER than or equal to number of active MPs
+                       const int numMPs, // total number of MPs which is >= number of active MPs
                        int* mpsPerElm,
                        const int* mp2Elm,
                        const int* isMPActive) {
   checkMPMeshValid(p_mpmesh);
-
   //the mesh must be fixed/set before adding MPs
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
   PMT_ALWAYS_ASSERT(!p_mesh->meshEditable());
   PMT_ALWAYS_ASSERT(p_mesh->getNumElements() == numElms);
 
+  //Find the total no of MPs across all ranks
+  //And loop over all MPs and find the smallest element id associated across a MP
   int numActiveMPs = 0;
-  int minElmID = numElms+1;
+  int minElmID = INT_MAX;
   for(int i = 0; i < numMPs; i++) {
     if(isMPActive[i] == MP_ACTIVE) {
-      if(mp2Elm[i] < minElmID) {
+      numActiveMPs++;
+      if(mp2Elm[i] < minElmID) 
         minElmID = mp2Elm[i];
-        numActiveMPs++;
-      }
     }
   }
-  //TODO do we care about empty ranks? check just in case...
-  PMT_ALWAYS_ASSERT(numActiveMPs>0);
-
-  int firstElmWithMPs=-1;
+  int globalNumActiveMPs = 0;
+  int globalMinElmID;
+  MPI_Allreduce(&numActiveMPs, &globalNumActiveMPs, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); 
+  MPI_Allreduce(&minElmID, &globalMinElmID, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD); 
+  PMT_ALWAYS_ASSERT(globalNumActiveMPs>0);
+  
+  //Loop over all mesh elements 0,1,... and find the first element that has an associated MP
+  int firstElmWithMPs=INT_MAX;
   for (int i=0; i<numElms; i++) {
     if(mpsPerElm[i]) {
       firstElmWithMPs = i;
       break;
     }
   }
+  int globalFirstElmWithMPs;
+  MPI_Allreduce(&firstElmWithMPs, &globalFirstElmWithMPs, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD);
+  
   int offset = -1;
-  if(minElmID-firstElmWithMPs==1) {
+  if(globalMinElmID-globalFirstElmWithMPs==1) {
     offset = 1;
-  }else if (minElmID-firstElmWithMPs==0){
+  }else if (globalMinElmID-globalFirstElmWithMPs==0){
     offset = 0;
   }else {
     fprintf(stderr,"The minElmID is incorrect! Offset is wrong!\n");
@@ -112,20 +119,22 @@ void polympo_createMPs_f(MPMesh_ptr p_mpmesh,
       numActiveMPs++;
     }
   }
-
+  auto elm2global = p_mesh->getElmGlobal();
   auto mpsPerElm_d = create_mirror_view_and_copy(mpsPerElm, numElms);
   auto active_mp2Elm_d = create_mirror_view_and_copy(active_mp2Elm.data(), numActiveMPs);
   auto active_mpIDs_d = create_mirror_view_and_copy(active_mpIDs.data(), numActiveMPs);
 
   delete ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   ((polyMPO::MPMesh*)p_mpmesh)->p_MPs =
-     new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d);
+     new polyMPO::MaterialPoints(numElms, numActiveMPs, mpsPerElm_d, active_mp2Elm_d, active_mpIDs_d, elm2global);
+
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   p_MPs->setElmIDoffset(offset);
+  
 }
 
 void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh,
-                         const int numMPs, // total number of MPs which is GREATER than or equal to number of active MPs
+                         const int numMPs, // Total MPs which is GREATER than or equal to number of active MPs
                          const int* allMP2Elm,
                          const int* addedMPMask) {
   checkMPMeshValid(p_mpmesh);
@@ -186,11 +195,52 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh,
   }
 }
 
-void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh)
-{
+void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh,
+                         const int sizeMP2elm,
+                         const int* elem_ids,
+                         const int nMPs_delete,
+                         const int nMPs_add,
+                         int* recvMPs_elm,
+                         int* recvMPs_ids) {
+  
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  int offset = p_MPs->getElmIDoffset();
+
+  for (int k=0; k<nMPs_add; k++){
+    recvMPs_elm[k]=recvMPs_elm[k]-offset;
+    recvMPs_ids[k]=recvMPs_ids[k]-offset;
+  }
+  
+  auto elem_ids_d = create_mirror_view_and_copy(elem_ids, sizeMP2elm);
+  auto recvMPs_elm_d = create_mirror_view_and_copy(recvMPs_elm, nMPs_add);
+  auto recvMPs_ids_d = create_mirror_view_and_copy(recvMPs_ids, nMPs_add);
+  
+  Kokkos::View<int*> mp2Elm("mp2Elm", p_MPs->getCapacity());
+  Kokkos::View<int*> numDeletedMPs_d("numDeletedMPs", 1);
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  auto setMP2Elm = PS_LAMBDA(const int& e, const int& mp, const int& mask) {
+    if(mask) {
+      mp2Elm(mp) = elem_ids_d(mpAppID(mp))-offset;
+      if (mp2Elm(mp) == MP_DELETE){
+          Kokkos::atomic_increment(&numDeletedMPs_d(0));
+      }
+    }
+  };
+  p_MPs->parallel_for(setMP2Elm, "setMP2Elm");
+  int numDeletedMPs = pumipic::getLastValue(numDeletedMPs_d);
+  assert(nMPs_delete==numDeletedMPs);
+  p_MPs->startRebuild(mp2Elm, nMPs_add, recvMPs_elm_d, recvMPs_ids_d);
+  pumipic::RecordTime("polympo_startRebuildMPs2_f", timer.seconds());
+}
+
+void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh){
+  Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   p_MPs->finishRebuild();
+  pumipic::RecordTime("polympo_finishRebuildMPs_f", timer.seconds());
 }
 
 void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs) {
@@ -200,13 +250,39 @@ void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appI
   p_MPs->setAppIDFunc(getNextAppID);
 }
 
+void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh,
+                            const int numMPs,
+                            int* elmIDs){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  auto mpTgtElmID = p_MPs->getData<polyMPO::MPF_Tgt_Elm_ID>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  auto elmIDoffset = p_MPs->getElmIDoffset();
+
+  kkIntViewHostU arrayHost(elmIDs,numMPs);
+  polyMPO::IntView mpTgtElmIDCopy("mpTgtElmIDNewValue",numMPs);
+
+  auto setTgtElmId = PS_LAMBDA(const int&, const int& mp, const int& mask){
+    if(mask){
+        mpTgtElmIDCopy(mpAppID(mp)) = mpTgtElmID(mp)+elmIDoffset;
+    }
+  };
+  p_MPs->parallel_for(setTgtElmId, "set mpTgtElmID");
+  Kokkos::deep_copy( arrayHost, mpTgtElmIDCopy);
+  pumipic::RecordTime("PolyMPO_getMPTgtElmID", timer.seconds());
+}
+
 void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh,
                            const int numMPs,
                            int* elmIDs){
+   
   checkMPMeshValid(p_mpmesh);
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
   auto mpCurElmID = p_MPs->getData<polyMPO::MPF_Cur_Elm_ID>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
   auto elmIDoffset = p_MPs->getElmIDoffset();
@@ -221,6 +297,7 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh,
   };
   p_MPs->parallel_for(getElmId, "get mpCurElmID");
   Kokkos::deep_copy( arrayHost, mpCurElmIDCopy);
+  
 }
 
 void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag){
@@ -239,7 +316,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh,
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries);
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
   kkViewHostU<const double**> mpPositionsIn_h(mpPositionsIn,nComps,numMPs);
 
   if (p_MPs->rebuildOngoing()) {
@@ -249,6 +326,7 @@ void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh,
 
   auto mpPositions = p_MPs->getData<polyMPO::MPF_Cur_Pos_XYZ>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  
   Kokkos::View<double**> mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs);
   Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h);
   auto setPos = PS_LAMBDA(const int&, const int& mp, const int& mask){
@@ -270,7 +348,7 @@ void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh,
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries);
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpPositions = p_MPs->getData<polyMPO::MPF_Cur_Pos_XYZ>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -287,10 +365,72 @@ void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh,
   Kokkos::deep_copy(arrayHost, mpPositionsCopy);
 }
 
+
+void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh,
+                            const int nComps,
+                            const int numMPs,
+                            const double* mpPositionsIn){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  kkViewHostU<const double**> mpPositionsIn_h(mpPositionsIn,nComps,numMPs);
+
+  if (p_MPs->rebuildOngoing()) {
+    p_MPs->setRebuildMPSlice<polyMPO::MPF_Cur_Pos_XYZ>(mpPositionsIn_h);
+    return;
+  }
+
+  auto  mpPositions = p_MPs->getData<polyMPO::MPF_Tgt_Pos_XYZ>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  Kokkos::View<double**> mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs);
+  Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h);
+  auto setPos = PS_LAMBDA(const int&, const int& mp, const int& mask){
+    if(mask){
+      mpPositions(mp,0) = mpPositionsIn_d(0, mpAppID(mp));
+      mpPositions(mp,1) = mpPositionsIn_d(1, mpAppID(mp));
+      mpPositions(mp,2) = mpPositionsIn_d(2, mpAppID(mp));
+    }
+  };
+  p_MPs->parallel_for(setPos, "setMPPositions");
+  pumipic::RecordTime("PolyMPO_setMPTgtPositions", timer.seconds());
+}
+
+void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh,
+                            const int nComps,
+                            const int numMPs,
+                            double* mpPositionsHost){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(nComps == vec3d_nEntries);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+
+  auto mpPositions = p_MPs->getData<polyMPO::MPF_Tgt_Pos_XYZ>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  Kokkos::View<double**> mpPositionsCopy("mpPositionsCopy",vec3d_nEntries,numMPs);
+  auto getPos = PS_LAMBDA(const int&, const int& mp, const int& mask){
+    if(mask){
+      mpPositionsCopy(0,mpAppID(mp)) = mpPositions(mp,0);
+      mpPositionsCopy(1,mpAppID(mp)) = mpPositions(mp,1);
+      mpPositionsCopy(2,mpAppID(mp)) = mpPositions(mp,2);
+    }
+  };
+  p_MPs->parallel_for(getPos, "getMPPositions");
+  kkDbl2dViewHostU arrayHost(mpPositionsHost,nComps,numMPs);
+  Kokkos::deep_copy(arrayHost, mpPositionsCopy);
+  pumipic::RecordTime("PolyMPO_getMPTgtPositions", timer.seconds());
+}
+
+
 void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh,
                          const int nComps,
                          const int numMPs,
                          const double* mpRotLatLonIn){
+  Kokkos::Timer timer;
   static int callCount = 0;
   PMT_ALWAYS_ASSERT(callCount == 0);
   checkMPMeshValid(p_mpmesh);
@@ -312,6 +452,7 @@ void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh,
   };
   p_MPs->parallel_for(setPos, "setMPRotLatLon");
   callCount++;
+  pumipic::RecordTime("PolyMPO_setMPRotLatLon", timer.seconds());
 }
 
 void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh,
@@ -338,13 +479,67 @@ void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh,
   Kokkos::deep_copy(arrayHost, mpRotLatLonCopy);
 }
 
+
+void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh,
+                         const int nComps,
+                         const int numMPs,
+                         const double* mpRotLatLonIn){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+
+  auto mpRotLatLon = p_MPs->getData<polyMPO::MPF_Tgt_Pos_Rot_Lat_Lon>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  kkViewHostU<const double**> mpRotLatLonIn_h(mpRotLatLonIn,nComps,numMPs);
+  Kokkos::View<double**> mpRotLatLonIn_d("mpRotLatLonDevice",vec2d_nEntries,numMPs);
+  Kokkos::deep_copy(mpRotLatLonIn_d, mpRotLatLonIn_h);
+  auto setPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+    if(mask){
+      mpRotLatLon(mp,0) = mpRotLatLonIn_d(0, mpAppID(mp));
+      mpRotLatLon(mp,1) = mpRotLatLonIn_d(1, mpAppID(mp));
+    }
+  };
+  p_MPs->parallel_for(setPos, "setMPTgtRotLatLon");
+  pumipic::RecordTime("PolyMPO_setMPTgtRotLatLon", timer.seconds());
+}
+
+void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh,
+                         const int nComps,
+                         const int numMPs,
+                         double* mpRotLatLonHost){
+  Kokkos::Timer timer;
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+
+  auto mpRotLatLon = p_MPs->getData<polyMPO::MPF_Tgt_Pos_Rot_Lat_Lon>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  Kokkos::View<double**> mpRotLatLonCopy("mpRotLatLonCopy",vec2d_nEntries,numMPs);
+  auto getPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+    if(mask){
+      mpRotLatLonCopy(0,mpAppID(mp)) = mpRotLatLon(mp,0);
+      mpRotLatLonCopy(1,mpAppID(mp)) = mpRotLatLon(mp,1);
+    }
+  };
+  p_MPs->parallel_for(getPos, "getMPRotLatLon");
+  kkDbl2dViewHostU arrayHost(mpRotLatLonHost,nComps,numMPs);
+  Kokkos::deep_copy(arrayHost, mpRotLatLonCopy);
+  pumipic::RecordTime("PolyMPO_getMPTgtRotLatLon", timer.seconds());
+}
+
+
 void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn) {
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == 1); //TODO mp_sclr_t
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpMass = p_MPs->getData<polyMPO::MPF_Mass>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -366,7 +561,7 @@ void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == 1); //TODO mp_sclr_t
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpMass = p_MPs->getData<polyMPO::MPF_Mass>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -388,7 +583,7 @@ void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs,
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpVel = p_MPs->getData<polyMPO::MPF_Vel>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -411,7 +606,7 @@ void polympo_getMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs,
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
   PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
-  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+  //PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpVel = p_MPs->getData<polyMPO::MPF_Vel>();
   auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
@@ -743,6 +938,38 @@ void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* x
   }
 }
 
+void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle){
+
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+  //copy the host array to the device
+  auto dualArea = p_mesh->getMeshField<polyMPO::MeshF_DualTriangleArea>();
+  auto h_dualArea = Kokkos::create_mirror_view(dualArea);
+  for(int i=0; i<nVertices; i++)
+    h_dualArea(i,0) = areaTriangle[i];
+  Kokkos::deep_copy(dualArea, h_dualArea);
+
+}
+
+void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle){
+  
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices()==nVertices);
+  //copy the device to host 
+  auto dualArea = p_mesh->getMeshField<polyMPO::MeshF_DualTriangleArea>();
+  auto h_dualArea = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), dualArea);
+  for(int i=0; i<nVertices; i++)
+    areaTriangle[i] = h_dualArea(i,0);
+
+}
+
+
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
@@ -851,11 +1078,16 @@ void polympo_getMeshElmMass_f(MPMesh_ptr p_mpmesh, const int nCells, double* elm
   pumipic::RecordTime("PolyMPO_getMeshElmMass", timer.seconds());
 }
 
+//Increments in vertex velcoity and displacement
 void polympo_setMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
+  
+  Kokkos::Timer timer;
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
   auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
-  kkViewHostU<const double*[vec2d_nEntries]> arrayHost(array,nVertices);
+  kkViewHostU<const double**> arrayHost(array,nComps,nVertices);
+  Kokkos::View<double**> array_d("meshVelIncrDevice",nComps,nVertices);
+  Kokkos::deep_copy(array_d, arrayHost);
 
   auto vtxField = p_mesh->getMeshField<polyMPO::MeshF_OnSurfVeloIncr>();
 
@@ -864,24 +1096,11 @@ void polympo_setMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
   PMT_ALWAYS_ASSERT(static_cast<size_t>(nVertices*vec2d_nEntries)==vtxField.size());
 
   //copy the host array to the device
-  Kokkos::deep_copy(vtxField,arrayHost);
-}
-
-void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) {
-  //check mpMesh is valid
-  checkMPMeshValid(p_mpmesh);
-  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
-  kkVec2dViewHostU arrayHost(array,nVertices);
-
-  auto vtxField = p_mesh->getMeshField<polyMPO::MeshF_OnSurfVeloIncr>();
-
-  //check the size
-  PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
-  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); 
-  PMT_ALWAYS_ASSERT(static_cast<size_t>(nVertices*vec2d_nEntries)==vtxField.size());
-
-  //copy the device array to the host
-  Kokkos::deep_copy(arrayHost, vtxField);
+  Kokkos::parallel_for("set mesh dispIncr", nVertices, KOKKOS_LAMBDA(const int iVtx){
+    vtxField(iVtx,0) = array_d(0,iVtx);
+    vtxField(iVtx,1) = array_d(1,iVtx);
+  });
+  pumipic::RecordTime("PolyMPO_setMeshVtxOnSurfVelIncr", timer.seconds());
 }
 
 void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array) {
@@ -907,6 +1126,29 @@ void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
   pumipic::RecordTime("PolyMPO_setMeshVtxOnSurfDispIncr", timer.seconds());
 }
 
+
+void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) {
+  //check mpMesh is valid
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  kkDbl2dViewHostU arrayHost(array,nComps,nVertices);
+  Kokkos::View<double**> array_d("meshVelIncrDevice",nComps,nVertices);
+
+  auto vtxField = p_mesh->getMeshField<polyMPO::MeshF_OnSurfVeloIncr>();
+
+  //check the size
+  PMT_ALWAYS_ASSERT(nComps == vec2d_nEntries);
+  PMT_ALWAYS_ASSERT(p_mesh->getNumVertices() == nVertices); 
+  PMT_ALWAYS_ASSERT(static_cast<size_t>(nVertices*vec2d_nEntries)==vtxField.size());
+
+  //copy the device array to the host
+  Kokkos::parallel_for("get mesh dispIncr", nVertices, KOKKOS_LAMBDA(const int iVtx){
+    array_d(0,iVtx) = vtxField(iVtx,0);
+    array_d(1,iVtx) = vtxField(iVtx,1);
+  });
+  Kokkos::deep_copy(arrayHost, array_d);
+}
+
 void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array) {
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);
@@ -929,9 +1171,30 @@ void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
   Kokkos::deep_copy(arrayHost, array_d);
 }
 
+bool polympo_push1P_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  bool is_migrating=((polyMPO::MPMesh*)p_mpmesh)->push1P();
+  return is_migrating;
+}
+
+void polympo_push_ahead_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  ((polyMPO::MPMesh*)p_mpmesh)->push_ahead();
+}
+
+void polympo_push_swap_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  ((polyMPO::MPMesh*)p_mpmesh)->push_swap();
+}
+
+void polympo_push_swap_pos_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  ((polyMPO::MPMesh*)p_mpmesh)->push_swap_pos();
+}
+
 void polympo_push_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
-  ((polyMPO::MPMesh*)p_mpmesh) ->push();
+  ((polyMPO::MPMesh*)p_mpmesh)->push();
 }
 
 //TODO skeleton of reconstruction functions
@@ -973,6 +1236,50 @@ void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, c
   (void)meshEntType;
 }
 
+//With MPI communication done via MPAS
+void polympo_vtxSubAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array){
+  checkMPMeshValid(p_mpmesh);
+  auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm);
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+  PMT_ALWAYS_ASSERT(comp == 0 || comp== 1);  //either first or second component
+  mpmesh->subAssemblyVtx1<polyMPO::MeshF_VtxMass>(vtxPerElm, nCells, comp, array);
+}
+
+void polympo_vtxSubAssemblyVelocity_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array){
+  checkMPMeshValid(p_mpmesh);
+  auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm);
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+  PMT_ALWAYS_ASSERT(comp == 0 || comp== 1);  //either first or second component
+  mpmesh->subAssemblyVtx1<polyMPO::MeshF_Vel>(vtxPerElm, nCells, comp, array);
+}
+
+void polympo_subAssemblyCoeffs_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, double* m11, double* m12, double* m13, double* m14,
+                                                      double* m22, double* m23, double* m24,
+                                                      double* m33, double* m34,
+                                                      double* m44){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  PMT_ALWAYS_ASSERT(vtxPerElm <= maxVtxsPerElm);
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+  auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  mpmesh->subAssemblyCoeffs(vtxPerElm, nCells, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44); 
+}
+
+void polympo_regularize_and_solve_matrix_f(MPMesh_ptr p_mpmesh, int nVertices, double* m11, double* m12, double* m13, double* m14,
+                                                                double* m22, double* m23, double* m24,
+                                                                double* m33, double* m34,
+                                                                double* m44){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+  PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices());
+  auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  mpmesh->solveMatrixAndRegularize(nVertices, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44);
+}
+
 void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
@@ -993,6 +1300,40 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a
   p_mesh->setOwningProc(owningProc);
 }
 
+void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 
+  Kokkos::View<int*, Kokkos::HostSpace> arrayHost("arrayHost", nCells);
+  for (int i = 0; i < nCells; i++) {
+    arrayHost(i) = array[i] - 1;  // TODO right now elmID offset is set after MPs initialized
+  }
+  //check the size
+  PMT_ALWAYS_ASSERT(nCells == p_mesh->getNumElements());
+
+  Kokkos::View<int*> elmGlobal("elmGlobal",nCells);
+  Kokkos::deep_copy(elmGlobal, arrayHost);
+  p_mesh->setElmGlobal(elmGlobal);
+}
+
+void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array){
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh; 
+  PMT_ALWAYS_ASSERT(nVertices==p_mesh->getNumVertices());
+  Kokkos::View<int*, Kokkos::HostSpace> arrayHost("arrayHost", nVertices);
+  for (int i = 0; i < nVertices; i++) {
+    arrayHost(i) = array[i] - 1;  // TODO right now elmID offset is set after MPs initialized
+  }
+
+  Kokkos::View<int*> vtxGlobal("vtxGlobal",nVertices);
+  Kokkos::deep_copy(vtxGlobal, arrayHost);
+  p_mesh->setVtxGlobal(vtxGlobal);
+}
+
+int polympo_getMPCount_f(MPMesh_ptr p_mpmesh) {
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  return p_MPs->getCount();
+}
+
 void polympo_enableTiming_f(){
   pumipic::EnableTiming();
 }
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 6c38106..d39c195 100644
--- a/src/pmpo_c.h
+++ b/src/pmpo_c.h
@@ -21,16 +21,30 @@ void polympo_setMPICommunicator_f(MPMesh_ptr p_mpmesh, MPI_Fint fcomm);
 //MP info
 void polympo_createMPs_f(MPMesh_ptr p_mpmesh, const int numElms, const int numMPs, int* mpsPerElm, const int* mp2Elm, const int* isMPActive);
 void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh, const int numMPs, const int* allTgtMpElmIn, const int* addedMPMask);
+void polympo_startRebuildMPs2_f(MPMesh_ptr p_mpmesh, const int size1, const int* arg1, const int size2, const int size3, int* arg2, int* arg3);
+
+int polympo_getMPCount_f(MPMesh_ptr p_mpmesh);
+
 void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh);
+
 void polympo_setAppIDFunc_f(MPMesh_ptr p_mpmesh, IntVoidFunc getNext, void* appIDs);
+
+void polympo_getMPTgtElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs);
 void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh, const int numMPs, int* elmIDs);
 void polympo_setMPLatLonRotatedFlag_f(MPMesh_ptr p_mpmesh, const int isRotateFlag);
 
 //MP slices
+//Positions
 void polympo_setMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn);
 void polympo_getMPPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsIn);
+void polympo_setMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpPositionsIn);
+void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpPositionsIn);
+//LatLon
 void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn);
 void polympo_getMPRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost);
+void polympo_setMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpRotLatLonIn);
+void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpRotLatLonHost);
+
 void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn);
 void polympo_getMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, double* mpMassHost);
 void polympo_setMPVel_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpVelIn);
@@ -57,6 +71,9 @@ void polympo_setMeshNumEdgesPerElm_f(MPMesh_ptr p_mpmesh, const int nCells, cons
 void polympo_setMeshElm2VtxConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
 void polympo_setMeshElm2ElmConn_f(MPMesh_ptr p_mpmesh, const int maxEdges, const int nCells, const int* array);
 void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
+void polympo_setElmGlobal_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array);
+void polympo_setVtxGlobal_f(MPMesh_ptr p_mpmesh, const int nVertices, const int* array);
+
 
 //Mesh fields
 int polympo_getMeshFVtxType_f();
@@ -77,9 +94,16 @@ void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
 void polympo_getMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray);
 void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* xArray, double* yArray, double* zArray);
-
+//Area Triangle
+void polympo_setMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* areaTriangle);
+void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices, double* areaTriangle);
+  
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);
+void polympo_push_ahead_f(MPMesh_ptr p_mpmesh);
+bool polympo_push1P_f(MPMesh_ptr p_mpmesh);
+void polympo_push_swap_f(MPMesh_ptr p_mpmesh);
+void polympo_push_swap_pos_f(MPMesh_ptr p_mpmesh);
 
 // Reconstruction of variables from MPs to mesh vertices
 void polympo_setReconstructionOfMass_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType);
@@ -88,6 +112,18 @@ void polympo_setReconstructionOfStrainRate_f(MPMesh_ptr p_mpmesh, const int orde
 void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, const int meshEntType);
 void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh);
 
+//Reconstruction using MPAS
+void polympo_vtxSubAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array);
+void polympo_vtxSubAssemblyVelocity_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int comp, double* array);
+void polympo_subAssemblyCoeffs_f(MPMesh_ptr p_mpmesh, 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 polympo_regularize_and_solve_matrix_f(MPMesh_ptr p_mpmesh, int nVetices, double* m11, double* m12, double* m13, double* m14,
+                                                                double* m22, double* m23, double* m24,
+                                                                double* m33, double* m34,
+                                                                double* m44);
+ 
 // Timing
 void polympo_enableTiming_f();
 void polympo_summarizeTime_f();
diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h
index 34857b9..4539a63 100644
--- a/src/pmpo_defines.h
+++ b/src/pmpo_defines.h
@@ -4,10 +4,11 @@
 
 #define MP_ACTIVE 1
 #define MP_DELETE -1
-
+#define INVALID_ELM_ID -1
 typedef void* MPMesh_ptr;
 //Function that receives void* and returns an int
 typedef int (*IntVoidFunc)(void*);
+typedef void (*VoidVoidFunc)(void*, const int, const int, const int, int&);
 
 using space_t = Kokkos::DefaultExecutionSpace::memory_space;
 
diff --git a/src/pmpo_fortran.f90 b/src/pmpo_fortran.f90
index f16f57a..10d04f2 100644
--- a/src/pmpo_fortran.f90
+++ b/src/pmpo_fortran.f90
@@ -87,6 +87,20 @@ subroutine polympo_startRebuildMPs(mpMesh, numMPs, allMP2Elm, addedMPMask) &
     type(c_ptr), intent(in), value :: allMP2Elm
     type(c_ptr), intent(in), value :: addedMPMask
   end subroutine
+
+
+  subroutine polympo_startRebuildMPs2(mpMesh, size1, arg1, size2, size3, arg2, arg3) &
+    bind(C, NAME='polympo_startRebuildMPs2_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: size1
+    type(c_ptr), intent(in), value :: arg1
+    integer(c_int), value :: size2
+    integer(c_int), value :: size3
+    type(c_ptr), value :: arg2
+    type(c_ptr), value :: arg3
+  end subroutine
+  !
   !---------------------------------------------------------------------------
   !> @brief called after startRebuild()
   !> @brief called after initializing MP fields
@@ -110,6 +124,7 @@ subroutine polympo_setAppIDFunc(mpMesh, getNext, appIDs) &
     type(c_funptr), value :: getNext
     type(c_ptr), value :: appIDs
   end subroutine
+  
   !---------------------------------------------------------------------------
   !> @brief get the current element ID MP array from a polympo array
   !> @param mpmesh(in/out) MPMesh object
@@ -123,6 +138,16 @@ subroutine polympo_getMPCurElmID(mpMesh, numMPs, array) &
     integer(c_int), value :: numMPs
     type(c_ptr), value :: array
   end subroutine
+
+
+  subroutine polympo_getMPTgtElmID(mpMesh, numMPs, array) &
+             bind(C, NAME='polympo_getMPTgtElmID_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: numMPs
+    type(c_ptr), value :: array
+  end subroutine
+  !
   !---------------------------------------------------------------------------
   !> @brief set the mp lat lon is rotational or normal
   !> @param mpmesh(in/out) MPMesh object
@@ -134,6 +159,8 @@ subroutine polympo_setMPLatLonRotatedFlag(mpMesh, isRotateFlag) &
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: isRotateFlag
   end subroutine
+  
+  
   !---------------------------------------------------------------------------
   !> @brief set the MP positions array from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -163,6 +190,37 @@ subroutine polympo_getMPPositions(mpMesh, nComps, numMPs, array) &
     integer(c_int), value :: nComps, numMPs
     type(c_ptr), value :: array
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set the MP positions array from a host array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nComps(in) number of components, should always be 3
+  !> @param numMPs(in) number of the MPs
+  !> @param array(in) MP current position 2D array (3,numMPs), allocated by user on host
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMPTgtPositions(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_setMPTgtPositions_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the MP positions array from a polympo array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nComps(in) number of components, should always be 3
+  !> @param numMPs(in) number of the MPs
+  !> @param array(in/out) output MP current position 2D array (3,numMPs),
+  !>                      allocated by user
+  !---------------------------------------------------------------------------
+  subroutine polympo_getMPTgtPositions(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_getMPTgtPositions_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+
+
   !---------------------------------------------------------------------------
   !> @brief set the MP latitude and longtitude array from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -193,6 +251,40 @@ subroutine polympo_getMPRotLatLon(mpMesh, nComps, numMPs, array) &
     integer(c_int), value :: nComps, numMPs
     type(c_ptr), value :: array
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set the MP latitude and longtitude array from a host array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nComps(in) number of components, should always be 2
+  !> @param numMPs(in) number of the MPs
+  !> @param array(in)  input MP current lat and lon 2D array (2,numMPs),
+  !>                   allocated by user
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMPTgtRotLatLon(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_setMPTgtRotLatLon_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), intent(in), value :: array
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the MP latitude and longtitude array from a polympo array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nComps(in) number of components, should always be 2
+  !> @param numMPs(in) number of the MPs
+  !> @param array(in/out) output MP current lat and lon 2D array (2,numMPs),
+  !>                      allocated by user
+  !---------------------------------------------------------------------------
+  subroutine polympo_getMPTgtRotLatLon(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_getMPTgtRotLatLon_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+
+  
+  
+  
   !---------------------------------------------------------------------------
   !> @brief set the Mass MP array from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -435,6 +527,7 @@ integer function polympo_getMeshFElmType() &
              bind(C, NAME='polympo_getMeshFElmType_f')
     use :: iso_c_binding
   end function
+
   !---------------------------------------------------------------------------
   !> @brief set the polympo mesh vertices coordinates
   !> @param mpmesh(in/out) MPMesh object
@@ -513,6 +606,34 @@ subroutine polympo_getMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) &
     integer(c_int), value :: nCells
     type(c_ptr), value :: xArray, yArray, zArray
   end subroutine
+
+  !---------------------------------------------------------------------------
+  !> @brief set the polympo dual mesh triangle area
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVertices(in) length of array in, use for assertion 
+  !> @param Array(in) 1D array of area of dual triangle elements
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
+             bind(C, NAME='polympo_setMeshDualTriangleArea_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: areaTriangle
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the polympo mesh dual mesh triangle area
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nVetices(in) length of array in, use for assertion
+  !> @param Array(in/out) 1D array of area of dual triangle elements
+  !---------------------------------------------------------------------------
+  subroutine polympo_getMeshDualTriangleArea(mpMesh, nVertices, areaTriangle) &
+             bind(C, NAME='polympo_getMeshDualTriangleArea_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: areaTriangle
+  end subroutine
+
   !---------------------------------------------------------------------------
   !> @brief set the vertices velocity from a host array
   !> @param mpmesh(in/out) MPMesh object
@@ -670,6 +791,29 @@ subroutine polympo_setOwningProc(mpMesh, nCells, array) &
     integer(c_int), value :: nCells
     type(c_ptr), intent(in), value :: array
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set the owning process array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nCells(in) number of cells
+  !> @param array(in) input mesh cell to process array
+  !---------------------------------------------------------------------------
+  subroutine polympo_setElmGlobal(mpMesh, nCells, array) &
+    bind(C, NAME='polympo_setElmGlobal_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nCells
+    type(c_ptr), intent(in), value :: array
+  end subroutine
+  subroutine polympo_setVtxGlobal(mpMesh, nVertices, array) &
+    bind(C, NAME='polympo_setVtxGlobal_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), intent(in), value :: array
+  end subroutine
+  
+
+  
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices
@@ -681,6 +825,37 @@ subroutine polympo_push(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
+
+
+  !---------------------------------------------------------------------------
+  !> @brief calculate the MPs from given mesh vertices rotational latitude
+  logical function polympo_push1P(mpMesh) &
+             bind(C, NAME='polympo_push1P_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end function
+
+ subroutine polympo_push_ahead(mpMesh) &
+             bind(C, NAME='polympo_push_ahead_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+ 
+  subroutine polympo_push_swap(mpMesh) &
+             bind(C, NAME='polympo_push_swap_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+  
+  subroutine polympo_push_swap_pos(mpMesh) &
+             bind(C, NAME='polympo_push_swap_pos_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
+
+
+
+
   !---------------------------------------------------------------------------
   !> @brief start the reconstruction of MP Mass to Mesh Vertices
   !> @param mpmesh(in/out) MPMesh object
@@ -733,6 +908,40 @@ subroutine polympo_setReconstructionOfStress(mpMesh, order, meshEntType) &
     type(c_ptr), value :: mpMesh
     integer(c_int), value :: order, meshEntType
   end subroutine
+
+
+  subroutine polympo_vtxSubAssemblyIceArea(mpMesh, vtxPerElm, nCells, comp, array) &
+             bind(C, NAME='polympo_vtxSubAssemblyIceArea_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: vtxPerElm, nCells, comp
+    type(c_ptr), value :: array
+  end subroutine
+
+  subroutine polympo_vtxSubAssemblyVelocity(mpMesh, vtxPerElm, nCells, comp, array) &
+             bind(C, NAME='polympo_vtxSubAssemblyVelocity_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: vtxPerElm, nCells, comp
+    type(c_ptr), value :: array
+  end subroutine
+  
+  subroutine polympo_subAssemblyCoeffs(mpMesh, vtxPerElm, nCells, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44) &
+             bind(C, NAME='polympo_subAssemblyCoeffs_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: vtxPerElm, nCells
+    type(c_ptr), value :: m11, m12, m13, m14, m22, m23, m24, m33, m34, m44
+  end subroutine
+  
+  subroutine polympo_regularize_and_solve_matrix(mpMesh, nVertices, m11, m12, m13, m14, m22, m23, m24, m33, m34, m44) &
+             bind(C, NAME='polympo_regularize_and_solve_matrix_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nVertices
+    type(c_ptr), value :: m11, m12, m13, m14, m22, m23, m24, m33, m34, m44
+  end subroutine
+ 
   !---------------------------------------------------------------------------
   !> @brief directly call the reconstruct of the MP fields to mesh fields
   !> @param mpmesh(in/out) MPMesh object
@@ -761,6 +970,13 @@ subroutine polympo_setTimingVerbosity(v) bind(C, NAME='polympo_setTimingVerbosit
     use :: iso_c_binding
     integer(c_int), value :: v
   end subroutine
+
+  integer function polympo_getMPCount(mpMesh) bind(C, name="polympo_getMPCount_f")
+    use :: iso_c_binding
+    implicit none
+    type(c_ptr), value :: mpMesh
+  end function
+
   end interface
   contains
   !---------------------------------------------------------------------------
@@ -776,4 +992,5 @@ subroutine polympo_checkPrecisionForRealKind(APP_RKIND)
         call exit(1)
     end if
   end subroutine
+
 end module
diff --git a/src/pmpo_materialPoints.cpp b/src/pmpo_materialPoints.cpp
index c20fe00..0f6b6c1 100644
--- a/src/pmpo_materialPoints.cpp
+++ b/src/pmpo_materialPoints.cpp
@@ -8,11 +8,13 @@ template<typename MemSpace = defaultSpace, typename View>
 pumipic::MemberTypeViews createInternalMemberViews(int numMPs, View mp2elm, View mpAppID){
   auto mpInfo = ps::createMemberViews<MaterialPointTypes, MemSpace>(numMPs);
   auto mpCurElmPos_m = ps::getMemberView<MaterialPointTypes, MPF_Cur_Elm_ID, MemSpace>(mpInfo);
+  auto mpTgtElmPos_m = ps::getMemberView<MaterialPointTypes, MPF_Tgt_Elm_ID, MemSpace>(mpInfo);
   auto mpAppID_m = ps::getMemberView<MaterialPointTypes, MPF_MP_APP_ID, MemSpace>(mpInfo);
   auto mpStatus_m = ps::getMemberView<MaterialPointTypes, MPF_Status, MemSpace>(mpInfo);
   auto policy = Kokkos::RangePolicy<typename MemSpace::execution_space>(typename MemSpace::execution_space(), 0, numMPs);
   Kokkos::parallel_for("setMPinfo", policy, KOKKOS_LAMBDA(int i) {
     mpCurElmPos_m(i) = mp2elm(i);
+    mpTgtElmPos_m(i) = INVALID_ELM_ID;
     mpStatus_m(i) = MP_ACTIVE;
     mpAppID_m(i) = mpAppID(i);
   });
@@ -40,8 +42,12 @@ PS* createDPS(int numElms, int numMPs, MPSView<MPF_Cur_Pos_XYZ> positions, IntVi
   return dps;
 }
 
-PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID) {
+PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global) {
   PS::kkGidView elmGids("elementGlobalIds", numElms); //TODO - initialize this to [0..numElms)
+  Kokkos::parallel_for("setGids", numElms, KOKKOS_LAMBDA(const int elm){
+    elmGids(elm) = elm2global(elm);
+  });
+  
   auto mpInfo = createInternalMemberViews(numMPs, mp2elm, mpAppID);
   Kokkos::TeamPolicy<Kokkos::DefaultExecutionSpace> policy(numElms,Kokkos::AUTO);
   auto dps = new DPS<MaterialPointTypes>(policy, numElms, numMPs, mpsPerElm, elmGids, mp2elm, mpInfo);
@@ -57,8 +63,8 @@ MaterialPoints::MaterialPoints(int numElms, int numMPs, MPSView<MPF_Cur_Pos_XYZ>
   operating_mode = MP_RELEASE;
 };
 
-MaterialPoints::MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID) {
-  MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, mpAppID);
+MaterialPoints::MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global){
+  MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, mpAppID, elm2global);
   updateMaxAppID();
   operating_mode = MP_RELEASE;
 };
@@ -91,6 +97,16 @@ void MaterialPoints::startRebuild(IntView tgtElm, int addedNumMPs, IntView added
   rebuildFields.addedSlices_h = createInternalMemberViews<hostSpace>(addedNumMPs, addedMP2elm_h, addedMPAppID_h);
 }
 
+void MaterialPoints::startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID) {
+  rebuildFields.ongoing = true;
+  rebuildFields.addedNumMPs = addedNumMPs;
+  rebuildFields.addedMP2elm = addedMP2elm;
+  rebuildFields.allTgtElm = tgtElm;
+  auto addedMP2elm_h = Kokkos::create_mirror_view_and_copy(hostSpace(), addedMP2elm);
+  auto addedMPAppID_h = Kokkos::create_mirror_view_and_copy(hostSpace(), addedMPAppID);
+  rebuildFields.addedSlices_h = createInternalMemberViews<hostSpace>(addedNumMPs, addedMP2elm_h, addedMPAppID_h);
+}
+
 void MaterialPoints::finishRebuild() {
   auto addedSlices_d = ps::createMemberViews<MaterialPointTypes, defaultSpace>(rebuildFields.addedNumMPs);
   ps::CopyMemSpaceToMemSpace<defaultSpace, hostSpace, MaterialPointTypes>(addedSlices_d, rebuildFields.addedSlices_h);
@@ -109,14 +125,34 @@ void MaterialPoints::setMPIComm(MPI_Comm comm) {
   mpi_comm = comm;
 }
 
-bool MaterialPoints::migrate() {
+bool MaterialPoints::check_migrate(){
+  
+  Kokkos::Timer timer; 
+  auto MPs2Proc = getData<MPF_Tgt_Proc_ID>();
+  IntView isMigrating("isMigrating", 1);
+  int rank;
+  MPI_Comm_rank(mpi_comm, &rank);
+
+  auto setMigrationFields = PS_LAMBDA(const int& e, const int& mp, const bool& mask) {
+    if (mask) {
+      if (MPs2Proc(mp) != rank) isMigrating(0) = 1;
+    }
+  };
+  parallel_for(setMigrationFields, "setMigrationFields");
+  
+  if (getOpMode() == polyMPO::MP_DEBUG)
+    printf("Material point check migration: %f\n", timer.seconds());
+  pumipic::RecordTime("PolyMPO_check_migrate", timer.seconds());
+  return pumipic::getLastValue(isMigrating) > 0;
+}
+
+void MaterialPoints::migrate() {
   Kokkos::Timer timer;
   auto MPs2Elm = getData<MPF_Tgt_Elm_ID>();
   auto MPs2Proc = getData<MPF_Tgt_Proc_ID>();
 
   IntView new_elem("new_elem", MPs->capacity());
   IntView new_process("new_process", MPs->capacity());
-  IntView isMigrating("isMigrating", 1);
 
   int rank;
   MPI_Comm_rank(mpi_comm, &rank);
@@ -124,7 +160,6 @@ bool MaterialPoints::migrate() {
     if (mask) {
       new_elem(mp) = MPs2Elm(mp);
       new_process(mp) = MPs2Proc(mp);
-      if (new_process(mp) != rank) isMigrating(0) = 1;
     }
   };
   parallel_for(setMigrationFields, "setMigrationFields");
@@ -133,7 +168,6 @@ bool MaterialPoints::migrate() {
   if (getOpMode() == polyMPO::MP_DEBUG)
     printf("Material point migration: %f\n", timer.seconds());
   pumipic::RecordTime("PolyMPO_migrate", timer.seconds());
-  return pumipic::getLastValue(isMigrating) > 0;
 }
 
 bool MaterialPoints::rebuildOngoing() { return rebuildFields.ongoing; }
diff --git a/src/pmpo_materialPoints.hpp b/src/pmpo_materialPoints.hpp
index 99d2964..f8391d1 100644
--- a/src/pmpo_materialPoints.hpp
+++ b/src/pmpo_materialPoints.hpp
@@ -18,6 +18,7 @@ using hostSpace = Kokkos::HostSpace;
 using defaultSpace = Kokkos::DefaultExecutionSpace::memory_space;
 
 typedef std::function<int()> IntFunc;
+typedef std::function<int(int)> IntIntFunc;
 
 enum MaterialPointSlice {
   MPF_Status = 0,
@@ -33,6 +34,7 @@ enum MaterialPointSlice {
   MPF_Mass,
   MPF_Vel,
   MPF_Rot_Lat_Lon_Incr,
+  MPF_Vel_Incr,
   MPF_Strain_Rate,
   MPF_Stress,
   MPF_Stress_Div,
@@ -61,6 +63,7 @@ template <> struct mpSliceToMeshField < MPF_Basis_Grad_Vals     > { using type =
 template <> struct mpSliceToMeshField < MPF_Mass                > { using type = doubleSclr_t; };
 template <> struct mpSliceToMeshField < MPF_Vel                 > { using type = vec2d_t; };
 template <> struct mpSliceToMeshField < MPF_Rot_Lat_Lon_Incr    > { using type = vec2d_t; };
+template <> struct mpSliceToMeshField < MPF_Vel_Incr            > { using type = vec2d_t; };
 template <> struct mpSliceToMeshField < MPF_Strain_Rate         > { using type = double[6]; };
 template <> struct mpSliceToMeshField < MPF_Stress              > { using type = double[6]; };
 template <> struct mpSliceToMeshField < MPF_Stress_Div          > { using type = vec3d_t; };
@@ -97,6 +100,7 @@ typedef MemberTypes<mpSliceToMeshField < MPF_Status              >::type,
                     mpSliceToMeshField < MPF_Mass                >::type,
                     mpSliceToMeshField < MPF_Vel                 >::type,
                     mpSliceToMeshField < MPF_Rot_Lat_Lon_Incr    >::type,
+                    mpSliceToMeshField < MPF_Vel_Incr            >::type,
                     mpSliceToMeshField < MPF_Strain_Rate         >::type,
                     mpSliceToMeshField < MPF_Stress              >::type,
                     mpSliceToMeshField < MPF_Stress_Div          >::type,
@@ -130,15 +134,18 @@ class MaterialPoints {
   public:
     MaterialPoints() : MPs(nullptr) {};
     MaterialPoints(int numElms, int numMPs, MPSView<MPF_Cur_Pos_XYZ> positions, IntView mpsPerElm, IntView mp2elm);
-    MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID);
+    MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID, IntView elm2global);
     ~MaterialPoints();
 
     void rebuild(IntView addedMP2elm, IntView addedMPAppID);
     void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID, Kokkos::View<const int*> addedMPMask);
+    void startRebuild(IntView tgtElm, int addedNumMPs, IntView addedMP2elm, IntView addedMPAppID);
+    
     void finishRebuild();
     bool rebuildOngoing();
 
-    bool migrate();
+    bool check_migrate();
+    void migrate();
     MPI_Comm getMPIComm();
     void setMPIComm(MPI_Comm comm);
 
@@ -156,7 +163,7 @@ class MaterialPoints {
     void rebuild() {
       IntView tgtElm("tgtElm", MPs->capacity());
       auto tgtMpElm = MPs->get<MPF_Tgt_Elm_ID>();
-      auto setTgtElm = PS_LAMBDA(const int&, const int& mp, const int& mask) {
+      auto setTgtElm = PS_LAMBDA(const int& e, const int& mp, const int& mask) {
         if(mask) {
           tgtElm(mp) = tgtMpElm(mp);
         }
@@ -214,7 +221,10 @@ class MaterialPoints {
         auto tgtPosRotLatLon = MPs->get<MPF_Tgt_Pos_Rot_Lat_Lon>();
         auto tgtPosXYZ = MPs->get<MPF_Tgt_Pos_XYZ>();
         auto rotLatLonIncr = MPs->get<MPF_Rot_Lat_Lon_Incr>();
-        
+        //Velocity   
+        auto velMPs = MPs->get<MPF_Vel>();
+        auto velIncr = MPs->get<MPF_Vel_Incr>();
+      
         auto is_rotated = getRotatedFlag(); 
         auto updateRotLatLon = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
             if(mask){
@@ -229,10 +239,13 @@ class MaterialPoints {
                   auto xyz_geo = grid_rotation_backward(xyz_rot);
                   lat_lon_from_xyz(geoLat, geoLon, xyz_geo, radius);
                 }	
-                // x = cosLon cosLat, y = sinLon cosLat, z = sinLat
+                // x=cosLon cosLat, y=sinLon cosLat, z= sinLat
                 tgtPosXYZ(mp,0) = radius * std::cos(geoLon) * std::cos(geoLat);
                 tgtPosXYZ(mp,1) = radius * std::sin(geoLon) * std::cos(geoLat);
                 tgtPosXYZ(mp,2) = radius * std::sin(geoLat);
+                
+                velMPs(mp,0) = velMPs(mp,0) + velIncr(mp,0);
+                velMPs(mp,1) = velMPs(mp,1) + velIncr(mp,1);
             } 
         };
         ps::parallel_for(MPs, updateRotLatLon,"updateRotationalLatitudeLongitude"); 
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index cdf9b28..57b838f 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -41,6 +41,10 @@ namespace polyMPO{
         auto vtxRotLatLonIncrMapEntry = meshFields2TypeAndString.at(MeshF_RotLatLonIncr);
         PMT_ALWAYS_ASSERT(vtxRotLatLonIncrMapEntry.first == MeshFType_VtxBased);
         vtxRotLatLonIncr_ = MeshFView<MeshF_RotLatLonIncr>(vtxRotLatLonIncrMapEntry.second,numVtxs_);
+
+        auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea);
+        PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased);
+        dualTriangleArea_ = MeshFView<MeshF_DualTriangleArea>(dualTriangleAreaEntry.second,numVtxs_);
     }
     
     void Mesh::setMeshElmBasedFieldSize(){
@@ -73,5 +77,11 @@ namespace polyMPO{
     IntView Mesh::getElm2Process() {
         return owningProc_;
     }
+    
+    IntView Mesh::getElmGlobal() {
+        return globalElm_;
+    }
+
+
 
 } // namespace polyMPO
diff --git a/src/pmpo_mesh.hpp b/src/pmpo_mesh.hpp
index 3b91dc3..97910fe 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -20,6 +20,7 @@ enum MeshFieldIndex{
     MeshF_VtxCoords,
     MeshF_VtxRotLat,
     MeshF_ElmCenterXYZ,
+    MeshF_DualTriangleArea,
     MeshF_Vel,
     MeshF_VtxMass,
     MeshF_ElmMass,
@@ -38,6 +39,7 @@ template <MeshFieldIndex> struct meshFieldToType;
 template <> struct meshFieldToType < MeshF_VtxCoords         > { using type = Kokkos::View<vec3d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxRotLat         > { using type = DoubleView; };
 template <> struct meshFieldToType < MeshF_ElmCenterXYZ      > { using type = Kokkos::View<vec3d_t*>; };
+template <> struct meshFieldToType < MeshF_DualTriangleArea  > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_Vel               > { using type = Kokkos::View<vec2d_t*>; };
 template <> struct meshFieldToType < MeshF_VtxMass           > { using type = Kokkos::View<doubleSclr_t*>; };
 template <> struct meshFieldToType < MeshF_ElmMass           > { using type = Kokkos::View<doubleSclr_t*>; };
@@ -55,6 +57,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
                {MeshF_VtxCoords,        {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
                {MeshF_VtxRotLat,        {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
                {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
+               {MeshF_DualTriangleArea, {MeshFType_VtxBased,"MeshField_DualTriangleArea"}},
                {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
                {MeshF_VtxMass,          {MeshFType_VtxBased,"MeshField_VerticesMass"}},
                {MeshF_ElmMass,          {MeshFType_ElmBased,"MeshField_ElementsMass"}},
@@ -85,11 +88,13 @@ class Mesh {
     IntVtx2ElmView elm2VtxConn_;
     IntElm2ElmView elm2ElmConn_;
     IntView owningProc_;
-
+    IntView globalElm_;
+    IntView globalVtx_;
     //start of meshFields
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
     MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
+    MeshFView<MeshF_DualTriangleArea> dualTriangleArea_;
     MeshFView<MeshF_Vel> vtxVel_;
     MeshFView<MeshF_VtxMass> vtxMass_;
     MeshFView<MeshF_ElmMass> elmMass_;
@@ -159,6 +164,11 @@ class Mesh {
     void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_);
                                             owningProc_ = owningProc; }
     
+    void setElmGlobal(IntView globalElm) {globalElm_ = globalElm;}
+    IntView getElmGlobal();
+    void setVtxGlobal(IntView globalVtx) {globalVtx_ = globalVtx;}
+    IntView getVtxGlobal() {return globalVtx_;}
+
     void computeRotLatLonIncr();
 };
 
@@ -181,6 +191,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_ElmCenterXYZ){
         return elmCenterXYZ_;
     }
+    else if constexpr (index==MeshF_DualTriangleArea){
+        return dualTriangleArea_;
+    }
     else if constexpr (index==MeshF_Vel){
         return vtxVel_;
     }
diff --git a/src/pmpo_utils.hpp b/src/pmpo_utils.hpp
index b9d5f7d..76c2a7d 100644
--- a/src/pmpo_utils.hpp
+++ b/src/pmpo_utils.hpp
@@ -51,36 +51,36 @@ class Vec2d {
     //constructors
     KOKKOS_INLINE_FUNCTION
     Vec2d():coords_{0.0, 0.0}{}
-    
+ 
     ~Vec2d() = default;
-    
+
     KOKKOS_INLINE_FUNCTION
     Vec2d(double x1, double x2):coords_{x1,x2}{}
-    
+
     KOKKOS_INLINE_FUNCTION
     Vec2d(double x[2]):coords_{x[0], x[1]}{}
 
     //operators 
     KOKKOS_INLINE_FUNCTION
     double operator[](int i) const { return coords_[i]; }
-    
+ 
     KOKKOS_INLINE_FUNCTION
     double &operator[](int i) { return coords_[i]; }
 
     KOKKOS_INLINE_FUNCTION
     Vec2d operator-() { return Vec2d(-coords_[0], -coords_[1]); }
-    
+ 
     KOKKOS_INLINE_FUNCTION
     Vec2d operator+(Vec2d v) { return Vec2d(coords_[0] + v[0],
-                                                 coords_[1] + v[1]); }
+                                            coords_[1] + v[1]); }
     
     KOKKOS_INLINE_FUNCTION
     Vec2d operator-(Vec2d v) { return Vec2d(coords_[0] - v[0],
-                                                 coords_[1] - v[1]); }
+                                            coords_[1] - v[1]); }
     
     KOKKOS_INLINE_FUNCTION
     Vec2d operator*(double scalar) const { return Vec2d(coords_[0]*scalar,
-                                                            coords_[1]*scalar); }
+                                                        coords_[1]*scalar); }
     
     KOKKOS_INLINE_FUNCTION
     double dot(Vec2d v) { return  coords_[0]*v[0]+coords_[1]*v[1]; }
@@ -119,12 +119,12 @@ class Vec3d{
     
     KOKKOS_INLINE_FUNCTION
     Vec3d operator+(const Vec3d& v) const {
-        return Vec3d(coords_[0] + v.coords_[0], coords_[1] + v.coords_[1], coords_[2] + v.coords_[2]);
+      return Vec3d(coords_[0] + v.coords_[0], coords_[1] + v.coords_[1], coords_[2] + v.coords_[2]);
     }
 
     KOKKOS_INLINE_FUNCTION
     Vec3d operator-(const Vec3d& v) const {
-        return Vec3d(coords_[0] - v.coords_[0], coords_[1] - v.coords_[1], coords_[2] - v.coords_[2]);
+      return Vec3d(coords_[0] - v.coords_[0], coords_[1] - v.coords_[1], coords_[2] - v.coords_[2]);
     }
     
     KOKKOS_INLINE_FUNCTION
@@ -132,24 +132,24 @@ class Vec3d{
 
     KOKKOS_INLINE_FUNCTION
     Vec3d operator*(double scalar) const {
-        return Vec3d(coords_[0] * scalar, coords_[1] * scalar, coords_[2] * scalar);
+      return Vec3d(coords_[0] * scalar, coords_[1] * scalar, coords_[2] * scalar);
     }
 
     KOKKOS_INLINE_FUNCTION
     double dot(const Vec3d& v) const {
-        return coords_[0] * v.coords_[0] + coords_[1] * v.coords_[1] + coords_[2] * v.coords_[2];
+      return coords_[0] * v.coords_[0] + coords_[1] * v.coords_[1] + coords_[2] * v.coords_[2];
     }
 
     KOKKOS_INLINE_FUNCTION
     Vec3d cross(const Vec3d& v) const {
-        return Vec3d(coords_[1] * v.coords_[2] - coords_[2] * v.coords_[1],
-                     coords_[2] * v.coords_[0] - coords_[0] * v.coords_[2],
-                     coords_[0] * v.coords_[1] - coords_[1] * v.coords_[0]);
+      return Vec3d(coords_[1] * v.coords_[2] - coords_[2] * v.coords_[1],
+                   coords_[2] * v.coords_[0] - coords_[0] * v.coords_[2],
+                   coords_[0] * v.coords_[1] - coords_[1] * v.coords_[0]);
     }
 
     KOKKOS_INLINE_FUNCTION
     double magnitude() const {
-        return std::sqrt(coords_[0] * coords_[0] + coords_[1] * coords_[1] + coords_[2] * coords_[2]);
+      return std::sqrt(coords_[0] * coords_[0] + coords_[1] * coords_[1] + coords_[2] * coords_[2]);
     }
 };
 
@@ -183,7 +183,7 @@ class Vec4d{
     Vec4d operator+(const Vec4d& v) const {
       return Vec4d(coords_[0] + v.coords_[0], coords_[1] + v.coords_[1], 
              coords_[2] + v.coords_[2], coords_[3] + v.coords_[3]);
-      }
+    }
       
     KOKKOS_INLINE_FUNCTION
     Vec4d operator-(const Vec4d& v) const {
@@ -274,6 +274,18 @@ class Matrix4d {
         data_[i][i]=data_[i][i]+eps;
       }
     }
+
+    KOKKOS_INLINE_FUNCTION
+    void scaleFirstRowAndColumn(double factor) {
+      // Scale the first row
+      for (int j = 0; j < 4; j++) {
+        data_[0][j] *= factor;
+      }
+      // Scale the first column
+      for (int i = 0; i < 4; i++) {
+        data_[i][0] *= factor;
+      }
+    }
 				        
 };
 
@@ -289,148 +301,144 @@ void initArray(Vec2d* arr, int n, Vec2d fill){
 KOKKOS_INLINE_FUNCTION
 void CholeskySolve4d_UnitRHS(Matrix4d& A, double* x){
 
-    double a_00=A(0,0);
-    if (A(0,0)==0){
-        x[0]=0;
-	x[1]=0;
-	x[2]=0;
-	x[3]=0;
-	return;
-    }
-    
-    A(0,0) = std::sqrt(A(0,0));
-    A(0,1) /= A(0,0);
-    A(0,2) /= A(0,0);
-    A(0,3) /= A(0,0);
-
-    double diag = A(1,1) - A(0,1)*A(0,1);
-    if(diag>EPSILON){
-         A(1,1)=std::sqrt(diag);
-	 A(1,2)=(A(1,2)-A(0,1)*A(0,2))/A(1,1);
-	 A(1,3)=(A(1,3)-A(0,1)*A(0,3))/A(1,1);
-    }
-    else{
-        A(1,1)=0.0;
-	A(1,2)=0.0;
-	A(1,3)=0.0;
-    }
-    
-    diag = A(2,2) - A(0,2)*A(0,2)-A(1,2)*A(1,2);
-    if(diag>EPSILON){
-        A(2,2)=std::sqrt(diag);
-        A(2,3)=(A(2,3)-A(0,2)*A(0,3)-A(1,2)*A(1,3))/A(2,2);
-    }
-    else{
-       A(2,2)=0.0;
-       A(2,3)=0.0;
-    }
+  double a_00=A(0,0);
+  if (A(0,0)==0){
+    x[0]=0;
+    x[1]=0;
+    x[2]=0;
+    x[3]=0;
+    return;
+  }
 
-    diag=A(3,3)-A(0,3)*A(0,3)-A(1,3)*A(1,3)-A(2,3)*A(2,3);
-    if(diag>EPSILON)
-       A(3,3)=std::sqrt(diag);
-    else 
-       A(3,3)=0.0;
-
-    if(abs(A(1,1))<EPSILON || abs(A(2,2))<EPSILON || abs(A(3,3))<EPSILON){
-       x[0]=1.0/a_00;
-       x[1]=0.0;
-       x[2]=0.0;
-       x[3]=0.0;
-    }
-    else{
-       x[0]= 1.0/A(0,0);
-       x[1]= -A(0,1)*x[0]/A(1,1);                
-       x[2]= -(A(0,2)*x[0]+A(1,2)*x[1])/A(2,2);  
-       x[3]= -(A(0,3)*x[0]+A(1,3)*x[1]+A(2,3)*x[2])/A(3,3); 
-
-       x[3] = x[3]/A(3,3);
-       x[2] = ( x[2] - A(2,3)*x[3] )/A(2,2);
-       x[1] = ( x[1] - A(1,2)*x[2] - A(1,3)*x[3])/A(1,1);
-       x[0] = ( x[0] - A(0,1)*x[1] - A(0,2)*x[2] - A(0,3)*x[3])/A(0,0);
-    } 
+  A(0,0) = std::sqrt(A(0,0));
+  A(0,1) /= A(0,0);
+  A(0,2) /= A(0,0);
+  A(0,3) /= A(0,0);
+
+  double diag = A(1,1) - A(0,1)*A(0,1);
+  if(diag>EPSILON){
+    A(1,1)=std::sqrt(diag);
+    A(1,2)=(A(1,2)-A(0,1)*A(0,2))/A(1,1);
+    A(1,3)=(A(1,3)-A(0,1)*A(0,3))/A(1,1);
+  }
+  else{
+    A(1,1)=0.0;
+    A(1,2)=0.0;
+    A(1,3)=0.0;
+  }
+
+  diag = A(2,2) - A(0,2)*A(0,2)-A(1,2)*A(1,2);
+  if(diag>EPSILON){
+    A(2,2)=std::sqrt(diag);
+    A(2,3)=(A(2,3)-A(0,2)*A(0,3)-A(1,2)*A(1,3))/A(2,2);
+  }
+  else{
+    A(2,2)=0.0;
+    A(2,3)=0.0;
+  }
+
+  diag=A(3,3)-A(0,3)*A(0,3)-A(1,3)*A(1,3)-A(2,3)*A(2,3);
+  if(diag>EPSILON)
+    A(3,3)=std::sqrt(diag);
+  else 
+    A(3,3)=0.0;
+
+  if(abs(A(1,1))<EPSILON || abs(A(2,2))<EPSILON || abs(A(3,3))<EPSILON){
+    x[0]=1.0/a_00;
+    x[1]=0.0;
+    x[2]=0.0;
+    x[3]=0.0;
+  }
+  else{
+    x[0]= 1.0/A(0,0);
+    x[1]= -A(0,1)*x[0]/A(1,1);                
+    x[2]= -(A(0,2)*x[0]+A(1,2)*x[1])/A(2,2);  
+    x[3]= -(A(0,3)*x[0]+A(1,3)*x[1]+A(2,3)*x[2])/A(3,3);
+
+    x[3] = x[3]/A(3,3);
+    x[2] = ( x[2] - A(2,3)*x[3] )/A(2,2);
+    x[1] = ( x[1] - A(1,2)*x[2] - A(1,3)*x[3])/A(1,1);
+    x[0] = ( x[0] - A(0,1)*x[1] - A(0,2)*x[2] - A(0,3)*x[3])/A(0,0);
+  } 
 }
 
 KOKKOS_INLINE_FUNCTION
 void QRDecomp4d(Matrix4d& A, Matrix4d& Q, Matrix4d& R){
   int m_size=4;
   for (int i=0; i<m_size; i++){
-      Vec4d A_column;
-      Vec4d u_column;
-      for (int l=0; l<m_size; l++){
-         A_column[l] = A(l,i);
-	 u_column[l] = A(l,i);
-      }
-
-      for (int j=0; j<i; j++){
-         Vec4d q_column;
-	 for (int l=0; l<m_size; l++)
-	   q_column[l] = Q(l,j);
-	 double r=q_column.dot(A_column);
-	 R(j,i) = r;
-         for (int k=0; k<m_size; k++) u_column[k] -= r*Q(k,j);
-      }
-      double uNorm = std::sqrt(u_column.dot(u_column));
-      R(i,i)=uNorm;
-      for (int k=0; k<m_size; k++)
-          Q(k,i)= u_column[k]/uNorm; 
+    Vec4d A_column;
+    Vec4d u_column;
+    for (int l=0; l<m_size; l++){
+      A_column[l] = A(l,i);
+      u_column[l] = A(l,i);
+    }
 
+    for (int j=0; j<i; j++){
+      Vec4d q_column;
+      for (int l=0; l<m_size; l++)
+        q_column[l] = Q(l,j);
+      double r=q_column.dot(A_column);
+      R(j,i) = r;
+      for (int k=0; k<m_size; k++) u_column[k] -= r*Q(k,j);
+    }
+    double uNorm = std::sqrt(u_column.dot(u_column));
+    R(i,i)=uNorm;
+    for (int k=0; k<m_size; k++)
+      Q(k,i)= u_column[k]/uNorm;
   }
-
 }
 
-
 KOKKOS_INLINE_FUNCTION
 void initArray(Vec3d* arr, int n, Vec3d fill){
-    for(int i=0; i<n; i++){
-        arr[i] = fill;
-    }
+  for(int i=0; i<n; i++){
+    arr[i] = fill;
+  }
 }
 
 KOKKOS_INLINE_FUNCTION
 void initArray(double* arr, int n, double fill){
-    for(int i=0; i<n; i++){
-        arr[i] = fill;
-    }
+  for(int i=0; i<n; i++){
+    arr[i] = fill;
+  }
 }
 
-
 KOKKOS_INLINE_FUNCTION
 double arcLength(Vec3d &a, Vec3d &b){
-    Vec3d c = b-a;
-    double r = a.magnitude();
-    return r * 2.0 * std::asin(c.magnitude() / (2.0 * r)) ;
+  Vec3d c = b-a;
+  double r = a.magnitude();
+  return r * 2.0 * std::asin(c.magnitude() / (2.0 * r)) ;
 }
 
 // implement from: https://github.com/MPAS-Seaice-MPM-Project/MPAS-Seaice-MPM/blob/9bb3b4f3a09ac59fbc93b8ab317ca0fa6a6f18bb/components/mpas-framework/src/operators/mpas_geometry_utils.F#L651-L671
 //TODO:check for more effcient calculation
 KOKKOS_INLINE_FUNCTION
 double sphericalTriangleArea(Vec3d& a, Vec3d& b, Vec3d& c, double radius){
-    double ab, bc, ca, semiperim, tanqe;
-    Vec3d ablen, aclen, dlen;
-
-    ab = arcLength(a,b) / radius;
-    bc = arcLength(b,c) / radius;
-    ca = arcLength(c,a) / radius;
-    semiperim = 0.5 * (ab + bc + ca);
-
-    tanqe = Kokkos::sqrt(Kokkos::fmax(0.0, Kokkos::tan(0.5 * semiperim) * 
-                                           Kokkos::tan(0.5 * (semiperim - ab)) *
-                                           Kokkos::tan(0.5 * (semiperim - bc)) * 
-                                           Kokkos::tan(0.5 * (semiperim - ca)))); 
-    return 4.0 * Kokkos::atan(tanqe);
+  double ab, bc, ca, semiperim, tanqe;
+  Vec3d ablen, aclen, dlen;
+
+  ab = arcLength(a,b) / radius;
+  bc = arcLength(b,c) / radius;
+  ca = arcLength(c,a) / radius;
+  semiperim = 0.5 * (ab + bc + ca);
+
+  tanqe = Kokkos::sqrt(Kokkos::fmax(0.0, Kokkos::tan(0.5 * semiperim) * 
+                                         Kokkos::tan(0.5 * (semiperim - ab)) *
+                                         Kokkos::tan(0.5 * (semiperim - bc)) * 
+                                         Kokkos::tan(0.5 * (semiperim - ca)))); 
+  return 4.0 * Kokkos::atan(tanqe);
 }
 
 //implement from: https://www.maa.org/sites/default/files/Eriksson14108673.pdf 
 KOKKOS_INLINE_FUNCTION
 double sphericalTriangleArea2(Vec3d &a, Vec3d &b, Vec3d &c, double radius){
-    double inv_radius = 1.0/radius;
-    Vec3d a_unit = a*inv_radius;
-    Vec3d b_unit = b*inv_radius;
-    Vec3d c_unit = c*inv_radius;
-    double tripleProduct = a_unit.dot(b_unit.cross(c_unit));
-    double tangent = tripleProduct / (1 + b_unit.dot(c_unit) + c_unit.dot(a_unit) + a_unit.dot(b_unit));
-
-    return 2.0 * Kokkos::atan(tangent);
+  double inv_radius = 1.0/radius;
+  Vec3d a_unit = a*inv_radius;
+  Vec3d b_unit = b*inv_radius;
+  Vec3d c_unit = c*inv_radius;
+  double tripleProduct = a_unit.dot(b_unit.cross(c_unit));
+  double tangent = tripleProduct / (1 + b_unit.dot(c_unit) + c_unit.dot(a_unit) + a_unit.dot(b_unit));
+
+  return 2.0 * Kokkos::atan(tangent);
 }
 
 //this is a lazy comparison and shouldn't be relied on beyond simple testing
diff --git a/src/pmpo_wachspressBasis.hpp b/src/pmpo_wachspressBasis.hpp
index cc3f49a..1dc262b 100644
--- a/src/pmpo_wachspressBasis.hpp
+++ b/src/pmpo_wachspressBasis.hpp
@@ -358,5 +358,74 @@ void sphericalInterpolation(MPMesh& mpMesh){
     pumipic::RecordTime("PolyMPO_sphericalInterpolation", timer.seconds());
 }
 
+
+inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
+    Kokkos::Timer timer;
+    auto p_mesh = mpMesh.p_mesh;
+    auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>();
+    int numVtxs = p_mesh->getNumVertices();
+    auto elm2VtxConn = p_mesh->getElm2VtxConn();
+    
+    auto p_MPs = mpMesh.p_MPs;
+    auto MPsPosition = p_MPs->getPositions();
+    double radius = p_mesh->getSphereRadius();
+    PMT_ALWAYS_ASSERT(radius > 0);
+ 
+    constexpr MeshFieldIndex meshFieldIndex1 = polyMPO::MeshF_RotLatLonIncr;
+    constexpr MeshFieldIndex meshFieldIndex2 = polyMPO::MeshF_OnSurfVeloIncr;
+    
+    auto meshField1 = p_mesh->getMeshField<meshFieldIndex1>();
+    auto meshField2 = p_mesh->getMeshField<meshFieldIndex2>();
+    
+    constexpr MaterialPointSlice mpfIndex1 = meshFieldIndexToMPSlice<meshFieldIndex1>;
+    constexpr MaterialPointSlice mpfIndex2 = meshFieldIndexToMPSlice<meshFieldIndex2>;
+
+    const int numEntries1 = mpSliceToNumEntries<mpfIndex1>();
+    const int numEntries2 = mpSliceToNumEntries<mpfIndex2>();
+
+    auto mpField1 = p_MPs->getData<mpfIndex1>();
+    auto mpField2 = p_MPs->getData<mpfIndex2>();
+
+    auto interpolation = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
+      if(mask) {
+        Vec3d position3d(MPsPosition(mp, 0), MPsPosition(mp, 1), MPsPosition(mp, 2));
+        Vec3d v3d[maxVtxsPerElm + 1];
+        int numVtx = elm2VtxConn(elm, 0);
+        for (int i = 1; i <= numVtx; i++) {
+          v3d[i-1][0] = vtxCoords(elm2VtxConn(elm, i) - 1, 0);
+          v3d[i-1][1] = vtxCoords(elm2VtxConn(elm, i) - 1, 1);
+          v3d[i-1][2] = vtxCoords(elm2VtxConn(elm, i) - 1, 2);
+        }
+        v3d[numVtx][0] = vtxCoords(elm2VtxConn(elm,1)-1,0);
+        v3d[numVtx][1] = vtxCoords(elm2VtxConn(elm,1)-1,1);
+        v3d[numVtx][2] = vtxCoords(elm2VtxConn(elm,1)-1,2);        
+
+        double basisByArea3d[maxVtxsPerElm] = {0.0};
+        initArray(basisByArea3d, maxVtxsPerElm, 0.0);
+
+        getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
+      
+        for(int entry=0; entry<numEntries1; entry++){
+          double mpValue = 0.0;
+          for(int i=1; i<= numVtx; i++){
+            mpValue += meshField1(elm2VtxConn(elm,i)-1,entry)*basisByArea3d[i-1];
+          }
+          mpField1(mp,entry) = mpValue;
+        }
+        
+        for(int entry=0; entry<numEntries2; entry++){
+          double mpValue = 0.0;
+          for(int i=1; i<= numVtx; i++){
+            mpValue += meshField2(elm2VtxConn(elm,i)-1,entry)*basisByArea3d[i-1];
+          }
+          mpField2(mp,entry) = mpValue;
+        }   
+      }
+    };
+    p_MPs->parallel_for(interpolation, "sphericalInterpolationMultiField");
+    pumipic::RecordTime("PolyMPO_sphericalInterpolationDispVelIncr", timer.seconds());
+  }
+
+
 } //namespace polyMPO end
 #endif
diff --git a/test/readMPAS.f90 b/test/readMPAS.f90
index ac7549d..c5a4116 100644
--- a/test/readMPAS.f90
+++ b/test/readMPAS.f90
@@ -27,7 +27,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
                         xVertex, yVertex, zVertex, &
                         latVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
     use :: netcdf
     use :: iso_c_binding
     implicit none
@@ -42,6 +42,8 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
     real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
+    real(kind=MPAS_RKIND), dimension(:), pointer, optional :: areaTriangle
+
 
     call polympo_checkPrecisionForRealKind(MPAS_RKIND)
     !check on maxEdges and vertexDegree
@@ -75,6 +77,11 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     
     !set mesh element center
     call polympo_setMeshElmCenter(mpMesh,nCells,c_loc(xCell),c_loc(yCell),c_loc(zCell))
+
+    if (present(areaTriangle)) then
+      call polympo_setMeshDualTriangleArea(mpMesh, nVertices, c_loc(areaTriangle)) 
+    end if
+
 end subroutine
 
 subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
@@ -83,7 +90,8 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, &
+                        areaTriangle)
     use :: netcdf
     use :: iso_c_binding
     implicit none
@@ -99,12 +107,14 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
     real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
+    
+    real(kind=MPAS_RKIND), dimension(:), pointer, optional :: areaTriangle(:)
 
     integer :: ncid, status, nCellsID, nVerticesID, maxEdgesID, vertexDegreeID, &
                nEdgesOnCellID, xVertexID, yVertexID, zVertexID, &
                latVertexID, lonVertexID, &
                xCellID, yCellID, zCellID, &
-               verticesOnCellID, cellsOnCellID
+               verticesOnCellID, cellsOnCellID, areaTriangleID
     
     status = nf90_open(path=trim(filename), mode=nf90_nowrite, ncid=ncid)
     if (status /= nf90_noerr) then
@@ -181,6 +191,10 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     allocate(verticesOnCell(maxEdges, nCells))
     allocate(cellsOnCell(maxEdges, nCells))
    
+    if (present(areaTriangle)) then
+      allocate(areaTriangle(nVertices))
+    end if
+   
     status = nf90_get_att(ncid, nf90_global, "on_a_sphere", onSphere)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on get attribute 'on_a_sphere'"
@@ -275,6 +289,16 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    if (present(areaTriangle)) then
+      status = nf90_inq_varid(ncid, 'areaTriangle', areaTriangleID)
+      if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'areaTriangleID'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+      end if
+    end if
+
+    !Reading the variables
     status = nf90_get_var(ncid, xVertexID, xVertex)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on get var of 'xVertex'"
@@ -352,6 +376,18 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    if (present(areaTriangle)) then
+      print *, "Reading area of dual triangles"
+      status = nf90_get_var(ncid, areaTriangleID, areaTriangle)
+      if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on get var of 'areaTriangle'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+      end if
+    else
+      print *, "Not reading area of dual cells"
+    end if
+
     status = nf90_close(ncid)   
 end subroutine readMPASMeshFromNCFile
 subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMPASMeshByFortran")
diff --git a/test/testFortranCreateRebuildMPs.f90 b/test/testFortranCreateRebuildMPs.f90
index 61ea24a..478708a 100644
--- a/test/testFortranCreateRebuildMPs.f90
+++ b/test/testFortranCreateRebuildMPs.f90
@@ -25,7 +25,7 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
     real(kind=MPAS_RKIND) :: ptOne = 0.1_MPAS_RKIND
     integer, parameter :: MP_ACTIVE = 1
     integer, parameter :: MP_INACTIVE = 0
-    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
     real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition
   
     isMPActive = MP_ACTIVE !no inactive MPs and some changed below
@@ -42,11 +42,13 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
     end do
   
     allocate(mpsPerElm(nCells))
+    allocate(globalElms(nCells))
     mpsPerElm = 1 !all elements have 1 MP and some changed below
     mpsPerElm(1) = 0 !1st element has 0 MPs
     mpsPerElm(2) = 2 !2nd element has 2 MPs
     mpsPerElm(3) = 2 !3rd element has 2 MPs 
 
+    call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
     call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
 
     !set mp positions
@@ -80,6 +82,7 @@ subroutine createMPsTest(mpMesh, nCells, numMPs, mp2Elm, isMPActive, mpPosition)
 
     !deallocate MP variables
     deallocate(mpsPerElm)
+    deallocate(globalElms)
 end subroutine
 
 subroutine rebuildMPsTests(mpMesh, numMPs, mp2Elm, isMPActive, mpPosition)
diff --git a/test/testFortranMPAdvection.f90 b/test/testFortranMPAdvection.f90
index abdedf7..292d806 100644
--- a/test/testFortranMPAdvection.f90
+++ b/test/testFortranMPAdvection.f90
@@ -218,7 +218,7 @@ program main
   real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
   integer :: numMPs, numMPsCount, numPush
-  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, mp2Elm_new
+  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, mp2Elm_new, globalElms
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon, mpPositions_new, mpLatLon_new
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
@@ -284,6 +284,7 @@ program main
 
   allocate(lonCell(nCells))
   allocate(mpsPerElm(nCells))
+  allocate(globalElms(nCells))
   allocate(mp2Elm(numMPs))
   allocate(mp2Elm_new(numMPs))
 
@@ -301,6 +302,7 @@ program main
     mp2Elm(numMPsCount+1:numMPsCount+localNumMPs) = i
     mpsPerElm(i) = localNumMPs
     numMPsCount = numMPsCount + localNumMPs
+    globalElms(i)=i
   end do
 
   call assert(numMPsCount == numMPs, "num mps miscounted")
@@ -356,24 +358,11 @@ program main
   end do
 
   call assert(numMPsCount == numMPs, "num mps miscounted")
-
+  call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
   call polympo_setMPPositions(mpMesh,3,numMPs,c_loc(mpPosition))
-
-  !Another advection test to test if material poins come back to the same position
-  call runAdvectionTest2(mpMesh, numPush, latVertex, lonVertex, nEdgesOnCell, verticesOnCell, nVertices, sphereRadius)
-  call polympo_getMPPositions(mpMesh, 3, numMPs, c_loc(mpPositions_new))
-  call polympo_getMPRotLatLon(mpMesh, 2, numMPs, c_loc(mpLatLon_new))
-    call polympo_getMPCurElmID(mpMesh, numMPS, c_loc(mp2Elm_new))
-
-  do i = 1, numMPs
-    if ( abs(mpLatLon_new(2,i)-mpLatLon(2,i)) > max_push_diff ) then
-      max_push_diff = abs(mpLatLon_new(2,i)-mpLatLon(2,i))
-    end if
-  end do
-  call assert(max_push_diff.le.TOLERANCE_PUSH , "MPs donot come back check push!")
-
+ 
   if (testType == "API") then
     call runApiTest(mpMesh, numMPs, nVertices, nCells, numPush, mpLatLon, mpPosition, xVertex, yVertex, zVertex, latVertex)
   else if (testType == "MIGRATION") then
@@ -412,7 +401,7 @@ program main
   deallocate(isMPActive)
   deallocate(mpPosition)
   deallocate(mpLatLon)
-
+  deallocate(globalElms)
   stop
 
 end program
diff --git a/test/testFortranMPAppIDs.f90 b/test/testFortranMPAppIDs.f90
index bcae9d7..a7cf7f2 100644
--- a/test/testFortranMPAppIDs.f90
+++ b/test/testFortranMPAppIDs.f90
@@ -40,7 +40,7 @@ subroutine testAppIDPointer(mpMesh) &
     integer :: setMeshOption, setMPOption
     integer :: ierr, self
     integer :: mpi_comm_handle = MPI_COMM_WORLD
-    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+    integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
     type(c_ptr) :: mpMesh
     integer :: nCells, numMPs, appID
     integer, parameter :: MP_ACTIVE = 1
@@ -56,15 +56,17 @@ subroutine testAppIDPointer(mpMesh) &
     setMeshOption = 1 !create a hard coded planar test mesh
     setMPOption = 0   !create an empty set of MPs
     mpMesh = polympo_createMPMesh(setMeshOption, setMPOption)
-
     numMPs = 1
     allocate(mpsPerElm(nCells))
+    allocate(globalElms(nCells))
     allocate(mp2Elm(numMPs))
     allocate(isMPActive(numMPs))
     call polympo_getMeshNumElms(mpMesh, nCells)
     mpsPerElm = 1
     mp2Elm = 1
     isMPActive = MP_ACTIVE
+    !This is a dummy array of global Cell IDs for each mesh element
+    call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
     call polympo_createMPs(mpMesh, nCells, numMPs, c_loc(mpsPerElm), c_loc(mp2Elm), c_loc(isMPActive))
     call polympo_setMPICommunicator(mpMesh, mpi_comm_handle)
     ! Set function and opaque data structure(list/queue) used to retrieve appIDS
@@ -74,9 +76,15 @@ subroutine testAppIDPointer(mpMesh) &
     call testAppIDPointer(mpMesh)
 
     ! Clean Up
+    deallocate(mpsPerElm)
+    deallocate(globalElms)
+    deallocate(mp2Elm)
+    deallocate(isMPActive)
+    
+    
     call polympo_deleteMPMesh(mpMesh)
     call queue_destroy(queue)
     call polympo_finalize()
     call mpi_finalize(ierr)
 
-end program
\ No newline at end of file
+end program
diff --git a/test/testFortranMPReconstruction.f90 b/test/testFortranMPReconstruction.f90
index e88720e..4a9e910 100644
--- a/test/testFortranMPReconstruction.f90
+++ b/test/testFortranMPReconstruction.f90
@@ -19,7 +19,7 @@ program main
   real(kind=MPAS_RKIND) :: pi = 4.0_MPAS_RKIND*atan(1.0_MPAS_RKIND)
   real(kind=MPAS_RKIND) :: TEST_VAL = 1.1_MPAS_RKIND
   real(kind=MPAS_RKIND) :: TOLERANCE = 0.0001_MPAS_RKIND
-  real(kind=MPAS_RKIND) :: TOLERANCE1 = 0.0000000001_MPAS_RKIND
+  real(kind=MPAS_RKIND) :: TOLERANCE1 = 0.000000001_MPAS_RKIND
   character (len=2048) :: filename
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: dispIncr
   character (len=64) :: onSphere
@@ -28,12 +28,15 @@ program main
   real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
   real(kind=MPAS_RKIND), dimension(:), pointer :: latVertex, lonVertex
   real(kind=MPAS_RKIND), dimension(:), pointer :: xCell, yCell, zCell
+  real(kind=MPAS_RKIND), dimension(:), pointer :: areaTriangle
+  
   integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
   integer :: numMPs, vID
-  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
+  integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive, globalElms
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition, mpLatLon
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpMass, mpVel
   real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxMass, meshElmMass, meshVtxMass1
+  real(kind=MPAS_RKIND), dimension(:), pointer :: meshVtxVelu, meshVtxVelv
   logical :: inBound
   integer, parameter :: MP_ACTIVE = 1
   integer, parameter :: MP_INACTIVE = 0
@@ -58,7 +61,7 @@ program main
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
   if (onSphere .ne. 'YES') then
     write (*,*) "The mesh is not spherical!"
     call exit(1)
@@ -74,13 +77,14 @@ program main
                         xVertex, yVertex, zVertex, &
                         latVertex, &
                         xCell, yCell, zCell, &
-                        verticesOnCell, cellsOnCell)
+                        verticesOnCell, cellsOnCell, areaTriangle)
  
   nCompsDisp = 2
   allocate(dispIncr(nCompsDisp,nVertices))
   !createMPs
   numMPs = nCells
   allocate(mpsPerElm(nCells))
+  allocate(globalElms(nCells))
   allocate(mp2Elm(numMPs))
   allocate(isMPActive(numMPs))
   allocate(mpPosition(3,numMPs))
@@ -90,6 +94,8 @@ program main
   allocate(meshVtxMass(nVertices))
   allocate(meshVtxMass1(nVertices))
   allocate(meshElmMass(nCells))
+  allocate(meshVtxVelu(nVertices))
+  allocate(meshVtxVelv(nVertices))
 
   isMPActive = MP_ACTIVE !all active MPs and some changed below
   mpsPerElm = 1 !all elements have 1 MP and some changed below
@@ -104,7 +110,8 @@ program main
     mpPosition(2,i) = yCell(i)
     mpPosition(3,i) = zCell(i)
   end do
-
+  
+  call polympo_setElmGlobal(mpMesh, nCells, c_loc(globalElms))
   call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
   call polympo_setMPICommunicator(mpMesh, mpi_comm_handle)
   call polympo_setMPRotLatLon(mpMesh,2,numMPs,c_loc(mpLatLon))
@@ -120,14 +127,20 @@ program main
 
   do i = 1, nVertices
     call assert(meshVtxMass(i) < TEST_VAL+TOLERANCE .and. meshVtxMass(i) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
+    !write(*, *) 'The value of LRV is:', meshVtxMass(i)
   end do
   
   !Test vtx order 1 reconstruction
   call polympo_setReconstructionOfMass(mpMesh,1,polympo_getMeshFVtxType())
+  call polympo_setReconstructionOfVel(mpMesh, 1, polympo_getMeshFVtxType())
   call polympo_applyReconstruction(mpMesh)
   call polympo_getMeshVtxMass(mpMesh,nVertices,c_loc(meshVtxMass1))
+  call polympo_getMeshVtxVel(mpMesh, nVertices, c_loc(meshVtxVelu), c_loc(meshVtxVelv))
   do i = 1, nVertices
-    call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass LINEAR")
+    call assert(meshVtxMass1(i) < TEST_VAL+TOLERANCE1 .and. meshVtxMass1(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass order 1")
+    call assert(meshVtxVelu(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelu(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velU order 1")
+    call assert(meshVtxVelv(i) < TEST_VAL+TOLERANCE1 .and. meshVtxVelv(i) > TEST_VAL-TOLERANCE1, "Error: wrong vtx velV order 1")
+    write(*, *) 'The value of LRV is:', meshVtxVelu(i)-TEST_VAL, meshVtxVelv(i)-TEST_VAL
   end do
 
 
@@ -143,7 +156,7 @@ program main
 
     do i = 1, numMPs
       vID = verticesOnCell(1,mp2Elm(i))
-      call assert(meshVtxMass(vID) < TEST_VAL+TOLERANCE .and. meshVtxMass(vID) > TEST_VAL-TOLERANCE, "Error: wrong vtx mass")
+      call assert(meshVtxMass(vID) < TEST_VAL+TOLERANCE1 .and. meshVtxMass(vID) > TEST_VAL-TOLERANCE1, "Error: wrong vtx mass push")
 
       call assert(meshElmMass(mp2Elm(i)) < TEST_VAL+TOLERANCE &
             .and. meshElmMass(mp2Elm(i)) > TEST_VAL-TOLERANCE, "Error: wrong elm mass")
@@ -175,7 +188,9 @@ program main
   deallocate(meshVtxMass)
   deallocate(meshVtxMass1)
   deallocate(meshElmMass)
-
+  deallocate(meshVtxVelu)
+  deallocate(meshVtxVelv)
+  deallocate(globalElms)
   stop
 
   contains
diff --git a/test/testMPAppIDs.cpp b/test/testMPAppIDs.cpp
index 4d33658..c3064cc 100644
--- a/test/testMPAppIDs.cpp
+++ b/test/testMPAppIDs.cpp
@@ -67,4 +67,4 @@ void testAppIDPointer(MPMesh_ptr p_mpmesh) {
   PMT_ALWAYS_ASSERT(numAddedMPsAfter_h == numAddedMPs);
 }
 
-#endif
\ No newline at end of file
+#endif