Author: eudoxos
Date: 2009-03-22 15:01:25 +0100 (Sun, 22 Mar 2009)
New Revision: 1726
Added:
trunk/scripts/test/shear.py
trunk/scripts/test/triax-identical-results.py
Modified:
trunk/SConstruct
trunk/core/BexContainer.hpp
trunk/core/Omega.cpp
trunk/core/yade.cpp
trunk/gui/py/PythonUI_rc.py
trunk/lib/factory/ClassFactory.cpp
trunk/lib/factory/ClassFactory.hpp
trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.hpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.hpp
Log:
1. Remove logger from ClassFactory
2. Cleanup logging stuff in main (use constructor function, among other)
3. Cleanup system exit from python in PythonUI_rc.py
4. Add SpheresContactGeometry::updateShear that can be optionally used with
ElasticContactLaw to update shear displacement instead of updating shearForce.
triax-identical-results.py show quite large difference between both
implementations, but I am not able to tell which one is correct.
scripts/test/shear.py shown almost no difference for 2-sphere scenarios, modulo
differences at 15th decimal place or so.
5. Remove debug output from BexContainer
6. Remove warning about meniscus data from CapillaryCohesiveLaw (warning is
given in postProcessAttributes now, i.e. iff the class is actually used)
7. Add logger to snow.
8. Removed shear computation code in ElasticContactLaw, use
SpheresContactGeometry::updateShearForce.
9. Fix scons deprecation warnings
Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/SConstruct 2009-03-22 14:01:25 UTC (rev 1726)
@@ -66,8 +66,8 @@
env=Environment(tools=['default'])
profileFile='scons.current-profile'
-profOpts=Options(profileFile)
-profOpts.AddOptions(('profile','Config profile to use (predefined: default or
"", opt); append ! to use it but not save for next build (in
scons.current-profile)','default'))
+profOpts=Variables(profileFile)
+profOpts.Add(('profile','Config profile to use (predefined: default or "",
opt); append ! to use it but not save for next build (in
scons.current-profile)','default'))
profOpts.Update(env)
# multiple profiles - run them all at the same time
# take care not to save current profile for those parallel builds
@@ -107,7 +107,7 @@
else: defOptions={'debug':0,'optimize':0,'variant':profile,'openmp':True}
-opts=Options(optsFile)
+opts=Variables(optsFile)
#
# The convention now is, that
# 1. CAPITALIZED options are
@@ -115,19 +115,19 @@
# (b) c preprocessor macros available to the program source (like PREFIX and
SUFFIX)
# 2. lowercase options influence the building process, compiler options and
the like.
#
-opts.AddOptions(
+opts.AddVariables(
### OLD: use PathOption with PathOption.PathIsDirCreate, but that
doesn't exist in 0.96.1!
('PREFIX','Install path prefix','/usr/local'),
('runtimePREFIX','Runtime path prefix; DO NOT USE, inteded for
packaging only.','$PREFIX'),
('variant','Build variant, will be suffixed to all files, along with
version (beware: if PREFIX is the same, headers of the older version will still
be overwritten',defOptions['variant'],None,lambda x:x),
- BoolOption('debug', 'Enable debugging information and disable
optimizations',defOptions['debug']),
- BoolOption('gprof','Enable profiling information for gprof',0),
- BoolOption('optimize','Turn on heavy
optimizations',defOptions['optimize']),
- BoolOption('openmp','Compile with openMP parallelization
support',defOptions['openmp']),
- ListOption('exclude','Yade components that will not be
built','none',names=['qt3','gui','extra','common','dem','fem','lattice','mass-spring','realtime-rigidbody','snow']),
- EnumOption('arcs','Whether to generate or use branch
probabilities','',['','gen','use'],{'no':'','0':'','false':''},1),
+ BoolVariable('debug', 'Enable debugging information and disable
optimizations',defOptions['debug']),
+ BoolVariable('gprof','Enable profiling information for gprof',0),
+ BoolVariable('optimize','Turn on heavy
optimizations',defOptions['optimize']),
+ BoolVariable('openmp','Compile with openMP parallelization
support',defOptions['openmp']),
+ ListVariable('exclude','Yade components that will not be
built','none',names=['qt3','gui','extra','common','dem','fem','lattice','mass-spring','realtime-rigidbody','snow']),
+ EnumVariable('arcs','Whether to generate or use branch
probabilities','',['','gen','use'],{'no':'','0':'','false':''},1),
# OK, dummy prevents bug in scons: if one selects all, it says all in
scons.config, but without quotes, which generates error.
- ListOption('features','Optional features that are turned
on','python,log4cxx',names=['python','log4cxx','binfmt','dummy']),
+ ListVariable('features','Optional features that are turned
on','python,log4cxx',names=['python','log4cxx','binfmt','dummy']),
('jobs','Number of jobs to run at the same time (same as -j, but
saved)',4,None,int),
('extraModules', 'Extra directories with their own SConscript files
(must be in-tree) (whitespace separated)',None,None,Split),
('buildPrefix','Where to create build-[version][variant] directory for
intermediary files','..'),
@@ -140,10 +140,10 @@
('march','Architecture to use with -march=... when
optimizing','native',None,None),
#('SHLINK','Linker for shared objects','g++'),
('SHCCFLAGS','Additional compiler flags for linking (for
plugins).',None,None,Split),
- BoolOption('QUAD_PRECISION','typedef Real as long double (=quad)',0),
- BoolOption('pretty',"Don't show compiler command line (like the Linux
kernel)",1),
- BoolOption('useMiniWm3','use local miniWm3 library instead of
Wm3Foundation',1),
- #BoolOption('useLocalQGLViewer','use in-tree QGLViewer library instead
of the one installed in system',1),
+ BoolVariable('QUAD_PRECISION','typedef Real as long double (=quad)',0),
+ BoolVariable('pretty',"Don't show compiler command line (like the Linux
kernel)",1),
+ BoolVariable('useMiniWm3','use local miniWm3 library instead of
Wm3Foundation',1),
+ #BoolVariable('useLocalQGLViewer','use in-tree QGLViewer library
instead of the one installed in system',1),
)
opts.Update(env)
opts.Save(optsFile,env)
Modified: trunk/core/BexContainer.hpp
===================================================================
--- trunk/core/BexContainer.hpp 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/core/BexContainer.hpp 2009-03-22 14:01:25 UTC (rev 1726)
@@ -122,7 +122,7 @@
_force.resize(newSize);
_torque.resize(newSize);
size=newSize;
- std::cerr<<"[DEBUG] BexContainer: Resized to
"<<size<<std::endl;
+ // std::cerr<<"[DEBUG] BexContainer: Resized to
"<<size<<std::endl;
}
};
Modified: trunk/core/Omega.cpp
===================================================================
--- trunk/core/Omega.cpp 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/core/Omega.cpp 2009-03-22 14:01:25 UTC (rev 1726)
@@ -46,7 +46,7 @@
CREATE_LOGGER(Omega);
Omega::Omega(){
- if(getenv("YADE_DEBUG")) cerr<<"Constructing Omega; _must_ be only
once, otherwise linking is broken (missing -rdynamic?)\n";
+ if(getenv("YADE_DEBUG")) cerr<<"Constructing Omega; _must_ be only
once, otherwise linking is broken (missing -rdynamic?)"<<endl;
}
Omega::~Omega(){LOG_INFO("Shuting down; duration
"<<(microsec_clock::local_time()-msStartingSimulationTime)/1000<<" s");}
Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/core/yade.cpp 2009-03-22 14:01:25 UTC (rev 1726)
@@ -34,11 +34,29 @@
#ifdef LOG4CXX
// provides parent logger for everybody
log4cxx::LoggerPtr logger=log4cxx::Logger::getLogger("yade");
+
+ #ifdef LOG4CXX_TRACE
+ log4cxx::LevelPtr debugLevel=log4cxx::Level::getDebug(),
infoLevel=log4cxx::Level::getInfo(), warnLevel=log4cxx::Level::getWarn();
+ #else
+ log4cxx::LevelPtr debugLevel=log4cxx::Level::DEBUG,
infoLevel=log4cxx::Level::INFO, warnLevel=log4cxx::Level::WARN;
+ #endif
+
+ /* initialization of log4cxx for early logging */
+ __attribute__((constructor(1000))) void initLog4cxx(){
+ log4cxx::BasicConfigurator::configure();
+ log4cxx::LoggerPtr
localLogger=log4cxx::Logger::getLogger("yade");
+ if(getenv("YADE_DEBUG")){
+ LOG_INFO("YADE_DEBUG environment variable is defined,
logging level is DEBUG.");
+ localLogger->setLevel(debugLevel);
+ }
+ else localLogger->setLevel(warnLevel);
+ }
#endif
void nullHandler(int sig){}
+void termHandler(int sig){cerr<<"terminating..."<<endl; raise(SIGTERM);}
void warnOnceHandler(int sig){
- cerr<<"WARN: nullHandler (probably log4cxx error), signal
"<<(sig==SIGSEGV?"SEGV":"[other]. Signal will be ignored since now.")<<endl;
+ cerr<<"WARN: nullHandler (probably log4cxx error), signal
"<<(sig==SIGSEGV?"SEGV":"[other]")<<". Signal will be ignored since now."<<endl;
signal(sig,nullHandler);
}
@@ -185,22 +203,19 @@
optind=0; opterr=0;
#ifdef LOG4CXX
- #ifdef LOG4CXX_TRACE
- log4cxx::LevelPtr
debugLevel=log4cxx::Level::getDebug(), infoLevel=log4cxx::Level::getInfo(),
warnLevel=log4cxx::Level::getWarn();
- #else
- log4cxx::LevelPtr debugLevel=log4cxx::Level::DEBUG,
infoLevel=log4cxx::Level::INFO, warnLevel=log4cxx::Level::WARN;
- #endif
// read logging configuration from file and watch it (creates a
separate thread)
std::string logConf=configPath+"/logging.conf";
if(filesystem::exists(logConf)){
log4cxx::PropertyConfigurator::configure(logConf);
LOG_INFO("Loading "<<logConf);
} else { // otherwise use simple console-directed logging
- log4cxx::BasicConfigurator::configure();
logger->setLevel(warnLevel);
LOG_INFO("Logger uses basic (console) configuration
since `"<<logConf<<"' was not found. INFO and DEBUG messages will be omitted.");
LOG_INFO("Look at the file doc/logging.conf.sample in
the source distribution as an example on how to customize logging.");
}
+ // command-line switches override levels
+ if(verbose==1) logger->setLevel(infoLevel);
+ else if (verbose>=2) logger->setLevel(debugLevel);
#endif
Omega::instance().yadeVersionName = "Yet Another Dynamic Engine 0.12.x,
beta, SVN snapshot.";
@@ -209,16 +224,7 @@
filesystem::path yadeConfigPath =
filesystem::path(Omega::instance().yadeConfigPath, filesystem::native);
filesystem::path yadeConfigFile =
filesystem::path(Omega::instance().yadeConfigPath + "/preferences.xml",
filesystem::native);
- #ifdef LOG4CXX
- if(verbose==1) logger->setLevel(infoLevel);
- else if (verbose>=2) logger->setLevel(debugLevel);
- if(getenv("YADE_DEBUG")){
- LOG_INFO("YADE_DEBUG environment variable is defined,
setting logging level to DEBUG.");
- logger->setLevel(debugLevel);
- }
- #endif
-
#ifdef EMBED_PYTHON
/* see http://www.python.org/dev/peps/pep-0311 for threading
with Python embedded */
LOG_DEBUG("Initializing Python...");
@@ -288,8 +294,8 @@
#ifdef YADE_DEBUG
if(useGdb){
signal(SIGABRT,SIG_DFL); signal(SIGHUP,SIG_DFL); //
default handlers
- signal(SIGSEGV,warnOnceHandler); // FIXME: this is to
cover up crash that occurs in log4cxx on i386 sometimes
unlink(Omega::instance().gdbCrashBatch.c_str());
+ signal(SIGSEGV,termHandler);
}
#endif
Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/gui/py/PythonUI_rc.py 2009-03-22 14:01:25 UTC (rev 1726)
@@ -23,6 +23,14 @@
def _quit(): import sys; sys.exit(0)
__builtins__.quit=_quit
+def cleanup():
+ try:
+ import yade.qt, time
+ yade.qt.close()
+ while True: time.sleep(.1) # wait to be killed
+ except ImportError: pass
+
+
## run the TCP server
import yade.PythonTCPServer
srv=yade.PythonTCPServer.PythonTCPServer(minPort=9000)
@@ -31,50 +39,54 @@
runtime.argv=[runtime.script]+runtime.argv
sys.argv=runtime.argv # could be [] as well
-## run simulation if requested from the command line
-if runtime.simulation:
- print "Running simulation "+runtime.simulation
- o=Omega(); o.load(runtime.simulation); o.run();
+try:
+ ## run simulation if requested from the command line
+ if runtime.simulation:
+ print "Running simulation "+runtime.simulation
+ o=Omega(); o.load(runtime.simulation); o.run();
-## run script if requested from the command line
-if runtime.script:
- print "Running script "+runtime.script
- # an exception from python would propagate to c++ unhandled and cause
crash
- try: execfile(runtime.script)
- except SystemExit: raise # re-raise sys.exit
- except: # all other exceptions
- import traceback
- traceback.print_exc()
- if(runtime.nonInteractive or runtime.stopAfter): sys.exit(1)
-if runtime.stopAfter:
- print "Finished, bye."
- sys.exit(0)
+ ## run script if requested from the command line
+ if runtime.script:
+ print "Running script "+runtime.script
+ # an exception from python would propagate to c++ unhandled and
cause crash
+ try:
+ execfile(runtime.script)
+ except SystemExit: raise
+ except: # all other exceptions
+ import traceback
+ traceback.print_exc()
+ if(runtime.nonInteractive or runtime.stopAfter):
sys.exit(1)
+ if runtime.stopAfter:
+ print "Finished, bye."
+ sys.exit(0)
-# run commands if requested from the command line
-#if yadeRunCommands:
-# print "Running commands from commandline: "+yadeRunCommands
-# exec(yadeRunCommands)
+ # run commands if requested from the command line
+ #if yadeRunCommands:
+ # print "Running commands from commandline: "+yadeRunCommands
+ # exec(yadeRunCommands)
-if runtime.nonInteractive:
- import time;
- while True: time.sleep(1)
-else:
- sys.argv[0]='<embedded python interpreter>'
- from IPython.Shell import IPShellEmbed
- ipshell = IPShellEmbed(banner=r"""__ __ ____ ____
_
-\ \ / /_ _| _ \ ___ / ___|___ _ __ ___ ___ | | ___
- \ V / _` | | | |/ _ \ | | / _ \| '_ \/ __|/ _ \| |/ _ \
- | | (_| | |_| | __/ | |__| (_) | | | \__ \ (_) | | __/
- |_|\__,_|____/ \___| \____\___/|_| |_|___/\___/|_|\___|
-""",exit_msg='Bye.'
-
,rc_override={'execfile':[runtime.prefix+'/lib/yade'+runtime.suffix+'/gui/yade/ipython.py']})
+ if runtime.nonInteractive:
+ import time;
+ while True: time.sleep(1)
+ else:
+ sys.argv[0]='<embedded python interpreter>'
+ from IPython.Shell import IPShellEmbed
+ ipshell = IPShellEmbed(banner=r"""__ __ ____ ____
_
+ \ \ / /_ _| _ \ ___ / ___|___ _ __ ___ ___ | | ___
+ \ V / _` | | | |/ _ \ | | / _ \| '_ \/ __|/ _ \| |/ _ \
+ | | (_| | |_| | __/ | |__| (_) | | | \__ \ (_) | | __/
+ |_|\__,_|____/ \___| \____\___/|_| |_|___/\___/|_|\___|
+ """,exit_msg='Bye.'
+
,rc_override={'execfile':[runtime.prefix+'/lib/yade'+runtime.suffix+'/gui/yade/ipython.py']})
- ipshell()
- # save history -- a workaround for atexit handlers not being run (why?)
- #
http://lists.ipython.scipy.org/pipermail/ipython-user/2008-September/005839.html
- import IPython.ipapi
- IPython.ipapi.get().IP.atexit_operations()
- try:
- import yade.qt
- yade.qt.close()
- except ImportError: pass
+ ipshell()
+ # save history -- a workaround for atexit handlers not being
run (why?)
+ #
http://lists.ipython.scipy.org/pipermail/ipython-user/2008-September/005839.html
+ import IPython.ipapi
+ IPython.ipapi.get().IP.atexit_operations()
+except SystemExit:
+ cleanup()
+finally:
+ cleanup()
+
+
Modified: trunk/lib/factory/ClassFactory.cpp
===================================================================
--- trunk/lib/factory/ClassFactory.cpp 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/lib/factory/ClassFactory.cpp 2009-03-22 14:01:25 UTC (rev 1726)
@@ -12,8 +12,6 @@
#include<boost/algorithm/string/regex.hpp>
-CREATE_LOGGER(ClassFactory);
-
class Factorable;
void ClassFactory::addBaseDirectory(const string& dir)
@@ -128,16 +126,17 @@
void ClassFactory::registerPluginClasses(const char* fileAndClasses[]){
assert(fileAndClasses[0]!=NULL); // must be file name
+ bool log=getenv("YADE_DEBUG");
// only filename given, no classes names explicitly
if(fileAndClasses[1]==NULL){
/* strip leading path (if any; using / as path separator) and
strip one suffix (if any) to get the contained class name */
string
heldClass=boost::algorithm::replace_regex_copy(string(fileAndClasses[0]),boost::regex("^(.*/)?(.*?)(\\.[^.]*)?$"),string("\\2"));
- LOG_DEBUG("Plugin "<<fileAndClasses[0]<<", class
"<<heldClass<<" (deduced)");
+ if(log) cerr<<"Plugin "<<fileAndClasses[0]<<", class
"<<heldClass<<" (deduced)"<<endl;
pluginClasses.push_back(heldClass); // last item with
everything up to last / take off and .suffix strip
}
else {
for(int i=1; fileAndClasses[i]!=NULL; i++){
- LOG_DEBUG("Plugin "<<fileAndClasses[0]<<", class
"<<fileAndClasses[i]);
+ if(log) cerr<<"Plugin "<<fileAndClasses[0]<<", class
"<<fileAndClasses[i]<<endl;
pluginClasses.push_back(fileAndClasses[i]);
}
}
Modified: trunk/lib/factory/ClassFactory.hpp
===================================================================
--- trunk/lib/factory/ClassFactory.hpp 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/lib/factory/ClassFactory.hpp 2009-03-22 14:01:25 UTC (rev 1726)
@@ -26,7 +26,6 @@
#include<yade/lib-loki/Singleton.hpp>
-#include<yade/lib-base/Logging.hpp>
#include "FactoryExceptions.hpp"
@@ -105,7 +104,7 @@
/// Map that contains the name of the registered class and
their description
FactorableCreatorsMap map;
- ClassFactory() { if(getenv("YADE_DEBUG")) cerr<<"Constructing
ClassFactory; _must_ be only once, otherwise linking is broken (missing
-rdynamic?)\n"; };
+ ClassFactory() { if(getenv("YADE_DEBUG")) cerr<<"Constructing
ClassFactory; _must_ be only once, otherwise linking is broken (missing
-rdynamic?)\n"<<endl; };
ClassFactory(const ClassFactory&);
ClassFactory& operator=(const ClassFactory&);
virtual ~ClassFactory() {};
@@ -152,8 +151,6 @@
virtual string getClassName() const { return "Factorable"; };
virtual string getBaseClassName(int ) const { return "";};
- DECLARE_LOGGER;
-
FRIEND_SINGLETON(ClassFactory);
};
Modified:
trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/pkg/common/DataClass/InteractionPhysics/NormalShearInteractions.hpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -11,15 +11,10 @@
Real kn;
//! normal force
Vector3r normalForce;
- NormalInteraction(){createIndex(); }
+ NormalInteraction(): normalForce(Vector3r::ZERO)
{createIndex(); }
virtual ~NormalInteraction();
- protected:
- virtual void registerAttributes(){
- REGISTER_ATTRIBUTE(kn);
- REGISTER_ATTRIBUTE(normalForce);
- }
- REGISTER_CLASS_NAME(NormalInteraction);
- REGISTER_BASE_CLASS_NAME(InteractionPhysics);
+ REGISTER_ATTRIBUTES(/*no base class attributes*/,(kn)(normalForce));
+ REGISTER_CLASS_AND_BASE(NormalInteraction,InteractionPhysics);
REGISTER_CLASS_INDEX(NormalInteraction,InteractionPhysics);
};
REGISTER_SERIALIZABLE(NormalInteraction);
@@ -33,16 +28,10 @@
Real ks;
//! shear force
Vector3r shearForce;
- NormalShearInteraction(){ createIndex(); }
+ NormalShearInteraction(): shearForce(Vector3r::ZERO){
createIndex(); }
virtual ~NormalShearInteraction();
- protected:
- virtual void registerAttributes(){
- NormalInteraction::registerAttributes();
- REGISTER_ATTRIBUTE(ks);
- REGISTER_ATTRIBUTE(shearForce);
- }
- REGISTER_CLASS_NAME(NormalShearInteraction);
- REGISTER_BASE_CLASS_NAME(NormalInteraction);
+ REGISTER_ATTRIBUTES(NormalInteraction,(ks)(shearForce));
+ REGISTER_CLASS_AND_BASE(NormalShearInteraction,NormalInteraction);
REGISTER_CLASS_INDEX(NormalShearInteraction,NormalInteraction);
};
REGISTER_SERIALIZABLE(NormalShearInteraction);
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -9,6 +9,63 @@
// At least one virtual method must be in the .cpp file (!!!)
SpheresContactGeometry::~SpheresContactGeometry(){};
+#ifdef SCG_SHEAR
+void SpheresContactGeometry::updateShear(const RigidBodyParameters* rbp1,
const RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
+
+ Vector3r axis;
+ Real angle;
+
+ shearIncrement=Vector3r::ZERO;
+
+ // approximated rotations
+ axis = prevNormal.Cross(normal);
+ shearIncrement -= shear.Cross(axis);
+ angle = dt*0.5*normal.Dot(rbp1->angularVelocity +
rbp2->angularVelocity);
+ axis = angle*normal;
+ shearIncrement -= (shear+shearIncrement).Cross(axis);
+
+ // exact rotations (not adapted to shear/shearIncrement!)
+ #if 0
+ Quaternionr q;
+ axis =
prevNormal.Cross(normal);
+ angle =
acos(normal.Dot(prevNormal));
+ q.FromAngleAxis(angle,axis);
+ shearForce = shearForce*q;
+ angle =
dt*0.5*normal.dot(rbp1->angularVelocity+rbp2->angularVelocity);
+ axis = normal;
+ q.FromAngleAxis(angle,axis);
+ shearForce = q*shearForce;
+ #endif
+
+ Vector3r& x = contactPoint;
+ Vector3r c1x, c2x;
+
+ if(avoidGranularRatcheting){
+ /* The following definition of c1x and c2x is to avoid
"granular ratcheting"
+ * (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
+ * "Micro-mechanical investigation of granular ratcheting, in
Cyclic Behaviour of Soils and Liquefaction Phenomena",
+ * ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 -
and a lot more papers from the same authors) */
+
+ // FIXME: For sphere-facet contact this will give an erroneous
value of relative velocity...
+ c1x = radius1*normal;
+ c2x = -radius2*normal;
+ }
+ else {
+ // FIXME: It is correct for sphere-sphere and sphere-facet
contact
+ c1x = (x - rbp1->se3.position);
+ c2x = (x - rbp2->se3.position);
+ }
+
+ Vector3r relativeVelocity =
(rbp2->velocity+rbp2->angularVelocity.Cross(c2x)) -
(rbp1->velocity+rbp1->angularVelocity.Cross(c1x));
+ Vector3r shearVelocity =
relativeVelocity-normal.Dot(relativeVelocity)*normal;
+ Vector3r shearDisplacement = shearVelocity*dt;
+ shearIncrement -= shearDisplacement;
+
+ shear+=shearIncrement;
+ shearUpdateIter=Omega::instance().getCurrentIteration();
+}
+#endif
+
void SpheresContactGeometry::updateShearForce(Vector3r& shearForce, Real ks,
const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting){
Vector3r axis;
Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -8,31 +8,31 @@
#include<yade/pkg-common/RigidBodyParameters.hpp>
/*! Class representing geometry of two spheres in contact.
*
- * exactRot code
- * =============
- * At initial contact, each of the two spheres fixes a point on its surface
relative to its own orientation.
- * It is therefore possible to derive at any later point how much slipping
between the two spheres has taken
- * place since the first contact.
- *
- * The surface point is stored as quaternion.
- *
- * Function is provided to manipulate those points:
- *
- * (a) plastic slip, when we want to limit their maximum distance (in the
projected plane)
- * (b) rolling, where those points must be relocated to not flip over the π
angle from the current contact point
- *
- * TODO: add torsion code, that will (if a flag is on) compute torsion angle
on the contact axis.
- *
- * TODO: add bending code, that will compute bending angle of the contact axis
- *
- *
+ * The code under SCG_SHEAR is experimental and is used only if
ElasticContactLaw::useShear is explicitly true
*/
+
+#define SCG_SHEAR
+
class SpheresContactGeometry: public InteractionGeometry{
public:
Vector3r normal, // unit vector in the direction from sphere1
center to sphere2 center
contactPoint;
Real radius1,radius2,penetrationDepth;
+ #ifdef SCG_SHEAR
+ //! Total value of the current shear. Update the value
using updateShear.
+ Vector3r shear;
+ //! Normal of the contact in the previous step
+ Vector3r prevNormal;
+ //! Increment of shear displacement in last step (is
already added to shear); perhaps useful for viscosity or something similar
+ Vector3r shearIncrement;
+ long shearUpdateIter; // debugging only, to check when
shear was updated last time
+ //! update shear on this contact given dynamic
parameters of bodies. Should be called from constitutive law, exactly once per
iteration
+ void updateShear(const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
+ //const Vector3r& getShear(){
if(Omega::instance().getCurrentIteration()>shearUpdateIter) throw
runtime_error("SpheresContactGeometry::updateShear must be called prior to
getShear()."); return shear; }
+ #endif
+
+ // all the rest here will hopefully disappear at some point...
// begin abusive storage
bool hasNormalViscosity;
@@ -97,7 +97,11 @@
Vector3r relRotVector() const;
-
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();}
+
SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();
+ #ifdef SCG_SHEAR
+ shear=Vector3r::ZERO; prevNormal=Vector3r::ZERO
/*initialized to aproper value by geom functor*/;
shearIncrement=Vector3r::ZERO; shearUpdateIter=-1 /* proper value from geom
functor */;
+ #endif
+ }
virtual ~SpheresContactGeometry();
void updateShearForce(Vector3r& shearForce, Real ks, const
Vector3r& prevNormal, const RigidBodyParameters* rbp1, const
RigidBodyParameters* rbp2, Real dt, bool avoidGranularRatcheting=true);
@@ -107,6 +111,12 @@
(contactPoint)
(radius1)
(radius2)
+ #ifdef SCG_SHEAR
+ (shear)
+ (prevNormal)
+ (shearIncrement)
+ (shearUpdateIter)
+ #endif
(facetContactFace)
// hasShear
(hasShear)
Modified:
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
---
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
2009-03-19 17:55:37 UTC (rev 1725)
+++
trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -37,6 +37,17 @@
if(c->interactionGeometry)
scm=YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
else { scm=shared_ptr<SpheresContactGeometry>(new
SpheresContactGeometry()); c->interactionGeometry=scm; }
+ #ifdef SCG_SHEAR
+ if(c->isNew){
+ scm->prevNormal=normal;
+
scm->shearUpdateIter=Omega::instance().getCurrentIteration(); /* no shear at
the very beginning; shear initialized to zero vector in SCG ctor */
+ } else {
+ scm->prevNormal=scm->normal;
+ // make sure updateShear was properly called at
last iteration; debugging only
+
//assert(scm->shearUpdateIter==Omega::instance().getCurrentIteration()-1);
+ }
+ #endif
+
Real penetrationDepth=s1->radius+s2->radius-normal.Normalize();
/* Normalize() works in-place and returns length before normalization; from
here, normal is unit vector */
scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
scm->normal=normal;
Modified: trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.cpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -39,26 +39,30 @@
{
sdecGroupMask=1;
- capillary = shared_ptr<capillarylaw>(new capillarylaw); // ????????
-
- capillary->fill("M(r=1)");
- capillary->fill("M(r=1.1)");
- capillary->fill("M(r=1.25)");
- capillary->fill("M(r=1.5)");
- capillary->fill("M(r=1.75)");
- capillary->fill("M(r=2)");
- capillary->fill("M(r=3)");
- capillary->fill("M(r=4)");
- capillary->fill("M(r=5)");
- capillary->fill("M(r=10)");
-
CapillaryPressure=0;
fusionDetection = false;
binaryFusion = true;
+ // capillary setup moved to postProcessAttributes
+
}
+void CapillaryCohesiveLaw::postProcessAttributes(bool deserializing){
+ if(!deserializing) return;
+ capillary = shared_ptr<capillarylaw>(new capillarylaw); // ????????
+ capillary->fill("M(r=1)");
+ capillary->fill("M(r=1.1)");
+ capillary->fill("M(r=1.25)");
+ capillary->fill("M(r=1.5)");
+ capillary->fill("M(r=1.75)");
+ capillary->fill("M(r=2)");
+ capillary->fill("M(r=3)");
+ capillary->fill("M(r=4)");
+ capillary->fill("M(r=5)");
+ capillary->fill("M(r=10)");
+}
+
void CapillaryCohesiveLaw::registerAttributes()
{
InteractionSolver::registerAttributes();
Modified: trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.hpp
2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/Engine/StandAloneEngine/CapillaryCohesiveLaw.hpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -89,6 +89,7 @@
protected :
void registerAttributes();
+ virtual void postProcessAttributes(bool deserializing);
NEEDS_BEX("Force","Momentum");
REGISTER_CLASS_NAME(CapillaryCohesiveLaw);
REGISTER_BASE_CLASS_NAME(InteractionSolver);
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-19
17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp 2009-03-22
14:01:25 UTC (rev 1726)
@@ -64,6 +64,9 @@
momentRotationLaw = true;
actionForceIndex = actionForce->getClassIndex();
actionMomentumIndex = actionMomentum->getClassIndex();
+ #ifdef SCG_SHEAR
+ useShear=false;
+ #endif
}
@@ -72,6 +75,9 @@
InteractionSolver::registerAttributes();
REGISTER_ATTRIBUTE(sdecGroupMask);
REGISTER_ATTRIBUTE(momentRotationLaw);
+ #ifdef SCG_SHEAR
+ REGISTER_ATTRIBUTE(useShear);
+ #endif
}
@@ -111,62 +117,30 @@
Real un=currentContactGeometry->penetrationDepth;
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real)
0)*currentContactGeometry->normal;
- #if 0
- // the core under #else, refactored
-
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,isDynamic1,isDynamic2,dt);
- #else
- Vector3r axis;
- Real angle;
-
- /// Here is the code with approximated rotations ///
- axis =
currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal);
- shearForce -= shearForce.Cross(axis);
- angle =
dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
- axis = angle*currentContactGeometry->normal;
- shearForce -= shearForce.Cross(axis);
-
- /// Here is the code with exact rotations ///
-
- // Quaternionr q;
- //
- // axis =
currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
- // angle =
acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
- // q.fromAngleAxis(angle,axis);
- //
- // currentContactPhysics->shearForce =
currentContactPhysics->shearForce*q;
- //
- // angle =
dt*0.5*currentContactGeometry->normal.dot(de1->angularVelocity+de2->angularVelocity);
- // axis =
currentContactGeometry->normal;
- // q.fromAngleAxis(angle,axis);
- // currentContactPhysics->shearForce =
q*currentContactPhysics->shearForce;
-
- /// ///
-
- Vector3r x =
currentContactGeometry->contactPoint;
- Vector3r c1x = (x -
de1->se3.position);
- Vector3r c2x = (x -
de2->se3.position);
- /// The following definition of c1x and c2x is to
avoid "granular ratcheting"
- /// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J.
HERRMANN,
- /// "Micro-mechanical investigation of granular
ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- /// ed. T. Triantafyllidis (Balklema, London, 2004),
p. 3-10 - and a lot more papers from the same authors)
- //Vector3r _c1x_ =
currentContactGeometry->radius1*currentContactGeometry->normal;
- //Vector3r _c2x_ =
-currentContactGeometry->radius2*currentContactGeometry->normal;
- //Vector3r relativeVelocity
= (de2->velocity+de2->angularVelocity.Cross(_c2x_)) -
(de1->velocity+de1->angularVelocity.Cross(_c1x_));
- Vector3r relativeVelocity =
(de2->velocity+de2->angularVelocity.Cross(c2x)) -
(de1->velocity+de1->angularVelocity.Cross(c1x));
- Vector3r shearVelocity =
relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
- Vector3r shearDisplacement =
shearVelocity*dt;
- shearForce -=
currentContactPhysics->ks*shearDisplacement;
-
- #endif
-
- // PFC3d SlipModel, is using friction angle. CoulombCriterion
+ // the same as under #else, but refactored
+ #ifdef SCG_SHEAR
+ if(useShear){
+
currentContactGeometry->updateShear(de1,de2,dt);
+
shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
+ } else {
+
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+ }
+ #else
+
currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
+ #endif
+
+ // PFC3d SlipModel, is using friction angle.
CoulombCriterion
Real maxFs =
currentContactPhysics->normalForce.SquaredLength() *
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.SquaredLength() > maxFs )
{
- maxFs = Mathr::Sqrt(maxFs) /
shearForce.Length();
- shearForce *= maxFs;
+ Real ratio = Mathr::Sqrt(maxFs) /
shearForce.Length();
+ shearForce *= ratio;
+ #ifdef SCG_SHEAR
+ // in this case, total shear
displacement must be updated as well
+ if(useShear)
currentContactGeometry->shear*=ratio;
+ #endif
}
- ////////// PFC3d SlipModel
+ ////////// PFC3d SlipModel
Vector3r f=currentContactPhysics->normalForce +
shearForce;
Vector3r
_c1x(currentContactGeometry->contactPoint-de1->se3.position),
Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-03-19
17:55:37 UTC (rev 1725)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.hpp 2009-03-22
14:01:25 UTC (rev 1726)
@@ -11,6 +11,9 @@
#include<yade/core/InteractionSolver.hpp>
#include<yade/core/PhysicalAction.hpp>
+// only to see whether SCG_SHEAR is defined, may be removed in the future
+#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+
#include <set>
#include <boost/tuple/tuple.hpp>
@@ -45,6 +48,9 @@
public :
int sdecGroupMask;
bool momentRotationLaw;
+ #ifdef SCG_SHEAR
+ bool useShear;
+ #endif
ElasticContactLaw();
void action(MetaBody*);
Modified:
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
===================================================================
---
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
2009-03-19 17:55:37 UTC (rev 1725)
+++
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.cpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -22,6 +22,8 @@
#include <Wm3Plane3.h>
#endif
+CREATE_LOGGER(Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact);
+
bool is_point_inside_box(InteractingBox* b, Vector3r P)
{
return std::abs(P[0]) < std::abs(b->extents[0])
@@ -231,7 +233,7 @@
//return true;
#else
- LOG_FATAL("Using miniWm3; recompile with full Wm3 support to
make snow folly functional.");
+ LOG_FATAL("Using miniWm3; recompile with full Wm3 support to
make snow fully functional.");
throw runtime_error("full wm3 required (message above).");
#endif
}
Modified:
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.hpp
===================================================================
---
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.hpp
2009-03-19 17:55:37 UTC (rev 1725)
+++
trunk/pkg/snow/Engine/Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact.hpp
2009-03-22 14:01:25 UTC (rev 1726)
@@ -30,6 +30,8 @@
const Se3r& se32,
const shared_ptr<Interaction>& c);
+ DECLARE_LOGGER;
+
REGISTER_CLASS_NAME(Ef2_InteractingBox_BssSnowGrain_makeIstSnowLayersContact);
REGISTER_BASE_CLASS_NAME(InteractionGeometryEngineUnit);
Added: trunk/scripts/test/shear.py
===================================================================
--- trunk/scripts/test/shear.py 2009-03-19 17:55:37 UTC (rev 1725)
+++ trunk/scripts/test/shear.py 2009-03-22 14:01:25 UTC (rev 1726)
@@ -0,0 +1,53 @@
+#
+# script for testing SpheresContactGeometry shear computation
+# Runs the same 2-sphere setup for useShear=False and for useShear=True
+#
+
+O.bodies.append([
+ utils.sphere([0,0,0],.5000001,dynamic=False,color=(1,0,0)),
+ utils.sphere([0,0,1],.5000001,dynamic=True,color=(0,0,1))
+])
+O.engines=[
+ MetaEngine('BoundingVolumeMetaEngine',[
+ EngineUnit('InteractingSphere2AABB'),
+ EngineUnit('InteractingFacet2AABB'),
+ ]),
+ StandAloneEngine('PersistentSAPCollider',{'haveDistantTransient':True}),
+ MetaEngine('InteractionGeometryMetaEngine',[
+
EngineUnit('InteractingSphere2InteractingSphere4SpheresContactGeometry'),
+ ]),
+ MetaEngine('InteractionPhysicsMetaEngine',[
+ EngineUnit('SimpleElasticRelationships')
+ ]),
+
DeusExMachina('RotationEngine',{'rotationAxis':[1,1,0],'angularVelocity':.001,'subscribedBodies':[1]}),
+ StandAloneEngine('ElasticContactLaw',{'useShear':False}),
+
StandAloneEngine('PeriodicPythonRunner',{'iterPeriod':10000,'command':'interInfo()'}),
+ #DeusExMachina('NewtonsDampedLaw',{'damping':0})
+]
+O.miscParams=[
+ Generic('GLDrawSphere',{'glutUse':True})
+]
+
+O.dt=1e-8
+O.saveTmp('init')
+
+def interInfo():
+ i=O.interactions[0,1]
+ if 1:
+ print O.time,i.phys['shearForce']
+ else:
+ r1,r2=O.bodies[0].shape['radius'],O.bodies[1].shape['radius']
+ theta=[e['angularVelocity'] for e in O.engines if
e.name=='RotationEngine'][0]
+ f=.5*(r1+r2)*theta*O.time*i.phys['ks']
+ print O.time,i.phys['shearForce'],f,i.phys['shearForce'][0]/f
+
+
+print '=========== no shear ============'
+O.run(100000,True)
+nIter=O.iter
+print '============= shear ============='
+O.loadTmp('init')
+ecl=[e for e in O.engines if e.name=='ElasticContactLaw'][0]
+ecl['useShear']=True
+O.run(nIter,True)
+quit()
Added: trunk/scripts/test/triax-identical-results.py
===================================================================
--- trunk/scripts/test/triax-identical-results.py 2009-03-19 17:55:37 UTC
(rev 1725)
+++ trunk/scripts/test/triax-identical-results.py 2009-03-22 14:01:25 UTC
(rev 1726)
@@ -0,0 +1,32 @@
+# Test if different algorithm version give same results for TriaxialTest.
+#
+# The first run creates initial sphere packing, and saves resulting positions
+# after 2000 steps. Subsequent runs only save resulting positions.
+#
+# Compare output files with diff to see if the are 100% identical
+#
+from os.path import exists
+sph='triax-identical-results'
+i=0; outSph=''
+while True:
+ outSph='%s-out%02d.spheres'%(sph,i)
+ if not exists(outSph): break
+ i+=1
+inSph='%s-in.spheres'%sph
+if exists(inSph):
+ print "Using existing initial configuration",inSph
+Preprocessor('TriaxialTest',{'importFilename':(inSph if exists(inSph) else
''),'numberOfGrains':400}).load()
+if not exists(inSph):
+ print "Saving new initial configuration to",inSph
+ utils.spheresToFile(inSph)
+O.usesTimeStepper=False
+O.dt=utils.PWaveTimeStep()
+#
+# uncomment this line to enable shear computation in SpheresContactGeometry
and then compare results with this line commented
+#
+#[e for e in O.engines if e.name=='ElasticContactLaw'][0]['useShear']=True
+if 1:
+ O.run(2000,True)
+ utils.spheresToFile(outSph)
+ print "Results saved to",outSph
+ quit()
_______________________________________________
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
