ESyS-Particle  4.0.1
RotParticleVi.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 __ROTPARTICLEVI_H
14 #define __ROTPARTICLEVI_H
15 
16 // -- project includes --
17 #include "Foundation/vec3.h"
18 #include "Foundation/Matrix3.h"
19 #include "Model/Particle.h"
20 #include "Foundation/Quaternion.h"
21 
22 template <class T> class ParallelParticleArray;
23 class AMPISGBufferRoot;
24 class AMPIBuffer;
25 class AField;
26 
27 //--- MPIincludes ---
28 #include <mpi.h>
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 
52 class CRotParticleVi : public CParticle
53 {
54  public: // types
55 
57  {
58  public:
59  exchangeType()
60  : m_pos(),
61  m_initPos(),
62  m_vel(),
63  m_angVel(),
64  m_angVel_t(),
65  m_quat()
66  {
67  }
68 
70  const Vec3 &pos,
71  const Vec3 &initPos,
72  const Vec3 &vel,
73  const Vec3 &AngVel,
74  const Vec3 &currAngVel,
75  const Quaternion &quat
76  )
77  : m_pos(pos),
78  m_initPos(initPos),
79  m_vel(vel),
80  m_angVel(AngVel),
81  m_angVel_t(currAngVel),
82  m_quat(quat)
83  {
84  }
85  public:
86  Vec3 m_pos;
87  Vec3 m_initPos;
88  Vec3 m_vel;
89  Vec3 m_angVel;
90  Vec3 m_angVel_t ;
91  Quaternion m_quat;
92 
93  friend class TML_PackedMessageInterface;
94  };
95  typedef double (CRotParticleVi::* ScalarFieldFunction)() const;
96  typedef Vec3 (CRotParticleVi::* VectorFieldFunction)() const;
97 
98  protected:
99 
100  Quaternion m_quat;
101  Quaternion m_initquat;
102  Vec3 m_angVel; // ! Angular velocity at time t -0.5*dt
103  Vec3 m_angVel_t ;
105  double m_inertRot;
106  double m_div_inertRot;
107  bool m_is_dynamic;
108 
109  public:
110  CRotParticleVi();
111  CRotParticleVi(const esys::lsm::SimpleParticleData &particleData);
113  double rad,
114  double mass,
115  const Vec3& pos,
116  const Vec3& vel,
117  const Vec3& force,
118  int id,
119  bool is_dyn
120  );
122  double rad,
123  double mass,
124  const Vec3& pos,
125  const Vec3& vel,
126  const Vec3& force,
127  int id,
128  Quaternion& quat,
129  double inertRot,
130  const Vec3& moment,
131  const Vec3& angvel,
132  const Vec3& angvel_t
133  );
135  double rad,
136  double mass,
137  const Vec3& pos,
138  const Vec3& oldpos,
139  const Vec3& initpos,
140  const Vec3& vel,
141  const Vec3& force,
142  int id,
143  const Quaternion& quat,
144  const Quaternion& initquat,
145  double inertRot,
146  const Vec3& moment,
147  const Vec3& angvel,
148  const Vec3& angvel_t
149  );
150 
151  CRotParticleVi(const CParticle &p);
152 
153  virtual ~CRotParticleVi(){};
154 
155  static int getPackSize();
156  static ScalarFieldFunction getScalarFieldFunction(const string&);
157  static VectorFieldFunction getVectorFieldFunction(const string&);
158 
159  // Need to define this for template class using forAllParticles call in Parallel/SubLattice.hpp.
160  Vec3 getDisplacement() const {return CParticle::getDisplacement();};
161  void resetDisplacement() {CParticle::resetDisplacement();}
162 
163  inline const Vec3 &getAngVel () const { return m_angVel;}
164  inline const Vec3 getAngVel_t () const { return m_angVel_t ;} ;
165  inline void setAngVel_t (const Vec3 &v) { m_angVel_t = v;} ;
166  inline Vec3 getAngVelNR () const { return m_angVel_t;} // for field functions
167  inline void setAngVel(const Vec3 &V) { m_angVel = V;}
168  inline Quaternion getInitQuat() const { return m_initquat;}
169  inline Quaternion getQuat() const { return m_quat;}
170  inline void setQuat(const Quaternion &q) {m_quat = q;}
171  inline double getInertRot () const { return m_inertRot; }
172  inline void setInertRot (double inertRot)
173  {
174  m_inertRot = inertRot;
175  m_div_inertRot = 1.0/m_inertRot;
176  }
177  inline double getInvInertRot () const { return m_div_inertRot; }
178  inline Vec3 getMoment() const {return m_moment;}
179  inline void setMoment(const Vec3 &moment) {m_moment = moment;}
180  Vec3 getAngVector()const ;
181  void applyMoment( const Vec3&);
182  void integrate(double);
183  void zeroForce();
184 
185 //wycnewadded
186  virtual void zeroHeat(){} ;
187  virtual void integrateTherm(double){} ;
188  virtual void setTemperature(double){} ;
189  virtual void setCp(double){} ;
190  virtual void setThermExpansion0(double){} ;
191  virtual void setThermExpansion1(double){} ;
192  virtual void setThermExpansion2(double){} ;
193  virtual void thermExpansion(){} ;
194  virtual double get_y() {return m_pos.Y(); } ;
195 
196 
197  void rescale();
198  void setCircular(const Vec3& cv);
199  double getAngularKineticEnergy() const {return (0.2*m_mass*m_rad*m_rad)*(m_angVel_t*m_angVel_t);}
200  double getLinearKineticEnergy() const {return (0.5*m_mass)*(m_vel*m_vel);}
201  double getKineticEnergy() const {return getLinearKineticEnergy() + getAngularKineticEnergy();}
202  void writeAsDXLine(ostream&,int slid=0);
203  virtual void setNonRot() {m_inertRot=0.0;m_div_inertRot=0.0;};
204 
205  inline Quaternion getQuatFromRotVec(const Vec3 &vec) const
206  {
207  const double angle = vec.norm();
208  const double halfAngle = angle/2.0;
209  if (halfAngle > 0.0) {
210  return Quaternion(cos(halfAngle), (vec)*(sin(halfAngle)/angle));
211  }
212  return Quaternion();
213  }
214  void rotateBy(const Vec3 &vec) {m_quat = getQuatFromRotVec(vec)*m_quat;}
215  void rotateTo(const Vec3 &vec) {m_quat = getQuatFromRotVec(vec);}
216 
217  static map<string,AField*> generateFields(ParallelParticleArray<CRotParticleVi>*);
218 
219  friend ostream& operator<<(ostream&, const CRotParticleVi&);
220  void print(){cout << *this << endl << flush;};
221 
222  // -- checkpointing --
223  virtual void saveSnapShotData(std::ostream& oStream);
224  virtual void saveCheckPointData(std::ostream& oStream);
225  virtual void loadCheckPointData(std::istream& iStream);
226 
229 
230  // stress
231  inline double sigma_xx_2D() const {return m_sigma(0,0)/(M_PI*m_rad*m_rad);};
232  inline double sigma_xy_2D() const {return m_sigma(0,1)/(M_PI*m_rad*m_rad);};
233  inline double sigma_yy_2D() const {return m_sigma(1,1)/(M_PI*m_rad*m_rad);};
234 // inline double sigma_d() const;
235  static void get_type() {cout <<" CRotParticleVi" ;};
236  friend class TML_PackedMessageInterface;
237 
238  template <typename TmplVisitor>
239  void visit(TmplVisitor &visitor)
240  {
241  visitor.visitRotParticleVi(*this);
242  }
243 
244 };
245 
246 #endif //__ROTPARTICLEVI_H