yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-05-24 20:22:30 +0200 (Sun, 24 May 2009)
New Revision: 1778

Added:
   trunk/extra/Shop.cpp
   trunk/extra/Shop.hpp
   trunk/pkg/dem/ConcretePM.cpp
   trunk/pkg/dem/ConcretePM.hpp
   trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp
   trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.hpp
   trunk/pkg/dem/PreProcessor/SimpleScene.cpp
   trunk/pkg/dem/PreProcessor/SimpleScene.hpp
Removed:
   trunk/core/yadeExceptions.cpp
   trunk/core/yadeExceptions.hpp
   trunk/extra/Brefcom.cpp
   trunk/extra/Brefcom.hpp
   trunk/extra/BrefcomTestGen.cpp
   trunk/extra/BrefcomTestGen.hpp
   trunk/extra/SimpleScene.cpp
   trunk/extra/SimpleScene.hpp
   trunk/extra/clump/
   trunk/extra/usct/
Modified:
   trunk/SConstruct
   trunk/core/InteractionContainer.cpp
   trunk/core/Omega.cpp
   trunk/core/Omega.hpp
   trunk/core/Preferences.cpp
   trunk/core/Preferences.hpp
   trunk/core/SConscript
   trunk/core/yade.cpp
   trunk/extra/SConscript
   trunk/gui/SConscript
   trunk/gui/py/PythonUI_rc.py
   trunk/gui/py/_eudoxos.cpp
   trunk/gui/qt3/SimulationController.cpp
   trunk/lib/factory/ClassFactory.cpp
   trunk/lib/factory/ClassFactory.hpp
   trunk/lib/factory/DynLibManager.cpp
   trunk/lib/factory/DynLibManager.hpp
   trunk/pkg/dem/PreProcessor/HydraulicTest.cpp
   trunk/pkg/dem/PreProcessor/ThreePointBending.cpp
   trunk/pkg/dem/SConscript
Log:
1. Add logic to scons to delete files that will not be installed in this run but are on-disk. This should make it unnecessary to remove files by hand. In the same way, disable implicit target cache which was making problems if files shuffle around.
2. Remove Brefcom* files, move to pkg/dem/ConcretePM.?pp. (not sure what other place, since it contains Engines, DataClasses etc in one file that I would refuse to break.).
3. Add detailed comments on ConcretePM to the ConcretePM.hpp
4. Remove extra/usct, moveUniaxialStrainer and put to pkg/dem; move SimpleScene from extra to pkg/dem
5. Remove core/yadeExceptions.*
6. Remove Preferences::dynLibDirectores and Preferences::version. Omega now recursively searches the lib directory for plugins and loads them. preferences.xml has only defaultGUILibName now. Finally. Remove baseDirs and other cruft from DynLibManager.



Modified: trunk/SConstruct
===================================================================
--- trunk/SConstruct	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/SConstruct	2009-05-24 18:22:30 UTC (rev 1778)
@@ -346,7 +346,7 @@
 SCons.Defaults.DefaultEnvironment(tools = [])
 env.Decider('MD5-timestamp')
 env.SetOption('max_drift',5) # cache md5sums of files older than 5 seconds
-SetOption('implicit_cache',1) # cache #include files etc.
+SetOption('implicit_cache',0) # cache #include files etc.
 env.SourceCode(".",None) # skip dotted directories
 SetOption('num_jobs',env['jobs'])
 
@@ -561,4 +561,14 @@
 # read top-level SConscript file. It is used only so that build_dir is set. This file reads all SConscripts from in yadeModules
 env.SConscript(dirs=['.'],build_dir=buildDir,duplicate=0)
 
+#################################################################################
+## remove plugins that are in the target dir but will not be installed now
+toInstall=[str(node) for node in env.FindInstalledFiles()]
+for root,dirs,files in os.walk(env.subst('$PREFIX/lib/yade${SUFFIX}')):
+	for f in files:
+		ff=os.path.join(root,f)
+		if ff not in toInstall and not ff.endswith('.pyo'):
+			print "Deleting extra plugin", ff
+			os.remove(ff)
+
 #Progress('.', interval=100, file=sys.stderr)

Modified: trunk/core/InteractionContainer.cpp
===================================================================
--- trunk/core/InteractionContainer.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/InteractionContainer.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -23,7 +23,7 @@
 	}
 	else
 	{
-		interaction.clear();
+		interaction.clear(); interaction.reserve(this->size());
 		InteractionContainer::iterator i    = this->begin();
 		InteractionContainer::iterator iEnd = this->end();
 		for( ; i!=iEnd ; ++i )

Modified: trunk/core/Omega.cpp
===================================================================
--- trunk/core/Omega.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/Omega.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -9,11 +9,9 @@
 *************************************************************************/
 
 #include"Omega.hpp"
-#include"yadeExceptions.hpp"
 #include"MetaBody.hpp"
 #include"TimeStepper.hpp"
 #include"ThreadRunner.hpp"
-#include"Preferences.hpp"
 #include<Wm3Vector3.h>
 #include<yade/lib-base/yadeWm3.hpp>
 #include<yade/lib-serialization/IOFormatManager.hpp>
@@ -23,6 +21,7 @@
 #include<cstdlib>
 #include<boost/filesystem/operations.hpp>
 #include<boost/filesystem/convenience.hpp>
+#include<boost/filesystem/exception.hpp>
 #include<boost/algorithm/string.hpp>
 #include<boost/thread/mutex.hpp>
 #include<boost/version.hpp>
@@ -164,36 +163,33 @@
 	return (dynlibs[className].baseClasses.find(baseClassName)!=dynlibs[className].baseClasses.end());
 }
 
