yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #01163
[svn] r1726 - in trunk: . core gui/py lib/factory pkg/common/DataClass/InteractionPhysics pkg/dem/DataClass/InteractionGeometry pkg/dem/Engine/EngineUnit pkg/dem/Engine/StandAloneEngine pkg/snow/Engine scripts/test
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()
