Dn/rec elm center

src/pmpo_MPMesh.cpp

@@ -116,7 +116,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
     int numElms = p_mesh->getNumElements();
     auto numMPs = p_MPs->getCount();
 
-    const auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>(); 
+    const auto elmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+
     auto elm2VtxConn = p_mesh->getElm2VtxConn();
     auto elm2ElmConn = p_mesh->getElm2ElmConn();
 
@@ -128,20 +129,6 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
       printVTP_mesh(printVTPIndex);
     }
 
-    Vec3dView elmCenter("elementCenter",numElms);
-    Kokkos::parallel_for("calcElementCenter", numElms, KOKKOS_LAMBDA(const int elm){  
-        int numVtx = elm2VtxConn(elm,0);
-        double sum_x = 0.0, sum_y = 0.0, sum_z = 0.0;
-        for(int i=1; i<= numVtx; i++){
-            sum_x += vtxCoords(elm2VtxConn(elm,i)-1,0);
-            sum_y += vtxCoords(elm2VtxConn(elm,i)-1,1);
-            sum_z += vtxCoords(elm2VtxConn(elm,i)-1,2);
-        }
-        elmCenter(elm)[0] = sum_x/numVtx;
-        elmCenter(elm)[1] = sum_y/numVtx;
-        elmCenter(elm)[2] = sum_z/numVtx;
-    });
-
     Vec3dView history("positionHistory",numMPs);
     Vec3dView resultLeft("positionResult",numMPs);
     Vec3dView resultRight("positionResult",numMPs);
@@ -155,19 +142,27 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
             Vec3d dx = MPnew-MP;
             while(true){
                 int numConnElms = elm2ElmConn(iElm,0);
-                Vec3d delta = MPnew - elmCenter(iElm);
+
+                Vec3d center(elmCenter(iElm, 0), elmCenter(iElm, 1), elmCenter(iElm, 2));
+                Vec3d delta = MPnew - center;
+
                 double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
                 int closestElm = -1;
                 //go through all the connected elm, calc distance
                 for(int i=1; i<=numConnElms; i++){
                     int elmID = elm2ElmConn(iElm,i)-1;
-                    delta = MPnew - elmCenter(elmID);
+
+	            //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];
                     if(neighborDistSq < minDistSq){
                         closestElm = elmID;
                         minDistSq = neighborDistSq;
                     }
                 }
+
                 if(closestElm<0){
                     MPs2Elm(mp) = iElm;
                     break;

src/pmpo_c.cpp

@@ -696,6 +696,43 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double
   }
 }
 
+void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, const double* xArray, const double* yArray, const double* zArray){
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); 
+
+  //copy the host array to the device
+  auto elmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+  auto h_elmCenter = Kokkos::create_mirror_view(elmCenter);
+  for(int i=0; i<nCells; i++){
+    h_elmCenter(i,0) = xArray[i];
+    h_elmCenter(i,1) = yArray[i];
+    h_elmCenter(i,2) = zArray[i];
+  }
+  Kokkos::deep_copy(elmCenter, h_elmCenter);
+}
+
+void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nCells, double* xArray, double* yArray, double* zArray){
+  //chech validity
+  checkMPMeshValid(p_mpmesh);
+  auto p_mesh = ((polyMPO::MPMesh*)p_mpmesh)->p_mesh;
+
+  //check the size
+  PMT_ALWAYS_ASSERT(p_mesh->getNumElements()==nCells); 
+
+  //copy the device to host 
+  auto elmCenter = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+  auto h_elmCenter = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), elmCenter);
+  for(int i=0; i<nCells; i++){
+    xArray[i] = h_elmCenter(i,0);
+    yArray[i] = h_elmCenter(i,1);
+    zArray[i] = h_elmCenter(i,2);
+  }
+}
+
 void polympo_setMeshVtxVel_f(MPMesh_ptr p_mpmesh, const int nVertices, const double* uVelIn, const double* vVelIn){
   //check mpMesh is valid
   checkMPMeshValid(p_mpmesh);

src/pmpo_c.h

@@ -74,6 +74,8 @@ void polympo_setMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, c
 void polympo_getMeshVtxOnSurfVeloIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, double* array);//vec2d
 void polympo_setMeshVtxOnSurfDispIncr_f(MPMesh_ptr p_mpmesh, const int nComps, const int nVertices, const double* array);//vec2d
 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);
 
 // calculations
 void polympo_push_f(MPMesh_ptr p_mpmesh);

src/pmpo_fortran.f90

@@ -487,6 +487,32 @@ subroutine polympo_getMeshVtxRotLat(mpMesh, nVertices, latitude) &
     type(c_ptr), value :: latitude
   end subroutine
   !---------------------------------------------------------------------------
