Removed old routines that does reconstruction assembly in MPAS

Author: dhyan1272

src/pmpo_MPMesh.cpp

@@ -18,7 +18,8 @@ void MPMesh::calcBasis(bool use3DArea) {
     //For Gnomonic Projection
     auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
     auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
-     
+
+    bool isRotated = true;    
     auto calcbasis = PS_LAMBDA(const int& elm, const int& mp, const int& mask) {
         if(mask) { //if material point is 'active'/'enabled'
             Vec3d position3d(MPsPosition(mp,0),MPsPosition(mp,1),MPsPosition(mp,2));
@@ -40,6 +41,10 @@ void MPMesh::calcBasis(bool use3DArea) {
             if(!use3DArea){
               double mpProjX, mpProjY;
               auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
+              if(isRotated){
+                position3d[0] = -MPsPosition(mp, 2);
+                position3d[2] = MPsPosition(mp, 0);
+              }
               computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
               auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL); 
               compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);

src/pmpo_MPMesh_assembly.hpp

@@ -1249,6 +1249,7 @@ void polympo_setReconstructionOfStress_f(MPMesh_ptr p_mpmesh, const int order, c
   (void)meshEntType;
 }
 
+//TO DO DELETE
 //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);
@@ -1299,7 +1300,6 @@ void polympo_applyReconstruction_f(MPMesh_ptr p_mpmesh){
   mpmesh->reconstructSlices();
 }
 
-//Method 1
 void polympo_fullAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCells, int nVerticesSolve, int nVertices, 
                                    double* array_sub, double* array_full){
   checkMPMeshValid(p_mpmesh);
@@ -1308,18 +1308,27 @@ void polympo_fullAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCell
   PMT_ALWAYS_ASSERT(nVertices == p_mesh->getNumVertices());
   mpmesh->assembleField(vtxPerElm, nCells, nVerticesSolve, nVertices, array_sub, array_full);
 }
+//
 
+//Only these are needed
 void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
   mpmesh->reconstruct_full<polyMPO::MeshF_VtxMass>();
 }
 
+void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh){
+  checkMPMeshValid(p_mpmesh);
+  auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
+  mpmesh->reconstruct_full<polyMPO::MeshF_Vel>();
+}
+
 void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh){
   checkMPMeshValid(p_mpmesh);
   auto mpmesh = ((polyMPO::MPMesh*)p_mpmesh);
   mpmesh->reconstruct_coeff_full();
 }
+//Till here
 
 void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* array){
   checkMPMeshValid(p_mpmesh);

src/pmpo_c.cpp

@@ -130,8 +130,8 @@ void polympo_fullAssemblyIceArea_f(MPMesh_ptr p_mpmesh, int vtxPerElm, int nCell
                                    double* array_sub, double* array_full);
 
 void polympo_reconstruct_coeff_with_MPI_f(MPMesh_ptr p_mpmesh);
-
 void polympo_reconstruct_iceArea_with_MPI_f(MPMesh_ptr p_mpmesh);
+void polympo_reconstruct_velocity_with_MPI_f(MPMesh_ptr p_mpmesh);
 
 // Timing
 void polympo_enableTiming_f();

src/pmpo_c.h

@@ -977,7 +977,14 @@ subroutine polympo_reconstruct_iceArea_with_MPI(mpMesh) &
     use :: iso_c_binding
     type(c_ptr), value :: mpMesh
   end subroutine
+  
+  subroutine polympo_reconstruct_velocity_with_MPI(mpMesh) &
+             bind(C, NAME='polympo_reconstruct_velocity_with_MPI_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+  end subroutine
 
+
   !---------------------------------------------------------------------------
   !> @brief directly call the reconstruct of the MP fields to mesh fields
   !> @param mpmesh(in/out) MPMesh object

src/pmpo_fortran.f90

@@ -79,7 +79,7 @@ namespace polyMPO{
         pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds());
     }
 
-    void Mesh::setGnomonicProjection(bool isRotated, double radius){
+    void Mesh::setGnomonicProjection(bool isRotated){
       std::cout<<__FUNCTION__<<std::endl;
       auto gnomProjVtx = getMeshField<MeshF_VtxGnomProj>();
       auto gnomProjElmCenter = getMeshField<MeshF_ElmCenterGnomProj>();
@@ -93,7 +93,8 @@ namespace polyMPO{
 
         Vec3d elmCenter(elmCenters(iElm, 0), elmCenters(iElm, 1), elmCenters(iElm, 2));
         if(isRotated){
-        
+          elmCenter[0] = - elmCenters(iElm, 2);
+          elmCenter[2] =   elmCenters(iElm, 0);
         }
         auto cos2LatR = elmCenter[0]*elmCenter[0] + elmCenter[1]*elmCenter[1];
         auto invR = 1.0/ sqrt(cos2LatR + elmCenter[2]*elmCenter[2]);
@@ -108,8 +109,9 @@ namespace polyMPO{
         for(int i=0; i<nVtxE; i++){
           int vID = elm2VtxConn(iElm, i+1) - 1;
           Vec3d vtxCord(vtxCoords(vID, 0), vtxCoords(vID, 1), vtxCoords(vID, 2));
-          if (isRotated){
-
+          if(isRotated){
+            vtxCord[0] = - vtxCoords(vID, 2);
+            vtxCord[2] = vtxCoords(vID, 0);
           }
 
           double outX, outY;

src/pmpo_mesh.cpp

@@ -186,7 +186,7 @@ class Mesh {
     IntView getElmGlobal() {return globalElm_;}
     IntView getVtxGlobal() {return globalVtx_;}
 
-    void setGnomonicProjection(bool isRotated, double radius);
+    void setGnomonicProjection(bool isRotated);
 
     void computeRotLatLonIncr();
 };

src/pmpo_mesh.hpp

@@ -450,7 +450,8 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
     auto mpField2 = p_MPs->getData<mpfIndex2>();
 
     // Field required for calculting gnomonic projection of MPs
-    bool use3DArea=false;
+    bool use3DArea = false;
+    bool isRotated = true;
     auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
     auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
 
@@ -477,6 +478,10 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
         else{ //if using gnomonic Projection for weights
           double mpProjX, mpProjY;
           auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
+          if(isRotated){
+            position3d[0] = -MPsPosition(mp, 2);
+            position3d[2] = MPsPosition(mp, 0);
+          }
           computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
           auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL); 
           compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);