 |
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
|
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
00001 /************************************************************************ 00002 * MechSys - Open Library for Mechanical Systems * 00003 * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand * 00004 * Copyright (C) 2009 Sergio Galindo * 00005 * * 00006 * This program is free software: you can redistribute it and/or modify * 00007 * it under the terms of the GNU General Public License as published by * 00008 * the Free Software Foundation, either version 3 of the License, or * 00009 * any later version. * 00010 * * 00011 * This program is distributed in the hope that it will be useful, * 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00014 * GNU General Public License for more details. * 00015 * * 00016 * You should have received a copy of the GNU General Public License * 00017 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00018 ************************************************************************/ 00019 00022 #ifndef MECHSYS_DEM_DOMAIN_H 00023 #define MECHSYS_DEM_DOMAIN_H 00024 00025 // Std lib 00026 #include <cmath> 00027 #include <stdlib.h> // for M_PI 00028 #include <iostream> 00029 #include <fstream> 00030 #include <set> 00031 #include <list> 00032 00033 // Hdf5 00034 #ifdef USE_HDF5 00035 #include <hdf5.h> 00036 #include <hdf5_hl.h> 00037 #endif 00038 00039 // Voro++ 00040 #include "src/voro++.cc" 00041 00042 // MechSys 00043 #include <mechsys/dem/interacton.h> 00044 #include <mechsys/util/array.h> 00045 #include <mechsys/util/util.h> 00046 #include <mechsys/util/numstreams.h> 00047 #include <mechsys/mesh/mesh.h> 00048 #include <mechsys/util/maps.h> 00049 #include <mechsys/util/stopwatch.h> 00050 #include <mechsys/util/tree.h> 00051 00052 namespace DEM 00053 { 00054 00055 #ifdef USE_THREAD 00056 void GlobalForce(Array<Interacton *> * I, double dt, size_t n, size_t Np) 00057 { 00058 size_t Ni = I->Size()/Np; 00059 size_t In = n*Ni; 00060 size_t Fn; 00061 n == Np-1 ? Fn = I->Size() : Fn = (n+1)*Ni; 00062 for (size_t i=In;i<Fn;i++) 00063 { 00064 (*I)[i]->CalcForce(dt); 00065 } 00066 } 00067 00068 void GlobalMove(Array<Particle *> * P, double dt, double & MaxD, size_t n, size_t Np) 00069 { 00070 size_t Ni = P->Size()/Np; 00071 size_t In = n*Ni; 00072 size_t Fn; 00073 n == Np-1 ? Fn = P->Size() : Fn = (n+1)*Ni; 00074 MaxD = 0.0; 00075 for (size_t i=In;i<Fn;i++) 00076 { 00077 (*P)[i]->Translate(dt); 00078 (*P)[i]->Rotate(dt); 00079 //double temp = (*P)[i]->MaxDisplacement(); 00080 //if (temp > MaxD) MaxD = temp; 00081 00082 } 00083 //std::cout << MaxD << " mu " << n << std::endl; 00084 } 00085 00086 #endif 00087 00088 00089 00090 class Domain 00091 { 00092 public: 00093 // typedefs 00094 typedef void (*ptFun_t) (Domain & Dom, void * UserData); 00095 00096 // Constructor 00097 Domain(void * UserData=NULL); 00098 00099 // Destructor 00100 ~Domain(); 00101 00102 // Particle generation 00103 void GenSpheres (int Tag, double L, size_t N, double rho, char const * Type, 00104 size_t Randomseed, double fraction, double RminFraction = 1.0); 00105 void GenRice (int Tag, double L, size_t N, double R, double rho, size_t Randomseed, double fraction); 00106 void GenBox (int InitialTag, double Lx, double Ly, double Lz, double R, double Cf, bool Cohesion=false); 00107 void GenOpenBox (int InitialTag, double Lx, double Ly, double Lz, double R, double Cf); 00108 void GenBoundingBox (int InitialTag, double R, double Cf,bool Cohesion=false); 00109 void GenBoundingPlane(int InitialTag, double R, double Cf,bool Cohesion=false); 00110 void GenFromMesh (Mesh::Generic & M, double R, double rho, bool cohesion=false, bool MC=true, double thickness = 0.0); 00111 void AddVoroPack (int Tag, double R, double Lx, double Ly, double Lz, size_t nx, size_t ny, size_t nz, 00112 double rho, bool Cohesion, bool Periodic,size_t Randomseed, double fraction, Vec3_t q = OrthoSys::O); 00113 // Single particle addition 00114 void AddSphere (int Tag, Vec3_t const & X, double R, double rho); 00115 void AddCube (int Tag, Vec3_t const & X, double R, double L, double rho, double Angle=0, Vec3_t * Axis=NULL); 00116 void AddTetra (int Tag, Vec3_t const & X, double R, double L, double rho, double Angle=0, Vec3_t * Axis=NULL); 00117 void AddDrill (int Tag, Vec3_t const & X, double R, double Lt, double Ll, double rho); 00118 void AddRice (int Tag, Vec3_t const & X, double R, double L, double rho, double Angle=0, Vec3_t * Axis=NULL); 00119 void AddPlane (int Tag, Vec3_t const & X, double R, double Lx,double Ly, double rho, double Angle=0, Vec3_t * Axis=NULL); 00120 void AddVoroCell (int Tag, voronoicell_neighbor & VC, double R, double rho, bool Erode, iVec3_t nv = iVec3_t(1,1,1)); 00121 void AddTorus (int Tag, Vec3_t const & X, Vec3_t const & N, double Rmax, double R, double rho); 00122 void AddCylinder (int Tag, Vec3_t const & X0, double R0, Vec3_t const & X1, double R1, double R, double rho); 00123 00124 // Methods 00125 void SetProps (Dict & D); 00126 void Initialize (double dt=0.0); 00127 void Solve (double tf, double dt, double dtOut, ptFun_t ptSetup=NULL, ptFun_t ptReport=NULL, 00128 char const * FileKey=NULL, bool RenderVideo=true, size_t Nproc=1); 00129 void WritePOV (char const * FileKey); 00130 void WriteBPY (char const * FileKey); 00131 #ifdef USE_HDF5 00132 void Save (char const * FileKey); 00133 void Load (char const * FileKey); 00134 #endif 00135 void BoundingBox (Vec3_t & minX, Vec3_t & maxX); 00136 void Center (Vec3_t C = Vec3_t(0.0,0.0,0.0)); 00137 void ResetInteractons (); 00138 void ResetDisplacements(); 00139 double MaxDisplacement (); 00140 void ResetContacts (); 00141 void EnergyOutput (size_t IdxOut, std::ostream & OutFile); 00142 void GetGSD (Array<double> & X, Array<double> & Y, Array<double> & D, size_t NDiv=10) const; 00143 void Clusters (); 00144 void DelParticles (Array<int> const & Tags); 00145 00146 // Access methods 00147 Particle * GetParticle (int Tag, bool Check=true); 00148 Particle const & GetParticle (int Tag, bool Check=true) const; 00149 void GetParticles (int Tag, Array<Particle*> & P); 00150 00151 // Auxiliar methods 00152 void LinearMomentum (Vec3_t & L); 00153 void AngularMomentum (Vec3_t & L); 00154 double CalcEnergy (double & Ekin, double & Epot); 00155 00156 // Data 00157 bool Initialized; 00158 bool Finished; 00159 Array<Particle*> Particles; 00160 Array<Interacton*> Interactons; 00161 Array<CInteracton*> CInteractons; 00162 Array<BInteracton*> BInteractons; 00163 Vec3_t CamPos; 00164 double Time; 00165 double Evis; 00166 double Efric; 00167 double Wext; 00168 double Vs; 00169 double Ms; 00170 double Alpha; 00171 void * UserData; 00172 String FileKey; 00173 size_t idx_out; 00174 set<pair<Particle *, Particle *> > Listofpairs; 00175 Array<Array <int> > Listofclusters; 00176 00177 #ifdef USE_BOOST_PYTHON 00178 void PyAddSphere (int Tag, BPy::tuple const & X, double R, double rho) { AddSphere (Tag,Tup2Vec3(X),R,rho); } 00179 void PyAddCube (int Tag, BPy::tuple const & X, double R, double L, double rho, double Ang, BPy::tuple const & Ax) { Vec3_t a(Tup2Vec3(Ax)); AddCube (Tag,Tup2Vec3(X),R,L,rho,Ang,&a); } 00180 void PyAddTetra (int Tag, BPy::tuple const & X, double R, double L, double rho, double Ang, BPy::tuple const & Ax) { Vec3_t a(Tup2Vec3(Ax)); AddTetra (Tag,Tup2Vec3(X),R,L,rho,Ang,&a); } 00181 void PyAddRice (int Tag, BPy::tuple const & X, double R, double L, double rho, double Ang, BPy::tuple const & Ax) { Vec3_t a(Tup2Vec3(Ax)); AddRice (Tag,Tup2Vec3(X),R,L,rho,Ang,&a); } 00182 void PyAddPlane (int Tag, BPy::tuple const & X, double R, double Lx,double Ly, double rho, double Ang, BPy::tuple const & Ax) { Vec3_t a(Tup2Vec3(Ax)); AddPlane (Tag,Tup2Vec3(X),R,Lx,Ly,rho,Ang,&a); } 00183 void PySetCamPos (BPy::tuple const & PyCamPos) { CamPos = Tup2Vec3(PyCamPos); } 00184 void PyGetParticles(BPy::list & P) 00185 { 00186 for (size_t i=0; i<Particles.Size(); ++i) 00187 { 00188 BPy::list p,V,E,F; 00189 double radius = Particles[i]->PyGetFeatures (V, E, F); 00190 p.append (radius); 00191 p.append (V); 00192 p.append (E); 00193 p.append (F); 00194 P.append (p); 00195 } 00196 } 00197 void PyGetGSD (BPy::list & X, BPy::list & Y, BPy::list & D, int NDiv=10) 00198 { 00199 Array<double> x, y, d; 00200 GetGSD (x, y, d, NDiv); 00201 for (size_t i=0; i<x.Size(); ++i) 00202 { 00203 X.append (x[i]); 00204 Y.append (y[i]); 00205 } 00206 } 00207 #endif 00208 00209 }; 00210 00212 00213 00214 // Constructor & Destructor 00215 00216 inline Domain::Domain (void * UD) 00217 : Initialized(false), Time(0.0), Alpha(0.05), UserData(UD) 00218 { 00219 CamPos = 1.0, 2.0, 3.0; 00220 } 00221 00222 inline Domain::~Domain () 00223 { 00224 for (size_t i=0; i<Particles.Size(); ++i) if (Particles [i]!=NULL) delete Particles [i]; 00225 for (size_t i=0; i<Interactons.Size(); ++i) if (Interactons[i]!=NULL) delete Interactons[i]; 00226 } 00227 00228 // Particle generation 00229 00230 inline void Domain::GenSpheres (int Tag, double L, size_t N, double rho,char const * Type, size_t Randomseed, double fraction, double RminFraction) 00231 { 00232 // find radius from the edge's length 00233 Util::Stopwatch stopwatch; 00234 printf("\n%s--- Generating packing of spheres -----------------------------------------------%s\n",TERM_CLR1,TERM_RST); 00235 srand(Randomseed); 00236 double R = L/(2.0*N); 00237 if (strcmp(Type,"Normal")==0) 00238 { 00239 for (size_t n=0; n<N*N*N; ++n) 00240 { 00241 Vec3_t pos(-L/2.0+R, -L/2.0+R, -L/2.0+R); 00242 size_t i = (n%N); 00243 size_t j = (n/N)%N; 00244 size_t k = (n/(N*N)); 00245 pos += Vec3_t(2.0*i*R, 2.0*j*R, 2.0*k*R); 00246 if (rand()<fraction*RAND_MAX) AddSphere (Tag,pos,R*RminFraction+(1.0*rand())/RAND_MAX*(R-R*RminFraction),rho); 00247 } 00248 } 00249 else if (strcmp(Type,"HCP")==0) 00250 { 00251 size_t nx = N; 00252 size_t ny = int(L/(sqrt(3.0)*R)); 00253 size_t nz = int(L/(sqrt(8.0/3.0)*R)); 00254 for (size_t k = 0; k < nz; k++) 00255 { 00256 for (size_t j = 0; j < ny; j++) 00257 { 00258 Vec3_t X; 00259 if (k%2==0) X = Vec3_t(-2*R-L/2.0,R-L/2.0,2*R-L/2.0+k*sqrt(8.0/3.0)*R); 00260 else X = Vec3_t(-R-L/2.0,R+sqrt(1.0/3.0)*R-L/2.0,2*R-L/2.0+k*sqrt(8.0/3.0)*R); 00261 if (j%2==0) X += Vec3_t(R,j*sqrt(3.0)*R,0.0); 00262 else X += Vec3_t(0.0,j*sqrt(3.0)*R,0.0); 00263 for (size_t i = 0; i < nx; i++) 00264 { 00265 X += Vec3_t(2*R,0.0,0.0); 00266 if (rand()<fraction*RAND_MAX) AddSphere(Tag,X,R*RminFraction+(1.0*rand())/RAND_MAX*(R-R*RminFraction),rho); 00267 } 00268 } 00269 } 00270 } 00271 else throw new Fatal ("Right now there are only two possible packings available the Normal and the HCP, packing %s is not implemented yet",Type); 00272 printf("%s Num of particles = %zd%s\n",TERM_CLR2,Particles.Size(),TERM_RST); 00273 } 00274 00275 inline void Domain::GenRice (int Tag, double L, size_t N, double R, double rho, size_t Randomseed, double fraction) 00276 { 00277 Util::Stopwatch stopwatch; 00278 printf("\n%s--- Generating packing of 'rices' -----------------------------------------------%s\n",TERM_CLR1,TERM_RST); 00279 srand(Randomseed); 00280 double dL = L/N; 00281 for (size_t n=0; n<N*N*N; ++n) 00282 { 00283 Vec3_t pos(-L/2.0+dL, -L/2.0+dL, -L/2.0+dL); 00284 size_t i = (n%N); 00285 size_t j = (n/N)%N; 00286 size_t k = (n/(N*N)); 00287 pos += Vec3_t(2.0*i*dL, 2.0*j*dL, 2.0*k*dL); 00288 if (rand()<fraction*RAND_MAX) AddRice (Tag, pos, R, dL-2*R, rho); 00289 } 00290 printf("%s Num of particles = %zd%s\n",TERM_CLR2,Particles.Size(),TERM_RST); 00291 } 00292 00293 inline void Domain::GenBox (int InitialTag, double Lx, double Ly, double Lz, double R, double Cf, bool Cohesion) 00294 { 00295 /* +----------------+ 00296 * ,'| ,'| 00297 * ,' | ___ ,' | 00298 * z ,' |,'4,' [1],' | 00299 * | ,' |~~~ ,' | 00300 * ,+--y +'===============+' ,'| | 00301 * x' | ,'| | | |2| | 00302 * | |3| | | |,' | 00303 * | |,' +- - - | +- - - -+ 00304 * | ,' | ,' 00305 * | ,' [0] ___| ,' 00306 * | ,' ,'5,'| ,' 00307 * | ,' ~~~ | ,' 00308 * +----------------+' 00309 */ 00310 00311 00312 // add faces of box 00313 Vec3_t axis0(OrthoSys::e0); // rotation of face 00314 Vec3_t axis1(OrthoSys::e1); // rotation of face 00315 size_t IIndex = Particles.Size(); // First face index 00316 AddPlane (InitialTag, Vec3_t(Lx/2.0,0.0,0.0), R, Cf*Lz, Cf*Ly, 1.0, M_PI/2.0, &axis1); 00317 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00318 AddPlane (InitialTag-1, Vec3_t(-Lx/2.0,0.0,0.0), R, Cf*Lz, Cf*Ly, 1.0, 3.0*M_PI/2.0, &axis1); 00319 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00320 AddPlane (InitialTag-2, Vec3_t(0.0,Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, 3.0*M_PI/2.0, &axis0); 00321 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00322 AddPlane (InitialTag-3, Vec3_t(0.0,-Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, M_PI/2.0, &axis0); 00323 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00324 AddPlane (InitialTag-4, Vec3_t(0.0,0.0,Lz/2.0), R, Cf*Lx, Cf*Ly, 1.0); 00325 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00326 AddPlane (InitialTag-5, Vec3_t(0.0,0.0,-Lz/2.0), R, Cf*Lx, Cf*Ly, 1.0, M_PI, &axis0); 00327 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00328 00329 // define some tolerance for comparissions 00330 if (Cohesion) 00331 { 00332 double tol1 = 1.0e-8; 00333 double tol2 = 1.0e-3; 00334 for (size_t i=0;i<IIndex;i++) 00335 { 00336 Particle * P1 = Particles[i]; 00337 for (size_t j=IIndex;j<Particles.Size();j++) 00338 { 00339 Particle * P2 = Particles[j]; 00340 for (size_t k=0;k<P1->Faces.Size();k++) 00341 { 00342 Face * F1 = P1->Faces[k]; 00343 Vec3_t n1,c1; 00344 F1->Normal (n1); 00345 F1->Centroid(c1); 00346 Face * F2 = P2->Faces[0]; 00347 Vec3_t n2,c2; 00348 F2->Normal (n2); 00349 F2->Centroid(c2); 00350 Vec3_t n = 0.5*(n1-n2); 00351 n/=norm(n); 00352 if ((fabs(dot(n1,n2)+1.0)<tol1) 00353 &&(fabs(dot(c2-c1,n)-2*R)<tol2)) 00354 { 00355 BInteractons.Push(new BInteracton(P1,P2,k,1)); 00356 break; 00357 } 00358 } 00359 } 00360 } 00361 } 00362 } 00363 00364 inline void Domain::GenOpenBox (int InitialTag, double Lx, double Ly, double Lz, double R, double Cf) 00365 { 00366 /* +----------------+ 00367 * ,'| ,'| 00368 * ,' | ___ ,' | 00369 * z ,' |,'N,' [1],' | 00370 * | ,' |~~~ ,' | 00371 * ,+--y +'===============+' ,'| | 00372 * x' | ,'| | | |2| | 00373 * | |3| | | |,' | 00374 * | |,' +- - - | +- - - -+ 00375 * | ,' | ,' 00376 * | ,' [0] ___| ,' 00377 * | ,' ,'4,'| ,' 00378 * | ,' ~~~ | ,' 00379 * +----------------+' 00380 */ 00381 00382 00383 // Creates an open box without the top lid, acts as a container 00384 00385 // add faces of box 00386 Vec3_t axis0(OrthoSys::e0); // rotation of face 00387 Vec3_t axis1(OrthoSys::e1); // rotation of face 00388 AddPlane (InitialTag, Vec3_t(Lx/2.0,0.0,0.0), R, Cf*Lz, Cf*Ly, 1.0, M_PI/2.0, &axis1); 00389 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00390 AddPlane (InitialTag-1, Vec3_t(-Lx/2.0,0.0,0.0), R, Cf*Lz, Cf*Ly, 1.0, 3.0*M_PI/2.0, &axis1); 00391 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00392 AddPlane (InitialTag-2, Vec3_t(0.0,Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, 3.0*M_PI/2.0, &axis0); 00393 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00394 AddPlane (InitialTag-3, Vec3_t(0.0,-Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, M_PI/2.0, &axis0); 00395 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00396 AddPlane (InitialTag-4, Vec3_t(0.0,0.0,-Lz/2.0), R, Cf*Lx, Cf*Ly, 1.0, M_PI, &axis0); 00397 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00398 } 00399 00400 inline void Domain::GenBoundingBox (int InitialTag, double R, double Cf,bool Cohesion) 00401 { 00402 Center(); 00403 Vec3_t minX,maxX; 00404 BoundingBox(minX,maxX); 00405 GenBox(InitialTag, maxX(0)-minX(0)+2*R, maxX(1)-minX(1)+2*R, maxX(2)-minX(2)+2*R, R, Cf,Cohesion); 00406 } 00407 00408 inline void Domain::GenBoundingPlane (int InitialTag, double R, double Cf,bool Cohesion) 00409 { 00410 Center(); 00411 Vec3_t minX,maxX; 00412 BoundingBox(minX,maxX); 00413 Vec3_t axis0(OrthoSys::e0); // rotation of face 00414 Vec3_t axis1(OrthoSys::e1); // rotation of face 00415 size_t IIndex = Particles.Size(); // First face index 00416 double Lx = maxX(0)-minX(0)+2*R; 00417 double Ly = maxX(1)-minX(1)+2*R; 00418 double Lz = maxX(2)-minX(2)+2*R; 00419 AddPlane (InitialTag , Vec3_t(0.0,Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, 3.0*M_PI/2.0, &axis0); 00420 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00421 AddPlane (InitialTag-1, Vec3_t(0.0,-Ly/2.0,0.0), R, Cf*Lx, Cf*Lz, 1.0, M_PI/2.0, &axis0); 00422 Particles[Particles.Size()-1]->Initialize(Particles.Size()-1); 00423 00424 // define some tolerance for comparissions 00425 if (Cohesion) 00426 { 00427 double tol1 = 1.0e-8; 00428 double tol2 = 1.0e-3; 00429 for (size_t i=0;i<IIndex;i++) 00430 { 00431 Particle * P1 = Particles[i]; 00432 for (size_t j=IIndex;j<Particles.Size();j++) 00433 { 00434 Particle * P2 = Particles[j]; 00435 for (size_t k=0;k<P1->Faces.Size();k++) 00436 { 00437 Face * F1 = P1->Faces[k]; 00438 Vec3_t n1,c1; 00439 F1->Normal (n1); 00440 F1->Centroid(c1); 00441 Face * F2 = P2->Faces[0]; 00442 Vec3_t n2,c2; 00443 F2->Normal (n2); 00444 F2->Centroid(c2); 00445 Vec3_t n = 0.5*(n1-n2); 00446 n/=norm(n); 00447 if ((fabs(dot(n1,n2)+1.0)<tol1) 00448 &&(fabs(dot(c2-c1,n)-2*R)<tol2)) 00449 { 00450 BInteractons.Push(new BInteracton(P1,P2,k,1)); 00451 break; 00452 } 00453 } 00454 } 00455 } 00456 } 00457 } 00458 00459 inline void Domain::GenFromMesh (Mesh::Generic & M, double R, double rho, bool Cohesion, bool MC, double thickness) 00460 { 00461 // info 00462 Util::Stopwatch stopwatch; 00463 printf("\n%s--- Generating particles from mesh ----------------------------------------------%s\n",TERM_CLR1,TERM_RST); 00464 00465 size_t IIndex = Particles.Size(); 00466 00467 for (size_t i=0; i<M.Cells.Size(); ++i) 00468 { 00469 00470 Array<Vec3_t> V; // Array of vertices 00471 Array<Array <int> > E; // Array of edges 00472 Array<Array <int> > F; // array of faces 00473 if (M.NDim==3) 00474 { 00475 if (thickness > 0.0) throw new Fatal("Domain::GenFromMesh: Thickness should not be used in a 3D mesh"); 00476 Array<Mesh::Vertex*> const & verts = M.Cells[i]->V; 00477 size_t nverts = verts.Size(); 00478 00479 // verts 00480 V.Resize(nverts); 00481 for (size_t j=0; j<nverts; ++j) 00482 { 00483 V[j] = verts[j]->C; 00484 } 00485 00486 // edges 00487 size_t nedges = Mesh::NVertsToNEdges3D[nverts]; 00488 E.Resize(nedges); 00489 for (size_t j=0; j<nedges; ++j) 00490 { 00491 E[j].Push (Mesh::NVertsToEdge3D[nverts][j][0]); 00492 E[j].Push (Mesh::NVertsToEdge3D[nverts][j][1]); 00493 } 00494 00495 size_t nfaces = Mesh::NVertsToNFaces3D[nverts]; 00496 size_t nvperf = Mesh::NVertsToNVertsPerFace3D[nverts]; 00497 F.Resize(nfaces); 00498 for (size_t j=0; j<nfaces; ++j) 00499 { 00500 for (size_t k=0; k<nvperf; ++k) 00501 { 00502 // TODO: check if face is planar or not 00503 F[j].Push(Mesh::NVertsToFace3D[nverts][j][k]); 00504 } 00505 } 00506 } 00507 else if (M.NDim==2) 00508 { 00509 if (thickness <= 0.0) throw new Fatal("Domain::GenFromMesh: Thickness should be positive in a 2D mesh"); 00510 Array<Mesh::Vertex*> const & verts = M.Cells[i]->V; 00511 size_t nverts = verts.Size(); 00512 V.Resize(2*nverts); 00513 for (size_t j=0; j<nverts; ++j) 00514 { 00515 V[j] = verts[j]->C; 00516 V[j+nverts] = verts[j]->C + Vec3_t(0.0,0.0,thickness); 00517 } 00518 size_t nedges = 3*nverts; 00519 E.Resize(nedges); 00520 for (size_t j=0; j<nverts; ++j) 00521 { 00522 E[j].Push (Mesh::NVertsToEdge2D[nverts][j][0]); 00523 E[j].Push (Mesh::NVertsToEdge2D[nverts][j][1]); 00524 E[j+nverts].Push (Mesh::NVertsToEdge2D[nverts][j][0]+nverts); 00525 E[j+nverts].Push (Mesh::NVertsToEdge2D[nverts][j][1]+nverts); 00526 E[j+2*nverts].Push(j); 00527 E[j+2*nverts].Push(j+nverts); 00528 } 00529 size_t nfaces = nverts+2; 00530 F.Resize(nfaces); 00531 for (size_t j=0; j<nverts; ++j) 00532 { 00533 F[j].Push (Mesh::NVertsToEdge2D[nverts][j][0]); 00534 F[j].Push (Mesh::NVertsToEdge2D[nverts][j][1]); 00535 F[j].Push (Mesh::NVertsToEdge2D[nverts][j][1]+nverts); 00536 F[j].Push (Mesh::NVertsToEdge2D[nverts][j][0]+nverts); 00537 F[nverts].Push(nverts-1-j); 00538 F[nverts+1].Push(j+nverts); 00539 } 00540 } 00541 00542 double vol; // volume of the polyhedron 00543 Vec3_t CM; // Center of mass of the polyhedron 00544 Mat3_t It; // Inertia tensor of the polyhedron 00545 PolyhedraMP(V,F,vol,CM,It); 00546 Erosion(V,E,F,R); 00547 00548 // add particle 00549 Particles.Push (new Particle(M.Cells[i]->Tag, V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 00550 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 00551 if (!MC) 00552 { 00553 Particles[Particles.Size()-1]->x = CM; 00554 Particles[Particles.Size()-1]->Props.V = vol; 00555 Particles[Particles.Size()-1]->Props.m = vol*rho; 00556 Vec3_t I; 00557 Quaternion_t Q; 00558 Vec3_t xp,yp,zp; 00559 Eig(It,I,xp,yp,zp); 00560 CheckDestroGiro(xp,yp,zp); 00561 I *= rho; 00562 Q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2)); 00563 Q(1) = (yp(2)-zp(1))/(4*Q(0)); 00564 Q(2) = (zp(0)-xp(2))/(4*Q(0)); 00565 Q(3) = (xp(1)-yp(0))/(4*Q(0)); 00566 Q = Q/norm(Q); 00567 Particles[Particles.Size()-1]->I = I; 00568 Particles[Particles.Size()-1]->Q = Q; 00569 double Dmax = Distance(CM,V[0])+R; 00570 for (size_t i=1; i<V.Size(); ++i) 00571 { 00572 if (Distance(CM,V[i])+R > Dmax) Dmax = Distance(CM,V[i])+R; 00573 } 00574 Particles[Particles.Size()-1]->Ekin = 0.0; 00575 Particles[Particles.Size()-1]->Erot = 0.0; 00576 Particles[Particles.Size()-1]->Dmax = Dmax; 00577 Particles[Particles.Size()-1]->PropsReady = true; 00578 } 00579 } 00580 00581 Array<Array <int> > Neigh(Particles.Size()-IIndex); 00582 Array<Array <int> > FNeigh(Particles.Size()-IIndex); 00583 if(Cohesion) 00584 { 00585 M.FindNeigh(); 00586 //std::cout << M; 00587 for (size_t i=0; i<M.Cells.Size(); ++i) 00588 { 00589 for (Mesh::Neighs_t::const_iterator p=M.Cells[i]->Neighs.begin(); p!=M.Cells[i]->Neighs.end(); ++p) 00590 { 00591 Neigh[i].Push(p->second.second->ID); 00592 FNeigh[i].Push(p->second.first); 00593 } 00594 } 00595 for (size_t i=0; i<Neigh.Size(); ++i) 00596 { 00597 for (size_t j=0; j<Neigh[i].Size(); ++j) 00598 { 00599 size_t index = Neigh[Neigh[i][j]].Find(i); 00600 if ((size_t)Neigh[i][j]>i) BInteractons.Push(new BInteracton(Particles[i+IIndex],Particles[Neigh[i][j]+IIndex],FNeigh[i][j],FNeigh[Neigh[i][j]][index])); 00601 } 00602 } 00603 00604 } 00605 00606 // info 00607 printf("%s Num of particles = %zd%s\n",TERM_CLR2,Particles.Size(),TERM_RST); 00608 } 00609 00610 inline void Domain::AddVoroPack (int Tag, double R, double Lx, double Ly, double Lz, size_t nx, size_t ny, size_t nz, double rho 00611 , bool Cohesion, bool Periodic,size_t Randomseed, double fraction, Vec3_t qin) 00612 { 00613 // info 00614 Util::Stopwatch stopwatch; 00615 printf("\n%s--- Adding Voronoi particles packing --------------------------------------------%s\n",TERM_CLR1,TERM_RST); 00616 00617 srand(Randomseed); 00618 const double x_min=-(nx*Lx/2.0), x_max=nx*Lx/2.0; 00619 const double y_min=-(ny*Ly/2.0), y_max=ny*Ly/2.0; 00620 const double z_min=-(nz*Lz/2.0), z_max=nz*Lz/2.0; 00621 container con(x_min,x_max,y_min,y_max,z_min,z_max,nx,ny,nz, Periodic,Periodic,Periodic,8); 00622 int n = 0; 00623 for (size_t i=0; i<nx; i++) 00624 { 00625 for (size_t j=0; j<ny; j++) 00626 { 00627 for (size_t k=0; k<nz; k++) 00628 { 00629 double x = x_min+(i+0.5*qin(0)+(1-qin(0))*double(rand())/RAND_MAX)*(x_max-x_min)/nx; 00630 double y = y_min+(j+0.5*qin(1)+(1-qin(1))*double(rand())/RAND_MAX)*(y_max-y_min)/ny; 00631 double z = z_min+(k+0.5*qin(2)+(1-qin(2))*double(rand())/RAND_MAX)*(z_max-z_min)/nz; 00632 con.put (n,x,y,z); 00633 n++; 00634 } 00635 } 00636 } 00637 00638 fpoint x,y,z,px,py,pz; 00639 container *cp = & con; 00640 voropp_loop l1(cp); 00641 int q,s; 00642 voronoicell_neighbor c; 00643 s=l1.init(con.ax,con.bx,con.ay,con.by,con.az,con.bz,px,py,pz); 00644 00645 Array<Array <size_t> > ListBpairs(n); 00646 size_t IIndex = Particles.Size(); 00647 do 00648 { 00649 for(q=0;q<con.co[s];q++) 00650 { 00651 x=con.p[s][con.sz*q]+px;y=con.p[s][con.sz*q+1]+py;z=con.p[s][con.sz*q+2]+pz; 00652 if(x>con.ax&&x<con.bx&&y>con.ay&&y<con.by&&z>con.az&&z<con.bz) 00653 { 00654 if(con.compute_cell(c,l1.ip,l1.jp,l1.kp,s,q,x,y,z)) 00655 { 00656 00657 if (rand()<fraction*RAND_MAX) 00658 { 00659 AddVoroCell(Tag,c,R,rho,true,iVec3_t(nx,ny,nz)); 00660 Vec3_t trans(x/nx,y/ny,z/nz); 00661 Particle * P = Particles[Particles.Size()-1]; 00662 P->Translate(trans); 00663 } 00664 } 00665 } 00666 } 00667 } while((s=l1.inc(px,py,pz))!=-1); 00668 00669 // info 00670 printf("%s Num of particles = %zd%s\n",TERM_CLR2,Particles.Size(),TERM_RST); 00671 00672 00673 if (Cohesion) 00674 { 00675 //if (fraction<1.0) throw new Fatal("Domain::AddVoroPack: With the Cohesion all particles should be considered, plese change the fraction to 1.0"); 00676 00677 // define some tolerance for comparissions 00678 double tol1 = 1.0e-8; 00679 double tol2 = 1.0e-3; 00680 for (size_t i=IIndex;i<Particles.Size()-1;i++) 00681 { 00682 Particle * P1 = Particles[i]; 00683 for (size_t j=i+1;j<Particles.Size();j++) 00684 { 00685 Particle * P2 = Particles[j]; 00686 if (Distance(P1->x,P2->x)<P1->Dmax+P2->Dmax) 00687 { 00688 for (size_t k=0;k<P1->Faces.Size();k++) 00689 { 00690 Face * F1 = P1->Faces[k]; 00691 Vec3_t n1,c1; 00692 F1->Normal (n1); 00693 F1->Centroid(c1); 00694 bool found = false; 00695 for (size_t l=0;l<P2->Faces.Size();l++) 00696 { 00697 Face * F2 = P2->Faces[l]; 00698 Vec3_t n2,c2; 00699 F2->Normal (n2); 00700 F2->Centroid(c2); 00701 Vec3_t n = 0.5*(n1-n2); 00702 n/=norm(n); 00703 if ((fabs(dot(n1,n2)+1.0)<tol1) 00704 &&(fabs(Distance(c1,*F2)-2*R)<tol2) 00705 &&(fabs(Distance(c2,*F1)-2*R)<tol2)) 00706 { 00707 BInteractons.Push(new BInteracton(P1,P2,k,l)); 00708 found = true; 00709 break; 00710 } 00711 } 00712 if (found) break; 00713 } 00714 } 00715 } 00716 } 00717 } 00718 } 00719 00720 // Single particle addition 00721 00722 inline void Domain::AddSphere (int Tag,Vec3_t const & X, double R, double rho) 00723 { 00724 // vertices 00725 Array<Vec3_t> V(1); 00726 V[0] = X; 00727 00728 // edges 00729 Array<Array <int> > E(0); // no edges 00730 00731 // faces 00732 Array<Array <int> > F(0); // no faces 00733 00734 // add particle 00735 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 00736 } 00737 00738 inline void Domain::AddCube (int Tag, Vec3_t const & X, double R, double L, double rho, double Angle, Vec3_t * Axis) 00739 { 00740 // vertices 00741 Array<Vec3_t> V(8); 00742 double l = L/2.0; 00743 V[0] = -l, -l, -l; 00744 V[1] = l, -l, -l; 00745 V[2] = l, l, -l; 00746 V[3] = -l, l, -l; 00747 V[4] = -l, -l, l; 00748 V[5] = l, -l, l; 00749 V[6] = l, l, l; 00750 V[7] = -l, l, l; 00751 00752 // edges 00753 Array<Array <int> > E(12); 00754 for (size_t i=0; i<12; ++i) E[i].Resize(2); 00755 E[ 0] = 0, 1; 00756 E[ 1] = 1, 2; 00757 E[ 2] = 2, 3; 00758 E[ 3] = 3, 0; 00759 E[ 4] = 4, 5; 00760 E[ 5] = 5, 6; 00761 E[ 6] = 6, 7; 00762 E[ 7] = 7, 4; 00763 E[ 8] = 0, 4; 00764 E[ 9] = 1, 5; 00765 E[10] = 2, 6; 00766 E[11] = 3, 7; 00767 00768 // faces 00769 Array<Array <int> > F(6); 00770 for (size_t i=0; i<6; i++) F[i].Resize(4); 00771 F[0] = 4, 7, 3, 0; 00772 F[1] = 1, 2, 6, 5; 00773 F[2] = 0, 1, 5, 4; 00774 F[3] = 2, 3, 7, 6; 00775 F[4] = 0, 3, 2, 1; 00776 F[5] = 4, 5, 6, 7; 00777 00778 // calculate the rotation 00779 bool ThereisanAxis = true; 00780 if (Axis==NULL) 00781 { 00782 Angle = (1.0*rand())/RAND_MAX*2*M_PI; 00783 Axis = new Vec3_t((1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX); 00784 ThereisanAxis = false; 00785 } 00786 Quaternion_t q; 00787 NormalizeRotation (Angle,(*Axis),q); 00788 for (size_t i=0; i<V.Size(); i++) 00789 { 00790 Vec3_t t; 00791 Rotation (V[i],q,t); 00792 V[i] = t+X; 00793 } 00794 00795 // add particle 00796 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 00797 00798 // clean up 00799 if (!ThereisanAxis) delete Axis; 00800 q(0) = 1.0; 00801 q(1) = 0.0; 00802 q(2) = 0.0; 00803 q(3) = 0.0; 00804 q = q/norm(q); 00805 00806 Particles[Particles.Size()-1]->Q = q; 00807 Particles[Particles.Size()-1]->Props.V = L*L*L; 00808 Particles[Particles.Size()-1]->Props.m = rho*L*L*L; 00809 Particles[Particles.Size()-1]->I = L*L, L*L, L*L; 00810 Particles[Particles.Size()-1]->I *= Particles[Particles.Size()-1]->Props.m/6.0; 00811 Particles[Particles.Size()-1]->x = X; 00812 Particles[Particles.Size()-1]->Ekin = 0.0; 00813 Particles[Particles.Size()-1]->Erot = 0.0; 00814 Particles[Particles.Size()-1]->Dmax = sqrt(3.0*L*L/4.0)+R; 00815 Particles[Particles.Size()-1]->PropsReady = true; 00816 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 00817 00818 00819 } 00820 00821 inline void Domain::AddTetra (int Tag, Vec3_t const & X, double R, double L, double rho, double Angle, Vec3_t * Axis) 00822 { 00823 // vertices 00824 double sq8 = sqrt(8.0); 00825 Array<Vec3_t> V(4); 00826 V[0] = L/sq8, L/sq8, L/sq8; 00827 V[1] = -L/sq8, -L/sq8, L/sq8; 00828 V[2] = -L/sq8, L/sq8,-L/sq8; 00829 V[3] = L/sq8, -L/sq8,-L/sq8; 00830 00831 // edges 00832 Array<Array <int> > E(6); 00833 for (size_t i=0; i<6; ++i) E[i].Resize(2); 00834 E[0] = 0, 1; 00835 E[1] = 1, 2; 00836 E[2] = 2, 0; 00837 E[3] = 0, 3; 00838 E[4] = 1, 3; 00839 E[5] = 2, 3; 00840 00841 // face 00842 Array<Array <int> > F; 00843 F.Resize(4); 00844 for (size_t i=0; i<4; ++i) F[i].Resize(3); 00845 F[0] = 0, 3, 2; 00846 F[1] = 0, 1, 3; 00847 F[2] = 0, 2, 1; 00848 F[3] = 1, 2, 3; 00849 00850 // calculate the rotation 00851 bool ThereisanAxis = true; 00852 if (Axis==NULL) 00853 { 00854 Angle = (1.0*rand())/RAND_MAX*2*M_PI; 00855 Axis = new Vec3_t((1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX); 00856 ThereisanAxis = false; 00857 } 00858 Quaternion_t q; 00859 NormalizeRotation (Angle,(*Axis),q); 00860 for (size_t i=0; i<V.Size(); i++) 00861 { 00862 Vec3_t t; 00863 Rotation (V[i],q,t); 00864 V[i] = t+X; 00865 } 00866 00867 // add particle 00868 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 00869 00870 // clean up 00871 if (!ThereisanAxis) delete Axis; 00872 q(0) = 1.0; 00873 q(1) = 0.0; 00874 q(2) = 0.0; 00875 q(3) = 0.0; 00876 q = q/norm(q); 00877 00878 Particles[Particles.Size()-1]->Q = q; 00879 Particles[Particles.Size()-1]->Props.V = sqrt(2.0)*L*L*L/12.0; 00880 Particles[Particles.Size()-1]->Props.m = rho*sqrt(2.0)*L*L*L/12.0; 00881 Particles[Particles.Size()-1]->I = L*L, L*L, L*L; 00882 Particles[Particles.Size()-1]->I *= Particles[Particles.Size()-1]->Props.m/20.0; 00883 Particles[Particles.Size()-1]->x = X; 00884 Particles[Particles.Size()-1]->Ekin = 0.0; 00885 Particles[Particles.Size()-1]->Erot = 0.0; 00886 Particles[Particles.Size()-1]->Dmax = sqrt(3.0*L*L/8.0)+R; 00887 Particles[Particles.Size()-1]->PropsReady = true; 00888 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 00889 } 00890 00891 inline void Domain::AddDrill (int Tag, const Vec3_t & X, double R, double Lt, double Ll, double rho) 00892 { 00893 Array<Vec3_t> V(9); 00894 V[0] = Lt/2.0, Lt/2.0, Ll/2.0; 00895 V[1] = -Lt/2.0, Lt/2.0, Ll/2.0; 00896 V[2] = -Lt/2.0, -Lt/2.0, Ll/2.0; 00897 V[3] = Lt/2.0, -Lt/2.0, Ll/2.0; 00898 V[4] = Lt/2.0, Lt/2.0, 0.0; 00899 V[5] = -Lt/2.0, Lt/2.0, 0.0; 00900 V[6] = -Lt/2.0, -Lt/2.0, 0.0; 00901 V[7] = Lt/2.0, -Lt/2.0, 0.0; 00902 V[8] = 0.0, 0.0,-Ll/2.0; 00903 00904 Array<Array <int> > E(12); 00905 for (size_t i=0; i<12; ++i) E[i].Resize(2); 00906 for (size_t i=0;i<4;i++) 00907 { 00908 E[i] = i , (i+1)%4 ; 00909 E[i+4] = i+4, (i+5)%4 + 4; 00910 E[i+8] = i+4, 8; 00911 } 00912 00913 Array<Array <int> > F; 00914 F.Resize(9); 00915 F[0].Resize(4); 00916 F[0] = 0, 1, 2, 3; 00917 F[1].Resize(4); 00918 F[1] = 0, 4, 5, 1; 00919 F[2].Resize(4); 00920 F[2] = 1, 5, 6, 2; 00921 F[3].Resize(4); 00922 F[3] = 2, 6, 7, 3; 00923 F[4].Resize(4); 00924 F[4] = 3, 7, 4, 0; 00925 F[5].Resize(3); 00926 F[5] = 4, 8, 5; 00927 F[6].Resize(3); 00928 F[6] = 5, 8, 6; 00929 F[7].Resize(3); 00930 F[7] = 6, 8, 7; 00931 F[8].Resize(3); 00932 F[8] = 7, 8, 4; 00933 00934 00935 double vol; // volume of the polyhedron 00936 Vec3_t CM; // Center of mass of the polyhedron 00937 Mat3_t It; // Inertia tensor of the polyhedron 00938 PolyhedraMP(V,F,vol,CM,It); 00939 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 00940 Particles[Particles.Size()-1]->x = CM; 00941 Particles[Particles.Size()-1]->Props.V = vol; 00942 Particles[Particles.Size()-1]->Props.m = vol*rho; 00943 Vec3_t I; 00944 Quaternion_t Q; 00945 Vec3_t xp,yp,zp; 00946 Eig(It,I,xp,yp,zp); 00947 CheckDestroGiro(xp,yp,zp); 00948 I *= rho; 00949 Q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2)); 00950 Q(1) = (yp(2)-zp(1))/(4*Q(0)); 00951 Q(2) = (zp(0)-xp(2))/(4*Q(0)); 00952 Q(3) = (xp(1)-yp(0))/(4*Q(0)); 00953 Q = Q/norm(Q); 00954 Particles[Particles.Size()-1]->I = I; 00955 Particles[Particles.Size()-1]->Q = Q; 00956 double Dmax = Distance(CM,V[0])+R; 00957 for (size_t i=1; i<V.Size(); ++i) 00958 { 00959 if (Distance(CM,V[i])+R > Dmax) Dmax = Distance(CM,V[i])+R; 00960 } 00961 Particles[Particles.Size()-1]->Ekin = 0.0; 00962 Particles[Particles.Size()-1]->Erot = 0.0; 00963 Particles[Particles.Size()-1]->Dmax = Dmax; 00964 Particles[Particles.Size()-1]->PropsReady = true; 00965 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 00966 00967 Vec3_t Y = X; 00968 Particles[Particles.Size()-1]->Translate(Y); 00969 00970 } 00971 00972 inline void Domain::AddRice (int Tag, const Vec3_t & X, double R, double L, double rho, double Angle, Vec3_t * Axis) 00973 { 00974 // vertices 00975 Array<Vec3_t> V(2); 00976 V[0] = 0.0, 0.0, L/2; 00977 V[1] = 0.0, 0.0, -L/2; 00978 00979 // edges 00980 Array<Array <int> > E(1); 00981 E[0].Resize(2); 00982 E[0] = 0, 1; 00983 00984 // faces 00985 Array<Array <int> > F(0); // no faces 00986 00987 // calculate the rotation 00988 bool ThereisanAxis = true; 00989 if (Axis==NULL) 00990 { 00991 Angle = (1.0*rand())/RAND_MAX*2*M_PI; 00992 Axis = new Vec3_t((1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX); 00993 ThereisanAxis = false; 00994 } 00995 Quaternion_t q; 00996 NormalizeRotation (Angle,(*Axis),q); 00997 for (size_t i=0; i<V.Size(); i++) 00998 { 00999 Vec3_t t; 01000 Rotation (V[i],q,t); 01001 V[i] = t+X; 01002 } 01003 01004 // add particle 01005 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 01006 01007 Particles[Particles.Size()-1]->Q = q; 01008 Particles[Particles.Size()-1]->Props.V = (4.0/3.0)*M_PI*R*R*R + M_PI*L*R*R; 01009 Particles[Particles.Size()-1]->Props.m = rho*Particles[Particles.Size()-1]->Props.V; 01010 Particles[Particles.Size()-1]->I = (1.0/3.0)*M_PI*rho*R*R*R*L*L + (1.0/12.0)*M_PI*rho*R*R*L*L*L + (3.0/4.0)*M_PI*rho*R*R*R*R*L + (8.0/15.0)*M_PI*rho*R*R*R*R*R, 01011 (1.0/3.0)*M_PI*rho*R*R*R*L*L + (1.0/12.0)*M_PI*rho*R*R*L*L*L + (3.0/4.0)*M_PI*rho*R*R*R*R*L + (8.0/15.0)*M_PI*rho*R*R*R*R*R, 01012 0.5*M_PI*rho*R*R*R*R*L + (8.0/15.0)*M_PI*rho*R*R*R*R*R; 01013 Particles[Particles.Size()-1]->x = X; 01014 Particles[Particles.Size()-1]->Ekin = 0.0; 01015 Particles[Particles.Size()-1]->Erot = 0.0; 01016 Particles[Particles.Size()-1]->Dmax = 0.5*L + R; 01017 Particles[Particles.Size()-1]->PropsReady = true; 01018 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 01019 01020 // clean up 01021 if (!ThereisanAxis) delete Axis; 01022 } 01023 01024 inline void Domain::AddPlane (int Tag, const Vec3_t & X, double R, double Lx, double Ly, double rho, double Angle, Vec3_t * Axis) 01025 { 01026 // vertices 01027 Array<Vec3_t> V(4); 01028 double lx = Lx/2.0, ly = Ly/2.0; 01029 V[0] = -lx, -ly, 0.0; 01030 V[1] = lx, -ly, 0.0; 01031 V[2] = lx, ly, 0.0; 01032 V[3] = -lx, ly, 0.0; 01033 01034 // edges 01035 Array<Array <int> > E(4); 01036 for (size_t i=0; i<4; ++i) E[i].Resize(2); 01037 E[ 0] = 0, 1; 01038 E[ 1] = 1, 2; 01039 E[ 2] = 2, 3; 01040 E[ 3] = 3, 0; 01041 01042 // faces 01043 Array<Array <int> > F(1); 01044 F[0].Resize(4); 01045 F[0] = 0, 3, 2, 1; 01046 01047 bool ThereisanAxis = true; 01048 if (Axis==NULL) 01049 { 01050 Angle = 0.; 01051 Axis = new Vec3_t((1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX, (1.0*rand())/RAND_MAX); 01052 ThereisanAxis = false; 01053 } 01054 Quaternion_t q; 01055 NormalizeRotation (Angle,(*Axis),q); 01056 for (size_t i=0; i<V.Size(); i++) 01057 { 01058 Vec3_t t; 01059 Rotation (V[i],q,t); 01060 V[i] = t+X; 01061 } 01062 01063 // add particle 01064 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 01065 Particles[Particles.Size()-1]->Q = q; 01066 Particles[Particles.Size()-1]->Props.V = Lx*Ly*2*R; 01067 Particles[Particles.Size()-1]->Props.m = rho*Lx*Ly*2*R; 01068 Particles[Particles.Size()-1]->I = (1.0/12.0)*(Ly*Ly+4*R*R),(1.0/12.0)*(Lx*Lx+4*R*R),(1.0/12.0)*(Lx*Lx+Ly*Ly); 01069 Particles[Particles.Size()-1]->I *= Particles[Particles.Size()-1]->Props.m; 01070 Particles[Particles.Size()-1]->x = X; 01071 Particles[Particles.Size()-1]->Ekin = 0.0; 01072 Particles[Particles.Size()-1]->Erot = 0.0; 01073 Particles[Particles.Size()-1]->Dmax = sqrt(Lx*Lx+Ly*Ly)+R; 01074 Particles[Particles.Size()-1]->PropsReady = true; 01075 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 01076 // clean up 01077 if (!ThereisanAxis) delete Axis; 01078 } 01079 01080 inline void Domain::AddVoroCell (int Tag, voronoicell_neighbor & VC, double R, double rho, bool Erode, iVec3_t nv) 01081 { 01082 Array<Vec3_t> V(VC.p); 01083 Array<Array <int> > E; 01084 Array<int> Eaux(2); 01085 for(int i=0;i<VC.p;i++) 01086 { 01087 V[i] = Vec3_t(0.5*VC.pts[3*i]/nv(0),0.5*VC.pts[3*i+1]/nv(1),0.5*VC.pts[3*i+2]/nv(2)); 01088 for(int j=0;j<VC.nu[i];j++) 01089 { 01090 int k=VC.ed[i][j]; 01091 if (VC.ed[i][j]<i) 01092 { 01093 Eaux[0] = i; 01094 Eaux[1] = k; 01095 E.Push(Eaux); 01096 } 01097 } 01098 } 01099 Array<Array <int> > F; 01100 Array<int> Faux; 01101 for(int i=0;i<VC.p;i++) 01102 { 01103 for(int j=0;j<VC.nu[i];j++) 01104 { 01105 int k=VC.ed[i][j]; 01106 if (k>=0) 01107 { 01108 Faux.Push(i); 01109 VC.ed[i][j]=-1-k; 01110 int l=VC.cycle_up(VC.ed[i][VC.nu[i]+j],k); 01111 do 01112 { 01113 Faux.Push(k); 01114 int m=VC.ed[k][l]; 01115 VC.ed[k][l]=-1-m; 01116 l=VC.cycle_up(VC.ed[k][VC.nu[k]+l],m); 01117 k=m; 01118 } while (k!=i); 01119 Array<int> Faux2(Faux.Size()); 01120 for (size_t l = 0; l < Faux.Size();l++) 01121 { 01122 Faux2[l] = Faux[Faux.Size()-1-l]; 01123 } 01124 01125 F.Push(Faux2); 01126 Faux.Clear(); 01127 Faux2.Clear(); 01128 } 01129 } 01130 } 01131 VC.reset_edges(); 01132 double vol; // volume of the polyhedron 01133 Vec3_t CM; // Center of mass of the polyhedron 01134 Mat3_t It; // Inertia tensor of the polyhedron 01135 PolyhedraMP(V,F,vol,CM,It); 01136 if (Erode) Erosion(V,E,F,R); 01137 // add particle 01138 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 01139 Particles[Particles.Size()-1]->x = CM; 01140 Particles[Particles.Size()-1]->Props.V = vol; 01141 Particles[Particles.Size()-1]->Props.m = vol*rho; 01142 Vec3_t I; 01143 Quaternion_t Q; 01144 Vec3_t xp,yp,zp; 01145 Eig(It,I,xp,yp,zp); 01146 CheckDestroGiro(xp,yp,zp); 01147 I *= rho; 01148 Q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2)); 01149 Q(1) = (yp(2)-zp(1))/(4*Q(0)); 01150 Q(2) = (zp(0)-xp(2))/(4*Q(0)); 01151 Q(3) = (xp(1)-yp(0))/(4*Q(0)); 01152 Q = Q/norm(Q); 01153 Particles[Particles.Size()-1]->I = I; 01154 Particles[Particles.Size()-1]->Q = Q; 01155 double Dmax = Distance(CM,V[0])+R; 01156 for (size_t i=1; i<V.Size(); ++i) 01157 { 01158 if (Distance(CM,V[i])+R > Dmax) Dmax = Distance(CM,V[i])+R; 01159 } 01160 Particles[Particles.Size()-1]->Ekin = 0.0; 01161 Particles[Particles.Size()-1]->Erot = 0.0; 01162 Particles[Particles.Size()-1]->Dmax = Dmax; 01163 Particles[Particles.Size()-1]->PropsReady = true; 01164 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 01165 } 01166 01167 inline void Domain::AddTorus (int Tag, Vec3_t const & X, Vec3_t const & N, double Rmax, double R, double rho) 01168 { 01169 // Normalize normal vector 01170 Vec3_t n = N/norm(N); 01171 01172 // Create the 2 vertices that define the torus 01173 Vec3_t P1 = OrthoSys::e0 - dot(OrthoSys::e0,n)*n; 01174 if (norm(P1)<1.0e-12) P1 = OrthoSys::e1 - dot(OrthoSys::e1,n)*n; 01175 P1 /= norm(P1); 01176 Vec3_t P2 = cross(n,P1); 01177 P2 /= norm(P2); 01178 Array<Vec3_t > V(2); 01179 V[0] = X + Rmax*P1; 01180 V[1] = X + Rmax*P2; 01181 01182 // the torus has no edges or faces 01183 Array<Array <int> > E(0); 01184 Array<Array <int> > F(0); 01185 01186 //Add the particle just with two vertices 01187 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 01188 01189 //Input all the mass properties 01190 Vec3_t xp,yp,zp; 01191 xp = P1; 01192 zp = P2; 01193 yp = cross(zp,xp); 01194 CheckDestroGiro(xp,yp,zp); 01195 Quaternion_t q; 01196 q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2)); 01197 q(1) = (yp(2)-zp(1))/(4*q(0)); 01198 q(2) = (zp(0)-xp(2))/(4*q(0)); 01199 q(3) = (xp(1)-yp(0))/(4*q(0)); 01200 q = q/norm(q); 01201 01202 Particles[Particles.Size()-1]->Q = q; 01203 Particles[Particles.Size()-1]->Props.V = 2*M_PI*M_PI*R*R*Rmax; 01204 Particles[Particles.Size()-1]->Props.m = rho*2*M_PI*M_PI*R*R*Rmax; 01205 Particles[Particles.Size()-1]->I = 0.125*(4*Rmax*Rmax + 5*R*R), (Rmax*Rmax+0.75*R*R), 0.125*(4*Rmax*Rmax + 5*R*R); 01206 Particles[Particles.Size()-1]->I *= Particles[Particles.Size()-1]->Props.m; 01207 Particles[Particles.Size()-1]->x = X; 01208 Particles[Particles.Size()-1]->Ekin = 0.0; 01209 Particles[Particles.Size()-1]->Erot = 0.0; 01210 Particles[Particles.Size()-1]->Dmax = Rmax + R; 01211 Particles[Particles.Size()-1]->PropsReady = true; 01212 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 01213 01214 Particles[Particles.Size()-1]->Tori.Push(new Torus(&Particles[Particles.Size()-1]->x,Particles[Particles.Size()-1]->Verts[0],Particles[Particles.Size()-1]->Verts[1])); 01215 01216 } 01217 01218 inline void Domain::AddCylinder (int Tag, Vec3_t const & X0, double R0, Vec3_t const & X1, double R1, double R, double rho) 01219 { 01220 // Define all key variables of the cylinder 01221 Vec3_t n = X1 - X0; 01222 n /= norm(n); 01223 Vec3_t P1 = OrthoSys::e0 - dot(OrthoSys::e0,n)*n; 01224 if (norm(P1)<1.0e-12) P1 = OrthoSys::e1 - dot(OrthoSys::e1,n)*n; 01225 P1 /= norm(P1); 01226 Vec3_t P2 = cross(n,P1); 01227 P2 /= norm(P2); 01228 01229 //The cylinder is defined by 6 vertices 01230 Array<Vec3_t > V(6); 01231 V[0] = X0 + R0*P1; 01232 V[1] = X0 + R0*P2; 01233 V[2] = X0 - R0*P1; 01234 V[3] = X1 + R1*P1; 01235 V[4] = X1 + R1*P2; 01236 V[5] = X1 - R1*P1; 01237 01238 // It has no edges or faces 01239 Array<Array <int> > E(0); 01240 Array<Array <int> > F(0); 01241 01242 //Add the particle just with the vertices 01243 Particles.Push (new Particle(Tag,V,E,F,OrthoSys::O,OrthoSys::O,R,rho)); 01244 01245 //Input all the mass properties 01246 Vec3_t xp,yp,zp; 01247 xp = P1; 01248 zp = P2; 01249 yp = cross(zp,xp); 01250 CheckDestroGiro(xp,yp,zp); 01251 Quaternion_t q; 01252 q(0) = 0.5*sqrt(1+xp(0)+yp(1)+zp(2)); 01253 q(1) = (yp(2)-zp(1))/(4*q(0)); 01254 q(2) = (zp(0)-xp(2))/(4*q(0)); 01255 q(3) = (xp(1)-yp(0))/(4*q(0)); 01256 q = q/norm(q); 01257 01258 Particles[Particles.Size()-1]->Q = q; 01259 Particles[Particles.Size()-1]->Props.V = M_PI*R0*R0*norm(X1-X0); 01260 Particles[Particles.Size()-1]->Props.m = rho*M_PI*R0*R0*norm(X1-X0); 01261 Particles[Particles.Size()-1]->I = 1.0, 1.0, 1.0; 01262 Particles[Particles.Size()-1]->x = 0.5*(X0 + X1); 01263 Particles[Particles.Size()-1]->Ekin = 0.0; 01264 Particles[Particles.Size()-1]->Erot = 0.0; 01265 Particles[Particles.Size()-1]->Dmax = sqrt(0.25*dot(X1-X0,X1-X0)+max(R0,R1)*max(R0,R1)) + R; 01266 Particles[Particles.Size()-1]->PropsReady = true; 01267 Particles[Particles.Size()-1]->Index = Particles.Size()-1; 01268 Particles[Particles.Size()-1]->Tori.Push (new Torus(&X0,Particles[Particles.Size()-1]->Verts[0],Particles[Particles.Size()-1]->Verts[1])); 01269 Particles[Particles.Size()-1]->Tori.Push (new Torus(&X1,Particles[Particles.Size()-1]->Verts[3],Particles[Particles.Size()-1]->Verts[4])); 01270 Particles[Particles.Size()-1]->Cylinders.Push(new Cylinder(Particles[Particles.Size()-1]->Tori[0],Particles[Particles.Size()-1]->Tori[1],Particles[Particles.Size()-1]->Verts[2],Particles[Particles.Size()-1]->Verts[5])); 01271 } 01272 01273 // Methods 01274 01275 inline void Domain::SetProps (Dict & D) 01276 { 01277 for (size_t i =0 ; i<Particles.Size(); i++) 01278 { 01279 for (size_t j=0; j<D.Keys.Size(); ++j) 01280 { 01281 int tag = D.Keys[j]; 01282 if (tag==Particles[i]->Tag) 01283 { 01284 SDPair const & p = D(tag); 01285 if (p.HasKey("Gn")) 01286 { 01287 Particles[i]->Props.Gn = p("Gn"); 01288 } 01289 if (p.HasKey("Gt")) 01290 { 01291 Particles[i]->Props.Gt = p("Gt"); 01292 } 01293 if (p.HasKey("Kn")) 01294 { 01295 Particles[i]->Props.Kn = p("Kn"); 01296 } 01297 if (p.HasKey("Kt")) 01298 { 01299 Particles[i]->Props.Kt = p("Kt"); 01300 } 01301 if (p.HasKey("Bn")) 01302 { 01303 Particles[i]->Props.Bn = p("Bn"); 01304 } 01305 if (p.HasKey("Bt")) 01306 { 01307 Particles[i]->Props.Bt = p("Bt"); 01308 } 01309 if (p.HasKey("Bm")) 01310 { 01311 Particles[i]->Props.Bm = p("Bm"); 01312 } 01313 if (p.HasKey("Mu")) 01314 { 01315 Particles[i]->Props.Mu = p("Mu"); 01316 } 01317 if (p.HasKey("Eps")) 01318 { 01319 Particles[i]->Props.eps = p("Eps"); 01320 } 01321 if (p.HasKey("Beta")) 01322 { 01323 Particles[i]->Props.Beta = p("Beta"); 01324 } 01325 if (p.HasKey("Eta")) 01326 { 01327 Particles[i]->Props.Eta = p("Eta"); 01328 } 01329 } 01330 } 01331 } 01332 for (size_t i=0; i<BInteractons.Size(); i++) 01333 { 01334 BInteractons[i]->UpdateParameters(); 01335 } 01336 for (size_t i=0; i<CInteractons.Size(); i++) 01337 { 01338 CInteractons[i]->UpdateParameters(); 01339 } 01340 } 01341 01342 inline void Domain::Initialize (double dt) 01343 { 01344 if (!Initialized) 01345 { 01346 // initialize all particles 01347 for (size_t i=0; i<Particles.Size(); i++) 01348 { 01349 Particles[i]->Initialize(i); 01350 Particles[i]->InitializeVelocity(dt); 01351 } 01352 //Initializing the energies 01353 Evis = 0.0; 01354 Efric = 0.0; 01355 Wext = 0.0; 01356 01357 // initialize 01358 ResetInteractons(); 01359 // info 01360 Util::Stopwatch stopwatch; 01361 printf("\n%s--- Initializing particles ------------------------------------------------------%s\n",TERM_CLR1,TERM_RST); 01362 // set flag 01363 Initialized = true; 01364 01365 // info 01366 double Ekin, Epot, Etot; 01367 Etot = CalcEnergy (Ekin, Epot); 01368 printf("%s Kinematic energy = %g%s\n",TERM_CLR4, Ekin, TERM_RST); 01369 printf("%s Potential energy = %g%s\n",TERM_CLR4, Epot, TERM_RST); 01370 printf("%s Total energy = %g%s\n",TERM_CLR2, Etot, TERM_RST); 01371 } 01372 else 01373 { 01374 for (size_t i=0; i<Particles.Size(); i++) 01375 { 01376 if (Particles[i]->vxf) Particles[i]->xb(0) = Particles[i]->x(0) - Particles[i]->v(0)*dt; 01377 if (Particles[i]->vyf) Particles[i]->xb(1) = Particles[i]->x(1) - Particles[i]->v(1)*dt; 01378 if (Particles[i]->vzf) Particles[i]->xb(2) = Particles[i]->x(2) - Particles[i]->v(2)*dt; 01379 } 01380 } 01381 01382 } 01383 01384 inline void Domain::Solve (double tf, double dt, double dtOut, ptFun_t ptSetup, ptFun_t ptReport, char const * TheFileKey, bool RenderVideo, size_t Nproc) 01385 { 01386 // Assigning some domain particles especifically to the output 01387 FileKey.Printf("%s",TheFileKey); 01388 idx_out = 0; 01389 01390 // initialize particles 01391 Initialize (dt); 01392 01393 // set the displacement of the particles to zero (for the Halo) 01394 ResetDisplacements(); 01395 01396 // build the map of possible contacts (for the Halo) 01397 ResetContacts(); 01398 01399 // calc the total volume of particles (solids) 01400 Vs = 0.0; 01401 Ms = 0.0; 01402 for (size_t i=0; i<Particles.Size(); i++) 01403 { 01404 if (Particles[i]->IsFree()) 01405 { 01406 Vs += Particles[i]->Props.V; 01407 Ms += Particles[i]->Props.m; 01408 } 01409 } 01410 01411 // info 01412 Util::Stopwatch stopwatch; 01413 printf("\n%s--- Solving ---------------------------------------------------------------------%s\n",TERM_CLR1,TERM_RST); 01414 01415 // solve 01416 double t0 = Time; // initial time 01417 double tout = t0+dtOut; // time position for output 01418 01419 // report 01420 Finished = false; 01421 if (ptReport!=NULL) (*ptReport) ((*this), UserData); 01422 01423 // string to output energy data, if user gives the FileKey 01424 std::ostringstream oss_energy; 01425 EnergyOutput (idx_out, oss_energy); 01426 01427 // run 01428 while (Time<tf) 01429 { 01430 01431 // initialize forces and torques 01432 for (size_t i=0; i<Particles.Size(); i++) 01433 { 01434 // set the force and torque to the fixed values 01435 Particles[i]->F = Particles[i]->Ff; 01436 Particles[i]->T = Particles[i]->Tf; 01437 for (size_t n=0;n<3;n++) 01438 { 01439 for (size_t m=0;m<3;m++) 01440 { 01441 Particles[i]->M(n,m)=0.0; 01442 Particles[i]->B(n,m)=0.0; 01443 } 01444 } 01445 01446 // initialize the coordination (number of contacts per particle) number 01447 Particles[i]->Cn = 0.0; 01448 01449 // external work added to the system by the fixed forces Ff 01450 Wext += dot(Particles[i]->Ff,Particles[i]->v)*dt; 01451 } 01452 01453 // calc contact forces: collision and bonding (cohesion) 01454 #ifdef USE_THREAD 01455 // declare a thread array for the interactons 01456 vector<std::thread> IT; 01457 for (size_t n=0;n<Nproc;n++) 01458 { 01459 IT.push_back(std::thread(std::bind(GlobalForce,&Interactons,dt,n,Nproc))); 01460 } 01461 for (size_t n=0;n<Nproc;n++) 01462 { 01463 IT[n].join(); 01464 } 01465 #else 01466 for (size_t i=0; i<Interactons.Size(); i++) 01467 { 01468 Interactons[i]->CalcForce (dt); 01469 } 01470 #endif 01471 01472 // calculate the collision energy 01473 for (size_t i=0; i<CInteractons.Size(); i++) 01474 { 01475 Evis += CInteractons[i]-> dEvis; 01476 Efric += CInteractons[i]-> dEfric; 01477 } 01478 01479 // tell the user function to update its data 01480 if (ptSetup!=NULL) (*ptSetup) ((*this), UserData); 01481 01482 01483 // move particles 01484 #ifdef USE_THREAD 01485 // declare a thread array for the particles 01486 vector<std::thread> MT; 01487 Array<double> TMD; 01488 TMD.Resize(Nproc); 01489 for (size_t n=0;n<Nproc;n++) 01490 { 01491 //call the function move for each thread 01492 MT.push_back(std::thread(std::bind(GlobalMove,&Particles,dt,TMD[n],n,Nproc))); 01493 } 01494 for (size_t n=0;n<Nproc;n++) 01495 { 01496 // wait for all the threads to finish 01497 MT[n].join(); 01498 } 01499 #else 01500 for (size_t i=0; i<Particles.Size(); i++) 01501 { 01502 Particles[i]->Translate (dt); 01503 Particles[i]->Rotate (dt); 01504 } 01505 #endif 01506 01507 // output 01508 if (Time>=tout) 01509 { 01510 idx_out++; 01511 if (BInteractons.Size()>0) Clusters(); 01512 if (ptReport!=NULL) (*ptReport) ((*this), UserData); 01513 if (TheFileKey!=NULL) 01514 { 01515 String fn; 01516 fn.Printf ("%s_%04d", TheFileKey, idx_out); 01517 if(RenderVideo) WritePOV (fn.CStr()); 01518 EnergyOutput (idx_out, oss_energy); 01519 } 01520 tout += dtOut; 01521 } 01522 01523 01524 // next time position 01525 Time += dt; 01526 01527 01528 #ifdef USE_THREAD 01529 double maxdis = TMD.TheMax(); 01530 //std::cout << TMD.TheMax() << std::endl; 01531 maxdis = MaxDisplacement(); 01532 #else 01533 double maxdis = MaxDisplacement(); 01534 #endif 01535 // update the Halos 01536 if (maxdis>Alpha) 01537 { 01538 ResetDisplacements(); 01539 ResetContacts(); 01540 } 01541 01542 //if (fabs(Time-7.15013)<1.0e-5) WriteBPY("test"); 01543 } 01544 01545 // last output 01546 Finished = true; 01547 if (ptReport!=NULL) (*ptReport) ((*this), UserData); 01548 01549 // save energy data 01550 if (TheFileKey!=NULL) 01551 { 01552 String fn; 01553 fn.Printf("%s_energy.res",TheFileKey); 01554 std::ofstream fe(fn.CStr()); 01555 fe << oss_energy.str(); 01556 fe.close(); 01557 } 01558 01559 // info 01560 double Ekin, Epot, Etot; 01561 Etot = CalcEnergy (Ekin, Epot); 01562 printf("%s Kinematic energy = %g%s\n",TERM_CLR4, Ekin, TERM_RST); 01563 printf("%s Potential energy = %g%s\n",TERM_CLR4, Epot, TERM_RST); 01564 printf("%s Total energy = %g%s\n",TERM_CLR2, Etot, TERM_RST); 01565 } 01566 01567 inline void Domain::WritePOV (char const * FileKey) 01568 { 01569 String fn(FileKey); 01570 fn.append(".pov"); 01571 std::ofstream of(fn.CStr(), std::ios::out); 01572 POVHeader (of); 01573 POVSetCam (of, CamPos, OrthoSys::O); 01574 Array <String> Colors(10); 01575 Colors = "Gray","Blue","Yellow","Gold","Green","Blue","Orange","Salmon","Copper","Aquamarine"; 01576 for (size_t i=0; i<Particles.Size(); i++) 01577 { 01578 if (!Particles[i]->IsFree()) Particles[i]->Draw(of,"Col_Glass_Bluish"); 01579 else 01580 { 01581 bool found = false; 01582 for (size_t j=0;j<Listofclusters.Size();j++) 01583 { 01584 if (Listofclusters[j].Has(i)&&BInteractons.Size()>0) 01585 { 01586 Particles[i]->Draw(of,Colors[j%10].CStr()); 01587 found = true; 01588 break; 01589 } 01590 } 01591 if (!found) Particles[i]->Draw(of,"Red"); 01592 } 01593 } 01594 01595 } 01596 01597 inline void Domain::WriteBPY (char const * FileKey) 01598 { 01599 String fn(FileKey); 01600 fn.append(".bpy"); 01601 std::ofstream of(fn.CStr(), std::ios::out); 01602 BPYHeader(of); 01603 for (size_t i=0; i<Particles.Size(); i++) Particles[i]->Draw (of,"",true); 01604 of.close(); 01605 } 01606 01607 #ifdef USE_HDF5 01608 inline void Domain::Save (char const * FileKey) 01609 { 01610 01611 // Opening the file for writing 01612 String fn(FileKey); 01613 fn.append(".hdf5"); 01614 hid_t file_id; 01615 file_id = H5Fcreate(fn.CStr(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); 01616 01617 // Storing the number of particles in the domain 01618 int data[1]; 01619 data[0]=Particles.Size(); 01620 hsize_t dims[1]; 01621 dims[0]=1; 01622 H5LTmake_dataset_int(file_id,"/NP",1,dims,data); 01623 01624 for (size_t i=0; i<Particles.Size(); i++) 01625 { 01626 // Creating the string and the group for each particle 01627 hid_t group_id; 01628 String par; 01629 par.Printf("/Particle_%08d",i); 01630 group_id = H5Gcreate(file_id, par.CStr(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 01631 01632 01633 // Storing some scalar variables 01634 double dat[1]; 01635 dat[0] = Particles[i]->Props.R; 01636 H5LTmake_dataset_double(group_id,"SR",1,dims,dat); 01637 dat[0] = Particles[i]->Props.rho; 01638 H5LTmake_dataset_double(group_id,"Rho",1,dims,dat); 01639 dat[0] = Particles[i]->Props.m; 01640 H5LTmake_dataset_double(group_id,"m",1,dims,dat); 01641 dat[0] = Particles[i]->Props.V; 01642 H5LTmake_dataset_double(group_id,"V",1,dims,dat); 01643 dat[0] = Particles[i]->Diam; 01644 H5LTmake_dataset_double(group_id,"Diam",1,dims,dat); 01645 dat[0] = Particles[i]->Dmax; 01646 H5LTmake_dataset_double(group_id,"Dmax",1,dims,dat); 01647 int datint[1]; 01648 datint[0] = Particles[i]->Index; 01649 H5LTmake_dataset_int(group_id,"Index",1,dims,datint); 01650 01651 01652 int tag[1]; 01653 tag[0] = Particles[i]->Tag; 01654 H5LTmake_dataset_int(group_id,"Tag",1,dims,tag); 01655 01656 // Storing vectorial variables 01657 double cd[3]; 01658 hsize_t dd[1]; 01659 dd[0] = 3; 01660 01661 cd[0]=Particles[i]->x(0); 01662 cd[1]=Particles[i]->x(1); 01663 cd[2]=Particles[i]->x(2); 01664 H5LTmake_dataset_double(group_id,"x",1,dd,cd); 01665 01666 cd[0]=Particles[i]->xb(0); 01667 cd[1]=Particles[i]->xb(1); 01668 cd[2]=Particles[i]->xb(2); 01669 H5LTmake_dataset_double(group_id,"xb",1,dd,cd); 01670 01671 cd[0]=Particles[i]->v(0); 01672 cd[1]=Particles[i]->v(1); 01673 cd[2]=Particles[i]->v(2); 01674 H5LTmake_dataset_double(group_id,"v",1,dd,cd); 01675 01676 cd[0]=Particles[i]->w(0); 01677 cd[1]=Particles[i]->w(1); 01678 cd[2]=Particles[i]->w(2); 01679 H5LTmake_dataset_double(group_id,"w",1,dd,cd); 01680 01681 cd[0]=Particles[i]->wb(0); 01682 cd[1]=Particles[i]->wb(1); 01683 cd[2]=Particles[i]->wb(2); 01684 H5LTmake_dataset_double(group_id,"wb",1,dd,cd); 01685 01686 cd[0]=Particles[i]->I(0); 01687 cd[1]=Particles[i]->I(1); 01688 cd[2]=Particles[i]->I(2); 01689 H5LTmake_dataset_double(group_id,"I",1,dd,cd); 01690 01691 double cq[4]; 01692 dd[0] = 4; 01693 cq[0]=Particles[i]->Q(0); 01694 cq[1]=Particles[i]->Q(1); 01695 cq[2]=Particles[i]->Q(2); 01696 cq[3]=Particles[i]->Q(3); 01697 H5LTmake_dataset_double(group_id,"Q",1,dd,cq); 01698 01699 01700 01701 01702 // Storing the number of vertices of each particle 01703 data[0] = Particles[i]->Verts.Size(); 01704 H5LTmake_dataset_int(group_id,"n_vertices",1,dims,data); 01705 hid_t gv_id; 01706 gv_id = H5Gcreate(group_id,"Verts", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 01707 01708 // Storing each vertex 01709 for (size_t j=0;j<Particles[i]->Verts.Size();j++) 01710 { 01711 String parv; 01712 parv.Printf("Verts_%08d",j); 01713 double cod[3]; 01714 cod[0]=(*Particles[i]->Verts[j])(0); 01715 cod[1]=(*Particles[i]->Verts[j])(1); 01716 cod[2]=(*Particles[i]->Verts[j])(2); 01717 hsize_t dim[1]; 01718 dim[0]=3; 01719 H5LTmake_dataset_double(gv_id,parv.CStr(),1,dim,cod); 01720 } 01721 01722 // Number of edges of the particle 01723 data[0] = Particles[i]->Edges.Size(); 01724 H5LTmake_dataset_int(group_id,"n_edges",1,dims,data); 01725 gv_id = H5Gcreate(group_id,"Edges", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 01726 01727 // Edges 01728 for (size_t j=0;j<Particles[i]->Edges.Size();j++) 01729 { 01730 String parv; 01731 parv.Printf("Edges_%08d",j); 01732 int co[2]; 01733 co[0] = Particles[i]->EdgeCon[j][0]; 01734 co[1] = Particles[i]->EdgeCon[j][1]; 01735 hsize_t dim[1]; 01736 dim[0] =2; 01737 H5LTmake_dataset_int(gv_id,parv.CStr(),1,dim,co); 01738 } 01739 01740 // Number of faces of the particle 01741 data[0] = Particles[i]->Faces.Size(); 01742 H5LTmake_dataset_int(group_id,"n_faces",1,dims,data); 01743 gv_id = H5Gcreate(group_id,"Faces", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); 01744 01745 // Faces 01746 for (size_t j=0;j<Particles[i]->Faces.Size();j++) 01747 { 01748 String parv; 01749 parv.Printf("Faces_%08d",j); 01750 int co[Particles[i]->FaceCon[j].Size()]; 01751 hsize_t dim[1]; 01752 dim[0]= Particles[i]->FaceCon[j].Size(); 01753 for (size_t k=0;k<Particles[i]->FaceCon[j].Size();k++) 01754 { 01755 co[k]=Particles[i]->FaceCon[j][k]; 01756 } 01757 H5LTmake_dataset_int(gv_id,parv.CStr(),1,dim,co); 01758 } 01759 01760 } 01761 01762 H5Fclose(file_id); 01763 01764 01765 } 01766 01767 inline void Domain::Load (char const * FileKey) 01768 { 01769 01770 // Opening the file for reading 01771 String fn(FileKey); 01772 fn.append(".hdf5"); 01773 hid_t file_id; 01774 file_id = H5Fopen(fn.CStr(), H5F_ACC_RDONLY, H5P_DEFAULT); 01775 01776 // Number of particles in the domain 01777 int data[1]; 01778 H5LTread_dataset_int(file_id,"/NP",data); 01779 size_t NP = data[0]; 01780 01781 // Loading the particles 01782 for (size_t i=0; i<NP; i++) 01783 { 01784 01785 // Creating the string and the group for each particle 01786 hid_t group_id; 01787 String par; 01788 par.Printf("/Particle_%08d",i); 01789 group_id = H5Gopen(file_id, par.CStr(),H5P_DEFAULT); 01790 01791 // Finding the particle's position for the domain decomposition 01792 double X[3]; 01793 H5LTread_dataset_double(group_id,"x",X); 01794 01795 01796 // Loading the Vertices 01797 H5LTread_dataset_int(group_id,"n_vertices",data); 01798 size_t nv = data[0]; 01799 hid_t gv_id; 01800 gv_id = H5Gopen(group_id,"Verts", H5P_DEFAULT); 01801 Array<Vec3_t> V; 01802 01803 for (size_t j=0;j<nv;j++) 01804 { 01805 String parv; 01806 parv.Printf("Verts_%08d",j); 01807 double cod[3]; 01808 H5LTread_dataset_double(gv_id,parv.CStr(),cod); 01809 V.Push(Vec3_t(cod[0],cod[1],cod[2])); 01810 } 01811 01812 // Loading the edges 01813 H5LTread_dataset_int(group_id,"n_edges",data); 01814 size_t ne = data[0]; 01815 gv_id = H5Gopen(group_id,"Edges", H5P_DEFAULT); 01816 Array<Array <int> > E; 01817 01818 for (size_t j=0;j<ne;j++) 01819 { 01820 String parv; 01821 parv.Printf("Edges_%08d",j); 01822 int cod[2]; 01823 H5LTread_dataset_int(gv_id,parv.CStr(),cod); 01824 Array<int> Ep(2); 01825 Ep[0]=cod[0]; 01826 Ep[1]=cod[1]; 01827 E.Push(Ep); 01828 } 01829 01830 // Loading the faces 01831 01832 // Number of faces of the particle 01833 H5LTread_dataset_int(group_id,"n_faces",data); 01834 size_t nf = data[0]; 01835 gv_id = H5Gopen(group_id,"Faces", H5P_DEFAULT); 01836 Array<Array <int> > F; 01837 01838 // Faces 01839 for (size_t j=0;j<nf;j++) 01840 { 01841 String parv; 01842 parv.Printf("Faces_%08d",j); 01843 hsize_t dim[1]; 01844 H5LTget_dataset_info(gv_id,parv.CStr(),dim,NULL,NULL); 01845 size_t ns = (size_t)dim[0]; 01846 int co[ns]; 01847 Array<int> Fp(ns); 01848 01849 H5LTread_dataset_int(gv_id,parv.CStr(),co); 01850 01851 for (size_t k=0;k<ns;k++) 01852 { 01853 Fp[k] = co[k]; 01854 } 01855 01856 F.Push(Fp); 01857 01858 } 01859 01860 Particles.Push (new Particle(-1,V,E,F,OrthoSys::O,OrthoSys::O,0.1,1.0)); 01861 01862 // Loading vectorial variables 01863 Particles[Particles.Size()-1]->x = Vec3_t(X[0],X[1],X[2]); 01864 double cd[3]; 01865 H5LTread_dataset_double(group_id,"xb",cd); 01866 Particles[Particles.Size()-1]->xb = Vec3_t(cd[0],cd[1],cd[2]); 01867 H5LTread_dataset_double(group_id,"v",cd); 01868 Particles[Particles.Size()-1]->v = Vec3_t(cd[0],cd[1],cd[2]); 01869 H5LTread_dataset_double(group_id,"w",cd); 01870 Particles[Particles.Size()-1]->w = Vec3_t(cd[0],cd[1],cd[2]); 01871 H5LTread_dataset_double(group_id,"wb",cd); 01872 Particles[Particles.Size()-1]->wb = Vec3_t(cd[0],cd[1],cd[2]); 01873 H5LTread_dataset_double(group_id,"I",cd); 01874 Particles[Particles.Size()-1]->I = Vec3_t(cd[0],cd[1],cd[2]); 01875 01876 double cq[4]; 01877 H5LTread_dataset_double(group_id,"Q",cq); 01878 Particles[Particles.Size()-1]->Q = Quaternion_t(cq[0],cq[1],cq[2],cq[3]); 01879 01880 // Loading the scalar quantities of the particle 01881 double dat[1]; 01882 H5LTread_dataset_double(group_id,"SR",dat); 01883 Particles[Particles.Size()-1]->Props.R = dat[0]; 01884 H5LTread_dataset_double(group_id,"Rho",dat); 01885 Particles[Particles.Size()-1]->Props.rho = dat[0]; 01886 H5LTread_dataset_double(group_id,"m",dat); 01887 Particles[Particles.Size()-1]->Props.m = dat[0]; 01888 H5LTread_dataset_double(group_id,"V",dat); 01889 Particles[Particles.Size()-1]->Props.V = dat[0]; 01890 H5LTread_dataset_double(group_id,"Diam",dat); 01891 Particles[Particles.Size()-1]->Diam = dat[0]; 01892 H5LTread_dataset_double(group_id,"Dmax",dat); 01893 Particles[Particles.Size()-1]->Dmax = dat[0]; 01894 int datint[1]; 01895 H5LTread_dataset_int(group_id,"Index",datint); 01896 Particles[Particles.Size()-1]->Index = datint[0]; 01897 int tag[1]; 01898 H5LTread_dataset_int(group_id,"Tag",tag); 01899 Particles[Particles.Size()-1]->Tag = tag[0]; 01900 Particles[Particles.Size()-1]->PropsReady = true; 01901 01902 } 01903 01904 01905 H5Fclose(file_id); 01906 01907 } 01908 #endif 01909 01910 inline void Domain::BoundingBox(Vec3_t & minX, Vec3_t & maxX) 01911 { 01912 minX = Vec3_t(Particles[0]->MinX(), Particles[0]->MinY(), Particles[0]->MinZ()); 01913 maxX = Vec3_t(Particles[0]->MaxX(), Particles[0]->MaxY(), Particles[0]->MaxZ()); 01914 for (size_t i=1; i<Particles.Size(); i++) 01915 { 01916 if (minX(0)>Particles[i]->MinX()&&Particles[i]->IsFree()) minX(0) = Particles[i]->MinX(); 01917 if (minX(1)>Particles[i]->MinY()&&Particles[i]->IsFree()) minX(1) = Particles[i]->MinY(); 01918 if (minX(2)>Particles[i]->MinZ()&&Particles[i]->IsFree()) minX(2) = Particles[i]->MinZ(); 01919 if (maxX(0)<Particles[i]->MaxX()&&Particles[i]->IsFree()) maxX(0) = Particles[i]->MaxX(); 01920 if (maxX(1)<Particles[i]->MaxY()&&Particles[i]->IsFree()) maxX(1) = Particles[i]->MaxY(); 01921 if (maxX(2)<Particles[i]->MaxZ()&&Particles[i]->IsFree()) maxX(2) = Particles[i]->MaxZ(); 01922 } 01923 } 01924 01925 inline void Domain::Center(Vec3_t C) 01926 { 01927 Vec3_t minX,maxX; 01928 BoundingBox(minX,maxX); 01929 Vec3_t Transport(-0.5*(maxX+minX)); 01930 Transport += C; 01931 for (size_t i=0; i<Particles.Size(); i++) Particles[i]->Translate(Transport); 01932 } 01933 01934 inline void Domain::ResetInteractons() 01935 { 01936 // delete old interactors 01937 for (size_t i=0; i<CInteractons.Size(); ++i) 01938 { 01939 if (CInteractons[i]!=NULL) delete CInteractons[i]; 01940 } 01941 01942 // new interactors 01943 CInteractons.Resize(0); 01944 for (size_t i=0; i<Particles.Size()-1; i++) 01945 { 01946 bool pi_has_vf = !Particles[i]->IsFree(); 01947 for (size_t j=i+1; j<Particles.Size(); j++) 01948 { 01949 bool pj_has_vf = !Particles[j]->IsFree(); 01950 01951 01952 // if both particles have any component specified or they are far away, don't create any intereactor 01953 bool close = (Distance(Particles[i]->x,Particles[j]->x)<=Particles[i]->Dmax+Particles[j]->Dmax+2*Alpha); 01954 if ((pi_has_vf && pj_has_vf) || !close ) continue; 01955 Listofpairs.insert(make_pair(Particles[i],Particles[j])); 01956 01957 // if both particles are spheres (just one vertex) 01958 if (Particles[i]->Verts.Size()==1 && Particles[j]->Verts.Size()==1) 01959 { 01960 CInteractons.Push (new CInteractonSphere(Particles[i],Particles[j])); 01961 } 01962 01963 // normal particles 01964 else 01965 { 01966 CInteractons.Push (new CInteracton(Particles[i],Particles[j])); 01967 } 01968 } 01969 } 01970 } 01971 01972 inline void Domain::ResetDisplacements() 01973 { 01974 for (size_t i=0; i<Particles.Size(); i++) 01975 { 01976 Particles[i]->ResetDisplacements(); 01977 } 01978 } 01979 01980 inline double Domain::MaxDisplacement() 01981 { 01982 double md = 0.0; 01983 for (size_t i=0; i<Particles.Size(); i++) 01984 { 01985 double mpd = Particles[i]->MaxDisplacement(); 01986 if (mpd > md) md = mpd; 01987 } 01988 return md; 01989 } 01990 01991 inline void Domain::ResetContacts() 01992 { 01993 for (size_t i=0; i<Particles.Size()-1; i++) 01994 { 01995 bool pi_has_vf = !Particles[i]->IsFree(); 01996 for (size_t j=i+1; j<Particles.Size(); j++) 01997 { 01998 bool pj_has_vf = !Particles[j]->IsFree(); 01999 02000 bool close = (Distance(Particles[i]->x,Particles[j]->x)<=Particles[i]->Dmax+Particles[j]->Dmax+2*Alpha); 02001 if ((pi_has_vf && pj_has_vf) || !close) continue; 02002 02003 // checking if the interacton exist for that pair of particles 02004 set<pair<Particle *, Particle *> >::iterator it = Listofpairs.find(make_pair(Particles[i],Particles[j])); 02005 if (it != Listofpairs.end()) 02006 { 02007 continue; 02008 } 02009 Listofpairs.insert(make_pair(Particles[i],Particles[j])); 02010 02011 // if both particles are spheres (just one vertex) 02012 if (Particles[i]->Verts.Size()==1 && Particles[j]->Verts.Size()==1) 02013 { 02014 CInteractons.Push (new CInteractonSphere(Particles[i],Particles[j])); 02015 } 02016 02017 // normal particles 02018 else 02019 { 02020 CInteractons.Push (new CInteracton(Particles[i],Particles[j])); 02021 } 02022 } 02023 } 02024 02025 Interactons.Resize(0); 02026 for (size_t i=0; i<CInteractons.Size(); i++) 02027 { 02028 if(CInteractons[i]->UpdateContacts(Alpha)) Interactons.Push(CInteractons[i]); 02029 } 02030 for (size_t i=0; i<BInteractons.Size(); i++) 02031 { 02032 if(BInteractons[i]->UpdateContacts(Alpha)) Interactons.Push(BInteractons[i]); 02033 } 02034 } 02035 02036 // Auxiliar methods 02037 02038 inline void Domain::LinearMomentum (Vec3_t & L) 02039 { 02040 L = 0.,0.,0.; 02041 for (size_t i=0; i<Particles.Size(); i++) 02042 { 02043 L += Particles[i]->Props.m*Particles[i]->v; 02044 } 02045 } 02046 02047 inline void Domain::AngularMomentum (Vec3_t & L) 02048 { 02049 L = 0.,0.,0.; 02050 for (size_t i=0; i<Particles.Size(); i++) 02051 { 02052 Vec3_t t1,t2; 02053 t1 = Particles[i]->I(0)*Particles[i]->w(0),Particles[i]->I(1)*Particles[i]->w(1),Particles[i]->I(2)*Particles[i]->w(2); 02054 Rotation (t1,Particles[i]->Q,t2); 02055 L += Particles[i]->Props.m*cross(Particles[i]->x,Particles[i]->v)+t2; 02056 } 02057 } 02058 02059 inline double Domain::CalcEnergy (double & Ekin, double & Epot) 02060 { 02061 // kinematic energy 02062 Ekin = 0.0; 02063 for (size_t i=0; i<Particles.Size(); i++) 02064 { 02065 Ekin += Particles[i]->Ekin + Particles[i]->Erot; 02066 } 02067 02068 // potential energy 02069 Epot = 0.0; 02070 for (size_t i=0; i<CInteractons.Size(); i++) 02071 { 02072 Epot += CInteractons[i]->Epot; 02073 } 02074 02075 // total energy 02076 return Ekin + Epot; 02077 } 02078 02079 inline void Domain::EnergyOutput (size_t IdxOut, std::ostream & OF) 02080 { 02081 // header 02082 if (IdxOut==0) 02083 { 02084 OF << Util::_10_6 << "Time" << Util::_8s << "Ekin" << Util::_8s << "Epot" << Util::_8s << "Evis" << Util::_8s << "Efric" << Util::_8s << "Wext" << std::endl; 02085 } 02086 double Ekin,Epot; 02087 CalcEnergy(Ekin,Epot); 02088 OF << Util::_10_6 << Time << Util::_8s << Ekin << Util::_8s << Epot << Util::_8s << Evis << Util::_8s << Efric << Util::_8s << Wext << std::endl; 02089 } 02090 02091 inline void Domain::GetGSD (Array<double> & X, Array<double> & Y, Array<double> & D, size_t NDiv) const 02092 { 02093 // calc GSD information 02094 Array<double> Vg; 02095 double Vs = 0.0; 02096 02097 for (size_t i=0; i<Particles.Size(); i++) 02098 { 02099 Particle * P = Particles[i]; 02100 double Diam = sqrt((P->MaxX()-P->MinX())*(P->MaxX()-P->MinX())+(P->MaxY()-P->MinY())*(P->MaxY()-P->MinY())+(P->MaxZ()-P->MinZ())*(P->MaxX()-P->MinX())); 02101 Vs += Particles[i]->Props.V; 02102 Vg.Push(Particles[i]->Props.V); 02103 D.Push(Diam); 02104 } 02105 double Dmin = D[D.TheMin()]; // minimum diameter 02106 double Dmax = D[D.TheMax()]; // maximum diameter 02107 double Dspan = (Dmax-Dmin)/NDiv; 02108 for (size_t i=0; i<=NDiv; i++) 02109 { 02110 X.Push (i*Dspan+Dmin); 02111 double cumsum = 0; 02112 for (size_t j=0; j<D.Size(); j++) 02113 { 02114 if (D[j]<=i*Dspan+Dmin) cumsum++; 02115 } 02116 Y.Push (cumsum/Particles.Size()); 02117 } 02118 } 02119 02120 inline void Domain::Clusters () 02121 { 02122 Array<int> connections; 02123 for (size_t i=0;i<BInteractons.Size();i++) 02124 { 02125 if (BInteractons[i]->valid) 02126 { 02127 connections.Push(BInteractons[i]->P1->Index); 02128 connections.Push(BInteractons[i]->P2->Index); 02129 } 02130 } 02131 02132 Util::Tree tree(connections); 02133 tree.GetClusters(Listofclusters); 02134 } 02135 02136 inline void Domain::DelParticles (Array<int> const & Tags) 02137 { 02138 Array<int> idxs; // indices to be deleted 02139 for (size_t i=0; i<Particles.Size(); ++i) 02140 { 02141 for (size_t j=0; j<Tags.Size(); ++j) 02142 { 02143 if (Particles[i]->Tag==Tags[j]) idxs.Push(i); 02144 } 02145 } 02146 if (idxs.Size()<1) throw new Fatal("Domain::DelParticles: Could not find any particle to be deleted"); 02147 Particles.DelItems (idxs); 02148 } 02149 02150 inline Particle * Domain::GetParticle (int Tag, bool Check) 02151 { 02152 size_t idx = 0; 02153 size_t count = 0; 02154 for (size_t i=0; i<Particles.Size(); ++i) 02155 { 02156 if (Particles[i]->Tag==Tag) 02157 { 02158 if (!Check) return Particles[i]; 02159 idx = i; 02160 count++; 02161 } 02162 } 02163 if (count==0) throw new Fatal("Domain::GetParticle: Could not find Particle with Tag==%d",Tag); 02164 else if (count>1) throw new Fatal("Domain::GetParticle: There are more than one particle with Tag==%d",Tag); 02165 return Particles[idx]; 02166 } 02167 02168 inline Particle const & Domain::GetParticle (int Tag, bool Check) const 02169 { 02170 size_t idx = 0; 02171 size_t count = 0; 02172 for (size_t i=0; i<Particles.Size(); ++i) 02173 { 02174 if (Particles[i]->Tag==Tag) 02175 { 02176 if (!Check) return (*Particles[i]); 02177 idx = i; 02178 count++; 02179 } 02180 } 02181 if (count==0) throw new Fatal("Domain::GetParticle: Could not find Particle with Tag==%d",Tag); 02182 else if (count>1) throw new Fatal("Domain::GetParticle: There are more than one particle with Tag==%d",Tag); 02183 return (*Particles[idx]); 02184 } 02185 02186 inline void Domain::GetParticles (int Tag, Array<Particle*> & P) 02187 { 02188 P.Resize(0); 02189 for (size_t i=0; i<Particles.Size(); ++i) 02190 { 02191 if (Particles[i]->Tag==Tag) P.Push(Particles[i]); 02192 } 02193 if (P.Size()==0) throw new Fatal("Domain::GetParticles: Could not find any Particle with Tag==%d",Tag); 02194 } 02195 02196 }; // namespace DEM 02197 02198 #endif // MECHSYS_DEM_DOMAIN_H
1.7.6.1
