support getting and setting the positions of material points

src/pmpo_c.cpp

@@ -77,7 +77,7 @@ typedef Kokkos::View<
           Kokkos::DefaultHostExecutionSpace,
           Kokkos::MemoryTraits<Kokkos::Unmanaged>
         > kkVec2dViewHostU;//TODO:put it to mesh.hpp
-                         
+
 typedef Kokkos::View<
           double**,
           Kokkos::LayoutLeft,
@@ -191,26 +191,56 @@ void polympo_getMPCurElmID(MPMesh_ptr p_mpmesh,
   Kokkos::deep_copy( arrayHost, mpCurElmIDCopy);
 }
 
+void polympo_setMPPositions(MPMesh_ptr p_mpmesh,
+                            int numComps,
+                            int numMPs,
+                            double* mpPositionsIn){
+  static int callCount = 0;
+  PMT_ALWAYS_ASSERT(callCount == 0);
+  checkMPMeshValid(p_mpmesh);
+  auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
+  PMT_ALWAYS_ASSERT(numComps == vec3d_nEntries);
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
+
+  auto mpPositions = p_MPs->getData<polyMPO::MPF_Cur_Pos_XYZ>();
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
+  kkDbl2dViewHostU mpPositionsIn_h(mpPositionsIn,numComps,numMPs);
+  Kokkos::View<double**> mpPositionsIn_d("mpPositionsDevice",vec3d_nEntries,numMPs);
+  Kokkos::deep_copy(mpPositionsIn_d, mpPositionsIn_h);
+  auto setPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+    if(mask){
+      mpPositions(mp,0) = mpPositionsIn_d(0, mpAppID(mp));
+      mpPositions(mp,1) = mpPositionsIn_d(1, mpAppID(mp));
+      mpPositions(mp,2) = mpPositionsIn_d(2, mpAppID(mp));
+    }
+  };
+  p_MPs->parallel_for(setPos, "setMPPositions");
+  callCount++;
+}
+
 void polympo_getMPPositions(MPMesh_ptr p_mpmesh,
-                           int numComps,
-                           int numMPs,
-                           double* mpPositionsIn){ 
+                            int numComps,
+                            int numMPs,
+                            double* mpPositionsHost){
   checkMPMeshValid(p_mpmesh);
   auto p_MPs = ((polyMPO::MPMesh*)p_mpmesh)->p_MPs;
   PMT_ALWAYS_ASSERT(numComps == vec3d_nEntries);
-  PMT_ALWAYS_ASSERT(numMPs == p_MPs->getCount());
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getCount());
+  PMT_ALWAYS_ASSERT(numMPs >= p_MPs->getMaxAppID());
 
   auto mpPositions = p_MPs->getData<polyMPO::MPF_Cur_Pos_XYZ>();
-  kkDbl2dViewHostU arrayHost(mpPositionsIn,numComps,numMPs);
+  auto mpAppID = p_MPs->getData<polyMPO::MPF_MP_APP_ID>();
   Kokkos::View<double**> mpPositionsCopy("mpPositionsCopy",vec3d_nEntries,numMPs);
-  auto setVel = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
+  auto getPos = PS_LAMBDA(const int& elm, const int& mp, const int& mask){
     if(mask){
-        mpPositionsCopy(0,mp) = mpPositions(mp,0);
-        mpPositionsCopy(1,mp) = mpPositions(mp,1);
-        mpPositionsCopy(2,mp) = mpPositions(mp,2);
+      mpPositionsCopy(0,mpAppID(mp)) = mpPositions(mp,0);
+      mpPositionsCopy(1,mpAppID(mp)) = mpPositions(mp,1);
+      mpPositionsCopy(2,mpAppID(mp)) = mpPositions(mp,2);
     }
   };
-  p_MPs->parallel_for(setVel, "get mpCurElmID");
+  p_MPs->parallel_for(getPos, "getMPPositions");
+  kkDbl2dViewHostU arrayHost(mpPositionsHost,numComps,numMPs);
   Kokkos::deep_copy(arrayHost, mpPositionsCopy);
 }
 

src/pmpo_c.h

@@ -22,9 +22,10 @@ void polympo_setMPICommunicator(MPI_Fint fcomm);//TODO:is MPI_Fint best? or some
 //MP info
 void polympo_createMPs(MPMesh_ptr p_mpmesh, int numElms, int numMPs, int* mpsPerElm, int* mp2Elm, int* isMPActive);
 void polympo_getMPCurElmID(MPMesh_ptr p_mpmesh, int numMPs, int* elmIDs);
