MercuryDPM: MD Class Reference

References add_ene_ela(), BaseParticle::add_Force(), BaseParticle::add_Torque(), CSpecies::AdhesionForceType, computeShortRangeForceWithParticle(), Cross(), CTangentialSpring::delta, CSpecies::disp, CSpecies::disprolling, CSpecies::dispt, CSpecies::disptorsion, STD_save::fstat_file, gather_statistics_collision(), BaseParticle::get_AngularVelocity(), get_do_stat_always(), get_dt(), CSpecies::get_ForceType(), BaseParticle::get_Index(), BaseParticle::get_IndSpecies(), BaseParticle::get_InteractionRadius(), BaseParticle::get_Mass(), get_MixedSpecies(), BaseParticle::get_PeriodicFromParticle(), BaseParticle::get_Position(), BaseParticle::get_Radius(), get_save_data_ene(), get_save_data_fstat(), get_save_data_stat(), get_t(), BaseParticle::get_Velocity(), Vec3D::GetLength(), Vec3D::GetLength2, getTangentialSpring(), Hertzian, HM, HMD, BaseParticle::isFixed(), CSpecies::k, CSpecies::krolling, CSpecies::kt, CSpecies::ktorsion, CSpecies::mu, CSpecies::murolling, CSpecies::mus, CSpecies::musrolling, CSpecies::mustorsion, CSpecies::mutorsion, None, R, CTangentialSpring::RollingSpring, Vec3D::set_zero(), CTangentialSpring::sliding, CTangentialSpring::SlidingForce, CTangentialSpring::slidingRolling, CTangentialSpring::slidingTorsion, Species, sqr, CTangentialSpring::TorsionSpring, CSpecies::update_disp(), and CSpecies::useRestitution.

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

MD::add_Species

void add_Species(void)

Definition: MD.h:505

virtual void MD::broad_phase ( BaseParticle * i )

inlineprotectedvirtual

Broad phase of contact detection goes here. Default check all contacts.

Reimplemented in HGRID_base.

Definition at line 605 of file MD.h.

References BaseHandler< T >::begin(), compute_internal_forces(), and particleHandler.

Referenced by compute_internal_forces().

     {
         for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); (*it)!=i; ++it)
         {
             compute_internal_forces(i,*it);
         }
     }

int MD::Check_and_Duplicate_Periodic_Particles ( )

private

Calls Check_and_Duplicate_Periodic_Particle for all Particles curently in Particles[] and returns the number of particles created.

Definition at line 2719 of file MD.cc.

References boundaryHandler, BaseBoundary::createPeriodicParticles(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), and particleHandler.

Referenced by solve(), and statistics_from_restart_data().

     {
 //      std::cout<<"Largest R="<<particleHandler.getLargestParticle()->get_Radius()<<std::endl;
         int C=0;
         for (unsigned int k=0; k<boundaryHandler.getNumberOfObjects(); k++)
         {
 //          std::cout<<"Checking wall number: "<<k<<std::endl;
             int N=particleHandler.getNumberOfObjects(); //Required because the number of particles increaeses
             for (int i=0; i<N; i++)
             {
 //              std::cout<<"Particle number: "<<i<<std::endl;
                 C+=boundaryHandler.getObject(k)->createPeriodicParticles(particleHandler.getObject(i), particleHandler);
             }
         }
         return C;           
     }

void MD::checkInteractionWithBoundaries ( )

protectedvirtual

Definition at line 1794 of file MD.cc.

References BaseHandler< T >::begin(), boundaryHandler, BaseHandler< T >::end(), and particleHandler.

Referenced by solve().

 {
     for(std::vector<BaseBoundary*>::iterator B = boundaryHandler.begin(); B!=boundaryHandler.end(); ++B)
     {
         for (std::vector<BaseParticle*>::iterator P = particleHandler.begin(); P!=particleHandler.end(); ++P)
         {
             if((*B)->checkBoundaryAfterParticleMoved(*P,particleHandler)) //Returns true if the particle is deleted
                 --P;
         }
     }
 }

void MD::compute_all_forces ( )

protectedvirtual

This does the force computation.

Reset all forces to zero

Now loop over all particles contacts computing force contributions

Now loop over all other particles looking for contacts

Definition at line 1935 of file MD.cc.

References WallHandler::begin(), BaseHandler< T >::begin(), compute_external_forces(), compute_internal_forces(), WallHandler::end(), BaseHandler< T >::end(), PossibleContact::get_Next(), PossibleContact::get_P1(), PossibleContact::get_P2(), particleHandler, and wallHandler.

Referenced by solve(), and statistics_from_restart_data().

 {
     for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
     {
         (*it)->set_Force(Vec3D(0.0,0.0,0.0));
         (*it)->set_Torque(Vec3D(0.0,0.0,0.0));
     } 
     for (std::vector<BaseWall*>::iterator it = wallHandler.begin(); it!=wallHandler.end(); ++it)
     {
         (*it)->set_Force(Vec3D(0.0,0.0,0.0));
         (*it)->set_Torque(Vec3D(0.0,0.0,0.0));
     } //end reset forces loop
 
 
     #ifdef DEBUG_OUTPUT
         std::cerr << "Have all forces set to zero " << std::endl;
     #endif
 
 
     for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
     {
         
         compute_internal_forces(*it);
         //end inner loop over contacts.
             
         compute_external_forces(*it);
 
     }
 
 #ifdef ContactListHgrid     
     PossibleContact* Curr=possibleContactList.getFirstPossibleContact();
     while(Curr)
     {
         compute_internal_forces(Curr->get_P1(),Curr->get_P2());
         Curr=Curr->get_Next();
     }
 #endif
     //end outer loop over contacts.
 
 }

void MD::compute_external_forces ( BaseParticle * PI )

protectedvirtual

This is were the computation of external forces takes place (e.g. gravity)

This computes the external forces e.g. here it is gravity and walls.

Now add on gravity

Finally walls

Definition at line 1410 of file MD.cc.

References BaseParticle::add_Force(), compute_walls(), get_gravity(), BaseParticle::get_Mass(), and BaseParticle::isFixed().

Referenced by compute_all_forces().

 {
     if (!CI->isFixed()) {
         CI->add_Force(get_gravity() * CI->get_Mass());
         compute_walls(CI);
     }
 }

void MD::compute_internal_forces ( BaseParticle * i )

protectedvirtual

Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces)

Definition at line 1980 of file MD.cc.

References broad_phase().

Referenced by broad_phase(), HGRID_2D::CheckCell(), HGRID_3D::CheckCell(), HGRID_2D::CheckCell_current(), HGRID_3D::CheckCell_current(), and compute_all_forces().

                                                 {
     broad_phase(i);
 }

void MD::compute_internal_forces	(	BaseParticle *	P1,
		BaseParticle *	P2
	)

protectedvirtual

Computes the forces between particles (internal in the sence that the sum over all these forces is zero i.e. fully modelled forces)

This computer the internal forces (internal in the sence that they sum to zero) i.e. the fully modelled forces.

Tangential spring information is always store in the real particle with highest index When a Periodic contact is encountered it is always encoutered twice (for normal periodic walls), but the force is only applied to the real particle The tangential spring information for periodic particles is stored in the normal particle (and thus twice for normal periodic interactions) When a Particle is removed the tangential spring information has to be moved

Both options are up to first order the same (the first one is nicer becaus it always keeps the spring tangential, whereas the second one is in a nicer intergration form)

Todo:: {The spring should be cut back such that fdott=mu*fdotn. This is simple for dispt=0; we have to think about what happens in the sliding case with tang. dissipation; same for walls; Update Dec-2-2011: fixed}

Todo:: TW: the following 13 lines concern only the sliding spring and could be moved into the if statement above

Todo:: Define the center this way or are radii involved? Or maybe just use middle of overlap region? Thomas: Yes, we should correct that for polydispersed problems; however, then we also have to correct it in StatisticsVector::gather_statistics_collision.

The second particle (i.e. the particle the force acts on) is always a flow particle

Todo:: {Is it the first particle the force acts on?}

Definition at line 756 of file MD.cc.

 {
     //this is because the rough bottom allows overlapping fixed particles
     if (P1->isFixed()&&P2->isFixed()) return;
     
     
     //Dificult cases for tangential springs in comination with periodic walls are:
     //
     // A new contact over a periodic wall:
     // Starting situation: There are no tangential springs stored at all
     // Creating periodic particles: There are no tangential springs so nothing happens
     // Contact detection: There are two contacts detected, one (CA) for P1 with P2per and one (CB) for P2 with P1per
     // Switch to the 4 cases:
     //   CA: PI=P1, PJ=P2per PJreal=P2
     //   CB: PI=P2, PJ=P1per PJreal=P1
     // Reconnecting springs: 
     //   CA: PI=P1 and PJ=P2per do not have a spring, so a new spring is created in either PI or PJ (could be in periodic particle)
     //   CB: PI=P2 and PJ=P1Per do not have a spring, so a new spring is created in either PI or PJ (could be in periodic particle)
     // Removing periodic particles: One of the springs will be stored in a periodic particle and thus removed, the other spring is kept and used (this is the real particle with the kighest index)
     
     // Reconnect to a contact over a periodic wall:
     // Starting situation: There is a tangential spring stored in the particle with the highest index, lets assume this is P1
     // Creating periodic particles: P1 has a tangential spring, so P1per has a reversed copy of this spring
     // Contact detection: There are two contacts detected, one (CA) for P1 with P2per and one (CB) for P2 with P1per
     // Switch to the 4 cases:
     //   CA: PI=P1, PJ=P2per PJreal=P2
     //   CB: PI=P2, PJ=P1per PJreal=P1
     // Reconnecting springs: 
     //   CA: PI=P1 and PJ=P2per have a spring in P1, so we reconnect to this spring in the normal way and integrate it.
     //   CB: PI=P2 and PJ=P1Per have a spring in P1per, however this spring has to be a reversed spring since it is in PJ.
     // Removing periodic particles: The spring in P1per is removed
     
     //Cases (start from P1 and P2 and go to PI and PJ
     //1 Normal-Normal       ->PI=max(P1,P2), PJ=min(P1,P2)
     //2 Periodic-Normal     ->(PI=P2 PJ=real(P1))
     //3 Normal-Periodic     ->(PI=P1 PJ=real(P2))
     //4 Periodic-Periodic   ->do nothing
     
     //Just some statements to handle the 4 cases
     BaseParticle *PI, *PJ, *PJreal;
     bool isPeriodic;
     BaseParticle *P1Per=P1->get_PeriodicFromParticle();
     BaseParticle *P2Per=P2->get_PeriodicFromParticle();
     if(P1Per==NULL)
     {
         if(P2Per==NULL)//N-N
             if(P1->get_Index()<P2->get_Index())
                 {PI=P2;PJ=P1;PJreal=P1;isPeriodic=false;}
             else
                 {PI=P1;PJ=P2;PJreal=P2;isPeriodic=false;}
         else //N-P
             {PI=P1;PJ=P2;PJreal=P2Per;isPeriodic=true;}
     } else {
         if(P2Per==NULL) //P-N
             {PI=P2;PJ=P1;PJreal=P1Per;isPeriodic=true;}
         else return;
     }
     
 #ifdef DEBUG_OUTPUT
     std::cerr << "In computing internal forces between particle "<<PI->get_Position()<<" and "<<PJ->get_Position()<<std::endl;
 #endif
     
     //Get the square of the distance between particle i and particle j
     Mdouble dist_squared=GetDistance2(PI->get_Position(), PJ->get_Position());
     Mdouble interactionRadii_sum=PI->get_InteractionRadius()+PJ->get_InteractionRadius();
     
 #ifdef DEBUG_OUTPUT_FULL
     std::cerr << "Square of distance between " << dist_squared << " square sum of radii " << radii_sum*radii_sum <<std::endl;
 #endif
     
     // True if the particles are in contact
     if (dist_squared<(interactionRadii_sum*interactionRadii_sum))
     {
         // For particles of the same species, set species vector to Species(PI);
         // for particles of different species, set species vector to MixedSpecies(PI,PJ)
         CSpecies* pSpecies = (PI->get_IndSpecies()==PJ->get_IndSpecies())?&Species[PI->get_IndSpecies()]:get_MixedSpecies(PI->get_IndSpecies(),PJ->get_IndSpecies());
         
         // Calculate distance between the particles
         Mdouble dist=sqrt(dist_squared);
         
         // Compute normal vector
         
         Vec3D normal=(PI->get_Position()-PJ->get_Position())/dist;
         
         // Compute the overlap between the particles
         Mdouble radii_sum=PI->get_Radius()+PJ->get_Radius();
         Mdouble deltan = std::max(0.0,radii_sum-dist);
     
         Vec3D force= Vec3D(0,0,0);;
         Vec3D forcet; forcet.set_zero();
         Vec3D forcerolling; forcerolling.set_zero();
         Vec3D forcetorsion; forcetorsion.set_zero();
         Mdouble fdotn = 0;
         CTangentialSpring* TangentialSpring = NULL;
 
         //evaluate shortrange non-contact forces
         if (pSpecies->AdhesionForceType!=None)
             fdotn += computeShortRangeForceWithParticle(PI, PJ, PJreal, pSpecies, dist);
 
         if (deltan>0) //if contact forces
         {
 
             // Compute the relative velocity vector v_ij
             Vec3D vrel;
             if (!pSpecies->mu) {
                 vrel=(PI->get_Velocity()-PJ->get_Velocity());
             } else {
                 vrel=(PI->get_Velocity()-PJ->get_Velocity()) + Cross(normal, PI->get_AngularVelocity() * (PI->get_Radius() - .5 * deltan) + PJ->get_AngularVelocity() * (PJ->get_Radius() - .5 * deltan) );
             }
             
             // Compute the projection of vrel onto the normal (can be negative)
             Mdouble vdotn=-Dot(vrel,normal);
             
             //update restitution coeff
             if (pSpecies->useRestitution) pSpecies->update_disp(PI->get_Mass(),PJ->get_Mass());
             Mdouble a=0, R=0;
             
             // Compute normal force on particle i due to contact
             if (pSpecies->get_ForceType()==HM||pSpecies->get_ForceType()==HMD) 
             {
                 //R is twice the effective radius
                 R = 2.0*PI->get_Radius()*PJ->get_Radius()/(PI->get_Radius()+PJ->get_Radius());
                 a = sqrt(R*deltan);
                 //pSpecies->k stores the elastic modulus
                 Mdouble kn = 4./3. * pSpecies->k * a;
                 fdotn += kn * deltan + pSpecies->disp*vdotn;
             } else {
                 fdotn += pSpecies->k*deltan+pSpecies->disp*vdotn;
             }       
             force += normal * fdotn;
 
             //If tangential forces are present
             if (pSpecies->mu || pSpecies->murolling || pSpecies->mutorsion) { 
                 //call tangential spring
                 if (pSpecies->kt || pSpecies->krolling || pSpecies->ktorsion) 
                     TangentialSpring = getTangentialSpring(PI, PJ, PJreal);
                 
                 //Compute norm of normal force
                 Mdouble norm_fn = fabs(fdotn);
                 
                 //calculate sliding friction
                 if (pSpecies->mu) {
                     //Compute the tangential component of vrel
                     Vec3D vrelt=vrel+normal*vdotn;
                     //Compute norm of vrelt
                     Mdouble vdott=vrelt.GetLength();
                     
                     if (pSpecies->kt) {
                         Vec3D* delta = &(TangentialSpring->delta);
                         
                         //Integrate the spring
                         //(*delta) += vrelt * dt - Dot(*delta,normal)*normal;
                         Vec3D ddelta = (vrelt - Dot(*delta,PI->get_Velocity()-PJ->get_Velocity())*normal/dist) * get_dt();
                         (*delta) += ddelta;
                         
                         //Calculate test force including viscous force
                         if (pSpecies->get_ForceType()==HM) {
                             //pSpecies->kt stores the elastic shear modulus
                             Mdouble kt = 8. * pSpecies->kt * a;
                             TangentialSpring->SlidingForce += - kt * ddelta;
                             forcet = TangentialSpring->SlidingForce - pSpecies->dispt * vrelt;
                         } else if (pSpecies->get_ForceType()==HMD) {
                             //pSpecies->kt stores the elastic shear modulus
                             forcet = TangentialSpring->SlidingForce - pSpecies->dispt * vrelt;
                             Mdouble kt = 8. * pSpecies->kt * a * pow(1- GetLength(forcet)/pSpecies->mu/fdotn,0.33);
                             TangentialSpring->SlidingForce += - kt * ddelta;
                             forcet = TangentialSpring->SlidingForce - pSpecies->dispt * vrelt;
                         } else {
                             forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * (*delta);
                         }
 
                         Mdouble forcet2 = forcet.GetLength2();
                         
                         //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                         //but the force is limited by Coulomb friction (sliding):
                         //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                         double muact=(TangentialSpring->sliding)?(pSpecies->mu):(pSpecies->mus); // select mu from previous sliding mode
                         if( forcet2 <= sqr(muact*norm_fn) ) { 
                             //sticking++;
                             TangentialSpring->sliding=false;
                         } else { 
                             //sliding++;
                             TangentialSpring->sliding=true;
                             Mdouble norm_forcet = sqrt(forcet2);
                             forcet *= pSpecies->mu * norm_fn / norm_forcet;
                             TangentialSpring->SlidingForce = forcet + pSpecies->dispt * vrelt;
                             (*delta) = TangentialSpring->SlidingForce/(-pSpecies->kt);
                         }
                         
                         //Add tangential force to total force
                         force += forcet;
                         
                     } else { //if no tangential spring
                         //tangential forces are modelled by a damper of viscosity dispt (sticking), 
                         //but the force is limited by Coulomb friction (sliding):
                         //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                         if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
                             forcet = -pSpecies->dispt * vrelt; 
                         } else { //sliding++;
                             //set force to Coulomb limit
                             forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
                         }               
                         //Add tangential force to total force
                         force += forcet;
                     }
                 }
                 
                 //calculate rolling friction
                 if (pSpecies->murolling) {
                     //From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
                     Mdouble reducedRadiusI = PI->get_Radius() - .5 * deltan;
                     Mdouble reducedRadiusJ = PJ->get_Radius() - .5 * deltan;
                     Mdouble reducedRadiusIJ = 2.0*reducedRadiusI*reducedRadiusJ/(reducedRadiusI+reducedRadiusJ);
                     Vec3D vrolling=-reducedRadiusIJ * Cross(normal, PI->get_AngularVelocity() - PJ->get_AngularVelocity());
                     if (pSpecies->krolling) {
                         Vec3D* RollingSpring = &(TangentialSpring->RollingSpring);
                         
                         //Integrate the spring
                         (*RollingSpring) += vrolling * get_dt();// - Dot(*RollingSpring,normal)*normal;
                         
                         //Calculate test force including viscous force
                         forcerolling = (-pSpecies->disprolling) * vrolling - pSpecies->krolling * (*RollingSpring);
                         Mdouble forcerolling2 = forcerolling.GetLength2();
                         
                         //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                         //but the force is limited by Coulomb friction (sliding):
                         //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                         double muact=(TangentialSpring->slidingRolling)?(pSpecies->murolling):(pSpecies->musrolling); // select mu from previous sliding mode
                         if( forcerolling2 <= sqr(muact*norm_fn) ) { 
                             //sticking++;
                             TangentialSpring->slidingRolling=false;
                         } else { 
                             //sliding++;
                             TangentialSpring->slidingRolling=true;
                             forcerolling *= pSpecies->murolling * norm_fn / sqrt(forcerolling2);
                             (*RollingSpring) = (forcerolling + pSpecies->disprolling * vrolling)/(-pSpecies->krolling);
                         }
                         
                         //Add tangential force to torque
                         Vec3D Torque = reducedRadiusIJ * Cross(normal, forcerolling);
                         PI->add_Torque(Torque);
                         PJreal->add_Torque(-Torque);
                     }
                 }
 
                 
                 //calculate torsive friction
                 if (pSpecies->mutorsion) {
                     //From Luding2008, spin velocity (eq 16) w/o 2.0!
                     Mdouble RadiusIJ = 2.0*PI->get_Radius()*PJ->get_Radius()/(PI->get_Radius()+PJ->get_Radius());
                     Vec3D vtorsion=RadiusIJ*Dot(normal,PI->get_AngularVelocity() - PJ->get_AngularVelocity())*normal;
                     if (pSpecies->ktorsion) {
                         //~ std::cout << "Error; not yet implemented" << std::endl;
                         //~ exit(-1);
                         Vec3D* TorsionSpring = &(TangentialSpring->TorsionSpring);
                         
                         //Integrate the spring
                         (*TorsionSpring) = Dot((*TorsionSpring) + vtorsion * get_dt(),normal)*normal;
                         
                         //Calculate test force including viscous force
                         forcetorsion = (-pSpecies->disptorsion) * vtorsion - pSpecies->ktorsion * (*TorsionSpring);
                         Mdouble forcetorsion2 = forcetorsion.GetLength2();
                         
                         //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                         //but the force is limited by Coulomb friction (sliding):
                         //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                         double muact=(TangentialSpring->slidingTorsion)?(pSpecies->mutorsion):(pSpecies->mustorsion); // select mu from previous sliding mode
                         if( forcetorsion2 <= sqr(muact*norm_fn) ) { 
                             //sticking++;
                             TangentialSpring->slidingTorsion=false;
                         } else { 
                             //sliding++;
                             TangentialSpring->slidingTorsion=true;
                             //~ std::cout << "sliding " << std::endl;
                             forcetorsion *= pSpecies->mutorsion * norm_fn / sqrt(forcetorsion2);
                             (*TorsionSpring) = (forcetorsion + pSpecies->disptorsion * vtorsion)/(-pSpecies->ktorsion);
                         }
                         
                         //Add tangential force to torque
                         Vec3D torque = RadiusIJ * forcetorsion;
                         PI->add_Torque(torque);
                         PJreal->add_Torque(-torque);
                         
                     }
                 }           
             } //end if tangential forces
             
             //Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
             if (pSpecies->get_ForceType()==Hertzian) force *= sqrt(deltan/(PI->get_Radius()+PJ->get_Radius()));
         
         } else {//end if contact forces
             force += fdotn*normal;
         }
 
         //Add forces to total force
         PI->add_Force(force);
         if(!isPeriodic)
             PJreal->add_Force(-force);
         
         // Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
         if (pSpecies->mu) {
             Vec3D cross = Cross(normal, force);
             PI    ->add_Torque(-cross * (PI->get_Radius() - .5 * deltan));
             if(!isPeriodic)
                 PJreal->add_Torque(-cross * (PJ->get_Radius() - .5 * deltan));
         }
         
         // output for ene and stat files:
         if (get_save_data_ene()) {
             if(!isPeriodic)
                 add_ene_ela(0.5 * (pSpecies->k * sqr(deltan) + (TangentialSpring?
                                                              (pSpecies->kt       * TangentialSpring->delta.GetLength2()
                                                              +pSpecies->krolling * TangentialSpring->RollingSpring.GetLength2()
                                                              +pSpecies->ktorsion * TangentialSpring->TorsionSpring.GetLength2()):0.0)));
             else
                 add_ene_ela(0.25* (pSpecies->k * sqr(deltan) + (TangentialSpring?
                                                              (pSpecies->kt       * TangentialSpring->delta.GetLength2()
                                                              +pSpecies->krolling * TangentialSpring->RollingSpring.GetLength2()
                                                              +pSpecies->ktorsion * TangentialSpring->TorsionSpring.GetLength2()):0.0)));
         }
         if (get_save_data_fstat()||get_save_data_stat()||get_do_stat_always()) {            
             
             Mdouble fdott = forcet.GetLength();
             Mdouble deltat_norm = TangentialSpring?(-TangentialSpring->delta.GetLength()):0.0;
             
             Vec3D centre = 0.5 * (PI->get_Position() + PJ->get_Position());
             
             
             if (!PI->isFixed())
             {
                 if(get_save_data_stat()||get_do_stat_always()) gather_statistics_collision(PJreal->get_Index(),PI->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,-normal,-(fdott?forcet/fdott:forcet));
                 if(get_save_data_fstat()) fstat_file 
                     << get_t() << " "
                     << PJreal->get_Index() << " "
                     << PI->get_Index() << " "               
                     << centre << " "
                     << radii_sum-dist << " "
                     << deltat_norm << " "
                     << fdotn << " "
                     << fdott << " "
                     << -normal << " "
                     << -(fdott?forcet/fdott:forcet) << std::endl;
             }
             if (!PJreal->isFixed()&&!isPeriodic)
             {
                 if(get_save_data_stat()||get_do_stat_always()) gather_statistics_collision(PI->get_Index(),PJreal->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,normal,(fdott?forcet/fdott:forcet));
                 if(get_save_data_fstat()) fstat_file 
                     << get_t() << " "
                     << PI->get_Index() << " "
                     << PJreal->get_Index() << " "
                     << centre << " "
                     << radii_sum-dist << " "
                     << deltat_norm << " "
                     << fdotn << " "
                     << fdott << " "
                     << normal << " "
                     << (fdott?forcet/fdott:forcet) << std::endl;
             }
         }
 #ifdef FOLLOWPARTICLE
     if(PI->get_Index()==FPID||PJ->get_Index()==FPID)
         std::cout<<"Particle collission at time="<<get_t()<<" between "<<PI->get_Index()<<" and "<<PJ->get_Index()<<std::endl;
 #endif
         
     } // end if particle i and j are overlapping
 }

