Merge lp:~jonmcc/esys-particle/remote-force-periodicity into lp:~llaniewski/esys-particle/remote-force

[Jon McCullough]

Hopeful merge of Jon's and Lukasz's esys-particle codes for coupling with TCLB

- Introduction of full periodicity to esys in the remote force periodicity framework.
- should ensure it works when periodicity is not required. Anecdotally this works but I am skeptical.

=== modified file 'Geometry/CircularNeighbourTable.h'
--- Geometry/CircularNeighbourTable.h	2017-01-04 08:26:43 +0000
+++ Geometry/CircularNeighbourTable.h	2018-11-27 02:58:52 +0000
@@ -109,6 +109,10 @@
       ParticleSet     m_clonedParticleSet;
       int             m_circGridWidth;
       int             m_periodicDimIndex;
+      
+      //JM
+      int            m_periodicDimList[3];// = {-1,-1,-1};
+      int            m_circGridList[3];// = {0,0,0};
     };
   }
 }

=== modified file 'Geometry/CircularNeighbourTable.hpp'
--- Geometry/CircularNeighbourTable.hpp	2017-01-04 08:26:43 +0000
+++ Geometry/CircularNeighbourTable.hpp	2018-11-27 02:58:52 +0000
@@ -48,6 +48,21 @@
         circBorderWidth = this->getGridSpacing();
       }
       setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
+      
+      //JM
+      for (int i = 0; i < 3; i++) {
+      double checkcircBorderWidth =
+                   min(
+                     (bBox.getSizes()[i]),
+                     circBorderWidth
+                   );
+                   
+                   m_circGridList[i] = int(ceil(checkcircBorderWidth/gridSpacing));
+        //m_circGridList[i] = m_circGridWidth;
+      }
+      
+      //JM
+      cerr << "init1 circGridlist = " << m_circGridList[0] << "," << m_circGridList[1] << "," << m_circGridList[2] << endl;
     }
 
     template <class TmplParticle>
@@ -70,6 +85,14 @@
         circBorderWidth = this->getGridSpacing();
       }
       setCircularBorderWidth(circBorderWidth, this->getGridSpacing());
+      
+      //JM
+     for (int i = 0; i < 3; i++) {
+        m_circGridList[i] = m_circGridWidth;
+      }
+  
+      //JM
+      cerr << "init2 circGridlist = " << m_circGridList[0] << "," << m_circGridList[1] << "," << m_circGridList[2] << endl;
     }
 
     template <class TmplParticle>
@@ -86,20 +109,26 @@
       }
       int numPeriodic = 0;
       for (int i = 0; i < 3; i++) {
+         m_periodicDimList[i] = -1;
         if (m_periodicDimensions[i])
         {
-          m_periodicDimIndex = i;
+          m_periodicDimIndex = i; // JM need to select multiple indices here
           numPeriodic++;
+          
+          //JM
+          m_periodicDimList[i] = i;
         }
       }
-
+      cerr << "dimlist = " << m_periodicDimList[0] << "," << m_periodicDimList[1] << "," << m_periodicDimList[2] << endl;   
+      
       if (numPeriodic > 1) {
         std::stringstream msg;
         msg
           << "CircularNeighbourTable::CircularNeighbourTable -"
           << " only a single dimension may be periodic.";
-        throw std::runtime_error(msg.str());
-      }
+        //throw std::runtime_error(msg.str());
+      } 
+           
     }
 
     template <class TmplParticle>
@@ -133,13 +162,30 @@
     {
       if (havePeriodicDimensions())
       {
-        circBorderWidth =
+      //JM repeat for each periodic direction
+        /* circBorderWidth =
           min(
             (bBox.getSizes()[m_periodicDimIndex])/2.0,
             circBorderWidth
-          );
+          ); */ 
+       // end loop
+       
+          for (int i = 0; i < 3; i++) {
+             if (m_periodicDimList[i] != -1) {
+                circBorderWidth =
+                   min(
+                     (bBox.getSizes()[m_periodicDimList[i]])/2.0,
+                     circBorderWidth
+                   );
+                   
+                   m_circGridList[i] = int(ceil(circBorderWidth/gridSpacing)); //JM direct use of setCircularBorderWidth(circBorderWidth, gridSpacing);
+             }
+         }
+      } 
+      else { //JM
+      setCircularBorderWidth(circBorderWidth, gridSpacing); //JM now included in loop above when periodic in place.
       }
-      setCircularBorderWidth(circBorderWidth, gridSpacing);
+      //JM setCircularBorderWidth(circBorderWidth, gridSpacing);
       ParticleVector particles = getNonClonedParticles();
       clearClonedParticles();
       this->clearAndRecomputeGrid(bBox, gridSpacing);
@@ -151,6 +197,9 @@
       {
         this->insert(*it);
       }
+      
+      //JM
+      cerr << "Resize circGridlist = " << m_circGridList[0] << "," << m_circGridList[1] << "," << m_circGridList[2] << endl;
     }
 
     template <class TmplParticle>
@@ -189,33 +238,45 @@
         havePeriodicDimensions()
       )
       {
-        const int &dimIdx = m_periodicDimIndex;
-        if (
-          (posn[dimIdx] < this->getBBox().getMinPt()[dimIdx])
-          ||
-          (posn[dimIdx] > this->getBBox().getMaxPt()[dimIdx])
-        )
-        {
-          Vec3 moddedPosn = posn;
-          const double dimSize = this->getBBox().getSizes()[dimIdx];
-          moddedPosn[dimIdx] -= this->getBBox().getMinPt()[dimIdx];
-          moddedPosn[dimIdx] =
-            (moddedPosn[dimIdx] > 0.0)
-            ?
-            (
-              this->getBBox().getMinPt()[dimIdx]
-              +
-              moddedPosn[dimIdx] - (floor(moddedPosn[dimIdx]/dimSize)*dimSize)
-            )
-            :
-            (
-              this->getBBox().getMaxPt()[dimIdx]
-              -
-              (fabs(moddedPosn[dimIdx]) - (floor(fabs(moddedPosn[dimIdx])/dimSize)*dimSize))
-            );
-
-          return moddedPosn;
-        }
+         
+         Vec3 moddedPosn = posn;
+              
+         for (int i = 0; i < 3; i++) {
+             if (m_periodicDimList[i] != -1) {
+        
+                 //JM repeat this loop over each periodic direction
+                 //const int &dimIdx = m_periodicDimIndex;
+                 const int &dimIdx = m_periodicDimList[i];
+                                  
+                 if (
+                   (posn[dimIdx] < this->getBBox().getMinPt()[dimIdx])
+                   ||
+                   (posn[dimIdx] > this->getBBox().getMaxPt()[dimIdx])
+                 )
+                 {
+                   //Vec3 moddedPosn = posn;
+                   const double dimSize = this->getBBox().getSizes()[dimIdx];
+                   moddedPosn[dimIdx] -= this->getBBox().getMinPt()[dimIdx];
+                   moddedPosn[dimIdx] =
+                     (moddedPosn[dimIdx] > 0.0)
+                     ?
+                     (
+                       this->getBBox().getMinPt()[dimIdx]
+                       +
+                       moddedPosn[dimIdx] - (floor(moddedPosn[dimIdx]/dimSize)*dimSize)
+                     )
+                     :
+                     (
+                       this->getBBox().getMaxPt()[dimIdx]
+                       -
+                       (fabs(moddedPosn[dimIdx]) - (floor(fabs(moddedPosn[dimIdx])/dimSize)*dimSize))
+                     );
+                }
+            }
+         }
+         //JM end periodic repeat here
+             return moddedPosn;
+           //}
       }
       return posn;
     }
