Author: eudoxos
Date: 2009-04-08 22:27:04 +0200 (Wed, 08 Apr 2009)
New Revision: 1750

Added:
   trunk/scripts/test/facet-topo.py
Modified:
   trunk/lib/QGLViewer/VRender/BSPSortMethod.cpp
   trunk/lib/QGLViewer/VRender/Exporter.cpp
   trunk/lib/QGLViewer/VRender/FIGExporter.cpp
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
   trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
   trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
   trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp
   trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
   trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
Log:
1. Parallelize initial bound filling in PersistentSAPCollider. Gives almost 3x 
speedup for the first step.
2. Fix missing headers so that it compiles with g++-4.4 (in QGLViewer)
3. FacetTopologyAnalyzer works (not tested extensively), questions of where to 
put topology data for InteractingFacet (will be raised on yade-dev)
4. Test script for FacetTopologyAnalyzer.



Modified: trunk/lib/QGLViewer/VRender/BSPSortMethod.cpp
===================================================================
--- trunk/lib/QGLViewer/VRender/BSPSortMethod.cpp       2009-04-06 21:40:30 UTC 
(rev 1749)
+++ trunk/lib/QGLViewer/VRender/BSPSortMethod.cpp       2009-04-08 20:27:04 UTC 
(rev 1750)
@@ -49,6 +49,7 @@
 #include "Primitive.h"
 #include "SortMethod.h"
 #include "math.h" // fabs
+#include <cstdio>
 
 using namespace vrender ;
 using namespace std;

Modified: trunk/lib/QGLViewer/VRender/Exporter.cpp
===================================================================
--- trunk/lib/QGLViewer/VRender/Exporter.cpp    2009-04-06 21:40:30 UTC (rev 
1749)
+++ trunk/lib/QGLViewer/VRender/Exporter.cpp    2009-04-08 20:27:04 UTC (rev 
1750)
@@ -48,6 +48,7 @@
 #include <stdexcept>
 #include "VRender.h"
 #include "Exporter.h"
+#include <cstdio>
 
 using namespace vrender ;
 using namespace std ;

Modified: trunk/lib/QGLViewer/VRender/FIGExporter.cpp
===================================================================
--- trunk/lib/QGLViewer/VRender/FIGExporter.cpp 2009-04-06 21:40:30 UTC (rev 
1749)
+++ trunk/lib/QGLViewer/VRender/FIGExporter.cpp 2009-04-08 20:27:04 UTC (rev 
1750)
@@ -47,6 +47,7 @@
 
 #include "Exporter.h"
 #include "math.h"
+#include <cstdio>
 
 using namespace vrender ;
 using namespace std ;

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp 
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.cpp 
2009-04-08 20:27:04 UTC (rev 1750)
@@ -9,7 +9,10 @@
 
 InteractingFacet::InteractingFacet() : InteractingGeometry()
 {
-    createIndex();
+       createIndex();
+       #ifdef FACET_TOPO
+               edgeAdjIds.resize(3,Body::ID_NONE);     
+       #endif
 }
 
 InteractingFacet::~InteractingFacet()
@@ -20,6 +23,9 @@
 {
     InteractingGeometry::registerAttributes();
     REGISTER_ATTRIBUTE(vertices);
+       #ifdef FACET_TOPO
+               REGISTER_ATTRIBUTE(edgeAdjIds);
+       #endif
 }
 
 void InteractingFacet::postProcessAttributes(bool deserializing)

Modified: trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp
===================================================================
--- trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/common/DataClass/InteractingGeometry/InteractingFacet.hpp 
2009-04-08 20:27:04 UTC (rev 1750)
@@ -9,7 +9,11 @@
 
 
 #include <yade/core/InteractingGeometry.hpp>
+#include<yade/core/Body.hpp>
 