void MD::compute_particle_masses ( )

inlineprotected

Computes the mass of each particle.

Definition at line 618 of file MD.h.

References BaseParticle::compute_particle_mass(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), particleHandler, and Species.

Referenced by solve(), and statistics_from_restart_data().

618 {for (unsigned int i=0;i<particleHandler.getNumberOfObjects();i++) particleHandler.getObject(i)->compute_particle_mass(Species);}

BaseParticle::compute_particle_mass

void compute_particle_mass(std::vector< CSpecies > &Species)

Compute Particle mass function, which required a reference to the Species vector. It copmuters the Pa...

Definition: BaseParticle.cc:159

BaseHandler::getObject

T * getObject(const unsigned int id) const

Gets a pointer to the Object at the specified index in the BaseHandler.

Definition: BaseHandler.h:176

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

BaseHandler::getNumberOfObjects

ParticleHandler particleHandler

Definition: MD.h:707

unsigned int getNumberOfObjects() const

Gets the number of Object in this BaseHandler.

Definition: BaseHandler.h:199

void MD::compute_plastic_internal_forces	(	BaseParticle *	P1,
		BaseParticle *	P2
	)

protected

Computes plastic forces between particles.

This computer the internal forces for the plastic model.

Tangential spring information is always store in the real particle with highest index When a Periodic contact is encountered it is always enstd::coutered twice, but only applied when the real particle has the highest index of both real indices When a Particle is removed the tangential spring information has to be moved

newcode begin

Both options are up to first order the same (the first one is nicer becaus it always keeps the spring tangential, whereas the second one is in a nicer intergration form)

Todo:: {The spring should be cut back such that fdott=mu*fdotn. This is simple for dispt=0; we have to think about what happens in the sliding case with tang. dissipation; same for walls}

The second particle (i.e. the particle the force acts on) is always a flow particle

Definition at line 1139 of file MD.cc.

 {
     //this is because the rough bottom allows overlapping fixed particles
     if (P1->isFixed()&&P2->isFixed()) return;
 
 
     //Cases (start from P1 and P2 and go to PI and PJ
     //1 Normal-Normal       ->PI=max(P1,P2), PJ=min(P1,P2)
     //2 Periodic-Normal     ->if(P2>Real(P1)) (PI=P2 PJ=real(P1)) otherwise do nothing
     //3 Normal-Periodic     ->if(P1>Real(P2)) (PI=P1 PJ=real(P2)) otherwise do nothing
     //4 Periodic-Periodic   ->do nothing
     
     //Just some statements to handle the 4 cases
     BaseParticle *PI, *PJ,*PJreal;
     BaseParticle *P1Per=P1->get_PeriodicFromParticle();
     BaseParticle *P2Per=P2->get_PeriodicFromParticle();
     if(P1Per==NULL)
     {
         if(P2Per==NULL)
             //N-N
             if(P1->get_Index()<P2->get_Index())
                 {PI=P2;PJ=P1;PJreal=PJ;}
             else
                 {PI=P1;PJ=P2;PJreal=PJ;}
         else
             //N-P
             if(P1->get_Index()>P2Per->get_Index())
                 {PI=P1;PJ=P2;PJreal=P2Per;}
             else return;
     } else {
         if(P2Per==NULL)
             //P-N
             if(P2->get_Index()>P1Per->get_Index())
                 {PI=P2;PJ=P1;PJreal=P1Per;}
             else return;
         else return;
     }
 
     #ifdef DEBUG_OUTPUT
         std::cerr << "In computing interal forces between particle "<<PI->get_Position()<<" and "<<PJ->get_Position()<<std::endl;
     #endif
 
     
     //Get the square of the distance between particle i and particle j
     Mdouble dist_squared=GetDistance2(PI->get_Position(), PJ->get_Position());
     Mdouble interactionRadii_sum=PI->get_InteractionRadius()+PJ->get_InteractionRadius();
     
     #ifdef DEBUG_OUTPUT_FULL
         std::cerr << "Square of distance between " << dist_squared << " square sum of radii " << radii_sum*radii_sum <<std::endl;
     #endif
     
     // True if the particles are in contact
     if (dist_squared<(interactionRadii_sum*interactionRadii_sum))
     {
         // For particles of the same species, set species vector to Species(PI);
         // for particles of different species, set species vector to MixedSpecies(PI,PJ)
         CSpecies* pSpecies = (PI->get_IndSpecies()==PJ->get_IndSpecies())?&Species[PI->get_IndSpecies()]:get_MixedSpecies(PI->get_IndSpecies(),PJ->get_IndSpecies());
         
         // Calculate distance between the particles
         Mdouble dist=sqrt(dist_squared);
         
         // Compute normal vector
         Vec3D normal=(PI->get_Position()-PJ->get_Position())/dist;
         
         // Compute the overlap between the particles
         Mdouble radii_sum=PI->get_Radius()+PJ->get_Radius();
         Mdouble deltan = std::max(0.0,radii_sum-dist);
         
         Vec3D force= Vec3D(0,0,0);;
         Vec3D forcet; forcet.set_zero();
         Vec3D deltat; deltat.set_zero();
         Mdouble fdotn = 0;
         
         //evaluate shortrange non-contact forces
         if (pSpecies->AdhesionForceType!=None)
             fdotn += computeShortRangeForceWithParticle(PI, PJ, PJreal, pSpecies, dist);
         
         if (deltan>0) //if contact forces
         {
             
             // Compute the relative velocity vector v_ij
             Vec3D vrel;
             if (!pSpecies->mu) {
                 vrel=(PI->get_Velocity()-PJ->get_Velocity());
             } else {
                 vrel=(PI->get_Velocity()-PJ->get_Velocity()) + Cross(normal, PI->get_AngularVelocity() * (PI->get_Radius() - .5 * deltan) + PJ->get_AngularVelocity() * (PJ->get_Radius() - .5 * deltan) );
             }
             
             // BEGIN add viscous normal force
 
             // Compute the projection of vrel onto the normal (can be negative)
             Mdouble vdotn=-Dot(vrel,normal);
 
             // Compute normal force on particle i due to contact
             if (pSpecies->useRestitution) pSpecies->update_disp(PI->get_Mass(),PJ->get_Mass());
             fdotn += pSpecies->disp*vdotn;
 
             // END add viscous normal force
 
 
             // BEGIN add elastic force
             
             //calculate deltastar, the overlap above which k2max becomes active (the 'fluid branch')
             Mdouble deltastar = (pSpecies->k2max/(pSpecies->k2max-pSpecies->k))*pSpecies->depth*((2*PI->get_Radius()*PJ->get_Radius())/(PI->get_Radius()+PJ->get_Radius()));
             //2*depth*r_eff is the overlap at which fn=0 on the k2max branch
             //deltastar is the overlap at which the k1 and the k2max branch meet
 
             //retrieve history parameter deltamax, the max. achieved overlap
             DeltaMaxsParticle* DMParticle=dynamic_cast<DeltaMaxsParticle*>(PI);
             Mdouble *deltamax = DMParticle->get_DeltaMaxs().select_particle(PJreal->get_Index(), t, dt);
             //update deltastar
             *deltamax = std::min(deltastar,std::max(*deltamax,deltan)); 
             
             //calculate k2(deltamax) (only use first case for Walton-Braun spring)
             Mdouble k2;     
             if (deltan>deltastar) {
                 k2 = pSpecies->k2max;   
             } else {
                 k2 = pSpecies->k+(pSpecies->k2max-pSpecies->k)*((*deltamax)/deltastar);
             }       
                     
             //calculate delta0(deltamax), the overlap where the force is zero
             Mdouble delta0 = (k2-pSpecies->k)/k2*(*deltamax); 
             
             //add elastic force
             //std::cout << k2*(deltan-(delta0))-pSpecies->k*deltan << std::endl;
             if (deltan>deltastar) {
                 fdotn += pSpecies->k2max*(deltan-(delta0));
             } else if (k2*(deltan-(delta0))>=pSpecies->k*deltan){
                 fdotn += pSpecies->k*deltan;
             } else if (k2*(deltan-delta0)>=-pSpecies->kc*deltan){
                 fdotn += k2*(deltan-delta0);
             } else {
                 fdotn += -pSpecies->kc*deltan;
                 *deltamax=(k2+pSpecies->kc)/(k2-pSpecies->k)*deltan;
             }
             
             force = normal * fdotn;
             // END add elastic force
 
             
             //If tangential forces are present
             if (pSpecies->mu) {
                 
                 //Compute the tangential component of vrel
                 Vec3D vrelt=vrel+normal*vdotn;
 
                 //Compute norm of vrelt
                 Mdouble vdott=vrelt.GetLength();
 
                 //Compute norm of normal force
                 Mdouble norm_fn = fabs(fdotn);
                 
                 if (!pSpecies->kt) { //if no tangential spring
                     //tangential forces are modelled by a damper of viscosity dispt (sticking), 
                     //but the force is limited by Coulomb friction (sliding):
                     //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                     if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
                         forcet = -pSpecies->dispt * vrelt; 
                     } else { //sliding++;
                         //set force to Coulomb limit
                         forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
                     }
                     
                 } else { //if with tangential spring    
                     //Retrieve the spring
                     CTangentialSpring* TangentialSpring = getTangentialSpring(PI, PJ, PJreal);
                     Vec3D* delta = &(TangentialSpring->delta);
                     
                     //Integrate the spring
                     //(*delta) += vrelt * dt - Dot(*delta,normal)*normal;
                     (*delta) += (vrelt - Dot(*delta,PI->get_Velocity()-PJ->get_Velocity())*normal/dist) * dt;
                     deltat = (*delta);
                     
                     //Calculate test force including viscous force
                     forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * deltat;
                     Mdouble forcet2 = forcet.GetLength2();
 
                     //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                     //but the force is limited by Coulomb friction (sliding):
                     //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                     if( forcet2 <= sqr(pSpecies->mus*norm_fn) ) { 
                         //sticking++;
                     } else { 
                         //sliding++;
                         Mdouble norm_forcet = sqrt(forcet2);
                         forcet *= pSpecies->mu * norm_fn / norm_forcet;
                         (*delta) = -(forcet + pSpecies->dispt * vrelt)/pSpecies->kt;
                         //This line should be removed before release it is the old tangential model (the new is shown to be better).
                         //(*delta) = forcet/(-pSpecies->kt);
                     }
                 } //end if tangential spring
                 
                 //Add tangential force to total force
                 force += forcet;
 
             } //end if tangential forces
         } else {//end if contact forces
             force = fdotn*normal;
         }
         
         //Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
         if (pSpecies->get_ForceType()==Hertzian) force *= sqrt(deltan/(PI->get_Radius()+PJ->get_Radius()));
         
         PI    ->add_Force( force);
         PJreal->add_Force(-force);
                 
         // Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
         if (pSpecies->mu) {
             Vec3D cross = -Cross(normal, force);
             PI    ->add_Torque(cross * (PI->get_Radius() - .5 * deltan));
             PJreal->add_Torque(cross * (PI->get_Radius() - .5 * deltan));
         }
 
         // output for ene and stat files:
         if (save_data_ene) {
             ene_ela += 0.5 * (pSpecies->k * sqr(deltan) + pSpecies->kt * deltat.GetLength2());
         }
         if (save_data_fstat||save_data_stat||do_stat_always) {          
             Mdouble fdott = forcet.GetLength();
             Mdouble deltat_norm = -deltat.GetLength();
             
             Vec3D centre = 0.5 * (PI->get_Position() + PJ->get_Position());
             
             if (!PI->isFixed())
             {
                  if(save_data_stat||do_stat_always) gather_statistics_collision(PJreal->get_Index(),PI->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,-normal,-(fdott?forcet/fdott:forcet));
                  if(save_data_fstat) fstat_file 
                     << t << " "
                     << PJreal->get_Index() << " "
                     << PI->get_Index() << " "               
                     << centre << " "
                     << radii_sum-dist << " "
                     << deltat_norm << " "
                     << fdotn << " "
                     << fdott << " "
                     << -normal << " "
                     << -(fdott?forcet/fdott:forcet) << std::endl;
             }
             if (!PJreal->isFixed())
             {
                  if(save_data_stat||do_stat_always) gather_statistics_collision(PI->get_Index(),PJreal->get_Index(), centre, deltan, deltat_norm,fdotn,fdott,normal,(fdott?forcet/fdott:forcet));
                  if(save_data_fstat){ fstat_file 
                     << t << " "
                     << PI->get_Index() << " "
                     << PJreal->get_Index() << " "
                     << centre << " "
                     << radii_sum-dist << " "
                     << deltat_norm << " "
                     << fdotn << " "
                     << fdott << " "
                     << normal << " "
                     << (fdott?forcet/fdott:forcet) << std::endl;
                 }
             }
         }
         
     } // end if particle i and j are overlapping
 }

void MD::compute_walls ( BaseParticle * PI )

protectedvirtual

This is were the walls are.

This is were the walls are implemented - normals are outward normals.

Todo:: HMD still has to be fully implemented; HM & HMD requires thinking about how to keep it tangential; you dont need to update delta AND the force

Todo:: : reducedRadiusIJ should have a factor of 2.0, but then the rolling frequency differs from the normal frequency when krolling=2/5*k;

Todo:: : RadiusIJ should have a factor of 2.0, but then the rolling frequency differs from the normal frequency when krolling=2/5*k;

Definition at line 1423 of file MD.cc.

