yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-04-10 13:36:17 +0200 (Fri, 10 Apr 2009)
New Revision: 1754

Modified:
   trunk/extra/Brefcom.cpp
   trunk/extra/Brefcom.hpp
   trunk/gui/py/eudoxos.py
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
   trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
   trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
   trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
   trunk/scripts/test/facet-topo.py
Log:
1. Fix bug in FacetTopologyAnalyzer algorithm
2. Add angle info to InteractingFacet (not yet used)
3. Add triangulated sphere test for FacetTopoloyAnalyzer to facet-topo.py
4. Fixes in rate-dep in brefcom


Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/extra/Brefcom.cpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -172,8 +172,9 @@
 
 Real BrefcomContact::computeViscoplScalingFactor(Real sigmaTNorm, Real sigmaTYield,Real dt){
 	if(sigmaTNorm<sigmaTYield) return 1.;
-	Real c=sigmaTNorm*pow(plTau/(G*dt),plRateExp)*pow(sigmaTNorm-sigmaTYield,plRateExp-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);
 }
 

Modified: trunk/extra/Brefcom.hpp
===================================================================
--- trunk/extra/Brefcom.hpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/extra/Brefcom.hpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -13,6 +13,7 @@
 #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.

Modified: trunk/gui/py/eudoxos.py
===================================================================
--- trunk/gui/py/eudoxos.py	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/gui/py/eudoxos.py	2009-04-10 11:36:17 UTC (rev 1754)
@@ -174,7 +174,7 @@
 	f.write("SimpleCS 1 thick 1.0 width 1.0\n")
 	# material
 	ph=inters[0].phys
-	f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g damchartime %g damrateexp %g d 1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle'],ph['dmgTau'],ph['dmgRateExp']))
+	f.write("CohInt 1 kn %g ks %g e0 %g ef %g c 0. ksi %g coh %g tanphi %g damchartime %g damrateexp %g plchartime %g plrateexp %g d 1.0\n"%(ph['E'],ph['G'],ph['epsCrackOnset'],ph['epsFracture'],ph['xiShear'],ph['undamagedCohesion'],ph['tanFrictionAngle'],ph['dmgTau'],ph['dmgRateExp'],ph['plTau'],ph['plRateExp']))
 	# boundary conditions
 	f.write('BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0\n')
 	f.write('BoundaryCondition 2 loadTimeFunction 1 prescribedvalue 1.e-4\n')

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -12,6 +12,7 @@
 	createIndex();
 	#ifdef FACET_TOPO
 		edgeAdjIds.resize(3,Body::ID_NONE);	
+		edgeAdjAngle.resize(3,0);
 	#endif
 }
 
@@ -25,6 +26,7 @@
     REGISTER_ATTRIBUTE(vertices);
 	#ifdef FACET_TOPO
 		REGISTER_ATTRIBUTE(edgeAdjIds);
+		REGISTER_ATTRIBUTE(edgeAdjAngle);
 	#endif
 }
 

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -41,6 +41,8 @@
 	#ifdef FACET_TOPO
 		//! facet id's that are adjacent to respective edges
 		vector<body_id_t> edgeAdjIds;
+		//! angle between normals of this facet and the adjacent facet
+		vector<body_id_t> edgeAdjAngle;
 	#endif
 
 	protected:

Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -251,8 +251,6 @@
 
 /* Note that this function is called only for bodies that actually overlap along some axis */
 void PersistentSAPCollider::updateOverlapingBBSet(int id1,int id2){
-	#pragma omp critical
-	{
 		// look if the pair (id1,id2) already exists in the overlappingBB collection
 		const shared_ptr<Interaction>& interaction=transientInteractions->find(body_id_t(id1),body_id_t(id2));
 		bool found=(bool)interaction;
@@ -279,7 +277,6 @@
 			//LOG_DEBUG("Erasing interaction #"<<id1<<"=#"<<id2<<" (isReal="<<interaction->isReal<<")");
 			transientInteractions->erase(body_id_t(id1),body_id_t(id2));
 		}
-	}
 }
 
 