@@ -223,6 +284,9 @@
     template <class TmplParticle>
     void CircularNeighbourTable<TmplParticle>::insert(Particle *pParticle)
     {
+      //JM
+      //cerr << "insert circGridlist = " << m_circGridList[0] << "," << m_circGridList[1] << "," << m_circGridList[2] << endl;
+         
       pParticle->moveTo(getModdedPosn(pParticle->getPos()));
       const Vec3L minIdx = this->getVecIndex(pParticle->getPos() - pParticle->getRad());
       const Vec3L maxIdx = this->getVecIndex(pParticle->getPos() + pParticle->getRad());
@@ -230,6 +294,7 @@
       this->insertInTable(pParticle, minIdx, maxIdx);
       this->addInserted(pParticle);
 
+       //original
       if (havePeriodicDimensions())
       {
         for (int dimIdx = 0; dimIdx < 3; dimIdx++) {
@@ -247,6 +312,24 @@
           }
         }
       }
+        
+      /* //modified
+      if (havePeriodicDimensions())
+      {
+        Vec3 shift = Vec3::ZERO;
+        for (int dimIdx = 0; dimIdx < 3; dimIdx++) {
+          if (m_periodicDimensions[dimIdx]) {
+            if (minIdx[dimIdx] < (this->getMinVecIndex()[dimIdx] + m_circGridList[dimIdx])) {  //JM
+              shift[dimIdx] = this->getBBox().getSizes()[dimIdx];
+            }
+            if (maxIdx[dimIdx] > (this->getMaxVecIndex()[dimIdx] - m_circGridList[dimIdx])) { //JM
+              shift[dimIdx] = -1.0d*this->getBBox().getSizes()[dimIdx];
+            }
+          }
+        }
+        this->insertClone(pParticle, pParticle->getPos() + shift);
+      }
+      */
     }
 
     template <class TmplParticle>

=== modified file 'Geometry/FaultedBlock2d.cpp'
--- Geometry/FaultedBlock2d.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/FaultedBlock2d.cpp	2018-11-27 02:58:52 +0000
@@ -29,6 +29,10 @@
   \param rmax maximum particle radius
   \param pad thickness of the padding region at each diving edge
   \param circ_x circular boudary condition in x-direction
+  
+  JM
+  \param circ_y circular or open boundary conditions in y-direction
+  \param circ_z circular or open boundary conditions in z-direction
 */
 FaultedBlock2D::FaultedBlock2D(double xmin,double xmax,double ymin,double ymax,double rmin,double rmax,double pad,bool circ_x) :
   CRandomBlock2D(xmin,xmax,ymin,ymax,rmin,rmax,circ_x)

=== modified file 'Geometry/GougeBlock3D.cpp'
--- Geometry/GougeBlock3D.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/GougeBlock3D.cpp	2018-11-27 02:58:52 +0000
@@ -454,7 +454,10 @@
           m_prms.getPeriodicDimensions()[2] ? 1 : 0
         )*m_prms.getMaxRadius();
       if (m_prms.getBBox().getSizes().Z() >= 4*m_prms.getMaxRadius()) {
-        ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        //ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        
+        //JM For y/z circular bounds
+        ntableAdjust += Vec3(m_prms.getMaxRadius(), m_prms.getMaxRadius(), m_prms.getMaxRadius());
       }
       const BoundingBox nTableBBox(bBox.getMinPt(), bBox.getMaxPt() - ntableAdjust);
       m_nTablePtr =

=== modified file 'Geometry/GougeConfig.hpp'
--- Geometry/GougeConfig.hpp	2017-01-04 08:26:43 +0000
+++ Geometry/GougeConfig.hpp	2018-11-27 02:58:52 +0000
@@ -652,7 +652,10 @@
           m_prms.getPeriodicDimensions()[2] ? 1 : 0
         )*m_prms.getMaxRadius();
       if (m_prms.getBBox().getSizes().Z() >= 4*m_prms.getMaxRadius()) {
-        ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        //ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        
+        //JM For y/z circular bounds
+        ntableAdjust += Vec3(m_prms.getMaxRadius(), m_prms.getMaxRadius(), m_prms.getMaxRadius());
       }
       const BoundingBox nTableBBox(bBox.getMinPt(), bBox.getMaxPt() - ntableAdjust);
       m_nTablePtr =

=== modified file 'Geometry/SimpleNTable.cpp'
--- Geometry/SimpleNTable.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/SimpleNTable.cpp	2018-11-27 02:58:52 +0000
@@ -142,6 +142,23 @@
       }
     }
   }
