yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-06-15 19:34:09 +0200 (Mon, 15 Jun 2009)
New Revision: 1799

Modified:
   trunk/core/Interaction.cpp
   trunk/core/Interaction.hpp
   trunk/lib/smoothing/WeightedAverage2d.cpp
   trunk/lib/smoothing/WeightedAverage2d.hpp
   trunk/pkg/dem/ConcretePM.cpp
   trunk/pkg/dem/ConcretePM.hpp
Log:
1. Add resetting interactionGeometry and interactionPhysics to Interaction::reset(), as it should be...
2. Small things in CPM
3. Add GaussAverage.data method to python (for debugging and plotting)


Modified: trunk/core/Interaction.cpp
===================================================================
--- trunk/core/Interaction.cpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/core/Interaction.cpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -12,20 +12,25 @@
 
 #include<yade/core/MetaBody.hpp>
 
-Interaction::Interaction(): id1(0), id2(0){ reset(); }
+Interaction::Interaction(): id1(0), id2(0){ init(); }
 Interaction::Interaction(body_id_t newId1,body_id_t newId2): id1(newId1), id2(newId2){ reset(); }
 
 bool Interaction::isFresh(MetaBody* rb){ return iterMadeReal==rb->currentIteration;}
 
-
-void Interaction::reset(){
+void Interaction::init(){
 	isNeighbor = true;//NOTE : TriangulationCollider needs that
 	iterMadeReal=-1;
 	functorCache.geomExists=true;
 	//functorCache.geom=shared_ptr<InteractionGeometryEngineUnit>(); functorCache.phys=shared_ptr<InteractionPhysicsEngineUnit>(); functorCache.constLaw=shared_ptr<ConstitutiveLaw>();
 }
 
+void Interaction::reset(){
+	interactionGeometry=shared_ptr<InteractionGeometry>();
+	interactionPhysics=shared_ptr<InteractionPhysics>();
+	init();
+}
 
+
 void Interaction::swapOrder(){
 	if(interactionGeometry || interactionPhysics){
 		throw std::logic_error("Bodies in interaction cannot be swapped if they have interactionGeometry or interactionPhysics.");

Modified: trunk/core/Interaction.hpp
===================================================================
--- trunk/core/Interaction.hpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/core/Interaction.hpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -63,6 +63,8 @@
 
 		//! Reset interaction to the intial state (keep only body ids)
 		void reset();
+		//! common initialization called from both constructor and reset()
+		void init();
 			
 	REGISTER_ATTRIBUTES(/*no base*/,
 		(id1)

Modified: trunk/lib/smoothing/WeightedAverage2d.cpp
===================================================================
--- trunk/lib/smoothing/WeightedAverage2d.cpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/lib/smoothing/WeightedAverage2d.cpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -22,6 +22,7 @@
 		.add_property("stDev",&pyGaussAverage::stDev_get,&pyGaussAverage::stDev_set)
 		.add_property("relThreshold",&pyGaussAverage::relThreshold_get,&pyGaussAverage::relThreshold_set)
 		.add_property("clips",&pyGaussAverage::clips_get,&pyGaussAverage::clips_set)
+		.add_property("data",&pyGaussAverage::data_get)
 	;
 };
 

Modified: trunk/lib/smoothing/WeightedAverage2d.hpp
===================================================================
--- trunk/lib/smoothing/WeightedAverage2d.hpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/lib/smoothing/WeightedAverage2d.hpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -47,6 +47,7 @@
 		return ret;
 	}
 	Vector2r cell2xyMid(Vector2i cxy) const { return Vector2r(lo[0]+cellSizes[0]*(.5+cxy[0]),lo[1]+cellSizes[1]*(.5+cxy[1])); }
+	const Vector2i& getSize() const{ return nCells;}
 
 	// add new element to the right cell
 	void add(const T& t, Vector2r xy){Vector2i cxy=xy2cell(xy); grid[cxy[0]][cxy[1]].push_back(t);}
@@ -196,5 +197,17 @@
 			clips.push_back(poly);
 		}
 	}
+	python::tuple data_get(){
+		python::list x,y,val;
+		const Vector2i& dim=sgka->grid->getSize();
+		for(int i=0; i<dim[0]; i++){
+			for(int j=0; j<dim[1]; j++){
+				FOREACH(const Scalar2d& element, sgka->grid->grid[i][j]){
+					x.append(element.pos[0]); y.append(element.pos[1]); val.append(element.val);
+				}
+			}
+		}
+		return python::make_tuple(x,y,val);
+	}
 };
 

Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/pkg/dem/ConcretePM.cpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -49,7 +49,7 @@
 	contPhys->ks=contPhys->G*contPhys->crossSection;
 
 	if(neverDamage) contPhys->neverDamage=true;
-	if(cohesiveThresholdIter<0 || Omega::instance().getCurrentIteration()<cohesiveThresholdIter) contPhys->isCohesive=true;
+	if(cohesiveThresholdIter<0 || (Omega::instance().getCurrentIteration()<cohesiveThresholdIter)) contPhys->isCohesive=true;
 	else contPhys->isCohesive=false;
 	contPhys->dmgTau=dmgTau;
 	contPhys->dmgRateExp=dmgRateExp;
@@ -192,13 +192,12 @@
 
 	// handle broken contacts
 	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.");
+			rbp1->numBrokenCohesive+=1; rbp2->numBrokenCohesive+=1; rbp1->epsPlBroken+=epsPlSum; rbp2->epsPlBroken+=epsPlSum;
 		}
+		rootBody->interactions->requestErase(I->getId1(),I->getId2());
 		return;
 	}
 
@@ -333,6 +332,7 @@
 
 
 /********************** CpmPhysDamageColorizer ****************************/
+CREATE_LOGGER(CpmPhysDamageColorizer);
 void CpmPhysDamageColorizer::action(MetaBody* rootBody){
 	//vector<pair<short,Real> > bodyDamage; /* number of cohesive interactions per body; cummulative damage of interactions */
 	//vector<pair<short,
@@ -344,19 +344,20 @@
 		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();
+		const 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;
+		int 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;
+			if(bpp->normDmg>1){
+				LOG_WARN("#"<<id<<" normDmg="<<bpp->normDmg<<" nCohLinks="<<bodyStats[id].nCohLinks<<", numBrokenCohesive="<<bpp->numBrokenCohesive<<", dmgSum="<<bodyStats[id].dmgSum<<", numAllCohLinks"<<cohLinksWhenever);
+			}
 		}
 		else { bpp->normDmg=0; bpp->normEpsPl=0;}
 		B->geometricalModel->diffuseColor=Vector3r(bpp->normDmg,1-bpp->normDmg,B->isDynamic?0:1);

Modified: trunk/pkg/dem/ConcretePM.hpp
===================================================================
--- trunk/pkg/dem/ConcretePM.hpp	2009-06-15 15:39:57 UTC (rev 1798)
+++ trunk/pkg/dem/ConcretePM.hpp	2009-06-15 17:34:09 UTC (rev 1799)
@@ -306,12 +306,13 @@
 #endif
 
 class CpmPhysDamageColorizer: public PeriodicEngine {
-	struct BodyStats{ short nCohLinks; Real dmgSum; Real epsPlSum; BodyStats(): nCohLinks(0), dmgSum(0), epsPlSum(0.){} };
+	struct BodyStats{ int 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*);
+	DECLARE_LOGGER;
 	REGISTER_ATTRIBUTES(PeriodicEngine,(maxOmega));
 	REGISTER_CLASS_AND_BASE(CpmPhysDamageColorizer,PeriodicEngine);
 };