Checked gnomonic projection of MPs, vertices

Author: dhyan1272

src/pmpo_MPMesh_assembly.hpp

@@ -687,7 +687,7 @@ void MPMesh::startCommunication(){
   //Temporary
   double radius = p_mesh->getSphereRadius();
   bool isRotated=false;
-  //p_mesh->setGnomonicProjection(isRotated, radius);
+  p_mesh->setGnomonicProjection(isRotated, radius);
 
 
   if (p_MPs->getOpMode() != polyMPO::MP_DEBUG) 

src/pmpo_c.h

@@ -75,7 +75,6 @@ void polympo_setOwningProc_f(MPMesh_ptr p_mpmesh, const int nCells, const int* a
 void polympo_setOwningProcVertex_f(MPMesh_ptr p_mpmesh, const int nVertices, 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);
-void polympo_setGnomonicProjection(MPMesh_ptr p_mpmesh);
 
 //Mesh fields
 int polympo_getMeshFVtxType_f();

src/pmpo_mesh.cpp

@@ -44,6 +44,7 @@ namespace polyMPO{
 
         auto dualTriangleAreaEntry = meshFields2TypeAndString.at(MeshF_DualTriangleArea);
         PMT_ALWAYS_ASSERT(dualTriangleAreaEntry.first == MeshFType_VtxBased);
+        dualTriangleArea_ = MeshFView<MeshF_DualTriangleArea>(dualTriangleAreaEntry.second,numVtxs_);
     }
 
     void Mesh::setMeshElmBasedFieldSize(){
@@ -78,17 +79,6 @@ namespace polyMPO{
         pumipic::RecordTime("PolyMPO_computeRotLatLonIncr", timer.seconds());
     }
 
-    KOKKOS_INLINE_FUNCTION
-    void computeGnomonicProjectionAtPoint( const Vec3d& vtxCoord, const double gnomProjElmCenter[4], double& outX, double& outY) {
-      const double iDen = 1.0 / ( gnomProjElmCenter[1] * gnomProjElmCenter[3] * vtxCoord[0] +
-                                  gnomProjElmCenter[0] * gnomProjElmCenter[3] * vtxCoord[1] +
-                                  gnomProjElmCenter[2] * vtxCoord[2]);
-      outX = iDen * (vtxCoord[1] * gnomProjElmCenter[1] -
-                    vtxCoord[1] * gnomProjElmCenter[0]);
-      outY = iDen * (vtxCoord[2] * gnomProjElmCenter[3] - vtxCoord[1] * gnomProjElmCenter[2] * gnomProjElmCenter[0] -
-                     vtxCoord[0] * gnomProjElmCenter[1] * gnomProjElmCenter[2]);
-    }
-
     void Mesh::setGnomonicProjection(bool isRotated, double radius){
       std::cout<<__FUNCTION__<<std::endl;
       auto gnomProjVtx = getMeshField<MeshF_VtxGnomProj>();
@@ -114,6 +104,10 @@ namespace polyMPO{
         gnomProjElmCenter(iElm, 2) = invR*elmCenter[2];
         gnomProjElmCenter(iElm, 3) = invR*cosLatR;
 
+        if(iElm < 10)
+          printf("Elm %d, Centres %.15e %.15e %.15e %.15e \n", iElm, gnomProjElmCenter(iElm, 0), gnomProjElmCenter(iElm, 1),
+                                                                     gnomProjElmCenter(iElm, 2), gnomProjElmCenter(iElm, 3) );
+
         int nVtxE = elm2VtxConn(iElm,0);
         for(int i=0; i<nVtxE; i++){
           int vID = elm2VtxConn(iElm, i+1) - 1;
@@ -123,10 +117,14 @@ namespace polyMPO{
           }
 
           double outX, outY;
-          computeGnomonicProjectionAtPoint(vtxCord, &gnomProjElmCenter(iElm, 0), outX, outY);
+          auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, iElm, Kokkos::ALL);
+          computeGnomonicProjectionAtPoint(vtxCord, gnomProjElmCenter_sub, outX, outY);
 
           gnomProjVtx(iElm, i, 0) = outX;
           gnomProjVtx(iElm, i, 1) = outY;
+          
+          if(iElm < 10)
+            printf("Elm %d, Vtx %d -> outX: %.15e, outY: %.15e \n", iElm, i, outX, outY);
         }
       });
     }

src/pmpo_mesh.hpp

@@ -250,6 +250,23 @@ void Mesh::fillMeshField(int size, int numEntries, double val){
     });
 }
 