+  
+  // JM Circular y bound addition
+  if(m_ycirc){
+    if (int((cbp.getPos().Y()-m_p0.Y())/m_dim)==1) {
+      cbp.moveTo(cbp.getPos()+m_yshift);
+      const vector<int> idy=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idy.begin();citer!=idy.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    } else if (int((cbp.getPos().Y()-m_p0.Y())/m_dim)==m_ysize-2) {
+      cbp.moveTo(cbp.getPos()-m_yshift);
+      const vector<int> idy=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idy.begin();citer!=idy.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    }
+  }
 }
 
 /*!

=== modified file 'Geometry/SimpleNTable3D.cpp'
--- Geometry/SimpleNTable3D.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/SimpleNTable3D.cpp	2018-11-27 02:58:52 +0000
@@ -152,6 +152,37 @@
       }
     }
   }
+  // JM Circular y/z bound addition
+  if(m_ycirc){
+    if (int((cbp.getPos().Y()-m_p0.Y())/m_dim)==1) {
+      cbp.moveTo(cbp.getPos()+m_yshift);
+      const vector<int> idy=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idy.begin();citer!=idy.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    } else if (int((cbp.getPos().Y()-m_p0.Y())/m_dim)==m_ysize-2) {
+      cbp.moveTo(cbp.getPos()-m_yshift);
+      const vector<int> idy=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idy.begin();citer!=idy.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    }
+  }
+  if(m_zcirc){
+    if (int((cbp.getPos().Z()-m_p0.Z())/m_dim)==1) {
+      cbp.moveTo(cbp.getPos()+m_zshift);
+      const vector<int> idz=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idz.begin();citer!=idz.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    } else if (int((cbp.getPos().Z()-m_p0.Z())/m_dim)==m_zsize-2) {
+      cbp.moveTo(cbp.getPos()-m_zshift);
+      const vector<int> idz=allidx(cbp.getPos());
+      for(vector<int>::const_iterator citer=idz.begin();citer!=idz.end();citer++){
+        m_data[*citer].push_back(cbp);
+      }
+    }
+  }
 }
 
 
@@ -187,7 +218,7 @@
   if (m_zcirc) {
     m_zsize+=2;
     m_p0-=Vec3(0.0,0.0,r);
-    m_yshift=Vec3(0.0,0.0,dim.Z());
+    m_zshift=Vec3(0.0,0.0,dim.Z());
   }
   //cout << "NT size : " << m_xsize << " , " << m_ysize << " , " << m_zsize << endl;
   m_data=new vector<SimpleParticle>[m_xsize*m_ysize*m_zsize];

=== modified file 'Geometry/SphAggGougeBlock.cpp'
--- Geometry/SphAggGougeBlock.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/SphAggGougeBlock.cpp	2018-11-27 02:58:52 +0000
@@ -162,7 +162,10 @@
       const BoundingBox bBox = m_prms.getBBox();
       Vec3 ntableAdjust = Vec3(0.0,0.0,0.0);
       if (m_prms.getBBox().getSizes().Z() >= 4*m_prms.getMaxRadius()) {
-        ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        //ntableAdjust += Vec3(m_prms.getMaxRadius(), 0, 0);
+        
+        //JM For y/z circular bounds
+        ntableAdjust += Vec3(m_prms.getMaxRadius(), m_prms.getMaxRadius(), m_prms.getMaxRadius());
       }
       const BoundingBox nTableBBox(bBox.getMinPt(), bBox.getMaxPt() - ntableAdjust);
       m_nTablePtr2 =

=== modified file 'Geometry/SplitBlock3D.cpp'
--- Geometry/SplitBlock3D.cpp	2017-01-04 08:26:43 +0000
+++ Geometry/SplitBlock3D.cpp	2018-11-27 02:58:52 +0000
@@ -27,6 +27,10 @@
   \param ysplit the position of the split plane 
   \param dir the direction of the split plane (2=y,3=z)
   \param circ_x circular or open boundary conditions in x-direction
+  
+  JM
+  \param circ_y circular or open boundary conditions in y-direction
+  \param circ_z circular or open boundary conditions in z-direction
 */
 CSplitBlock3D::CSplitBlock3D(double xmin,double xmax,double ymin,double ymax,double zmin,double zmax,double rmin,double rmax,double ysplit,int dir, bool circ_x,bool rough):CRandomBlock3D(xmin,xmax,ymin,ymax,zmin,zmax,rmin,rmax,1.05,circ_x)
 {

=== modified file 'Model/Particle.cpp'
--- Model/Particle.cpp	2017-01-04 08:26:43 +0000
+++ Model/Particle.cpp	2018-11-27 02:58:52 +0000
@@ -34,7 +34,8 @@
     m_mass(0.0),
     m_div_mass(0.0),
     flag(false),
-  m_is_dynamic(true)
+  m_is_dynamic(true),
+  m_is_acc(true)
 {
 }
 
@@ -49,7 +50,8 @@
     m_mass(data.getMass()),
     m_div_mass((data.getMass() != 0) ? 1.0/data.getMass() : 0.0),
     flag(false),
-    m_is_dynamic(true) //! SimpleParticleData has no is_dynamic info -> default to true 
+    m_is_dynamic(true), //! SimpleParticleData has no is_dynamic info -> default to true 
+    m_is_acc(true) //! SimpleParticleData has no is_dynamic info -> default to true 
 {
 }
 
@@ -63,7 +65,7 @@
   \param force currently applied force
   \param id particle id
 */
-CParticle::CParticle(double rad,double mass,const Vec3& pos,const Vec3& vel,const Vec3& force,int id,bool is_dyn) : CBasicParticle(pos,rad)
+CParticle::CParticle(double rad,double mass,const Vec3& pos,const Vec3& vel,const Vec3& force,int id,bool is_dyn,bool is_acc) : CBasicParticle(pos,rad)
 {
   m_oldpos=pos;
   m_initpos=pos;
@@ -79,9 +81,10 @@
   m_circular_shift = zeroVec3;
   flag=false;
   m_is_dynamic=is_dyn;
+  m_is_acc=is_acc;
 }
 
-CParticle::CParticle(double rad,double mass,const Vec3& pos,const Vec3& oldpos,const Vec3& initpos,const Vec3& vel,const Vec3& force,int id,bool is_dyn) : CBasicParticle(pos,rad)
+CParticle::CParticle(double rad,double mass,const Vec3& pos,const Vec3& oldpos,const Vec3& initpos,const Vec3& vel,const Vec3& force,int id,bool is_dyn,bool is_acc) : CBasicParticle(pos,rad)
 {
   m_oldpos=oldpos;
   m_initpos=initpos;
@@ -97,6 +100,7 @@
   m_circular_shift = zeroVec3;
   flag=false;
   m_is_dynamic=is_dyn;
+  m_is_acc=is_acc;
 }
 
 /*!
@@ -211,6 +215,7 @@
   append(p.m_rad);
   append(p.m_mass);
   append(p.m_is_dynamic);
+  append(p.m_is_acc);
 }
 
 /*!
@@ -240,6 +245,7 @@
     p.m_div_mass=0.0;
   }
   p.m_is_dynamic=pop_bool();
+  p.m_is_acc=pop_bool();
 }
 
 /*!
@@ -253,6 +259,9 @@
     if (getDo2dCalculations()) {
       m_force = Vec3(m_force.X(), m_force.Y(), 0);
     }
+    if (!m_is_acc) {
+      m_force = Vec3(0.0,0.0,0.0);
+    }
 #if 1
     // 1st order scheme,  13 Flops
     m_vel+=m_force*m_div_mass*dt; // 7 Flops 
@@ -303,7 +312,8 @@
       m_initpos - m_circular_shift,
       m_oldpos  - m_circular_shift,
       m_vel,
-      m_is_dynamic
+      m_is_dynamic,
+      m_is_acc
     );
 }
 
@@ -319,6 +329,7 @@
   m_oldpos  = E.m_oldPos  + m_circular_shift;
   m_vel     = E.m_vel;
   m_is_dynamic = E.m_is_dynamic;
+  m_is_acc = E.m_is_acc;
 }
 
 /*!
@@ -355,6 +366,7 @@
   append(p.m_vel.Y());
   append(p.m_vel.Z());
   append(p.m_is_dynamic);
+  append(p.m_is_acc);
 }
 
 /*!
@@ -372,6 +384,7 @@
   p.m_oldPos  = Vec3(db[6],  db[7],  db[8]);
   p.m_vel     = Vec3(db[9], db[10], db[11]);
   p.m_is_dynamic=pop_bool();
+  p.m_is_acc=pop_bool();
 }
 
 
@@ -448,7 +461,8 @@
     << getVel()                        << delim // 13,14,15
     << getForce()                      << delim // 16,17,18
     << m_circular_shift << delim                // 19,20,21
-    << m_is_dynamic;
+    << m_is_dynamic << delim
+    << m_is_acc;
 }
 
 /*!
@@ -467,6 +481,7 @@
   int id;
   int tag;
   bool is_dyn;
+  bool is_acc;
 
   iStream
     >> pos
@@ -479,7 +494,8 @@
     >> vel
     >> force
     >> circular_shift
-    >> is_dyn;
+    >> is_dyn
+    >> is_acc;
 
   CParticle 
     particle
@@ -491,7 +507,7 @@
       initpos + circular_shift,
       vel,
       force,
-      id, is_dyn
+      id, is_dyn, is_acc
     );
   particle.setTag(tag);
   

=== modified file 'Model/Particle.h'
--- Model/Particle.h	2017-01-04 08:26:43 +0000
+++ Model/Particle.h	2018-11-27 02:58:52 +0000
@@ -60,15 +60,17 @@
       m_vel()
     {
       m_is_dynamic=true;
+      m_is_acc=true;
     }
 
-      exchangeType(const Vec3 &pos, const Vec3 &initPos, const Vec3 &oldPos, const Vec3 &vel,bool is_dyn)
+      exchangeType(const Vec3 &pos, const Vec3 &initPos, const Vec3 &oldPos, const Vec3 &vel,bool is_dyn,bool is_acc)
       : m_pos(pos),
         m_initPos(initPos),
         m_oldPos(oldPos),
         m_vel(vel)
     {
       m_is_dynamic=is_dyn;
+      m_is_acc=is_acc;
     }
     
     Vec3 m_pos;
@@ -76,6 +78,7 @@
     Vec3 m_oldPos;    
     Vec3 m_vel;
     bool m_is_dynamic;
+    bool m_is_acc;
   };
   
   typedef double (CParticle::* ScalarFieldFunction)() const; 
@@ -92,13 +95,14 @@
 
   bool flag;
   bool m_is_dynamic; 
+  bool m_is_acc;
 
   void setForce(const Vec3 &force) {m_force = force;}
 
  public:
   CParticle();
-  CParticle(double,double,const Vec3&,const Vec3&,const Vec3&,int,bool);
-  CParticle(double,double,const Vec3&,const Vec3&,const Vec3&,const Vec3&,const Vec3&,int,bool); // including oldpos
+  CParticle(double,double,const Vec3&,const Vec3&,const Vec3&,int,bool,bool);
+  CParticle(double,double,const Vec3&,const Vec3&,const Vec3&,const Vec3&,const Vec3&,int,bool,bool); // including oldpos
   CParticle(const esys::lsm::SimpleParticleData &particleData);
   virtual ~CParticle(){};
 
@@ -136,6 +140,7 @@
   virtual void setNonDynamic() {m_is_dynamic=false;};
   virtual void setNonDynamicLinear() {m_is_dynamic=false;};
   virtual void setNonDynamicRot(){}; // do nothing
+  virtual void setNonDynamicAcc(){m_is_dynamic=true; m_is_acc=false;}; // JM
 
   void setFlag(bool b=true){flag=b;};
   bool isFlagged() const {return flag;};

=== modified file 'Model/RotParticle.cpp'
--- Model/RotParticle.cpp	2017-01-04 08:26:43 +0000
+++ Model/RotParticle.cpp	2018-11-27 02:58:52 +0000
@@ -66,8 +66,9 @@
   const Vec3& force,
   int id,
   bool is_dyn,
-  bool is_rot
-  ) : CParticle(rad, mass, pos, vel, force, id, is_dyn)
+  bool is_rot,
+  bool is_acc
+  ) : CParticle(rad, mass, pos, vel, force, id, is_dyn, is_acc)
 {
   m_quat = Quaternion(1.0,Vec3(0.0,0.0,0.0));
   m_initquat = m_quat;
@@ -109,8 +110,9 @@
   const Vec3& moment,
   const Vec3& angvel,
   bool is_dyn,
-  bool is_rot
-  ) :CParticle(rad, mass,pos,vel,force,id,is_dyn)
+  bool is_rot,
+  bool is_acc
+  ) :CParticle(rad, mass,pos,vel,force,id,is_dyn,is_acc)
 {
   m_circular_shift=Vec3(0.0,0.0,0.0);
   flag=false;
@@ -139,8 +141,9 @@
   const Vec3& moment,
   const Vec3& angvel,
   bool is_dyn,
-  bool is_rot
-  ) :CParticle(rad, mass,pos,oldpos, initpos,vel,force,id,is_dyn)
+  bool is_rot,
+  bool is_acc
+  ) :CParticle(rad, mass,pos,oldpos, initpos,vel,force,id,is_dyn,is_acc)
 {
   m_circular_shift=Vec3(0.0,0.0,0.0);
   flag=false;
@@ -228,6 +231,7 @@
   
   append(p.m_is_dynamic);
   append(p.m_is_rot);
+  append(p.m_is_acc);
 }
 
 
@@ -265,6 +269,7 @@
   p.m_global_id=pop_int();
   p.m_is_dynamic=pop_bool();
   p.m_is_rot=pop_bool();
+  p.m_is_acc=pop_bool();
 }
 
 /*!
@@ -339,7 +344,8 @@
       m_angVel,
       m_quat,
       m_is_dynamic,
-      m_is_rot
+      m_is_rot,
+      m_is_acc
     );
 }
 
@@ -357,6 +363,7 @@
   m_quat       = E.m_quat;
   m_is_dynamic = E.m_is_dynamic;
   m_is_rot = E.m_is_rot;
+  m_is_acc = E.m_is_acc;
 }
 
 template<>
@@ -380,6 +387,7 @@
   append(p.m_quat.return_vec().Z());
   append(p.m_is_dynamic);
   append(p.m_is_rot);
+  append(p.m_is_acc);
 }
 
 /*!
@@ -399,6 +407,7 @@
   p.m_quat       = Quaternion(db[12],Vec3(db[13],db[14],db[15]));
   p.m_is_dynamic = pop_bool();
   p.m_is_rot     = pop_bool();
+  p.m_is_acc     = pop_bool();
 }
 
 CRotParticle::ScalarFieldFunction CRotParticle::getScalarFieldFunction(const string& name)

=== modified file 'Model/RotParticle.h'
--- Model/RotParticle.h	2017-01-04 08:26:43 +0000
+++ Model/RotParticle.h	2018-11-27 02:58:52 +0000
@@ -66,6 +66,7 @@
     {
       m_is_dynamic=true;
       m_is_rot=true;
+      m_is_acc=true;
     }
 
     exchangeType(
@@ -75,7 +76,8 @@
       const Vec3 &currAngVel,
       const Quaternion &quat,
       const bool is_dyn,
-      const bool is_rot
+      const bool is_rot,
+      const bool is_acc
     )
       : m_pos(pos),
         m_initPos(initPos),
@@ -85,6 +87,7 @@
     {
       m_is_dynamic=is_dyn;
       m_is_rot=is_rot;
+      m_is_acc=is_acc;
     }
   public:
     Vec3       m_pos;
@@ -94,6 +97,7 @@
     Quaternion m_quat;
     bool m_is_dynamic;
     bool m_is_rot;
+    bool m_is_acc;
 
     friend class TML_PackedMessageInterface;    
   };
@@ -124,7 +128,8 @@
     const Vec3& force,
     int         id,
     bool is_dyn,
-    bool is_rot
+    bool is_rot,
+    bool is_acc
   );
   CRotParticle(
     double      rad,
@@ -138,7 +143,8 @@
     const Vec3& moment,
     const Vec3& angvel, 
     bool is_dyn, 
-    bool is_rot
+    bool is_rot,
+    bool is_acc
   );
   CRotParticle(
     double            rad,
@@ -155,7 +161,8 @@
     const Vec3&       moment,
     const Vec3&       angvel,  
     bool is_dyn,
-    bool is_rot
+    bool is_rot,
+    bool is_acc
   );
   
   virtual ~CRotParticle(){};
@@ -200,6 +207,7 @@
   virtual void setNonDynamic() {m_is_dynamic=false;m_is_rot=false;};
   virtual void setNonDynamicLinear() {m_is_dynamic=false;};
   virtual void setNonDynamicRot() {m_is_rot=false;};
+  virtual void setNonDynamicAcc() {m_is_dynamic=true; m_is_acc=false;};
 
   inline Quaternion getQuatFromRotVec(const Vec3 &vec) const
   {

=== modified file 'Model/RotParticleVi.cpp'
--- Model/RotParticleVi.cpp	2017-01-04 08:26:43 +0000
+++ Model/RotParticleVi.cpp	2018-11-27 02:58:52 +0000
@@ -42,8 +42,9 @@
   const Vec3& vel,
   const Vec3& force,
   int id,
-  bool is_dyn
-) : CParticle(rad, mass, pos, vel, force, id, is_dyn)
+  bool is_dyn,
+  bool is_acc
+) : CParticle(rad, mass, pos, vel, force, id, is_dyn,is_acc)
 {
   m_quat = Quaternion(1.0,Vec3(0.0,0.0,0.0));
   m_initquat = m_quat;
@@ -53,6 +54,7 @@
   m_angVel_t = Vec3(0.0,0.0,0.0) ;
   m_moment     = Vec3(0.0,0.0,0.0);
   m_is_dynamic=is_dyn;
+  m_is_acc=is_acc;
 }
 
 CRotParticleVi::CRotParticleVi(const CParticle &particle) : CParticle(particle)
@@ -91,7 +93,7 @@
   const Vec3& moment,
   const Vec3& angvel,
   const Vec3& angvel_t
-) :CParticle(rad, mass,pos,vel,force,id,true)
+) :CParticle(rad, mass,pos,vel,force,id,true,true)
 {
   m_circular_shift=Vec3(0.0,0.0,0.0);
   flag=false;
@@ -120,7 +122,7 @@
   const Vec3& moment,
   const Vec3& angvel,
   const Vec3& angvel_t
-) :CParticle(rad, mass,pos,oldpos, initpos,vel,force,id,true)
+) :CParticle(rad, mass,pos,oldpos, initpos,vel,force,id,true,true)
 {
   m_circular_shift=Vec3(0.0,0.0,0.0);
   flag=false;

=== modified file 'Model/RotParticleVi.h'
--- Model/RotParticleVi.h	2017-01-04 08:26:43 +0000
+++ Model/RotParticleVi.h	2018-11-27 02:58:52 +0000
@@ -105,6 +105,7 @@
   double      m_inertRot;
   double      m_div_inertRot;
   bool        m_is_dynamic; 
+  bool        m_is_acc; 
 
  public:
   CRotParticleVi();
@@ -116,7 +117,8 @@
     const Vec3& vel,
     const Vec3& force,
     int         id,
-    bool        is_dyn
+    bool        is_dyn,
+    bool	is_acc
   );
   CRotParticleVi(
     double      rad,

=== modified file 'Model/RotThermParticle.cpp'
--- Model/RotThermParticle.cpp	2017-01-04 08:26:43 +0000
+++ Model/RotThermParticle.cpp	2018-11-27 02:58:52 +0000
@@ -52,9 +52,10 @@
   const Vec3& vel,
   const Vec3& force,
   int id,
-  bool is_dyn
+  bool is_dyn,
+  bool is_acc
 ) :
-    CRotParticleVi(rad, mass, pos, vel, force, id, is_dyn),
+    CRotParticleVi(rad, mass, pos, vel, force, id, is_dyn, is_acc),
     CThermParticle(rad)
 {
 }

=== modified file 'Model/RotThermParticle.h'
--- Model/RotThermParticle.h	2017-01-04 08:26:43 +0000
+++ Model/RotThermParticle.h	2018-11-27 02:58:52 +0000
@@ -136,7 +136,8 @@
     const Vec3& vel,
     const Vec3& force,
     int         id,
-    bool        is_dyn
+    bool        is_dyn,
+    bool 	is_acc
   );
 
   CRotThermParticle(

=== modified file 'Parallel/ASubLattice.h'
--- Parallel/ASubLattice.h	2017-09-07 08:07:19 +0000
+++ Parallel/ASubLattice.h	2018-11-27 02:58:52 +0000
@@ -131,6 +131,7 @@
   virtual void setParticleAngularVelocity(){};
   virtual void setParticleNonDynamic()=0;
   virtual void setParticleNonRot()=0;
+  virtual void setParticleNonAcc()=0;
   virtual void tagParticleNearestTo()=0;
   virtual void moveSingleNode()=0;
   virtual void moveTaggedNodes()=0;

=== modified file 'Parallel/LatticeMaster.cpp'
--- Parallel/LatticeMaster.cpp	2017-12-11 00:23:07 +0000
+++ Parallel/LatticeMaster.cpp	2018-11-27 02:58:52 +0000
@@ -1091,6 +1091,22 @@
 }
 
 /*!
+  Set all particles with a given tag to be linear non-dynamic, i.e. switch off velocity updates JM version
+
+  \param tag the tag
+*/
+void CLatticeMaster::setParticleNonAcc(int tag)
+{
+  CMPILCmdBuffer cmd_buffer(m_global_comm,m_global_rank);
+  CVarMPIBuffer buffer(m_global_comm);
+  CMPIBarrier barrier(m_global_comm);
+
+  cmd_buffer.broadcast(CMD_PSETNA);
+  buffer.append(tag);
+  buffer.broadcast(m_global_rank);
+  barrier.wait("setParticleNonAcc");
+}
+/*!
   Set all particles with a given tag to be linear non-dynamic, i.e. switch off velocity updates
 
   \param tag the tag

=== modified file 'Parallel/LatticeMaster.h'
--- Parallel/LatticeMaster.h	2017-12-11 00:23:07 +0000
+++ Parallel/LatticeMaster.h	2018-11-27 02:58:52 +0000
@@ -303,6 +303,7 @@
     void setParticleNonDynamic(int);
     void setParticleNonRot(int);
     void setParticleNonTrans(int);
+    void setParticleNonAcc(int);
     void setParticleVel(int,const Vec3&);
     void setParticleAngVel(int,const Vec3&);
     void setParticleDensity(int tag,int mask,double rho);

=== modified file 'Parallel/SubLattice.h'
--- Parallel/SubLattice.h	2017-09-07 08:07:19 +0000
+++ Parallel/SubLattice.h	2018-11-27 02:58:52 +0000
@@ -226,6 +226,7 @@
   virtual void setParticleNonDynamic();
   virtual void setParticleNonRot();
   virtual void setParticleNonTrans();
+  virtual void setParticleNonAcc();
   virtual void setParticleDensity();
   virtual void setTaggedParticleVel();
   virtual void tagParticleNearestTo();

=== modified file 'Parallel/SubLattice.hpp'
--- Parallel/SubLattice.hpp	2017-12-11 00:23:07 +0000
+++ Parallel/SubLattice.hpp	2018-11-27 02:58:52 +0000
@@ -293,9 +293,13 @@
     ysize=max.Y()-min.Y();
     if(fabs(ysize/m_nrange-lrint(ysize/m_nrange))>1e-6)
     {
-      console.Critical() << "circular y-dimension doesn't fit range !\n";
+      //JM as per x-dim
+      //console.Critical() << "circular y-dimension doesn't fit range !\n";
+      console.Info() << "Circular y-size incompatible with range, adjusting...\n";
+      m_nrange = ysize/floor(ysize/m_nrange);
+      console.Info() << "New range = " << m_nrange << "\n";
     }
-    ysize+=2.0*m_nrange; // padding on the circular ends
+    //ysize+=2.0*m_nrange; // padding on the circular ends
   }
   else
   {
@@ -308,9 +312,13 @@
     zsize=max.Z()-min.Z();
     if(fabs(zsize/m_nrange-lrint(zsize/m_nrange))>1e-6)
     {
-      console.Critical() << "circular z-dimension doesn't fit range !\n";
+      //JM as per x-dim
+      //console.Critical() << "circular z-dimension doesn't fit range !\n";
+      console.Info() << "Circular z-size incompatible with range, adjusting...\n";
+      m_nrange = zsize/floor(zsize/m_nrange);
+      console.Info() << "New range = " << m_nrange << "\n";
     }
-    zsize+=2.0*m_nrange; // padding on the circular ends
+    //zsize+=2.0*m_nrange; // padding on the circular ends
   }
   else
   {
@@ -2537,6 +2545,22 @@
 }
 
 /*!
+  Make tagged particles non-accelerating, i.e. don't change velocity based on forces JM
+  Parameters (tag) are received from master. 
+*/
+template <class T>
+void TSubLattice<T>::setParticleNonAcc()
+{
+  console.Debug() << "TSubLattice<T>::setParticleNonAcc()\n";
+  CVarMPIBuffer buffer(m_comm);
+  
+  buffer.receiveBroadcast(0); // get data from master
+  int tag=buffer.pop_int();
+  m_ppa->forParticleTag(tag,(void (T::*)())(&T::setNonDynamicAcc));
+  console.Debug() << "end TSubLattice<T>::setParticleNonAcc()\n";
+}
+
+/*!
   Make tagged particles lieanr non-dynamic, i.e. don't update linear velocity 
   Parameters (tag) are received from master. 
 */

=== modified file 'Parallel/SubLatticeControler.cpp'
--- Parallel/SubLatticeControler.cpp	2017-12-11 00:23:07 +0000
+++ Parallel/SubLatticeControler.cpp	2018-11-27 02:58:52 +0000
@@ -551,6 +551,7 @@
     case CMD_PMOVETAGGEDBY : m_lattice->moveTaggedParticlesBy(); break;
     case CMD_PSETND : m_lattice->setParticleNonDynamic(); break;
     case CMD_PSETNR : m_lattice->setParticleNonRot(); break;
+    case CMD_PSETNA : m_lattice->setParticleNonAcc(); break;
     case CMD_PVEL : m_lattice->setParticleVelocity(); break;
     case CMD_PANGVEL : m_lattice->setParticleAngularVelocity(); break;
     case CMD_PDENS : m_lattice->setParticleDensity(); break;     

=== modified file 'Parallel/sublattice_cmd.h'
--- Parallel/sublattice_cmd.h	2017-09-07 08:07:19 +0000
+++ Parallel/sublattice_cmd.h	2018-11-27 02:58:52 +0000
@@ -101,7 +101,7 @@
 const int CMD_SETCONSOLEFNAME=1113;
 const int CMD_SETCONSOLEBUFF=1114;
 const int CMD_SETINTERACTIONPARAMS=1115;
-
+const int CMD_PSETNA = 1116;
 
 
 // ---------------------------------

=== modified file 'Python/esys/lsm/LsmMpiPy.cpp'
--- Python/esys/lsm/LsmMpiPy.cpp	2017-12-11 00:49:38 +0000
+++ Python/esys/lsm/LsmMpiPy.cpp	2018-11-27 02:58:52 +0000
@@ -1134,6 +1134,12 @@
     {
       getLatticeMaster().setParticleNonTrans(id);
     }
+    
+    //JM non-accelerative flag for fluid coupling
+    void LsmMpiPy::setParticleNonAcc(int id)
+    {
+      getLatticeMaster().setParticleNonAcc(id);
+    }
 
 	// --- modifying interaction parameters ---
 	void LsmMpiPy::setInteractionParameter(const std::string& igname,const std::string& pname,double val)
@@ -2879,6 +2885,17 @@
         "@type tag: int\n"
         "@kwarg tag: the tag of the particle."
       )