@@ -322,7 +319,14 @@
 		while (i<2*nbElements && !bounds[i]->lower) i++;
 		j=i+1;
 		while (j<2*nbElements && bounds[j]->id!=bounds[i]->id){
-			if (bounds[j]->lower) updateOverlapingBBSet(bounds[i]->id,bounds[j]->id);
+			if (bounds[j]->lower){
+				/* findOverlappingBB is called in parallel for different data at initialization,
+				 * but updateOverlapingBBSet touches shared global data, therefore must be protected by critical section;
+				 * normally updateOverlapingBBSet is also called sequentially from sortBounds, where the critical section
+				 * would penalize performance; therefore we have to protect it from here */
+				#pragma omp critical 
+					updateOverlapingBBSet(bounds[i]->id,bounds[j]->id);
+			}
 			j++;
 		}
 		i++;

Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -61,6 +61,7 @@
 		if(v->vertexId>=0) continue; // already assigned
 	}
 	LOG_DEBUG("Found "<<maxVertexId<<" unique vertices.");
+	commonVerticesFound=maxVertexId;
 	// add FacetTopology for all facets; index within the topo array is the body id
 	vector<shared_ptr<FacetTopology> > topo(rb->bodies->size()); // initialized with the default ctor
 	FOREACH(shared_ptr<VertexData>& v, vv){
@@ -80,21 +81,18 @@
 		}
 		topo1.push_back(t);
 	}
-	sort(topo1.begin(),topo1.end(),TopologyIndexComparator());
 	size_t nTopo=topo1.size();
 	for(size_t i=0; i<nTopo; i++){
 		size_t j=i;
 		while(++j<nTopo){
 			const shared_ptr<FacetTopology>& ti(topo1[i]), &tj(topo1[j]);
-			LOG_TRACE("Analyzing possibly-adjacent facets #"<<ti->id<<" #"<<tj->id<<" (vertices "<<ti->vertices[0]<<","<<ti->vertices[1]<<","<<ti->vertices[2]<<"; "<<tj->vertices[0]<<","<<ti->vertices[1]<<","<<ti->vertices[2]<<")");
 			vector<size_t> vvv; // array of common vertices
 			for(size_t k=0; k<3; k++){
 				if     (ti->vertices[k]==tj->vertices[0]) vvv.push_back(ti->vertices[k]);
 				else if(ti->vertices[k]==tj->vertices[1]) vvv.push_back(ti->vertices[k]);
 				else if(ti->vertices[k]==tj->vertices[2]) vvv.push_back(ti->vertices[k]);
 			}
-			if(vvv.size()==0) break; // reached end of those that have the lowest-id vertex in common
-			if(vvv.size()==1) continue; // only one edge in common
+			if(vvv.size()<2) continue;
 			assert(vvv.size()!=3); // same coords? nonsense
 			assert(vvv.size()==2);
 			vector<int> edge(2,0);
@@ -112,7 +110,7 @@
 			}
 			// add adjacency information to the facet itself
 			YADE_PTR_CAST<InteractingFacet>((*rb->bodies)[ti->id]->interactingGeometry)->edgeAdjIds[edge[0]]=tj->id;
-			YADE_PTR_CAST<InteractingFacet>((*rb->bodies)[tj->id]->interactingGeometry)->edgeAdjIds[edge[1]]=tj->id;
+			YADE_PTR_CAST<InteractingFacet>((*rb->bodies)[tj->id]->interactingGeometry)->edgeAdjIds[edge[1]]=ti->id;
 			commonEdgesFound++;
 			LOG_TRACE("Added adjacency information for #"<<ti->id<<"+#"<<tj->id<<" (common edges "<<edge[0]<<"+"<<edge[1]<<")");
 		}

Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp	2009-04-10 11:36:17 UTC (rev 1754)
@@ -40,9 +40,6 @@
 		//! facet id, for back reference
 		body_id_t id;
 	};
