Author: eudoxos
Date: 2009-05-31 12:09:14 +0200 (Sun, 31 May 2009)
New Revision: 1787
Added:
trunk/lib/computational-geometry/Hull2d.hpp
Modified:
trunk/core/InteractionContainer.cpp
trunk/examples/concrete/pack/mk-pack.py
trunk/gui/py/_utils.cpp
trunk/gui/py/eudoxos.py
trunk/gui/py/utils.py
trunk/gui/py/yadeControl.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
trunk/pkg/dem/ConcretePM.cpp
trunk/pkg/dem/ConcretePM.hpp
trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp
Log:
1. Syntax fixes in ConcretePM model
2. Interactions will be deleted also if the ye4xist, but geometry functor fails
(returns false)
3. Add utils.approxSectionArea, utils.spheresFromFileUnixial (get some
parameters useful for uniaxial tests from the sphere packing, such as
cross-section, principal axis, etc)
4. Add lib/computational-geometry/Hull2d.hpp for computing convex hull in 2d
using Graham scan algorithm (I know it is in CGAL already, but I don't want to
depend on that commercial stuff)
5. Remove some unused stuff from utils.
Modified: trunk/core/InteractionContainer.cpp
===================================================================
--- trunk/core/InteractionContainer.cpp 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/core/InteractionContainer.cpp 2009-05-31 10:09:14 UTC (rev 1787)
@@ -12,7 +12,9 @@
void InteractionContainer::requestErase(body_id_t id1, body_id_t id2){
find(id1,id2)->reset(); bodyIdPair v(0,2); v.push_back(id1);
v.push_back(id2);
- boost::mutex::scoped_lock lock(pendingEraseMutex);
+ #ifdef YADE_OPENMP
+ boost::mutex::scoped_lock lock(pendingEraseMutex);
+ #endif
pendingErase.push_back(v);
}
Modified: trunk/examples/concrete/pack/mk-pack.py
===================================================================
--- trunk/examples/concrete/pack/mk-pack.py 2009-05-29 14:35:22 UTC (rev
1786)
+++ trunk/examples/concrete/pack/mk-pack.py 2009-05-31 10:09:14 UTC (rev
1787)
@@ -12,7 +12,7 @@
cutoffRound=2
tt=TriaxialTest(
- numberOfGrains=500,
+ numberOfGrains=5000,
radiusMean=3e-4,
# this is just size ratio if radiusMean is specified
# if you comment out the line above, it will be the corner (before
compaction) and radiusMean will be set accordingly
@@ -44,7 +44,7 @@
ext=utils.aabbExtrema()
rSphere=tt['radiusMean']
-outFile=tempfile.NamedTemporaryFile(delete=False)
+outFile=tempfile.NamedTemporaryFile(delete=False).name
if cylAxis<0: # box-shaped packing
aabbMin,aabbMax=tuple([ext[0][i]+rSphere*cutoffFlat for i in
0,1,2]),tuple([ext[1][i]-rSphere*cutoffFlat for i in 0,1,2])
@@ -67,7 +67,7 @@
return True
nSpheres=utils.spheresToFile(outFile,consider=isInCyl)
-outFile2='pack-%d-%s.sphere'%(nSpheres,'box' if cylAxis<0 else 'cyl')
+outFile2='pack-%d-%s.spheres'%(nSpheres,'box' if cylAxis<0 else 'cyl')
shutil.move(outFile,outFile2)
print nSpheres,'spheres written to',outFile2
quit()
Modified: trunk/gui/py/_utils.cpp
===================================================================
--- trunk/gui/py/_utils.cpp 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/gui/py/_utils.cpp 2009-05-31 10:09:14 UTC (rev 1787)
@@ -8,6 +8,7 @@
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/SimpleViscoelasticBodyParameters.hpp>
#include<yade/pkg-common/NormalShearInteractions.hpp>
+#include<yade/lib-computational-geometry/Hull2d.hpp>
#include<cmath>
#include<numpy/ndarrayobject.h>
@@ -284,7 +285,39 @@
return inside;
}
+/* Compute area of convex hull when when taking (swept) spheres crossing the
plane at coord, perpendicular to axis.
+ All spheres that touch the plane are projected as hexagons on their
circumference to the plane.
+ Convex hull from this cloud is computed.
+ The area of the hull is returned.
+
+*/
+Real approxSectionArea(Real coord, int axis){
+ std::list<Vector2r> cloud;
+ if(axis<0 || axis>2) throw invalid_argument("Axis must be ∈ {0,1,2}");
+ const int ax1=(axis+1)%3, ax2=(axis+2)%3;
+ const Real sqrt3=sqrt(3);
+ Vector2r mm(-10.,0),mx(0,0);
+ FOREACH(const shared_ptr<Body>& b,
*Omega::instance().getRootBody()->bodies){
+ Sphere* s=dynamic_cast<Sphere*>(b->geometricalModel.get());
+ if(!s) continue;
+ const Vector3r& pos(b->physicalParameters->se3.position); const
Real r(s->radius);
+ if((pos[axis]>coord && (pos[axis]-r)>coord) || (pos[axis]<coord
&& (pos[axis]+r)<coord)) continue;
+ Vector2r c(pos[ax1],pos[ax2]);
+ cloud.push_back(c+Vector2r(r,0));
cloud.push_back(c+Vector2r(-r,0));
+ cloud.push_back(c+Vector2r( r/2., sqrt3*r));
cloud.push_back(c+Vector2r( r/2.,-sqrt3*r));
+ cloud.push_back(c+Vector2r(-r/2., sqrt3*r));
cloud.push_back(c+Vector2r(-r/2.,-sqrt3*r));
+ if(mm[0]==-10.){ mm=c, mx=c;}
+ mm=Vector2r(min(c[0]-r,mm[0]),min(c[1]-r,mm[1]));
+ mx=Vector2r(max(c[0]+r,mx[0]),max(c[1]+r,mx[1]));
+ }
+ if(cloud.size()==0) return 0;
+ ConvexHull2d ch2d(cloud);
+ vector<Vector2r> hull=ch2d();
+ return simplePolygonArea2d(hull);
+}
+
+
/* Project 3d point into 2d using spiral projection along given axis;
* the returned tuple is
*
@@ -338,6 +371,7 @@
def("aabbExtrema",aabbExtrema,aabbExtrema_overloads(args("cutoff","centers"),"Return
coordinates of box enclosing all bodies\n centers: do not take sphere radii in
account, only their centroids (default=False)\n cutoff: 0-1 number by which the
box will be scaled around its center (default=0)"));
def("ptInAABB",ptInAABB,"Return True/False whether the point (3-tuple)
p is within box given by its min (3-tuple) and max (3-tuple) corners");
def("negPosExtremeIds",negPosExtremeIds,negPosExtremeIds_overloads(args("axis","distFactor"),"Return
list of ids for spheres (only) that are on extremal ends of the specimen along
given axis; distFactor multiplies their radius so that sphere that do not touch
the boundary coordinate can also be returned."));
+ def("approxSectionArea",approxSectionArea,"Compute area of convex hull
when when taking (swept) spheres crossing the plane at coord, perpendicular to
axis.");
def("coordsAndDisplacements",coordsAndDisplacements,coordsAndDisplacements_overloads(args("AABB"),"Return
tuple of 2 same-length lists for coordinates and displacements (coordinate
minus reference coordinate) along given axis (1st arg); if the
AABB=((x_min,y_min,z_min),(x_max,y_max,z_max)) box is given, only bodies within
this box will be considered."));
def("setRefSe3",setRefSe3,"Set reference positions and orientation of
all bodies equal to their current ones.");
def("interactionAnglesHistogram",interactionAnglesHistogram,interactionAnglesHistogram_overloads(args("axis","mask","bins","aabb")));
Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/gui/py/eudoxos.py 2009-05-31 10:09:14 UTC (rev 1787)
@@ -74,7 +74,7 @@
ph=o.interactions.nth(0).phys # some params are the same everywhere
material.append("%g %g"%(ph['E'],ph['G']))
- material.append("%g %g %g
%g"%(ph['epsCrackOnset'],ph['epsFracture'],1e50,ph['xiShear']))
+ material.append("%g %g %g
%g"%(ph['epsCrackOnset'],ph['epsFracture'],1e50,0.0))
material.append("%g
%g"%(ph['undamagedCohesion'],ph['tanFrictionAngle']))
# need strainer for getting bodies in positive/negative boundary
@@ -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
damchartime %g damrateexp %g plchartime %g plrateexp %g d
1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle'],ph['dmgTau'],ph['dmgRateExp'],ph['plTau'],ph['plRateExp']))
+ f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g
damchartime %g damrateexp %g plchartime %g plrateexp %g d
1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],0.0,ph['undamagedCohesion'],ph['tanFrictionAngle'],ph['dmgTau'],ph['dmgRateExp'],ph['plTau'],ph['plRateExp']))
# 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/utils.py
===================================================================
--- trunk/gui/py/utils.py 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/gui/py/utils.py 2009-05-31 10:09:14 UTC (rev 1787)
@@ -4,8 +4,6 @@
#
# 2008 © Václav Šmilauer <[email protected]>
-#from yade._utils import *
-
import math,random
from yade.wrapper import *
try: # use psyco if available
@@ -16,16 +14,6 @@
# c++ implementations for performance reasons
from yade._utils import *
-def resetEngineClocks():
- for e in O.engines: e.execTime=0
-def engineClockStats():
- tSum=sum([e.execTime for e in O.engines])
- for e in O.engines:
- print
e.name.ljust(30),(str(e.execTime/1000)+'ms').rjust(15),'%6.2f%%'%(e.execTime*100./tSum)
- print '='*53
- print 'TOTAL'.ljust(30),(str(tSum/1000)+'ms').rjust(15),'100.00%'
-
-
def saveVars(mark='',loadNow=False,**kw):
"""Save passed variables into the simulation so that it can be
recovered when the simulation is loaded again.
@@ -381,44 +369,6 @@
dictParams.update(dictDefaults); saveVars('table',**dictParams)
return len(tagsParams)
-
-def
basicDEMEngines(interPhysics='SimpleElasticRelationships',constitutiveLaw='ElasticContactLaw',gravity=None,damping=.4):
- """Set basic set of DEM engines and initializers.
-
- interPhysics and constitutiveLaw specify class of respective engines to
use instead of defaults.
-
- Gravity can be list or tuple to specify numeric value, it can also be
an object that will be inserted into
- engines, however. By default, no gravity is applied.
- """
- O.initializers=[
-
MetaEngine('BoundingVolumeMetaEngine',[EngineUnit('InteractingSphere2AABB'),EngineUnit('InteractingBox2AABB'),EngineUnit('InteractingFacet2AABB'),EngineUnit('MetaInteractingGeometry2AABB')])
- ]
- O.engines=[
- StandAloneEngine('PhysicalActionContainerReseter'),
- MetaEngine('BoundingVolumeMetaEngine',[
- EngineUnit('InteractingSphere2AABB'),
- EngineUnit('InteractingBox2AABB'),
- EngineUnit('InteractingFacet2AABB'),
- EngineUnit('MetaInteractingGeometry2AABB')
- ]),
- StandAloneEngine('PersistentSAPCollider'),
- MetaEngine('InteractionGeometryMetaEngine',[
-
EngineUnit('InteractingSphere2InteractingSphere4SpheresContactGeometry'),
-
EngineUnit('InteractingFacet2InteractingSphere4SpheresContactGeometry'),
-
EngineUnit('InteractingBox2InteractingSphere4SpheresContactGeometry')
- ]),
-
MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleElasticRelationships')]),
- StandAloneEngine('ElasticContactLaw'),
- ]
- if gravity:
- if islist(gravity) or istuple(gravity):
-
O.engines=O.engines+[DeusExMachina('GravityEngine',{'gravity':gravity}),]
- else:
- O.engines=O.engines+[gravity]
-
O.engines=O.engines+[DeusExMachina('NewtonsDampedLaw',{'damping':damping}),]
-
-
-
def
ColorizedVelocityFilter(isFilterActivated=True,autoScale=True,minValue=0,maxValue=0,posX=0,posY=0.2,width=0.05,height=0.5,title='Velocity,
m/s'):
f =
DeusExMachina('ColorizedVelocityFilter',{'isFilterActivated':isFilterActivated,'autoScale':autoScale,'minValue':minValue,'maxValue':maxValue,'posX':posX,'posY':posY,'width':width,'height':height,'title':title})
O.engines+=[f]
@@ -494,4 +444,24 @@
ret+=[sphere(xyz,radius=radius,**kw)]
return ret
+def spheresFromFileUniaxial(filename,areaSections=10,**kw):
+ """Load spheres from file, but do some additional work useful for
uniaxial test:
+
+ 1. Find the dimensions that is the longest (uniaxial loading axis)
+ 2. Find the minimum cross-section area of the speciment by examining
several (areaSections)
+ sections perpendicular to axis, computing area of the convex
hull for each one. This will
+ work also for non-prismatic specimen.
+ 3. Find the bodies that are on the negative/positive boundary, to which
the straining condition
+ should be applied.
+ Returns dictionary with keys 'negIds', 'posIds', 'axis', 'area'.
+ """
+ ids=spheresFromFile(filename,**kw)
+ mm,mx=aabbExtrema()
+ dim=aabbDim(); axis=dim.index(max(dim))
+ import numpy
+ areas=[approxSectionArea(coord,axis) for coord in
numpy.linspace(mm[axis],mx[axis],num=10)[1:-1]]
+ negIds,posIds=negPosExtremeIds(axis=axis)
+ return {'negIds':negIds,'posIds':posIds,'axis':axis,'area':min(areas)}
+
+
Modified: trunk/gui/py/yadeControl.cpp
===================================================================
--- trunk/gui/py/yadeControl.cpp 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/gui/py/yadeControl.cpp 2009-05-31 10:09:14 UTC (rev 1787)
@@ -316,6 +316,7 @@
NONPOD_ATTRIBUTE_ACCESS(geom,pyInteractionGeometry,interactionGeometry);
NONPOD_ATTRIBUTE_ACCESS(phys,pyInteractionPhysics,interactionPhysics);
/* shorthands */ unsigned id1_get(void){ensureAcc(); return
proxee->getId1();} unsigned id2_get(void){ensureAcc(); return proxee->getId2();}
+ bool isReal_get(void){ensureAcc(); return proxee->isReal(); }
BASIC_PY_PROXY_TAIL;
BASIC_PY_PROXY_HEAD(pyBody,Body)
@@ -807,7 +808,8 @@
.add_property("phys",&pyInteraction::phys_get,&pyInteraction::phys_set)
.add_property("geom",&pyInteraction::geom_get,&pyInteraction::geom_set)
.add_property("id1",&pyInteraction::id1_get)
- .add_property("id2",&pyInteraction::id2_get);
+ .add_property("id2",&pyInteraction::id2_get)
+ .add_property("isReal",&pyInteraction::isReal_get);
BASIC_PY_PROXY_WRAPPER(pyFileGenerator,"Preprocessor")
.def("generate",&pyFileGenerator::generate)
Added: trunk/lib/computational-geometry/Hull2d.hpp
===================================================================
--- trunk/lib/computational-geometry/Hull2d.hpp 2009-05-29 14:35:22 UTC (rev
1786)
+++ trunk/lib/computational-geometry/Hull2d.hpp 2009-05-31 10:09:14 UTC (rev
1787)
@@ -0,0 +1,71 @@
+// 2009 © Václav Šmilauer <[email protected]>
+
+
+/*! Computing convex hull of a 2d cloud of points passed to the constructor,
+ using Graham scan algorithm.
+
+ Use the operator() to launch computation and get the hull as
std::list<Vector2r>.
+
+ The source http://marknelson.us/2007/08/22/convex/ is gratefully
acknowledged.
+ Look there for detailed description and more information.
+*/
+class ConvexHull2d{
+ list<Vector2r> raw_points, lower_partition_points,
upper_partition_points, hull;
+ Vector2r left, right;
+ static Real direction(const Vector2r& p0, const Vector2r& p1, const
Vector2r& p2) {
+ return
((p0[0]-p1[0])*(p2[1]-p1[1]))-((p2[0]-p1[0])*(p0[1]-p1[1]));
+ }
+ struct Vector2r_xComparator{
+ bool operator()(const Vector2r& p1, const Vector2r& p2){return
p1[0]<p2[0];}
+ };
+
+ void partition_points(){
+ raw_points.sort(Vector2r_xComparator());
+ left=raw_points.front(); raw_points.pop_front();
+ right=raw_points.back(); raw_points.pop_back();
+ FOREACH(const Vector2r& p, raw_points){
+ if(direction(left,right,p)<0)
upper_partition_points.push_back(p);
+ else lower_partition_points.push_back(p);
+ }
+ }
+ vector<Vector2r> build_half_hull(list<Vector2r>& in, int factor){
+ vector<Vector2r> out;
+ in.push_back(right); out.push_back(left);
+ while(in.size()>0){
+ out.push_back(in.front()); in.pop_front();
+ while(out.size()>=3){
+ size_t end=out.size()-1;
+
if(factor*direction(out[end-2],out[end],out[end-1])<=0)
out.erase(out.begin()+end-1);
+ else break;
+ }
+ }
+ return out;
+ }
+ public:
+ ConvexHull2d(const list<Vector2r>&
pts){raw_points.assign(pts.begin(),pts.end());};
+ ConvexHull2d(const vector<Vector2r>&
pts){raw_points.assign(pts.begin(),pts.end());};
+ vector<Vector2r> operator()(void){
+ partition_points();
+ vector<Vector2r>
lower_hull=build_half_hull(lower_partition_points,1);
+ vector<Vector2r>
upper_hull=build_half_hull(upper_partition_points,-1);
+ vector<Vector2r> ret;
ret.reserve(lower_hull.size()+upper_hull.size()-2);
+ for(size_t i=upper_hull.size()-1; i>0; i--)
ret.push_back(upper_hull[i]);
+ size_t lsize=lower_hull.size();
+ for(size_t i=0; i<lsize-1; i++) ret.push_back(lower_hull[i]);
+ return ret;
+ }
+};
+
+/*! Compute area of a simple 2d polygon, using the Surveyor's formula.
+ http://en.wikipedia.org/wiki/Polygon
+
+ The first and last points shouldn't be given twice.
+*/
+Real simplePolygonArea2d(vector<Vector2r> P){
+ Real ret=0.; size_t n=P.size();
+ for(size_t i=0; i<n-1; i++) { ret+=P[i][0]*P[i+1][1]-P[i+1][0]*P[i][1];}
+ ret+=P[n-1][0]*P[0][1]-P[0][0]*P[n-1][1];
+ return abs(ret/2.);
+}
+
+
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionDispatchers.cpp
2009-05-31 10:09:14 UTC (rev 1787)
@@ -53,8 +53,12 @@
const shared_ptr<Body>&
b2=Body::byId(I->getId2(),rootBody);
assert(I->functorCache.geom);
+ bool wasReal=I->isReal();
bool
geomCreated=I->functorCache.geom->go(b1->interactingGeometry,b2->interactingGeometry,b1->physicalParameters->se3,
b2->physicalParameters->se3,I);
- if(!geomCreated) continue;
+ if(!geomCreated){
+ if(wasReal)
rootBody->interactions->requestErase(I->getId1(),I->getId2()); // fully created
interaction without geometry is reset and perhaps erased in the next step
+ continue; // in any case don't care about this
one anymore
+ }
// InteractionPhysicsMetaEngine
if(!I->functorCache.phys){
@@ -80,10 +84,14 @@
const shared_ptr<Body>&
b1=Body::byId(I->getId1(),rootBody);
const shared_ptr<Body>&
b2=Body::byId(I->getId2(),rootBody);
// InteractionGeometryMetaEngine
+ bool wasReal=I->isReal();
bool geomCreated =
b1->interactingGeometry &&
b2->interactingGeometry && // some bodies do not have interactingGeometry
geomDispatcher->operator()(b1->interactingGeometry, b2->interactingGeometry,
b1->physicalParameters->se3, b2->physicalParameters->se3,I);
- if(!geomCreated) continue;
+ if(!geomCreated){
+ if(wasReal)
*rootBody->interactions->requestErase(I->getId1(),I->getId2());
+ continue;
+ }
// InteractionPhysicsMetaEngine
// geom may have swapped bodies, get bodies again
physDispatcher->operator()(Body::byId(I->getId1(),rootBody)->physicalParameters,
Body::byId(I->getId2(),rootBody)->physicalParameters,I);
Modified: trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/pkg/common/Engine/MetaEngine/InteractionGeometryMetaEngine.cpp
2009-05-31 10:09:14 UTC (rev 1787)
@@ -58,14 +58,17 @@
const long size=ncb->transientInteractions->size();
#pragma omp parallel for
for(long i=0; i<size; i++){
- const shared_ptr<Interaction>&
interaction=(*ncb->transientInteractions)[i];
+ const shared_ptr<Interaction>&
I=(*ncb->transientInteractions)[i];
#else
- FOREACH(const shared_ptr<Interaction>& interaction,
*ncb->interactions){
+ FOREACH(const shared_ptr<Interaction>& I, *ncb->interactions){
#endif
- const shared_ptr<Body>&
b1=(*bodies)[interaction->getId1()];
- const shared_ptr<Body>&
b2=(*bodies)[interaction->getId2()];
- b1->interactingGeometry && b2->interactingGeometry &&
// some bodies do not have interactingGeometry
- operator()(b1->interactingGeometry,
b2->interactingGeometry, b1->physicalParameters->se3,
b2->physicalParameters->se3, interaction);
+ const shared_ptr<Body>& b1=(*bodies)[I->getId1()];
+ const shared_ptr<Body>& b2=(*bodies)[I->getId2()];
+ bool wasReal=I->isReal();
+ if (!b1->interactingGeometry ||
!b2->interactingGeometry) { assert(!wasReal); continue; } // some bodies do not
have interactingGeometry
+ bool geomCreated=operator()(b1->interactingGeometry,
b2->interactingGeometry, b1->physicalParameters->se3,
b2->physicalParameters->se3, I);
+ // reset && erase interaction that existed but now has
no geometry anymore
+ if(wasReal && !geomCreated){
ncb->interactions->requestErase(I->getId1(),I->getId2()); }
}
}
Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/pkg/dem/ConcretePM.cpp 2009-05-31 10:09:14 UTC (rev 1787)
@@ -164,14 +164,14 @@
CPM_MATERIAL_MODEL
#else
// very simplified version of the constitutive law
- kappaD=max(max(0,epsN),kappaD); // internal variable, max
positive strain (non-decreasing)
+ kappaD=max(max(0.,epsN),kappaD); // internal variable, max
positive strain (non-decreasing)
omega=isCohesive?funcG(kappaD,epsCrackOnset,epsFracture,neverDamage):1.; //
damage variable (non-decreasing, as funcG is also non-decreasing)
sigmaN=(1-(epsN>0?omega:0))*E*epsN; // damage taken in account
in tension only
sigmaT=G*epsT; // trial stress
Real
yieldSigmaT=max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle);
// Mohr-Coulomb law with damage
if(sigmaT.SquaredLength()>yieldSigmaT*yieldSigmaT){
sigmaT*=yieldSigmaT/sigmaT.Length(); // stress return
- epsPlSum+=rT*contGeom->slipToStrainTMax(rT/G); //
adjust strain
+
epsPlSum+=yieldSigmaT*contGeom->slipToStrainTMax(yieldSigmaT/G); // adjust
strain
}
relResidualStrength=isCohesive?(kappaD<epsCrackOnset?1.:(1-omega)*(kappaD)/epsCrackOnset):0;
#endif
Modified: trunk/pkg/dem/ConcretePM.hpp
===================================================================
--- trunk/pkg/dem/ConcretePM.hpp 2009-05-29 14:35:22 UTC (rev 1786)
+++ trunk/pkg/dem/ConcretePM.hpp 2009-05-31 10:09:14 UTC (rev 1787)
@@ -210,7 +210,7 @@
// init to signaling_NaN to force crash if not
initialized (better than unknowingly using garbage values)
sigmaT=epsCrackOnset=relDuctility=G_over_E=std::numeric_limits<Real>::signaling_NaN();
neverDamage=false;
- cohesiveThresholdIter=-1;
+ cohesiveThresholdIter=10;
dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
isoPrestress=0;
}
Modified: trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp 2009-05-29
14:35:22 UTC (rev 1786)
+++ trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp 2009-05-31
10:09:14 UTC (rev 1787)
@@ -81,7 +81,7 @@
YADE_CAST<ParticleParameters*>(b->physicalParameters.get())->velocity[axis]=pNormalized*(v1-v0)+v0;
}
}
-
stressUpdateInterval=max(1,(int)(2e-5/(abs(strainRate)*Omega::instance().getTimeStep())));
+
stressUpdateInterval=max(1,(int)(1e-5/(abs(strainRate)*Omega::instance().getTimeStep())));
LOG_INFO("Stress will be updated every "<<stressUpdateInterval<<"
steps.");
/* if we have default (<0) crossSectionArea, try to get it from root's
AABB;
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