-void Omega::scanPlugins(){
-	FOREACH(string dld,preferences->dynlibDirectories) ClassFactory::instance().addBaseDirectory(dld);
-	vector<string> dynlibsList;
-	FOREACH(string si, preferences->dynlibDirectories){
-		filesystem::path directory(si);
-		if(!filesystem::exists(directory)){LOG_ERROR("Nonexistent plugin directory: "<<directory.native_directory_string()<<".");continue; }
-		filesystem::directory_iterator di(directory),diEnd;
-		FOREACH(filesystem::path pth,std::make_pair(di,diEnd)){
-			// node is not a directory and is either regular file or non-dangling symlink; and extension is not ".a"; AND moreover transforming it to library name and back to filename is identity; otherwise the file wouldn't be loaded by the DynLibManager anyway
-			if (!filesystem::is_directory(*di) && filesystem::exists(*di) && filesystem::extension(*di)!=".a" &&
-				ClassFactory::instance().libNameToSystemName(ClassFactory::instance().systemNameToLibName(filesystem::basename(pth)))==(pth.leaf())){
-				filesystem::path name(filesystem::basename(pth));
-				if(name.leaf().length()<1) continue; // filter out 0-length filenames
-				string plugin=name.leaf();
-				if(!ClassFactory::instance().load(ClassFactory::instance().systemNameToLibName(plugin))){
-					string err=ClassFactory::instance().lastError();
-					if(err.find(": undefined symbol: ")!=std::string::npos){
-						size_t pos=err.rfind(":");	assert(pos!=std::string::npos);
-						std::string sym(err,pos+2); //2 removes ": " from the beginning
-						int status=0; char* demangled_sym=abi::__cxa_demangle(sym.c_str(),0,0,&status);
-						LOG_FATAL(plugin<<": undefined symbol `"<<demangled_sym<<"'"); LOG_FATAL(plugin<<": "<<err); LOG_FATAL("Bailing out.");
-					}
-					else {
-						LOG_FATAL(plugin<<": "<<err<<" ."); /* leave space to not to confuse c++filt */ LOG_FATAL("Bailing out.");
-					}
-					abort();
+void Omega::scanPlugins(string baseDir){
+	try{
+		filesystem::recursive_directory_iterator Iend;
+		for(filesystem::recursive_directory_iterator I(baseDir); I!=Iend; ++I){ 
+			filesystem::path pth=I->path();
+			if(filesystem::is_directory(pth) || ClassFactory::instance().libNameToSystemName(ClassFactory::instance().systemNameToLibName(filesystem::basename(pth)))!=(pth.leaf())){ LOG_DEBUG("File not considered a plugin: "<<pth.leaf()<<"."); continue; }
+			LOG_DEBUG("Trying "<<pth.leaf());
+			filesystem::path name(filesystem::basename(pth));
+			if(name.leaf().length()<1) continue; // filter out 0-length filenames
+			string plugin=name.leaf();
+			if(!ClassFactory::instance().load(ClassFactory::instance().systemNameToLibName(plugin))){
+				string err=ClassFactory::instance().lastError();
+				if(err.find(": undefined symbol: ")!=std::string::npos){
+					size_t pos=err.rfind(":");	assert(pos!=std::string::npos);
+					std::string sym(err,pos+2); //2 removes ": " from the beginning
+					int status=0; char* demangled_sym=abi::__cxa_demangle(sym.c_str(),0,0,&status);
+					LOG_FATAL(plugin<<": undefined symbol `"<<demangled_sym<<"'"); LOG_FATAL(plugin<<": "<<err); LOG_FATAL("Bailing out.");
 				}
+				else {
+					LOG_FATAL(plugin<<": "<<err<<" ."); /* leave space to not to confuse c++filt */ LOG_FATAL("Bailing out.");
+				}
+				abort();
 			}
-			else LOG_DEBUG("File not considered a plugin: "<<pth.leaf()<<".");
 		}
+	} catch(filesystem::basic_filesystem_error<filesystem::path>& e) {
+		LOG_FATAL("Error from recursive plugin directory scan (unreadable directory?): "<<e.what());
+		throw;
 	}
 	list<string>& plugins(ClassFactory::instance().pluginClasses);
 	plugins.sort(); plugins.unique();
@@ -205,6 +201,7 @@
 	resetRootBody();
 	IOFormatManager::loadFromStream("XMLFormatManager",stream,"rootBody",rootBody);
 }
+
 void Omega::saveSimulationToStream(std::ostream& stream){
 	LOG_DEBUG("Saving simulation to stream.");
 	IOFormatManager::saveToStream("XMLFormatManager",stream,"rootBody",rootBody);
@@ -212,8 +209,8 @@
 
 void Omega::loadSimulation(){
 
-	if(Omega::instance().getSimulationFileName().size()==0) throw yadeBadFile("Simulation filename to load has zero length");
-	if(!filesystem::exists(simulationFileName) && !algorithm::starts_with(simulationFileName,":memory")) throw yadeBadFile((std::string("Simulation file to load doesn't exist")+simulationFileName).c_str());
+	if(Omega::instance().getSimulationFileName().size()==0) throw runtime_error("Empty simulation filename to load.");
+	if(!filesystem::exists(simulationFileName) && !algorithm::starts_with(simulationFileName,":memory")) throw runtime_error("Simulation file to load doesn't exist: "+simulationFileName);
 	
 	LOG_INFO("Loading file " + simulationFileName);
 
@@ -229,16 +226,16 @@
 			RenderMutexLock lock; IOFormatManager::loadFromFile("BINFormatManager",simulationFileName,"rootBody",rootBody);
 		}
 		else if(algorithm::starts_with(simulationFileName,":memory:")){
-			if(memSavedSimulations.count(simulationFileName)==0) throw yadeBadFile(("Cannot load nonexistent memory-saved simulation "+simulationFileName).c_str());
+			if(memSavedSimulations.count(simulationFileName)==0) throw runtime_error("Cannot load nonexistent memory-saved simulation "+simulationFileName);
 			istringstream iss(memSavedSimulations[simulationFileName]);
 			joinSimulationLoop();
 			resetRootBody();
 			RenderMutexLock lock; IOFormatManager::loadFromStream("XMLFormatManager",iss,"rootBody",rootBody);
 		}
-		else throw (yadeBadFile("Extension of file not recognized."));
+		else throw runtime_error("Extension of file to load not recognized "+simulationFileName);
 	}
 
-	if( rootBody->getClassName() != "MetaBody") throw yadeBadFile("Wrong file format (rootBody is not a MetaBody!) ??");
+	if( rootBody->getClassName() != "MetaBody") throw runtime_error("Wrong file format (rootBody is not a MetaBody!?) in "+simulationFileName);
 
 	timeInit();
 
@@ -249,7 +246,7 @@
 
 void Omega::saveSimulation(const string name)
 {
-	if(name.size()==0) throw yadeBadFile("Filename with zero length.");
+	if(name.size()==0) throw runtime_error("Name of file to save has zero length.");
 	LOG_INFO("Saving file " << name);
 
 	if(algorithm::ends_with(name,".xml") || algorithm::ends_with(name,".xml.bz2")){
@@ -267,7 +264,7 @@
 		memSavedSimulations[name]=oss.str();
 	}
 	else {
-		throw(yadeBadFile(("Filename extension not recognized in `"+name+"'").c_str()));
+		throw runtime_error("Filename extension not recognized in `"+name+"'");
 	}
 }
 

Modified: trunk/core/Omega.hpp
===================================================================
--- trunk/core/Omega.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/Omega.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -121,7 +121,7 @@
 		bool		isRunning();
 
 		const		map<string,DynlibDescriptor>& getDynlibsDescriptor();
-		void		scanPlugins();
+		void		scanPlugins(string baseDir);
 		bool		isInheritingFrom(const string& className, const string& baseClassName );
 
 		void		setTimeStep(const Real);

Modified: trunk/core/Preferences.cpp
===================================================================
--- trunk/core/Preferences.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/Preferences.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -8,9 +8,5 @@
 
 #include "Preferences.hpp"
 
+Preferences::Preferences (){}
 
-Preferences::Preferences ()
-{
-	version = 1;
-}
-

Modified: trunk/core/Preferences.hpp
===================================================================
--- trunk/core/Preferences.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/Preferences.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -15,19 +15,13 @@
 
 using namespace std;
 
-class Preferences : public Serializable
-{
+class Preferences : public Serializable{
 	public :
-		int version;	
-		vector<string>	 dynlibDirectories;
-
 		string		 defaultGUILibName;
 		Preferences ();
-
-	REGISTER_ATTRIBUTES(,(version)(dynlibDirectories)(defaultGUILibName));
+	REGISTER_ATTRIBUTES(/* no base class*/,(defaultGUILibName));
 	REGISTER_CLASS_AND_BASE(Preferences,Serializable);
 };
-
 REGISTER_SERIALIZABLE(Preferences);
 
 

Modified: trunk/core/SConscript
===================================================================
--- trunk/core/SConscript	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/SConscript	2009-05-24 18:22:30 UTC (rev 1778)
@@ -25,7 +25,6 @@
 			'ThreadWorker.cpp',
 			'TimeStepper.cpp',
 			'yade.cpp',
-			'yadeExceptions.cpp',
 			'containers/BodyRedirectionVector.cpp',
 			'containers/BodyAssocVector.cpp',
 			'containers/InteractionHashMap.cpp',

Modified: trunk/core/yade.cpp
===================================================================
--- trunk/core/yade.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/yade.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -99,19 +99,11 @@
 
 void firstRunSetup(shared_ptr<Preferences>& pref)
 {
-	const char* libDirs[]={"extra","gui","lib","pkg-common","pkg-dem","pkg-fem","pkg-snow","pkg-lattice","pkg-mass-spring","pkg-realtime-rigidbody",NULL /* sentinel */};
 	string cfgFile=Omega::instance().yadeConfigPath+"/preferences.xml";
-	LOG_INFO("Creating default configuration file: "<<cfgFile<<". Please tune by hand if needed.");
-	string expLibDir;
-	for(int i=0; libDirs[i]!=NULL; i++) {
-		expLibDir=string(PREFIX "/lib/yade" SUFFIX "/") + libDirs[i];
-		if(!filesystem::exists(expLibDir)) continue; // skip modules that were not built/installed
-		LOG_INFO("Adding plugin directory "<<expLibDir<<".");
-		pref->dynlibDirectories.push_back(expLibDir);
-	}
+	LOG_INFO("Creating default configuration file: "<<cfgFile<<". Tune by hand if needed.");
 	LOG_INFO("Setting GUI: QtGUI.");
 	pref->defaultGUILibName="QtGUI";
-	IOFormatManager::saveToFile("XMLFormatManager",Omega::instance().yadeConfigPath+"/preferences.xml","preferences",pref);
+	IOFormatManager::saveToFile("XMLFormatManager",cfgFile,"preferences",pref);
 }
 
 string findRecoveryCandidate(filesystem::path dir, string start){
@@ -257,7 +249,7 @@
 	
 	LOG_INFO("Loading "<<yadeConfigFile.string()); IOFormatManager::loadFromFile("XMLFormatManager",yadeConfigFile.string(),"preferences",Omega::instance().preferences);
 
-	LOG_INFO("Loading plugins"); Omega::instance().scanPlugins();
+	LOG_INFO("Loading plugins"); Omega::instance().scanPlugins(string(PREFIX "/lib/yade" SUFFIX));
 	Omega::instance().init();
 
 	Omega::instance().setSimulationFileName(simulationFileName); //init() resets to "";

Deleted: trunk/core/yadeExceptions.cpp
===================================================================
--- trunk/core/yadeExceptions.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/yadeExceptions.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,12 +0,0 @@
-/*************************************************************************
-*  Copyright (C) 2004 by Janek Kozicki                                   *
-*  cosurgi@xxxxxxxxxx                                                    *
-*                                                                        *
-*  This program is free software; it is licensed under the terms of the  *
-*  GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#include "yadeExceptions.hpp"
-
-const char* yadeExceptions::BadFile	= "Please specify correct filename using your frontend. Following filename is not correct: ";
-

Deleted: trunk/core/yadeExceptions.hpp
===================================================================
--- trunk/core/yadeExceptions.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/core/yadeExceptions.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,31 +0,0 @@
-/*************************************************************************
-*  Copyright (C) 2004 by Janek Kozicki                                   *
-*  cosurgi@xxxxxxxxxx                                                    *
-*                                                                        *
-*  This program is free software; it is licensed under the terms of the  *
-*  GNU General Public License v2 or later. See file LICENSE for details. *
-*************************************************************************/
-
-#pragma once
-
-#include <string>
-#include <stdexcept>
-
-using namespace std;
-
-struct yadeError : public std::runtime_error
-{
-	explicit yadeError(const char* msg) : std::runtime_error(msg) {};
-};
-
-struct yadeBadFile : public yadeError
-{
-	explicit yadeBadFile(const char* msg) : yadeError(msg) {};
-};
-
-struct yadeExceptions
-{
-	static const char* BadFile;
-};
-
-

Deleted: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/Brefcom.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,362 +0,0 @@
-// 2007,2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
-#include"Brefcom.hpp"
-#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-dem/BodyMacroParameters.hpp>
-#include<yade/pkg-common/Sphere.hpp>
-#include<yade/lib-QGLViewer/qglviewer.h>
-#include<yade/lib-opengl/GLUtils.hpp>
-#include<yade/pkg-dem/DemXDofGeom.hpp>
-
-YADE_PLUGIN("Ip2_CpmMat_CpmMat_CpmPhys","CpmPhys","GLDrawCpmPhys","CpmPhysDamageColorizer", "CpmMat", "CpmGlobalCharacteristics", "Law2_Dem3DofGeom_CpmPhys_Cpm");
-
-CREATE_LOGGER(CpmGlobalCharacteristics);
-
-void CpmGlobalCharacteristics::compute(MetaBody* rb, bool useMaxForce){
-	rb->bex.sync();
-
-	// 1. reset volumetric strain (cummulative in the next loop)
-	// 2. get maximum force on a body and sum of all forces (for averaging)
-	Real sumF=0,maxF=0,currF;
-	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
-	CpmMat* bpp(YADE_CAST<CpmMat*>(b->physicalParameters.get()));
-		bpp->epsVolumetric=0;
-		bpp->numContacts=0;
-		currF=rb->bex.getForce(b->id).Length(); maxF=max(currF,maxF); sumF+=currF;
-	}
-	Real meanF=sumF/rb->bodies->size(); 
-
-	// commulate normal strains from contacts
-	// get max force on contacts
-	Real maxContactF=0;
-	FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
-		if(!I->isReal) continue;
-		shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
-		maxContactF=max(maxContactF,max(BC->Fn,BC->Fs.Length()));
-		CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
-		CpmMat* bpp2(YADE_CAST<CpmMat*>(Body::byId(I->getId2())->physicalParameters.get()));
-		bpp1->epsVolumetric+=BC->epsN; bpp1->numContacts+=1;
-		bpp2->epsVolumetric+=BC->epsN; bpp2->numContacts+=1;
-	}
-	unbalancedForce=(useMaxForce?maxF:meanF)/maxContactF;
-
-	FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
-		if(!I->isReal) continue;
-		shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
-		CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
-		CpmMat* bpp2(YADE_CAST<CpmMat*>(Body::byId(I->getId2())->physicalParameters.get()));
-		Real epsVolAvg=.5*((3./bpp1->numContacts)*bpp1->epsVolumetric+(3./bpp2->numContacts)*bpp2->epsVolumetric);
-		BC->epsTrans=(epsVolAvg-BC->epsN)/2.;
-		//TRVAR5(I->getId1(),I->getId2(),BC->epsTrans,(3./bpp1->numContacts)*bpp1->epsVolumetric,(3./bpp2->numContacts)*bpp2->epsVolumetric);
-		//TRVAR4(bpp1->numContacts,bpp1->epsVolumetric,bpp2->numContacts,bpp2->epsVolumetric);
-	}
-	#if 0
-		FOREACH(const shared_ptr<Body>& b, *rb->bodies){
-			CpmMat* bpp(YADE_PTR_CAST<CpmMat>(b->physicalParameters.get()));
-			bpp->epsVolumeric*=3/bpp->numContacts;
-		}
-	#endif
-
-
-}
-
-
-/********************** Ip2_CpmMat_CpmMat_CpmPhys ****************************/
-CREATE_LOGGER(Ip2_CpmMat_CpmMat_CpmPhys);
-
-
-void Ip2_CpmMat_CpmMat_CpmPhys::go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction){
-	Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
-
-	assert(contGeom); // for now, don't handle anything other than SpheresContactGeometry and Dem3DofGeom
-
-	if(!interaction->isNew && interaction->interactionPhysics){ /* relax */ } 
-	else {
-		interaction->isNew=false; // just in case
-
-		const shared_ptr<BodyMacroParameters>& elast1=static_pointer_cast<BodyMacroParameters>(pp1);
-		const shared_ptr<BodyMacroParameters>& elast2=static_pointer_cast<BodyMacroParameters>(pp2);
-
-		Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young); // harmonic Young's modulus average
-		//Real nu12=2*elast1->poisson*elast2->poisson/(elast1->poisson+elast2->poisson); // dtto for Poisson ratio 
-		Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
-		Real S12=Mathr::PI*pow(minRad,2); // "surface" of interaction
-		//Real E=(E12 /* was here for Kn:  *S12/d0  */)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
-		//Real E=E12; // apply alpha, beta, gamma: garbage values of E !?
-
-		if(!neverDamage) { assert(!isnan(sigmaT)); }
-
-		shared_ptr<CpmPhys> contPhys(new CpmPhys());
-
-		contPhys->E=E12;
-		contPhys->G=E12*G_over_E;
-		contPhys->tanFrictionAngle=tan(.5*(elast1->frictionAngle+elast2->frictionAngle));
-		contPhys->undamagedCohesion=sigmaT;
-		contPhys->crossSection=S12;
-		contPhys->epsCrackOnset=epsCrackOnset;
-		contPhys->epsFracture=relDuctility*epsCrackOnset;
-		contPhys->omegaThreshold=omegaThreshold;
-		// inherited from NormalShearInteracion, used in the timestepper
-		contPhys->kn=contPhys->E*contPhys->crossSection;
-		contPhys->ks=contPhys->G*contPhys->crossSection;
-
-		if(neverDamage) contPhys->neverDamage=true;
-		if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
-		else contPhys->isCohesive=false;
-		contPhys->dmgTau=dmgTau;
-		contPhys->dmgRateExp=dmgRateExp;
-		contPhys->plTau=plTau;
-		contPhys->plRateExp=plRateExp;
-		contPhys->isoPrestress=isoPrestress;
-
-		interaction->interactionPhysics=contPhys;
-	}
-}
-
-
-
-
-/********************** CpmPhys ****************************/
-CREATE_LOGGER(CpmPhys);
-
-// !! at least one virtual function in the .cpp file
-CpmPhys::~CpmPhys(){};
-
-CREATE_LOGGER(Law2_Dem3DofGeom_CpmPhys_Cpm);
-
-long CpmPhys::cummBetaIter=0, CpmPhys::cummBetaCount=0;
-
-Real CpmPhys::solveBeta(const Real c, const Real N){
-	#ifdef YADE_DEBUG
-		cummBetaCount++;
-	#endif
-	const int maxIter=20;
-	const Real maxError=1e-12;
-	Real f, ret=0.;
-	for(int i=0; i<maxIter; i++){
-		#ifdef YADE_DEBUG
-			cummBetaIter++;
-		#endif
-		Real aux=c*exp(N*ret)+exp(ret);
-		f=log(aux);
-		if(fabs(f)<maxError) return ret;
-		Real df=(c*N*exp(N*ret)+exp(ret))/aux;
-		ret-=f/df;
-	}
-	LOG_FATAL("No convergence after "<<maxIter<<" iters; c="<<c<<", N="<<N<<", ret="<<ret<<", f="<<f);
-	throw runtime_error("Law2_Dem3DofGeom_CpmPhys_Cpm::solveBeta failed to converge.");
-}
-
-Real CpmPhys::computeDmgOverstress(Real dt){
-	if(dmgStrain>=epsN*omega){ // unloading, no viscous stress
-		dmgStrain=epsN*omega;
-		LOG_TRACE("Elastic/unloading, no viscous overstress");
-		return 0.;
-	}
-	Real c=epsCrackOnset*(1-omega)*pow(dmgTau/dt,dmgRateExp)*pow(epsN*omega-dmgStrain,dmgRateExp-1.);
-	Real beta=solveBeta(c,dmgRateExp);
-	Real deltaDmgStrain=(epsN*omega-dmgStrain)*exp(beta);
-	dmgStrain+=deltaDmgStrain;
-	LOG_TRACE("deltaDmgStrain="<<deltaDmgStrain<<", viscous overstress "<<(epsN*omega-dmgStrain)*E);
-	/* σN=Kn(εN-εd); dmgOverstress=σN-(1-ω)*Kn*εN=…=Kn(ω*εN-εd) */
-	return (epsN*omega-dmgStrain)*E;
-}
-
-Real CpmPhys::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt){
-	if(sigmaTNorm<sigmaTYield) return 1.;
-	Real c=undamagedCohesion*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
-	Real beta=solveBeta(c,plRateExp);
-	//LOG_DEBUG("scaling factor "<<1.-exp(beta)*(1-sigmaTYield/sigmaTNorm));
-	return 1.-exp(beta)*(1-sigmaTYield/sigmaTNorm);
-}
-
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2=1.; /* deactivated if > 0 */
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed=1.;
-Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift=0.;
-
-void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
-	//timingDeltas->start();
-	Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
-	CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
-
-	/* kept fully damaged contacts; note that normally the contact is deleted _after_ the CPM_MATERIAL_MODEL,
-	 * i.e. if it is 1.0 here, omegaThreshold is >= 1.0 for sure.
-	 * &&'ing that just to make sure anyway ...
-	 */
-	// if(BC->omega>=1.0 && BC->omegaThreshold>=1.0) return;
-
-	// shorthands
-	Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(BC->omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep();  const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
-	/* const Real& transStrainCoeff(BC->transStrainCoeff); const Real& epsTrans(BC->epsTrans); const Real& xiShear(BC->xiShear); */
-	Real& omega(BC->omega); Real& sigmaN(BC->sigmaN);  Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
-	const Real& yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed); const int& yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType);
-	const Real& yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift); 
-
-	#define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); throw;}
-	
-	#define NNAN(a) YADE_VERIFY(!isnan(a));
-	#define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
-
-	//timingDeltas->checkpoint("setup");
-	// if(contGeom->refR1<0) contGeom->refLength=contGeom->refR2; // make facet-sphere contact always at equilibrium when touching exactly (and not the initial distance)
-	epsN=contGeom->strainN(); epsT=contGeom->strainT();
-	#ifdef YADE_DEBUG
-		if(isnan(epsN)){
-			LOG_FATAL("refLength="<<contGeom->refLength<<"; pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<"; displacementN="<<contGeom->displacementN());
-			throw runtime_error("!! epsN==NaN !!");
-		}
-		NNAN(epsN); NNANV(epsT);
-	#endif
-	// already in SpheresContactGeometry:
-	// contGeom->relocateContactPoints(); // allow very large mutual rotations
-	if(logStrain && epsN<0){
-		Real epsN0=epsN;
-		epsN=log(epsN0+1); epsT*=epsN/epsN0;
-	}
-	NNAN(epsN); NNANV(epsT);
-	//timingDeltas->checkpoint("geom");
-
-	epsN+=BC->isoPrestress/E;
-	//TRVAR1(epsN);
-	#ifdef CPM_MATERIAL_MODEL
-		CPM_MATERIAL_MODEL
-	#else
-		sigmaN=E*epsN;
-		sigmaT=G*epsT;
-	#endif
-	sigmaN-=BC->isoPrestress;
-	if(contGeom->refR1<0 && Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2<=0 && epsN<Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2){
-		/* move Body2 (the sphere) so that minStrain is satisfied */
-		rootBody->bex.addMove(I->getId2(),contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
-		LOG_TRACE("Moving by "<<contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
-	}
-	NNAN(kappaD); NNAN(epsCrackOnset); NNAN(epsFracture); NNAN(omega);
-	NNAN(sigmaN); NNANV(sigmaT); NNAN(crossSection);
-	//timingDeltas->checkpoint("material");
-
-	//const int watch1=6300, watch2=6299;
-	//#define SHOW(a) if((I->getId1()==watch1 && I->getId2()==watch2) || (I->getId2()==watch1 && I->getId1()==watch2)) cerr<<__FILE__<<":"<<__LINE__<<" "<<a<<endl;
-	//SHOW("epsN"<<epsN);
-	if(epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)){
-		rootBody->interactions->requestErase(I->getId1(),I->getId2());
-		if(isCohesive){
-			const shared_ptr<Body>& body1=Body::byId(I->getId1(),rootBody), body2=Body::byId(I->getId2(),rootBody); assert(body1); assert(body2);
-			const shared_ptr<CpmMat>& rbp1=YADE_PTR_CAST<CpmMat>(body1->physicalParameters), rbp2=YADE_PTR_CAST<CpmMat>(body2->physicalParameters);
-			if(BC->isCohesive){rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum;}
-			LOG_DEBUG("Contact #"<<I->getId1()<<"=#"<<I->getId2()<<" is damaged over thershold ("<<omega<<">"<<omegaThreshold<<") and will be deleted.");
-		}
-		return;
-	}
-
-	Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
-	Fs=sigmaT*crossSection; BC->shearForce=Fs;
-
-	applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, rootBody);
-	//timingDeltas->checkpoint("rest");
-}
-
-
-/********************** GLDrawCpmPhys ****************************/
-
-#include<yade/lib-opengl/OpenGLWrapper.hpp>
-
-CREATE_LOGGER(GLDrawCpmPhys);
-
-bool GLDrawCpmPhys::contactLine=true;
-bool GLDrawCpmPhys::dmgLabel=true;
-bool GLDrawCpmPhys::dmgPlane=false;
-bool GLDrawCpmPhys::epsNLabel=true;
-bool GLDrawCpmPhys::epsT=false;
-bool GLDrawCpmPhys::epsTAxes=false;
-bool GLDrawCpmPhys::normal=false;
-bool GLDrawCpmPhys::colorStrain=false;
-
-
-void GLDrawCpmPhys::go(const shared_ptr<InteractionPhysics>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
-	const shared_ptr<CpmPhys>& BC=static_pointer_cast<CpmPhys>(ip);
-	const shared_ptr<Dem3DofGeom>& geom=YADE_PTR_CAST<Dem3DofGeom>(i->interactionGeometry);
-
-	//Vector3r lineColor(BC->omega,1-BC->omega,0.0); /* damaged links red, undamaged green */
-	Vector3r lineColor=Shop::scalarOnColorScale(1.-BC->relResidualStrength);
-
-	if(colorStrain) lineColor=Vector3r(
-		min((Real)1.,max((Real)0.,-BC->epsTrans/BC->epsCrackOnset)),
-		min((Real)1.,max((Real)0.,BC->epsTrans/BC->epsCrackOnset)),
-		min((Real)1.,max((Real)0.,abs(BC->epsTrans)/BC->epsCrackOnset-1)));
-
-	if(contactLine) GLUtils::GLDrawLine(b1->physicalParameters->dispSe3.position,b2->physicalParameters->dispSe3.position,lineColor);
-	if(dmgLabel){ GLUtils::GLDrawNum(BC->omega,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
-	else if(epsNLabel){ GLUtils::GLDrawNum(BC->epsN,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
-	if(BC->omega>0 && dmgPlane){
-		Real halfSize=sqrt(1-BC->relResidualStrength)*.5*.705*sqrt(BC->crossSection);
-		Vector3r midPt=.5*Vector3r(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position);
-		glDisable(GL_CULL_FACE);
-		glPushMatrix();
-			glTranslatev(midPt);
-			Quaternionr q; q.Align(Vector3r::UNIT_Z,geom->normal);
-			Vector3r axis; Real angle; q.ToAxisAngle(axis,angle);
-			glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
-			glBegin(GL_POLYGON);
-				glColor3v(lineColor); 
-				glVertex3d(halfSize,0.,0.);
-				glVertex3d(.5*halfSize,.866*halfSize,0.);
-				glVertex3d(-.5*halfSize,.866*halfSize,0.);
-				glVertex3d(-halfSize,0.,0.);
-				glVertex3d(-.5*halfSize,-.866*halfSize,0.);
-				glVertex3d(.5*halfSize,-.866*halfSize,0.);
-			glEnd();
-		glPopMatrix();
-	}
-
-	const Vector3r& cp=static_pointer_cast<Dem3DofGeom>(i->interactionGeometry)->contactPoint;
-	if(epsT){
-		Real maxShear=(BC->undamagedCohesion-BC->sigmaN*BC->tanFrictionAngle)/BC->G;
-		Real relShear=BC->epsT.Length()/maxShear;
-		Real scale=.5*geom->refLength;
-		Vector3r dirShear=BC->epsT; dirShear.Normalize();
-		if(epsTAxes){
-			GLUtils::GLDrawLine(cp-Vector3r(scale,0,0),cp+Vector3r(scale,0,0));
-			GLUtils::GLDrawLine(cp-Vector3r(0,scale,0),cp+Vector3r(0,scale,0));
-			GLUtils::GLDrawLine(cp-Vector3r(0,0,scale),cp+Vector3r(0,0,scale));
-		}
-		GLUtils::GLDrawArrow(cp,cp+dirShear*relShear*scale,Vector3r(1.,0.,0.));
-		GLUtils::GLDrawLine(cp+dirShear*relShear*scale,cp+dirShear*scale,Vector3r(.3,.3,.3));
-
-		/* normal strain */ GLUtils::GLDrawArrow(cp,cp+geom->normal*(BC->epsN/maxShear),Vector3r(0.,1.,0.));
-	}
-	//if(normal) GLUtils::GLDrawArrow(cp,cp+geom->normal*.5*BC->equilibriumDist,Vector3r(0.,1.,0.));
-}
-
-struct BodyStats{ short nCohLinks; Real dmgSum; Real epsPlSum; BodyStats(): nCohLinks(0), dmgSum(0), epsPlSum(0.){} };
-
-/********************** CpmPhysDamageColorizer ****************************/
-void CpmPhysDamageColorizer::action(MetaBody* rootBody){
-	//vector<pair<short,Real> > bodyDamage; /* number of cohesive interactions per body; cummulative damage of interactions */
-	//vector<pair<short,
-	vector<BodyStats> bodyStats; bodyStats.resize(rootBody->bodies->size());
-	assert(bodyStats[0].nCohLinks==0); // should be initialized by dfault ctor
-	FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
-		shared_ptr<CpmPhys> BC=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
-		if(!BC || !BC->isCohesive) continue;
-		const body_id_t id1=I->getId1(), id2=I->getId2();
-		bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-BC->relResidualStrength); bodyStats[id1].epsPlSum+=BC->epsPlSum;
-		bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-BC->relResidualStrength); bodyStats[id2].epsPlSum+=BC->epsPlSum;
-		//bodyDamage[id1].first++; bodyDamage[id2].first++;
-		//bodyDamage[id1].second+=(1-BC->relResidualStrength); bodyDamage[id2].second+=(1-BC->relResidualStrength);
-		maxOmega=max(maxOmega,BC->omega);
-	}
-	FOREACH(shared_ptr<Body> B, *rootBody->bodies){
-		body_id_t id=B->getId();
-		// add damaged contacts that have already been deleted
-		CpmMat* bpp=dynamic_cast<CpmMat*>(B->physicalParameters.get());
-		if(!bpp) continue;
-		short cohLinksWhenever=bodyStats[id].nCohLinks+bpp->numBrokenCohesive;
-		if(cohLinksWhenever>0){
-			bpp->normDmg=(bodyStats[id].dmgSum+bpp->numBrokenCohesive)/cohLinksWhenever;
-			bpp->normEpsPl=(bodyStats[id].epsPlSum+bpp->epsPlBroken)/cohLinksWhenever;
-		}
-		else { bpp->normDmg=0; bpp->normEpsPl=0;}
-		B->geometricalModel->diffuseColor=Vector3r(bpp->normDmg,1-bpp->normDmg,B->isDynamic?0:1);
-	}
-}
-
-

Deleted: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/Brefcom.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,272 +0,0 @@
-// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
-#pragma once
-#include<yade/extra/Shop.hpp>
-
-#include<yade/core/InteractionSolver.hpp>
-#include<yade/core/FileGenerator.hpp>
-#include<yade/pkg-common/RigidBodyParameters.hpp>
-#include<yade/pkg-dem/BodyMacroParameters.hpp>
-#include<yade/pkg-common/InteractionPhysicsEngineUnit.hpp>
-#include<yade/pkg-dem/SpheresContactGeometry.hpp>
-#include<yade/pkg-common/GLDrawFunctors.hpp>
-#include<yade/pkg-common/PeriodicEngines.hpp>
-#include<yade/pkg-common/NormalShearInteractions.hpp>
-#include<yade/pkg-common/ConstitutiveLaw.hpp>
-
-
-/* Engine encompassing several computations looping over all bodies/interactions
- *
- * * Compute and store unbalanced force over the whole simulation.
- * * Compute and store volumetric strain for every body.
- *
- * May be extended in the future to compute global stiffness etc as well.
- */
-class CpmGlobalCharacteristics: public PeriodicEngine{
-	public:
-		bool useMaxForce; // use maximum unbalanced force instead of mean unbalanced force
-		Real unbalancedForce;
-		void compute(MetaBody* rb, bool useMax=false);
-		virtual void action(MetaBody* rb){compute(rb,useMaxForce);}
-		CpmGlobalCharacteristics(){};
-	REGISTER_ATTRIBUTES(PeriodicEngine,
-		(unbalancedForce)
-		(useMaxForce)
-	);
-	DECLARE_LOGGER;
-	REGISTER_CLASS_AND_BASE(CpmGlobalCharacteristics,PeriodicEngine);
-};
-REGISTER_SERIALIZABLE(CpmGlobalCharacteristics);
-
-/*! @brief representation of a single interaction of the CPM type: storage for relevant parameters.
- *
- * Evolution of the contact is governed by Law2_Dem3DofGeom_CpmPhys_Cpm:
- * that includes damage effects and chages of parameters inside CpmPhys.
- *
- */
-class CpmPhys: public NormalShearInteraction {
-	private:
-	public:
-		/*! Fundamental parameters (constants) */
-		Real
-			//! normal modulus (stiffness / crossSection)
-			E,
-			//! shear modulus
-			G,
-			//! tangens of internal friction angle
-			tanFrictionAngle, 
-			//! virgin material cohesion
-			undamagedCohesion,
-			//! equivalent cross-section associated with this contact
-			crossSection,
-			//! strain at which the material starts to behave non-linearly
-			epsCrackOnset,
-			//! strain where the damage-evolution law tangent from the top (epsCrackOnset) touches the axis;
-			/// since the softening law is exponential, this doesn't mean that the contact is fully damaged at this point,
-			/// that happens only asymptotically 
-			epsFracture,
-			//! damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞
-			omegaThreshold,
-			//! weight coefficient for shear strain when computing the strain semi-norm kappaD
-			xiShear,
-			//! characteristic time for damage (if non-positive, the law without rate-dependence is used)
-			dmgTau,
-			//! exponent in the rate-dependent damage evolution
-			dmgRateExp,
-			//! damage strain (at previous or current step)
-			dmgStrain,
-			//! damage viscous overstress (at previous step or at current step)
-			dmgOverstress,
-			//! characteristic time for viscoplasticity (if non-positive, no rate-dependence for shear)
-			plTau,
-			//! exponent in the rate-dependent viscoplasticity
-			plRateExp,
-			//! "prestress" of this link (used to simulate isotropic stress)
-			isoPrestress;
-		/*! Up to now maximum normal strain (semi-norm), non-decreasing in time. */
-		Real kappaD;
-		/*! Transversal strain (perpendicular to the contact axis) */
-		Real epsTrans;
-		/*! if not cohesive, interaction is deleted when distance is greater than zero. */
-		bool isCohesive;
-		/*! the damage evlution function will always return virgin state */
-		bool neverDamage;
-		/*! cummulative plastic strain measure (scalar) on this contact */
-		Real epsPlSum;
-		//! debugging, to see convergence rate
-		static long cummBetaIter, cummBetaCount;
-
-		/*! auxiliary variable for visualization, recalculated in Law2_Dem3DofGeom_CpmPhys_Cpm at every iteration */
-		// FIXME: Fn and Fs are stored as Vector3r normalForce, shearForce in NormalShearInteraction 
-		Real omega, Fn, sigmaN, epsN, relResidualStrength; Vector3r epsT, sigmaT, Fs;
-
-
-		static Real solveBeta(const Real c, const Real N);
-		Real computeDmgOverstress(Real dt);
-		Real computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt);
-
-
-
-		CpmPhys(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), isoPrestress(0.), kappaD(0.), epsTrans(0.), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; dmgOverstress=0; }
-		virtual ~CpmPhys();
-
-		REGISTER_ATTRIBUTES(NormalShearInteraction,
-			(E)
-			(G)
-			(tanFrictionAngle)
-			(undamagedCohesion)
-			(crossSection)
-			(epsCrackOnset)
-			(epsFracture)
-			(omegaThreshold)
-			(xiShear)
-			(dmgTau)
-			(dmgRateExp)
-			(dmgStrain)
-			(dmgOverstress)
-			(plTau)
-			(plRateExp)
-			(isoPrestress)
-
-			(cummBetaIter)
-			(cummBetaCount)
-
-			(kappaD)
-			(neverDamage)
-			(epsT)
-			(epsTrans)
-			(epsPlSum)
-
-			(isCohesive)
-
-			// auxiliary params to make them accessible from python
-			(omega)
-			(Fn)
-			(Fs)
-			(epsN)
-			(sigmaN)
-			(sigmaT)
-			(relResidualStrength)
-		);
-	REGISTER_CLASS_AND_BASE(CpmPhys,NormalShearInteraction);
-	DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(CpmPhys);
-
-/* This class holds information associated with each body */
-class CpmMat: public BodyMacroParameters {
-	public:
-		//! volumetric strain around this body
-		Real epsVolumetric;
-		//! number of (cohesive) contacts that damaged completely
-		int numBrokenCohesive;
-		//! number of contacts with this body
-		int numContacts;
-		//! average damage including already deleted contacts (it is really not damage, but 1-relResidualStrength now)
-		Real normDmg;
-		//! plastic strain on contacts already deleted
-		Real epsPlBroken;
-		//! sum of plastic strains normalized by number of contacts
-		Real normEpsPl;
-		CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.) {createIndex();};
-		REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl));
-		REGISTER_CLASS_AND_BASE(CpmMat,BodyMacroParameters);
-};
-REGISTER_SERIALIZABLE(CpmMat);
-
-class Law2_Dem3DofGeom_CpmPhys_Cpm: public ConstitutiveLaw{
-	public:
-	/*! Damage evolution law */
-	static Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
-		if(kappaD<epsCrackOnset || neverDamage) return 0;
-		return 1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
-	}
-		bool logStrain;
-		//! yield function: 0: mohr-coulomb (original); 1: parabolic; 2: logarithmic, 3: log+lin_tension, 4: elliptic, 5: elliptic+log
-		int yieldSurfType;
-		//! scaling in the logarithmic yield surface (should be <1 for realistic results; >=0 for meaningful results)
-		static Real yieldLogSpeed;
-		static Real yieldEllipseShift;
-		//! HACK: limit strain on some contacts by moving body #2 in the contact; only if refR1<0 (facet); deactivated if > 0
-		static Real minStrain_moveBody2;
-		Law2_Dem3DofGeom_CpmPhys_Cpm(): logStrain(false), yieldSurfType(0) { /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
-		void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
-	FUNCTOR2D(Dem3DofGeom,CpmPhys);
-	REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_CpmPhys_Cpm,ConstitutiveLaw);
-	REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain)(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(minStrain_moveBody2));
-	DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
-
-/*! @brief Convert macroscopic properties to CpmPhys with corresponding parameters.
- *
- * */
-class Ip2_CpmMat_CpmMat_CpmPhys: public InteractionPhysicsEngineUnit{
-	private:
-	public:
-		/* nonelastic material parameters */
-		/* alternatively (and more cleanly), we would have subclass of ElasticBodyParameters,
-		 * which would define just those in addition to the elastic ones.
-		 * This might be done later, for now hardcode that here. */
-		/* uniaxial tension resistance, bending parameter of the damage evolution law, whear weighting constant for epsT in the strain seminorm (kappa) calculation. Default to NaN so that user gets loudly notified it was not set.
-		
-		*/
-		Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, dmgTau, dmgRateExp, plTau, plRateExp, isoPrestress;
-		//! Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If 0, they will never be. If negative, they will always be created as cohesive.
-		long cohesiveThresholdIter;
-		//! Create contacts that don't receive any damage (CpmPhys::neverDamage=true); defaults to false
-		bool neverDamage;
-
-		Ip2_CpmMat_CpmMat_CpmPhys(){
-			// 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;
-			dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
-			omegaThreshold=0.999;
-			isoPrestress=0;
-		}
-
-		virtual void go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction);
-		REGISTER_ATTRIBUTES(InteractionPhysicsEngineUnit,
-			(cohesiveThresholdIter)
-			(G_over_E)
-			(sigmaT)
-			(neverDamage)
-			(epsCrackOnset)
-			(relDuctility)
-			(dmgTau)
-			(dmgRateExp)
-			(plTau)
-			(plRateExp)
-			(omegaThreshold)
-			(isoPrestress)
-		);
-
-		FUNCTOR2D(CpmMat,CpmMat);
-		REGISTER_CLASS_AND_BASE(Ip2_CpmMat_CpmMat_CpmPhys,InteractionPhysicsEngineUnit);
-		DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(Ip2_CpmMat_CpmMat_CpmPhys);
-
-class GLDrawCpmPhys: public GLDrawInteractionPhysicsFunctor {
-	public: virtual void go(const shared_ptr<InteractionPhysics>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
-	virtual ~GLDrawCpmPhys() {};
-	REGISTER_ATTRIBUTES(/*no base*/,(contactLine)(dmgLabel)(dmgPlane)(epsT)(epsTAxes)(normal)(colorStrain)(epsNLabel));
-	RENDERS(CpmPhys);
-	REGISTER_CLASS_AND_BASE(GLDrawCpmPhys,GLDrawInteractionPhysicsFunctor);
-	DECLARE_LOGGER;
-	static bool contactLine,dmgLabel,dmgPlane,epsT,epsTAxes,normal,colorStrain,epsNLabel;
-};
-REGISTER_SERIALIZABLE(GLDrawCpmPhys);
-
-class CpmPhysDamageColorizer: public PeriodicEngine {
-	public:
-		//! maximum damage over all contacts
-		Real maxOmega;
-		CpmPhysDamageColorizer(){maxOmega=0; /* run at the very beginning */ initRun=true;}
-		virtual void action(MetaBody*);
-	REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega));
-	REGISTER_CLASS_AND_BASE(CpmPhysDamageColorizer,PeriodicEngine);
-};
-REGISTER_SERIALIZABLE(CpmPhysDamageColorizer);
-