+	//JM
+       .def(
+        "setParticleNonAccelerating",
+        &LsmMpiPy::setParticleNonAcc,
+        (arg("tag")),
+        "Set the particle with tag C{tag} to be non-accelerating,"
+        " this is to allow for fluid coupled particles to mainain\n"
+        " a constant velocity e.g. for shear platen tests"
+        "@type tag: int\n"
+        "@kwarg tag: the tag of the particle."
+      )
       .def(
         "addPreTimeStepRunnable",
         &LsmMpiPy::addPreTimeStepRunnable,

=== modified file 'Python/esys/lsm/LsmMpiPy.h'
--- Python/esys/lsm/LsmMpiPy.h	2017-12-11 00:23:07 +0000
+++ Python/esys/lsm/LsmMpiPy.h	2018-11-27 02:58:52 +0000
@@ -205,6 +205,7 @@
       void setParticleNonDynamic(int);
       void setParticleNonRot(int);
       void setParticleNonTrans(int);
+      void setParticleNonAcc(int); //JM non-accelerating to allow constant velocity in fluid
       
       // --- modifying interaction parameters ---
       void setInteractionParameter(const std::string&,const std::string&,double);

=== modified file 'Python/esys/lsm/ParticlePy.cpp'
--- Python/esys/lsm/ParticlePy.cpp	2017-12-04 00:00:38 +0000
+++ Python/esys/lsm/ParticlePy.cpp	2018-11-27 02:58:52 +0000
@@ -23,7 +23,7 @@
     }
 
     ParticlePy::ParticlePy(int id, const Vec3Py &posn, double radius, double mass)
