Clean up and option to use 2D/3D basis functions

Author: dhyan1272

src/pmpo_mesh.cpp

@@ -104,10 +104,6 @@ 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;
@@ -122,9 +118,6 @@ namespace polyMPO{
 
           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_wachspressBasis.hpp

@@ -450,60 +450,56 @@ inline void sphericalInterpolationDispVelIncr(MPMesh& mpMesh){
     auto mpField2 = p_MPs->getData<mpfIndex2>();
 
     // Field required for calculting gnomonic projection of MPs
+    bool use3DArea=false;
     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));
-        Vec3d v3d[maxVtxsPerElm + 1];
+
+        double basisByArea[maxVtxsPerElm] = {0.0};
+        initArray(basisByArea, maxVtxsPerElm, 0.0);
         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);
+        Vec3d position3d(MPsPosition(mp, 0), MPsPosition(mp, 1), MPsPosition(mp, 2));
+
+        if(use3DArea){
+          Vec3d v3d[maxVtxsPerElm + 1];
+          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);        
+          getBasisByAreaGblFormSpherical(position3d, numVtx, v3d, radius, basisByArea);
+        }
+        else{ //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); 
+          compute2DplanarTriangleArea(numVtx, gnom_vtx_subview, mpProjX, mpProjY, basisByArea);
         }
-        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);
-         
-        //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++){
-            //mpValue += meshField1(elm2VtxConn(elm,i)-1,entry)*basisByArea3d[i-1];
-            mpValue += meshField1(elm2VtxConn(elm,i)-1,entry)*basisByArea2D[i-1];
-          }
+          for(int i=1; i<= numVtx; i++)
+            mpValue += meshField1(elm2VtxConn(elm,i)-1,entry)*basisByArea[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];
-            mpValue += meshField2(elm2VtxConn(elm,i)-1,entry)*basisByArea2D[i-1];
-          }
+          for(int i=1; i<= numVtx; i++)
+            mpValue += meshField2(elm2VtxConn(elm,i)-1,entry)*basisByArea[i-1];
           mpField2(mp,entry) = mpValue;
-        }   
-      }
-    };
+        } 
+      }//mask
+    };//lambda
     p_MPs->parallel_for(interpolation, "sphericalInterpolationMultiField");
     pumipic::RecordTime("PolyMPO_sphericalInterpolationDispVelIncr", timer.seconds());
-  }
-
+}
 
 } //namespace polyMPO end
 #endif