+KOKKOS_INLINE_FUNCTION
+void computeGnomonicProjectionAtPoint(const Vec3d& Coord, 
+     const Kokkos::View<double[4], Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnomProjElmCenter_sub, 
+     double& outX, double& outY){
+   /*
+   printf("Inline %.15e %.15e %.15e %.15e \n", gnomProjElmCenter_sub(0), gnomProjElmCenter_sub(1), gnomProjElmCenter_sub(2),
+          gnomProjElmCenter_sub(3));
+   */
+   const double iDen = 1.0 / (gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(3) * Coord[0] +
+                              gnomProjElmCenter_sub(0) * gnomProjElmCenter_sub(3) * Coord[1] +
+                              gnomProjElmCenter_sub(2) * Coord[2]);
+   outX = iDen * (Coord[1] * gnomProjElmCenter_sub(1) -
+                  Coord[0] * gnomProjElmCenter_sub(0));
+   outY = iDen * (Coord[2] * gnomProjElmCenter_sub(3) - Coord[1] * gnomProjElmCenter_sub(2) * gnomProjElmCenter_sub(0) -
+                  Coord[0] * gnomProjElmCenter_sub(1) * gnomProjElmCenter_sub(2));
+  }
+
 }
 
 #endif

src/pmpo_wachspressBasis.hpp

@@ -272,7 +272,6 @@ void getBasisByAreaGblFormSpherical2(Vec3d MP, int numVtxs, Vec3d* v,
     calcBasis(numVtxs, a, c, basis);
 }
 
-
 /*
 KOKKOS_INLINE_FUNCTION
 void getBasisByAreaGblForm_1(Vec2d MP, int numVtxs, Vec2d* vtxCoords, double* basis) {
@@ -359,6 +358,13 @@ void sphericalInterpolation(MPMesh& mpMesh){
 }
 
 
+KOKKOS_INLINE_FUNCTION
+void compute2DplanarTriangleArea(int numVtx, 
+     const Kokkos::View<double[maxVtxsPerElm][2], Kokkos::LayoutStride, Kokkos::MemoryTraits<Kokkos::Unmanaged>>& gnom_vtx_subview, 
+     double mpProjX, double mpProjY, double* basis){
+
+}
+
 inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
     Kokkos::Timer timer;
     auto p_mesh = mpMesh.p_mesh;
@@ -386,6 +392,10 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
     auto mpField1 = p_MPs->getData<mpfIndex1>();
     auto mpField2 = p_MPs->getData<mpfIndex2>();
 
+    // Field required for calculting gnomonic projection of MPs
+    auto gnomProjVtx = p_mesh->getMeshField<polyMPO::MeshF_VtxGnomProj>();
+    auto gnomProjElmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterGnomProj>();
+      
     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));
@@ -404,7 +414,15 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
         initArray(basisByArea3d, maxVtxsPerElm, 0.0);
 
         getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea3d);
-      
+         
+        //if using gnomonic Projection for weights
+        double mpProjX, mpProjY;
+        auto gnomProjElmCenter_sub = Kokkos::subview(gnomProjElmCenter, elm, Kokkos::ALL);
+        computeGnomonicProjectionAtPoint(position3d, gnomProjElmCenter_sub, mpProjX, mpProjY);
+        auto gnom_vtx_subview = Kokkos::subview(gnomProjVtx, elm, Kokkos::ALL, Kokkos::ALL); 
+        double basisbyArea2D[maxVtxsPerElm] = {0.0};
+        compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisbyArea2D);
+
         for(int entry=0; entry<numEntries1; entry++){
           double mpValue = 0.0;
           for(int i=1; i<= numVtx; i++){