-      : CParticle(radius, mass, posn, Vec3(), Vec3(), id, true)
+      : CParticle(radius, mass, posn, Vec3(), Vec3(), id, true,true)
     {
     }
 

=== modified file 'Python/esys/lsm/RotParticlePy.cpp'
--- Python/esys/lsm/RotParticlePy.cpp	2017-01-04 08:26:43 +0000
+++ Python/esys/lsm/RotParticlePy.cpp	2018-11-27 02:58:52 +0000
@@ -34,7 +34,7 @@
 
     // dynamic & rotational default to true
     RotParticlePy::RotParticlePy(int id, const Vec3Py &posn, double radius, double mass)
-      : CRotParticle(radius, mass, posn, Vec3(), Vec3(), id,true,true)
+      : CRotParticle(radius, mass, posn, Vec3(), Vec3(), id,true,true,true)
     {
     }
 

=== modified file 'Python/esys/lsm/RotParticleViPy.cpp'
--- Python/esys/lsm/RotParticleViPy.cpp	2017-01-04 08:26:43 +0000
+++ Python/esys/lsm/RotParticleViPy.cpp	2018-11-27 02:58:52 +0000
@@ -24,7 +24,7 @@
     }
 
     RotParticleViPy::RotParticleViPy(int id, const Vec3Py &posn, double radius, double mass)