-	struct TopologyIndexComparator{
-		bool operator()(const shared_ptr<FacetTopology>& t1, const shared_ptr<FacetTopology>& t2){ return min(t1->vertices[0],min(t1->vertices[1],t1->vertices[2]))<min(t2->vertices[0],min(t2->vertices[1],t2->vertices[2])); }
-	};
 	public:
 		//! Axis along which to do the initial vertex sort
 		Vector3r projectionAxis;
@@ -50,11 +47,13 @@
 		Real relTolerance;
 		//! how many common edges were identified during last run
 		long commonEdgesFound;
+		//! how many common vertices were identified during last run
+		long commonVerticesFound;
 	void action(MetaBody*); 
-	FacetTopologyAnalyzer(): projectionAxis(Vector3r::UNIT_X), relTolerance(1e-4), commonEdgesFound(0) {}
+	FacetTopologyAnalyzer(): projectionAxis(Vector3r::UNIT_X), relTolerance(1e-4), commonEdgesFound(0), commonVerticesFound(0) {}
 	DECLARE_LOGGER;
 	REGISTER_CLASS_AND_BASE(FacetTopologyAnalyzer,StandAloneEngine);
-	REGISTER_ATTRIBUTES(StandAloneEngine, (projectionAxis)(relTolerance));
+	REGISTER_ATTRIBUTES(StandAloneEngine, (projectionAxis)(relTolerance)(commonEdgesFound)(commonVerticesFound));
 };
 REGISTER_SERIALIZABLE(FacetTopologyAnalyzer);
 

Modified: trunk/scripts/test/facet-topo.py
===================================================================
--- trunk/scripts/test/facet-topo.py	2009-04-10 08:49:36 UTC (rev 1753)
+++ trunk/scripts/test/facet-topo.py	2009-04-10 11:36:17 UTC (rev 1754)
@@ -1,9 +1,9 @@
 import yade.log
-yade.log.setLevel('FacetTopologyAnalyzer',yade.log.TRACE)
+#yade.log.setLevel('FacetTopologyAnalyzer',yade.log.TRACE)
 
 # Note: FacetTopologyAnalyzer is normally run as an initializer;
 # it is only for testing sake that it is in O.engines here.
-O.engines=[FacetTopologyAnalyzer(projectionAxis=(0,0,1),label='topo'),]
+O.engines=[FacetTopologyAnalyzer(projectionAxis=(1,1,1),label='topo'),]
 
 # most simple case: no touch at all
 if 1:
@@ -22,3 +22,26 @@
 	O.step()
 	assert(O.bodies[0].phys['edgeAdjIds'][1]==1 and O.bodies[1].phys['edgeAdjIds'][0]==1)
 	assert(topo['commonEdgesFound']==1)
+if 1:
+	O.bodies.clear()
+	r=.5 # radius of the sphere
+	nPoly=12; # try 128, it is still quite fast
+	def sphPt(i,j):
+		if i==0: return (0,0,-r)
+		if i==nPoly/2: return (0,0,r)
+		assert(i>0 and i<nPoly/2)
+		assert(j>=0 and j<=nPoly)
+		theta=i*2*pi/nPoly
+		rr=r*sin(theta)
+		phi=j*2*pi/nPoly
+		return rr*cos(phi),rr*sin(phi),-r*cos(theta)
+	for i in range(0,nPoly/2):
+		for j in range(0,nPoly):
+			if i!=0: O.bodies.append(utils.facet([sphPt(i,j),sphPt(i,j+1),sphPt(i+1,j)]))
+			if i!=nPoly/2-1: O.bodies.append(utils.facet([sphPt(i+1,j),sphPt(i,j+1),sphPt(i+1,j+1)]))
+	print 'Sphere created, has',len(O.bodies),'facets'
+	O.step()
+	assert(topo['commonVerticesFound']==nPoly*(nPoly/2-1)+2)
+	assert(topo['commonEdgesFound']==nPoly*((nPoly/2-1)+(nPoly/2-2)*2+2))
+	print topo['commonVerticesFound'],'vertices; ',topo['commonEdgesFound'],'edges'
+#quit()