yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2430
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: trunk
timestamp: Thu 2010-09-09 22:47:40 +0200
message:
  1. Changed some Mathr:: functions to standard analogues, discussed here http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg05026.html
modified:
  lib/base/Math.hpp
  pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
  pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
  pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp
  pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp
  pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp
  pkg/dem/Engine/PartialEngine/KinemSimpleShearBox.cpp
  pkg/dem/Engine/PartialEngine/TriaxialStressController.cpp
  pkg/dem/PreProcessor/SimpleShear.cpp
  pkg/dem/meta/Shop.cpp
  pkg/dem/meta/ViscoelasticPM.cpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp	2010-07-24 18:10:24 +0000
+++ lib/base/Math.hpp	2010-09-09 20:47:40 +0000
@@ -95,13 +95,6 @@
 	static const Scalar EPSILON;
 	static const Scalar ZERO_TOLERANCE;
 	static Scalar Sign(Scalar f){ if(f<0) return -1; if(f>0) return 1; return 0; }
-	static Scalar FAbs(Scalar f){ return abs(f); }
-	static Scalar Sqrt(Scalar f){ return sqrt(f); }
-	static Scalar Log(Scalar f){ return log(f); }
-	static Scalar Exp(Scalar f){ return exp(f); }
-	static Scalar ATan(Scalar f){ return atan(f); }
-	static Scalar Tan(Scalar f){ return tan(f); }
-	static Scalar Pow(Scalar base,Scalar exponent){ return pow(base,exponent); }
 
 	static Scalar UnitRandom(){ return ((double)rand()/((double)(RAND_MAX))); }
 	static Scalar SymmetricRandom(){ return 2.*(((double)rand())/((double)(RAND_MAX)))-1.; }

=== modified file 'pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp'
--- pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp	2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp	2010-09-09 20:47:40 +0000
@@ -62,7 +62,7 @@
 	Real maxFs = phys->normalForce.squaredNorm() * std::pow(phys->tangensOfFrictionAngle,2);
 	if( shearForce.squaredNorm() > maxFs )
 	{
-		maxFs = Mathr::Sqrt(maxFs) / shearForce.norm();
+		maxFs = sqrt(maxFs) / shearForce.norm();
 		shearForce *= maxFs;
 	}
 