Deleted: trunk/extra/BrefcomTestGen.cpp
===================================================================
--- trunk/extra/BrefcomTestGen.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/BrefcomTestGen.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,130 +0,0 @@
-// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
-#include"BrefcomTestGen.hpp"
-#include<yade/core/MetaBody.hpp>
-#include<yade/pkg-dem/BodyMacroParameters.hpp>
-
-YADE_PLUGIN("BrefcomTestGen");
-
-
-
-
-#include<yade/pkg-common/PhysicalActionContainerReseter.hpp>
-#include<yade/pkg-common/InteractingSphere.hpp>
-#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
-#include<yade/pkg-common/AABB.hpp>
-#include<yade/pkg-common/InteractingSphere2AABB.hpp>
-#include<yade/pkg-common/MetaInteractingGeometry.hpp>
-#include<yade/pkg-common/MetaInteractingGeometry2AABB.hpp>
-#include<yade/pkg-common/InteractionGeometryMetaEngine.hpp>
-#include<yade/pkg-common/InteractionPhysicsMetaEngine.hpp>
-#include<yade/pkg-common/PhysicalActionApplier.hpp>
-#include<yade/pkg-common/PhysicalParametersMetaEngine.hpp>
-#include<yade/pkg-common/NewtonsForceLaw.hpp>
-#include<yade/pkg-common/NewtonsMomentumLaw.hpp>
-#include<yade/pkg-common/LeapFrogPositionIntegrator.hpp>
-#include<yade/pkg-common/LeapFrogOrientationIntegrator.hpp>
-#include<yade/pkg-common/PersistentSAPCollider.hpp>
-#include<yade/pkg-common/ConstitutiveLawDispatcher.hpp>
-#include<yade/pkg-dem/PositionOrientationRecorder.hpp>
-#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
-#include<yade/pkg-dem/Dem3DofGeom_SphereSphere.hpp>
-#include<yade/extra/UniaxialStrainControlledTest.hpp>
-
-
-CREATE_LOGGER(BrefcomTestGen);
-
-void BrefcomTestGen::createEngines(){
-	rootBody->initializers.clear();
-
-	shared_ptr<BoundingVolumeMetaEngine> boundingVolumeDispatcher	= shared_ptr<BoundingVolumeMetaEngine>(new BoundingVolumeMetaEngine);
-	boundingVolumeDispatcher->add(new InteractingSphere2AABB);
-	boundingVolumeDispatcher->add(new MetaInteractingGeometry2AABB);
-	rootBody->initializers.push_back(boundingVolumeDispatcher);
-
-	rootBody->engines.clear();
-
-	rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
-	rootBody->engines.push_back(boundingVolumeDispatcher);
-
-	shared_ptr<PersistentSAPCollider> collider(new PersistentSAPCollider);
-	collider->haveDistantTransient=true;
-	rootBody->engines.push_back(collider);
-
-	shared_ptr<InteractionGeometryMetaEngine> igeomDispatcher(new InteractionGeometryMetaEngine);
-	shared_ptr<ef2_Sphere_Sphere_Dem3DofGeom> ef2ssd3d(new ef2_Sphere_Sphere_Dem3DofGeom);
-	igeomDispatcher->add(ef2ssd3d);
-	rootBody->engines.push_back(igeomDispatcher);
-
-	shared_ptr<InteractionPhysicsMetaEngine> iphysDispatcher(new InteractionPhysicsMetaEngine);
-		shared_ptr<Ip2_CpmMat_CpmMat_CpmPhys> bmc(new Ip2_CpmMat_CpmMat_CpmPhys);
-		bmc->cohesiveThresholdIter=-1; bmc->G_over_E=1; bmc->sigmaT=3e9; bmc->neverDamage=true; bmc->epsCrackOnset=1e-4; bmc->relDuctility=5;
-		//bmc->calibratedEpsFracture=.5; /* arbitrary, but large enough */
-		iphysDispatcher->add(bmc);
-	rootBody->engines.push_back(iphysDispatcher);
-
-	shared_ptr<ConstitutiveLawDispatcher> clDisp(new ConstitutiveLawDispatcher);
-		clDisp->add(shared_ptr<ConstitutiveLaw>(new Law2_Dem3DofGeom_CpmPhys_Cpm));
-	rootBody->engines.push_back(clDisp);
-
-	shared_ptr<PhysicalActionApplier> applyActionDispatcher(new PhysicalActionApplier);
-	applyActionDispatcher->add(new NewtonsForceLaw);
-	applyActionDispatcher->add(new NewtonsMomentumLaw);
-	rootBody->engines.push_back(applyActionDispatcher);
-	
-	shared_ptr<PhysicalParametersMetaEngine> positionIntegrator(new PhysicalParametersMetaEngine);
-	positionIntegrator->add(new LeapFrogPositionIntegrator);
-	rootBody->engines.push_back(positionIntegrator);
-
-	shared_ptr<PhysicalParametersMetaEngine> orientationIntegrator(new PhysicalParametersMetaEngine);
-	orientationIntegrator->add(new LeapFrogOrientationIntegrator);
-	rootBody->engines.push_back(orientationIntegrator);
-
-	shared_ptr<CpmPhysDamageColorizer> dmg(new CpmPhysDamageColorizer);
-	rootBody->engines.push_back(dmg);
-
-#if 0
-	shared_ptr<BrefcomStiffnessCounter> bsc(new BrefcomStiffnessCounter);
-	bsc->interval=100;
-	rootBody->engines.push_back(bsc);
-
-	shared_ptr<GlobalStiffnessTimeStepper> gsts(new GlobalStiffnessTimeStepper);
-	gsts->sdecGroupMask=1023;
-	gsts->timeStepUpdateInterval=100;
-	gsts->defaultDt=1e-4;
-	rootBody->engines.push_back(gsts);
-#endif
-}
-
-bool BrefcomTestGen::generate(){
-	message="";
-	rootBody=Shop::rootBody();
-
-	createEngines();
-
-	shared_ptr<UniaxialStrainer> strainer(new UniaxialStrainer);
-	strainer->strainRate=strainRate;
-	strainer->axis=2; // z-oriented straining
-	strainer->limitStrain=-4;
-	rootBody->engines.push_back(strainer);
-	
-	// control normal/shear ratio
-	//Real zCoord=.1; Real yCoord=sqrt(1-zCoord*zCoord); // distance is always 2, with contact at origin
-	Real zCoord=.9, yCoord=0;
-	shared_ptr<Body>
-		s1=Shop::sphere(Vector3r(0,-yCoord,-zCoord),.5),
-		s2=Shop::sphere(Vector3r(0,yCoord,zCoord),.5),
-		sMid=Shop::sphere(Vector3r(0,0,0.01),.5);
-	body_id_t id1=rootBody->bodies->insert(s1), id2=rootBody->bodies->insert(s2); //, id3=rootBody->bodies->insert(sMid);
-	
-	//  engines should take care of the rest of interaction; this is what collider would do normally
-	/*
-	rootBody->transientInteractions->insert(id1,id2);
-	rootBody->transientInteractions->find(id1,id2)->isReal=1;
-	rootBody->transientInteractions->find(id1,id2)->isNew=1;
-	*/
-
-	strainer->negIds.push_back(id1); strainer->negCoords.push_back(-zCoord);
-	strainer->posIds.push_back(id2); strainer->posCoords.push_back(zCoord);
-
-	return true;
-}