-void polympo_getMPPositions(MPMesh_ptr p_mpmesh, int numComps, int numMPs, double* mpPositionsIn);
 
 //MP slices
+void polympo_getMPPositions(MPMesh_ptr p_mpmesh, int numComps, int numMPs, double* mpPositionsIn);
+void polympo_setMPPositions(MPMesh_ptr p_mpmesh, int numComps, int numMPs, double* mpPositionsIn);
 void polympo_setMPVel(MPMesh_ptr p_mpmesh, int size, double* array);
 void polympo_getMPVel(MPMesh_ptr p_mpmesh, int size, double* array);
 

src/pmpo_fortran.f90

@@ -98,6 +98,21 @@ subroutine polympo_getMPPositions(mpMesh, nComps, numMPs, array) &
     integer(c_int), value :: nComps, numMPs
     type(c_ptr), value :: array
   end subroutine
+  !---------------------------------------------------------------------------
+  !> @brief set the MP positions array from a host array
+  !> @param mpmesh(in/out) MPMesh object
+  !> @param nComps(in) number of components, should always be 3
+  !> @param numMPs(in) number of the MPs
+  !> @param array(in) MP current position 2D array (3,numMPs), allocated by user on host
+  !---------------------------------------------------------------------------
+  subroutine polympo_setMPPositions(mpMesh, nComps, numMPs, array) &
+             bind(C, NAME='polympo_setMPPositions')
+    use :: iso_c_binding
+    type(c_ptr), value :: mpMesh
+    integer(c_int), value :: nComps, numMPs
+    type(c_ptr), value :: array
+  end subroutine
+
   !---------------------------------------------------------------------------
   !> @brief set the velocity MP array from a host array
   !> @warning THIS IS NOT SUPPORTED YET 

src/pmpo_materialPoints.cpp

@@ -8,12 +8,14 @@ PS* createDPS(int numElms, int numMPs, DoubleVec3dView positions, IntView mpsPer
   auto mpPositions = ps::getMemberView<MaterialPointTypes, MPF_Cur_Pos_XYZ>(mpInfo);
   auto mpCurElmPos = ps::getMemberView<MaterialPointTypes, MPF_Cur_Elm_ID>(mpInfo);
   auto mpStatus = ps::getMemberView<MaterialPointTypes, MPF_Status>(mpInfo);
+  auto mpAppID_m = ps::getMemberView<MaterialPointTypes, MPF_MP_APP_ID>(mpInfo);
   Kokkos::parallel_for("setMPinfo", numMPs, KOKKOS_LAMBDA(int i) {
     mpPositions(i,0) = positions(i,0);
     mpPositions(i,1) = positions(i,1);
     mpPositions(i,2) = positions(i,2);
     mpCurElmPos(i) = mp2elm(i);
     mpStatus(i) = 1; 
+    mpAppID_m(i) = i;
   });
   Kokkos::TeamPolicy<Kokkos::DefaultExecutionSpace> policy(numElms,Kokkos::AUTO);
   auto dps = new DPS<MaterialPointTypes>(policy, numElms, numMPs, mpsPerElm, elmGids, mp2elm, mpInfo);
@@ -38,4 +40,15 @@ PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntVie
   return dps;
 }
 
+int getMaxAppID(IntView mpAppID) {
+  int maxAppID = 0;
+  Kokkos::parallel_reduce("getMax" , mpAppID.size(),
+    KOKKOS_LAMBDA(const int i, int & valueToUpdate) {
+      if ( mpAppID(i) > valueToUpdate ) valueToUpdate = mpAppID(i) ;
+    },
+    Kokkos::Max<int>(maxAppID)
+  );
+  return maxAppID;
+}
+
 }

src/pmpo_materialPoints.hpp

@@ -95,20 +95,24 @@ typedef ps::ParticleStructure<MaterialPointTypes> PS;
 
 
 PS* createDPS(int numElms, int numMPs, DoubleVec3dView positions, IntView mpsPerElm, IntView mp2elm);
-PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView isMPActive);
+PS* createDPS(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID);
+int getMaxAppID(IntView mpAppID);
 
 class MaterialPoints {
   private:
     PS* MPs;
     int elmIDoffset = -1;
+    int maxAppID = -1;
 
   public:
     MaterialPoints() : MPs(nullptr) {};
     MaterialPoints(int numElms, int numMPs, DoubleVec3dView positions, IntView mpsPerElm, IntView mp2elm) {
       MPs = createDPS(numElms, numMPs, positions, mpsPerElm, mp2elm);
+      maxAppID = numMPs; //this ctor does not support inactive MPs
     };
-    MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView isMPActive) {
-      MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, isMPActive);
+    MaterialPoints(int numElms, int numMPs, IntView mpsPerElm, IntView mp2elm, IntView mpAppID) {
+      MPs = createDPS(numElms, numMPs, mpsPerElm, mp2elm, mpAppID);
+      maxAppID = polyMPO::getMaxAppID(mpAppID);
     };
     ~MaterialPoints() {
       if(MPs != nullptr)
@@ -179,6 +183,10 @@ class MaterialPoints {
       PMT_ALWAYS_ASSERT(elmIDoffset == 0 || elmIDoffset == 1);
       return elmIDoffset;
     }
+    int getMaxAppID() {
+      PMT_ALWAYS_ASSERT(maxAppID != -1);
+      return maxAppID;
+    }
 
 //MUTATOR  
     template <MaterialPointSlice index> void fillData(double value);//use PS_LAMBDA fill up to 1

test/readMPAS.f90

@@ -8,13 +8,29 @@ subroutine assert(condition,message)
   endif
 end subroutine
 
+
+
 module readMPAS
     use :: polympo
     use iso_c_binding
     implicit none
     integer, parameter :: MPAS_RKIND = selected_real_kind(12)
 
 contains
+
+function epsilonDiff(a,b) result(isSame)
+  implicit none
+  real(kind=MPAS_RKIND) :: a,b,delta
+  parameter (delta=1.0e-8)
+  logical :: isSame
+  !delta=1.0e-8
+  if (abs(a-b) < delta) then
+    isSame = .true.
+  else
+    isSame = .false.
+  endif
+end function
+
 !---------------------------------------------------------------------------
 !> @brief get the MP positions array from a polympo array
 !> @param mpmesh(in/out) MPMesh object to fill, allocated by users
@@ -30,13 +46,16 @@ subroutine loadMPASMesh(mpMesh, filename)
     character (len=64) :: onSphere, stringYes = "YES"
     integer :: i
     integer :: maxEdges, vertexDegree, nCells, nVertices
+    integer, parameter :: nDims = 3
     integer :: numMPs
+    real(kind=MPAS_RKIND) :: ptOne = 0.100000000000000000
     real(kind=MPAS_RKIND) :: sphereRadius
     integer, dimension(:), pointer :: nEdgesOnCell
     real(kind=MPAS_RKIND), dimension(:), pointer :: xVertex, yVertex, zVertex
     integer, dimension(:,:), pointer :: verticesOnCell, cellsOnCell
     integer, dimension(:), pointer :: mpsPerElm, mp2Elm, isMPActive
-
+    real(kind=MPAS_RKIND), dimension(:,:), pointer :: mpPosition
+
     call readMPASMesh(trim(filename), maxEdges, vertexDegree, &
                               nCells, nVertices, nEdgesOnCell, &
                               onSphere, sphereRadius, &
@@ -94,6 +113,25 @@ subroutine loadMPASMesh(mpMesh, filename)
                       !i=numMPs leads to mp2Elm(numMPs=nCells+2)=numMPs-2=nCells
     end do
     call polympo_createMPs(mpMesh,nCells,numMPs,c_loc(mpsPerElm),c_loc(mp2Elm),c_loc(isMPActive))
+
+    !set mp positions
+    allocate(mpPosition(nDims,numMPs))
+    do i = 1,numMPs
+      mpPosition(1,i) = i+ptOne
+      mpPosition(2,i) = numMPs+i+ptOne
+      mpPosition(3,i) = (2*numMPs)+i+ptOne
+    end do
+
+    call polympo_setMPPositions(mpMesh,nDims,numMPs,c_loc(mpPosition))
+    mpPosition = 42
+    call polympo_getMPPositions(mpMesh,nDims,numMPs,c_loc(mpPosition))
+    do i = 1,numMPs
+      if(isMPActive(i) .eq. 1) then
+        call assert(epsilonDiff(mpPosition(1,i),i+ptOne), "x position of MP does not match")
+        call assert(epsilonDiff(mpPosition(2,i),numMPs+i+ptOne), "y position of MP does not match")
+        call assert(epsilonDiff(mpPosition(3,i),(2*numMPs)+i+ptOne), "z position of MP does not match")
+      endif
+    end do
 
     mp2Elm = -99 !override values and then use get function below
     call polympo_getMPCurElmID(mpMesh,numMPs,c_loc(mp2Elm))
@@ -106,6 +144,7 @@ subroutine loadMPASMesh(mpMesh, filename)
     end do
     !test end
 
+    deallocate(mpPosition)
     deallocate(mpsPerElm)
     deallocate(mp2Elm)
     deallocate(isMPActive)