Referenced by compute_external_forces().

 {
     //No need to compute interactions between periodic particle images and walls
     if(pI->get_PeriodicFromParticle()!=NULL)
         return;
     
     for (unsigned int i=0; i<wallHandler.getNumberOfWalls(); i++) 
     {
         
         // note: normal points away from particle i
         Vec3D normal;
         Mdouble dist;
         bool touch = wallHandler.getWall(i)->get_distance_and_normal(*pI, dist, normal);
         
         //If the wall is being touched (I think :Ant)
         if (touch) 
         {
             
             CSpecies* pSpecies = (pI->get_IndSpecies()==wallHandler.getWall(i)->indSpecies)?&Species[pI->get_IndSpecies()]:get_MixedSpecies(pI->get_IndSpecies(),wallHandler.getWall(i)->indSpecies);
             
             Mdouble deltan = std::max(0.0,pI->get_Radius()-dist);
 
             Vec3D force= Vec3D(0,0,0);;
             Vec3D forcet; forcet.set_zero();
             Vec3D forcerolling; forcerolling.set_zero();
             Vec3D forcetorsion; forcetorsion.set_zero();
             Mdouble fdotn = 0;
             CTangentialSpring* TangentialSpring = NULL;
 
             //evaluate shortrange non-contact forces
             if (pSpecies->AdhesionForceType!=None)
                 fdotn += computeShortRangeForceWithWall(pI, i, pSpecies, dist);
             
             if (deltan>0) //if contact forces
             {
 
             
                 // Compute the relative velocity vector v_ij
                 Vec3D vrel;
                 if (!pSpecies->mu) {
                     vrel = pI->get_Velocity() - wallHandler.getWall(i)->get_Velocity();
                 } else {
                     vrel = pI->get_Velocity() - wallHandler.getWall(i)->get_Velocity() - Cross(normal, pI->get_AngularVelocity()) * dist;
                 }
                 
                 // Compute the projection of vrel onto the normal (can be negative)
                 Mdouble vdotn = -Dot(vrel, normal);
                 
                 //possibly update restitution coefficient
                 if (pSpecies->useRestitution) pSpecies->update_disp(pI->get_Mass(),1e20);
 
                 // Compute normal force on particle i due to contact
                 Mdouble a=0, R=0;
                 if (pSpecies->get_ForceType()==HM||pSpecies->get_ForceType()==HMD) {
                     //R is twice the effective radius
                     R = pI->get_Radius();
                     a = sqrt(deltan*R);
                     //pSpecies->k stores the elastic modulus
                     Mdouble kn = 4./3. * pSpecies->k * a; //it's not really kn
                     fdotn += kn*deltan - 2*pSpecies->disp*sqrt(pI->get_Mass()*2.*pSpecies->k * a)*vdotn;
                 } else {
                     fdotn += pSpecies->k*deltan - pSpecies->disp*vdotn;
                 }
                 force = normal * (-fdotn);
                 
                 //If tangential forces are present
                 if (pSpecies->mu || pSpecies->murolling || pSpecies->mutorsion) {
                     //call tangential spring
                     if (pSpecies->kt || pSpecies->krolling || pSpecies->ktorsion)
                         TangentialSpring = getTangentialSpringWall(pI, i);
                     
                     //Compute norm of normal force
                     Mdouble norm_fn = fabs(fdotn);
                     
                     //calculate sliding friction
                     if (pSpecies->mu) {
                         //Compute the tangential component of vrel
                         Vec3D vrelt=vrel+normal*vdotn;
                         //Compute norm of vrelt
                         Mdouble vdott=vrelt.GetLength();
                         
                         if (pSpecies->kt) {
                             Vec3D* delta = &(TangentialSpring->delta);
                             
                             //Integrate the spring
                             (*delta) += vrelt * get_dt(); //no correction because the normal direction is constant
                             //std::cout << (vrel) << std::endl;
                             //Calculate test force including viscous force
                             if (pSpecies->get_ForceType()==HM) {
                                 //pSpecies->kt stores the elastic shear modulus
                                 Mdouble kt = 8. * pSpecies->kt * a;
                                 //std::cout << "kt" << kt << std::endl;
                                 TangentialSpring->SlidingForce += - kt * (vrelt * get_dt());
                                 a = sqrt(deltan*R);
                                 forcet = TangentialSpring->SlidingForce - 2*pSpecies->dispt * sqrt(pI->get_Mass()*kt) * vrelt;
                             } else if (pSpecies->get_ForceType()==HMD) {  
                                 //pSpecies->kt stores the elastic shear modulus
                                 forcet = TangentialSpring->SlidingForce - pSpecies->dispt * vrelt;
                                 Mdouble kt;
                                 if (get_t()<0.034/10) {
                                     // loading
                                     kt = 8. * pSpecies->kt * a * pow(1.001- GetLength(forcet)/pSpecies->mus/fdotn,1./3.);
                                 } else if (Dot(vrelt,forcet)<=0) { 
                                     //unloading in opposite direction
                                     kt = 8. * pSpecies->kt * a * pow(.501 - .5 * GetLength(forcet)/pSpecies->mus/fdotn,1./3.);
                                 } else {
                                     //unloading
                                     kt = 8. * pSpecies->kt * a * pow(.501 + .5* GetLength(forcet)/pSpecies->mus/fdotn,1./3.);
                                 }
                                 TangentialSpring->SlidingForce += - kt * (vrelt * get_dt());
                                 forcet = TangentialSpring->SlidingForce - pSpecies->dispt * vrelt;
                             } else {
                                 forcet = (-pSpecies->dispt) * vrelt - pSpecies->kt * (*delta);
                             }
                             
                             Mdouble forcet2 = forcet.GetLength2();
                             
                             //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                             //but the force is limited by Coulomb friction (sliding):
                             //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                             double muact=(TangentialSpring->sliding)?(pSpecies->mu):(pSpecies->mus); // select mu from previous sliding mode
                             if( forcet2 <= sqr(muact*norm_fn) ) { 
                                 //sticking++;
                                 TangentialSpring->sliding=false;
                             } else { 
                                 //sliding++;
                                 TangentialSpring->sliding=true;
                                 Mdouble norm_forcet = sqrt(forcet2);
                                 forcet *= pSpecies->mu * norm_fn / norm_forcet;
                                 TangentialSpring->SlidingForce = forcet + pSpecies->dispt * vrelt;
                                 (*delta) = TangentialSpring->SlidingForce/(-pSpecies->kt);
                             }
                             
                             //Add tangential force to total force
                             force += forcet;
                             
                         } else { //if no tangential spring
                             //tangential forces are modelled by a damper of viscosity dispt (sticking), 
                             //but the force is limited by Coulomb friction (sliding):
                             //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                             if (vdott*pSpecies->dispt <= pSpecies->mus*norm_fn) { //sticking++;
                                 forcet = -pSpecies->dispt * vrelt; 
                             } else { //sliding++;
                                 //set force to Coulomb limit
                                 forcet = -(pSpecies->mu * norm_fn / vdott) * vrelt;
                                 //std::cout << "sliding " << std::endl;
                             }               
                             //Add tangential force to total force
                             force += forcet;
                         }
                     }
                     
                     //calculate rolling friction
                     if (pSpecies->murolling) {
                         //From Luding2008, objective rolling velocity (eq 15) w/o 2.0!
                         Mdouble reducedRadiusI = pI->get_Radius() - .5 * deltan;
                         Mdouble reducedRadiusIJ = reducedRadiusI; 
                         Vec3D vrolling=-reducedRadiusIJ * Cross(normal, pI->get_AngularVelocity());
                         if (pSpecies->krolling) {
                             Vec3D* RollingSpring = &(TangentialSpring->RollingSpring);
                             
                             //Integrate the spring
                             (*RollingSpring) += vrolling * get_dt();// - Dot(*RollingSpring,normal)*normal;
                             
                             //Calculate test force including viscous force
                             forcerolling = (-pSpecies->disprolling) * vrolling - pSpecies->krolling * (*RollingSpring);
                             Mdouble forcerolling2 = forcerolling.GetLength2();
                             
                             //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                             //but the force is limited by Coulomb friction (sliding):
                             //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                             double muact=(TangentialSpring->slidingRolling)?(pSpecies->murolling):(pSpecies->musrolling); // select mu from previous sliding mode
                             if( forcerolling2 <= sqr(muact*norm_fn) ) { 
                                 TangentialSpring->slidingRolling=false;
                             } else { 
                                 TangentialSpring->slidingRolling=true;
                                 forcerolling *= pSpecies->murolling * norm_fn / sqrt(forcerolling2);
                                 (*RollingSpring) = (forcerolling + pSpecies->disprolling * vrolling)/(-pSpecies->krolling);
                             }
                             
                             //Add tangential force to torque
                             Vec3D Torque = reducedRadiusIJ * Cross(normal, forcerolling);
                             pI->add_Torque(Torque);
                             wallHandler.getWall(i)->add_Torque(Torque);
 
                         }
                     }
                     
                     //calculate torsive friction
                     if (pSpecies->mutorsion) {
                         //From Luding2008, spin velocity (eq 16) w/o 2.0!
                         Mdouble RadiusIJ = pI->get_Radius();
                         Vec3D vtorsion=RadiusIJ*Dot(normal,pI->get_AngularVelocity())*normal;
                         if (pSpecies->ktorsion) {
                             //~ std::cout << "Error; not yet implemented" << std::endl;
                             //~ exit(-1);
                             Vec3D* TorsionSpring = &(TangentialSpring->TorsionSpring);
                             
                             //Integrate the spring
                             (*TorsionSpring) = Dot((*TorsionSpring) + vtorsion * get_dt(),normal)*normal;
                             
                             //Calculate test force including viscous force
                             forcetorsion = (-pSpecies->disptorsion) * vtorsion - pSpecies->ktorsion * (*TorsionSpring);
                             Mdouble forcetorsion2 = forcetorsion.GetLength2();
                             
                             //tangential forces are modelled by a spring-damper of elastisity kt and viscosity dispt (sticking), 
                             //but the force is limited by Coulomb friction (sliding):
                             //f_t = -dispt*vrelt, if dispt*vrelt<=mu_s*fdotn, f_t=mu+s*fdotn*t, else
                             double muact=(TangentialSpring->slidingTorsion)?(pSpecies->mutorsion):(pSpecies->mustorsion); // select mu from previous 
                             if( forcetorsion2 <= sqr(muact*norm_fn) ) { 
                                 //sticking++;
                                 TangentialSpring->slidingTorsion=false;
                             } else { 
                                 //sliding++;
                                 TangentialSpring->slidingTorsion=true;
                                 // std::cout << "torsion sliding " << std::endl;
                                 forcetorsion *= pSpecies->mutorsion * norm_fn / sqrt(forcetorsion2);
                                 (*TorsionSpring) = (forcetorsion + pSpecies->disptorsion * vtorsion)/(-pSpecies->ktorsion);
                             }
                             
                             //Add tangential force to torque
                             Vec3D Torque = RadiusIJ * forcetorsion;
                             
                             if (pSpecies->get_ForceType()==HM||pSpecies->get_ForceType()==HMD) {
                                 //R is twice the effective radius
                                 R = pI->get_Radius();
                                 a = sqrt(deltan*R);
                                 Torque = 20 * a * forcetorsion;
                             }
                             
                             pI->add_Torque(Torque);
                             wallHandler.getWall(i)->add_Torque(Torque);
                             //std::cout << " " << Torque.Z << std::endl;
                         }
                     }           
                     
                 } //end if tangential forces
             } else {//end if contact forces
                 force += (-fdotn)*normal;
             }
             
             //Make force Hertzian (note: deltan is normalized by the normal distanct of two particles in contact, as in Silbert)
             if (pSpecies->get_ForceType()==Hertzian) force *= sqrt(deltan/(2.0*pI->get_Radius())); 
             
             //Add forces to total force
             pI->add_Force(force);
             wallHandler.getWall(i)->add_Force(-force);
             
             if (pSpecies->mu) pI->set_Torque(pI->get_Torque()+Cross(normal, force) * dist);
             // Add torque due to tangential forces: t = Cross(l,f), l=dist*Wall.normal
             
             if (get_save_data_ene()) {
                 add_ene_ela(0.5 * (pSpecies->k * sqr(deltan) + (TangentialSpring?
                                                                 (pSpecies->kt * TangentialSpring->delta.GetLength2()
                                                                  +pSpecies->krolling * TangentialSpring->RollingSpring.GetLength2()
                                                                  +pSpecies->ktorsion * TangentialSpring->TorsionSpring.GetLength2()):0.0)));
             }
             if (get_save_data_fstat()||get_save_data_stat()||get_do_stat_always()) 
             {
                 
                 deltan = pI->get_Radius()-dist;
                 Mdouble fdott = forcet.GetLength();
                 Mdouble deltat_norm = TangentialSpring?(-TangentialSpring->delta.GetLength()):0.0;
                 
                 if (!pI->isFixed())
                 {
                     if(get_save_data_stat()||get_do_stat_always()) gather_statistics_collision(pI->get_Index(),-(i+1), pI->get_Position() + normal*(pI->get_Radius()-deltan), deltan, deltat_norm,fdotn,fdott,normal,-(fdott?forcet/fdott:forcet));
                     if(get_save_data_fstat()) fstat_file 
                         << get_t() << " "
                         << pI->get_Index() << " "
                         << -(static_cast<int>(i)+1) << " "
                         << pI->get_Position() + normal*(pI->get_Radius()-deltan) << " "
                         << deltan << " "
                         << deltat_norm << " "
                         << fdotn << " "
                         << fdott << " "
                         << normal << " "
                         << -(fdott?forcet/fdott:forcet) << std::endl;
                 }
             } 
             
         } // end if particle i touches the wall
         
     }//end if for wall i
 }

Mdouble MD::computeShortRangeForceWithParticle	(	BaseParticle *	PI,
		BaseParticle *	PJ,
		BaseParticle *	PJreal,
		CSpecies *	pSpecies,
		Mdouble	dist
	)

protected

Definition at line 1711 of file MD.cc.

References CSpecies::AdhesionForceType, CSpecies::f0, get_dt(), BaseParticle::get_Index(), BaseParticle::get_Radius(), get_t(), TangentialSpringParticle::get_TangentialSprings(), CSpecies::k0, LinearIrreversible, LinearReversible, and CTangentialSprings::select_particle_spring().

Referenced by compute_internal_forces(), and compute_plastic_internal_forces().

                                                                                                                                          {
     Mdouble fdotn;
     if (pSpecies->AdhesionForceType==LinearReversible)
     {
         fdotn = std::min(0.0,pSpecies->k0*std::max(dist-PI->get_Radius()-PJ->get_Radius(),0.0)-pSpecies->f0);
     }
     else if (pSpecies->AdhesionForceType==LinearIrreversible)
     {
         TangentialSpringParticle* TSParticleI=dynamic_cast<TangentialSpringParticle*>(PI);
         TangentialSpringParticle* TSParticleJ=dynamic_cast<TangentialSpringParticle*>(PJ);
         //First check if particle I has a spring stores that is connected to particle J
         CTangentialSpring* TangentialSpring = TSParticleI->get_TangentialSprings().select_particle_spring(PJreal->get_Index(), get_t(), get_dt());
         if(TangentialSpring==NULL)
         {
             //Then check if particle J has a spring stored that is connected to particle I
             TangentialSpring = TSParticleJ->get_TangentialSprings().select_particle_spring(PI->get_Index(), get_t(), get_dt());
         }
         if(TangentialSpring==NULL)
         {
             fdotn = (PI->get_Radius()+PJ->get_Radius()>dist)?-pSpecies->f0:0.0;
         }
         else
         {
             fdotn = std::min(0.0,pSpecies->k0*std::max(dist-PI->get_Radius()-PJ->get_Radius(),0.0)-pSpecies->f0);
         }
     }
     else
     {
         fdotn=0;
     }
     return fdotn;
 }

Mdouble MD::computeShortRangeForceWithWall	(	BaseParticle *	pI,
		int	wall,
		CSpecies *	pSpecies,
		Mdouble	dist
	)

protected

Definition at line 1744 of file MD.cc.

References CSpecies::AdhesionForceType, CSpecies::f0, get_dt(), BaseParticle::get_Radius(), get_t(), TangentialSpringParticle::get_TangentialSprings(), CSpecies::k0, LinearIrreversible, LinearReversible, and CTangentialSprings::select_wall_spring().

Referenced by compute_walls().

                                                                                                        {
     Mdouble fdotn;
     if (pSpecies->AdhesionForceType==LinearReversible)
     {
         fdotn = std::min(0.0,pSpecies->k0*std::max(dist-pI->get_Radius(),0.0)-pSpecies->f0);
     }
     else if (pSpecies->AdhesionForceType==LinearIrreversible)
     {
         TangentialSpringParticle* TSParticle=dynamic_cast<TangentialSpringParticle*>(pI);
         CTangentialSpring* TangentialSpring = TSParticle->get_TangentialSprings().select_wall_spring(wall, get_t(), get_dt());
         if(TangentialSpring==NULL)
         {
             fdotn = (pI->get_Radius()>dist)?-pSpecies->f0:0.0;
         }
         else
         {
             fdotn = std::min(0.0,pSpecies->k0*std::max(dist-pI->get_Radius(),0.0)-pSpecies->f0);
         }
     }
     else
     {
         fdotn=0;
     }
     return fdotn;
 }

void MD::constructor ( )

A public constructor, which sets defaults so the problem can be solved off the shelf.

todo{to incorporate changes in icc a make clean is required. Why?}

This is where the constructor is defines setup a basic two dimensional problem which can be solved off the shelf///

Definition at line 50 of file MD.cc.

References data_FixedParticles, dim, dt, format, gravity, particleHandler, PeriodicCreated, STD_save::problem_name, random, save_restart_data_counter, set_append(), set_dim_particle(), set_do_stat_always(), set_k(), RNG::set_RandomSeed(), set_restarted(), set_rho(), set_save_count_all(), ParticleHandler::setSpecies(), BaseHandler< T >::setStorageCapacity(), Species, tmax, xballs_additional_arguments, xballs_cmode, xballs_scale, xballs_vscale, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by MD().

 {
     Species.resize(1);
     set_restarted(false);
 
     //This sets the maximum number of particles
     particleHandler.setSpecies(&Species);
     particleHandler.setStorageCapacity(2);
     dim = 2;
     format = 0;
     
     //These are the particle parameters like dissipation etc..
     set_k(1e4);
     set_rho(2000);
     //\todo{Why shouldn't we set dim_particle (which defines the mass calculations) ALWAYS to three, even for 2D problems?}
     set_dim_particle(2); 
     
     // Set gravitational acceleration
     gravity = Vec3D(0.0, -9.8, 0.0);
     
     //This is the parameter of the numerical part
     dt=0e-5; // if dt is not user-specified, this is set in actions_before_time_loop()
     set_save_count_all(20);
     set_do_stat_always(false);
     tmax=1.0;
     
     //This sets the default xballs domain
     xmin=0.0;
     xmax=0.01;
     ymin=0.0;
     ymax=0.01;
     zmin=0.0;
     zmax=0.0;
 
     //This set the default position of walls
     /*Wall w0;
     w0.set(Vec3D(-1,0,0), -xmin);
     wallHandler.copyAndAddWall(w0);
     w0.set(Vec3D( 1,0,0),  xmax);
     wallHandler.copyAndAddWall(w0);
     w0.set(Vec3D(0,-1,0), -ymin);
     wallHandler.copyAndAddWall(w0);
     w0.set(Vec3D(0, 1,0),  ymax);
     wallHandler.copyAndAddWall(w0);*/
     
     problem_name << "Problem_1";
     data_FixedParticles = 0;
     
     //Default mode is energy with no scale of the vectors
     xballs_cmode=0;
     xballs_vscale=-1;
     xballs_scale=-1;
     xballs_additional_arguments = "";
     set_append(false);
     
     save_restart_data_counter = 0;
 
     PeriodicCreated=0;
 
     //The default random seed is 0
     random.set_RandomSeed(0);
     #ifdef DEBUG_OUTPUT
         std::cerr << "MD problem constructor finished " << std::endl;
     #endif
 }

virtual bool MD::continue_solve ( )

inlineprotectedvirtual

Definition at line 630 of file MD.h.

Referenced by solve().

630 {return true;}

virtual void MD::cout_time ( )

inlineprotectedvirtual

std::couts time

Reimplemented in Chute.

Definition at line 621 of file MD.h.

References t, and tmax.

Referenced by solve().

                              {
         std::cout << "\rt=" << std::setprecision(3) << std::left << std::setw(6) << t 
             << ", tmax=" << std::setprecision(3) << std::left << std::setw(6) << tmax << "        \r";
         std::cout.flush();
         #ifdef FOLLOWPARTICLE
             std::cout<<std::endl;
         #endif
     }

void MD::create_xballs_script ( )

virtual

This creates a scipt which can be used to load the xballs problem to display the data just generated.

This automatically creates an xballs script to plot the data you have just generated///.

First put in all the script lines. All these lines do is move you to the correct directory from any location

Todo:: {thomas:why does vscale have to be integer?}

Definition at line 84 of file MD_xballs.icc.

References STD_save::data_filename, dim, format, STD_save::get_options_data(), STD_save::problem_name, xballs_additional_arguments, xballs_cmode, xballs_scale, xballs_vscale, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by solve(), and StatisticsVector< T >::statistics_from_p3().

 {
     
     std::stringstream file_name;
     std::ofstream script_file;
     file_name << problem_name.str() <<".xballs";
     script_file.open((file_name.str()).c_str());
 
     script_file << "#!/bin/bash" << std::endl;
     script_file << "x=$(echo $0 | cut -c2-)" << std::endl;
     script_file << "file=$PWD$x" << std::endl;
     script_file << "dirname=`dirname \"$file\"`" << std::endl;
     script_file << "cd $dirname" << std::endl;
     
     Mdouble scale;
     int format;
 
     if (dim<3) 
         { // dim = 1 or 2
         format = 8;
         if (xballs_scale<0)
             {
             scale = 1.0 / std::max( ymax-ymin, xmax-xmin );
             }
         else
             {
             scale=xballs_scale;
             }
         } 
     else 
         { //dim==3
         format = 14;
         if (xballs_scale<0)
             {
                 scale = 1.2 / std::max( zmax-zmin,  xmax-xmin );
             }
         else
             {
                 scale=xballs_scale;
             }
         
         }
     
     script_file << "../../xballs/xballs -format " << format 
         << " -f " << data_filename.str() << ((get_options_data()==2)?".0000":"")
         << " -s " << scale
         << " -cmode " << xballs_cmode 
         << " -cmax -sort " 
         << xballs_additional_arguments 
         << " $*";
     if (xballs_vscale>-1)
         {
             script_file << " -vscale " << xballs_vscale;
         }
     script_file.close();
     
     //This line changes teh file permision and give the owener (i.e. you) read, write and excute permission to the file.
     #ifdef UNIX
         chmod((file_name.str().c_str()),S_IRWXU);
     #endif
 
 }

void MD::do_integration_after_force_computation ( BaseParticle * pI )

protectedvirtual

This is were the integration is done.

This is were the integration is done, at the moment it is Velocity Verlet integration.

Definition at line 1843 of file MD.cc.

References BaseParticle::accelerate(), BaseParticle::angularAccelerate(), boundaryHandler, PeriodicBoundary::distance(), dt, BaseParticle::get_AngularVelocity(), BaseParticle::get_Force(), get_HGRID_UpdateEachTimeStep(), BaseParticle::get_Index(), BaseParticle::get_InvInertia(), BaseParticle::get_InvMass(), BaseParticle::get_Position(), get_t(), BaseParticle::get_Torque(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), HGRID_RemoveParticleFromHgrid(), HGRID_UpdateParticleInHgrid(), BaseParticle::rotate(), rotation, BaseParticle::set_Position(), BaseParticle::set_Velocity(), and PeriodicBoundary::shift_position().