+// define this to have topology information about facets enabled;
+// it is necessary for FacetTopologyAnalyzer
+#define FACET_TOPO
 
 class InteractingFacet : public InteractingGeometry {
     public:
@@ -34,6 +38,10 @@
        Real vl[3];
        /// Unit vertice vectors
        Vector3r vu[3];
+       #ifdef FACET_TOPO
+               //! facet id's that are adjacent to respective edges
+               vector<body_id_t> edgeAdjIds;
+       #endif
 
        protected:
 

Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp  
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.cpp  
2009-04-08 20:27:04 UTC (rev 1750)
@@ -19,6 +19,7 @@
 PersistentSAPCollider::PersistentSAPCollider() : Collider()
 {
        haveDistantTransient=false;
+       ompBodiesMin=0;
 
        nbObjects=0;
        xBounds.clear();
@@ -27,7 +28,7 @@
        minima.clear();
        maxima.clear();
 
-       //timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
+//     timingDeltas=shared_ptr<TimingDeltas>(new TimingDeltas);
 }
 
 PersistentSAPCollider::~PersistentSAPCollider()
@@ -108,15 +109,22 @@
        nbObjects=bodies->size();
 
        // permutation sort of the AABBBounds along the 3 axis performed in a 
independant manner
-       //#pragma omp parallel sections
-       {
-       //#pragma omp section
-               sortBounds(xBounds, nbObjects);
-       //#pragma omp section
-               sortBounds(yBounds, nbObjects);
-       //#pragma omp section
-               sortBounds(zBounds, nbObjects);
-       }
+       // serial version
+       //if(nbObjects>ompBodiesMin || ompBodiesMin==0){ … }
+       sortBounds(xBounds,nbObjects); sortBounds(yBounds,nbObjects); 
sortBounds(zBounds,nbObjects);
+       #if 0
+               else {
+                       #pragma omp parallel sections
+                       {
+                       #pragma omp section
+                               sortBounds(xBounds, nbObjects);
+                       #pragma omp section
+                               sortBounds(yBounds, nbObjects);
+                       #pragma omp section
+                               sortBounds(zBounds, nbObjects);
+                       }
+               }
+       #endif
 
 //     timingDeltas->checkpoint("sortBounds");
 }
@@ -167,26 +175,60 @@
                // initialization if the field "value" of the xBounds, yBounds, 
zBounds arrays
                updateBounds(nbElements);
 
-               // modified quick sort of the xBounds, yBounds, zBounds arrays
+               /* Performance note: such was the timing result on initial step 
of 8k sphere in triaxial test.
+                       the findX, findY, findZ take almost the totality of the 
time.
+                       Parallelizing those is vastly beneficial (almost 3x 
speed increase, which can be quite sensible as the initial
+                       findOverlappingBB is really slow 
http://yade.wikia.com/wiki/Colliders_performace and is done in 3
+                       orthogonal directions. Therefore, it is enabled by 
default.
+                       
+                       Now sortX is right before findX etc, in the same openMP 
section. Beware that timingDeltas will give garbage
+                       results if used in parallelized code.
+
+                       ===
+
+                       8k spheres:
+                       Name                                                    
Count                 Time            Rel. time
+                       
-------------------------------------------------------------------------------------------------------
+                       PersistentSAPCollider                                 1 
           3568091us              100.00%      
+                         init                                                  
1              21178us                0.59%    
+                         sortX                                                 
1              33225us                0.93%    
+                         sortY                                                 
1              29300us                0.82%    
+                         sortZ                                                 
1              28334us                0.79%    
+                         findX                                                 
1            1708426us               47.88%    
+                         findY                                                 
1             869150us               24.36%    
+                         findZ                                                 
1             867378us               24.31%    
+                         TOTAL                                                 
             3556994us               99.69%    
+
+               */
+
                // The first time these arrays are not sorted so it is faster 
to use such a sort instead
                // of the permutation we are going to use next
