diff --git a/src/pmpo_MPMesh.cpp b/src/pmpo_MPMesh.cpp
index a5f709a..d525b04 100644
--- a/src/pmpo_MPMesh.cpp
+++ b/src/pmpo_MPMesh.cpp
@@ -2,7 +2,7 @@
 #include "pmpo_utils.hpp"
 #include "pmpo_MPMesh.hpp"
 #include "pmpo_wachspressBasis.hpp"
-
+#include <unistd.h>
 namespace polyMPO{
 
 void printVTP_mesh(MPMesh& mpMesh, int printVTPIndex=-1);
@@ -127,16 +127,18 @@ 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();
+
+    //Since Mesh is static print pnly for 1 time step 
     if(printVTPIndex>=0) {
       printVTP_mesh(printVTPIndex);
     }
-
+    
     Vec3dView history("positionHistory",numMPs);
     Vec3dView resultLeft("positionResult",numMPs);
     Vec3dView resultRight("positionResult",numMPs);
     Vec3dView mpTgtPosArray("positionTarget",numMPs);
-
+   
     auto CVTElmCalc = PS_LAMBDA(const int& elm, const int& mp, const int&mask){
         Vec3d MP(mpPositions(mp,0),mpPositions(mp,1),mpPositions(mp,2));
         if(mask){
@@ -155,7 +157,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
                 for(int i=1; i<=numConnElms; i++){
                     int elmID = elm2ElmConn(iElm,i)-1;
                     
-	            //New delta
+                    //New delta
 	            Vec3d center(elmCenter(elmID, 0), elmCenter(elmID, 1), elmCenter(elmID, 2));
                     delta = MPnew - center;
 
@@ -175,7 +177,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];
@@ -195,7 +197,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
         }
     };
     p_MPs->parallel_for(CVTElmCalc,"CVTTrackingElmCenterBasedCalc");
-
+    
     if(printVTPIndex>=0){
         Vec3dView::HostMirror h_history = Kokkos::create_mirror_view(history);
         Vec3dView::HostMirror h_resultLeft = Kokkos::create_mirror_view(resultLeft);
@@ -206,7 +208,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
         Kokkos::deep_copy(h_resultLeft, resultLeft);
         Kokkos::deep_copy(h_resultRight, resultRight);
         Kokkos::deep_copy(h_mpTgtPos, mpTgtPosArray);
-
         // printVTP file
         char* fileOutput = (char *)malloc(sizeof(char) * 256); 
         sprintf(fileOutput, "polyMPOCVTTrackingElmCenter_MPtracks_%d.vtp", printVTPIndex);
@@ -231,6 +232,7 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
         fprintf(pFile,"        </DataArray>\n      </Lines>\n    </Piece>\n  </PolyData>\n</VTKFile>\n");
         fclose(pFile);
     }
+
     pumipic::RecordTime("PolyMPO_CVTTrackingElmCenterBased", timer.seconds());
 }
 
@@ -323,11 +325,46 @@ bool getAnyIsMigrating(MaterialPoints* p_MPs, bool isMigrating) {
   return anyIsMigrating;
 }
 
+void MPMesh::push_ahead(){
+  //Latitude Longitude increment at mesh vertices and interpolate to particle position
+  p_mesh->computeRotLatLonIncr(); 
+  sphericalInterpolation<MeshF_RotLatLonIncr>(*this);
+  //Push the MPs
+  p_MPs->updateRotLatLonAndXYZ2Tgt(p_mesh->getSphereRadius());
+}
+
+bool MPMesh::push1P(){
+  //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());
+  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;
@@ -335,15 +372,19 @@ 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(); 
   } 
   while (anyIsMigrating);
-
+  
   pumipic::RecordTime("PolyMPO_push", timer.seconds());
 }
 
diff --git a/src/pmpo_MPMesh.hpp b/src/pmpo_MPMesh.hpp
index 33f0599..7b54afb 100644
--- a/src/pmpo_MPMesh.hpp
+++ b/src/pmpo_MPMesh.hpp
@@ -43,6 +43,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();
 