-      : CRotParticleVi(radius, mass, posn, Vec3(), Vec3(), id, true)
+      : CRotParticleVi(radius, mass, posn, Vec3(), Vec3(), id, true,true)
     {
     }
     

=== modified file 'Python/esys/lsm/RotThermalParticlePy.cpp'
--- Python/esys/lsm/RotThermalParticlePy.cpp	2017-01-04 08:26:43 +0000
+++ Python/esys/lsm/RotThermalParticlePy.cpp	2018-11-27 02:58:52 +0000
@@ -24,7 +24,7 @@
     }
 
     RotThermalParticlePy::RotThermalParticlePy(int id, const Vec3Py &posn, double radius, double mass)
-      : CRotThermParticle(radius, mass, posn, Vec3(), Vec3(), id, true)
+      : CRotThermParticle(radius, mass, posn, Vec3(), Vec3(), id, true, true)
     {
     }
     

=== modified file 'Tools/dump2vtk/frame_vtk.cpp'
--- Tools/dump2vtk/frame_vtk.cpp	2017-09-07 08:07:19 +0000
+++ Tools/dump2vtk/frame_vtk.cpp	2018-11-27 02:58:52 +0000
@@ -345,6 +345,16 @@
     vtkfile << (iter->second).angvel << endl;;
   } 
   vtkfile << "</DataArray>\n";