-               
std::sort(xBounds.begin(),xBounds.begin()+2*nbElements,AABBBoundComparator());
-               //timingDeltas->checkpoint("sortX");
-               
std::sort(yBounds.begin(),yBounds.begin()+2*nbElements,AABBBoundComparator());
-               //timingDeltas->checkpoint("sortY");
-               
std::sort(zBounds.begin(),zBounds.begin()+2*nbElements,AABBBoundComparator());
-               //timingDeltas->checkpoint("sortZ");
 
-               // initialize the overlappingBB collection
-               //for(unsigned int j=0;j<nbElements;j++)
-               //      overlappingBB[j].clear(); //attention memoire
-
-               findOverlappingBB(xBounds, nbElements);
-               //timingDeltas->checkpoint("findX");
-               findOverlappingBB(yBounds, nbElements);
-               //timingDeltas->checkpoint("findY");
-               findOverlappingBB(zBounds, nbElements);
-               //timingDeltas->checkpoint("findZ");
+               // do not juse timingDeltas with openMP enabled, results will 
be garbage
+               #pragma omp parallel sections
+               {
+                       #pragma omp section
+                       {
+                               
std::sort(xBounds.begin(),xBounds.begin()+2*nbElements,AABBBoundComparator());
+                               //timingDeltas->checkpoint("sortX");
+                               findOverlappingBB(xBounds, nbElements);
+                               //timingDeltas->checkpoint("findX");
+                       }
+                       #pragma omp section
+                       {
+                               
std::sort(yBounds.begin(),yBounds.begin()+2*nbElements,AABBBoundComparator());
+                               //timingDeltas->checkpoint("sortY");
+                               findOverlappingBB(yBounds, nbElements);
+                               //timingDeltas->checkpoint("findY");
+                       }
+                       #pragma omp section
+                       {
+                               
std::sort(zBounds.begin(),zBounds.begin()+2*nbElements,AABBBoundComparator());
+                               //timingDeltas->checkpoint("sortZ");
+                               findOverlappingBB(zBounds, nbElements);
+                               //timingDeltas->checkpoint("findZ");
+                       }
+               }
        }
        else updateBounds(nbElements);
 }