Deleted: trunk/extra/BrefcomTestGen.hpp
===================================================================
--- trunk/extra/BrefcomTestGen.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/BrefcomTestGen.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,22 +0,0 @@
-#pragma once
-#include<yade/extra/Shop.hpp>
-#include<yade/core/FileGenerator.hpp>
-#include<yade/extra/Brefcom.hpp>
-
-class BrefcomTestGen: public FileGenerator {
-	private:
-		void createEngines();
-	public:
-		Real strainRate;
-		BrefcomTestGen(){strainRate=1e-2;};
-		bool generate();
-	protected :
-		void registerAttributes(){
-			FileGenerator::registerAttributes();
-			REGISTER_ATTRIBUTE(strainRate);
-		}
-	REGISTER_CLASS_NAME(BrefcomTestGen);
-	REGISTER_BASE_CLASS_NAME(FileGenerator);
-	DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(BrefcomTestGen);

Modified: trunk/extra/SConscript
===================================================================
--- trunk/extra/SConscript	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/SConscript	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,13 +1,7 @@
+# vim: set filetype=python :
 Import('*')
 
 import os.path, os
-brefcomMaterialModel='../brefcom-mm.hh'
-if os.path.exists('../'+brefcomMaterialModel):
-	print "Will include local file "+brefcomMaterialModel
-	brefcomInclude=['-include',brefcomMaterialModel]
-	Depends('Brefcom.cpp','../../brefcom-mm.hh')
-else:
-	brefcomInclude=['']
 
 env.Install('$PREFIX/lib/yade$SUFFIX/extra',[
 	env.SharedLibrary('SphericalDEMSimulator',
@@ -37,16 +31,8 @@
 	
 	env.SharedLibrary('TetraTestGen',['tetra/TetraTestGen.cpp'],LIBS=env['LIBS']+['Shop','Tetra']),
 
-	env.SharedLibrary('UniaxialStrainControlledTest',['usct/UniaxialStrainControlledTest.cpp'],LIBS=env['LIBS']+['Shop','ConstitutiveLawDispatcher','GlobalStiffnessTimeStepper','Brefcom','PositionOrientationRecorder']),
-
-	env.SharedLibrary('Brefcom',['Brefcom.cpp'],CXXFLAGS=env['CXXFLAGS']+brefcomInclude,LIBS=env['LIBS']+['Shop','InteractingSphere2InteractingSphere4SpheresContactGeometry','DemXDofGeom']),
-
-	env.SharedLibrary('BrefcomTestGen',['BrefcomTestGen.cpp'],LIBS=env['LIBS']+['Shop','UniaxialStrainControlledTest','PositionOrientationRecorder','Brefcom','ConstitutiveLawDispatcher','Dem3DofGeom_SphereSphere']),
-
-	env.SharedLibrary('SimpleScene',['SimpleScene.cpp'],LIBS=env['LIBS']+['Shop','SimpleElasticRelationships']),
-
 	env.SharedLibrary('Shop',
-		['clump/Shop.cpp'],
+		['Shop.cpp'],
 		LIBS=env['LIBS']+[
 			'ElasticContactLaw',
 			'ElasticCohesiveLaw',
@@ -74,7 +60,7 @@
 			'ParticleParameters',
 			'MetaInteractingGeometry',
 			'MetaInteractingGeometry2AABB',
-			'InteractingSphere2AABB',\
+			'InteractingSphere2AABB',
 			'InteractingBox2AABB',
 			'NewtonsMomentumLaw',
 			'NewtonsForceLaw',

Copied: trunk/extra/Shop.cpp (from rev 1771, trunk/extra/clump/Shop.cpp)


Property changes on: trunk/extra/Shop.cpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: trunk/extra/Shop.hpp (from rev 1767, trunk/extra/clump/Shop.hpp)


Property changes on: trunk/extra/Shop.hpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Deleted: trunk/extra/SimpleScene.cpp
===================================================================
--- trunk/extra/SimpleScene.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/SimpleScene.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,157 +0,0 @@
-#include"SimpleScene.hpp"
-#include<yade/extra/Shop.hpp>
-#include<yade/pkg-dem/BodyMacroParameters.hpp>
-#include<yade/pkg-common/InteractingSphere2AABB.hpp>
-#include<yade/pkg-common/InteractingBox2AABB.hpp>
-#include<yade/pkg-common/MetaInteractingGeometry.hpp>
-#include<yade/pkg-common/MetaInteractingGeometry2AABB.hpp>
-#include<yade/pkg-common/PhysicalActionContainerReseter.hpp>
-#include<yade/pkg-common/PhysicalParametersMetaEngine.hpp>
-#include<yade/pkg-common/AABB.hpp>
-#include<yade/pkg-common/Box.hpp>
-#include<yade/pkg-common/Sphere.hpp>
-#include<yade/pkg-common/InteractingBox.hpp>
-#include<yade/pkg-common/NewtonsForceLaw.hpp>
-#include<yade/pkg-common/NewtonsMomentumLaw.hpp>
-#include<yade/pkg-common/LeapFrogPositionIntegrator.hpp>
-#include<yade/pkg-common/LeapFrogOrientationIntegrator.hpp>
-#include<yade/pkg-dem/InteractingSphere2InteractingSphere4SpheresContactGeometry.hpp>
-#include<yade/pkg-dem/InteractingBox2InteractingSphere4SpheresContactGeometry.hpp>
-#include<yade/pkg-common/PhysicalParametersMetaEngine.hpp>
-#include<yade/pkg-common/InteractionGeometryMetaEngine.hpp>
-#include<yade/pkg-common/InteractionPhysicsMetaEngine.hpp>
-#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
-#include<yade/pkg-common/PhysicalActionDamper.hpp>
-#include<yade/pkg-common/PhysicalActionApplier.hpp>
-#include<yade/pkg-common/CundallNonViscousDamping.hpp>
-#include<yade/pkg-common/CundallNonViscousDamping.hpp>
-#include<yade/pkg-common/GravityEngines.hpp>
-#include<yade/pkg-dem/ElasticContactLaw.hpp>
-#include<yade/pkg-dem/SpheresContactGeometry.hpp>
-#include<yade/pkg-dem/SimpleElasticRelationships.hpp>
-#include<yade/pkg-common/PersistentSAPCollider.hpp>
-
-
-
-YADE_PLUGIN("SimpleScene");
-
-bool SimpleScene::generate(){
-	message="";
-	//@
-	rootBody=Shop::rootBody();
-	//@
-	/* initializers */
-		rootBody->initializers.clear();
-		//@
-		shared_ptr<BoundingVolumeMetaEngine> boundingVolumeDispatcher	= shared_ptr<BoundingVolumeMetaEngine>(new BoundingVolumeMetaEngine);
-			boundingVolumeDispatcher->add(new InteractingSphere2AABB);
-			boundingVolumeDispatcher->add(new InteractingBox2AABB);
-			boundingVolumeDispatcher->add(new MetaInteractingGeometry2AABB);
-			rootBody->initializers.push_back(boundingVolumeDispatcher);
-	//@
-	/* engines */
-		rootBody->engines.clear();
-		//@
-		rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
-		//@
-		// use boundingVolumeDispatcher that we defined above
-		rootBody->engines.push_back(boundingVolumeDispatcher);
-		//@
-		shared_ptr<PersistentSAPCollider> collider(new PersistentSAPCollider);
-			rootBody->engines.push_back(collider);
-		//@
-		shared_ptr<InteractionGeometryMetaEngine> igeomDispatcher(new InteractionGeometryMetaEngine);
-			igeomDispatcher->add(new InteractingSphere2InteractingSphere4SpheresContactGeometry);
-			igeomDispatcher->add(new InteractingBox2InteractingSphere4SpheresContactGeometry);
-			rootBody->engines.push_back(igeomDispatcher);
-		//@
-		shared_ptr<InteractionPhysicsMetaEngine> iphysDispatcher(new InteractionPhysicsMetaEngine);
-			iphysDispatcher->add(new SimpleElasticRelationships);
-			rootBody->engines.push_back(iphysDispatcher);
-		//@
-		shared_ptr<ElasticContactLaw> ecl(new ElasticContactLaw);
-			rootBody->engines.push_back(ecl);
-		//@
-		shared_ptr<GravityEngine> ge(new GravityEngine);
-			ge->gravity=Vector3r(0,0,-9.81);
-			rootBody->engines.push_back(ge);
-		//@
-		shared_ptr<PhysicalActionDamper> dampingDispatcher(new PhysicalActionDamper);
-			shared_ptr<CundallNonViscousForceDamping> forceDamper(new CundallNonViscousForceDamping);
-			forceDamper->damping=0.2;
-			dampingDispatcher->add(forceDamper);
-			shared_ptr<CundallNonViscousMomentumDamping> momentumDamper(new CundallNonViscousMomentumDamping);
-			momentumDamper->damping=0.2;
-			dampingDispatcher->add(momentumDamper);
-			rootBody->engines.push_back(dampingDispatcher);
-		//@
-		shared_ptr<PhysicalActionApplier> applyActionDispatcher(new PhysicalActionApplier);
-			applyActionDispatcher->add(new NewtonsForceLaw);
-			applyActionDispatcher->add(new NewtonsMomentumLaw);
-			rootBody->engines.push_back(applyActionDispatcher);
-		//@
-		shared_ptr<PhysicalParametersMetaEngine> positionIntegrator(new PhysicalParametersMetaEngine);
-			positionIntegrator->add(new LeapFrogPositionIntegrator);
-			rootBody->engines.push_back(positionIntegrator);
-		//@
-		shared_ptr<PhysicalParametersMetaEngine> orientationIntegrator(new PhysicalParametersMetaEngine);
-			orientationIntegrator->add(new LeapFrogOrientationIntegrator);
-			rootBody->engines.push_back(orientationIntegrator);
-	//@
-	// set default values for Shop	
-	Shop::setDefault("body_sdecGroupMask",1);
-	Shop::setDefault("phys_density",2400.); Shop::setDefault("phys_young",30e9); Shop::setDefault("phys_poisson",.3);
-	Shop::setDefault("aabb_randomColor",false);Shop::setDefault("shape_randomColor",false); Shop::setDefault("mold_randomColor",false);
-	Shop::setDefault("aabb_color",Vector3r(0,1,0)); Shop::setDefault("shape_color",Vector3r(1,0,0)); Shop::setDefault("mold_color",Vector3r(1,0,0));
-
-	//@
-	shared_ptr<Body> box=Shop::box(Vector3r(0,0,0),Vector3r(.5,.5,.5));
-	box->isDynamic=false;
-	rootBody->bodies->insert(box);
-	
-	//@
-	if(false){
-		Vector3r extents(.5,.5,.5);
-		shared_ptr<Body> b=shared_ptr<Body>(new Body(body_id_t(0),0));
-		b->isDynamic=true;
-		
-		// phys
-		shared_ptr<BodyMacroParameters> physics(new BodyMacroParameters);
-		physics->mass=8*extents[0]*extents[1]*extents[2]*2400;
-		physics->inertia=Vector3r(physics->mass*(4*extents[1]*extents[1]+4*extents[2]*extents[2])/12.,physics->mass*(4*extents[0]*extents[0]+4*extents[2]*extents[2])/12.,physics->mass*(4*extents[0]*extents[0]+4*extents[1]*extents[1])/12.);
-		physics->se3=Se3r(Vector3r(0,0,0),Quaternionr::IDENTITY);
-		physics->young=30e9;
-		physics->poisson=.3;
-		b->physicalParameters=physics;
-
-		// aabb
-		shared_ptr<AABB> aabb(new AABB);
-		aabb->diffuseColor=Vector3r(0,1,0);
-		b->boundingVolume=aabb;
-
-		//shape
-		shared_ptr<Box> shape(new Box);
-		shape->extents=extents;
-		shape->diffuseColor=Vector3r(1,0,0);
-		b->geometricalModel=shape;
-
-		// mold
-		shared_ptr<InteractingBox> mold(new InteractingBox);
-		mold->extents=extents;
-		mold->diffuseColor=Vector3r(1,0,0);
-		b->interactingGeometry=mold;
-
-		rootBody->bodies->insert(b);
-	}
-
-	//@
-	Shop::setDefault("shape_color",Vector3r(0,1,0)); Shop::setDefault("mold_color",Vector3r(0,1,0));
-	shared_ptr<Body> sphere(Shop::sphere(Vector3r(0,0,2),1));
-	rootBody->bodies->insert(sphere);
-	
-	//@
-	rootBody->dt=.2*Shop::PWaveTimeStep(rootBody);
-
-	//@
-	return true;
-}

Deleted: trunk/extra/SimpleScene.hpp
===================================================================
--- trunk/extra/SimpleScene.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/extra/SimpleScene.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,18 +0,0 @@
-#pragma once
-#include<yade/core/FileGenerator.hpp>
-#include<yade/extra/Shop.hpp>
-
-
-class SimpleScene: public FileGenerator {
-	public:
-		SimpleScene(){};
-		~SimpleScene (){};
-		virtual bool generate();
-	protected :
-		void registerAttributes(){ FileGenerator::registerAttributes(); }
-	REGISTER_CLASS_NAME(SimpleScene);
-	REGISTER_BASE_CLASS_NAME(FileGenerator);
-	DECLARE_LOGGER;
-};
-REGISTER_SERIALIZABLE(SimpleScene);
-

Modified: trunk/gui/SConscript
===================================================================
--- trunk/gui/SConscript	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/gui/SConscript	2009-05-24 18:22:30 UTC (rev 1778)
@@ -74,8 +74,8 @@
 			],
 			),
 		env.SharedLibrary('_utils',['py/_utils.cpp'],SHLIBPREFIX='',
-			LIBS=env['LIBS']+['Shop','Brefcom']),
-		env.SharedLibrary('_eudoxos',['py/_eudoxos.cpp'],SHLIBPREFIX='',CXXFLAGS=env['CXXFLAGS']+([] if not os.path.exists('../../brefcom-mm.hh') else ['-include','../brefcom-mm.hh']),LIBS=env['LIBS']+['Shop','Brefcom']),
+			LIBS=env['LIBS']+['Shop','ConcretePM']),
+		env.SharedLibrary('_eudoxos',['py/_eudoxos.cpp'],SHLIBPREFIX='',CXXFLAGS=env['CXXFLAGS']+([] if not os.path.exists('../../brefcom-mm.hh') else ['-include','../brefcom-mm.hh']),LIBS=env['LIBS']+['Shop','ConcretePM']),
 		env.SharedLibrary('log',['py/log.cpp'],SHLIBPREFIX=''),
 		env.File('__init__.py','py'),
 		env.File('utils.py','py'),

Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/gui/py/PythonUI_rc.py	2009-05-24 18:22:30 UTC (rev 1778)
@@ -34,20 +34,22 @@
 	for p in listChildClassesRecursive(root):
 		class argStorage:
 			def __init__(self,_root,_class): self._root,self._class=_root,_class
-			def wrap(self,*args): return yade.wrapper.__dict__[self._root](self._class,*args)
-		if root=='MetaEngine': _dd[p]=argStorage(root2,p).wrap
+			def __call__(self,*args): return yade.wrapper.__dict__[self._root](self._class,*args)
+		if root=='MetaEngine': _dd[p]=argStorage(root2,p)
 		else:                  _dd[p]=lambda __r_=root2,__p_=p,**kw : yade.wrapper.__dict__[__r_](__p_,kw) # eval(root2)(p,kw)
 ### end wrappers
 
 #### HANDLE RENAMED CLASSES ####
+# if old class name is used, the new object is constructed and a warning is issued about old name being used
 renamed={
+	# renamed 23.5.2009, may be removed in a few months
 	'BrefcomMakeContact':'Ip2_CpmMat_CpmMat_CpmPhys',
 	'BrefcomContact':'CpmPhys',
 	'BrefcomPhysParams':'CpmMat',
 	'ef2_Spheres_Brefcom_BrefcomLaw':'Law2_Dem3DofGeom_CpmPhys_Cpm',
 	'GLDrawBrefcomContact':'GLDrawCpmPhys',
 	'BrefcomDamageColorizer':'CpmPhysDamageColorizer',
-	'BrefcomGlobalCharacteristics':'CpmGlobalCharacteristics'
+	'BrefcomGlobalCharacteristics':'CpmGlobalCharacteristics',
 }
 
 for oldName in renamed:

Modified: trunk/gui/py/_eudoxos.cpp
===================================================================
--- trunk/gui/py/_eudoxos.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/gui/py/_eudoxos.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -1,4 +1,4 @@
-#include<yade/extra/Brefcom.hpp>
+#include<yade/pkg-dem/ConcretePM.hpp>
 #include<boost/python.hpp>
 #include<yade/extra/boost_python_len.hpp>
 using namespace boost::python;

Modified: trunk/gui/qt3/SimulationController.cpp
===================================================================
--- trunk/gui/qt3/SimulationController.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/gui/qt3/SimulationController.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -21,7 +21,6 @@
 #include <boost/filesystem/convenience.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
 #include <unistd.h>
-#include<yade/core/yadeExceptions.hpp>
 #include <Wm3Math.h>
 #include<yade/lib-base/yadeWm3.hpp>
 #include<boost/foreach.hpp>
@@ -242,7 +241,7 @@
 			pbResetSimulation->setDisabled(true);
 			pbOneSimulationStep->setDisabled(true);
 		}
-		catch(yadeError& e)
+		catch(std::runtime_error& e)
 		{
 			Omega::instance().resetRootBody();
 			Omega::instance().setSimulationFileName("");

Modified: trunk/lib/factory/ClassFactory.cpp
===================================================================
--- trunk/lib/factory/ClassFactory.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/lib/factory/ClassFactory.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -14,11 +14,6 @@
 
 class Factorable;
 
-void ClassFactory::addBaseDirectory(const string& dir)
-{
-	dlm.addBaseDirectory(dir);
-}
-
 bool ClassFactory::registerFactorable( std::string name 			   , CreateFactorableFnPtr create,
 					 CreateSharedFactorableFnPtr createShared, CreatePureCustomFnPtr createPureCustom)
 {
@@ -49,7 +44,7 @@
 	FactorableCreatorsMap::const_iterator i = map.find( name );
 	if( i == map.end() )
 	{
-		dlm.loadFromDirectoryList(name);
+		dlm.load(name);
 		if (dlm.isLoaded(name))
 		{
 			if( map.find( name ) == map.end() )
@@ -76,7 +71,7 @@
 	if( i == map.end() )
 	{
 		//cerr << "------------ going to load something" << endl;
-		dlm.loadFromDirectoryList(name);
+		dlm.load(name);
 		if (dlm.isLoaded(name))
 		{
 			if( map.find( name ) == map.end() )
@@ -105,7 +100,7 @@
 
 bool ClassFactory::load(const string& name)
 {
-        return dlm.loadFromDirectoryList(name);
+        return dlm.load(name);
 }
 
 string ClassFactory::lastError()

Modified: trunk/lib/factory/ClassFactory.hpp
===================================================================
--- trunk/lib/factory/ClassFactory.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/lib/factory/ClassFactory.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -137,8 +137,6 @@
 		*/
 		bool isFactorable(const type_info& tp,bool& fundamental);
 
-		void addBaseDirectory(const string& dir);
-
 		bool load(const string& name);
 		std::string lastError();
 

Modified: trunk/lib/factory/DynLibManager.cpp
===================================================================
--- trunk/lib/factory/DynLibManager.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/lib/factory/DynLibManager.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -26,41 +26,16 @@
 
 DynLibManager::DynLibManager ()
 {
-	baseDirs.clear();
 	autoUnload = true;
 }
 
 
 DynLibManager::~DynLibManager ()
 {
-	if (autoUnload)
-		unloadAll();
+	if(autoUnload) unloadAll();
 }
 
 
-void DynLibManager::addBaseDirectory(const string& dir)
-{
-	string tmpDir;
-	if ( dir[dir.size()-1]=='/' || dir[dir.size()-1]=='\\' )
-		tmpDir = dir.substr(0,dir.size()-1);
-	else
-		tmpDir = dir;
-
-	baseDirs.push_back(tmpDir);
-}
-
-
-bool DynLibManager::loadFromDirectoryList (const string& libName )
-{
-	if (libName.empty()) return false;
-	string libFileName = libNameToSystemName(libName);
-	string baseDir = findLibDir(libName);
-	string fullLibName;
-	if (baseDir.length()==0) return load(libFileName,libName);
-	else return load(baseDir+"/"+libFileName,libName);
-}
-
-
 bool DynLibManager::load (const string& fullLibName, const string& libName )
 {
 	if (libName.empty() || fullLibName.empty()){
@@ -159,8 +134,6 @@
 	
 		if (lastError != ERROR_SUCCESS)
 		{
-//			cerr << errMsg << endl;
-//			Omega::printErrorLog(errMsg);
 			lastError_ = errMsg;
 			return true;
 		}
@@ -170,8 +143,6 @@
  		char * error = dlerror();
 		if (error != NULL)  
 		{
-			//cerr << error << endl;
-			//Omega::printErrorLog(error);
 			lastError_ = error;
 		}
 		return (error!=NULL);
@@ -195,8 +166,8 @@
 {
 	string libName;
 	if(name.length()<=3){ // this arbitrary value may disappear once the logic below is dumped...
-		LOG_WARN("Filename `"<<name<<"' too short, returning empty string (cross thumbs).");
-		return "";
+		// LOG_WARN("Filename `"<<name<<"' too short, returning empty string (cross thumbs).");
+		return "[Garbage plugin file `"+name+"']";
 	}
 
 	#ifdef WIN32
@@ -208,26 +179,3 @@
 	return libName;
 }
 
-
-string DynLibManager::findLibDir(const string& name)
-{
-	string libFileName = libNameToSystemName(name);
-
-	string baseDir;
-	baseDir.clear();
-
-	vector<string>::iterator bdi    = baseDirs.begin();
-	vector<string>::iterator bdiEnd = baseDirs.end();
-	for( ; bdi != bdiEnd ; ++bdi)
-	{
-		filesystem::path name = filesystem::path((*bdi)+"/"+libFileName);
-		if ( filesystem::exists( name ) )
-		{
-			baseDir = (*bdi);
-			break;
-		}
-	}
-
-	return baseDir;
-}
-

Modified: trunk/lib/factory/DynLibManager.hpp
===================================================================
--- trunk/lib/factory/DynLibManager.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/lib/factory/DynLibManager.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -36,8 +36,6 @@
 		#else	
 		std::map<const std::string, void *> handles;
 		#endif
-
-		vector<string> baseDirs;
 		bool autoUnload;
 
 	public :
@@ -45,9 +43,10 @@
 		~DynLibManager ();
 		void addBaseDirectory(const std::string& dir);
 
-		bool load (const std::string& libName, const std::string& libName2);
-		bool loadFromDirectoryList (const std::string& fullLibName);
+		bool load(const std::string& libName, const std::string& libName2);
+		bool load(const std::string& pluginName){return load(libNameToSystemName(pluginName),pluginName);}
 
+
 		bool unload (const string libName);
 		bool isLoaded (const string libName);
 		bool unloadAll ();
@@ -55,7 +54,6 @@
 
 		string libNameToSystemName(const string& name);
 		string systemNameToLibName(const string& name);
-		string findLibDir(const string& name);
 		string lastError();
 		DECLARE_LOGGER;
 

Copied: trunk/pkg/dem/ConcretePM.cpp (from rev 1777, trunk/extra/Brefcom.cpp)
===================================================================
--- trunk/extra/Brefcom.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/ConcretePM.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -0,0 +1,364 @@
+// 2007,2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+#include"ConcretePM.hpp"
+#include<yade/core/MetaBody.hpp>
+#include<yade/pkg-dem/BodyMacroParameters.hpp>
+#include<yade/pkg-common/Sphere.hpp>
+#include<yade/lib-QGLViewer/qglviewer.h>
+#include<yade/lib-opengl/GLUtils.hpp>
+#include<yade/pkg-dem/DemXDofGeom.hpp>
+#include<yade/extra/Shop.hpp>
+
+YADE_PLUGIN("CpmMat","Ip2_CpmMat_CpmMat_CpmPhys","CpmPhys","Law2_Dem3DofGeom_CpmPhys_Cpm","CpmGlobalCharacteristics","GLDrawCpmPhys","CpmPhysDamageColorizer");
+
+
+/********************** Ip2_CpmMat_CpmMat_CpmPhys ****************************/
+
+CREATE_LOGGER(Ip2_CpmMat_CpmMat_CpmPhys);
+
+void Ip2_CpmMat_CpmMat_CpmPhys::go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction){
+	Dem3DofGeom* contGeom=YADE_CAST<Dem3DofGeom*>(interaction->interactionGeometry.get());
+
+	assert(contGeom); // for now, don't handle anything other than SpheresContactGeometry and Dem3DofGeom
+
+	if(!interaction->isNew && interaction->interactionPhysics){ /* relax */ } 
+	else {
+		interaction->isNew=false; // just in case
+
+		const shared_ptr<BodyMacroParameters>& elast1=static_pointer_cast<BodyMacroParameters>(pp1);
+		const shared_ptr<BodyMacroParameters>& elast2=static_pointer_cast<BodyMacroParameters>(pp2);
+
+		Real E12=2*elast1->young*elast2->young/(elast1->young+elast2->young); // harmonic Young's modulus average
+		//Real nu12=2*elast1->poisson*elast2->poisson/(elast1->poisson+elast2->poisson); // dtto for Poisson ratio 
+		Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
+		Real S12=Mathr::PI*pow(minRad,2); // "surface" of interaction
+		//Real E=(E12 /* was here for Kn:  *S12/d0  */)*((1+alpha)/(beta*(1+nu12)+gamma*(1-alpha*nu12)));
+		//Real E=E12; // apply alpha, beta, gamma: garbage values of E !?
+
+		if(!neverDamage) { assert(!isnan(sigmaT)); }
+
+		shared_ptr<CpmPhys> contPhys(new CpmPhys());
+
+		contPhys->E=E12;
+		contPhys->G=E12*G_over_E;
+		contPhys->tanFrictionAngle=tan(.5*(elast1->frictionAngle+elast2->frictionAngle));
+		contPhys->undamagedCohesion=sigmaT;
+		contPhys->crossSection=S12;
+		contPhys->epsCrackOnset=epsCrackOnset;
+		contPhys->epsFracture=relDuctility*epsCrackOnset;
+		contPhys->omegaThreshold=omegaThreshold;
+		// inherited from NormalShearInteracion, used in the timestepper
+		contPhys->kn=contPhys->E*contPhys->crossSection;
+		contPhys->ks=contPhys->G*contPhys->crossSection;
+
+		if(neverDamage) contPhys->neverDamage=true;
+		if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
+		else contPhys->isCohesive=false;
+		contPhys->dmgTau=dmgTau;
+		contPhys->dmgRateExp=dmgRateExp;
+		contPhys->plTau=plTau;
+		contPhys->plRateExp=plRateExp;
+		contPhys->isoPrestress=isoPrestress;
+
+		interaction->interactionPhysics=contPhys;
+	}
+}
+
+
+
+
+/********************** CpmPhys ****************************/
+CREATE_LOGGER(CpmPhys);
+
+// !! at least one virtual function in the .cpp file
+CpmPhys::~CpmPhys(){};
+
+CREATE_LOGGER(Law2_Dem3DofGeom_CpmPhys_Cpm);
+
+long CpmPhys::cummBetaIter=0, CpmPhys::cummBetaCount=0;
+
+Real CpmPhys::solveBeta(const Real c, const Real N){
+	#ifdef YADE_DEBUG
+		cummBetaCount++;
+	#endif
+	const int maxIter=20;
+	const Real maxError=1e-12;
+	Real f, ret=0.;
+	for(int i=0; i<maxIter; i++){
+		#ifdef YADE_DEBUG
+			cummBetaIter++;
+		#endif
+		Real aux=c*exp(N*ret)+exp(ret);
+		f=log(aux);
+		if(fabs(f)<maxError) return ret;
+		Real df=(c*N*exp(N*ret)+exp(ret))/aux;
+		ret-=f/df;
+	}
+	LOG_FATAL("No convergence after "<<maxIter<<" iters; c="<<c<<", N="<<N<<", ret="<<ret<<", f="<<f);
+	throw runtime_error("Law2_Dem3DofGeom_CpmPhys_Cpm::solveBeta failed to converge.");
+}
+
+Real CpmPhys::computeDmgOverstress(Real dt){
+	if(dmgStrain>=epsN*omega){ // unloading, no viscous stress
+		dmgStrain=epsN*omega;
+		LOG_TRACE("Elastic/unloading, no viscous overstress");
+		return 0.;
+	}
+	Real c=epsCrackOnset*(1-omega)*pow(dmgTau/dt,dmgRateExp)*pow(epsN*omega-dmgStrain,dmgRateExp-1.);
+	Real beta=solveBeta(c,dmgRateExp);
+	Real deltaDmgStrain=(epsN*omega-dmgStrain)*exp(beta);
+	dmgStrain+=deltaDmgStrain;
+	LOG_TRACE("deltaDmgStrain="<<deltaDmgStrain<<", viscous overstress "<<(epsN*omega-dmgStrain)*E);
+	/* σN=Kn(εN-εd); dmgOverstress=σN-(1-ω)*Kn*εN=…=Kn(ω*εN-εd) */
+	return (epsN*omega-dmgStrain)*E;
+}
+
+Real CpmPhys::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt){
+	if(sigmaTNorm<sigmaTYield) return 1.;
+	Real c=undamagedCohesion*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-1.);
+	Real beta=solveBeta(c,plRateExp);
+	//LOG_DEBUG("scaling factor "<<1.-exp(beta)*(1-sigmaTYield/sigmaTNorm));
+	return 1.-exp(beta)*(1-sigmaTYield/sigmaTNorm);
+}
+
+
+
+/********************** Law2_Dem3DofGeom_CpmPhys_Cpm ****************************/
+
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2=1.; /* deactivated if > 0 */
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed=1.;
+Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift=0.;
+
+void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody){
+	//timingDeltas->start();
+	Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
+	CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
+
+	/* kept fully damaged contacts; note that normally the contact is deleted _after_ the CPM_MATERIAL_MODEL,
+	 * i.e. if it is 1.0 here, omegaThreshold is >= 1.0 for sure.
+	 * &&'ing that just to make sure anyway ...
+	 */
+	// if(BC->omega>=1.0 && BC->omegaThreshold>=1.0) return;
+
+	// shorthands
+	Real& epsN(BC->epsN); Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(BC->omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& dt=Omega::instance().getTimeStep();  const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau); const bool& isCohesive(BC->isCohesive);
+	/* const Real& transStrainCoeff(BC->transStrainCoeff); const Real& epsTrans(BC->epsTrans); const Real& xiShear(BC->xiShear); */
+	Real& omega(BC->omega); Real& sigmaN(BC->sigmaN);  Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
+	const Real& yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed); const int& yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType);
+	const Real& yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift); 
+
+	#define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); throw;}
+	
+	#define NNAN(a) YADE_VERIFY(!isnan(a));
+	#define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
+
+	//timingDeltas->checkpoint("setup");
+	// if(contGeom->refR1<0) contGeom->refLength=contGeom->refR2; // make facet-sphere contact always at equilibrium when touching exactly (and not the initial distance)
+	epsN=contGeom->strainN(); epsT=contGeom->strainT();
+	#ifdef YADE_DEBUG
+		if(isnan(epsN)){
+			LOG_FATAL("refLength="<<contGeom->refLength<<"; pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<"; displacementN="<<contGeom->displacementN());
+			throw runtime_error("!! epsN==NaN !!");
+		}
+		NNAN(epsN); NNANV(epsT);
+	#endif
+	// already in SpheresContactGeometry:
+	// contGeom->relocateContactPoints(); // allow very large mutual rotations
+	if(logStrain && epsN<0){
+		Real epsN0=epsN;
+		epsN=log(epsN0+1); epsT*=epsN/epsN0;
+	}
+	NNAN(epsN); NNANV(epsT);
+	//timingDeltas->checkpoint("geom");
+
+	epsN+=BC->isoPrestress/E;
+	//TRVAR1(epsN);
+	#ifdef CPM_MATERIAL_MODEL
+		CPM_MATERIAL_MODEL
+	#else
+		sigmaN=E*epsN;
+		sigmaT=G*epsT;
+	#endif
+	sigmaN-=BC->isoPrestress;
+	if(contGeom->refR1<0 && Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2<=0 && epsN<Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2){
+		/* move Body2 (the sphere) so that minStrain is satisfied */
+		rootBody->bex.addMove(I->getId2(),contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
+		LOG_TRACE("Moving by "<<contGeom->normal*(Law2_Dem3DofGeom_CpmPhys_Cpm::minStrain_moveBody2-epsN)*contGeom->refLength);
+	}
+	NNAN(kappaD); NNAN(epsCrackOnset); NNAN(epsFracture); NNAN(omega);
+	NNAN(sigmaN); NNANV(sigmaT); NNAN(crossSection);
+	//timingDeltas->checkpoint("material");
+
+	//const int watch1=6300, watch2=6299;
+	//#define SHOW(a) if((I->getId1()==watch1 && I->getId2()==watch2) || (I->getId2()==watch1 && I->getId1()==watch2)) cerr<<__FILE__<<":"<<__LINE__<<" "<<a<<endl;
+	//SHOW("epsN"<<epsN);
+	if(epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)){
+		rootBody->interactions->requestErase(I->getId1(),I->getId2());
+		if(isCohesive){
+			const shared_ptr<Body>& body1=Body::byId(I->getId1(),rootBody), body2=Body::byId(I->getId2(),rootBody); assert(body1); assert(body2);
+			const shared_ptr<CpmMat>& rbp1=YADE_PTR_CAST<CpmMat>(body1->physicalParameters), rbp2=YADE_PTR_CAST<CpmMat>(body2->physicalParameters);
+			if(BC->isCohesive){rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum;}
+			LOG_DEBUG("Contact #"<<I->getId1()<<"=#"<<I->getId2()<<" is damaged over thershold ("<<omega<<">"<<omegaThreshold<<") and will be deleted.");
+		}
+		return;
+	}
+
+	Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
+	Fs=sigmaT*crossSection; BC->shearForce=Fs;
+
+	applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, rootBody);
+	//timingDeltas->checkpoint("rest");
+}
+
+
+/********************** GLDrawCpmPhys ****************************/
+
+#include<yade/lib-opengl/OpenGLWrapper.hpp>
+
+CREATE_LOGGER(GLDrawCpmPhys);
+
+bool GLDrawCpmPhys::contactLine=true;
+bool GLDrawCpmPhys::dmgLabel=true;
+bool GLDrawCpmPhys::dmgPlane=false;
+bool GLDrawCpmPhys::epsNLabel=true;
+bool GLDrawCpmPhys::epsT=false;
+bool GLDrawCpmPhys::epsTAxes=false;
+bool GLDrawCpmPhys::normal=false;
+bool GLDrawCpmPhys::colorStrain=false;
+
+
+void GLDrawCpmPhys::go(const shared_ptr<InteractionPhysics>& ip, const shared_ptr<Interaction>& i, const shared_ptr<Body>& b1, const shared_ptr<Body>& b2, bool wireFrame){
+	const shared_ptr<CpmPhys>& BC=static_pointer_cast<CpmPhys>(ip);
+	const shared_ptr<Dem3DofGeom>& geom=YADE_PTR_CAST<Dem3DofGeom>(i->interactionGeometry);
+
+	//Vector3r lineColor(BC->omega,1-BC->omega,0.0); /* damaged links red, undamaged green */
+	Vector3r lineColor=Shop::scalarOnColorScale(1.-BC->relResidualStrength);
+
+	if(colorStrain) lineColor=Vector3r(
+		min((Real)1.,max((Real)0.,-BC->epsTrans/BC->epsCrackOnset)),
+		min((Real)1.,max((Real)0.,BC->epsTrans/BC->epsCrackOnset)),
+		min((Real)1.,max((Real)0.,abs(BC->epsTrans)/BC->epsCrackOnset-1)));
+
+	if(contactLine) GLUtils::GLDrawLine(b1->physicalParameters->dispSe3.position,b2->physicalParameters->dispSe3.position,lineColor);
+	if(dmgLabel){ GLUtils::GLDrawNum(BC->omega,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
+	else if(epsNLabel){ GLUtils::GLDrawNum(BC->epsN,0.5*(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position),lineColor); }
+	if(BC->omega>0 && dmgPlane){
+		Real halfSize=sqrt(1-BC->relResidualStrength)*.5*.705*sqrt(BC->crossSection);
+		Vector3r midPt=.5*Vector3r(b1->physicalParameters->dispSe3.position+b2->physicalParameters->dispSe3.position);
+		glDisable(GL_CULL_FACE);
+		glPushMatrix();
+			glTranslatev(midPt);
+			Quaternionr q; q.Align(Vector3r::UNIT_Z,geom->normal);
+			Vector3r axis; Real angle; q.ToAxisAngle(axis,angle);
+			glRotatef(angle*Mathr::RAD_TO_DEG,axis[0],axis[1],axis[2]);
+			glBegin(GL_POLYGON);
+				glColor3v(lineColor); 
+				glVertex3d(halfSize,0.,0.);
+				glVertex3d(.5*halfSize,.866*halfSize,0.);
+				glVertex3d(-.5*halfSize,.866*halfSize,0.);
+				glVertex3d(-halfSize,0.,0.);
+				glVertex3d(-.5*halfSize,-.866*halfSize,0.);
+				glVertex3d(.5*halfSize,-.866*halfSize,0.);
+			glEnd();
+		glPopMatrix();
+	}
+
+	const Vector3r& cp=static_pointer_cast<Dem3DofGeom>(i->interactionGeometry)->contactPoint;
+	if(epsT){
+		Real maxShear=(BC->undamagedCohesion-BC->sigmaN*BC->tanFrictionAngle)/BC->G;
+		Real relShear=BC->epsT.Length()/maxShear;
+		Real scale=.5*geom->refLength;
+		Vector3r dirShear=BC->epsT; dirShear.Normalize();
+		if(epsTAxes){
+			GLUtils::GLDrawLine(cp-Vector3r(scale,0,0),cp+Vector3r(scale,0,0));
+			GLUtils::GLDrawLine(cp-Vector3r(0,scale,0),cp+Vector3r(0,scale,0));
+			GLUtils::GLDrawLine(cp-Vector3r(0,0,scale),cp+Vector3r(0,0,scale));
+		}
+		GLUtils::GLDrawArrow(cp,cp+dirShear*relShear*scale,Vector3r(1.,0.,0.));
+		GLUtils::GLDrawLine(cp+dirShear*relShear*scale,cp+dirShear*scale,Vector3r(.3,.3,.3));
+
+		/* normal strain */ GLUtils::GLDrawArrow(cp,cp+geom->normal*(BC->epsN/maxShear),Vector3r(0.,1.,0.));
+	}
+	//if(normal) GLUtils::GLDrawArrow(cp,cp+geom->normal*.5*BC->equilibriumDist,Vector3r(0.,1.,0.));
+}
+
+
+/********************** CpmGlobalCharacteristics ****************************/
+
+CREATE_LOGGER(CpmGlobalCharacteristics);
+void CpmGlobalCharacteristics::compute(MetaBody* rb, bool useMaxForce){
+	rb->bex.sync();
+
+	// 1. reset volumetric strain (cummulative in the next loop)
+	// 2. get maximum force on a body and sum of all forces (for averaging)
+	Real sumF=0,maxF=0,currF;
+	FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+	CpmMat* bpp(YADE_CAST<CpmMat*>(b->physicalParameters.get()));
+		bpp->epsVolumetric=0;
+		bpp->numContacts=0;
+		currF=rb->bex.getForce(b->id).Length(); maxF=max(currF,maxF); sumF+=currF;
+	}
+	Real meanF=sumF/rb->bodies->size(); 
+
+	// commulate normal strains from contacts
+	// get max force on contacts
+	Real maxContactF=0;
+	FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
+		if(!I->isReal) continue;
+		shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
+		maxContactF=max(maxContactF,max(BC->Fn,BC->Fs.Length()));
+		CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
+		CpmMat* bpp2(YADE_CAST<CpmMat*>(Body::byId(I->getId2())->physicalParameters.get()));
+		bpp1->epsVolumetric+=BC->epsN; bpp1->numContacts+=1;
+		bpp2->epsVolumetric+=BC->epsN; bpp2->numContacts+=1;
+	}
+	unbalancedForce=(useMaxForce?maxF:meanF)/maxContactF;
+
+	FOREACH(const shared_ptr<Interaction>& I, *rb->interactions){
+		if(!I->isReal) continue;
+		shared_ptr<CpmPhys> BC=YADE_PTR_CAST<CpmPhys>(I->interactionPhysics); assert(BC);
+		CpmMat* bpp1(YADE_CAST<CpmMat*>(Body::byId(I->getId1())->physicalParameters.get()));
+		CpmMat* bpp2(YADE_CAST<CpmMat*>(Body::byId(I->getId2())->physicalParameters.get()));
+		Real epsVolAvg=.5*((3./bpp1->numContacts)*bpp1->epsVolumetric+(3./bpp2->numContacts)*bpp2->epsVolumetric);
+		BC->epsTrans=(epsVolAvg-BC->epsN)/2.;
+		//TRVAR5(I->getId1(),I->getId2(),BC->epsTrans,(3./bpp1->numContacts)*bpp1->epsVolumetric,(3./bpp2->numContacts)*bpp2->epsVolumetric);
+		//TRVAR4(bpp1->numContacts,bpp1->epsVolumetric,bpp2->numContacts,bpp2->epsVolumetric);
+	}
+	#if 0
+		FOREACH(const shared_ptr<Body>& b, *rb->bodies){
+			CpmMat* bpp(YADE_PTR_CAST<CpmMat>(b->physicalParameters.get()));
+			bpp->epsVolumeric*=3/bpp->numContacts;
+		}
+	#endif
+}
+
+
+/********************** CpmPhysDamageColorizer ****************************/
+void CpmPhysDamageColorizer::action(MetaBody* rootBody){
+	//vector<pair<short,Real> > bodyDamage; /* number of cohesive interactions per body; cummulative damage of interactions */
+	//vector<pair<short,
+	vector<BodyStats> bodyStats; bodyStats.resize(rootBody->bodies->size());
+	assert(bodyStats[0].nCohLinks==0); // should be initialized by dfault ctor
+	FOREACH(const shared_ptr<Interaction>& I, *rootBody->interactions){
+		shared_ptr<CpmPhys> BC=dynamic_pointer_cast<CpmPhys>(I->interactionPhysics);
+		if(!BC || !BC->isCohesive) continue;
+		const body_id_t id1=I->getId1(), id2=I->getId2();
+		bodyStats[id1].nCohLinks++; bodyStats[id1].dmgSum+=(1-BC->relResidualStrength); bodyStats[id1].epsPlSum+=BC->epsPlSum;
+		bodyStats[id2].nCohLinks++; bodyStats[id2].dmgSum+=(1-BC->relResidualStrength); bodyStats[id2].epsPlSum+=BC->epsPlSum;
+		//bodyDamage[id1].first++; bodyDamage[id2].first++;
+		//bodyDamage[id1].second+=(1-BC->relResidualStrength); bodyDamage[id2].second+=(1-BC->relResidualStrength);
+		maxOmega=max(maxOmega,BC->omega);
+	}
+	FOREACH(shared_ptr<Body> B, *rootBody->bodies){
+		body_id_t id=B->getId();
+		// add damaged contacts that have already been deleted
+		CpmMat* bpp=dynamic_cast<CpmMat*>(B->physicalParameters.get());
+		if(!bpp) continue;
+		short cohLinksWhenever=bodyStats[id].nCohLinks+bpp->numBrokenCohesive;
+		if(cohLinksWhenever>0){
+			bpp->normDmg=(bodyStats[id].dmgSum+bpp->numBrokenCohesive)/cohLinksWhenever;
+			bpp->normEpsPl=(bodyStats[id].epsPlSum+bpp->epsPlBroken)/cohLinksWhenever;
+		}
+		else { bpp->normDmg=0; bpp->normEpsPl=0;}
+		B->geometricalModel->diffuseColor=Vector3r(bpp->normDmg,1-bpp->normDmg,B->isDynamic?0:1);
+	}
+}


Property changes on: trunk/pkg/dem/ConcretePM.cpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: trunk/pkg/dem/ConcretePM.hpp (from rev 1777, trunk/extra/Brefcom.hpp)
===================================================================
--- trunk/extra/Brefcom.hpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/ConcretePM.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -0,0 +1,318 @@
+// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+
+/*
+
+=== HIGH LEVEL OVERVIEW OF CPM ===
+
+Concrete Particle Model (ConcretePM, Cpm) is a set of classes for modelling
+mechanical behavior of concrete. Several classes are needed for Cpm.
+
+1. CpmMat (Cpm material) deriving from BodyMacroParameters, which additionally has
+   some information about damage on the body, cummulative plastic strain etc.
+
+2.	Ip2_CpmMat_CpmMat_CpmPhys is 2-ary functor for creating CpmPhys from CpmMat's of
+	2 bodies that collide. Some parameters of the CpmPhys created are computed from
+	CpmMat's, others are passed as parameters of the functor.
+
+3. CpmPhys (Cpm (interaction)Physics) holds various parameters as well as internal
+   variables of the contact that can change as result of plasticity, damage, viscosity.
+
+4. Law2_Dem3Dof_CpmPhys_Cpm is constitutive law that takes geometry of the interaction
+	(Dem3Dof, which can be either Dem3Dof_SphereSphere or Dem3Dof_FacetSphere) and
+	CpmPhys, computing forces on both bodies and updating contact variables.
+
+	The model itself is defined in the macro CPM_MATERIAL_MODEL, but due to 
+	commercial reasons, those about 30 lines of code cannot be disclosed now and the macro
+	is defined in an external file. The model will be, however, described in enough detail
+	in my thesis (once it is written), along
+	with calibration procedures; it features damage, plasticity and viscosity
+	and is quite tunable (rigidity, poisson's	ratio, compressive/tensile strength
+	ratio, fracture energy, behavior under confinement, rate-dependence).
+
+There are other classes, which are not strictly necessary:
+
+ * CpmGlobalCharacteristics computes a few information about individual bodies based on
+   interactions they are involved in. It is probably quite useless now since volumetricStrain
+	is not used in the constitutive law anymore.
+
+ * GLDrawCpmPhys draws interaction physics (color for damage and a few other); rarely used, though.
+
+ * CpmPhysDamageColorizer changes bodies' colors depending on average damage of their interactions
+   and number of interactions that were already fully broken and have disappeared. This engine
+	contains its own loop (2 loops, more precisely) over all bodies and is run periodically
+	to update colors.
+
+*/
+
+#pragma once
+
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<yade/pkg-dem/BodyMacroParameters.hpp>
+#include<yade/pkg-common/InteractionPhysicsEngineUnit.hpp>
+#include<yade/pkg-dem/SpheresContactGeometry.hpp>
+#include<yade/pkg-common/GLDrawFunctors.hpp>
+#include<yade/pkg-common/PeriodicEngines.hpp>
+#include<yade/pkg-common/NormalShearInteractions.hpp>
+#include<yade/pkg-common/ConstitutiveLaw.hpp>
+
+/* This class holds information associated with each body */
+class CpmMat: public BodyMacroParameters {
+	public:
+		//! volumetric strain around this body
+		Real epsVolumetric;
+		//! number of (cohesive) contacts that damaged completely
+		int numBrokenCohesive;
+		//! number of contacts with this body
+		int numContacts;
+		//! average damage including already deleted contacts (it is really not damage, but 1-relResidualStrength now)
+		Real normDmg;
+		//! plastic strain on contacts already deleted
+		Real epsPlBroken;
+		//! sum of plastic strains normalized by number of contacts
+		Real normEpsPl;
+		CpmMat(): epsVolumetric(0.), numBrokenCohesive(0), numContacts(0), normDmg(0.), epsPlBroken(0.), normEpsPl(0.) {createIndex();};
+		REGISTER_ATTRIBUTES(BodyMacroParameters, (epsVolumetric) (numBrokenCohesive) (numContacts) (normDmg) (epsPlBroken) (normEpsPl));
+		REGISTER_CLASS_AND_BASE(CpmMat,BodyMacroParameters);
+};
+REGISTER_SERIALIZABLE(CpmMat);
+
+
+/*! @brief representation of a single interaction of the CPM type: storage for relevant parameters.
+ *
+ * Evolution of the contact is governed by Law2_Dem3DofGeom_CpmPhys_Cpm:
+ * that includes damage effects and chages of parameters inside CpmPhys.
+ *
+ */
+class CpmPhys: public NormalShearInteraction {
+	private:
+	public:
+		/*! Fundamental parameters (constants) */
+		Real
+			//! normal modulus (stiffness / crossSection)
+			E,
+			//! shear modulus
+			G,
+			//! tangens of internal friction angle
+			tanFrictionAngle, 
+			//! virgin material cohesion
+			undamagedCohesion,
+			//! equivalent cross-section associated with this contact
+			crossSection,
+			//! strain at which the material starts to behave non-linearly
+			epsCrackOnset,
+			//! strain where the damage-evolution law tangent from the top (epsCrackOnset) touches the axis;
+			/// since the softening law is exponential, this doesn't mean that the contact is fully damaged at this point,
+			/// that happens only asymptotically 
+			epsFracture,
+			//! damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞
+			omegaThreshold,
+			//! weight coefficient for shear strain when computing the strain semi-norm kappaD
+			xiShear,
+			//! characteristic time for damage (if non-positive, the law without rate-dependence is used)
+			dmgTau,
+			//! exponent in the rate-dependent damage evolution
+			dmgRateExp,
+			//! damage strain (at previous or current step)
+			dmgStrain,
+			//! damage viscous overstress (at previous step or at current step)
+			dmgOverstress,
+			//! characteristic time for viscoplasticity (if non-positive, no rate-dependence for shear)
+			plTau,
+			//! exponent in the rate-dependent viscoplasticity
+			plRateExp,
+			//! "prestress" of this link (used to simulate isotropic stress)
+			isoPrestress;
+		/*! Up to now maximum normal strain (semi-norm), non-decreasing in time. */
+		Real kappaD;
+		/*! Transversal strain (perpendicular to the contact axis) */
+		Real epsTrans;
+		/*! if not cohesive, interaction is deleted when distance is greater than zero. */
+		bool isCohesive;
+		/*! the damage evlution function will always return virgin state */
+		bool neverDamage;
+		/*! cummulative plastic strain measure (scalar) on this contact */
+		Real epsPlSum;
+		//! debugging, to see convergence rate
+		static long cummBetaIter, cummBetaCount;
+
+		/*! auxiliary variable for visualization, recalculated in Law2_Dem3DofGeom_CpmPhys_Cpm at every iteration */
+		// Fn and Fs are also stored as Vector3r normalForce, shearForce in NormalShearInteraction 
+		Real omega, Fn, sigmaN, epsN, relResidualStrength; Vector3r epsT, sigmaT, Fs;
+
+
+		static Real solveBeta(const Real c, const Real N);
+		Real computeDmgOverstress(Real dt);
+		Real computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt);
+
+
+
+		CpmPhys(): NormalShearInteraction(),E(0), G(0), tanFrictionAngle(0), undamagedCohesion(0), crossSection(0), xiShear(0), dmgTau(-1), dmgRateExp(0), dmgStrain(0), plTau(-1), plRateExp(0), isoPrestress(0.), kappaD(0.), epsTrans(0.), epsPlSum(0.) { createIndex(); epsT=Vector3r::ZERO; isCohesive=false; neverDamage=false; omega=0; Fn=0; Fs=Vector3r::ZERO; epsPlSum=0; dmgOverstress=0; }
+		virtual ~CpmPhys();
+
+		REGISTER_ATTRIBUTES(NormalShearInteraction,
+			(E)
+			(G)
+			(tanFrictionAngle)
+			(undamagedCohesion)
+			(crossSection)
+			(epsCrackOnset)
+			(epsFracture)
+			(omegaThreshold)
+			(xiShear)
+			(dmgTau)
+			(dmgRateExp)
+			(dmgStrain)
+			(dmgOverstress)
+			(plTau)
+			(plRateExp)
+			(isoPrestress)
+
+			(cummBetaIter)
+			(cummBetaCount)
+
+			(kappaD)
+			(neverDamage)
+			(epsT)
+			(epsTrans)
+			(epsPlSum)
+
+			(isCohesive)
+
+			// auxiliary params to make them accessible from python
+			(omega)
+			(Fn)
+			(Fs)
+			(epsN)
+			(sigmaN)
+			(sigmaT)
+			(relResidualStrength)
+		);
+	REGISTER_CLASS_AND_BASE(CpmPhys,NormalShearInteraction);
+	DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(CpmPhys);
+
+
+/*! @brief Convert macroscopic properties to CpmPhys with corresponding parameters.
+ *
+ * */
+class Ip2_CpmMat_CpmMat_CpmPhys: public InteractionPhysicsEngineUnit{
+	private:
+	public:
+		/* nonelastic material parameters */
+		/* alternatively (and more cleanly), we would have subclass of ElasticBodyParameters,
+		 * which would define just those in addition to the elastic ones.
+		 * This might be done later, for now hardcode that here. */
+		/* uniaxial tension resistance, bending parameter of the damage evolution law, whear weighting constant for epsT in the strain seminorm (kappa) calculation. Default to NaN so that user gets loudly notified it was not set.
+		
+		*/
+		Real sigmaT, epsCrackOnset, relDuctility, G_over_E, tau, expDmgRate, omegaThreshold, dmgTau, dmgRateExp, plTau, plRateExp, isoPrestress;
+		//! Should new contacts be cohesive? They will before this iter#, they will not be afterwards. If 0, they will never be. If negative, they will always be created as cohesive.
+		long cohesiveThresholdIter;
+		//! Create contacts that don't receive any damage (CpmPhys::neverDamage=true); defaults to false
+		bool neverDamage;
+
+		Ip2_CpmMat_CpmMat_CpmPhys(){
+			// 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;
+			dmgTau=-1; dmgRateExp=0; plTau=-1; plRateExp=-1;
+			omegaThreshold=0.999;
+			isoPrestress=0;
+		}
+
+		virtual void go(const shared_ptr<PhysicalParameters>& pp1, const shared_ptr<PhysicalParameters>& pp2, const shared_ptr<Interaction>& interaction);
+		REGISTER_ATTRIBUTES(InteractionPhysicsEngineUnit,
+			(cohesiveThresholdIter)
+			(G_over_E)
+			(sigmaT)
+			(neverDamage)
+			(epsCrackOnset)
+			(relDuctility)
+			(dmgTau)
+			(dmgRateExp)
+			(plTau)
+			(plRateExp)
+			(omegaThreshold)
+			(isoPrestress)
+		);
+
+		FUNCTOR2D(CpmMat,CpmMat);
+		REGISTER_CLASS_AND_BASE(Ip2_CpmMat_CpmMat_CpmPhys,InteractionPhysicsEngineUnit);
+		DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Ip2_CpmMat_CpmMat_CpmPhys);
+
+
+
+class Law2_Dem3DofGeom_CpmPhys_Cpm: public ConstitutiveLaw{
+	public:
+	/*! Damage evolution law */
+	static Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
+		if(kappaD<epsCrackOnset || neverDamage) return 0;
+		return 1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
+	}
+		bool logStrain;
+		//! yield function: 0: mohr-coulomb (original); 1: parabolic; 2: logarithmic, 3: log+lin_tension, 4: elliptic, 5: elliptic+log
+		int yieldSurfType;
+		//! scaling in the logarithmic yield surface (should be <1 for realistic results; >=0 for meaningful results)
+		static Real yieldLogSpeed;
+		static Real yieldEllipseShift;
+		//! HACK: limit strain on some contacts by moving body #2 in the contact; only if refR1<0 (facet); deactivated if > 0
+		static Real minStrain_moveBody2;
+		Law2_Dem3DofGeom_CpmPhys_Cpm(): logStrain(false), yieldSurfType(0) { /*timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);*/ }
+		void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, MetaBody* rootBody);
+	FUNCTOR2D(Dem3DofGeom,CpmPhys);
+	REGISTER_CLASS_AND_BASE(Law2_Dem3DofGeom_CpmPhys_Cpm,ConstitutiveLaw);
+	REGISTER_ATTRIBUTES(ConstitutiveLaw,(logStrain)(yieldSurfType)(yieldLogSpeed)(yieldEllipseShift)(minStrain_moveBody2));
+	DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
+
+/* Engine encompassing several computations looping over all bodies/interactions
+ *
+ * * Compute and store unbalanced force over the whole simulation.
+ * * Compute and store volumetric strain for every body.
+ *
+ * May be extended in the future to compute global stiffness etc as well.
+ */
+class CpmGlobalCharacteristics: public PeriodicEngine{
+	public:
+		bool useMaxForce; // use maximum unbalanced force instead of mean unbalanced force
+		Real unbalancedForce;
+		void compute(MetaBody* rb, bool useMax=false);
+		virtual void action(MetaBody* rb){compute(rb,useMaxForce);}
+		CpmGlobalCharacteristics(){};
+	REGISTER_ATTRIBUTES(PeriodicEngine,
+		(unbalancedForce)
+		(useMaxForce)
+	);
+	DECLARE_LOGGER;
+	REGISTER_CLASS_AND_BASE(CpmGlobalCharacteristics,PeriodicEngine);
+};
+REGISTER_SERIALIZABLE(CpmGlobalCharacteristics);
+
+class GLDrawCpmPhys: public GLDrawInteractionPhysicsFunctor {
+	public: virtual void go(const shared_ptr<InteractionPhysics>&,const shared_ptr<Interaction>&,const shared_ptr<Body>&,const shared_ptr<Body>&,bool wireFrame);
+	virtual ~GLDrawCpmPhys() {};
+	REGISTER_ATTRIBUTES(/*no base*/,(contactLine)(dmgLabel)(dmgPlane)(epsT)(epsTAxes)(normal)(colorStrain)(epsNLabel));
+	RENDERS(CpmPhys);
+	REGISTER_CLASS_AND_BASE(GLDrawCpmPhys,GLDrawInteractionPhysicsFunctor);
+	DECLARE_LOGGER;
+	static bool contactLine,dmgLabel,dmgPlane,epsT,epsTAxes,normal,colorStrain,epsNLabel;
+};
+REGISTER_SERIALIZABLE(GLDrawCpmPhys);
+
+class CpmPhysDamageColorizer: public PeriodicEngine {
+	struct BodyStats{ short nCohLinks; Real dmgSum; Real epsPlSum; BodyStats(): nCohLinks(0), dmgSum(0), epsPlSum(0.){} };
+	public:
+		//! maximum damage over all contacts
+		Real maxOmega;
+		CpmPhysDamageColorizer(){maxOmega=0; /* run at the very beginning */ initRun=true;}
+		virtual void action(MetaBody*);
+	REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega));
+	REGISTER_CLASS_AND_BASE(CpmPhysDamageColorizer,PeriodicEngine);
+};
+REGISTER_SERIALIZABLE(CpmPhysDamageColorizer);
+