Referenced by solve().

 {
 #ifdef USE_SIMPLE_VERLET_INTEGRATION
 
     static Vec3D xtemp, atemp;
     xtemp=pI->get_Position();
     atemp=pI->get_Force()*pI->get_InvMass();
     pI->set_Position(xtemp*2.0-pI->PrevPosition+atemp*dt*dt);
     pI->set_Velocity((pI->get_Position()-pI->PrevPosition)/(2*dt)+atemp*dt);
     pI->PrevPosition=xtemp;
     if (rotation) { 
         pI->angularAccelerate(pI->get_Torque()*pI->get_InvInertia()*dt);
         pI->rotate(pI->get_AngularVelocity()*dt);
     }
     
     // This shifts particles that moved through periodic walls
     for (unsigned int j=0; j<boundaryHandler.getNumberOfObjects(); j++) 
     {
         PeriodicBoundary* b0=dynamic_cast<PeriodicBoundary*>(boundaryHandler.getObject(j));
         //Note, checking the dynamic cast allows you to detect if the wall is of a certain type (or inherited from this type).
         if(b0) //If the dyamic_cast succeeds 
             if (b0->distance(*iP)<0)
             {
                 b0->shift_position(iP->get_Position());
                 if(!get_HGRID_UpdateEachTimeStep())
                 {
                     HGRID_RemoveParticleFromHgrid(iP);
                     HGRID_UpdateParticleInHgrid(iP);
                 }
             }
     }   
     
 #else //use velocity verlet
     pI->accelerate(pI->get_Force()*pI->get_InvMass()*0.5*dt);
     if (rotation) pI->angularAccelerate(pI->get_Torque()*pI->get_InvInertia()*0.5*dt);;
 #endif
 #ifdef FOLLOWPARTICLE
     if(pI->get_Index()==FPID)
         std::cout<<"In integration after  time="<<get_t()<<", particle "<<FPID<<": "<<*pI<<std::endl;
 #endif
 }

void MD::do_integration_before_force_computation ( BaseParticle * pI )

protectedvirtual

This is were the integration is done.

This is were the integration is done, at the moment it is Velocity Verlet integration.

Definition at line 1773 of file MD.cc.

References BaseParticle::accelerate(), BaseParticle::angularAccelerate(), dt, BaseParticle::get_AngularVelocity(), BaseParticle::get_Displacement(), BaseParticle::get_Force(), BaseParticle::get_Index(), BaseParticle::get_InvInertia(), BaseParticle::get_InvMass(), get_t(), BaseParticle::get_Torque(), BaseParticle::get_Velocity(), Vec3D::GetLength(), HGRID_update_move(), BaseParticle::move(), BaseParticle::rotate(), rotation, and BaseParticle::set_Displacement().

Referenced by solve().

 {
 #ifdef USE_SIMPLE_VERLET_INTEGRATION
     
 #else //use velocity verlet
     iP->accelerate(iP->get_Force()*iP->get_InvMass()*0.5*dt);
     iP->set_Displacement(iP->get_Velocity()*dt);
     iP->move(iP->get_Displacement());
     HGRID_update_move(iP,iP->get_Displacement().GetLength());
     if (rotation)
     {   
         iP->angularAccelerate(iP->get_Torque()*iP->get_InvInertia()*0.5*dt);
         iP->rotate(iP->get_AngularVelocity()*dt);
     }
 #endif
 #ifdef FOLLOWPARTICLE
     if(iP->get_Index()==FPID)
         std::cout<<"In integration before time="<<get_t()<<", particle "<<FPID<<": "<<*iP<<std::endl;
 #endif
 }

bool MD::find_next_data_file	(	Mdouble	tmin,
		bool	verbose = `true`
	)

Definition at line 426 of file MD.cc.

References STD_save::data_file, STD_save::get_data_filename(), STD_save::get_file_counter(), STD_save::get_options_data(), STD_save::increase_counter_data(), STD_save::set_file_counter(), and t.

 {
     if (get_options_data()==2) {
         while(true) {
             increase_counter_data(std::fstream::in);
             //check if file exists and contains data
             int N;
             data_file >> N;     
             if (data_file.eof()||data_file.peek()==-1) {
                 std::cout << "file " << get_data_filename() << " not found" << std::endl;
                 return false;
             }
             //check if tmin condition is satisfied
             Mdouble t;
             data_file >> t;
             if (t>tmin) {
                 set_file_counter(get_file_counter()-1);
                 return true;
             }
             if (verbose) std::cout <<"Jumping file counter: "<<get_file_counter() << std::endl;
         }
     } else return true;
     //length = is.tellg();
     //is.seekg (0, ios::end);
 }

virtual void MD::finish_statistics ( )

inlineprotectedvirtual

Reimplemented in StatisticsVector< T >.

Definition at line 586 of file MD.h.

Referenced by solve().

586 {};

void MD::fstat_header ( )

protectedvirtual

Definition at line 145 of file MD.cc.

References STD_save::fstat_file, BaseParticle::get_Radius(), get_t(), get_xmax(), get_xmin(), get_ymax(), get_ymin(), get_zmax(), get_zmin(), getLargestParticle(), and getSmallestParticle().

Referenced by solve().

 {
     // line #1: time, volume fraction
     // line #2: wall box: wx0, wy0, wz0, wx1, wy1, wz1
     // line #3: radii-min-max & moments: rad_min, rad_max, r1, r2, r3, r4
     fstat_file << "#" 
         << " " << get_t()
         << " " << 0 
         << std::endl;
     fstat_file << "#" 
         << " " << get_xmin() 
         << " " << get_ymin() 
         << " " << get_zmin() 
         << " " << get_xmax()
         << " " << get_ymax() 
         << " " << get_zmax()
         << std::endl;
     fstat_file << "#" 
         << " " << getSmallestParticle()->get_Radius()
         << " " << getLargestParticle()->get_Radius()
         << " " << 0 
         << " " << 0 
         << " " << 0 
         << " " << 0 
         << std::endl;
 }

virtual void MD::gather_statistics_collision	(	int index1	UNUSED,
		int index2	UNUSED,
		Vec3D Contact	UNUSED,
		Mdouble delta	UNUSED,
		Mdouble ctheta	UNUSED,
		Mdouble fdotn	UNUSED,
		Mdouble fdott	UNUSED,
		Vec3D P1_P2_normal_	UNUSED,
		Vec3D P1_P2_tangential	UNUSED
	)

inlineprotectedvirtual

Reimplemented in StatisticsVector< T >.

Definition at line 584 of file MD.h.

Referenced by compute_internal_forces(), compute_plastic_internal_forces(), and compute_walls().

584 {};

bool MD::get_append ( )

inline

Gets restarted.

Definition at line 389 of file MD.h.

References append.

Referenced by solve(), and start_ene().

389 {return append;}

MD::append

bool append

Definition: MD.h:720

Mdouble MD::get_collision_time	(	Mdouble	mass,
		unsigned int	indSpecies = `0`
	)

inline

Calculates collision time for two copies of a particle of given disp, k, mass.

Definition at line 300 of file MD.h.

References Species.

Referenced by Chute::get_collision_time(), set_dt(), and set_dt_by_mass().

300 {return Species[indSpecies].get_collision_time(mass);}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_depth ( unsigned int indSpecies = 0 )

inline

Definition at line 198 of file MD.h.

References Species.