=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-09-09 14:02:49 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp	2010-09-09 20:47:40 +0000
@@ -71,12 +71,12 @@
 	if (!traceEnergy){//Update force but don't compute energy terms (see below))
 		// PFC3d SlipModel, is using friction angle. CoulombCriterion
 		if( shearForce.squaredNorm() > maxFs ){
-			Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
+			Real ratio = sqrt(maxFs) / shearForce.norm();
 			shearForce *= ratio;}
 	} else {
 		//almost the same with additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
 		if( shearForce.squaredNorm() > maxFs ){
-			Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
+			Real ratio = sqrt(maxFs) / shearForce.norm();
 			Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
 			//define the plastic work input and increment the total plastic energy dissipated
 			shearForce *= ratio;

=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-09-09 14:02:49 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp	2010-09-09 20:47:40 +0000
@@ -195,8 +195,8 @@
 	if (useDamping){
 		Real mbar = (!b1->isDynamic() && b2->isDynamic()) ? de2->mass : ((!b2->isDynamic() && b1->isDynamic()) ? de1->mass : (de1->mass*de2->mass / (de1->mass + de2->mass))); // get equivalent mass if both bodies are dynamic, if not set it equal to the one of the dynamic body
 		//Real mbar = de1->mass*de2->mass / (de1->mass + de2->mass); // equivalent mass
-		Real Cn_crit = 2.*Mathr::Sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)
-		Real Cs_crit = 2.*Mathr::Sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction)
+		Real Cn_crit = 2.*sqrt(mbar*phys->kn); // Critical damping coefficient (normal direction)
+		Real Cs_crit = 2.*sqrt(mbar*phys->ks); // Critical damping coefficient (shear direction)
 		// Note: to compare with the analytical solution you provide cn and cs directly (since here we used a different method to define c_crit)
 		cn = Cn_crit*betan; // Damping normal coefficient
 		cs = Cs_crit*betas; // Damping tangential coefficient
@@ -339,7 +339,7 @@
 	Real uN = scg->penetrationDepth; // get overlapping  
 	if (uN<0) {ncb->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
 	/*** Hertz-Mindlin's formulation (PFC) ***/
-	phys->kn = phys->kno*Mathr::Sqrt(uN); // normal stiffness
+	phys->kn = phys->kno*sqrt(uN); // normal stiffness
 	Real Fn = phys->kn*uN; // normal Force (scalar)
 	phys->normalForce = Fn*scg->normal; // normal Force (vector)
 
@@ -355,7 +355,7 @@
 	/*** MOHR-COULOMB LAW ***/
 	Real maxFs = phys->normalForce.squaredNorm();
 	if (trialFs.squaredNorm() > maxFs)
-	{Real ratio = Mathr::Sqrt(maxFs)/trialFs.norm(); trialFs *= ratio;}
+	{Real ratio = sqrt(maxFs)/trialFs.norm(); trialFs *= ratio;}
 	
 
 	/*** APPLY FORCES ***/

=== modified file 'pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp'
--- pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-09-09 09:51:23 +0000
+++ pkg/dem/Engine/GlobalEngine/PeriIsoCompressor.cpp	2010-09-09 20:47:40 +0000
@@ -97,7 +97,7 @@
 void PeriTriaxController::strainStressStiffUpdate(){
 	// update strain first
 	//"Natural" strain, correct for large deformations, only used for comparison with goals
-	for (int i=0;i<3;i++) strain[i]=Mathr::Log(scene->cell->trsf(i,i));
+	for (int i=0;i<3;i++) strain[i]=log(scene->cell->trsf(i,i));
 	//stress tensor and stiffness
 	
 	//Compute volume of the deformed cell
@@ -198,7 +198,7 @@
 			
 		} else {    // control strain, see "true strain" definition here http://en.wikipedia.org/wiki/Finite_strain_theory
 			///NOTE : everything could be generalized to 9 independant components by comparing F[i,i] vs. Matrix3r goal[i,i], but it would be simpler in that case to let the user set the prescribed loading rates velGrad[i,i] when [i,i] is not stress-controlled. This "else" would disappear.
-			strain_rate = (Mathr::Exp ( goal[axis]-strain[axis] ) -1)/scene->dt;
+			strain_rate = (exp ( goal[axis]-strain[axis] ) -1)/scene->dt;
 			LOG_TRACE ( axis<<": strain="<<strain[axis]<<", goal="<<goal[axis]<<", cellGrow="<<strain_rate*scene->dt);
 		}		
 		// steady evolution with fluctuations; see TriaxialStressController

=== modified file 'pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp'
--- pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp	2010-08-24 12:54:14 +0000
+++ pkg/dem/Engine/PartialEngine/Disp2DPropLoadEngine.cpp	2010-09-09 20:47:40 +0000
@@ -113,12 +113,12 @@
 	computeAlpha();
 	if (alpha == Mathr::PI/2.0)	// Case of the very beginning
 	{
-		dalpha = - Mathr::ATan( dgamma / (Ysup_mod -Ylat_mod) );
+		dalpha = - atan( dgamma / (Ysup_mod -Ylat_mod) );
 	}
 	else
 	{
-		Real A = (Ysup_mod - Ylat_mod) * 2.0*Mathr::Tan(alpha) / (2.0*(Ysup - Ylat) + dgamma*Mathr::Tan(alpha) );
-		dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
+		Real A = (Ysup_mod - Ylat_mod) * 2.0*tan(alpha) / (2.0*(Ysup - Ylat) + dgamma*tan(alpha) );
+		dalpha = atan( (A - tan(alpha))/(1.0 + A * tan(alpha)));
 	}
 
 	Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));

=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp'
--- pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp	2010-08-16 21:31:08 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNSEngine.cpp	2010-09-09 20:47:40 +0000
@@ -76,12 +76,12 @@
 //	Then computation of the angle of the rotation to be done :
 	if (alpha == Mathr::PI/2.0)	// Case of the very beginning
 	{
-		dalpha = - Mathr::ATan( dx / (Ysup_mod -Ylat_mod) );
+		dalpha = - atan( dx / (Ysup_mod -Ylat_mod) );
 	}
 	else
 	{
-		Real A = (Ysup_mod - Ylat_mod) * 2.0*Mathr::Tan(alpha) / (2.0*(Ysup - Ylat) + dx*Mathr::Tan(alpha) );
-		dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
+		Real A = (Ysup_mod - Ylat_mod) * 2.0*tan(alpha) / (2.0*(Ysup - Ylat) + dx*tan(alpha) );
+		dalpha = atan( (A - tan(alpha))/(1.0 + A * tan(alpha)));
 	}
 
 	Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));

=== modified file 'pkg/dem/Engine/PartialEngine/KinemSimpleShearBox.cpp'
--- pkg/dem/Engine/PartialEngine/KinemSimpleShearBox.cpp	2010-08-16 21:31:08 +0000
+++ pkg/dem/Engine/PartialEngine/KinemSimpleShearBox.cpp	2010-09-09 20:47:40 +0000
@@ -77,12 +77,12 @@
 	computeAlpha();
 	if (alpha == Mathr::PI/2.0)	// Case of the very beginning
 	{
-		dalpha = - Mathr::ATan( dX / (Ysup_mod -Ylat_mod) );
+		dalpha = - atan( dX / (Ysup_mod -Ylat_mod) );
 	}
 	else
 	{
-		Real A = (Ysup_mod - Ylat_mod) * 2.0*Mathr::Tan(alpha) / (2.0*(Ysup - Ylat) + dX*Mathr::Tan(alpha) );
-		dalpha = Mathr::ATan( (A - Mathr::Tan(alpha))/(1.0 + A * Mathr::Tan(alpha)));
+		Real A = (Ysup_mod - Ylat_mod) * 2.0*tan(alpha) / (2.0*(Ysup - Ylat) + dX*tan(alpha) );
+		dalpha = atan( (A - tan(alpha))/(1.0 + A * tan(alpha)));
 	}
 
 	Quaternionr qcorr(AngleAxisr(dalpha,Vector3r::UnitZ()));