+  !> @brief set the polympo mesh elements/cells center
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nElements(in) length of array in, use for assertion 
+  !> @param x/y/zArray(in) the 1D arrays of element centers coords 
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) &
+             bind(C, NAME='polympo_setMeshElmCenter_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nCells
+    type(c_ptr), intent(in), value :: xArray, yArray, zArray
+  end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief get the polympo mesh elements/cells center
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nElements(in) length of array in, use for assertion
+  !> @param x/y/zArray(in/out) the 1D arrays of element centers coords
+  !---------------------------------------------------------------------------
+  subroutine polympo_getMeshElmCenter(mpMesh, nCells, xArray, yArray, zArray) &
+             bind(C, NAME='polympo_getMeshElmCenter_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nCells
+    type(c_ptr), value :: xArray, yArray, zArray
+  end subroutine
+  !---------------------------------------------------------------------------
   !> @brief set the vertices velocity from a host array
   !> @param mpmesh(in/out) MPMesh object
   !> @param nVertices(in) numVertices

src/pmpo_mesh.cpp

@@ -48,9 +48,10 @@ namespace polyMPO{
 
         auto elmMassMapEntry = meshFields2TypeAndString.at(MeshF_ElmMass);
         PMT_ALWAYS_ASSERT(elmMassMapEntry.first == MeshFType_ElmBased);
-        elmMass_ = MeshFView<MeshF_ElmMass>(elmMassMapEntry.second,numElms_);
+        elmMass_      = MeshFView<MeshF_ElmMass>(elmMassMapEntry.second,numElms_);
+        elmCenterXYZ_ = MeshFView<MeshF_ElmCenterXYZ>(meshFields2TypeAndString.at(MeshF_ElmCenterXYZ).second,numElms_);
     }
-    
+
     void Mesh::computeRotLatLonIncr(){
         Kokkos::Timer timer;
         PMT_ALWAYS_ASSERT(geomType_ == geom_spherical_surf);

src/pmpo_mesh.hpp

@@ -19,6 +19,7 @@ enum MeshFieldIndex{
     MeshF_Unsupported,
     MeshF_VtxCoords,
     MeshF_VtxRotLat,
+    MeshF_ElmCenterXYZ,
     MeshF_Vel,
     MeshF_VtxMass,
     MeshF_ElmMass,
@@ -36,6 +37,7 @@ enum MeshFieldType{
 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_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*>; };
@@ -52,6 +54,7 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
                {MeshF_Unsupported,      {MeshFType_Unsupported,"MeshField_Unsupported"}},
                {MeshF_VtxCoords,        {MeshFType_VtxBased,"MeshField_VerticesCoords"}},
                {MeshF_VtxRotLat,        {MeshFType_VtxBased,"MeshField_VerticesLatitude"}},
+               {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}},
                {MeshF_Vel,              {MeshFType_VtxBased,"MeshField_Velocity"}},
                {MeshF_VtxMass,          {MeshFType_VtxBased,"MeshField_VerticesMass"}},
                {MeshF_ElmMass,          {MeshFType_ElmBased,"MeshField_ElementsMass"}},
@@ -85,6 +88,7 @@ class Mesh {
     //start of meshFields
     MeshFView<MeshF_VtxCoords> vtxCoords_;
     MeshFView<MeshF_VtxRotLat> vtxRotLat_;
+    MeshFView<MeshF_ElmCenterXYZ> elmCenterXYZ_;
     MeshFView<MeshF_Vel> vtxVel_;
     MeshFView<MeshF_VtxMass> vtxMass_;
     MeshFView<MeshF_ElmMass> elmMass_;
@@ -169,6 +173,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_VtxRotLat){
         return vtxRotLat_;
     }
+    else if constexpr (index==MeshF_ElmCenterXYZ){
+        return elmCenterXYZ_;
+    }
     else if constexpr (index==MeshF_Vel){
         return vtxVel_;
     }

test/readMPAS.f90

@@ -26,6 +26,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
                         onSphere, sphereRadius, &
                         xVertex, yVertex, zVertex, &
                         latVertex, &
+                        xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell)
     use :: netcdf
     use :: iso_c_binding
@@ -39,6 +40,7 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     integer, dimension(:), pointer :: nEdgesOnCell
     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
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
 
     call polympo_checkPrecisionForRealKind(MPAS_RKIND)
@@ -70,13 +72,17 @@ subroutine loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
     !set vtxCoords which is a mesh field 
     call polympo_setMeshVtxCoords(mpMesh,nVertices,c_loc(xVertex),c_loc(yVertex),c_loc(zVertex))
     call polympo_setMeshVtxRotLat(mpMesh,nVertices,c_loc(latVertex))
+
+    !set mesh element center
+    call polympo_setMeshElmCenter(mpMesh,nCells,c_loc(xCell),c_loc(yCell),c_loc(zCell))
 end subroutine
 
 subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
                         nCells, nVertices, nEdgesOnCell, &
                         onSphere, sphereRadius, &
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
+                        xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell)
     use :: netcdf
     use :: iso_c_binding
@@ -91,11 +97,13 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     integer, dimension(:), pointer :: nEdgesOnCell
     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
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
 
     integer :: ncid, status, nCellsID, nVerticesID, maxEdgesID, vertexDegreeID, &
                nEdgesOnCellID, xVertexID, yVertexID, zVertexID, &
                latVertexID, lonVertexID, &
+               xCellID, yCellID, zCellID, &
                verticesOnCellID, cellsOnCellID
 
     status = nf90_open(path=trim(filename), mode=nf90_nowrite, ncid=ncid)
@@ -166,6 +174,9 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
     allocate(zVertex(nVertices))
     allocate(latVertex(nVertices))
     allocate(lonVertex(nVertices))
+    allocate(xCell(nCells))
+    allocate(yCell(nCells))
+    allocate(zCell(nCells))
     allocate(nEdgesOnCell(nCells))
     allocate(verticesOnCell(maxEdges, nCells))
     allocate(cellsOnCell(maxEdges, nCells))
@@ -222,6 +233,27 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    status = nf90_inq_varid(ncid, 'xCell', xCellID)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'xCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
+    status = nf90_inq_varid(ncid, 'yCell', yCellID)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'yCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
+    status = nf90_inq_varid(ncid, 'zCell', zCellID)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'zCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
     status = nf90_inq_varid(ncid, 'nEdgesOnCell', nEdgesOnCellID)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on inquire varid of 'nEdgesOnCell'"
@@ -278,6 +310,27 @@ subroutine readMPASMeshFromNCFile(filename, maxEdges, vertexDegree, &
         stop
     end if
 
+    status = nf90_get_var(ncid, xCellID, xCell)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on get var of 'xCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
+    status = nf90_get_var(ncid, yCellID, yCell)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on get var of 'yCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
+    status = nf90_get_var(ncid, zCellID, zCell)
+    if (status /= nf90_noerr) then
+        write(0, *) "readMPASMeshFromNCFile: Error on get var of 'zCell'"
+        write(0, *) trim(nf90_strerror(status))
+        stop
+    end if
+
     status = nf90_get_var(ncid, nEdgesOnCellID, nEdgesOnCell)
     if (status /= nf90_noerr) then
         write(0, *) "readMPASMeshFromNCFile: Error on get var of 'nEdgesOnCell'"
@@ -315,6 +368,7 @@ subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMP
     integer, dimension(:), pointer :: nEdgesOnCell
     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
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
 
     fileNameFortran = transfer(fileName(1:n), fileNameFortran) 
@@ -325,18 +379,23 @@ subroutine setWithMPASMeshByFortran(mpMesh, fileName, n) bind(C, name="setWithMP
                         onSphere, sphereRadius, &
                         xVertex, yVertex, zVertex, &
                         latVertex, lonVertex, &
+                        xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell)
     call loadMPASMeshInPolyMPO(mpMesh, maxEdges, vertexDegree, &
                         nCells, nVertices, nEdgesOnCell, &
                         onSphere, sphereRadius, &
                         xVertex, yVertex, zVertex, &
                         latVertex, &
+                        xCell, yCell, zCell, &
                         verticesOnCell, cellsOnCell)
 
     deallocate(nEdgesOnCell)
     deallocate(xVertex)
     deallocate(yVertex)
     deallocate(zVertex)
+    deallocate(xCell)
+    deallocate(yCell)
+    deallocate(zCell)
     deallocate(verticesOnCell)
     deallocate(cellsOnCell)
 end subroutine

test/testFortran.f90

@@ -168,6 +168,26 @@ program main
   deallocate(yArray)
   deallocate(zArray)
 
+  ! test elmcenter
+  allocate(xArray(numElms))
+  allocate(yArray(numElms))
+  allocate(zArray(numElms))
+  xArray = value1
+  yArray = value2
+  zArray = value1 + value2 
+  call polympo_setMeshElmCenter(mpMesh, numElms, c_loc(xArray), c_loc(yArray), c_loc(zArray))
+  xArray = -1
+  yArray = -1
+  zArray = -1
+  call polympo_getMeshElmCenter(mpMesh, numElms, c_loc(xArray), c_loc(yArray), c_loc(zArray))
+  call assert(all(xArray .eq. value1), "Assert xArray == value1 Failed!")
+  call assert(all(yArray .eq. value2), "Assert yArray == value2 Failed!")
+  call assert(all(zArray .eq. value1 + value2), "Assert zArray == value1 + value2 Failed!")
+
+  deallocate(xArray)
+  deallocate(yArray)
+  deallocate(zArray)
+
   call polympo_deleteMPMesh(mpMesh)
   call polympo_finalize()