diff --git a/src/pmpo_MPMesh_assembly.hpp b/src/pmpo_MPMesh_assembly.hpp
index 25899b9..69e02e6 100644
--- a/src/pmpo_MPMesh_assembly.hpp
+++ b/src/pmpo_MPMesh_assembly.hpp
@@ -104,7 +104,7 @@ void MPMesh::resetPreComputeFlag(){
 }
 
 void MPMesh::computeMatricesAndSolve(){
-
+  Kokkos::Timer timer;
   //Mesh Information
   auto elm2VtxConn = p_mesh->getElm2VtxConn();  
   int numVtx = p_mesh->getNumVertices();
@@ -205,11 +205,12 @@ void MPMesh::computeMatricesAndSolve(){
       VtxCoeffs(vtx,i)=coeff[i];
   });
   this->precomputedVtxCoeffs = VtxCoeffs;
+  pumipic::RecordTime("PolyMPO_Calculate_MLS_Coeff", timer.seconds());
 }
 
 template <MeshFieldIndex meshFieldIndex>
 void MPMesh::assemblyVtx1() {
-   
+  Kokkos::Timer timer; 
   //If no reconstruction till now calculate the coeffs
   if (!isPreComputed) {
     computeMatricesAndSolve();
@@ -270,6 +271,7 @@ void MPMesh::assemblyVtx1() {
     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 67653b5..4ac7275 100644
--- a/src/pmpo_c.cpp
+++ b/src/pmpo_c.cpp
@@ -61,46 +61,56 @@ 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 nof of MPs which is >= no 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);
+  
+  //printf("With a MP, globally smallest mesh elm %d and first elm %d \n", globalMinElmID, globalFirstElmWithMPs);
+
   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");
     exit(1);
   }
+  printf("Offset %d \n", offset);
 
   std::vector<int> active_mpIDs(numMPs);
   std::vector<int> active_mp2Elm(numMPs);
@@ -112,20 +122,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,6 +198,45 @@ void polympo_startRebuildMPs_f(MPMesh_ptr p_mpmesh,
   }
 }
 
+void polympo_startRebuildMPs_f2(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) {
+  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);
+}
+
+
+
 void polympo_finishRebuildMPs_f(MPMesh_ptr p_mpmesh)
 {
   checkMPMeshValid(p_mpmesh);
@@ -200,13 +251,42 @@ 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){
+  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);
+
+  int rank;
+  auto mpi_comm=p_MPs->getMPIComm();
+  MPI_Comm_rank(mpi_comm, &rank);
+   
+  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);
+  
+}
+
 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();
@@ -214,13 +294,18 @@ void polympo_getMPCurElmID_f(MPMesh_ptr p_mpmesh,
   kkIntViewHostU arrayHost(elmIDs,numMPs);
   polyMPO::IntView mpCurElmIDCopy("mpCurElmIDNewValue",numMPs);
 
+  int rank;
+  auto mpi_comm=p_MPs->getMPIComm();
+  MPI_Comm_rank(mpi_comm, &rank);
+   
   auto getElmId = PS_LAMBDA(const int&, const int& mp, const int& mask){
     if(mask){
-        mpCurElmIDCopy(mpAppID(mp)) = mpCurElmID(mp)+elmIDoffset;
+        mpCurElmIDCopy(mpAppID(mp)) = mpCurElmID(mp)+elmIDoffset; 
     }
   };
   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 +324,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 +334,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 +356,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,6 +373,66 @@ 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());
+  //printf("NumMPs %d %d \n", 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_setMPPositions", timer.seconds());
+}
+
+void polympo_getMPTgtPositions_f(MPMesh_ptr p_mpmesh,
+                            const int nComps,
+                            const int numMPs,
+                            double* mpPositionsHost){
+  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);
+}
+
+
 void polympo_setMPRotLatLon_f(MPMesh_ptr p_mpmesh,
                          const int nComps,
                          const int numMPs,
@@ -338,6 +484,56 @@ 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){
+  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, "setMPRotLatLon");
+}
+
+void polympo_getMPTgtRotLatLon_f(MPMesh_ptr p_mpmesh,
+                         const int nComps,
+                         const int numMPs,
+                         double* mpRotLatLonHost){
+  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);
+}
+
+
 void polympo_setMPMass_f(MPMesh_ptr p_mpmesh, const int nComps, const int numMPs, const double* mpMassIn) {
   Kokkos::Timer timer;
   checkMPMeshValid(p_mpmesh);
@@ -965,6 +1161,27 @@ 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();
@@ -1029,6 +1246,27 @@ 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);
+}
+
+
+int polympo_getMPCount_f(MPMesh_ptr p_mpmesh) {
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  return p_MPs->getCount();  // This returns the int you're after
+}
+
 void polympo_enableTiming_f(){
   pumipic::EnableTiming();
 }