=== modified file 'pkg/dem/Engine/PartialEngine/TriaxialStressController.cpp'
--- pkg/dem/Engine/PartialEngine/TriaxialStressController.cpp	2010-08-16 21:31:08 +0000
+++ pkg/dem/Engine/PartialEngine/TriaxialStressController.cpp	2010-09-09 20:47:40 +0000
@@ -185,9 +185,9 @@
 	if (height0 == 0) height0 = height;
 	if (width0 == 0) width0 = width;
 	if (depth0 == 0) depth0 = depth;
-	strain[0] = Mathr::Log(width0/width);
-	strain[1] = Mathr::Log(height0/height);
-	strain[2] = Mathr::Log(depth0/depth);
+	strain[0] = log(width0/width);
+	strain[1] = log(height0/height);
+	strain[2] = log(depth0/depth);
 	volumetricStrain=strain[0]+strain[1]+strain[2];
 	
 	Real invXSurface = 1.f/(height*depth);
@@ -255,7 +255,7 @@
 			Real f = (static_cast<FrictPhys*> ((contact->interactionPhysics.get()))->normalForce+static_cast<FrictPhys*>(contact->interactionPhysics.get())->shearForce).squaredNorm();
 			if (f!=0)
 			{
-			MeanForce += Mathr::Sqrt(f);
+			MeanForce += sqrt(f);
 			++nForce;
 			}
 		}

=== modified file 'pkg/dem/PreProcessor/SimpleShear.cpp'
--- pkg/dem/PreProcessor/SimpleShear.cpp	2010-08-15 17:53:02 +0000
+++ pkg/dem/PreProcessor/SimpleShear.cpp	2010-09-09 20:47:40 +0000
@@ -232,7 +232,7 @@
 	sphere_list.clear();
 	long tries = 1000; //nb max of tries for positionning the next sphere
 	Vector3r dimensions = upperCorner - lowerCorner;
-	Real mean_radius = Mathr::Pow(dimensions.x()*dimensions.y()*dimensions.z()*(1-porosity)/(4.0/3.0*Mathr::PI*number),1.0/3.0);
+	Real mean_radius = pow(dimensions.x()*dimensions.y()*dimensions.z()*(1-porosity)/(4.0/3.0*Mathr::PI*number),1.0/3.0);
 	cerr << " mean radius " << mean_radius << endl;;
 
 // 	std::cerr << "generating aggregates ... ";

=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp	2010-09-09 09:51:23 +0000
+++ pkg/dem/meta/Shop.cpp	2010-09-09 20:47:40 +0000
@@ -1141,10 +1141,10 @@
 
 void Shop::getViscoelasticFromSpheresInteraction( Real m, Real tc, Real en, Real es, shared_ptr<ViscElMat> b)
 {
-    b->kn = m/tc/tc * ( Mathr::PI*Mathr::PI + Mathr::Pow(Mathr::Log(en),2) );
-    b->cn = -2.0*m/tc * Mathr::Log(en);
-    b->ks = 2.0/7.0 * m/tc/tc * ( Mathr::PI*Mathr::PI + Mathr::Pow(Mathr::Log(es),2) );
-    b->cs = -2.0/7.0 * m/tc * Mathr::Log(es);
+    b->kn = m/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(en),2) );
+    b->cn = -2.0*m/tc * log(en);
+    b->ks = 2.0/7.0 * m/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(es),2) );
+    b->cs = -2.0/7.0 * m/tc * log(es);
 
     if (abs(b->cn) <= Mathr::ZERO_TOLERANCE ) b->cn=0;
     if (abs(b->cs) <= Mathr::ZERO_TOLERANCE ) b->cs=0;

=== modified file 'pkg/dem/meta/ViscoelasticPM.cpp'
--- pkg/dem/meta/ViscoelasticPM.cpp	2010-07-17 17:37:54 +0000
+++ pkg/dem/meta/ViscoelasticPM.cpp	2010-09-09 20:47:40 +0000
@@ -89,7 +89,7 @@
 		// Then Mohr-Coulomb is violated (so, we slip), 
 		// we have the max value of the shear force, so 
 		// we consider only friction damping.
-		const Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
+		const Real ratio = sqrt(maxFs) / shearForce.norm();
 		shearForce *= ratio;
 	} 
 	else