|
HG-MD
1
|
|
HG-MD
1
|
00001 #ifndef PARTICLE_H 00002 #define PARTICLE_H 00003 00004 //#define USE_SIMPLE_VERLET_INTEGRATION 1 00005 00006 00007 //This is a header with some extra standard maths function and constants that are required by the code 00008 #include "ExtendedMath.h" 00009 00010 #define sqr(a) ((a)*(a)) 00011 #define cubic(a) ((a)*(a)*(a)) 00012 00013 #include "CTangentialSpring.h" 00014 #include "CDeltaMax.h" 00015 #include "CSpecies.h" 00016 00017 #include <stdio.h> 00018 #include <stdlib.h> 00019 #include <iostream> 00020 #include <iomanip> 00021 using namespace std; 00025 00026 class Particle 00027 { 00028 public: 00030 Particle() { 00031 position.set_zero(); 00032 velocity.set_zero(); 00033 displacement.set_zero(); 00034 force.set_zero(); 00035 radius=1.0; 00036 angle.set_zero(); 00037 angularVelocity.set_zero(); 00038 torque.set_zero(); 00039 mass = invMass = 1.0; 00040 inertia = invInertia = 1.0; 00041 HGRID_NextObject=NULL; 00042 00043 tangentialSprings.reset(); 00044 deltaMaxs.reset(); 00045 indSpecies = 0; 00046 periodicFromParticle=NULL; 00047 00048 #ifdef CONSTUCTOR_OUTPUT 00049 cerr << "Particle() finished"<<endl; 00050 #endif 00051 } 00052 00053 00054 00056 Particle(const Particle &p) { 00057 position = p.position; 00058 velocity = p.velocity; 00059 displacement = p.displacement; 00060 force = p.force; 00061 radius = p.radius; 00062 angle = p.angle; 00063 angularVelocity = p.angularVelocity; 00064 torque = p.torque; 00065 mass = p.get_Mass(); 00066 invMass = p.get_InvMass(); 00067 inertia = p.get_Inertia(); 00068 invInertia = p.get_InvInertia(); 00069 00070 HGRID_NextObject = p.HGRID_NextObject; 00071 HGRID_x = p.HGRID_x; 00072 HGRID_y = p.HGRID_y; 00073 HGRID_z = p.HGRID_z; 00074 HGRID_Level = p.HGRID_Level; 00075 00076 tangentialSprings = p.tangentialSprings; 00077 deltaMaxs = p.deltaMaxs; 00078 00079 indSpecies = p.indSpecies; 00080 periodicFromParticle=p.periodicFromParticle; 00081 00082 #ifdef CONSTUCTOR_OUTPUT 00083 cerr << "Particle(Particle &p) finished"<<endl; 00084 #endif 00085 } 00086 00088 virtual ~Particle() 00089 { 00090 #ifdef CONSTUCTOR_OUTPUT 00091 cerr << "virtual ~Particle finished"<<endl; 00092 #endif 00093 } 00094 00095 //~ ///This get_mass function ensures that the mass is set 00096 //~ Mdouble get_mass(vector<CSpecies>& Species) { 00097 //~ compute_particle_mass(Species); 00098 //~ return mass; 00099 //~ } 00100 00102 virtual Particle* copy() 00103 { 00104 return new Particle(*this); 00105 } 00106 00108 Mdouble get_Volume(vector<CSpecies>& Species) const 00109 { 00110 switch(Species[indSpecies].dim_particle) 00111 { 00112 case 3 : 00113 return (4.0/3.0 * constants::pi * radius * radius * radius); 00114 case 2 : 00115 return (constants::pi * radius * radius); 00116 case 1 : 00117 return (2.0 * radius); 00118 default : 00119 { 00120 cerr<<"In get_Volume(vector<CSpecies>& Species) the dimension of the particle is not set"<<endl; 00121 exit(-1); 00122 } 00123 } 00124 00125 } 00126 00128 void fixParticle(){ 00129 mass=1e20; invMass=0.0; velocity.set_zero(); inertia=1e20; invInertia=0.0; angularVelocity.set_zero(); 00130 #ifdef USE_SIMPLE_VERLET_INTEGRATION 00131 position = PrevPosition; 00132 #endif 00133 } 00134 00136 bool isFixed(){return (invMass==0.0);} 00137 00139 void unfix(vector<CSpecies>& Species){ 00140 invMass=1.0; 00141 compute_particle_mass(Species); 00142 } 00143 00145 void compute_particle_mass(vector<CSpecies>& Species) { 00146 if (!isFixed()) { 00147 switch(Species[indSpecies].dim_particle) 00148 { 00149 case 3 : 00150 { 00151 set_Mass(4.0/3.0 * constants::pi * radius * radius * radius * Species[indSpecies].rho); 00152 set_inertia(.4 * get_Mass() * sqr(radius)); 00153 break; 00154 } 00155 case 2 : 00156 { 00157 set_Mass(constants::pi * radius * radius * Species[indSpecies].rho); 00158 set_inertia(.5 * get_Mass() * sqr(radius)); 00159 break; 00160 } 00161 case 1 : 00162 { 00163 set_Mass(2.0 * radius * Species[indSpecies].rho); 00164 set_inertia(0.0); 00165 break; 00166 } 00167 default : 00168 { 00169 cerr<<"In compute_particle_mass(vector<CSpecies>& Species) the dimension of the particle is not set"<<endl; 00170 exit(-1); 00171 } 00172 } 00173 } 00174 } 00175 00177 friend inline std::ostream& operator<<(std::ostream& os, const Particle &p) 00178 { 00179 os << p.position << " " << p.velocity << " " << p.radius << " " << p.angle << " " << p.angularVelocity << " " << p.invMass << " " << p.invInertia << " " << p.tangentialSprings; 00181 return os; 00182 } 00183 00185 friend inline std::istream& operator>>(std::istream& is, Particle &p) 00186 { 00187 is >> p.position >> p.velocity >> p.radius >> p.angle >> p.angularVelocity >> p.invMass >> p.invInertia >> p.tangentialSprings; //add this for the new restart 00188 if (p.invMass) p.mass = 1.0/p.invMass; else p.mass = 1e20; 00189 if (p.invInertia) p.inertia = 1.0/p.invInertia; else p.inertia = 1e20; 00190 return is; 00191 } 00192 00193 00195 void print(std::ostream& os) { 00196 os << "Particle( Position:" << position 00197 << ", Velocity:" << velocity 00198 << ", Radius:" << radius 00199 << ", Mass:" << mass 00200 //<< ", InvMass:" << invMass 00201 //<< ", Angle:" << Angle 00202 << ", AngularVelocity:" << angularVelocity 00203 << ", Inertia:" << inertia 00204 << ")"; 00205 } 00206 00207 void print_HGRID(std::ostream& os) { 00208 os << "Particle( HGRID_Level:" << HGRID_Level 00209 << ", HGRID_x:" << HGRID_x 00210 << ", HGRID_y:" << HGRID_y 00211 << ", HGRID_z:" << HGRID_z 00212 << ")"; 00213 } 00214 00215 int get_HGRID_Level() {return HGRID_Level;} 00216 Particle* get_HGRID_NextObject() {return HGRID_NextObject;} 00217 int get_HGRID_x() {return HGRID_x;} 00218 int get_HGRID_y() {return HGRID_y;} 00219 int get_HGRID_z() {return HGRID_z;} 00220 int get_Index() {return index;} 00221 Mdouble get_Inertia() const {return inertia;} 00222 Mdouble get_InvInertia() const {return invInertia;} 00223 Mdouble get_InvMass() const {return invMass;} 00224 Mdouble get_KineticEnergy() {return 0.5 * mass * velocity.GetLength2();} 00225 Mdouble get_Mass() const {return mass;} 00226 Particle* get_PeriodicFromParticle() {return periodicFromParticle;} 00227 Vec3D& get_Position() {return position;} 00228 Mdouble get_Radius() {return radius;} 00229 int get_Species() {return indSpecies;} 00230 int get_IndSpecies() {return indSpecies;} 00231 Vec3D& get_Velocity() {return velocity;} 00232 Vec3D& get_Angle() {return angle;} 00233 Vec3D& get_AngularVelocity() {return angularVelocity;} 00234 Vec3D& get_Force() {return force;} 00235 Vec3D& get_Torque() {return torque;} 00236 Vec3D& get_Displacement() {return displacement;} 00237 00238 00239 void set_inertia (Mdouble _new) {if (_new>=0){inertia=_new; invInertia=1.0/_new;} else { cerr << "Error in set_inertia ("<<_new<<")"<< endl; exit(-1); }} 00240 void set_periodicFromParticle (Particle* _new) {periodicFromParticle=_new;} 00241 void set_species (Mdouble _new) {indSpecies=_new;} 00242 void set_Index (int _new) {index=_new;} 00243 void set_HGRID_x (int _new) {HGRID_x=_new;} 00244 void set_HGRID_y (int _new) {HGRID_y=_new;} 00245 void set_HGRID_z (int _new) {HGRID_z=_new;} 00246 void set_HGRID_Level (int _new) {HGRID_Level=_new;} 00247 void set_HGRID_NextObject (Particle* _new) {HGRID_NextObject=_new;} 00248 void set_Radius (Mdouble _new) {radius=_new;} 00249 void set_IndSpecies (int _new) {indSpecies=_new;} 00250 void set_Mass (Mdouble _new) {if(_new>=0.0) {if(invMass){mass=_new;invMass=1.0/_new;}} else { cerr << "Error in set_Mass(" << _new << ")" << endl; exit(-1); }} //InvMass=0 is a flag for a fixed particle 00251 void set_Velocity (Vec3D _new) {velocity=_new;} 00252 void set_Position (Vec3D _new) {position=_new;} 00253 void set_Angle (Vec3D _new) {angle=_new;} 00254 void set_AngularVelocity (Vec3D _new) {angularVelocity=_new;} 00255 void set_Force (Vec3D _new) {force=_new;} 00256 void set_Torque (Vec3D _new) {torque=_new;} 00257 00258 void move (Vec3D _new) {position+=_new;} 00259 void accelerate (Vec3D _new) {velocity+=_new;} 00260 void rotate (Vec3D _new) {angle+=_new;} 00261 void angularAccelerate (Vec3D _new) {angularVelocity+=_new;} 00262 void add_Force (Vec3D _new) {force+=_new;} 00263 void add_Torque (Vec3D _new) {torque+=_new;} 00264 00265 00266 00267 00268 #ifdef USE_SIMPLE_VERLET_INTEGRATION 00269 void set_PrevPosition(Vec3D _new) {prevPosition=_new;} 00270 void prevMove(Vec3D _new) {prevPosition+=_new;} 00271 Vec3D& get_PrevPosition() {return prevPosition;} 00272 #endif 00273 00274 CTangentialSprings& get_TangentialSprings() {return tangentialSprings;} 00275 CDeltaMaxs& get_DeltaMaxs() {return deltaMaxs;} 00276 00277 private: 00278 00280 int HGRID_x, HGRID_y, HGRID_z; 00281 int HGRID_Level; 00282 Particle* HGRID_NextObject; 00283 00285 int index; 00286 Mdouble mass; 00287 Mdouble invMass; 00288 Mdouble inertia; 00289 Mdouble invInertia; 00290 Mdouble radius; 00291 Particle* periodicFromParticle; 00292 int indSpecies; 00293 00294 Vec3D angle; 00295 Vec3D angularVelocity; 00296 Vec3D position; 00297 Vec3D velocity; 00298 00299 Vec3D force; 00300 Vec3D torque; 00301 00302 Vec3D displacement; 00303 #ifdef USE_SIMPLE_VERLET_INTEGRATION 00304 Vec3D prevPosition; 00305 #endif 00306 00307 CTangentialSprings tangentialSprings; 00308 CDeltaMaxs deltaMaxs; 00309 }; 00310 00311 #endif
1.7.6.1