diff --git a/src/pmpo_c.h b/src/pmpo_c.h
index 174f862..ead18d6 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_startRebuildMPs_f2(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,8 @@ 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);
+
 
 //Mesh fields
 int polympo_getMeshFVtxType_f();
@@ -85,6 +101,10 @@ void polympo_getMeshDualTriangleArea_f(MPMesh_ptr p_mpmesh, const int nVertices,
 
 // 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);
diff --git a/src/pmpo_defines.h b/src/pmpo_defines.h
index 34857b9..14435eb 100644
--- a/src/pmpo_defines.h
+++ b/src/pmpo_defines.h
@@ -8,6 +8,7 @@
 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 cf013ca..694b6e0 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_startRebuildMPs_f2')
+    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
@@ -698,6 +790,21 @@ 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
+  
+  
   !---------------------------------------------------------------------------
   !> @brief calculate the MPs from given mesh vertices rotational latitude
   !>        longitude, update the MP slices
@@ -709,6 +816,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
@@ -789,6 +927,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
   !---------------------------------------------------------------------------
@@ -804,4 +949,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..46fb41a 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) = -1;
     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);
@@ -98,7 +114,7 @@ void MaterialPoints::finishRebuild() {
   updateMaxAppID();
   ps::destroyViews<MaterialPointTypes>(rebuildFields.addedSlices_h);
   ps::destroyViews<MaterialPointTypes>(addedSlices_d);
-  rebuildFields.ongoing = false;
+  rebuildFields.ongoing = false;  
 }
 
 MPI_Comm MaterialPoints::getMPIComm() {
@@ -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..70c71dc 100644
--- a/src/pmpo_materialPoints.hpp
+++ b/src/pmpo_materialPoints.hpp
@@ -18,6 +18,9 @@ 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,
@@ -130,15 +133,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 +162,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);
         }
diff --git a/src/pmpo_mesh.cpp b/src/pmpo_mesh.cpp
index 9232ad1..57b838f 100644
--- a/src/pmpo_mesh.cpp
+++ b/src/pmpo_mesh.cpp
@@ -77,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 2058e13..22e6ddd 100644
--- a/src/pmpo_mesh.hpp
+++ b/src/pmpo_mesh.hpp
@@ -88,7 +88,7 @@ class Mesh {
     IntVtx2ElmView elm2VtxConn_;
     IntElm2ElmView elm2ElmConn_;
     IntView owningProc_;
-
+    IntView globalElm_;
     //start of meshFields
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
@@ -163,6 +163,9 @@ class Mesh {
     void setOwningProc(IntView owningProc) {PMT_ALWAYS_ASSERT(meshEdit_);
                                             owningProc_ = owningProc; }
     
+    void setElmGlobal(IntView globalElm) {globalElm_ = globalElm;}
+    IntView getElmGlobal();
+
     void computeRotLatLonIncr();
 };
 
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 03a4165..4a9e910 100644
--- a/test/testFortranMPReconstruction.f90
+++ b/test/testFortranMPReconstruction.f90
@@ -32,7 +32,7 @@ program main
   
   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
@@ -84,6 +84,7 @@ program main
   !createMPs
   numMPs = nCells
   allocate(mpsPerElm(nCells))
+  allocate(globalElms(nCells))
   allocate(mp2Elm(numMPs))
   allocate(isMPActive(numMPs))
   allocate(mpPosition(3,numMPs))
@@ -109,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))
@@ -188,7 +190,7 @@ program main
   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