Author: eudoxos
Date: 2009-04-02 16:31:23 +0200 (Thu, 02 Apr 2009)
New Revision: 1744
Modified:
trunk/extra/Brefcom.cpp
trunk/extra/Brefcom.hpp
trunk/extra/BrefcomTestGen.cpp
trunk/extra/usct/UniaxialStrainControlledTest.cpp
trunk/extra/usct/UniaxialStrainControlledTest.hpp
trunk/gui/py/eudoxos.py
trunk/gui/py/plot.py
trunk/gui/py/yadeControl.cpp
Log:
1. Remove cruft from brefcom, a few fixes there
2. Add gnuplot grid to gnuplot by default
3. Fix speed in USCT
4. Fix InteractionDispatchers traversal in pyOmega
Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/extra/Brefcom.cpp 2009-04-02 14:31:23 UTC (rev 1744)
@@ -82,7 +82,7 @@
//Real E=(E12 /* was here for Kn: *S12/d0
*/)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
//Real E=E12; // apply alpha, beta, gamma: garbage values of E
!?
- if(!neverDamage) { assert(!isnan(sigmaT));
assert(!isnan(xiShear));}
+ if(!neverDamage) { assert(!isnan(sigmaT)); }
shared_ptr<BrefcomContact> contPhys(new BrefcomContact());
@@ -93,9 +93,7 @@
contPhys->crossSection=S12;
contPhys->epsCrackOnset=epsCrackOnset;
contPhys->epsFracture=relDuctility*epsCrackOnset;
- contPhys->xiShear=xiShear;
contPhys->omegaThreshold=omegaThreshold;
- contPhys->transStrainCoeff=transStrainCoeff;
// inherited from NormalShearInteracion, used in the timestepper
contPhys->kn=contPhys->E*contPhys->crossSection;
contPhys->ks=contPhys->G*contPhys->crossSection;
@@ -104,7 +102,7 @@
if(cohesiveThresholdIter<0 ||
Omega::instance().getCurrentIteration()<cohesiveThresholdIter)
contPhys->isCohesive=true;
else contPhys->isCohesive=false;
contPhys->dmgTau=dmgTau;
- contPhys->dmgRateExp=plRateExp;
+ contPhys->dmgRateExp=dmgRateExp;
contPhys->plTau=plTau;
contPhys->plRateExp=plRateExp;
@@ -201,6 +199,9 @@
assert(contGeom->hasShear);
//timingDeltas->checkpoint("setup");
epsN=contGeom->epsN(); epsT=contGeom->epsT();
+ if(isnan(epsN)){
+ LOG_ERROR("d0="<<contGeom->d0<<",
d1,d2="<<contGeom->d1<<","<<contGeom->d2<<";
pos1,pos2="<<contGeom->pos1<<","<<contGeom->pos2);
+ }
NNAN(epsN); NNANV(epsT);
// already in SpheresContactGeometry:
// contGeom->relocateContactPoints(); // allow very large mutual
rotations
@@ -346,77 +347,3 @@
}
-/*****************************************************************************************************
- ********************* static versions of damage evolution law and fracture
energy integration *******
-
*****************************************************************************************************/
-
-#if 0
-/*! Damage evolution law (static version: all parameters are explicitly
passed).
- *
- * This function is zero for any eps<=epsCrackOnset (no damage before the
material starts cracking),
- * therefore it will be linear-elastic. Between epsCrackOnset and espFracture
(complete damage), it
- * is a non-decreasing function (an exp curve in this case, controlled by a
single parameter expBending).
- * For any exp>=expFracture, return 1 (complete damage).
- *
- * @param kappaD maximum positive strain so far
- * @param expBending determines whether the function is bent up or down (and
how much)
- * @param epsCrackOnsert strain at which material begins to behave non-linearly
- * @param epsFracture strain at which material is fully damaged
- */
-Real BrefcomLaw::damageEvolutionLaw_static(Real kappaD, Real expBending, Real
epsCrackOnset, Real epsFracture){
- //double g(double x, double c, double eps0, double eps1){
- if(kappaD<epsCrackOnset)return 0;
- if(kappaD>epsFracture)return 1;
- return
(1/(1-exp(-expBending)))*(1-exp(-expBending*(kappaD-epsCrackOnset)/(epsFracture-epsCrackOnset)));
-}
-/*! Compute fracture energy by numerically integrating
damageEvolutionLaw_static, using unitary stiffness.
- *
- * The value returned must be multiplied by E to obtain real fracture energy:
Fracture energy
- * Gf=integral(E*(1-damage(eps))*eps)=E*integral(...) and we compute just
the integral(...) part.
- * The integration uses trapezoidal rule. All parameters except steps have the
same
- * meaning as for damageEvolutionLaw_static.
- *
- * @param steps Number of subdivision intervals for integration.
- */
-
-Real BrefcomLaw::unitFractureEnergy(double expBending, double epsCrackOnset,
double epsFracture, int steps /*=50*/){
- assert(steps>=10);
- const double lo=0, hi=epsFracture;
- double sum=0,stepSize=(hi-lo)/steps,x1,x2;
- for(int i=0; i<steps; i++){
- x1=lo+i*stepSize; x2=x1+stepSize;
-
sum+=((1-damageEvolutionLaw_static(x1,expBending,epsCrackOnset,epsFracture))*x1+(1-damageEvolutionLaw_static(x2,expBending,epsCrackOnset,epsFracture))*x2)*.5*stepSize;
/* trapezoid rule */
- }
- return sum;
-}
-
-/*! Calibrate epsFracture so that the desired fracture energy is obtained,
while other parameters are keps constant.
- *
- * The iteration relies on the fact that fractureEnergy is increasing for
increasing epsFracture;
- * first we find some strain value eps at which
fractureEnergy(eps,...)<Gf<fractureEnergy(2*eps,...),
- * then the interval is bisected until fractureEnergy matches Gf with desired
precision.
- *
- * @param Gf the desired fracture energy.
- * @param relEpsilon relative (with regards to Gf) result precision.
- * @param maxIter number of iterations after which we throw exception, since
for a reason we don't converge or converge too slowly.
- *
- * Other params have the same meaning as for damageEvolutionLaw_static.
- */
-Real BrefcomLaw::calibrateEpsFracture(double Gf, double E, double expBending,
double epsCrackOnset, double epsFractureInit /*=1e-3*/, double relEpsilon
/*=1e-3*/, int maxIter /*=1000*/){
- double
E1,E2,Emid,epsLo=epsFractureInit,epsHi=epsFractureInit*2,epsMid,iter=0,Gf_div_E=Gf/E;
- bool goUp=unitFractureEnergy(expBending,epsCrackOnset,epsHi)<Gf_div_E;
// do we double up or down when finding margins?
- do { epsLo*=goUp?2:.5; epsHi=2*epsLo;
- E1=unitFractureEnergy(expBending,epsCrackOnset,epsLo);
E2=unitFractureEnergy(expBending,epsCrackOnset,epsHi);
- if((iter++)>maxIter) throw runtime_error("Convergence problem
when finding margin values for bisection.");
- } while(!(E1<Gf_div_E && E2>=Gf_div_E));
- // now E(epsLo)<Gf_div_E<=E(epsHi); go ahead using interval bisection
- do {
- epsMid=.5*(epsLo+epsHi);
Emid=unitFractureEnergy(expBending,epsCrackOnset,epsMid);
- if(Emid<Gf_div_E)epsLo=.5*(epsLo+epsHi);
- else if(Emid>Gf_div_E)epsHi=.5*(epsLo+epsHi);
- if((iter++)>maxIter) throw runtime_error("Convergence problem
during bisection (relEpsilon too low?).");
- } while (abs(Emid-Gf_div_E)>relEpsilon*Gf_div_E);
- return epsMid;
-}
-#endif
-
Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/extra/Brefcom.hpp 2009-04-02 14:31:23 UTC (rev 1744)
@@ -222,10 +222,8 @@
* This might be done later, for now hardcode that here. */
/* uniaxial tension resistance, bending parameter of the damage
evolution law, whear weighting constant for epsT in the strain seminorm (kappa)
calculation. Default to NaN so that user gets loudly notified it was not set.
- expBending is positive if the damage evolution function is
concave after fracture onset;
- reasonable value seems like 4.
*/
- Real sigmaT, expBending, xiShear, epsCrackOnset, relDuctility,
G_over_E, tau, expDmgRate, omegaThreshold, transStrainCoeff, dmgTau,
dmgRateExp, plTau, plRateExp;
+ Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau,
expDmgRate, omegaThreshold, dmgTau, dmgRateExp, plTau, plRateExp;
//! Should new contacts be cohesive? They will before this
iter#, they will not be afterwards. If 0, they will never be. If negative, they
will always be created as cohesive.
long cohesiveThresholdIter;
//! Create contacts that don't receive any damage
(BrefcomContact::neverDamage=true); defaults to false
@@ -233,8 +231,7 @@
BrefcomMakeContact(){
// init to signaling_NaN to force crash if not
initialized (better than unknowingly using garbage values)
-
sigmaT=epsCrackOnset=relDuctility=G_over_E=transStrainCoeff=std::numeric_limits<Real>::signaling_NaN();
- xiShear=0;
+
sigmaT=epsCrackOnset=relDuctility=G_over_E=std::numeric_limits<Real>::signaling_NaN();
neverDamage=false;
cohesiveThresholdIter=-1;
dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
@@ -245,8 +242,6 @@
REGISTER_ATTRIBUTES(InteractionPhysicsEngineUnit,
(cohesiveThresholdIter)
(G_over_E)
- (expBending)
- (xiShear)
(sigmaT)
(neverDamage)
(epsCrackOnset)
@@ -256,7 +251,6 @@
(plTau)
(plRateExp)
(omegaThreshold)
- (transStrainCoeff)
);
FUNCTOR2D(BrefcomPhysParams,BrefcomPhysParams);
Modified: trunk/extra/BrefcomTestGen.cpp
===================================================================
--- trunk/extra/BrefcomTestGen.cpp 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/extra/BrefcomTestGen.cpp 2009-04-02 14:31:23 UTC (rev 1744)
@@ -57,7 +57,7 @@
shared_ptr<InteractionPhysicsMetaEngine> iphysDispatcher(new
InteractionPhysicsMetaEngine);
shared_ptr<BrefcomMakeContact> bmc(new BrefcomMakeContact);
- bmc->cohesiveThresholdIter=-1; bmc->G_over_E=1;
bmc->expBending=1; bmc->xiShear=.8; bmc->sigmaT=3e9; bmc->neverDamage=true;
bmc->epsCrackOnset=1e-4; bmc->relDuctility=5; bmc->transStrainCoeff=.5;
+ bmc->cohesiveThresholdIter=-1; bmc->G_over_E=1;
bmc->sigmaT=3e9; bmc->neverDamage=true; bmc->epsCrackOnset=1e-4;
bmc->relDuctility=5;
//bmc->calibratedEpsFracture=.5; /* arbitrary, but large enough
*/
iphysDispatcher->add(bmc);
rootBody->engines.push_back(iphysDispatcher);
Modified: trunk/extra/usct/UniaxialStrainControlledTest.cpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-04-02 09:28:06 UTC
(rev 1743)
+++ trunk/extra/usct/UniaxialStrainControlledTest.cpp 2009-04-02 14:31:23 UTC
(rev 1744)
@@ -49,9 +49,10 @@
assert(originalLength>0 && !isnan(originalLength));
assert(!isnan(strainRate) || !isnan(absSpeed));
- if(isnan(strainRate)){ strainRate=absSpeed/originalLength; }
+ if(strainRate==0){ strainRate=absSpeed/originalLength;
LOG_INFO("COmputed new strainRate "<<strainRate); }
initAccelTime_s=initAccelTime>=0 ? initAccelTime :
Omega::instance().getTimeStep()*(-initAccelTime);
+ LOG_INFO("Strain speed will be "<<absSpeed<<", strain rate
"<<strainRate<<", will be reached after "<<initAccelTime_s<<"s
("<<initAccelTime_s/Omega::instance().getTimeStep()<<" steps).");
/* if we have default (<0) crossSectionArea, try to get it from root's
AABB;
* this will not work if there are foreign bodies in the simulation,
@@ -89,7 +90,8 @@
// linearly increase strain to the desired value
if(abs(currentStrainRate)<abs(strainRate)){
Real t=Omega::instance().getSimulationTime();
- currentStrainRate=(t/initAccelTime_s)*strainRate;
+ if(initAccelTime_s!=0)
currentStrainRate=(t/initAccelTime_s)*strainRate;
+ else currentStrainRate=strainRate;
} else currentStrainRate=strainRate;
// how much do we move (in total, symmetry handled below)
Real
dAX=currentStrainRate*originalLength*Omega::instance().getTimeStep();
@@ -250,7 +252,7 @@
shared_ptr<InteractionPhysicsMetaEngine> iphysDispatcher(new
InteractionPhysicsMetaEngine);
shared_ptr<BrefcomMakeContact> bmc(new BrefcomMakeContact);
bmc->cohesiveThresholdIter=cohesiveThresholdIter;
- bmc->cohesiveThresholdIter=-1; bmc->G_over_E=1;
bmc->expBending=1; bmc->xiShear=.8; bmc->sigmaT=3e9; bmc->neverDamage=true;
bmc->epsCrackOnset=1e-4; bmc->relDuctility=5; bmc->transStrainCoeff=.5;
+ bmc->cohesiveThresholdIter=-1; bmc->G_over_E=1;bmc->sigmaT=3e9;
bmc->neverDamage=true; bmc->epsCrackOnset=1e-4; bmc->relDuctility=5;
iphysDispatcher->add(bmc);
rootBody->engines.push_back(iphysDispatcher);
Modified: trunk/extra/usct/UniaxialStrainControlledTest.hpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.hpp 2009-04-02 09:28:06 UTC
(rev 1743)
+++ trunk/extra/usct/UniaxialStrainControlledTest.hpp 2009-04-02 14:31:23 UTC
(rev 1744)
@@ -110,7 +110,7 @@
Real strain, avgStress;
virtual void applyCondition(MetaBody* rootBody);
- UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0;
originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1;
needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */
recordFile=""; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false;
blockDisplacements=false;
absSpeed=strainRate=stopStrain=numeric_limits<Real>::quiet_NaN(); active=true;
idleIterations=0; };
+ UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0;
originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1;
needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */
recordFile=""; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false;
blockDisplacements=false; strainRate=0;
absSpeed=stopStrain=numeric_limits<Real>::quiet_NaN(); active=true;
idleIterations=0; initAccelTime=-200;};
virtual ~UniaxialStrainer(){};
REGISTER_ATTRIBUTES(DeusExMachina,
(strainRate)
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/gui/py/eudoxos.py 2009-04-02 14:31:23 UTC (rev 1744)
@@ -174,7 +174,7 @@
f.write("SimpleCS 1 thick 1.0 width 1.0\n")
# material
ph=inters[0].phys
- f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g
d
1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle']))
+ f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g
damchartime %g damrateexp %g d
1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle'],ph['dmgTau'],ph['dmgRateExp']))
# boundary conditions
f.write('BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0\n')
f.write('BoundaryCondition 2 loadTimeFunction 1 prescribedvalue
1.e-4\n')
Modified: trunk/gui/py/plot.py
===================================================================
--- trunk/gui/py/plot.py 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/gui/py/plot.py 2009-04-02 14:31:23 UTC (rev 1744)
@@ -147,7 +147,8 @@
if term in ['wxt','x11']: fPlot.write("set term %s %d
persist\n"%(term,i))
else: fPlot.write("set term %s; set output
'%s.%d.%s'\n"%(term,baseNameNoPath,i,extension))
fPlot.write("set xlabel '%s'\n"%p)
- fPlot.write("set datafile missing 'nan'\n"%p)
+ fPlot.write("set grid\n")
+ fPlot.write("set datafile missing 'nan'\n")
if title: fPlot.write("set title '%s'\n"%title)
fPlot.write("plot "+",".join([" %s using %d:%d title '%s(%s)'
with lines"%(dataFile,vars.index(p)+1,vars.index(pp[0])+1,pp[0],p) for pp in
plots_p])+"\n")
i+=1
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-04-02 09:28:06 UTC (rev 1743)
+++ trunk/gui/py/yadeControl.cpp 2009-04-02 14:31:23 UTC (rev 1744)
@@ -449,7 +449,7 @@
}
}
}
-
if(isChildClassOf(e->getClassName(),"InteractionDispatchers")){
+ if(e->getClassName()=="InteractionDispatchers"){
shared_ptr<InteractionDispatchers>
ee=dynamic_pointer_cast<InteractionDispatchers>(e);
list<shared_ptr<EngineUnit> > eus;
FOREACH(const shared_ptr<EngineUnit>&
eu,ee->geomDispatcher->functorArguments) eus.push_back(eu);
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/yade-dev