Property changes on: trunk/pkg/dem/ConcretePM.hpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp (from rev 1777, trunk/extra/usct/UniaxialStrainControlledTest.cpp)
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -0,0 +1,164 @@
+// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+#include"UniaxialStrainer.hpp"
+#include<boost/foreach.hpp>
+
+#include<yade/core/InteractionContainer.hpp>
+#include<yade/pkg-common/ParticleParameters.hpp>
+#include<yade/pkg-common/AABB.hpp>
+
+YADE_PLUGIN("UniaxialStrainer");
+
+/************************ UniaxialStrainer **********************/
+CREATE_LOGGER(UniaxialStrainer);
+
+void UniaxialStrainer::init(){
+	needsInit=false;
+
+	assert(posIds.size()>0);
+	assert(negIds.size()>0);
+	posCoords.clear(); negCoords.clear();
+	FOREACH(body_id_t id,posIds){ const shared_ptr<Body>& b=Body::byId(id,rootBody); posCoords.push_back(b->physicalParameters->se3.position[axis]);
+		if(blockDisplacements && blockRotations) b->isDynamic=false;
+		else{
+			shared_ptr<PhysicalParameters> &pp=b->physicalParameters;
+			if(!blockDisplacements)pp->blockedDOFs=PhysicalParameters::axisDOF(axis); else pp->blockedDOFs=PhysicalParameters::DOF_XYZ;
+			if(blockRotations) pp->blockedDOFs|=PhysicalParameters::DOF_RXRYRZ;
+		}
+	}
+	FOREACH(body_id_t id,negIds){ const shared_ptr<Body>& b=Body::byId(id,rootBody); negCoords.push_back(b->physicalParameters->se3.position[axis]);
+		if(blockDisplacements && blockRotations) b->isDynamic=false;
+		else{
+			shared_ptr<PhysicalParameters> &pp=b->physicalParameters;
+			if(!blockDisplacements)pp->blockedDOFs=PhysicalParameters::axisDOF(axis); else pp->blockedDOFs=PhysicalParameters::DOF_XYZ;
+			if(blockRotations) pp->blockedDOFs|=PhysicalParameters::DOF_RXRYRZ;
+		}
+	}
+
+	assert(posIds.size()==posCoords.size() && negIds.size()==negCoords.size());
+
+	originalLength=axisCoord(posIds[0])-axisCoord(negIds[0]);
+	LOG_DEBUG("Reference particles: positive #"<<posIds[0]<<" at "<<axisCoord(posIds[0])<<"; negative #"<<negIds[0]<<" at "<<axisCoord(negIds[0]));
+	LOG_INFO("Setting initial length to "<<originalLength<<" (between #"<<negIds[0]<<" and #"<<posIds[0]<<")");
+	if(originalLength<=0) LOG_FATAL("Initial length is negative or zero (swapped reference particles?)! "<<originalLength);
+	/* this happens is nan propagates from e.g. brefcom consitutive law in case 2 bodies have _exactly_ the same position
+	 * (the the normal strain is 0./0.=nan). That is an user's error, however and should not happen. */
+	if(isnan(originalLength)) LOG_FATAL("Initial length is NaN!");
+	assert(originalLength>0 && !isnan(originalLength));
+
+	assert(!isnan(strainRate) || !isnan(absSpeed));
+	if(!isnan(std::numeric_limits<Real>::quiet_NaN())){ LOG_FATAL("NaN's are not properly supported (compiled, with -ffast-math?), which is required."); throw; }
+
+	if(isnan(strainRate)){ strainRate=absSpeed/originalLength; LOG_INFO("Computed new strainRate "<<strainRate); }
+	else {absSpeed=strainRate*originalLength;}
+
+	if(!setSpeeds){
+		initAccelTime_s=initAccelTime>=0 ? initAccelTime : Omega::instance().getTimeStep()*(-initAccelTime);
+		LOG_INFO("Strain speed will be "<<absSpeed<<", strain rate "<<strainRate<<", will be reached after "<<initAccelTime_s<<"s ("<<initAccelTime_s/Omega::instance().getTimeStep()<<" steps).");
+	} else {
+		/* set speed such that it is linear on the strained axis; transversal speed is not set, which can perhaps create some problems.
+			Note: all bodies in the simulation will have their speed set, since there is no way to tell which ones are part of the specimen
+			and which are not.
+
+			Speeds will be linearly interpolated beween axis positions p0,p1 and velocities v0,v1.
+		*/
+		initAccelTime_s=0;
+		LOG_INFO("Strain speed will be "<<absSpeed<<", strain rate "<<strainRate<<"; velocities will be set directly at the beginning.");
+		Real p0=axisCoord(negIds[0]), p1=axisCoord(posIds[0]); // limit positions
+		Real v0,v1; // speeds at p0, p1
+		switch(asymmetry){
+			case -1: v0=-absSpeed; v1=0; break;
+			case  0: v0=-absSpeed/2; v1=absSpeed/2; break;
+			case  1: v0=0; v1=absSpeed; break;
+			default: LOG_FATAL("Unknown asymmetry value "<<asymmetry<<" (should be -1,0,1)"); throw;
+		}
+		assert(p1>p0);
+		// set speeds for particles on the boundary
+		FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
+			// skip bodies on the boundary, since those will have their positions updated directly
+			if(std::find(posIds.begin(),posIds.end(),b->id)!=posIds.end() || std::find(negIds.begin(),negIds.end(),b->id)!=negIds.end()) { continue; }
+			Real p=axisCoord(b->id);
+			Real pNormalized=(p-p0)/(p1-p0);
+			YADE_CAST<ParticleParameters*>(b->physicalParameters.get())->velocity[axis]=pNormalized*(v1-v0)+v0;
+		}
+	}
+	stressUpdateInterval=max(1,(int)(2e-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;
+	 * this will not work if there are foreign bodies in the simulation,
+	 * in which case you must give the value yourself as engine attribute.
+	 *
+	 * A TODO option is to get crossSectionArea as average area of bounding boxes' of ABBBs
+	 * of posIds and negIds perpendicular to axis. That might be better, except for cases where
+	 * reference particles on either end do not coincide with the specimen cross-section.
+	 *
+	 * */
+	if(crossSectionArea<=0){
+		shared_ptr<AABB> rbAABB;
+		if (Omega::instance().getRootBody()->boundingVolume && (rbAABB=dynamic_pointer_cast<AABB>(Omega::instance().getRootBody()->boundingVolume))){
+			int axis2=(axis+1)%3, axis3=(axis+2)%3; // perpendicular axes indices
+			crossSectionArea=4*rbAABB->halfSize[axis2]*rbAABB->halfSize[axis3];
+			LOG_INFO("Setting crossSectionArea="<<crossSectionArea<<", using axes #"<<axis2<<" and #"<<axis3<<".");
+		} else {
+			crossSectionArea=1.;
+			LOG_WARN("No Axis Aligned Bounding Box for rootBody, using garbage value ("<<crossSectionArea<<") for crossSectionArea!");
+		}
+	}
+	assert(crossSectionArea>0);
+}
+
+void UniaxialStrainer::action(MetaBody* _rootBody){
+	rootBody=_rootBody;
+	if(needsInit) init();
+	// postconditions for initParams
+	assert(posIds.size()==posCoords.size() && negIds.size()==negCoords.size() && originalLength>0 && crossSectionArea>0);
+	//nothing to do
+	if(posIds.size()==0 || negIds.size()==0) return;
+	// linearly increase strain to the desired value
+	if(abs(currentStrainRate)<abs(strainRate)){
+		Real t=Omega::instance().getSimulationTime();
+		if(initAccelTime_s!=0) currentStrainRate=(t/initAccelTime_s)*strainRate;
+		else currentStrainRate=strainRate;
+	} else currentStrainRate=strainRate;
+	// how much do we move (in total, symmetry handled below)
+	Real dAX=currentStrainRate*originalLength*Omega::instance().getTimeStep();
+	if(!isnan(stopStrain)){
+		Real axialLength=axisCoord(posIds[0])-axisCoord(negIds[0]);
+		Real newStrain=(axialLength+dAX)/originalLength-1;
+		if((newStrain*stopStrain>0) && abs(newStrain)>=stopStrain){ // same sign of newStrain and stopStrain && over the limit from below in abs values
+			dAX=originalLength*(stopStrain+1)-axialLength;
+			LOG_INFO("Reached stopStrain "<<stopStrain<<", deactivating self and stopping in "<<idleIterations+1<<" iterations.");
+			this->active=false;
+			rootBody->stopAtIteration=Omega::instance().getCurrentIteration()+1+idleIterations;
+		}
+	}
+	if(asymmetry==0) dAX*=.5; // apply half on both sides if straining symetrically
+	for(size_t i=0; i<negIds.size(); i++){
+		if(asymmetry==0 || asymmetry==-1 /* for +1, don't move*/) negCoords[i]-=dAX;
+		axisCoord(negIds[i])=negCoords[i]; // update current position
+	}
+	for(size_t i=0; i<posIds.size(); i++){
+		if(asymmetry==0 || asymmetry==1 /* for -1, don't move */) posCoords[i]+=dAX;
+		axisCoord(posIds[i])=posCoords[i];
+	}
+
+	Real axialLength=axisCoord(posIds[0])-axisCoord(negIds[0]);
+	strain=axialLength/originalLength-1;
+
+	// reverse if we're over the limit strain
+	if(notYetReversed && limitStrain!=0 && ((currentStrainRate>0 && strain>limitStrain) || (currentStrainRate<0 && strain<limitStrain))) { currentStrainRate*=-1; notYetReversed=false; LOG_INFO("Reversed strain rate to "<<currentStrainRate); }
+
+	// update forces and stresses
+	if(Omega::instance().getCurrentIteration()%stressUpdateInterval==0) {
+		computeAxialForce();
+		avgStress=(sumPosForces+sumNegForces)/(2*crossSectionArea); // average nominal stress
+	}
+}
+
+void UniaxialStrainer::computeAxialForce(){
+	sumPosForces=sumNegForces=0;
+	rootBody->bex.sync();
+	FOREACH(body_id_t id, negIds) sumNegForces+=rootBody->bex.getForce(id)[axis];
+	FOREACH(body_id_t id, posIds) sumPosForces-=rootBody->bex.getForce(id)[axis];
+}
+


Property changes on: trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.cpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.hpp (from rev 1767, trunk/extra/usct/UniaxialStrainControlledTest.hpp)
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.hpp	2009-05-05 06:53:11 UTC (rev 1767)
+++ trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.hpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -0,0 +1,100 @@
+// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+#pragma once
+#include<fstream>
+#include<limits>
+#include<yade/extra/Shop.hpp>
+#include<yade/core/FileGenerator.hpp>
+#include<yade/core/StandAloneEngine.hpp>
+
+#ifndef FOREACH
+#define FOREACH BOOST_FOREACH
+#endif
+
+/*! Axial displacing two groups of bodies in the opposite direction with given strain rate.
+ *
+ * Takes two groups of body IDs (in posIds and negIds) and displaces them at each timestep in the direction given by axis∈{0,1,2} (for axes x,y,z respectively). These bodies automatically have Body::isDynamic==false.
+ *
+ * This engine should be run once forces on particles have been computed.
+ */
+class UniaxialStrainer: public StandAloneEngine {
+	private:
+		MetaBody* rootBody;
+		bool needsInit;
+
+		void computeAxialForce();
+		Real& axisCoord(body_id_t id){ return Body::byId(id,rootBody)->physicalParameters->se3.position[axis]; };
+		void init();
+	public:
+		virtual bool isActivated(){return active;}
+		//! strain rate, starting at 0, linearly raising to strainRate
+		Real strainRate,currentStrainRate;
+		//! alternatively, absolute speed of boundary motion can be specified; this is effective only at the beginning and if strainRate is not set; changing absSpeed directly during simulation wil have no effect.
+		Real absSpeed;
+		//! strain at which we will pause simulation; inactive (nan) by default; must be reached from below (in absolute value)
+		Real stopStrain;
+		//! distance of reference bodies in the direction of axis before straining started
+		Real originalLength;
+		//! invert the sense of straining (sharply, without transition) one this value of strain is reached. Not effective if 0.
+		Real limitStrain;
+		//! Flag whether the sense of straining has already been reversed
+		bool notYetReversed;
+		Real sumPosForces,sumNegForces;
+		//! crossSection perpendicular to he strained axis, computed from AABB of MetaBody
+		Real crossSectionArea;		//! Apply strain along x (0), y (1) or z(2) axis
+		//! The axis which is strained (0,1,2 for x,y,z)
+		int axis;
+		//! If 0, straining is symmetric for negIds and posIds; for 1 (or -1), only posIds are strained and negIds don't move (or vice versa)
+		int asymmetry;
+		//! Whether displacement of boundary bodies perpendicular to the strained axis are blocked of are free
+		bool blockDisplacements;
+		//! Whether rotations of boundary bodies are blocked.
+		bool blockRotations;
+		//! Are we activated?
+		bool active;
+		//! Number of iterations that will pass without straining activity after stopStrain has been reached (default: 0)
+		long idleIterations;
+		//! Time for strain reaching the requested value (linear interpolation). If negative, the time is dt*(-initAccelTime), where dt is  the timestep at the first iteration.
+		Real initAccelTime, initAccelTime_s /* value always in s, computed from initAccelTime */;
+		//! should we set speeds at the beginning directly, instead of increasing strain rate progressively?
+		bool setSpeeds;
+		//! how often to update forces (initialized automatically)
+		int stressUpdateInterval;
+
+		/** bodies on which straining will be applied (on the positive and negative side of axis) */
+		vector<body_id_t> posIds, negIds;
+		/** coordinates of pos/neg bodies in the direction of axis */
+		vector<Real> posCoords,negCoords;
+		//! Auxiliary vars (serializable, for recording)
+		Real strain, avgStress;
+
+		virtual void action(MetaBody*);
+		UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0; originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1; needsInit=true; strain=avgStress=0; blockRotations=false; blockDisplacements=false; setSpeeds=false; strainRate=absSpeed=stopStrain=numeric_limits<Real>::quiet_NaN(); active=true; idleIterations=0; initAccelTime=-200;};
+		virtual ~UniaxialStrainer(){};
+		REGISTER_ATTRIBUTES(StandAloneEngine,
+				(strainRate) 
+				(absSpeed)
+				(initAccelTime)
+				(stopStrain) 
+				(active)
+				(idleIterations)
+				(currentStrainRate) 
+				(axis) 
+				(asymmetry) 
+				(posIds) 
+				(negIds) 
+				(originalLength) 
+				(limitStrain) 
+				(notYetReversed) 
+				(crossSectionArea) 
+				(strain) 
+				(avgStress) 
+				(blockDisplacements) 
+				(blockRotations) 
+				(setSpeeds)
+		);
+	REGISTER_CLASS_AND_BASE(UniaxialStrainer,StandAloneEngine);
+	DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(UniaxialStrainer);
+
+


Property changes on: trunk/pkg/dem/Engine/StandAloneEngine/UniaxialStrainer.hpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Modified: trunk/pkg/dem/PreProcessor/HydraulicTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/HydraulicTest.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/PreProcessor/HydraulicTest.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -19,7 +19,6 @@
 #include<yade/pkg-dem/ElasticCriterionTimeStepper.hpp>
 #include<yade/pkg-dem/PositionOrientationRecorder.hpp>
 
-#include<yade/core/yadeExceptions.hpp>
 #include<yade/pkg-common/Box.hpp>
 #include<yade/pkg-common/AABB.hpp>
 #include<yade/pkg-common/Sphere.hpp>
@@ -118,7 +117,7 @@
 		{
 			message="Error: cannot load the file that should contain spheres"; return false;
 		}
-		catch ( yadeError& e )
+		catch ( std::runtime_error& e )
 		{
 			message="Error: cannot load the file that should contain spheres"; return false;
 		}

Copied: trunk/pkg/dem/PreProcessor/SimpleScene.cpp (from rev 1767, trunk/extra/SimpleScene.cpp)


Property changes on: trunk/pkg/dem/PreProcessor/SimpleScene.cpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Copied: trunk/pkg/dem/PreProcessor/SimpleScene.hpp (from rev 1767, trunk/extra/SimpleScene.hpp)


Property changes on: trunk/pkg/dem/PreProcessor/SimpleScene.hpp
___________________________________________________________________
Name: svn:mergeinfo
   + 

Modified: trunk/pkg/dem/PreProcessor/ThreePointBending.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/ThreePointBending.cpp	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/PreProcessor/ThreePointBending.cpp	2009-05-24 18:22:30 UTC (rev 1778)
@@ -22,7 +22,6 @@
 #include<yade/pkg-dem/ElasticCriterionTimeStepper.hpp>
 
 
-#include<yade/core/yadeExceptions.hpp>
 #include<yade/pkg-common/Box.hpp>
 #include<yade/pkg-common/AABB.hpp>
 #include<yade/pkg-common/Sphere.hpp>

Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript	2009-05-24 10:53:29 UTC (rev 1777)
+++ trunk/pkg/dem/SConscript	2009-05-24 18:22:30 UTC (rev 1778)
@@ -25,6 +25,15 @@
 				'CGAL'])
 		])
 