+
+  // force - JM addition
+  vtkfile << "<DataArray type=\"Float64\" Name=\"force\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+  for(map<int,r_part>::iterator iter=datamap.begin();
+      iter!=datamap.end();
+      iter++){
+    vtkfile << (iter->second).force << endl;;
+  } 
+  vtkfile << "</DataArray>\n";
+
   // displacement
   vtkfile << "<DataArray type=\"Float64\" Name=\"displacement\" NumberOfComponents=\"3\" format=\"ascii\">\n";
   for(map<int,r_part>::iterator iter=datamap.begin();
@@ -440,6 +450,16 @@
     vtkfile << (iter->second).vel << endl;;
   } 
   vtkfile << "</DataArray>\n";
+
+  // force - JM addition
+  vtkfile << "<DataArray type=\"Float64\" Name=\"force\" NumberOfComponents=\"3\" format=\"ascii\">\n";
+  for(map<int,nr_part>::iterator iter=datamap.begin();
+      iter!=datamap.end();
+      iter++){
+    vtkfile << (iter->second).force << endl;;
+  } 
+  vtkfile << "</DataArray>\n";
+
   // displacement
   vtkfile << "<DataArray type=\"Float64\" Name=\"displacement\" NumberOfComponents=\"3\" format=\"ascii\">\n";
   for(map<int,nr_part>::iterator iter=datamap.begin();

=== modified file 'ppa/src/pp_array.h'
--- ppa/src/pp_array.h	2017-12-11 00:23:07 +0000
+++ ppa/src/pp_array.h	2018-11-27 02:58:52 +0000
@@ -84,7 +84,7 @@
   NeighborTable<T>* m_nt;
   Vec3 m_minpos,m_maxpos; //!< local minimum and maximum positions
   double m_xshift,m_yshift,m_zshift; //!< circular shift values
-  bool m_circ_edge_x_up,m_circ_edge_x_down; //!< circular edge flags
+  bool m_circ_edge_x_up,m_circ_edge_x_down,m_circ_edge_y_up,m_circ_edge_y_down,m_circ_edge_z_up,m_circ_edge_z_down; //!< circular edge flags JM addition for y/z
   static const int m_exchg_tag;
   
   // helper fnc

=== modified file 'ppa/src/pp_array.hpp'
--- ppa/src/pp_array.hpp	2017-01-04 08:26:43 +0000
+++ ppa/src/pp_array.hpp	2018-11-27 02:58:52 +0000
@@ -64,6 +64,12 @@
   //- set x-edges to non-circular
   m_circ_edge_x_down=false;
   m_circ_edge_x_up=false;
+  //JM - set y-edges to non-circular
+  m_circ_edge_y_down=false;
+  m_circ_edge_y_up=false;
+  //JM - set z-edges to non-circular
+  m_circ_edge_z_down=false;
+  m_circ_edge_z_up=false;
   //- get own process coords
   vector<int> pcoords=m_comm.get_coords();
   //- initialize ntable
@@ -159,6 +165,12 @@
     // setup circular edges
     m_circ_edge_x_down=(circ[0] && (pcoords[0]==0)) ? true : false ;
     m_circ_edge_x_up=(circ[0] && (pcoords[0]==m_comm.get_dim(0)-1)) ? true : false ;
+    //JM setup circular edges
+    m_circ_edge_y_down=(circ[1] && (pcoords[1]==0)) ? true : false ;
+    m_circ_edge_y_up=(circ[1] && (pcoords[1]==m_comm.get_dim(1)-1)) ? true : false ;
+    //JM setup circular edges
+    m_circ_edge_z_down=(circ[2] && (pcoords[2]==0)) ? true : false ;
+    m_circ_edge_z_up=(circ[2] && (pcoords[2]==m_comm.get_dim(2)-1)) ? true : false ;
   }
   // init time stamp
   m_timestamp=0; 