@@ -240,16 +282,34 @@
 
 void PersistentSAPCollider::updateBounds(int nbElements)
 {
-       for(int i=0; i < 2*nbElements; i++){
-               if (xBounds[i]->lower) xBounds[i]->value = 
minima[3*xBounds[i]->id+0];
-               else xBounds[i]->value = maxima[3*xBounds[i]->id+0];
-
-               if (yBounds[i]->lower) yBounds[i]->value = 
minima[3*yBounds[i]->id+1];
-               else yBounds[i]->value = maxima[3*yBounds[i]->id+1];
-
-               if (zBounds[i]->lower) zBounds[i]->value = 
minima[3*zBounds[i]->id+2];
-               else zBounds[i]->value = maxima[3*zBounds[i]->id+2];
+       #define _BOUND_UPDATE(bounds,offset) \
+               if (bounds[i]->lower) bounds[i]->value = 
minima[3*bounds[i]->id+offset]; \
+               else bounds[i]->value = maxima[3*bounds[i]->id+offset];
+       // for small number of bodies, run sequentially
+       #if 0
+       if(nbElements<ompBodiesMin || ompBodiesMin==0){
+       #endif
+               for(int i=0; i < 2*nbElements; i++){
+                       _BOUND_UPDATE(xBounds,0);
+                       _BOUND_UPDATE(yBounds,1);
+                       _BOUND_UPDATE(zBounds,2);
+               }
+       #if 0
        }
+       else{
+               // parallelize for large number of bodies (not used, 
updateBounds takes only about 5% of collider time
+               #pragma omp parallel sections
+               {
+                       #pragma omp section
+                       for(int i=0; i < 2*nbElements; i++){ 
_BOUND_UPDATE(xBounds,0); }
+                       #pragma omp section
+                       for(int i=0; i < 2*nbElements; i++){ 
_BOUND_UPDATE(yBounds,1); }
+                       #pragma omp section
+                       for(int i=0; i < 2*nbElements; i++){ 
_BOUND_UPDATE(zBounds,2); }
+               }
+       }
+       #endif
+       #undef _BOUND_UPDATE
 }
 
 

Modified: trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp  
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/common/Engine/StandAloneEngine/PersistentSAPCollider.hpp  
2009-04-08 20:27:04 UTC (rev 1750)
@@ -99,9 +99,13 @@
                //! Don't break transient interaction once bodies don't overlap 
anymore; material law will be responsible for breaking it.
                bool haveDistantTransient;
 
+               //! minimum number of bodies to run updateIds in parallel 
secions; if 0 (default for now), never run in parallel
+               long ompBodiesMin;
+
                void registerAttributes(){
                        Collider::registerAttributes();
                        REGISTER_ATTRIBUTE(haveDistantTransient);
+                       REGISTER_ATTRIBUTE(ompBodiesMin);
                }
 
        DECLARE_LOGGER;

Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp     
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.cpp     
2009-04-08 20:27:04 UTC (rev 1750)
@@ -2,11 +2,21 @@
 #include<yade/pkg-common/InteractingFacet.hpp>
 #include<yade/core/MetaBody.hpp>
 #include<yade/core/Body.hpp>
+
+CREATE_LOGGER(FacetTopologyAnalyzer);
+YADE_PLUGIN("FacetTopologyAnalyzer");
+
+#ifndef FACET_TOPO
 void FacetTopologyAnalyzer::action(MetaBody* rb){
+       throw runtime_error("FACET_TOPO was not enabled in InteractingFacet.hpp 
at compile-time. Do not use FacetTopologyAnalyzer or recompile.");
+}
+#else
+void FacetTopologyAnalyzer::action(MetaBody* rb){
+       commonEdgesFound=0;
        LOG_DEBUG("Projection axis for analysis is "<<projectionAxis);
        vector<shared_ptr<VertexData> > vv;
        // minimum facet edge length (tolerance scale)
-       Real minSqLen=0;
+       Real minSqLen=numeric_limits<Real>::infinity();
        FOREACH(const shared_ptr<Body>& b, *rb->bodies){
                shared_ptr<InteractingFacet> 
f=dynamic_pointer_cast<InteractingFacet>(b->interactingGeometry);
                if(!f) continue;
@@ -23,20 +33,89 @@
        size_t nVertices=vv.size(), j;
        for(size_t i=0; i<nVertices; i++){
                j=i;
-               while(++j<nVertices && (vv[j]->coord-vv[i]->coord)<tolerance){
+               while(++j<nVertices && (vv[j]->coord-vv[i]->coord)<=tolerance){
                        shared_ptr<VertexData> &vi=vv[i], &vj=vv[j];
-                       if(abs(vi->pos[0]-vj->pos[0])<tolerance &&
-                               abs(vi->pos[1]-vj->pos[1])<tolerance &&
-                               abs(vi->pos[2]-vj->pos[2])<tolerance &&
-                               
(vi->pos-vj->pos).SquaredLength()<tolerance*tolerance){
+                       if(abs(vi->pos[0]-vj->pos[0])<=tolerance &&
+                               abs(vi->pos[1]-vj->pos[1])<=tolerance &&
+                               abs(vi->pos[2]-vj->pos[2])<=tolerance &&
+                               
(vi->pos-vj->pos).SquaredLength()<=tolerance*tolerance){
                                // OK, these two vertices touch
-                               LOG_DEBUG("Vertices 
"<<vi->id<<"/"<<vi->vertexNo<<" and "<<vj->id<<"/"<<vj->vertexNo<<" close 
enough.");
+                               LOG_TRACE("Vertices 
"<<vi->id<<"/"<<vi->vertexNo<<" and "<<vj->id<<"/"<<vj->vertexNo<<" close 
enough.");
                                // add vertex to the nextIndetical of the one 
that has lower index; the one that is added will have isLowestIndex=false
                                if(vi->index<vj->index){ 
vi->nextIdentical.push_back(vj); vj->isLowestIndex=false; }
                                else{                    
vj->nextIdentical.push_back(vi); vi->isLowestIndex=false; }
                        }
                }
        }
-       /* TODO: create vertex clusters of touching vertices; find what facets 
they belong to; find whether facets that touch have common edges; compute 
mutual angle between faces on edge and save that to InteractingFacet */
+       // identity chains start always at lower indices, this way we get all 
of them
+       sort(vv.begin(),vv.end(),VertexIndexComparator());
+       int maxVertexId=0;
+       FOREACH(shared_ptr<VertexData>& v, vv){
+               if(v->vertexId<0){
+                       assert(v->isLowestIndex);
+                       v->vertexId=maxVertexId++;
+               }
+               FOREACH(shared_ptr<VertexData>& vNext, v->nextIdentical){
+                       vNext->vertexId=v->vertexId;
+               }
+               if(v->vertexId>=0) continue; // already assigned
+       }
+       LOG_DEBUG("Found "<<maxVertexId<<" unique vertices.");
+       // 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){
+               if(!topo[v->id]) topo[v->id]=shared_ptr<FacetTopology>(new 
FacetTopology(v->id));
+               topo[v->id]->vertices[v->vertexNo]=v->vertexId;
+       }
+       // make sure all facets have their vertex id's assigned
+       // add non-empty ones to topo1 that will be used for adjacency search 
afterwards
+       vector<shared_ptr<FacetTopology> > topo1;
+       for(size_t i=0; i<topo.size(); i++){
+               shared_ptr<FacetTopology> t=topo[i];
+               if(!t) continue;
+               if(t->vertices[0]<0 || t->vertices[1]<0 || t->vertices[2]<0){
+                       LOG_FATAL("Facet #"<<i<<": some vertex has no 
integrized vertexId assigned!!");
+                       LOG_FATAL("Vertices are: 
"<<t->vertices[0]<<","<<t->vertices[1]<<","<<t->vertices[2]);
+                       throw logic_error("Facet vertex has no integrized 
vertex assigned?!");
+               }
+               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
+                       assert(vvv.size()!=3); // same coords? nonsense
+                       assert(vvv.size()==2);
+                       vector<int> edge(2,0);
+                       // identify what edge are we at, for both facets
+                       for(int k=0; k<2; k++){
+                               for(edge[k]=0; edge[k]<3; edge[k]++){
+                                       const shared_ptr<FacetTopology>& tt( 
k==0 ? ti : tj);
+                                       size_t 
v1=tt->vertices[edge[k]],v2=tt->vertices[(edge[k]+1)%3];
+                                       if((vvv[0]==v1 && vvv[1]==v2) || 
(vvv[0]==v2 && vvv[1]==v1)) break;
+                                       if(edge[k]==2){
+                                               LOG_FATAL("No edge identified 
for 2 vertices "<<vvv[0]<<","<<vvv[1]<<" (facet #"<<tt->id<<" is: 
"<<tt->vertices[0]<<","<<tt->vertices[1]<<","<<tt->vertices[2]<<")");
+                                               throw logic_error("No edge 
identified for given vertices.");
+                                       }
+                               }
+                       }
+                       // 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;
+                       commonEdgesFound++;
+                       LOG_TRACE("Added adjacency information for 
#"<<ti->id<<"+#"<<tj->id<<" (common edges "<<edge[0]<<"+"<<edge[1]<<")");
+               }
+       }
 }
-
+#endif

Modified: trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp     
2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/pkg/dem/Engine/StandAloneEngine/FacetTopologyAnalyzer.hpp     
2009-04-08 20:27:04 UTC (rev 1750)
@@ -9,7 +9,7 @@
  */
 class FacetTopologyAnalyzer: public StandAloneEngine{
        struct VertexData{
-               VertexData(body_id_t _id, int _vertexNo, Vector3r _pos, Real 
_coord): id(_id), vertexNo(_vertexNo), coord(_coord), 
pos(_pos){index=3*id+vertexNo; isLowestIndex=true;}
+               VertexData(body_id_t _id, int _vertexNo, Vector3r _pos, Real 
_coord): id(_id), vertexNo(_vertexNo), coord(_coord), 
pos(_pos){index=3*id+vertexNo; isLowestIndex=true; vertexId=-1;}
                //! Facet (body id) that we represent
                body_id_t id;
                //! vertex number within this Facet
@@ -24,21 +24,37 @@
                bool isLowestIndex;
                //! vertices that are "identical" with this one, but have higer 
indices
                vector<shared_ptr<VertexData> > nextIdentical;
+               //! id of vertex, once they have id's assigned
+               long vertexId;
        };
        struct VertexComparator{
                bool operator()(const shared_ptr<VertexData>& v1, const 
shared_ptr<VertexData>& v2){return v1->coord<v2->coord;}
        };
-
+       struct VertexIndexComparator{
+               bool operator()(const shared_ptr<VertexData>& v1, const 
shared_ptr<VertexData>& v2){return v1->index<v2->index;}
+       };
+       struct FacetTopology{
+               FacetTopology(body_id_t _id): 
id(_id){vertices[0]=vertices[1]=vertices[2]=-1;}
+               //! integrized vertices
+               long vertices[3];
+               //! 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;
-               //! maximum distance of "identical" vertices, relative to 
maximum facet size
+               //! maximum distance of "identical" vertices, relative to 
minimum facet size
                Real relTolerance;
+               //! how many common edges were identified during last run
+               long commonEdgesFound;
        void action(MetaBody*); 
-       FacetTopologyAnalyzer(): projectionAxis(Vector3r::UNIT_X), 
relTolerance(1e-4) {}
+       FacetTopologyAnalyzer(): projectionAxis(Vector3r::UNIT_X), 
relTolerance(1e-4), commonEdgesFound(0) {}
        DECLARE_LOGGER;
        REGISTER_CLASS_AND_BASE(FacetTopologyAnalyzer,StandAloneEngine);
-       REGISTER_ATTRIBUTES(StandAloneEngine, /* none */);
+       REGISTER_ATTRIBUTES(StandAloneEngine, (projectionAxis)(relTolerance));
 };
 REGISTER_SERIALIZABLE(FacetTopologyAnalyzer);
 

Added: trunk/scripts/test/facet-topo.py
===================================================================
--- trunk/scripts/test/facet-topo.py    2009-04-06 21:40:30 UTC (rev 1749)
+++ trunk/scripts/test/facet-topo.py    2009-04-08 20:27:04 UTC (rev 1750)
@@ -0,0 +1,24 @@
+import yade.log
+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'),]
+
+# most simple case: no touch at all
+if 1:
+       O.bodies.append([
+               utils.facet([(0,0,0),(1,0,0),(0,1,0)]),
+               utils.facet([(0,0,1),(1,0,1),(0,1,1)]),
+       ])
+       O.step();
+       assert(topo['commonEdgesFound']==0)
+if 1:
+       O.bodies.clear()
+       O.bodies.append([
+               utils.facet([(0,0,0),(1,0,0),(0,1,0)]),
+               utils.facet([(1,1,0),(1,0,0),(0,1,0)]),
+       ])
+       O.step()
+       #assert(O.bodies[0].phys['edgeAdjIds'][1]==1 and 
O.bodies[1].phys['edgeAdjIds'][0]==1)
+       assert(topo['commonEdgesFound']==1)


_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to     : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help   : https://help.launchpad.net/ListHelp
_______________________________________________
yade-dev mailing list
[email protected]
https://lists.berlios.de/mailman/listinfo/yade-dev

Reply via email to