yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-07-08 11:24:58 +0200 (Wed, 08 Jul 2009)
New Revision: 1848

Added:
   trunk/examples/triax-perf/mkTextTable.sh
   trunk/examples/triax-perf/triax-perf.ods
Modified:
   trunk/examples/triax-perf/triax-perf.py
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
   trunk/pkg/dem/ConcretePM.cpp
   trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
   trunk/py/pack.py
Log:
1. Dump XML in ConcretePM in case of problems (nan's)
2. InteractingFacet::postProcessAttributes gives FATAL message if there are coincident vertices
3. Update triax-perf, add script and openoffice sheet to make table
4. InsertionSortCollider is used in TriaxialTest always (not only with fast) now.


Added: trunk/examples/triax-perf/mkTextTable.sh
===================================================================
--- trunk/examples/triax-perf/mkTextTable.sh	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/examples/triax-perf/mkTextTable.sh	2009-07-08 09:24:58 UTC (rev 1848)
@@ -0,0 +1,10 @@
+#!/bin/sh
+# run this to get nice text to paste into the openoffice sheet
+TMPS=/tmp/a$$ /tmp/b$$ /tmp/c$$ /tmp/d$$
+ls *.ser?.log
+awk '/TOTAL/ { print $2}' *.ser?.log > /tmp/a$$
+awk '/Collider/ { print $4}' *.ser?.log > /tmp/b$$
+awk '/TOTAL/ { print $2}' *.par?.log > /tmp/c$$
+awk '/Collider/ { print $4}' *.par?.log > /tmp/d$$
+paste /tmp/a$$ /tmp/b$$ /tmp/c$$ /tmp/d$$ | sed 's/us//g' | sed 's/%//g' | sed 's/\./,/g'
+rm -rf $TMPS

Added: trunk/examples/triax-perf/triax-perf.ods
===================================================================
(Binary files differ)


Property changes on: trunk/examples/triax-perf/triax-perf.ods
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Modified: trunk/examples/triax-perf/triax-perf.py
===================================================================
--- trunk/examples/triax-perf/triax-perf.py	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/examples/triax-perf/triax-perf.py	2009-07-08 09:24:58 UTC (rev 1848)
@@ -10,13 +10,14 @@
 # The -j1 ensures that only 1 job will run at time
 # (even if other cores are free, access to memory is limiting if running multiple jobs at time)
 #
-# You have to collect the results by hand from log files.
+# You have to collect the results by hand from log files, or run sh mkTextTable.sh and use
+# triax-perf.ods to get comparison
 #
 utils.readParamsFromTable(fast=False,noTableOk=True)
-p=Preprocessor('TriaxialTest',{'numberOfGrains':8000,'fast':fast}).load()
+TriaxialTest(numberOfGrains=50000,fast=fast,noFiles=True).load()
 O.run(10,True) # filter out initialization
 O.timingEnabled=True
-O.run(1000,True)
+O.run(200,True)
 from yade import timing
 timing.stats()
 print 'BexContainer synced %d times'%(O.bexSyncCount)

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp	2009-07-08 09:24:58 UTC (rev 1848)
@@ -7,6 +7,8 @@
 *************************************************************************/
 #include "InteractingFacet.hpp"
 
+CREATE_LOGGER(InteractingFacet);
+
 InteractingFacet::InteractingFacet() : InteractingGeometry()
 {
 	createIndex();
@@ -35,6 +37,9 @@
     if (deserializing)
     {
 		Vector3r e[3] = {vertices[1]-vertices[0] ,vertices[2]-vertices[1] ,vertices[0]-vertices[2]};
+		#define CHECK_EDGE(i) if(e[i].SquaredLength()==0){LOG_FATAL("InteractingFacet has coincident vertices "<<i<<" ("<<vertices[i]<<") and "<<(i+1)%3<<" ("<<vertices[(i+1)%3]<<")!");}
+			CHECK_EDGE(0); CHECK_EDGE(1);CHECK_EDGE(2);
+		#undef CHECK_EDGE
 		nf = e[0].UnitCross(e[1]);
 		for(int i=0; i<3; ++i) 
 		{

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp	2009-07-08 09:24:58 UTC (rev 1848)
@@ -45,6 +45,8 @@
 		vector<Real> edgeAdjHalfAngle;
 	#endif
 
+	DECLARE_LOGGER;
+
 	protected:
 
 	void registerAttributes(); void postProcessAttributes(bool deserializing);

Modified: trunk/pkg/dem/ConcretePM.cpp
===================================================================
--- trunk/pkg/dem/ConcretePM.cpp	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/pkg/dem/ConcretePM.cpp	2009-07-08 09:24:58 UTC (rev 1848)
@@ -147,7 +147,8 @@
 	epsN=contGeom->strainN(); epsT=contGeom->strainT();
 	
 	// debugging
-		#define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); throw;}
+		#define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); LOG_FATAL("in interaction #"<<I->getId1()<<"+#"<<I->getId2()); Omega::instance().saveSimulation("/tmp/verificationFailed.xml"); throw;}
+
 		#define NNAN(a) YADE_VERIFY(!isnan(a));
 		#define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
 		#ifdef YADE_DEBUG

Modified: trunk/pkg/dem/PreProcessor/TriaxialTest.cpp
===================================================================
--- trunk/pkg/dem/PreProcessor/TriaxialTest.cpp	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/pkg/dem/PreProcessor/TriaxialTest.cpp	2009-07-08 09:24:58 UTC (rev 1848)
@@ -630,8 +630,7 @@
 	rootBody->engines.clear();
 	rootBody->engines.push_back(shared_ptr<Engine>(new PhysicalActionContainerReseter));
 	rootBody->engines.push_back(boundingVolumeDispatcher);
-	if(!fast) rootBody->engines.push_back(shared_ptr<Engine>(new PersistentSAPCollider));
-	else rootBody->engines.push_back(shared_ptr<Engine>(new InsertionSortCollider));
+	rootBody->engines.push_back(shared_ptr<Engine>(new InsertionSortCollider));
 	if(fast){
 		shared_ptr<InteractionDispatchers> ids(new InteractionDispatchers);
 			ids->geomDispatcher=interactionGeometryDispatcher;

Modified: trunk/py/pack.py
===================================================================
--- trunk/py/pack.py	2009-07-08 08:41:42 UTC (rev 1847)
+++ trunk/py/pack.py	2009-07-08 09:24:58 UTC (rev 1848)
@@ -15,6 +15,8 @@
 # import SpherePack
 from _packSpheres import *
 
+from miniWm3Wrap import *
+
 class inGtsSurface_py(Predicate):
 	"""This class was re-implemented in c++, but should stay here to serve as reference for implementing
 	Predicates in pure python code. C++ allows us to play dirty tricks in GTS which are not accessible
@@ -73,9 +75,17 @@
 
 def sweptPolylines2gtsSurface(pts,threshold=0,capStart=False,capEnd=False):
 	"""Create swept suface (as GTS triangulation) given same-length sequences of points (as 3-tuples).
-	If threshold is given (>0), gts.Surface().cleanup(threshold) will be called before returning, which
-	removes vertices mutually closer than threshold. Can be used to create closed swept surface (revolved), as
-	we don't check for coincident vertices otherwise.
+
+	If threshold is given (>0), then
+
+	* degenerate faces (with edges shorter than threshold) will not be created
+	* gts.Surface().cleanup(threshold) will be called before returning, which merges vertices mutually
+		closer than threshold. In case your pts are closed (last point concident with the first one)
+		this will the surface strip of triangles. If you additionally have capStart==True and capEnd==True,
+		the surface will be closed.
+
+	Note: capStart and capEnd make the most naive polygon triangulation (diagonals) and will perhaps fail
+	for non-convex sections.
 	"""
 	if not len(set([len(pts1) for pts1 in pts]))==1: raise RuntimeError("Polylines must be all of the same length!")
 	vtxs=[[gts.Vertex(x,y,z) for x,y,z in pts1] for pts1 in pts]
@@ -85,12 +95,19 @@
 		for j in range(0,len(vtxs[i])):
 			interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j]))
 			if j<len(vtxs[i])-1: interSectEdges[i].append(gts.Edge(vtxs[i][j],vtxs[i+1][j+1]))
+	if threshold>0: # replace edges of zero length with None; their faces will be skipped
+		def fixEdges(edges):
+			for i,e in enumerate(edges):
+				if (Vector3(e.v1.x,e.v1.y,e.v1.z)-Vector3(e.v2.x,e.v2.y,e.v2.z)).Length()<threshold: edges[i]=None
+		for ee in sectEdges: fixEdges(ee)
+		for ee in interSectEdges: fixEdges(ee)
 	surf=gts.Surface()
 	for i in range(0,len(vtxs)-1):
 		for j in range(0,len(vtxs[i])-1):
-			newFaces=gts.Face(interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]),gts.Face(sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1])
-			for face in newFaces:
-				if face.is_ok: surf.add(face)
+			ee1=interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]
+			ee2=sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1]
+			if None not in ee1: surf.add(gts.Face(interSectEdges[i][2*j+1],sectEdges[i+1][j],interSectEdges[i][2*j]))
+			if None not in ee2: surf.add(gts.Face(sectEdges[i][j],interSectEdges[i][2*j+2],interSectEdges[i][2*j+1]))
 	def doCap(vtx,edg,start):
 		ret=[]
 		eFan=[edg[0]]+[gts.Edge(vtx[i],vtx[0]) for i in range(2,len(vtx))]