+cpmModel='../brefcom-mm.hh'
+if os.path.exists('../../'+cpmModel):
+	print "Will include local file "+cpmModel
+	cpmInclude=['-include',cpmModel]
+	Depends('ConcretePM.cpp','../../../brefcom-mm.hh')
+else:
+	cpmInclude=['']
+
+
 env.Install('$PREFIX/lib/yade$SUFFIX/pkg-dem',[
 	
 	env.SharedLibrary('DemXDofGeom',['DataClass/InteractionGeometry/DemXDofGeom.cpp']),
@@ -33,12 +42,20 @@
 
 	env.SharedLibrary('FacetTopologyAnalyzer',['Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp'],LIBS=env['LIBS']+['InteractingFacet']),
 
+	env.SharedLibrary('UniaxialStrainer',['Engine/StandAloneEngine/UniaxialStrainer.cpp'],LIBS=env['LIBS']+['ParticleParameters','AABB']),
+
+	env.SharedLibrary('ConcretePM',['ConcretePM.cpp'],CXXFLAGS=env['CXXFLAGS']+cpmInclude,LIBS=env['LIBS']+['Shop','DemXDofGeom']),
+
 	env.SharedLibrary('Clump',['DataClass/Clump.cpp'],LIBS=env['LIBS']+['Shop']),
 
 	env.SharedLibrary('SQLiteRecorder',
 		['Engine/StandAloneEngine/SQLiteRecorder.cpp'],
 		LIBS=env['LIBS']+['sqlite3x']),
 
+	env.SharedLibrary('SimpleScene',['PreProcessor/SimpleScene.cpp'],LIBS=env['LIBS']+['Shop','SimpleElasticRelationships']),
+
+	env.SharedLibrary('UniaxialStrainerGen',['PreProcessor/UniaxialStrainerGen.cpp'],LIBS=env['LIBS']+['Shop','ConstitutiveLawDispatcher','ConcretePM','Dem3DofGeom_SphereSphere','NewtonsDampedLaw','UniaxialStrainer']),
+
 	env.SharedLibrary('InteractingMyTetrahedron',
 		['DataClass/InteractingGeometry/InteractingMyTetrahedron.cpp']),