199 {return Species[indSpecies].get_depth();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

int MD::get_dim ( )

inline

Allows the dimension of the simulation to be accessed.

Definition at line 369 of file MD.h.

References dim.

Referenced by HGridOptimiser::Initialise(), output_xballs_data_particle(), and StatisticsVector< T >::set_nz().

369 {return dim;}

MD::dim

int dim

The dimension of the simulation.

Definition: MD.h:660

int MD::get_dim_particle ( unsigned int indSpecies = 0 )

inline

Allows the dimension of the particle (f.e. for mass) to be accessed.

Definition at line 263 of file MD.h.

References Species.

263 {return Species[indSpecies].get_dim_particle();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_disp ( unsigned int indSpecies = 0 )

inline

Allows the normal dissipation to be accessed.

Definition at line 239 of file MD.h.

References Species.

Referenced by Chute::readNextArgument(), and readNextArgument().

239 {return Species[indSpecies].get_dissipation();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_disprolling ( unsigned int indSpecies = 0 )

inline

Allows the tangential viscosity to be accessed.

Definition at line 230 of file MD.h.

References Species.

Referenced by readNextArgument().

230 {return Species[indSpecies].get_disprolling();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_dispt ( unsigned int indSpecies = 0 )

inline

Allows the tangential viscosity to be accessed.

Definition at line 226 of file MD.h.

References Species.

Referenced by Chute::readNextArgument(), and readNextArgument().

226 {return Species[indSpecies].get_dispt();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_disptorsion ( unsigned int indSpecies = 0 )

inline

Allows the tangential viscosity to be accessed.

Definition at line 234 of file MD.h.

References Species.

234 {return Species[indSpecies].get_disptorsion();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_dissipation ( unsigned int indSpecies = 0 )

inline

Allows the normal dissipation to be accessed.

Definition at line 243 of file MD.h.

References Species.

Referenced by Chute::set_collision_time_and_restitution_coefficient().

243 {return Species[indSpecies].get_dissipation();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

bool MD::get_do_stat_always ( )

inline

Definition at line 269 of file MD.h.

References do_stat_always.

Referenced by compute_internal_forces(), and compute_walls().

269 {return do_stat_always;}

MD::do_stat_always

bool do_stat_always

Definition: MD.h:683

Mdouble MD::get_dt ( )

inline

Allows the time step dt to be accessed.

Definition at line 339 of file MD.h.

References dt.

Referenced by StatisticsVector< T >::check_current_time_for_statistics(), compute_internal_forces(), compute_walls(), computeShortRangeForceWithParticle(), computeShortRangeForceWithWall(), getTangentialSpring(), getTangentialSpringWall(), load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().

339 {return dt;}

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

Mdouble MD::get_ene_ela ( )

inline

Gets ene_ela.

Definition at line 395 of file MD.h.

References ene_ela.

Referenced by output_ene().

395 {return ene_ela;}

MD::ene_ela

Mdouble ene_ela

used in force calculations

Definition: MD.h:690

int MD::get_format ( )

inline

Definition at line 508 of file MD.h.

References format.

Referenced by output_xballs_data_particle().

508 {return format;}

MD::format

int format

Definition: MD.h:699

Vec3D MD::get_gravity ( )

inline

Allows the gravitational acceleration to be accessed.

Definition at line 364 of file MD.h.

References gravity.

Referenced by compute_external_forces(), output_ene(), and Chute::set_ChuteAngle().

364 {return gravity;}

MD::gravity

Vec3D gravity

Gravitational acceleration.

Definition: MD.h:663

virtual bool MD::get_HGRID_UpdateEachTimeStep ( )

inlineprotectedvirtual

Definition at line 561 of file MD.h.

Referenced by do_integration_after_force_computation().

561 {return true;};

Mdouble MD::get_k ( int indSpecies = 0 )

inline

Allows the spring constant to be accessed.

Definition at line 206 of file MD.h.

References Species.

Referenced by Chute::readNextArgument(), and readNextArgument().

206 {return Species[indSpecies].get_k();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_k1 ( unsigned int indSpecies = 0 )

inline

Allows the plastic constants to be accessed.

Definition at line 192 of file MD.h.

References Species.

193 {return Species[indSpecies].get_k1();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_k2max ( unsigned int indSpecies = 0 )

inline

Definition at line 194 of file MD.h.

References Species.

Referenced by read_next_from_data_file().

195 {return Species[indSpecies].get_k2max();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_kc ( unsigned int indSpecies = 0 )

inline

Definition at line 196 of file MD.h.

References Species.

197 {return Species[indSpecies].get_kc();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_krolling ( int indSpecies = 0 )

inline

Allows the spring constant to be accessed.

Definition at line 214 of file MD.h.

References Species.

Referenced by readNextArgument().

214 {return Species[indSpecies].get_krolling();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_kt ( int indSpecies = 0 )

inline

Allows the spring constant to be accessed.

Definition at line 210 of file MD.h.

References Species.

Referenced by Chute::readNextArgument(), and readNextArgument().

210 {return Species[indSpecies].get_kt();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_ktorsion ( int indSpecies = 0 )

inline

Allows the spring constant to be accessed.

Definition at line 218 of file MD.h.

References Species.

218 {return Species[indSpecies].get_ktorsion();}

BaseParticle::get_IndSpecies

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_Mass_from_Radius	(	Mdouble	radius,
		int	indSpecies = `0`
	)

inline

Definition at line 407 of file MD.h.

References BaseParticle::compute_particle_mass(), BaseParticle::get_Mass(), BaseParticle::set_IndSpecies(), BaseParticle::set_Radius(), and Species.

                                                                    {
         BaseParticle P;
         P.set_Radius(radius);
         P.set_IndSpecies(indSpecies);
         P.compute_particle_mass(Species);
         return P.get_Mass();
     }

Mdouble MD::get_max_radius ( )

inline

Sets restarted.

Get maximum radius

Definition at line 381 of file MD.h.

References max_radius.

381 {return max_radius;}

MD::max_radius

Mdouble max_radius

Definition: MD.h:665

Mdouble MD::get_maximum_velocity ( BaseParticle & P )

inline

Calculates the maximum velocity allowed for a collision of two copies of P (for higher velocities particles could pass through each other)

Definition at line 416 of file MD.h.

References BaseParticle::get_IndSpecies(), BaseParticle::get_Mass(), BaseParticle::get_Radius(), and Species.

416 {return Species[P.get_IndSpecies()].get_maximum_velocity(P.get_Radius(),P.get_Mass());}

int get_IndSpecies() const

Definition: BaseParticle.cc:241

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

BaseParticle::get_Radius

Mdouble get_Radius() const

Definition: BaseParticle.cc:233

BaseParticle::get_Mass

Mdouble get_Mass() const

Definition: BaseParticle.cc:230

Mdouble MD::get_maximum_velocity ( )

inline

Todo:: this function is dangerous since it uses particle data, which the user might not intend Calculates the maximum velocity allowed for a collision of two copies of the smallest particle (for higher velocities particles could pass through each other)

Definition at line 431 of file MD.h.

References BaseParticle::get_IndSpecies(), BaseParticle::get_Mass(), BaseParticle::get_Radius(), getSmallestParticle(), and Species.

                                   {
         BaseParticle *P = getSmallestParticle();
         return P->get_Radius() * sqrt(Species[P->get_IndSpecies()].k/(.5*P->get_Mass()));
     }

CSpecies* MD::get_MixedSpecies	(	int	i,
		int	j
	)

inline

Allows the mixed species to be accessed.

Definition at line 131 of file MD.h.

References Species.

Referenced by compute_internal_forces(), compute_plastic_internal_forces(), and compute_walls().

                                              {
         if (i>j) return &Species[i].MixedSpecies[j];
         else return &Species[j].MixedSpecies[i];
     };

Mdouble MD::get_mu ( unsigned int indSpecies = 0 )

inline

Allows the Coulomb friction coefficient to be accessed.

Definition at line 247 of file MD.h.

References Species.

Referenced by read_next_from_data_file(), readNextArgument(), and solve().

247 {return Species[indSpecies].get_mu();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_murolling ( unsigned int indSpecies = 0 )

inline

Allows the Coulomb friction coefficient to be accessed.

Definition at line 251 of file MD.h.

References Species.

Referenced by readNextArgument().

251 {return Species[indSpecies].get_murolling();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_mutorsion ( unsigned int indSpecies = 0 )

inline

Allows the Coulomb friction coefficient to be accessed.

Definition at line 255 of file MD.h.

References Species.

255 {return Species[indSpecies].get_mutorsion();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

int MD::get_NSpecies ( )

inline

Allows the number of Species to be accessed.

Definition at line 123 of file MD.h.

References Species.

123 {return Species.size();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_plastic_dt	(	Mdouble	mass,
		unsigned int	indSpecies = `0`
	)

inline

Definition at line 200 of file MD.h.

References Species.

201 {return Species[indSpecies].get_plastic_dt(mass);}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

int MD::get_restart_version ( )

inline

Gets restart_version.

Definition at line 372 of file MD.h.

References restart_version.

Referenced by Chute::read().

372 {return restart_version;}

MD::restart_version

int restart_version

Definition: MD.h:700

bool MD::get_restarted ( )

inline

Gets restarted.

Definition at line 377 of file MD.h.

References restarted.

377 {return restarted;}

MD::restarted

bool restarted

Definition: MD.h:701

Mdouble MD::get_restitution_coefficient	(	Mdouble	mass,
		unsigned int	indSpecies = `0`
	)

inline

Calculates restitution coefficient for two copies of given disp, k, mass.

Definition at line 302 of file MD.h.

References Species.

Referenced by Chute::get_restitution_coefficient().

302 {return Species[indSpecies].get_restitution_coefficient(mass);}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_rho ( int indSpecies = 0 )

inline

Allows the density to be accessed.

Definition at line 222 of file MD.h.

References Species.

Referenced by Chute::set_collision_time_and_restitution_coefficient().

222 {return Species[indSpecies].get_rho();}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

bool MD::get_rotation ( )

inline

Definition at line 258 of file MD.h.

References rotation.

258 {return rotation;}

MD::rotation

bool rotation

Definition: MD.h:721

int MD::get_save_count ( )

inline

Definition at line 163 of file MD.h.

References get_save_count_data().

163 {return get_save_count_data ();}

MD::get_save_count_data

int get_save_count_data()

Definition: MD.h:164

int MD::get_save_count_data ( )

inline

Definition at line 164 of file MD.h.

References save_count_data.

Referenced by get_save_count(), and get_savecount().

164 {return save_count_data;}

int save_count_data

These are informations for saving.

Definition: MD.h:675

int MD::get_save_count_ene ( )

inline

Definition at line 165 of file MD.h.

References save_count_ene.

165 {return save_count_ene;}

int save_count_ene

Definition: MD.h:676

int MD::get_save_count_fstat ( )

inline

Definition at line 167 of file MD.h.

References save_count_fstat.

167 {return save_count_fstat;}

int save_count_fstat

Definition: MD.h:678

int MD::get_save_count_stat ( )

inline

Definition at line 166 of file MD.h.

References save_count_stat.

166 {return save_count_stat;}

int save_count_stat

Definition: MD.h:677

bool MD::get_save_data_data ( )

inline

Returns the data counter.

Definition at line 265 of file MD.h.

References save_data_data.

265 {return save_data_data;}

MD::save_data_data

bool save_data_data

Definition: MD.h:679

bool MD::get_save_data_ene ( )

inline

Definition at line 266 of file MD.h.

References save_data_ene.

Referenced by compute_internal_forces(), and compute_walls().

266 {return save_data_ene;}

MD::save_data_ene

bool save_data_ene

Definition: MD.h:680

bool MD::get_save_data_fstat ( )

inline

Definition at line 267 of file MD.h.

References save_data_fstat.

Referenced by compute_internal_forces(), and compute_walls().

267 {return save_data_fstat;}

MD::save_data_fstat

bool save_data_fstat

Definition: MD.h:681

bool MD::get_save_data_stat ( )

inline

Definition at line 268 of file MD.h.

References save_data_stat.

Referenced by compute_internal_forces(), and compute_walls().

268 {return save_data_stat;}

MD::save_data_stat

bool save_data_stat

Definition: MD.h:682

int MD::get_savecount ( )

inline

Allows the number of time steps between saves to be accessed.

Definition at line 162 of file MD.h.

References get_save_count_data().

162 {return get_save_count_data ();}

MD::get_save_count_data

int get_save_count_data()

Definition: MD.h:164

std::vector<CSpecies>& MD::get_Species ( void )

inline

Allows the species to be copied.

Definition at line 126 of file MD.h.

References Species.

126 {return Species;}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

CSpecies* MD::get_Species ( int i )

inline

Allows the species to be accessed.

Definition at line 128 of file MD.h.

References Species.

128 {return &Species[i];}

Referenced by StatisticsVector< T >::check_current_time_for_statistics(), compute_internal_forces(), compute_walls(), computeShortRangeForceWithParticle(), computeShortRangeForceWithWall(), Chute::cout_time(), do_integration_after_force_computation(), do_integration_before_force_computation(), StatisticsVector< T >::evaluate_particle_statistics(), fstat_header(), getTangentialSpring(), getTangentialSpringWall(), and output_ene().

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

Mdouble MD::get_t ( )

inline

Access function for the time.

Definition at line 108 of file MD.h.

References t.

108 {return t;}

MD::t

Mdouble t

This stores the current time.

Definition: MD.h:687

Mdouble MD::get_tmax ( )

inline

Allows the upper time limit to be accessed.

Definition at line 144 of file MD.h.

References tmax.

Referenced by Chute::cout_time(), StatisticsVector< T >::get_tmaxStat(), and readNextArgument().

144 {return tmax;}

MD::xballs_additional_arguments

Mdouble tmax

Definition: MD.h:672

std::string MD::get_xballs_additional_arguments ( )

inline

Definition at line 355 of file MD.h.

References xballs_additional_arguments.

355 {return xballs_additional_arguments;}

std::string xballs_additional_arguments

Definition: MD.h:698

int MD::get_xballs_cmode ( )

inline

Definition at line 347 of file MD.h.

References xballs_cmode.

347 {return xballs_cmode;}

MD::xballs_cmode

int xballs_cmode

Definition: MD.h:695

double MD::get_xballs_scale ( )

inline

Definition at line 359 of file MD.h.

References xballs_scale.

359 {return xballs_scale;}

MD::xballs_scale

Mdouble xballs_scale

Definition: MD.h:697

double MD::get_xballs_vscale ( )

inline

Definition at line 351 of file MD.h.

References xballs_vscale.

351 {return xballs_vscale;}

MD::xballs_vscale

Mdouble xballs_vscale

Definition: MD.h:696

Mdouble MD::get_xmax ( )

inline

Get xmax.

Definition at line 307 of file MD.h.

References xmax.

Referenced by Chute::clean_chute(), Chute::create_bottom(), fstat_header(), Chute::get_ChuteLength(), ChuteWithHopperAndInset::get_ChuteLength(), ChuteWithHopper::get_ChuteLength(), ChuteWithHopperAndInset::get_MaximumVelocityInducedByGravity(), ChuteWithHopper::get_MaximumVelocityInducedByGravity(), StatisticsVector< T >::get_xmaxStat(), HGridOptimiser::Initialise(), load_par_ini_file(), ChuteWithHopperAndInset::set_shift(), ChuteWithHopper::set_shift(), and ChuteBottom::setup_particles_initial_conditions().

307 {return xmax;}

MD::xmax

Mdouble xmax

Definition: MD.h:668

Mdouble MD::get_xmin ( )

inline

Get xmin.

Definition at line 305 of file MD.h.

References xmin.

Referenced by Chute::clean_chute(), Chute::create_bottom(), Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_xminStat(), HGridOptimiser::Initialise(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().

305 {return xmin;}

MD::xmin

Mdouble xmin

These store the size of the domain, assume walls at the ends.

Definition: MD.h:668

Mdouble MD::get_ymax ( )

inline

Gets ymax.

Definition at line 311 of file MD.h.

References ymax.

Referenced by ChuteWithHopperAndInset::add_hopper(), ChuteWithHopper::add_hopper(), Chute::create_bottom(), Chute::create_inflow_particle(), ChuteWithHopper::create_inflow_particle(), ChuteWithHopperAndInset::create_inflow_particle(), fstat_header(), Chute::get_ChuteWidth(), StatisticsVector< T >::get_ymaxStat(), HGridOptimiser::Initialise(), load_par_ini_file(), Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().

311 {return ymax;}

MD::ymax

Mdouble ymax

Definition: MD.h:668

Mdouble MD::get_ymin ( )

inline

Gets ymin.

Definition at line 309 of file MD.h.

References ymin.

Referenced by ChuteWithHopperAndInset::add_hopper(), ChuteWithHopper::add_hopper(), Chute::create_bottom(), Chute::create_inflow_particle(), ChuteWithHopper::create_inflow_particle(), ChuteWithHopperAndInset::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_yminStat(), HGridOptimiser::Initialise(), load_par_ini_file(), Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().

309 {return ymin;}

BaseParticle::get_Species

Mdouble ymin

Definition: MD.h:668

Mdouble MD::get_zmax ( )

inline

Gets zmax.

Definition at line 315 of file MD.h.

References zmax.

Referenced by Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_zmaxStat(), HGridOptimiser::Initialise(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().

315 {return zmax;}

MD::zmax

Mdouble zmax

Definition: MD.h:668

Mdouble MD::get_zmin ( )

inline

Gets zmin.

Definition at line 313 of file MD.h.

References zmin.

Referenced by Chute::create_bottom(), Chute::create_inflow_particle(), fstat_header(), StatisticsVector< T >::get_zminStat(), HGridOptimiser::Initialise(), load_par_ini_file(), and ChuteBottom::setup_particles_initial_conditions().

313 {return zmin;}

MD::zmin

Mdouble zmin

Definition: MD.h:668

BoundaryHandler& MD::getBoundaryHandler ( )

inline

Definition at line 149 of file MD.h.

References boundaryHandler.

Referenced by Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().

149 { return boundaryHandler;}

MD::boundaryHandler

BoundaryHandler boundaryHandler

Definition: MD.h:709

virtual double MD::getInfo ( BaseParticle & P )

inlinevirtual

Allows the user to set what is written into the info column in the data file. By default is store the Species ID number.

Definition at line 463 of file MD.h.

References BaseParticle::get_Species().

Referenced by output_xballs_data_particle().

463 {return P.get_Species();}

int get_Species() const

Definition: BaseParticle.cc:240

virtual BaseParticle* MD::getLargestParticle ( )

inlinevirtual

Reimplemented in Chute.

Definition at line 419 of file MD.h.

References ParticleHandler::getLargestParticle(), and particleHandler.

Referenced by fstat_header(), HGRID_base::InitBroadPhase(), and solve().

419 {return particleHandler.getLargestParticle();}

ParticleHandler::getLargestParticle

BaseParticle * getLargestParticle() const

Gets a pointer to the largest BaseParticle (by interactionRadius) in this ParticleHandler.

Definition: ParticleHandler.cc:213

Referenced by ChuteWithHopper::add_hopper(), Chute::add_particle(), Chute::clean_chute(), Chute::cout_time(), Chute::create_bottom(), Chute::get_LargestParticleInteractionRadius(), Chute::get_LightestParticleMass(), Chute::get_SmallestParticleInteractionRadius(), Chute::getLargestParticle(), Chute::getSmallestParticle(), HGRID_base::HGRID_actions_before_time_step(), HGRID_base::InitBroadPhase(), HGridOptimiser::Initialise(), Chute::IsInsertable(), ChuteBottom::make_rough_bottom(), output_ene(), output_xballs_data_particle(), HGRID_base::set_HGRID_num_buckets_to_power(), ChuteBottom::setup_particles_initial_conditions(), and solve().

ParticleHandler particleHandler

Definition: MD.h:707

ParticleHandler& MD::getParticleHandler ( )

inline

Definition at line 147 of file MD.h.

References particleHandler.

147 { return particleHandler;}

ParticleHandler particleHandler

Definition: MD.h:707

virtual BaseParticle* MD::getSmallestParticle ( )

inlinevirtual

Reimplemented in Chute.

Definition at line 418 of file MD.h.

References ParticleHandler::getSmallestParticle(), and particleHandler.

Referenced by fstat_header(), get_maximum_velocity(), HGRID_base::InitBroadPhase(), and set_dt().

418 {return particleHandler.getSmallestParticle();}

ParticleHandler::getSmallestParticle

ParticleHandler particleHandler

Definition: MD.h:707

BaseParticle * getSmallestParticle() const

Gets a pointer to the smallest BaseParticle (by interactionRadius) in this ParticleHandler.

Definition: ParticleHandler.cc:95

CTangentialSpring * MD::getTangentialSpring	(	BaseParticle *	PI,
		BaseParticle *	PJ,
		BaseParticle *	PJreal
	)

protected

Definition at line 706 of file MD.cc.

References CTangentialSprings::create_new(), dt, get_dt(), BaseParticle::get_Index(), get_t(), TangentialSpringParticle::get_TangentialSprings(), CTangentialSprings::select_particle_spring(), and t.

Referenced by compute_internal_forces(), and compute_plastic_internal_forces().

 {
     //std::cout<<"i="<<PI->get_Index()<<" j="<<PJ->get_Index()<<std::endl;
     TangentialSpringParticle* TSParticleI=dynamic_cast<TangentialSpringParticle*>(PI);
     TangentialSpringParticle* TSParticleJ=dynamic_cast<TangentialSpringParticle*>(PJ);
     
     //First check if particle I has a spring stores that is connected to particle J
     CTangentialSpring* TangentialSpring = TSParticleI->get_TangentialSprings().select_particle_spring(PJreal->get_Index(), get_t(), get_dt());
     if(TangentialSpring==NULL)
     {
         //Then check if particle J has a spring stored that is connected to particle I
         TangentialSpring = TSParticleJ->get_TangentialSprings().select_particle_spring(PI->get_Index(), get_t(), get_dt());
         if(TangentialSpring==NULL)
         {
             //A new spring has to be created, this is done in the real particle with the highest index
             if(PI->get_Index()>PJreal->get_Index())
             {
                 TangentialSpring = TSParticleI->get_TangentialSprings().create_new(PJreal->get_Index(), t, dt);
                 //std::cout<<"Created new spring in Particle "<<PI->get_Index()<<" between particle "<<PI->get_Index()<<" and "<<PJreal->get_Index()<<std::endl<<std::endl;
             }
             else
             {
                 TangentialSpring = TSParticleJ->get_TangentialSprings().create_new(PI->get_Index(), t, dt);
                 //std::cout<<"Created new spring in Particle "<<PJ->get_Index()<<" between particle "<<PI->get_Index()<<" and "<<PJreal->get_Index()<<std::endl<<std::endl;
             }
         }
         //else
         //std::cout<<"i="<<PI->get_Index()<<" j="<<PJ->get_Index()<<" Reconnected to tangential spring in j"<<std::endl;
     }
     //else
     //std::cout<<"i="<<PI->get_Index()<<" j="<<PJ->get_Index()<<" Reconnected to tangential spring in i"<<std::endl;
     return TangentialSpring;
 }

CTangentialSpring * MD::getTangentialSpringWall	(	BaseParticle *	pI,
		int	w
	)

protected

Definition at line 741 of file MD.cc.

References CTangentialSprings::create_new_wall(), dt, get_dt(), get_t(), TangentialSpringParticle::get_TangentialSprings(), CTangentialSprings::select_wall_spring(), and t.

Referenced by compute_walls().

 {
     TangentialSpringParticle* TSParticle=dynamic_cast<TangentialSpringParticle*>(pI);
     CTangentialSpring* TangentialSpring = TSParticle->get_TangentialSprings().select_wall_spring(w, get_t(), get_dt());
 
     if(TangentialSpring==NULL)
         TangentialSpring = TSParticle->get_TangentialSprings().create_new_wall(w, t, dt);
 
     return TangentialSpring;
 }

WallHandler& MD::getWallHandler ( )

inline

Definition at line 148 of file MD.h.

References wallHandler.

Referenced by ChuteWithHopper::add_hopper(), Chute::create_bottom(), Chute::setup_particles_initial_conditions(), and ChuteBottom::setup_particles_initial_conditions().

148 { return wallHandler;}

MD::wallHandler

WallHandler wallHandler

Definition: MD.h:708

virtual void MD::HGRID_actions_after_integration ( )

inlineprotectedvirtual

Reimplemented in HGRID_base.

Definition at line 596 of file MD.h.

Referenced by solve().

596 {};

virtual void MD::HGRID_actions_before_integration ( )

inlineprotectedvirtual

Reimplemented in HGRID_base.

Definition at line 595 of file MD.h.

Referenced by solve().

595 {};

virtual void MD::HGRID_actions_before_time_loop ( )

inlineprotectedvirtual

This is actions before the start of the main time loop.

Reimplemented in HGRID_base.

Definition at line 552 of file MD.h.

Referenced by solve(), and statistics_from_restart_data().

552 {};

virtual void MD::HGRID_actions_before_time_step ( )

inlineprotectedvirtual

This is action before the time step is started.

Reimplemented in HGRID_base.

Definition at line 555 of file MD.h.

Referenced by solve(), and statistics_from_restart_data().

555 {};

virtual void MD::HGRID_InsertParticleToHgrid ( BaseParticle *obj UNUSED )

inlineprotectedvirtual

This is action before the time step is started.

Definition at line 558 of file MD.h.

558 {};

virtual void MD::HGRID_RemoveParticleFromHgrid ( BaseParticle *obj UNUSED )

inlineprotectedvirtual

Definition at line 560 of file MD.h.

Referenced by do_integration_after_force_computation(), and removeParticle().

560 {};

virtual void MD::HGRID_update_move	(	BaseParticle *	,
		Mdouble
	)

inlineprotectedvirtual

Reimplemented in HGRID_base.

Definition at line 594 of file MD.h.

Referenced by do_integration_before_force_computation().

594 {};

virtual void MD::HGRID_UpdateParticleInHgrid ( BaseParticle *obj UNUSED )

inlineprotectedvirtual

Definition at line 559 of file MD.h.

Referenced by do_integration_after_force_computation(), HGRID_base::HGRID_actions_before_time_step(), and HGRID_base::InitBroadPhase().

559 {};

void MD::info ( )

inlinevirtual

Set up a virtual info this will be provided from the inhertiance.

Reimplemented from STD_save.

Definition at line 93 of file MD.h.

References print().

93 {print(std::cout);}

MD::print

virtual void print(std::ostream &os, bool print_all=false)

Outputs MD.

Definition: MD.cc:1806

virtual void MD::InitBroadPhase ( )

inlineprotectedvirtual

Initialisation of Broad Phase Information (Default no Broad Phase so empty)

Reimplemented in HGRID_base.

Definition at line 602 of file MD.h.

602 {};

virtual void MD::initialize_statistics ( )

inlineprotectedvirtual

Functions for statistics.

Reimplemented in StatisticsVector< T >.

Definition at line 582 of file MD.h.

Referenced by solve(), and statistics_from_restart_data().

582 {};

void MD::initialize_tangential_springs ( )

protected

bool MD::load_from_data_file	(	const char *	filename,
		unsigned int	format = `0`
	)

This allows particle data to be reloaded from data files.

This function loads data from the .data file i.e.

you get position, mass and velocity information only See also MD::load_restart_data Can read in format 14 - 8 or format 7 data This code saves in format 8 for 2D and format 14 for 3D so if no extra parameters are specified it will assume this note: many parameters, like density cannot be set using the data file

Todo:: change from deprecated const char* to std::string

Definition at line 250 of file MD.cc.

References STD_save::data_file, STD_save::get_options_data(), STD_save::options_data, read_next_from_data_file(), and STD_save::set_options_data().

Referenced by readNextArgument(), and statistics_from_restart_data().

 {
     //This opens the file the data will be recalled from
     data_file.open( filename, std::fstream::in);
     if (!data_file.is_open()||data_file.bad()) {
         std::cout << "Loading data file " << filename << " failed" << std::endl;
         return false;
     } 
 
     //options_data is ignored
     int options_data = get_options_data();
     set_options_data(1);
     read_next_from_data_file(format);
     set_options_data(options_data);
     
     data_file.close();
 
     //std::cout << "Loaded data file " << filename << " (t=" << get_t() << ")" << std::endl;
     return true;
 }

bool MD::load_par_ini_file ( const char * filename )

allows the user to read par.ini files (useful to read MDCLR files)

Definition at line 271 of file MD.cc.

References boundaryHandler, BaseHandler< T >::copyAndAddObject(), WallHandler::copyAndAddWall(), get_dt(), get_xmax(), get_xmin(), get_ymax(), get_ymin(), get_zmax(), get_zmin(), PeriodicBoundary::set(), InfiniteWall::set(), set_dim_particle(), set_disp(), set_dt(), set_gravity(), set_k(), set_rho(), set_save_count_data(), set_save_count_fstat(), set_savecount(), set_t(), set_tmax(), and wallHandler.

 {
     //Opens the par.ini file
     std::fstream file;
     file.open( filename, std::fstream::in);
     if (!file.is_open()||file.bad()) {
         //std::cout << "Loading par.ini file " << filename << " failed" << std::endl;
         return false;
     } 
 
     Mdouble doubleValue;
     int integerValue;
 
     // inputfile par.ini 
     // line 1 =============================================================
     // Example: 1 1 0
     //   1: integer (0|1) switches from non-periodic to periodic
     //      integer (5|6) does 2D integration only (y-coordinates fixed)
     //                    and switches from non-periodic to periodic
     //      integer (11) uses a quarter system with circular b.c.
     file >> integerValue;
     //~ std::cout << "11" << integerValue << std::endl;
     if (integerValue==0) 
     {
         InfiniteWall w0;
         w0.set(Vec3D(-1,0,0), -get_xmin());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D( 1,0,0),  get_xmax());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D( 0,-1,0), -get_ymin());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D(0, 1,0),  get_ymax());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D(0,0,-1), -get_zmin());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D(0,0,1), get_zmax());
         wallHandler.copyAndAddWall(w0);
     } 
     else if (integerValue==1) 
     {
         PeriodicBoundary b0;
         b0.set(Vec3D(1,0,0), get_xmin(), get_xmax());
         boundaryHandler.copyAndAddObject(b0);
         b0.set(Vec3D(0,1,0), get_ymin(), get_ymax());
         boundaryHandler.copyAndAddObject(b0);
         b0.set(Vec3D(0,0,1), get_zmin(), get_zmax());
         boundaryHandler.copyAndAddObject(b0);
     }
     else if (integerValue==5) 
     {
         InfiniteWall w0;
         w0.set(Vec3D(-1,0,0), -get_xmin());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D( 1,0,0),  get_xmax());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D( 0,-1,0), -get_ymin());
         wallHandler.copyAndAddWall(w0);
         w0.set(Vec3D(0, 1,0),  get_ymax());
         wallHandler.copyAndAddWall(w0);
         
         
     } 
     else if (integerValue==6) 
     {
         PeriodicBoundary b0;
         b0.set(Vec3D(1,0,0), get_xmin(), get_xmax());
         boundaryHandler.copyAndAddObject(b0);
         b0.set(Vec3D(0,1,0), get_ymin(), get_ymax());
         boundaryHandler.copyAndAddObject(b0);
         b0.set(Vec3D(0,0,1), get_zmin(), get_zmax());
         boundaryHandler.copyAndAddObject(b0);
     } 
     else 
     {
         std::cerr << "Error in par.ini: line 1, value 1 is " << integerValue << std::endl;
         exit(-1);
     }
     
     //   2: integer (0|1) dont use | use the search pattern for linked cells
     file >> integerValue; //ignore
     
     //   3: real - gravity in z direction: positive points upwards
     file >> doubleValue;
     set_gravity(Vec3D(0.0,0.0,doubleValue));
 
     // line 2 =============================================================
     // Example: -10000 .5e-2
     //   1: time end of simulation - (negative resets start time to zero
     //                                and uses -time as end-time)   
     file >> doubleValue;
     if (doubleValue<0) set_t(0);
     set_tmax(fabs(doubleValue));
     
     //   2: time-step of simulation
     file >> doubleValue;
     set_dt(doubleValue);
 
     // line 3 =============================================================
     // Example: 1e-1 100
     file >> doubleValue;
     if (doubleValue>=0) {
         //   1: time-step for output on time-series protocoll file  -> "ene"
         int savecount = round(doubleValue/get_dt());
         set_savecount(savecount);
 
         //   2: time-step for output on film (coordinate) file      -> "c3d"
         //      (fstat-output is coupled to c3d-output time-step)
         file >> doubleValue;
         savecount = round(doubleValue/get_dt());
         set_save_count_data(savecount);
         set_save_count_fstat(savecount);
     } else {
         //  or: ---------------------------------------------------------------
         //   1: negative number is multiplied to the previous log-output-time
         //   2: requires initial log-output time
         //   3: negative number is multiplied to the previous film-output-time
         //   4: requires initial film-output time
         std::cerr << "Error in par.ini: line 3, value 1 is " << doubleValue << std::endl;
         exit(-1);
     }
 
     // line 4 =============================================================
     // Example: 2000 1e5 1e3 1e2
     //   1: particle density (mass=4/3*constants::pi*density*rad^3)
     file >> doubleValue;
     set_dim_particle(3);
     set_rho(doubleValue);
 
     //   2: linear spring constant
     file >> doubleValue;
     set_k(doubleValue);
 
     //   3: linear dashpot constant
     file >> doubleValue;
     set_disp(doubleValue);
 
     //   4: background damping dashpot constant
     file >> doubleValue;
     if (doubleValue) std::cerr << "Warning in par.ini: ignored background damping " << doubleValue << std::endl;
     
     // line 5 =============================================================
     // Example: 0 0
     //   1: growth rate:  d(radius) = xgrow * dt
     file >> doubleValue;
     if (doubleValue) std::cerr << "Warning in par.ini: ignored growth rate " << doubleValue << std::endl;
 
     //   2: target volume_fraction
     file >> doubleValue;
     if (doubleValue) std::cerr << "Warning in par.ini: ignored target volume_fraction " << doubleValue << std::endl;
 
     file.close();
     //std::cout << "Loaded par.ini file " << filename << std::endl;
     return true;
 }

int MD::load_restart_data ( )

Loads all MD data.

This function loads all MD data - See also MD::save_restart_data, MD::load_from_data_file This function return 1 if sucessful else it returns -1.

Definition at line 685 of file MD.cc.

References STD_save::problem_name.

Referenced by readNextArgument(), and statistics_from_restart_data().

                            {
     std::stringstream restart_filename;
     restart_filename.str("");
     restart_filename << problem_name.str().c_str()  << ".restart";
     return load_restart_data(restart_filename.str());
 }

int MD::load_restart_data ( std::string filename )

Definition at line 692 of file MD.cc.

References read(), and set_restarted().

                                              {
     std::ifstream restart_data(filename.c_str());
     if( restart_data.good() )
     {
         read(restart_data);
         restart_data.close();
         set_restarted(true);
         return(1);
     } else {
         std::cerr << "restart_data " << filename << " could not be loaded from " << filename << std::endl; 
         exit(-1);
     }
 }

void MD::output_ene ( )

protectedvirtual

This function outputs statistical data - The default is to compute the rotational kinetic engergy, linear kinetic energy, and the centre of mass.

todo{Why is there a +6 here?}

Definition at line 175 of file MD.cc.

References BaseHandler< T >::end(), STD_save::ene_file, get_ene_ela(), get_gravity(), get_t(), getParticleHandler(), and set_ene_ela().

Referenced by solve(), and statistics_from_restart_data().

 {
     Mdouble ene_kin = 0, ene_rot = 0, ene_gra = 0, mass_sum= 0, x_masslength=0, y_masslength=0, z_masslength=0;
 
     for (std::vector<BaseParticle*>::iterator it = getParticleHandler().begin(); it!=getParticleHandler().end(); ++it) if (!(*it)->isFixed())
     {
         ene_kin += .5 * (*it)->get_Mass() * (*it)->get_Velocity().GetLength2();
         ene_rot += .5 * (*it)->get_Inertia() * (*it)->get_AngularVelocity().GetLength2();
         ene_gra -= (*it)->get_Mass() * Dot(get_gravity(),(*it)->get_Position());
         mass_sum+= (*it)->get_Mass();
         x_masslength +=(*it)->get_Mass()*(*it)->get_Position().X;
         y_masslength +=(*it)->get_Mass()*(*it)->get_Position().Y;
         z_masslength +=(*it)->get_Mass()*(*it)->get_Position().Z;
     } //end for loop over Particles
 
     static int width = ene_file.precision() + 6;
     ene_file  << std::setw(width) << get_t() 
         << " " << std::setw(width) << ene_gra
         << " " << std::setw(width) << ene_kin
         << " " << std::setw(width) << ene_rot
         << " " << std::setw(width) << get_ene_ela()
         << " " << std::setw(width) << (mass_sum?x_masslength/mass_sum:std::numeric_limits<double>::quiet_NaN())
         << " " << std::setw(width) << (mass_sum?y_masslength/mass_sum:std::numeric_limits<double>::quiet_NaN()) 
         << " " << std::setw(width) << (mass_sum?z_masslength/mass_sum:std::numeric_limits<double>::quiet_NaN())
         << std::endl;
     
     set_ene_ela(0);
     //sliding = sticking = 0;
 }

virtual void MD::output_statistics ( )

inlineprotectedvirtual

Reimplemented in StatisticsVector< T >.

Definition at line 583 of file MD.h.

Referenced by solve().

583 {};

void MD::output_xballs_data ( )

protectedvirtual

Output xball data for Particle i.

This function outputs the location and velocity of the particle in a format the xballs progream can read.

Definition at line 209 of file MD.cc.

References STD_save::data_file, dim, format, BaseHandler< T >::getNumberOfObjects(), output_xballs_data_particle(), particleHandler, t, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by solve(), and StatisticsVector< T >::statistics_from_p3().

 {
     //Set the correct formation based of dimension if the formation is not specified by the user
     if (format ==0) {
             switch (dim) {
                 case 1:
                 case 2:
                     format=8;
                     break;
                 case 3:
                     format=14;
                     break;
             }
         }
 
     // This outputs the location of walls and how many particles there are to file this is required by the xballs plotting
     if (format!=14) // dim = 1 or 2
         data_file << particleHandler.getNumberOfObjects() << " " <<t << " " 
         << xmin << " " << ymin << " " 
         << xmax << " " << ymax << " " << std::endl;
     else //dim==3
         data_file << particleHandler.getNumberOfObjects() << " " <<t << " " 
         << xmin << " " << ymin << " " << zmin << " "
         << xmax << " " << ymax << " " << zmax << " " << std::endl;
                 // This outputs the particle data
     for (unsigned int i = 0;i<particleHandler.getNumberOfObjects();i++) 
         output_xballs_data_particle(i);
     #ifdef DEBUG_OUTPUT
         std::cerr << "Have output the properties of the problem to disk " << std::endl;
     #endif
 }

void MD::output_xballs_data_particle ( int i )

protectedvirtual

This function outputs the location and velocity of the particle in a format the xballs progream can read.

Todo:: {changes in *.icc files are not immediately regognized by the makefile!}

Definition at line 30 of file MD_xballs.icc.

References STD_save::data_file, BaseParticle::get_Angle(), BaseParticle::get_AngularVelocity(), get_dim(), get_format(), BaseParticle::get_Position(), BaseParticle::get_Radius(), BaseParticle::get_Velocity(), getInfo(), BaseHandler< T >::getObject(), getParticleHandler(), Vec3D::X, Vec3D::Y, and Vec3D::Z.

Referenced by output_xballs_data().

 {
     //data_file.precision(14);
     //This outputs the data about particle i again to the file.
     switch(get_format())
     {
         case 8:
         {
             if (get_dim()==1) {
                 data_file << getParticleHandler().getObject(i)->get_Position().X<<" 0 "<<getParticleHandler().getObject(i)->get_Velocity().X<<" 0 " <<getParticleHandler().getObject(i)->get_Radius()<<" 0 0 0"<<std::endl;
                 break;
             } else {
                 data_file 
                     << getParticleHandler().getObject(i)->get_Position().X << " " 
                     << getParticleHandler().getObject(i)->get_Position().Y << " " 
                     << getParticleHandler().getObject(i)->get_Velocity().X << " " 
                     << getParticleHandler().getObject(i)->get_Velocity().Y << " "
                     << getParticleHandler().getObject(i)->get_Radius() << " "
                     << -getParticleHandler().getObject(i)->get_Angle().Z << " " // negative b/c we are plotting (x,y) coordinates on the xz-axis of xballs
                     << -getParticleHandler().getObject(i)->get_AngularVelocity().Z << " "
                     << getInfo(*getParticleHandler().getObject(i)) <<std::endl;
             }
             break;
         }                           
         case 14:
         {
             data_file 
                     << getParticleHandler().getObject(i)->get_Position().X << " " 
                     << getParticleHandler().getObject(i)->get_Position().Y << " " 
                     << getParticleHandler().getObject(i)->get_Position().Z << " " 
                     << getParticleHandler().getObject(i)->get_Velocity().X << " " 
                     << getParticleHandler().getObject(i)->get_Velocity().Y << " " 
                     << getParticleHandler().getObject(i)->get_Velocity().Z << " "
                     << getParticleHandler().getObject(i)->get_Radius() << " "
                     << getParticleHandler().getObject(i)->get_Angle().X << " " // negative b/c we are plotting (x,y) coordinates on the xz-axis of xballs
                     << getParticleHandler().getObject(i)->get_Angle().Y << " " // negative b/c we are plotting (x,y) coordinates on the xz-axis of xballs
                     << getParticleHandler().getObject(i)->get_Angle().Z << " " // negative b/c we are plotting (x,y) coordinates on the xz-axis of xballs
                     << getParticleHandler().getObject(i)->get_AngularVelocity().X << " "
                     << getParticleHandler().getObject(i)->get_AngularVelocity().Y << " "
                     << getParticleHandler().getObject(i)->get_AngularVelocity().Z << " "
                     << getInfo(*getParticleHandler().getObject(i)) <<std::endl;
                 break;
         } //end case 3
         default:
         {
             std::cerr << "format not found" << std::endl;
         }
     } //end switch statement
 }

void MD::print	(	std::ostream &	os,
		bool	print_all = `false`
	)

virtual

Outputs MD.

Reimplemented in ChuteWithHopperAndInset, Chute, and HGRID_base.

Definition at line 1806 of file MD.cc.

References boundaryHandler, dim, dt, BaseHandler< T >::getNumberOfObjects(), WallHandler::getNumberOfWalls(), BaseHandler< T >::getObject(), BaseHandler< T >::getStorageCapacity(), WallHandler::getWall(), gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, BaseBoundary::print(), BaseParticle::print(), BaseWall::print(), STD_save::problem_name, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, wallHandler, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by info(), HGRID_base::print(), StatisticsVector< T >::statistics_from_fstat_and_data(), and StatisticsVector< T >::statistics_from_p3().

                                              {
         os << "problem_name:" << problem_name.str().c_str() << std::endl;
         os << " t:" << t << " dt:" << dt << ", tmax:" << tmax << ", save_count_data:" << save_count_data << ", save_count_ene:" << save_count_ene  << ", save_count_stat:" << save_count_stat << ", save_count_fstat:" << save_count_fstat << std::endl;
         os << " x:[" << xmin << "," << xmax << "], y:[" << ymin << "," << ymax << "], z:[" << zmin << "," << zmax << "]" << std::endl;
         os << " dim:" << dim << ", gravity:" << gravity << std::endl;
         if (Species.size()==1) {
             os << " "; Species[0].print(os); os << std::endl;
         } else {
             os << " Species: size:" << Species.size() << std::endl;
             for (unsigned int i=0; i<Species.size(); i++) { 
                 os << "  "; Species[i].print(os); os << std::endl; 
                 for (unsigned int j=0; j<Species[i].MixedSpecies.size(); j++) { os << "    "; Species[i].MixedSpecies[j].print(os); os << std::endl; }
             }
         }
         os << " options_fstat=" << options_fstat 
            << " , options_data=" << options_data 
            << " , options_ene=" << options_ene 
            << " , options_restart=" << options_restart << std::endl; 
         os << " Particles: N:" << particleHandler.getNumberOfObjects() << ", Nmax:" << particleHandler.getStorageCapacity() << std::endl;
         if (print_all) for (unsigned int i=0; i<particleHandler.getNumberOfObjects(); i++) { os << "  "; particleHandler.getObject(i)->print(os); os << std::endl; }
         os << " Walls: NWall:" << wallHandler.getNumberOfWalls() << ", NBoundary:" << boundaryHandler.getNumberOfObjects() << std::endl;
         for (unsigned int i=0; i<wallHandler.getNumberOfWalls(); i++) { os << "  "; wallHandler.getWall(i)->print(os); os << std::endl; }
         for (unsigned int i=0; i<boundaryHandler.getNumberOfObjects(); i++) { os << "  "; boundaryHandler.getObject(i)->print(os); os << std::endl; }
     }

virtual void MD::process_statistics ( bool usethese UNUSED )

inlineprotectedvirtual

Reimplemented in StatisticsVector< T >.

Definition at line 585 of file MD.h.

Referenced by solve().

585 {};

void MD::read ( std::istream & is )

virtual

Reads all MD data.

Reimplemented in ChuteWithHopper, ChuteWithHopperAndInset, Chute, and HGRID_base.

Definition at line 2026 of file MD.cc.

References boundaryHandler, WallHandler::clear(), ParticleHandler::clear(), BaseHandler< T >::clear(), BaseParticle::compute_particle_mass(), data_FixedParticles, dim, dt, BaseHandler< T >::getLastObject(), getLineFromStringStream(), gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, STD_save::problem_name, read_v1(), read_v2(), BoundaryHandler::readObject(), ParticleHandler::readObject(), WallHandler::readWall(), restart_version, save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, wallHandler, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by load_restart_data(), and HGRID_base::read().

 {
     std::string dummy;
     is >> dummy;
     if (strcmp(dummy.c_str(),"restart_version"))
     {
         dim = atoi(dummy.c_str());
         restart_version = 1;
         read_v1(is);
     }
     else
     {
         is  >> restart_version 
             >> dummy >> dummy;
         problem_name.str(dummy);
         if(restart_version==2)
         {
             read_v2(is);
         }
         else
         {
             //std::cout << "version 3" << std::endl;
             is  >> dummy >> xmin
                 >> dummy >> xmax
                 >> dummy >> ymin
                 >> dummy >> ymax
                 >> dummy >> zmin
                 >> dummy >> zmax
                 >> dummy >> dt
                 >> dummy >> t
                 >> dummy >> tmax
                 >> dummy >> save_count_data
                 >> dummy >> save_count_ene
                 >> dummy >> save_count_stat
                 >> dummy >> save_count_fstat
                 >> dummy >> dim
                 >> dummy >> gravity
                 >> dummy >> options_fstat
                 >> dummy >> options_data
                 >> dummy >> options_ene  
                 >> dummy;
                 //this is optional to allow restart files with and without options_restart  
             if (!strcmp(dummy.c_str(),"options_restart")) is >> options_restart >> dummy;
             unsigned int N; 
             is >> N; 
             Species.resize(0); //to clear the Species vector
             Species.resize(N);
             //please don't change this code to iterators; the mixed species requires indexing
             //old code: for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) it->read(is);
             for (unsigned int i=0; i<N; i++) {
                 Species[i].read(is);
                 Species[i].MixedSpecies.resize(i);
                 for (unsigned int j=0; j<i; j++) {
                     Species[i].MixedSpecies[j].read(is);
                 }
             }
             
             std::stringstream line (std::stringstream::in | std::stringstream::out);
             
             is >> dummy >> N; 
             wallHandler.clear();
             getLineFromStringStream(is,line);
             for (unsigned int i=0;i<N;i++)
             {
                 getLineFromStringStream(is,line);
                 wallHandler.readWall(line);
             }
             
             is >> dummy >> N;
             boundaryHandler.clear();
             getLineFromStringStream(is,line);
             for (unsigned int i=0;i<N;i++)
             {
                 getLineFromStringStream(is,line);
                 boundaryHandler.readObject(line);
             }
             
             is >> dummy >> N;
             particleHandler.clear();
             getLineFromStringStream(is,line);
             for (unsigned int i=0;i<N;i++)
             { 
                 getLineFromStringStream(is,line);
                 particleHandler.readObject(line);
                 particleHandler.getLastObject()->compute_particle_mass(Species);
             }
             
             //After the lines storing particle data, an optional variable 'FixedParticles' can be read in, which fixes the first 'FixedParticles' particles
             if (is.peek()==70) is >> dummy >> data_FixedParticles; 
         }
     }
 }

int MD::read_dim_from_data_file ( )

Definition at line 628 of file MD.cc.

References STD_save::data_file, STD_save::data_filename, t, xmax, xmin, ymax, ymin, zmax, and zmin.

 {
     //This opens the file the data will be recalled from; note that data_filename needs to be set
     std::fstream file;
     file.open( data_filename.str().c_str() , std::fstream::in);
     
     //read first line and close file
     std::string header_string;
     getline(data_file,header_string);
     std::stringstream header (std::stringstream::in | std::stringstream::out);
     header << header_string;
     file.close();
     
     //extract data from first line
     Mdouble N, t, xmin, ymin, zmin, xmax, ymax, zmax;
     header >> N >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
 
     //Set the correct formation based of dimension if the formation is not specified by the user
     if (header.eof()) return 2;
     else return 3;
 }

bool MD::read_next_from_data_file ( unsigned int format = 0 )

by default format do not pass an argument; only specify format if you have to read a special format (f.e. dim=2, but format=14 (3d format))

Todo:: {TW: what if mu=0, but get_Species(1).mu isn't? We have to check all mu, right?}

Todo:: {TW: also plastic particles have to be checked.}

Definition at line 453 of file MD.cc.

References BaseHandler< T >::begin(), BaseHandler< T >::copyAndAddObject(), STD_save::data_file, data_FixedParticles, dim, BaseHandler< T >::end(), BaseParticle::fixParticle(), get_k2max(), get_mu(), STD_save::get_options_data(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), STD_save::increase_counter_data(), particleHandler, BaseHandler< T >::removeLastObject(), reset_TangentialSprings(), Species, t, Vec3D::X, xmax, xmin, Vec3D::Y, ymax, ymin, Vec3D::Z, zmax, and zmin.

Referenced by load_from_data_file(), and StatisticsVector< T >::read_next_from_data_file().

 {
     if (get_options_data()==2) increase_counter_data(std::fstream::in);
 
 
 
     //Set the correct formation based of dimension if the formation is not specified by the user
     if (format ==0)
         {
             switch (dim)
             {
                 case 1:
                 case 2:
                     format=8;
                     break;
                 case 3:
                     format=14;
                     break;
             }
         //end case
             
         }
     //end if
     
     static unsigned int N;
     static Mdouble dummy;
     
     //read first parameter and check if we reached the end of the file
     data_file >> N;
     BaseParticle* P0;
     if(get_k2max()) {
         P0 = new DeltaMaxsParticle();
         //std::cout << "in read_next_from_data_file(): using DeltaMaxsParticle" << std::endl;
     } else if(get_mu()) {
         P0 = new TangentialSpringParticle();
         //std::cout << "in read_next_from_data_file(): using TangentialSpringParticle" << std::endl;
     } else {
         P0 = new BaseParticle();
         //std::cout << "in read_next_from_data_file(): using BaseParticle" << std::endl;
     }
     
     if(particleHandler.getNumberOfObjects()<N)
         while(particleHandler.getNumberOfObjects()<N)
             particleHandler.copyAndAddObject(P0);
     else
         while(particleHandler.getNumberOfObjects()>N)
             particleHandler.removeLastObject();
     delete P0;
     
     #ifdef DEBUG_OUTPUT
         std::cout << "Number of particles read from file "<<N << std::endl;
     #endif
         
     if (data_file.eof()||data_file.peek()==-1) return false;
     //read all other data available for the time step
     switch(format)
     {
         //This is the format that Stefan's and Vitaley's code saves in - note the axis rotation
         case 7:
         {
             data_file >> t >> xmin >> zmin >> ymin >> xmax >> zmax >> ymax;
             //std::cout << " time " << t <<  std::endl;
             Mdouble radius;
             Vec3D position,velocity;
             for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
                 data_file
                     >> position.X
                     >> position.Z
                     >> position.Y
                     >> velocity.X
                     >> velocity.Z
                     >> velocity.Y
                     >> radius 
                     >> dummy;
                     (*it)->set_Position(position);
                     (*it)->set_Velocity(velocity);
                     (*it)->set_Angle(Vec3D(0.0,0.0,0.0));
                     (*it)->set_AngularVelocity(Vec3D(0.0,0.0,0.0));
                     (*it)->set_Radius(radius);
                     (*it)->compute_particle_mass(Species);
                 }
             //End the interator over all particles
             break;
         }
         //This is a 2D format
         case 8:
         {
             data_file >> t >> xmin >> ymin >> xmax >> ymax;
             zmin = 0.0;
             zmax = 0.0;
             Mdouble radius;
             Vec3D position,velocity,angle,angularVelocity;
             for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
                 data_file 
                     >> position.X 
                     >> position.Y 
                     >> velocity.X 
                     >> velocity.Y 
                     >> radius 
                     >> angle.Z
                     >> angularVelocity.Z 
                     >> dummy;
                 (*it)->set_Position(position);
                 (*it)->set_Velocity(velocity);
                 (*it)->set_Angle(-angle);
                 (*it)->set_AngularVelocity(-angularVelocity);
                 (*it)->set_Radius(radius);
                 (*it)->compute_particle_mass(Species);
             } //end for all particles
             break;
         }
         //This is a 3D format
         case 14:
         {   
             data_file >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
             Mdouble radius;
             Vec3D position,velocity,angle,angularVelocity;
             for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
                 data_file 
                     >> position
                     >> velocity
                     >> radius 
                     >> angle
                     >> angularVelocity
                     >> dummy;
                 (*it)->set_Position(position);
                 (*it)->set_Velocity(velocity);
                 (*it)->set_Angle(angle);
                 (*it)->set_AngularVelocity(angularVelocity);                    
                 (*it)->set_Radius(radius);
                 (*it)->compute_particle_mass(Species);
             } //end for all particles
             break;
         }
         //This is a 3D format
         case 15:
         {
             data_file >> t >> xmin >> ymin >> zmin >> xmax >> ymax >> zmax;
             Mdouble radius;
             Vec3D position,velocity,angle,angularVelocity;            
             for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it) {
                 data_file 
                     >> position
                     >> velocity
                     >> radius 
                     >> angle
                     >> angularVelocity
                     >> dummy >> dummy;
                 (*it)->set_Position(position);
                 (*it)->set_Velocity(velocity);
                 (*it)->set_Angle(angle);
                 (*it)->set_AngularVelocity(angularVelocity);                    
                 (*it)->set_Radius(radius);
                 (*it)->compute_particle_mass(Species);
             } //end for all particles
             break;
         }   
     } //end switch statement
     
     //fix particles, if data_FixedParticles!=0
     if (data_FixedParticles) {
         //std::cout << "fix " << min(data_FixedParticles,getParticleHandler().getNumberOfObjects()) << " Particles" << std::endl;
         for (int i=0; i<std::min(data_FixedParticles,(int) particleHandler.getNumberOfObjects()); i++)
             particleHandler.getObject(i)->fixParticle();
     }
         
     //clean up tangential springs
     reset_TangentialSprings();
     
     return true;
 }

void MD::read_v1 ( std::istream & is )

virtual

Reads all MD data.

Definition at line 2217 of file MD.cc.

References boundaryHandler, BaseParticle::compute_particle_mass(), dt, BaseHandler< T >::getLastObject(), getLineFromStringStream(), gravity, max_radius, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, particleHandler, STD_save::problem_name, BoundaryHandler::readObject(), ParticleHandler::readObject(), WallHandler::readWall(), save_count_data, save_count_ene, save_count_fstat, save_count_stat, WallHandler::setStorageCapacity(), BaseHandler< T >::setStorageCapacity(), Species, t, tmax, wallHandler, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by read().

 {
     std::cout << "reading restart data version 1" << std::endl;
     is  >> gravity 
         >> xmin >> xmax >> ymin >> ymax >> zmin >> zmax 
         >> dt >> t >> tmax 
         >> save_count_data >> max_radius;
     save_count_ene=save_count_data;
     save_count_fstat=save_count_data;
     save_count_stat=save_count_data;            
     std::string str;
     is >> str;
     problem_name.str(str);
     int N;
     unsigned int NParticles, NWalls, NBoundarys;
     is >> N; Species.resize(N);
     is >> NParticles;           particleHandler.setStorageCapacity(NParticles);
     is >> NWalls;               wallHandler.setStorageCapacity(NWalls);
     is >> NBoundarys;           boundaryHandler.setStorageCapacity(NBoundarys);
     is >> options_fstat >> options_data >> options_ene; 
     for (std::vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) is >> (*it);
     
     std::stringstream line (std::stringstream::in | std::stringstream::out);
     getLineFromStringStream(is,line);
     
     for (unsigned int i=0;i<NParticles;i++)
     {
         getLineFromStringStream(is,line);
         particleHandler.readObject(line);
         particleHandler.getLastObject()->compute_particle_mass(Species);
     }
     
     for (unsigned int i=0;i<NWalls;i++)
     {
         getLineFromStringStream(is,line);
         wallHandler.readWall(line);
     }
     
     for (unsigned int i=0;i<NBoundarys;i++)
     {
         getLineFromStringStream(is,line);
         boundaryHandler.readObject(line);
     }
 }

void MD::read_v2 ( std::istream & is )

virtual

Definition at line 2121 of file MD.cc.

References boundaryHandler, WallHandler::clear(), ParticleHandler::clear(), BaseHandler< T >::clear(), BaseParticle::compute_particle_mass(), data_FixedParticles, dim, dt, BaseHandler< T >::getLastObject(), getLineFromStringStream(), gravity, STD_save::options_data, STD_save::options_ene, STD_save::options_fstat, STD_save::options_restart, particleHandler, BoundaryHandler::readObject(), ParticleHandler::readObject(), WallHandler::readWall(), save_count_data, save_count_ene, save_count_fstat, save_count_stat, Species, t, tmax, wallHandler, xmax, xmin, ymax, ymin, zmax, and zmin.

Referenced by read().

 {
     std::string dummy;
     //std::cout << "version 2" << std::endl;
     is  >> dummy >> xmin
         >> dummy >> xmax
         >> dummy >> ymin
         >> dummy >> ymax
         >> dummy >> zmin
         >> dummy >> zmax
         >> dummy >> dt
         >> dummy >> t
         >> dummy >> tmax
         >> dummy >> save_count_data
         >> dummy >> dim
         >> dummy >> gravity
         >> dummy >> options_fstat
         >> dummy >> options_data
         >> dummy >> options_ene
         >> dummy;
     save_count_ene=save_count_data;
     save_count_fstat=save_count_data;
     save_count_stat=save_count_data;
 
     //this is optional to allow restart files with and without options_restart  
     if (!strcmp(dummy.c_str(),"options_restart")) is >> options_restart >> dummy;
 
     unsigned int N; 
 
     is >> N; 
     Species.resize(N);
     //please don't change this code to iterators; the mixed species requires indexing
     //old code: for (vector<CSpecies>::iterator it = Species.begin(); it!=Species.end(); ++it) it->read(is);
     for (unsigned int i=0; i<N; i++) {
         Species[i].read(is);
         Species[i].MixedSpecies.resize(i);
         for (unsigned int j=0; j<i; j++) {
             Species[i].MixedSpecies[j].read(is);
         }
     }
 
     std::stringstream line (std::stringstream::in | std::stringstream::out);
     
     is >> dummy >> N; 
     wallHandler.clear();
     getLineFromStringStream(is,line);   
     for (unsigned int i=0;i<N;i++)
     {
         getLineFromStringStream(is,line);
         wallHandler.readWall(line);
     }
 
     is >> dummy >> N;
     boundaryHandler.clear();
     getLineFromStringStream(is,line);   
     for (unsigned int i=0;i<N;i++)
     {
         getLineFromStringStream(is,line);
         boundaryHandler.readObject(line);   
     }
 
     is >> dummy >> N;
     particleHandler.clear();
     getLineFromStringStream(is,line);   
     for (unsigned int i=0;i<N;i++)
     { 
         getLineFromStringStream(is,line);
         particleHandler.readObject(line);
         particleHandler.getLastObject()->compute_particle_mass(Species);
     }
     
     //After the lines storing particle data, an optional variable 'FixedParticles' can be read in, which fixes the first 'FixedParticles' particles
     if (is.peek()==70) is >> dummy >> data_FixedParticles; 
 }

int MD::readArguments	(	unsigned int	argc,
		char *	argv[]
	)

Can interpret main function input arguments that are passed by the driver codes.

Definition at line 2535 of file MD.cc.

References readNextArgument().

Referenced by solve().

 {
     bool isRead = true;
     for (unsigned int i = 1; i<argc; i+=2) {
         std::cout << "interpreting input argument " << argv[i]; 
         for (unsigned int j = i+1; j<argc; j++) {
             if (argv[j][0]=='-') break;
             std::cout << " " << argv[j];
         }
         std::cout << std::endl;
         isRead &= readNextArgument(i, argc, argv);
         if (!isRead) {
             std::cerr << "Warning: not all arguments read correctly!" << std::endl;
             exit(-10);
         }
     }   
     return isRead;
 }

int MD::readNextArgument	(	unsigned int &	i,
		unsigned int	argc,
		char *	argv[]
	)

virtual

argv[i+1] interpreted as argument of type char*, Mdouble, integer or boolean unless noted

-gravity requires three arguments

-restart or -r loads a restart file. By default, it loads <get_name()>.restart. If an argument "arg" is given it loads the file "arg", or "arg".restart (if the ending is not given).

Reimplemented in ChuteWithHopper, ChuteWithHopperAndInset, Chute, and HGRID_base.

Definition at line 2554 of file MD.cc.

References STD_save::auto_number(), get_disp(), get_disprolling(), get_dispt(), get_dt(), get_k(), get_krolling(), get_kt(), get_mu(), get_murolling(), STD_save::get_name(), get_tmax(), Hertzian, load_from_data_file(), load_restart_data(), random, RNG::randomise(), set_append(), STD_save::set_counter(), set_dim(), set_dim_particle(), set_disp(), set_disprolling(), set_dispt(), set_dt(), set_FixedParticles(), set_gravity(), set_k(), set_krolling(), set_kt(), set_mu(), set_murolling(), set_name(), set_number_of_saves_all(), STD_save::set_options_data(), STD_save::set_options_ene(), STD_save::set_options_fstat(), STD_save::set_options_restart(), STD_save::set_options_stat(), set_rho(), set_save_count_data(), set_save_count_ene(), set_save_count_fstat(), set_save_count_stat(), set_savecount(), set_tmax(), set_xmax(), set_xmin(), set_ymax(), set_ymin(), set_zmax(), set_zmin(), and Species.

Referenced by readArguments(), and HGRID_base::readNextArgument().

 {
     if (!strcmp(argv[i],"-name")) {
         set_name(argv[i+1]);
     } else if (!strcmp(argv[i],"-xmin")) {
         set_xmin(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-ymin")) {
         set_ymin(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-zmin")) {
         set_zmin(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-xmax")) {
         set_xmax(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-ymax")) {
         set_ymax(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-zmax")) {
         set_zmax(atof(argv[i+1]));
     //} else if (!strcmp(argv[i],"-svn")) {
     //  std::cout << "svn version " << SVN_VERSION << std::endl;
     //  i--;
     } else if (!strcmp(argv[i],"-dt")) { 
         Mdouble old=get_dt();
         set_dt(atof(argv[i+1]));
         std::cout << "  reset dt from " << old << " to " << get_dt() << std::endl;
     } else if (!strcmp(argv[i],"-Hertzian")) { 
         Species[0].set_ForceType(Hertzian);
         i--;
     } else if (!strcmp(argv[i],"-tmax")) {
         Mdouble old = get_tmax();
         set_tmax(atof(argv[i+1]));
         std::cout << "  reset tmax from " << old << " to " << get_tmax() << std::endl;
     } else if (!strcmp(argv[i],"-save_count") || !strcmp(argv[i],"-savecount")) {
         set_savecount(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-save_count_data")) {
         set_save_count_data(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-save_count_fstat")) {
         set_save_count_fstat(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-save_count_stat")) {
         set_save_count_stat(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-save_count_ene")) {
         set_save_count_ene(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-dim")) {
         set_dim(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-gravity")) {
         set_gravity(Vec3D(atof(argv[i+1]),atof(argv[i+2]),atof(argv[i+3]))); 
         i+=2; 
     } else if (!strcmp(argv[i],"-options_fstat")) {
         set_options_fstat(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-options_restart")) {
         set_options_restart(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-options_data")) {
         set_options_data(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-options_stat")) {
         set_options_stat(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-options_ene")) {
         set_options_ene(atoi(argv[i+1]));   
     } else if (!strcmp(argv[i],"-auto_number")) {
         auto_number();  
         i--;
     } else if (!strcmp(argv[i],"-number_of_saves")) {
         set_number_of_saves_all(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-restart")||!strcmp(argv[i],"-r")) {
         std::string filename;
         
         //use default filename if no argument is given
         if (i+1>=argc||argv[i+1][0]=='-') {
             i--;
             filename = get_name();
             std::cout << get_name() << std::endl;
         } else filename = argv[i+1];
                 
         //add ".restart" if necessary
         const char* fileextension = (filename.c_str()+std::max(0,(int)filename.length()-8));
         if (strcmp(fileextension,".restart")) 
             filename=filename+".restart";
         
         std::cout << "restart from " << filename << std::endl;
         load_restart_data(filename);
     } else if (!strcmp(argv[i],"-data")) {
         std::string filename = argv[i+1];
         load_from_data_file(filename.c_str());
     } else if (!strcmp(argv[i],"-k")) {
         Mdouble old = get_k();
         set_k(atof(argv[i+1]));
         std::cout << "  reset k from " << old << " to " << get_k() << std::endl;
     } else if (!strcmp(argv[i],"-disp")) {
         Mdouble old = get_disp();
         set_disp(atof(argv[i+1]));
         std::cout << "  reset disp from " << old << " to " << get_disp() << std::endl;
     } else if (!strcmp(argv[i],"-kt")) {
         Mdouble old = get_kt();
         set_kt(atof(argv[i+1]));
         std::cout << "  reset kt from " << old << " to " << get_kt() << std::endl;
     } else if (!strcmp(argv[i],"-dispt")) {
         Mdouble old = get_dispt();
         set_dispt(atof(argv[i+1]));
         std::cout << "  reset dispt from " << old << " to " << get_dispt() << std::endl;
     } else if (!strcmp(argv[i],"-krolling")) {
         Mdouble old = get_krolling();
         set_krolling(atof(argv[i+1]));
         std::cout << "  reset krolling from " << old << " to " << get_krolling() << std::endl;
     } else if (!strcmp(argv[i],"-disprolling")) {
         Mdouble old = get_disprolling();
         set_disprolling(atof(argv[i+1]));
         std::cout << "  reset disprolling from " << old << " to " << get_disprolling() << std::endl;
     } else if (!strcmp(argv[i],"-mu")) {
         Mdouble old = get_mu();
         set_mu(atof(argv[i+1]));
         std::cout << "  reset mu from " << old << " to " << get_mu() << std::endl;
     } else if (!strcmp(argv[i],"-murolling")) {
         Mdouble old = get_mu();
         set_murolling(atof(argv[i+1]));
         std::cout << "  reset murolling from " << old << " to " << get_murolling() << std::endl;
     } else if (!strcmp(argv[i],"-randomise") || !strcmp(argv[i],"-randomize")) {
         random.randomise(); 
         i--;
     } else if (!strcmp(argv[i],"-k0")) {
         Mdouble old = Species[0].get_k0();
         Species[0].set_k0(atof(argv[i+1]));
         std::cout << "  reset k0 from " << old << " to " << Species[0].get_k0() << std::endl;
     } else if (!strcmp(argv[i],"-f0")) {
         Mdouble old = Species[0].get_f0();
         Species[0].set_f0(atof(argv[i+1]));
         std::cout << "  reset f0 from " << old << " to " << Species[0].get_f0() << std::endl;
     } else if (!strcmp(argv[i],"-AdhesionForceType")) {
         Mdouble old = Species[0].get_AdhesionForceType();
         Species[0].set_AdhesionForceType(argv[i+1]);
         std::cout << "  reset AdhesionForceType from " << old << " to " << Species[0].get_AdhesionForceType() << std::endl;
     } else if (!strcmp(argv[i],"-append")) {
         set_append(true); 
         i--;
     } else if (!strcmp(argv[i],"-fixedParticles")) {
         set_FixedParticles(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-rho")) {
         set_rho(atof(argv[i+1]));
     } else if (!strcmp(argv[i],"-dim_particle")) {
         set_dim_particle(atoi(argv[i+1]));
     } else if (!strcmp(argv[i],"-counter")) {
         set_counter(atoi(argv[i+1]));
     } else return false;
     return true; //returns true if argv is found
 }

void MD::Remove_Duplicate_Periodic_Particles ( )

private

Remove periodic duplicate Particles (i.e. removes particles created by Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)). Note that between these two functions it is not allowed to create additional functions.

Remove periodic duplicate Particles (i.e. removes particles created by Check_and_Duplicate_Periodic_Particle(int i, int nWallPeriodic)). Note that between these two functions it is not allowed to create additional functions. It returns the number of Particles removed.

Definition at line 2703 of file MD.cc.

References BaseParticle::get_PeriodicFromParticle(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), particleHandler, PeriodicCreated, R, and removeParticle().

Referenced by solve(), and statistics_from_restart_data().

     {
         int R=0;
         for(int i=particleHandler.getNumberOfObjects()-1; i>=0&&particleHandler.getObject(i)->get_PeriodicFromParticle()!=NULL; i--)
         {
             removeParticle(i);
             R++;
         }
         if(R!=PeriodicCreated)
         {
             std::cerr<<"ERROR :: Periodic Particles Removed not equal to Periodic Particles Created, it created "<<PeriodicCreated<<" Particles, but removed "<<R<<std::endl<<std::endl;
             exit(-1);
         }
     }

void MD::Remove_Particle ( int IP )

This function removes partice IP from the vector of particles by moving the last particle in the vector to the position if IP Also it checks if the moved Particle has any tangentialsspring-information, which needs to be moved to a different particle, because tangential spring information always needs to be stored in the real particle with highest particle index.

virtual void MD::removeParticle ( int iP )

inlinevirtual

Definition at line 420 of file MD.h.

References BaseHandler< T >::getObject(), HGRID_RemoveParticleFromHgrid(), particleHandler, and BaseHandler< T >::removeObject().

Referenced by Chute::clean_chute(), and Remove_Duplicate_Periodic_Particles().

     {
 #ifdef ContactListHgrid
         possibleContactList.remove_ParticlePosibleContacts(particleHandler.getObject(iP));
 #endif
         HGRID_RemoveParticleFromHgrid(particleHandler.getObject(iP));
         particleHandler.removeObject(iP);
     }

void MD::reset_DeltaMax ( )

inlineprotected

sets the history parameter DeltaMax of all particles to zero

Definition at line 633 of file MD.h.

References BaseHandler< T >::begin(), BaseHandler< T >::end(), DeltaMaxsParticle::get_DeltaMaxs(), and particleHandler.

     {
         for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
         {
             DeltaMaxsParticle* DMParticle=dynamic_cast<DeltaMaxsParticle*>(*it);
             if(DMParticle)
                 DMParticle->get_DeltaMaxs().resize(0);
         }
     }

void MD::reset_TangentialSprings ( )

inlineprotected

sets the history parameter TangentialSprings of all particles to zero

Definition at line 644 of file MD.h.

References BaseHandler< T >::begin(), BaseHandler< T >::end(), TangentialSpringParticle::get_TangentialSprings(), and particleHandler.

Referenced by read_next_from_data_file().

     {
         for (std::vector<BaseParticle*>::iterator it = particleHandler.begin(); it!=particleHandler.end(); ++it)
         {
             TangentialSpringParticle* TSParticle=dynamic_cast<TangentialSpringParticle*>(*it);
             if(TSParticle)
                 TSParticle->get_TangentialSprings().resize(0);
         }
     }

void MD::save_restart_data ( )

virtual

Stores all MD data.

This function stores all MD data - See also MD::load_restart_data.

todo{This can probably done a lot nicer}

Definition at line 654 of file MD.cc.

References STD_save::options_restart, STD_save::problem_name, save_restart_data_counter, and write().

Referenced by solve(), and StatisticsVector< T >::statistics_from_p3().

                             {   
     std::stringstream restart_filename;
     restart_filename.str("");
     restart_filename << problem_name.str().c_str()  << ".restart";
     
     if (options_restart==2) {
         restart_filename << ".";
         if (save_restart_data_counter<1000) restart_filename << "0";
         if (save_restart_data_counter<100) restart_filename << "0";
         if (save_restart_data_counter<10) restart_filename << "0";
         restart_filename << save_restart_data_counter;
         save_restart_data_counter++;
     }
     
     std::ofstream restart_data(restart_filename.str().c_str());
     restart_data.precision(8); //max. 16 for Mdouble precision
     if( restart_data.good() )
     {
         write(restart_data);
         restart_data.close();
     } else {
         std::cerr << "restart_data " << restart_filename.str() << " could not be saved" << std::endl; 
         exit(-1);
     }
 }

void MD::set_append ( bool new_ )

inline

Sets restarted.

Definition at line 391 of file MD.h.

References append.

Referenced by constructor(), readNextArgument(), and set_restarted().

391 {append=new_;}

MD::append

bool append

Definition: MD.h:720

void MD::set_collision_time_and_normal_and_tangential_restitution_coefficient	(	Mdouble	tc,
		Mdouble	eps,
		Mdouble	beta,
		Mdouble	mass1,
		Mdouble	mass2,
		unsigned int	indSpecies = `0`
	)

inline

See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient.

Definition at line 289 of file MD.h.

References Species.

                                                                                                                                                                               {
         if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_normal_and_tangential_restitution_coefficient(tc, eps, beta, mass1, mass2);} else {std::cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
     }

void MD::set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt	(	Mdouble	tc,
		Mdouble	eps,
		Mdouble	beta,
		Mdouble	mass1,
		Mdouble	mass2,
		unsigned int	indSpecies = `0`
	)

inline

See CSpecies::set_collision_time_and_normal_and_tangential_restitution_coefficient.

Definition at line 294 of file MD.h.

References Species.

                                                                                                                                                                                        {
         if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_normal_and_tangential_restitution_coefficient_nodispt(tc, eps, beta, mass1, mass2);} else {std::cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
 
     }

void MD::set_collision_time_and_restitution_coefficient	(	Mdouble	tc,
		Mdouble	eps,
		Mdouble	mass,
		unsigned int	indSpecies = `0`
	)

inline

Sets k, disp such that it matches a given tc and eps for a collision of two copies of P.

Definition at line 276 of file MD.h.

References Species.

                                                                                                                            {
         if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_restitution_coefficient(tc, eps, mass);} else {std::cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
     }

void MD::set_collision_time_and_restitution_coefficient	(	Mdouble	tc,
		Mdouble	eps,
		Mdouble	mass1,
		Mdouble	mass2,
		unsigned int	indSpecies = `0`
	)

inline

Set k, disp such that is matches a given tc and eps for a collision of two different masses.

Recall the resitution constant is a function of k, disp and the mass of each particle in the collision See also set_collision_time_and_restitution_coefficient(Mdouble tc, Mdouble eps, Mdouble mass)

Definition at line 283 of file MD.h.

References Species.

         {
                 if (indSpecies<Species.size()) {Species[indSpecies].set_collision_time_and_restitution_coefficient(tc, eps, mass1, mass2);} else {std::cerr << "Error in set_collision_time_and_restitution_coefficient: species does not exist"; exit(-1);}
         }

void MD::set_depth	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Definition at line 189 of file MD.h.

References Species.

190 {if (indSpecies<Species.size()) {Species[indSpecies].set_depth(new_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_dim ( int new_dim )

inline

Allows the dimension of the simulation to be changed.

Definition at line 367 of file MD.h.

References dim.

Referenced by HGRID_3D::constructor(), readNextArgument(), and StatisticsVector< T >::statistics_from_p3().

367 { if (new_dim>=1 && new_dim<=3) dim = new_dim; else { std::cerr << "Error in set_dim" << std::endl; exit(-1); }}

MD::dim

int dim

The dimension of the simulation.

Definition: MD.h:660

void MD::set_dim_particle	(	int	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the dimension of the particle (f.e. for mass) to be changed.

Definition at line 261 of file MD.h.

References Species.

Referenced by HGRID_3D::constructor(), constructor(), load_par_ini_file(), and readNextArgument().

261 {if (indSpecies<Species.size()) {Species[indSpecies].set_dim_particle(new_);} else {std::cerr << "Error in set_dim_particle: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_disp	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the normal dissipation to be changed.

Definition at line 237 of file MD.h.

References Species.

Referenced by load_par_ini_file(), and readNextArgument().

237 {if (indSpecies<Species.size()) {Species[indSpecies].set_dissipation(new_);} else {std::cerr << "Error in set_dissipation: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_disprolling	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the tangential viscosity to be changed.

Definition at line 228 of file MD.h.

References Species.

Referenced by readNextArgument().

228 {if (indSpecies<Species.size()) {Species[indSpecies].set_disprolling(new_);} else {std::cerr << "Error in set_disprolling: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_dispt	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the tangential viscosity to be changed.

Definition at line 224 of file MD.h.

References Species.

Referenced by readNextArgument().

224 {if (indSpecies<Species.size()) {Species[indSpecies].set_dispt(new_);} else {std::cerr << "Error in set_dispt: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_disptorsion	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the tangential viscosity to be changed.

Definition at line 232 of file MD.h.

References Species.

232 {if (indSpecies<Species.size()) {Species[indSpecies].set_disptorsion(new_);} else {std::cerr << "Error in set_disptorsion: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_dissipation	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the normal dissipation to be changed.

Definition at line 241 of file MD.h.

References Species.

Referenced by Chute::set_collision_time_and_restitution_coefficient().

241 {if (indSpecies<Species.size()) {Species[indSpecies].set_dissipation(new_);} else {std::cerr << "Error in set_dissipation: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_do_stat_always ( bool new_ )

inline

Sets how often the data is saved using the number of saves wanted, tmax, and dt. See also set_savecount.

Definition at line 170 of file MD.h.

References do_stat_always.

Referenced by constructor().

170 {do_stat_always=new_;}

MD::do_stat_always

bool do_stat_always

Definition: MD.h:683

void MD::set_dt ( Mdouble new_dt )

inline

Allows the time step dt to be changed.

Definition at line 337 of file MD.h.

References dt.

337 {if (dt>=0.0){dt=new_dt;} else { std::cerr << "Error in set_dt" << std::endl; exit(-1); }}

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

void MD::set_dt ( BaseParticle & P )

inline

Sets dt to 1/50-th of the collision time for two copies of P.

Definition at line 440 of file MD.h.

References dt, get_collision_time(), BaseParticle::get_IndSpecies(), BaseParticle::get_Mass(), and Species.

                                 {
         Mdouble mass = P.get_Mass(); 
         dt = .02 * get_collision_time(mass);
         if (Species[P.get_IndSpecies()].dispt) { dt = std::min(dt,0.125*mass/Species[P.get_IndSpecies()].dispt); std::cerr << "Warning: dispt is large, dt had to be lowered"; }
     }

void MD::set_dt ( )

inline

Sets dt to 1/50-th of the smallest possible collision time.

Todo:: should setup_... be used here or not (by Thomas)

Todo:: {Is it necesarry to reset initial conditions here and in solve? (i.e. should it be in constructor)?}

Definition at line 447 of file MD.h.

References BaseParticle::compute_particle_mass(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), getSmallestParticle(), particleHandler, setup_particles_initial_conditions(), and Species.

Referenced by actions_before_time_loop(), load_par_ini_file(), readNextArgument(), and Chute::set_dt().

                  {
         setup_particles_initial_conditions();
         for (unsigned int i=0;i<particleHandler.getNumberOfObjects();i++) particleHandler.getObject(i)->compute_particle_mass(Species);
         set_dt(getSmallestParticle()[0]);
     }

void MD::set_dt_by_mass ( Mdouble mass )

inline

Sets dt to 1/50-th of the collision time for two particles of mass P.

Definition at line 437 of file MD.h.

References dt, and get_collision_time().

437 {dt = .02 * get_collision_time(mass);}

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

MD::get_collision_time

Mdouble get_collision_time(Mdouble mass, unsigned int indSpecies=0)

Calculates collision time for two copies of a particle of given disp, k, mass.

Definition: MD.h:300

void MD::set_ene_ela ( Mdouble new_ )

inline

Sets ene_ela.

Definition at line 397 of file MD.h.

References ene_ela.

Referenced by output_ene().

397 {ene_ela=new_;}

MD::ene_ela

Mdouble ene_ela

used in force calculations

Definition: MD.h:690

void MD::set_FixedParticles ( unsigned int n )

inlineprotected

Definition at line 613 of file MD.h.

References BaseParticle::fixParticle(), BaseHandler< T >::getNumberOfObjects(), BaseHandler< T >::getObject(), and particleHandler.

Referenced by readNextArgument().

613 {for (unsigned int i=0; i<std::min((unsigned int)particleHandler.getNumberOfObjects(),n); i++) particleHandler.getObject(i)->fixParticle();}

BaseHandler::getObject

T * getObject(const unsigned int id) const

Gets a pointer to the Object at the specified index in the BaseHandler.

Definition: BaseHandler.h:176

BaseHandler::getNumberOfObjects

ParticleHandler particleHandler

Definition: MD.h:707

unsigned int getNumberOfObjects() const

Gets the number of Object in this BaseHandler.

Definition: BaseHandler.h:199

BaseParticle::fixParticle

void fixParticle()

Fix Particle function. It fixes a Particle by setting its inverse mass and inertia and velocities to ...

Definition: BaseParticle.cc:142

void MD::set_format ( int new_ )

inline

Definition at line 507 of file MD.h.

References format.

Referenced by StatisticsVector< T >::statistics_from_p3().

507 {format = new_;}

MD::format

int format

Definition: MD.h:699

void MD::set_gravity ( Vec3D new_gravity )

inline

Allows the gravitational acceleration to be changed.

Definition at line 362 of file MD.h.

References gravity.

Referenced by load_par_ini_file(), readNextArgument(), and Chute::set_ChuteAngle().

362 {gravity = new_gravity;}

MD::gravity

Vec3D gravity

Gravitational acceleration.

Definition: MD.h:663

virtual void MD::set_initial_pressures_for_pressure_controlled_walls ( )

inlineprotectedvirtual

Definition at line 587 of file MD.h.

Referenced by solve().

587 {};

void MD::set_k	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the spring constant to be changed.

Definition at line 204 of file MD.h.

References Species.

Referenced by constructor(), load_par_ini_file(), readNextArgument(), and Chute::set_collision_time_and_restitution_coefficient().

204 {if (indSpecies<Species.size()) {Species[indSpecies].set_k(new_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_k1	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Definition at line 183 of file MD.h.

References Species.

184 {if (indSpecies<Species.size()) {Species[indSpecies].set_k1(new_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_k2max	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Definition at line 185 of file MD.h.

References Species.

186 {if (indSpecies<Species.size()) {Species[indSpecies].set_k2max(new_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_k_and_restitution_coefficient	(	Mdouble	k_,
		Mdouble	eps,
		Mdouble	mass,
		unsigned int	indSpecies = `0`
	)

inline

Sets k, disp such that it matches a given tc and eps for a collision of two copies of P.

Definition at line 272 of file MD.h.

References Species.

                                                                                                               {
         if (indSpecies<Species.size()) {Species[indSpecies].set_k_and_restitution_coefficient(k_, eps, mass);} else {std::cerr << "Error in set_k_and_restitution_coefficient: species does not exist"; exit(-1);}
     }

void MD::set_kc	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Definition at line 187 of file MD.h.

References Species.

188 {if (indSpecies<Species.size()) {Species[indSpecies].set_kc(new_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_krolling	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the spring constant to be changed.

Definition at line 212 of file MD.h.

References Species.

Referenced by readNextArgument().

212 {if (indSpecies<Species.size()) {Species[indSpecies].set_krolling(new_);} else {std::cerr << "Error in set_krolling: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_kt	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the spring constant to be changed.

Definition at line 208 of file MD.h.

References Species.

Referenced by readNextArgument().

208 {if (indSpecies<Species.size()) {Species[indSpecies].set_kt(new_);} else {std::cerr << "Error in set_kt: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_ktorsion	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the spring constant to be changed.

Definition at line 216 of file MD.h.

References Species.

216 {if (indSpecies<Species.size()) {Species[indSpecies].set_ktorsion(new_);} else {std::cerr << "Error in set_ktorsion: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_MixedSpecies	(	int	i,
		int	j,
		CSpecies &	S
	)

inline

Allows the mixed species to be set.

Definition at line 136 of file MD.h.

References Species.

                                                      {
         if (i>j) Species[i].MixedSpecies[j] = S;
         else Species[j].MixedSpecies[i] = S;
     };

void MD::set_mu	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the Coulomb friction coefficient to be changed.

Definition at line 245 of file MD.h.

References Species.

Referenced by ChuteBottom::make_rough_bottom(), and readNextArgument().

245 {if (indSpecies<Species.size()) {Species[indSpecies].set_mu(new_);} else {std::cerr << "Error in set_mu: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_murolling	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the Coulomb friction coefficient to be changed.

Definition at line 249 of file MD.h.

References Species.

Referenced by readNextArgument().

249 {if (indSpecies<Species.size()) {Species[indSpecies].set_murolling(new_);} else {std::cerr << "Error in set_murolling: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_mutorsion	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the Coulomb friction coefficient to be changed.

Definition at line 253 of file MD.h.

References Species.

253 {if (indSpecies<Species.size()) {Species[indSpecies].set_mutorsion(new_);} else {std::cerr << "Error in set_mutorsion: species does not exist"; exit(-1);}}

MD::set_number_of_saves_all

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_name ( const char * name )

inline

Sets the name of the problem, used for the same data files.

Definition at line 342 of file MD.h.

References STD_save::problem_name.

Referenced by ChuteBottom::constructor(), and readNextArgument().

342 {problem_name.str(""); problem_name << name;}

STD_save::problem_name

std::stringstream problem_name

Stores the problem_name.

Definition: STD_save.h:242

void MD::set_number_of_saves ( Mdouble N )

inline

Definition at line 171 of file MD.h.

References set_number_of_saves_all().

171 {set_number_of_saves_all(N);}

void set_number_of_saves_all(Mdouble N)

Definition: MD.h:172

void MD::set_number_of_saves_all ( Mdouble N )

inline

Definition at line 172 of file MD.h.

References dt, set_save_count_all(), and tmax.

Referenced by readNextArgument(), and set_number_of_saves().

172 {if (dt) {set_save_count_all((N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000);} else {std::cerr << "Error in set_number_of_saves_all ; dt must be set"; exit(-1);} }

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

MD::set_save_count_all

void set_save_count_all(int new_)

Definition: MD.h:156

Mdouble tmax

Definition: MD.h:672

void MD::set_number_of_saves_data ( Mdouble N )

inline

Definition at line 173 of file MD.h.

References dt, save_count_data, and tmax.

173 {if (dt) {save_count_data = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {std::cerr << "Error in set_number_of_saves_data ; dt must be set"; exit(-1);} }

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

Mdouble tmax

Definition: MD.h:672

int save_count_data

These are informations for saving.

Definition: MD.h:675

void MD::set_number_of_saves_ene ( Mdouble N )

inline

Definition at line 174 of file MD.h.

References dt, save_count_ene, and tmax.

174 {if (dt) {save_count_ene = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {std::cerr << "Error in set_number_of_saves_ene ; dt must be set"; exit(-1);} }

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

int save_count_ene

Definition: MD.h:676

Mdouble tmax

Definition: MD.h:672

void MD::set_number_of_saves_fstat ( Mdouble N )

inline

Definition at line 176 of file MD.h.

References dt, save_count_fstat, and tmax.

176 {if (dt) {save_count_fstat = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {std::cerr << "Error in set_number_of_saves_fstat; dt must be set"; exit(-1);} }

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

int save_count_fstat

Definition: MD.h:678

Mdouble tmax

Definition: MD.h:672

void MD::set_number_of_saves_stat ( Mdouble N )

inline

Definition at line 175 of file MD.h.

References dt, save_count_stat, and tmax.

175 {if (dt) {save_count_stat = (N>1)?((int)ceil(tmax/dt/(N-1.0))):1000000; } else {std::cerr << "Error in set_number_of_saves_stat ; dt must be set"; exit(-1);} }

Mdouble dt

These are numerical constants like the time step size.

Definition: MD.h:671

int save_count_stat

Definition: MD.h:677

Mdouble tmax

Definition: MD.h:672

void MD::set_plastic_k1_k2max_kc_depth	(	Mdouble	k1_,
		Mdouble	k2max_,
		Mdouble	kc_,
		Mdouble	depth_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the plastic constants to be changed.

Todo:: {these functions should also update the mixed species}

Definition at line 181 of file MD.h.

References Species.

182 {if (indSpecies<Species.size()) {Species[indSpecies].set_plastic_k1_k2max_kc_depth(k1_, k2max_, kc_, depth_);} else {std::cerr << "Error in set_k: species does not exist"; exit(-1);} }

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_restart_version ( int new_ )

inline

Sets restart_version.

Definition at line 374 of file MD.h.

References restart_version.

374 {restart_version=new_;}

MD::restart_version

int restart_version

Definition: MD.h:700

void MD::set_restarted ( bool new_ )

inline

Definition at line 383 of file MD.h.

References restarted, and set_append().

Referenced by constructor(), and load_restart_data().

                                  {
         restarted=new_;
         set_append(new_);
     }

void MD::set_rho	(	Mdouble	new_,
		unsigned int	indSpecies = `0`
	)

inline

Allows the density to be changed.

Definition at line 220 of file MD.h.

References Species.

Referenced by constructor(), load_par_ini_file(), and readNextArgument().

220 {if (indSpecies<Species.size()) {Species[indSpecies].set_rho(new_);} else {std::cerr << "Error in set_rho: species does not exist"; exit(-1);}}

std::vector< CSpecies > Species

These are the particle parameters like dissipation etc.

Definition: MD.h:655

void MD::set_rotation ( bool new_ )

inline

Definition at line 257 of file MD.h.

References rotation.

257 {rotation=new_;}

MD::rotation

bool rotation

Definition: MD.h:721

void MD::set_save_count_all ( int new_ )

inline

Definition at line 156 of file MD.h.

References save_count_data, save_count_ene, save_count_fstat, and save_count_stat.

Referenced by constructor(), set_number_of_saves_all(), and set_savecount().

156 {if (new_>0) {save_count_data =new_;save_count_ene =new_;save_count_stat =new_;save_count_fstat=new_;} else {std::cerr << "Error in set_save_count_all (set_save_count_all("<<new_<<"))"<< std::endl; exit(-1);}}

int save_count_ene

Definition: MD.h:676

int save_count_stat

Definition: MD.h:677

int save_count_fstat

Definition: MD.h:678

int save_count_data

These are informations for saving.

Definition: MD.h:675

void MD::set_save_count_data ( int new_ )

inline

Definition at line 157 of file MD.h.

References save_count_data.

Referenced by load_par_ini_file(), and readNextArgument().

157 {if (new_>0) {save_count_data =new_;} else {std::cerr << "Error in set_save_count_data, (set_save_count_data ("<<new_<<"))"<< std::endl; exit(-1);}}

int save_count_data

These are informations for saving.

Definition: MD.h:675

void MD::set_save_count_ene ( int new_ )

inline

Definition at line 158 of file MD.h.

References save_count_ene.

Referenced by readNextArgument().

158 {if (new_>0) {save_count_ene =new_;} else {std::cerr << "Error in set_save_count_ene, (set_save_count_ene ("<<new_<<"))"<< std::endl; exit(-1);}}

int save_count_ene

Definition: MD.h:676

void MD::set_save_count_fstat ( int new_ )

inline

Definition at line 160 of file MD.h.

References save_count_fstat.

Referenced by load_par_ini_file(), and readNextArgument().

160 {if (new_>0) {save_count_fstat=new_;} else {std::cerr << "Error in set_save_count_fstat, (set_save_count_fstat("<<new_<<"))"<< std::endl; exit(-1);}}

int save_count_fstat

Definition: MD.h:678

void MD::set_save_count_stat ( int new_ )

inline

Definition at line 159 of file MD.h.

References save_count_stat.

Referenced by readNextArgument().

159 {if (new_>0) {save_count_stat =new_;} else {std::cerr << "Error in set_save_count_stat, (set_save_count_stat ("<<new_<<"))"<< std::endl; exit(-1);}}

int save_count_stat

Definition: MD.h:677

void MD::set_savecount ( int new_ )

inline

Allows the number of time steps between saves to be changed, see also set_number_of_saves.

Definition at line 155 of file MD.h.

References set_save_count_all().

Referenced by load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().

155 {if (new_>0) {set_save_count_all(new_);} else {std::cerr << "Error in set_savecount (set_savecount("<<new_<<"))"<< std::endl; exit(-1);}}

MD::set_save_count_all

void set_save_count_all(int new_)

Definition: MD.h:156

void MD::set_t ( Mdouble new_t )

inline

Access function for the time.

Definition at line 110 of file MD.h.

References t.

Referenced by load_par_ini_file().

110 {t=new_t;}

MD::t

Mdouble t

This stores the current time.

Definition: MD.h:687

void MD::set_tmax ( Mdouble new_tmax )

inline

Allows the upper time limit to be changed.

Definition at line 142 of file MD.h.

References tmax.

Referenced by load_par_ini_file(), ChuteBottom::make_rough_bottom(), and readNextArgument().

142 {if (new_tmax>=0){tmax=new_tmax;} else { std::cerr << "Error in set_tmax, new_tmax="<<new_tmax << std::endl; exit(-1); }}

MD::xballs_additional_arguments

Mdouble tmax

Definition: MD.h:672

void MD::set_xballs_additional_arguments ( std::string new_ )

inline

Set the additional arguments for xballs.

Definition at line 354 of file MD.h.

References xballs_additional_arguments.

354 {xballs_additional_arguments=new_.c_str();}

std::string xballs_additional_arguments

Definition: MD.h:698

void MD::set_xballs_cmode ( int new_cmode )

inline

Definition at line 346 of file MD.h.

References xballs_cmode.

346 {xballs_cmode=new_cmode;}

MD::xballs_cmode

int xballs_cmode

Definition: MD.h:695

void MD::set_xballs_colour_mode ( int new_cmode )

inline

Set the xball output mode.

Definition at line 345 of file MD.h.

References xballs_cmode.

345 {xballs_cmode=new_cmode;}

MD::xballs_cmode

int xballs_cmode

Definition: MD.h:695

void MD::set_xballs_scale ( Mdouble new_scale )

inline

Set the scale of the xballs problem. The default is fit to screen.

Definition at line 358 of file MD.h.

References xballs_scale.

358 {xballs_scale=new_scale;}

MD::xballs_scale

Mdouble xballs_scale

Definition: MD.h:697

void MD::set_xballs_vector_scale ( double new_vscale )

inline

Set the scale of vectors in xballs.

Definition at line 350 of file MD.h.

References xballs_vscale.

350 {xballs_vscale=new_vscale;}

MD::xballs_vscale

Mdouble xballs_vscale

Definition: MD.h:696

void MD::set_xmax ( Mdouble new_xmax )

inline

Sets xmax and walls, assuming the standard definition of walls as in the default constructor.

Definition at line 328 of file MD.h.

References xmax, and xmin.

Referenced by readNextArgument(), Chute::set_ChuteLength(), ChuteWithHopperAndInset::set_ChuteLength(), ChuteWithHopper::set_ChuteLength(), ChuteWithHopperAndInset::set_shift(), and ChuteWithHopper::set_shift().

328 {if (new_xmax>=xmin){xmax=new_xmax;} else { std::cerr << "Error in set_xmax(" << new_xmax << "): xmin=" << xmin << std::endl; exit(-1); }}

MD::xmax

Mdouble xmax

Definition: MD.h:668

MD::xmin

Mdouble xmin

These store the size of the domain, assume walls at the ends.

Definition: MD.h:668

void MD::set_xmin ( Mdouble new_xmin )

inline

Sets xmin and walls, assuming the standard definition of walls as in the default constructor.

Definition at line 318 of file MD.h.

References xmax, and xmin.

Referenced by readNextArgument(), ChuteWithHopperAndInset::set_ChuteLength(), and ChuteWithHopper::set_ChuteLength().

318 {if (new_xmin<=xmax){xmin=new_xmin;} else { std::cerr << "Warning in set_xmin(" << new_xmin << "): xmax=" << xmax << std::endl; }}

MD::xmax

Mdouble xmax

Definition: MD.h:668

MD::xmin

Mdouble xmin

These store the size of the domain, assume walls at the ends.

Definition: MD.h:668

void MD::set_ymax ( Mdouble new_ymax )

inline

Sets ymax and walls, assuming the standard definition of walls as in the default constructor.

Definition at line 331 of file MD.h.

References ymax, and ymin.

Referenced by readNextArgument(), and Chute::set_ChuteWidth().

331 {if (new_ymax>ymin){ymax=new_ymax;} else { std::cerr << "Warning in set_ymax(" << new_ymax << "): ymin=" << ymin << std::endl; }}

Mdouble ymin

Definition: MD.h:668

MD::ymax

Mdouble ymax

Definition: MD.h:668

void MD::set_ymin ( Mdouble new_ymin )

inline

Definition at line 320 of file MD.h.

References ymax, and ymin.

Referenced by readNextArgument().

320 {if (new_ymin<=ymax){ymin=new_ymin;} else { std::cerr << "Error in set_ymin(" << new_ymin << "): ymax=" << ymax << std::endl; exit(-1); }}

Mdouble ymin

Definition: MD.h:668

MD::ymax

Mdouble ymax

Definition: MD.h:668

void MD::set_zmax ( Mdouble new_zmax )

inline

Sets ymax and walls, assuming the standard definition of walls as in the default constructor.

Definition at line 334 of file MD.h.

References ymin, zmax, and zmin.

Referenced by ChuteWithHopperAndInset::add_hopper(), ChuteWithHopper::add_hopper(), Chute::readNextArgument(), readNextArgument(), and Chute::set_InflowHeight().

334 {if (new_zmax>ymin){zmax=new_zmax;} else { std::cerr << "Error in set_zmax(" << new_zmax << "): zmin=" << zmin << std::endl; exit(-1); }}

MD::zmax

Mdouble zmax

Definition: MD.h:668