@@ -325,13 +337,46 @@
     // define tranfer slabs
     up_slab=m_nt->xz_slab(m_nt->ysize()-2);
     down_slab=m_nt->xz_slab(1);
-    if(!circ_edge){ // normal boundary
-      // upstream
+    
+    // upstream
+    if(m_circ_edge_y_up){ // circ. bdry
+      // cout << "circular shift, node " << m_comm.rank() << " y-dim, up slab size " << up_slab.size() << endl;
+      // copy particles from transfer slab into buffer
+      vector<T> buffer;
+      for(typename NTSlab<T>::iterator iter=up_slab.begin();
+	  iter!=up_slab.end();
+	  iter++){
+	buffer.push_back(*iter);
+      }
+      // shift particles in buffer by circular shift value 
+      for(typename vector<T>::iterator iter=buffer.begin();
+	  iter!=buffer.end();
+	  iter++){
+	iter->setCircular(Vec3(0.0,-1.0*m_yshift,0.0));
+      }
+      m_comm.shift_cont_packed(buffer,*m_nt,1,1,m_exchg_tag);
+    } else {
       m_comm.shift_cont_packed(up_slab,*m_nt,1,1,m_exchg_tag);
-      // downstream
+    }
+    // downstream
+    if(m_circ_edge_y_down){// circ. bdry
+      // cout << "circular shift , node " << m_comm.rank() << " y-dim, down slab size " << down_slab.size() << endl;
+      // copy particles from transfer slab into buffer
+      vector<T> buffer;
+      for(typename NTSlab<T>::iterator iter=down_slab.begin();
+	  iter!=down_slab.end();
+	  iter++){
+	buffer.push_back(*iter);
+      }
+      // shift particles in buffer by circular shift value 
+      for(typename vector<T>::iterator iter=buffer.begin();
+	  iter!=buffer.end();
+	  iter++){
+	iter->setCircular(Vec3(0.0,m_yshift,0.0));
+      }
+      m_comm.shift_cont_packed(buffer,*m_nt,1,-1,m_exchg_tag);
+    } else {
       m_comm.shift_cont_packed(down_slab,*m_nt,1,-1,m_exchg_tag);
-    } else { // circular edge -> buffer & shift received particles
-         cout << "circular y shift not implemented" << endl;
     }
   }
   // z-dir
@@ -346,13 +391,46 @@
     // define tranfer slabs
     up_slab=m_nt->xy_slab(m_nt->zsize()-2);
     down_slab=m_nt->xy_slab(1);
-    if(!circ_edge){ // normal boundary
-      // upstream
+    
+    // upstream
+    if(m_circ_edge_z_up){ // circ. bdry
+      // cout << "circular shift, node " << m_comm.rank() << " z-dim, up slab size " << up_slab.size() << endl;
+      // copy particles from transfer slab into buffer
+      vector<T> buffer;
+      for(typename NTSlab<T>::iterator iter=up_slab.begin();
+	  iter!=up_slab.end();
+	  iter++){
+	buffer.push_back(*iter);
+      }
+      // shift particles in buffer by circular shift value 
+      for(typename vector<T>::iterator iter=buffer.begin();
+	  iter!=buffer.end();
+	  iter++){
+	iter->setCircular(Vec3(0.0,0.0,-1.0*m_zshift));
+      }
+      m_comm.shift_cont_packed(buffer,*m_nt,2,1,m_exchg_tag);
+    } else {
       m_comm.shift_cont_packed(up_slab,*m_nt,2,1,m_exchg_tag);
-      // downstream
+    }
+    // downstream
+    if(m_circ_edge_z_down){// circ. bdry
+      // cout << "circular shift , node " << m_comm.rank() << " z-dim, down slab size " << down_slab.size() << endl;
+      // copy particles from transfer slab into buffer
+      vector<T> buffer;
+      for(typename NTSlab<T>::iterator iter=down_slab.begin();
+	  iter!=down_slab.end();
+	  iter++){
+	buffer.push_back(*iter);
+      }
+      // shift particles in buffer by circular shift value 
+      for(typename vector<T>::iterator iter=buffer.begin();
+	  iter!=buffer.end();
+	  iter++){
+	iter->setCircular(Vec3(0.0,0.0,m_zshift));
+      }
+      m_comm.shift_cont_packed(buffer,*m_nt,2,-1,m_exchg_tag);
+    } else {
       m_comm.shift_cont_packed(down_slab,*m_nt,2,-1,m_exchg_tag);
-    } else { // circular edge -> buffer & shift received particles
-         cout << "circular x shift not implemented" << endl;
     }
   }
   // update timestamp
@@ -424,12 +502,14 @@
     send_data.push_back(((*iter).*rdf)());
 //     cout << "put exchange values from particle " << iter->getID() << " into buffer" << endl; 
   }
+  
   // exchange
   m_comm.shift_cont_packed(send_data,recv_data,dir,dist,m_exchg_tag);
-
+  
   // apply received data
   // check if sizes are correct
   if(recv_data.size()==recv_slab.size()){
+   
     int count=0;
     for(typename NTSlab<T>::iterator iter=recv_slab.begin();
 	iter!=recv_slab.end();
@@ -438,6 +518,7 @@
 //       cout << "wrote exchange value to particle " << iter->getID() << endl;
       count++;
     }
+      
   }else{
     cerr << "rank = " << m_comm.rank() << ":ParallelParticleArray<T>::exchange_single ERROR: size mismatch in received data, recv_data.size()!=recv_slab.size() - (" << recv_data.size() << "!=" << recv_slab.size() << "). dir = " << dir << ", dist = " << dist << endl;
   }

=== modified file 'tml/comm/cart_comm.hpp'
--- tml/comm/cart_comm.hpp	2017-01-04 08:26:43 +0000
+++ tml/comm/cart_comm.hpp	2018-11-27 02:58:52 +0000
@@ -129,7 +129,7 @@
 void TML_CartComm::shift_cont_packed(T send_data,P& recv_data,int dir,int dist,int tag)
 {
   int source,dest;
-
+  
   if(dir<m_ndims){
     MPI_Cart_shift(m_comm,dir,dist,&source,&dest);
     sendrecv_cont_packed(send_data,recv_data,dest,source,false,tag);

=== modified file 'tml/comm/comm.hpp'
--- tml/comm/comm.hpp	2017-01-04 08:26:43 +0000
+++ tml/comm/comm.hpp	2018-11-27 02:58:52 +0000
@@ -389,7 +389,7 @@
 
   //send/receive data
   sendrecv_array(send_msg->buffer(),send_msg_size,recv_msg->buffer(),recv_msg_size,dest,source,tag+1024);
-
+  
   if(source!=MPI_PROC_NULL){ // if the source exists
     // extract nuber of items
     recv_nb_data=recv_msg->pop_int();