ESyS-Particle  4.0.1
RotParticle.h
1 
2 // //
3 // Copyright (c) 2003-2011 by The University of Queensland //
4 // Earth Systems Science Computational Centre (ESSCC) //
5 // http://www.uq.edu.au/esscc //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 #ifndef __ROTPARTICLE_H
14 #define __ROTPARTICLE_H
15 
16 //--- MPIincludes ---
17 #include <mpi.h>
18 
19 // -- project includes --
20 #include "Foundation/vec3.h"
21 #include "Foundation/Matrix3.h"
22 #include "Model/Particle.h"
23 #include "Foundation/Quaternion.h"
24 
25 template <class T> class ParallelParticleArray;
26 class AMPISGBufferRoot;
27 class AMPIBuffer;
28 
29 
30 //--- STL includes ---
31 #include <map>
32 #include <vector>
33 #include <utility>
34 #include <string>
35 
36 using std::map;
37 using std::vector;
38 using std::pair;
39 using std::string;
40 
41 namespace esys
42 {
43  namespace lsm
44  {
45  class SimpleParticleData;
46  }
47 }
48 
53 class CRotParticle : public CParticle
54 {
55  public: // types
56 
58  {
59  public:
60  exchangeType()
61  : m_pos(),
62  m_initPos(),
63  m_vel(),
64  m_angVel(),
65  m_quat()
66  {
67  m_is_dynamic=true;
68  m_is_rot=true;
69  }
70 
72  const Vec3 &pos,
73  const Vec3 &initPos,
74  const Vec3 &vel,
75  const Vec3 &currAngVel,
76  const Quaternion &quat,
77  const bool is_dyn,
78  const bool is_rot
79  )
80  : m_pos(pos),
81  m_initPos(initPos),
82  m_vel(vel),
83  m_angVel(currAngVel),
84  m_quat(quat)
85  {
86  m_is_dynamic=is_dyn;
87  m_is_rot=is_rot;
88  }
89  public:
90  Vec3 m_pos;
91  Vec3 m_initPos;
92  Vec3 m_vel;
93  Vec3 m_angVel;
94  Quaternion m_quat;
95  bool m_is_dynamic;
96  bool m_is_rot;
97 
98  friend class TML_PackedMessageInterface;
99  };
100  typedef double (CRotParticle::* ScalarFieldFunction)() const;
101  typedef Vec3 (CRotParticle::* VectorFieldFunction)() const;
102 
103  protected:
104 
105  Quaternion m_quat;
106  Quaternion m_initquat;
107  Vec3 m_angVel;
109  double m_inertRot;
110  double m_div_inertRot;
111  bool m_is_rot;
112 
113  void setMoment(const Vec3 &moment) {m_moment = moment;}
114 
115  public:
116  CRotParticle();
117  CRotParticle(const esys::lsm::SimpleParticleData &particleData);
118  CRotParticle(const CParticle &particle);
119  CRotParticle(
120  double rad,
121  double mass,
122  const Vec3& pos,
123  const Vec3& vel,
124  const Vec3& force,
125  int id,
126  bool is_dyn,
127  bool is_rot
128  );
129  CRotParticle(
130  double rad,
131  double mass,
132  const Vec3& pos,
133  const Vec3& vel,
134  const Vec3& force,
135  int id,
136  Quaternion& quat,
137  double inertRot,
138  const Vec3& moment,
139  const Vec3& angvel,
140  bool is_dyn,
141  bool is_rot
142  );
143  CRotParticle(
144  double rad,
145  double mass,
146  const Vec3& pos,
147  const Vec3& oldpos,
148  const Vec3& initpos,
149  const Vec3& vel,
150  const Vec3& force,
151  int id,
152  const Quaternion& quat,
153  const Quaternion& initquat,
154  double inertRot,
155  const Vec3& moment,
156  const Vec3& angvel,
157  bool is_dyn,
158  bool is_rot
159  );
160 
161  virtual ~CRotParticle(){};
162 
163  static int getPackSize();
164  static ScalarFieldFunction getScalarFieldFunction(const string&);
165  static VectorFieldFunction getVectorFieldFunction(const string&);
166 
167  // Need to define this for template class using forAllParticles call in Parallel/SubLattice.hpp.
168  Vec3 getDisplacement() const {return CParticle::getDisplacement();};
169  void resetDisplacement() {CParticle::resetDisplacement();}
170 
171  inline const Vec3 &getAngVel () const { return m_angVel;}
172  inline Vec3 getAngVelNR () const { return m_angVel;} // for field functions
173  inline void setAngVel(const Vec3 &V) { m_angVel = V;}
174  inline Quaternion getInitQuat() const { return m_initquat;}
175  inline Quaternion getQuat() const { return m_quat;}
176  inline void setQuat(const Quaternion &quat) { m_quat = quat;}
177  inline double getInertRot () const { return m_inertRot; }
178  inline void setInertRot (double inertRot)
179  {
180  m_inertRot = inertRot;
181  m_div_inertRot = 1.0/m_inertRot;
182  }
183  inline double getInvInertRot () const { return m_div_inertRot; }
184  inline Vec3 getMoment() const {return m_moment;}
185  void applyMoment( const Vec3&);
186  void integrate(double);
187  void integrateTherm(double dt){}
188  virtual void thermExpansion() {}
189  void zeroForce();
190  virtual void zeroHeat(){}
191  void rescale();
192  void setCircular(const Vec3& cv);
193  double getAngularKineticEnergy() const {return (0.2*m_mass*m_rad*m_rad)*(m_angVel*m_angVel);}
194  double getLinearKineticEnergy() const {return (0.5*m_mass)*(m_vel*m_vel);}
195  double getKineticEnergy() const {return getLinearKineticEnergy() + getAngularKineticEnergy();}
196  void writeAsDXLine(ostream&,int slid=0);
197 
198  // switching on/off dynamic behaviour
199  virtual void setNonDynamic() {m_is_dynamic=false;m_is_rot=false;};
200  virtual void setNonDynamicLinear() {m_is_dynamic=false;};
201  virtual void setNonDynamicRot() {m_is_rot=false;};
202 
203  inline Quaternion getQuatFromRotVec(const Vec3 &vec) const
204  {
205  const double angle = vec.norm();
206  const double halfAngle = angle/2.0;
207  if (halfAngle > 0.0) {
208  return Quaternion(cos(halfAngle), (vec)*(sin(halfAngle)/angle));
209  }
210  return Quaternion();
211  }
212  void rotateBy(const Vec3 &vec) {m_quat = getQuatFromRotVec(vec)*m_quat;}
213  void rotateTo(const Vec3 &vec) {m_quat = getQuatFromRotVec(vec);}
214 
215  friend ostream& operator<<(ostream&, const CRotParticle&);
216  void print(){cout << *this << endl << flush;};
217 
218  // -- checkpointing --
219  virtual void saveSnapShotData(std::ostream& oStream);
220  virtual void saveCheckPointData(std::ostream& oStream);
221  virtual void loadCheckPointData(std::istream& iStream);
222 
225 
226  template <typename TmplVisitor>
227  void visit(TmplVisitor &visitor)
228  {
229  visitor.visitRotParticle(*this);
230  }
231 
232  // stress
233  inline double sigma_xx_2D() const {return m_sigma(0,0)/(M_PI*m_rad*m_rad);};
234  inline double sigma_xy_2D() const {return m_sigma(0,1)/(M_PI*m_rad*m_rad);};
235  inline double sigma_yy_2D() const {return m_sigma(1,1)/(M_PI*m_rad*m_rad);};
236 // inline double sigma_d() const;
237  static void get_type() {cout <<" CRotParticle" ;};
238  friend class TML_PackedMessageInterface;
239 };
240 
241 #endif //__ROTPARTICLE_H