Merge yg/elem_centre to ac/reconstrcuction

src/pmpo_MPMesh.cpp

@@ -117,6 +117,8 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
     auto numMPs = p_MPs->getCount();
 
     const auto vtxCoords = p_mesh->getMeshField<polyMPO::MeshF_VtxCoords>(); 
+    const auto elCenters = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+
     auto elm2VtxConn = p_mesh->getElm2VtxConn();
     auto elm2ElmConn = p_mesh->getElm2ElmConn();
 
@@ -156,18 +158,36 @@ void MPMesh::CVTTrackingElmCenterBased(const int printVTPIndex){
             while(true){
                 int numConnElms = elm2ElmConn(iElm,0);
                 Vec3d delta = MPnew - elmCenter(iElm);
-                double minDistSq = delta[0]*delta[0] + delta[1]*delta[1] + delta[2]*delta[2];
+                //New delta 
+	       	Vec3d center(elCenters(iElm, 0), elCenters(iElm, 1), elCenters(iElm, 2));
+		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(elCenters(elmID, 0), elCenters(elmID, 1), elCenters(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(elm==2095){
+                        printf("elm %d Loop %d Adj elm %d Mins %.15e %.15e=>%.15e %.15e %.15e C new: %.15e %.15e %.15e\n",
+	                  elm, i, elmID, neighborDistSq, minDistSq, mpTgtPos(mp,0), mpTgtPos(mp,1), mpTgtPos(mp,2), 
+			  elCenters(elmID, 0), elCenters(elmID, 1), elCenters(elmID, 2));
+                    }
+		    */
                 }
+
+		//if(elm==2095) printf("Starting %d Ending %d %d \n", elm, iElm, closestElm);
+
                 if(closestElm<0){
                     MPs2Elm(mp) = iElm;
                     break;

src/pmpo_c.cpp

@@ -696,6 +696,44 @@ void polympo_getMeshVtxRotLat_f(MPMesh_ptr p_mpmesh, const int nVertices, double
   }
 }
 
+void polympo_setMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, 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()==nElements); 
+
+  //copy the host array to the device
+  auto coordsArray = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+  auto h_coordsArray = Kokkos::create_mirror_view(coordsArray);
+  for(int i=0; i<nElements; i++){
+    h_coordsArray(i,0) = xArray[i];
+    h_coordsArray(i,1) = yArray[i];
+    h_coordsArray(i,2) = zArray[i];
+  }
+  Kokkos::deep_copy(coordsArray, h_coordsArray);
+}
+
+void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, 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()==nElements); 
+
+  //copy the device to host 
+  auto coordsArray = p_mesh->getMeshField<polyMPO::MeshF_ElmCenterXYZ>();
+  auto h_coordsArray = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(),
+                                                           coordsArray);
+  for(int i=0; i<nElements; i++){
+    xArray[i] = h_coordsArray(i,0);
+    yArray[i] = h_coordsArray(i,1);
+    zArray[i] = h_coordsArray(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 nElements, const double* xArray, const double* yArray, const double* zArray);
+void polympo_getMeshElmCenter_f(MPMesh_ptr p_mpmesh, const int nElements, 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, nElements, xArray, yArray, zArray) &
+             bind(C, NAME='polympo_setMeshElmCenter_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nElements
+    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, nElements, xArray, yArray, zArray) &
+             bind(C, NAME='polympo_getMeshElmCenter_f')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nElements
+    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

@@ -24,7 +24,8 @@ enum MeshFieldIndex{
     MeshF_ElmMass,
     MeshF_OnSurfVeloIncr,
     MeshF_OnSurfDispIncr,
-    MeshF_RotLatLonIncr
+    MeshF_RotLatLonIncr,
+    MeshF_ElmCenterXYZ
 };
 enum MeshFieldType{
     MeshFType_Invalid = -2,
@@ -42,6 +43,7 @@ template <> struct meshFieldToType < MeshF_ElmMass           > { using type = Ko
 template <> struct meshFieldToType < MeshF_OnSurfVeloIncr    > { using type = Kokkos::View<vec2d_t*>; };
 template <> struct meshFieldToType < MeshF_OnSurfDispIncr    > { using type = Kokkos::View<vec2d_t*>; };
 template <> struct meshFieldToType < MeshF_RotLatLonIncr     > { using type = Kokkos::View<vec2d_t*>; };
+template <> struct meshFieldToType < MeshF_ElmCenterXYZ      > { using type = Kokkos::View<vec3d_t*>; };
 
 template <MeshFieldIndex index>
 using MeshFView = typename meshFieldToType<index>::type;
@@ -57,7 +59,8 @@ const std::map<MeshFieldIndex, std::pair<MeshFieldType,
                {MeshF_ElmMass,          {MeshFType_ElmBased,"MeshField_ElementsMass"}},
                {MeshF_OnSurfVeloIncr,   {MeshFType_VtxBased,"MeshField_OnSurfaceVelocityIncrement"}},
                {MeshF_OnSurfDispIncr,   {MeshFType_VtxBased,"MeshField_OnSurfaceDisplacementIncrement"}},
-               {MeshF_RotLatLonIncr,    {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}}};
+               {MeshF_RotLatLonIncr,    {MeshFType_VtxBased,"MeshField_RotationalLatitudeLongitudeIncreasement"}},
+               {MeshF_ElmCenterXYZ,     {MeshFType_ElmBased,"MeshField_ElementCenterXYZ"}}};
 
 enum mesh_type {mesh_unrecognized_lower = -1,
                 mesh_general_polygonal, //other meshes
@@ -91,6 +94,7 @@ class Mesh {
     MeshFView<MeshF_OnSurfVeloIncr> vtxOnSurfVeloIncr_;
     MeshFView<MeshF_OnSurfDispIncr> vtxOnSurfDispIncr_;
     MeshFView<MeshF_RotLatLonIncr> vtxRotLatLonIncr_;
+    MeshFView<MeshF_ElmCenterXYZ>elmCenterXYZ_;
     //DoubleMat2DView vtxStress_;
 
   public:
@@ -187,6 +191,9 @@ auto Mesh::getMeshField(){
     else if constexpr (index==MeshF_RotLatLonIncr){
         return vtxRotLatLonIncr_;
     }
+    else if constexpr (index==MeshF_ElmCenterXYZ){
+        return elmCenterXYZ_;
+    }
     fprintf(stderr,"Mesh Field Index error!\n");
     exit(1);
 }

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()
 

test/testFortranCreateRebuildMPs.f90

@@ -254,6 +254,7 @@ program main
   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
   real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition
 
@@ -279,12 +280,14 @@ program main
                         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)
 
   !check for allocation
@@ -311,6 +314,9 @@ program main
   deallocate(xVertex)
   deallocate(yVertex)
   deallocate(zVertex)
+  deallocate(xCell)
+  deallocate(yCell)
+  deallocate(zCell)
   deallocate(verticesOnCell)
   deallocate(cellsOnCell)