yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #09939
[Branch ~yade-pkg/yade/git-trunk] Rev 3686: Tetra modification (Ig2_Tetra_Tetra_TTtraSimpleGeom, Law2_TTetraSimpleGeom_NormPhys_Simple, modif...
------------------------------------------------------------
revno: 3686
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Sun 2013-09-08 13:12:42 +0200
message:
Tetra modification (Ig2_Tetra_Tetra_TTtraSimpleGeom, Law2_TTetraSimpleGeom_NormPhys_Simple, modified GL functor, examples, utils.tetra)
added Ip2_ElastMat_... functors
fixed bug in utils.kineticEnergy for aspherical bodies
added:
examples/polyhedra/
examples/polyhedra/oneTetra.py
examples/polyhedra/twoTetras.py
pkg/dem/Ip2_ElastMat.cpp
pkg/dem/Ip2_ElastMat.hpp
modified:
lib/opengl/OpenGLWrapper.hpp
pkg/common/Facet.cpp
pkg/common/Facet.hpp
pkg/dem/Shop.cpp
pkg/dem/Shop.hpp
pkg/dem/Tetra.cpp
pkg/dem/Tetra.hpp
py/_utils.cpp
py/utils.py
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== added directory 'examples/polyhedra'
=== added file 'examples/polyhedra/oneTetra.py'
--- examples/polyhedra/oneTetra.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/oneTetra.py 2013-09-08 11:12:42 +0000
@@ -0,0 +1,42 @@
+################################################################################
+#
+# Script to test tetra gl functions and prescribe dmotion
+#
+################################################################################
+v1 = (0,0,0)
+v2 = (1,0,0)
+v3 = (0,1,0)
+v4 = (0,0,1)
+
+t1 = tetra((v1,v2,v3,v4),color=(0,1,0))
+O.bodies.append((t1))
+
+def changeVelocity():
+ if O.iter == 50000:
+ t1.state.vel = (-1,0,0)
+ if O.iter == 100000:
+ t1.state.vel = (0,1,-1)
+ if O.iter == 150000:
+ t1.state.vel = (0,0,0)
+ t1.state.angMom = (0,0,10)
+ if O.iter == 200000:
+ O.pause()
+
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Tetra_Aabb()]),
+ InteractionLoop([],[],[]),
+ NewtonIntegrator(),
+ PyRunner(iterPeriod=1,command="changeVelocity()"),
+]
+
+
+try:
+ from yade import qt
+ qt.View()
+except:
+ pass
+
+O.dt = 1e-5
+t1.state.vel = (1,0,0)
+O.run()
=== added file 'examples/polyhedra/twoTetras.py'
--- examples/polyhedra/twoTetras.py 1970-01-01 00:00:00 +0000
+++ examples/polyhedra/twoTetras.py 2013-09-08 11:12:42 +0000
@@ -0,0 +1,111 @@
+################################################################################
+#
+# Python script to test tetra-tetra contact detection for different possible
+# contact modes. Some tests are run twice to test the symmetry of the law.
+#
+# It runs several tests, making pause before each one. If run with GUI, you can
+# adjust the viewer
+#
+# During the test, momentum, angular momentum and kinetic energy is tracked in
+# plot.data. You can use plot.plot() to see the results (in sime modes the
+# energy is not conserved for instace).
+#
+################################################################################
+from yade import qt,plot
+
+O.materials.append(ElastMat(young=1e10))
+O.dt = 4e-5
+O.engines = [
+ ForceResetter(),
+ InsertionSortCollider([Bo1_Tetra_Aabb()]),
+ InteractionLoop(
+ [Ig2_Tetra_Tetra_TTetraSimpleGeom()],
+ [Ip2_ElastMat_ElastMat_NormPhys()],
+ [Law2_TTetraSimpleGeom_NormPhys_Simple()]
+ ),
+ NewtonIntegrator(damping=0),
+ PyRunner(iterPeriod=500,command="addPlotData()"),
+ PyRunner(iterPeriod=1000,command="printt()"),
+ PyRunner(iterPeriod=1,command="runExamples()"),
+]
+
+def addPlotData():
+ p = utils.momentum()
+ l = utils.angularMomentum()
+ plot.addData(
+ i = O.iter,
+ e = utils.kineticEnergy(),
+ px = p[0],
+ py = p[1],
+ pz = p[2],
+ lx = l[0],
+ ly = l[1],
+ lz = l[2],
+ )
+
+def prepareExample(vertices1,vertices2,vel1=(0,0,0),vel2=(0,0,0),amom1=(0,0,0),amom2=(0,0,0),label=''):
+ O.interactions.clear()
+ O.bodies.clear()
+ t1 = tetra(vertices1,color=(0,1,0),wire=False)
+ t2 = tetra(vertices2,color=(0,1,1),wire=False)
+ O.bodies.append((t1,t2))
+ t1.state.vel = vel1
+ t2.state.vel = vel2
+ t1.state.angMom = amom1
+ t2.state.angMom = amom2
+ if label: print label
+ O.pause()
+ O.step()
+
+def runExamples():
+ dt = 20000
+ if O.iter == 1:
+ vertices1 = ((0,0,0),(2,0,0),(0,2,0),(0,0,2))
+ vertices2 = ((1,1,1),(3,1,1),(1,3,1),(1,1,3))
+ prepareExample(vertices1,vertices2,vel2=(-1,-1,-1),label='\ntesting vertex-triangle contact...\n')
+ if O.iter == 1*dt:
+ plot.data = {}
+ vertices1 = ((1,1,1),(3,1,1),(1,3,1),(1,1,3))
+ vertices2 = ((0,0,0),(2,0,0),(0,2,0),(0,0,2))
+ prepareExample(vertices1,vertices2,vel1=(-1,-1,-1),label='\ntesting vertex-triangle contact 2...\n')
+ elif O.iter == 2*dt:
+ vertices1 = ((0,0,0),(0,0,1.1),(1,0,1),(0,1,.9))
+ vertices2 = ((0,.5,1.4),(-.5,.5,.6),(-1.6,0,1),(-1.6,1.1,1))
+ prepareExample(vertices1,vertices2,vel2=(1,0,0),amom2=(0,10,0),label='\ntesting edge-edge contact\n')
+ elif O.iter == 3*dt:
+ vertices1 = ((0,0,0),(0,0,1.1),(1,0,1),(0,1,.9))
+ vertices2 = ((-.5,.5,.6),(0,.5,1.4),(-1.6,1.1,1),(-1.6,0,1))
+ prepareExample(vertices1,vertices2,vel2=(1,0,0),amom2=(0,10,0),label='\ntesting edge-edge contact\n')
+ elif O.iter == 4*dt:
+ vertices1 = ((.1,-.4,-.3),(-.3,-.4,2),(3,-.2,2),(-.1,3,2))
+ vertices2 = ((.5,1.5,2.3),(1.5,.5,2.3),(2,2,3),(0,0,3))
+ prepareExample(vertices1,vertices2,vel2=(0,0,-1),amom2=(0,0,0),label='\ntesting edge-triangle contact\n')
+ elif O.iter == 5*dt:
+ vertices1 = ((.1,-.4,-.3),(-.3,-.4,2),(3,-.2,2),(-.1,3,2))
+ vertices2 = ((-.5,2.5,2.3),(2.5,-.5,2.3),(2,2,3),(0,0,3))
+ prepareExample(vertices1,vertices2,vel2=(0,0,-1),amom2=(0,0,0),label='\ntesting edge-triangle contact 2\n')
+ elif O.iter == 6*dt:
+ vertices1 = ((1,0,0),(0,1,0),(0,0,1),(1,1,1))
+ vertices2 = ((.5,.5,1.2),(0,.1,2),(2,0,2),(0,1,2))
+ prepareExample(vertices1,vertices2,vel2=(0,0,-1),label='\ntesting vertex-edge contact\n')
+ elif O.iter == 7*dt:
+ vertices1 = ((.5,.5,1.2),(0,.1,2),(2,0,2),(0,1,2))
+ vertices2 = ((1,0,0),(0,1,0),(0,0,1),(1,1,1))
+ prepareExample(vertices1,vertices2,vel1=(0,0,-1),label='\ntesting vertex-edge contact 2\n')
+ elif O.iter == 8*dt:
+ vertices1 = ((0,0,0),(1,0,0),(0,1,0),(0,0,1))
+ vertices2 = ((0,0,1.2),(.9,.8,2),(-.7,.61,2.3),(0,-.7,2.1))
+ prepareExample(vertices1,vertices2,vel2=(0,0,-1),label='\ntesting vertex-edge contact\n')
+ elif O.iter == 9*dt:
+ vertices1 = ((0,0,1.2),(.9,.8,2),(-.7,.61,2.3),(0,-.7,2.1))
+ vertices2 = ((0,0,0),(1,0,0),(0,1,0),(0,0,1))
+ prepareExample(vertices1,vertices2,vel1=(0,0,-1),label='\ntesting vertex-edge contact 2\n')
+ elif O.iter == 10*dt:
+ O.pause()
+
+viewer = qt.View()
+plot.plots = {'i':'e','i ':('px','py','pz'),'i ':('lx','ly','lz')}
+
+O.step()
+O.step()
+O.step()
=== modified file 'lib/opengl/OpenGLWrapper.hpp'
--- lib/opengl/OpenGLWrapper.hpp 2012-08-02 18:13:22 +0000
+++ lib/opengl/OpenGLWrapper.hpp 2013-09-08 11:12:42 +0000
@@ -26,71 +26,76 @@
/// Primary Templates
-template< typename Type > inline void glRotate ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glScale ( Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glRotate ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glScale ( Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glScalev ( const Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glTranslate ( Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glTranslate ( Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTranslatev ( const Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glVertex2 ( Type ,Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glVertex3 ( Type ,Type , Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glVertex4 ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glVertex2 ( Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glVertex3 ( Type, Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glVertex4 ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glVertex2v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glVertex3v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glVertex4v ( const Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glNormal3 ( Type ,Type ,Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glNormal3 ( Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glNormal3v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glIndex ( Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glIndexv ( Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glColor3 ( Type ,Type ,Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glColor4 ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glColor3 ( Type, Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glColor4 ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glColor3v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glColor4v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTexCoord1 ( Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glTexCoord2 ( Type ,Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glTexCoord3 ( Type ,Type , Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glTexCoord4 ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glTexCoord2 ( Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glTexCoord3 ( Type, Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glTexCoord4 ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTexCoord1v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTexCoord2v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTexCoord3v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glTexCoord4v ( const Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glRasterPos2 ( Type ,Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glRasterPos3 ( Type ,Type , Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glRasterPos4 ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glRasterPos2 ( Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glRasterPos3 ( Type, Type, Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glRasterPos4 ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glRasterPos2v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glRasterPos3v ( const Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glRasterPos4v ( const Type ) { STATIC_ASSERT(false); };
-template< typename Type > inline void glRect ( Type ,Type ,Type , Type ) { STATIC_ASSERT(false); };
+template< typename Type > inline void glRect ( Type, Type, Type, Type ) { STATIC_ASSERT(false); };
template< typename Type > inline void glMaterial ( GLenum face, GLenum pname, Type param ){ STATIC_ASSERT(false); };
template< typename Type > inline void glMaterialv ( GLenum face, GLenum pname, Type param ){ STATIC_ASSERT(false); };
template< typename Type > inline void glMultMatrix (const Type*){ STATIC_ASSERT(false); };
+#define LDOUBL long double
/// Template Specializations
template< > inline void glMultMatrix<double>(const double* m){glMultMatrixd(m); };
-template< > inline void glMultMatrix<long double>(const long double* m){double mm[16]; for(int i=0;i<16;i++)mm[i]=(double)m[i]; glMultMatrixd(mm);};
+template< > inline void glMultMatrix<LDOUBL>(const LDOUBL* m){double mm[16]; for(int i=0;i<16;i++)mm[i]=(double)m[i]; glMultMatrixd(mm);};
template< > inline void glRotate< double > (double angle,double x,double y, double z ) { glRotated(angle,x,y,z); };
+template< > inline void glRotate< LDOUBL > (LDOUBL angle, LDOUBL x, LDOUBL y, LDOUBL z ) { glRotated((double)angle,(double)x,(double)y,(double)z); };
template< > inline void glRotate< float > (float angle,float x,float y, float z ) { glRotatef(angle,x,y,z); };
template< > inline void glScale< double > ( double x,double y, double z ) { glScaled(x,y,z); };
-template< > inline void glScale< long double > ( long double x,long double y, long double z ) { glScaled(x,y,z); };
+template< > inline void glScale< LDOUBL > ( LDOUBL x,LDOUBL y, LDOUBL z ) { glScaled(x,y,z); };
template< > inline void glScale< float > ( float x,float y,float z ) { glScalef(x,y,z); };
template< > inline void glScalev< Vector3r > ( const Vector3r v ) { glScaled(v[0],v[1],v[2]);};
template< > inline void glTranslate< double > ( double x,double y, double z ) { glTranslated(x,y,z); };
-template< > inline void glTranslate< long double > ( long double x, long double y, long double z ) { glTranslated(x,y,z); };
+template< > inline void glTranslate< LDOUBL > ( LDOUBL x, LDOUBL y, LDOUBL z ) { glTranslated(x,y,z); };
template< > inline void glTranslate< float > ( float x,float y,float z ) { glTranslatef(x,y,z); };
template< > inline void glTranslatev< Vector3r > ( const Vector3r v ) { glTranslated(v[0],v[1],v[2]);};
template< > inline void glVertex2< double > ( double x,double y ) { glVertex2d(x,y); };
+template< > inline void glVertex2< LDOUBL > ( LDOUBL x, LDOUBL y ) { glVertex2d((double)x,(double)y); };
template< > inline void glVertex2< float > ( float x,float y ) { glVertex2f(x,y); };
template< > inline void glVertex2< int > ( int x,int y ) { glVertex2i(x,y); };
template< > inline void glVertex3< double > ( double x,double y, double z ) { glVertex3d(x,y,z); };
+template< > inline void glVertex3< LDOUBL > ( LDOUBL x, LDOUBL y, LDOUBL z ) { glVertex3d((double)x,(double)y,(double)z); };
template< > inline void glVertex3< float > ( float x,float y,float z ) { glVertex3f(x,y,z); };
template< > inline void glVertex3< int > ( int x, int y, int z ) { glVertex3i(x,y,z); };
template< > inline void glVertex4< double > ( double x,double y,double z, double w ){ glVertex4d(x,y,z,w); };
+template< > inline void glVertex4< LDOUBL > ( LDOUBL x, LDOUBL y, LDOUBL z, LDOUBL w ){ glVertex4d((double)x,(double)y,(double)z,(double)w); };
template< > inline void glVertex4< float > ( float x,float y,float z, float w ) { glVertex4f(x,y,z,w); };
template< > inline void glVertex4< int > ( int x,int y,int z, int w ) { glVertex4i(x,y,z,w); };
@@ -105,12 +110,14 @@
template< > inline void glVertex4v< Vector3i > ( const Vector3i v ) { glVertex4iv((int*)&v); };
template< > inline void glNormal3< double > (double nx,double ny,double nz ) { glNormal3d(nx,ny,nz); };
+template< > inline void glNormal3< LDOUBL > ( LDOUBL nx, LDOUBL ny, LDOUBL nz ) { glNormal3d((double)nx,(double)ny,(double)nz); };
template< > inline void glNormal3< int > (int nx,int ny,int nz ) { glNormal3i(nx,ny,nz); };
template< > inline void glNormal3v< Vector3r > ( const Vector3r v ) { glNormal3dv((double*)&v); };
template< > inline void glNormal3v< Vector3i > ( const Vector3i v ) { glNormal3iv((int*)&v); };
template< > inline void glIndex< double > ( double c ) { glIndexd(c); };
+template< > inline void glIndex< LDOUBL > ( LDOUBL c ) { glIndexd((double)c); };
template< > inline void glIndex< float > ( float c ) { glIndexf(c); };
template< > inline void glIndex< int > ( int c ) { glIndexi(c); };
template< > inline void glIndex< unsigned char > ( unsigned char c ) { glIndexub(c); };
@@ -119,10 +126,12 @@
template< > inline void glIndexv<const Vector3i > ( const Vector3i c) { glIndexiv((int*)&c); }
template< > inline void glColor3< double > (double red,double green,double blue ) { glColor3d(red,green,blue); };
+template< > inline void glColor3< LDOUBL > ( LDOUBL red, LDOUBL green, LDOUBL blue ) { glColor3d((double)red,(double)green,(double)blue); };
template< > inline void glColor3< float > (float red,float green,float blue ) { glColor3f(red,green,blue); };
template< > inline void glColor3< int > (int red,int green,int blue ) { glColor3i(red,green,blue); };
template< > inline void glColor4< double > (double red,double green,double blue, double alpha ) { glColor4d(red,green,blue,alpha); };
+template< > inline void glColor4< LDOUBL > (LDOUBL red,LDOUBL green,LDOUBL blue, LDOUBL alpha ) { glColor4d((double)red,(double)green,(double)blue,(double)alpha); };
template< > inline void glColor4< float > (float red,float green,float blue, float alpha ) { glColor4f(red,green,blue,alpha); };
template< > inline void glColor4< int > (int red,int green,int blue, int alpha ) { glColor4i(red,green,blue,alpha); };
@@ -135,18 +144,22 @@
template< > inline void glTexCoord1< double > ( double s ) { glTexCoord1d(s); };
+template< > inline void glTexCoord1< LDOUBL > ( LDOUBL s ) { glTexCoord1d((double)s); };
template< > inline void glTexCoord1< float > ( float s ) { glTexCoord1f(s); };
template< > inline void glTexCoord1< int > ( int s ) { glTexCoord1i(s); };
template< > inline void glTexCoord2< double > ( double s,double t ) { glTexCoord2d(s,t); };
+template< > inline void glTexCoord2< LDOUBL > ( LDOUBL s,LDOUBL t ) { glTexCoord2d((double)s,(double)t); };
template< > inline void glTexCoord2< float > ( float s,float t ) { glTexCoord2f(s,t); };
template< > inline void glTexCoord2< int > ( int s,int t ) { glTexCoord2i(s,t); };
template< > inline void glTexCoord3< double > ( double s,double t, double r ) { glTexCoord3d(s,t,r); };
+template< > inline void glTexCoord3< LDOUBL > ( LDOUBL s,LDOUBL t, LDOUBL r ) { glTexCoord3d((double)s,(double)t,(double)r); };
template< > inline void glTexCoord3< float > ( float s,float t,float r ) { glTexCoord3f(s,t,r); };
template< > inline void glTexCoord3< int > ( int s, int t, int r ) { glTexCoord3i(s,t,r); };
template< > inline void glTexCoord4< double > (double s,double t,double r, double q ) { glTexCoord4d(s,t,r,q); };
+template< > inline void glTexCoord4< LDOUBL > (LDOUBL s,LDOUBL t,LDOUBL r, LDOUBL q ) { glTexCoord4d((double)s,(double)t,(double)r,(double)q); };
template< > inline void glTexCoord4< float > (float s,float t,float r, float q ) { glTexCoord4f(s,t,r,q); };
template< > inline void glTexCoord4< int > (int s,int t,int r, int q ) { glTexCoord4i(s,t,r,q); };
@@ -164,14 +177,17 @@
template< > inline void glRasterPos2< double > ( double x,double y ) { glRasterPos2d(x,y); };
+template< > inline void glRasterPos2< LDOUBL > ( LDOUBL x,LDOUBL y ) { glRasterPos2d((double)x,(double)y); };
template< > inline void glRasterPos2< float > ( float x,float y ) { glRasterPos2f(x,y); };
template< > inline void glRasterPos2< int > ( int x,int y ) { glRasterPos2i(x,y); };
template< > inline void glRasterPos3< double > ( double x,double y, double z ) { glRasterPos3d(x,y,z); };
+template< > inline void glRasterPos3< LDOUBL > ( LDOUBL x,LDOUBL y, LDOUBL z ) { glRasterPos3d((double)x,(double)y,(double)z); };
template< > inline void glRasterPos3< float > ( float x,float y,float z ) { glRasterPos3f(x,y,z); };
template< > inline void glRasterPos3< int > ( int x, int y, int z ) { glRasterPos3i(x,y,z); };
template< > inline void glRasterPos4< double > (double x,double y,double z, double w ) { glRasterPos4d(x,y,z,w); };
+template< > inline void glRasterPos4< LDOUBL > (LDOUBL x,LDOUBL y,LDOUBL z, LDOUBL w ) { glRasterPos4d((double)x,(double)y,(double)z,(double)w); };
template< > inline void glRasterPos4< float > (float x,float y,float z, float w ) { glRasterPos4f(x,y,z,w); };
template< > inline void glRasterPos4< int > (int x,int y,int z, int w ) { glRasterPos4i(x,y,z,w); };
@@ -188,6 +204,7 @@
template< > inline void glRect< double > (double x1,double y1,double x2, double y2 ) { glRectd(x1,y1,x2,y2); };
+template< > inline void glRect< LDOUBL > (LDOUBL x1,LDOUBL y1,LDOUBL x2, LDOUBL y2 ) { glRectd((double)x1,(double)y1,(double)x2,(double)y2); };
template< > inline void glRect< float > (float x1,float y1,float x2,float y2 ) { glRectf(x1,y1,x2,y2); };
template< > inline void glRect< int > (int x1,int y1,int x2, int y2 ) { glRecti(x1,y1,x2,y2); };
@@ -197,4 +214,4 @@
template< > inline void glMaterialv< Vector3r > ( GLenum face, GLenum pname, const Vector3r params ) { const GLfloat _p[3]={(float) params[0], (float) params[1], (float) params[2]}; glMaterialfv(face,pname,_p); };
template< > inline void glMaterialv< Vector3i > ( GLenum face, GLenum pname, const Vector3i params ) { glMaterialiv(face,pname,(int*)¶ms); };
-
+#undef LDOUBL
=== modified file 'pkg/common/Facet.cpp'
--- pkg/common/Facet.cpp 2013-05-24 16:48:19 +0000
+++ pkg/common/Facet.cpp 2013-09-08 11:12:42 +0000
@@ -28,7 +28,6 @@
normal = e[0].cross(e[1]);
area = .5*normal.norm();
normal /= 2*area;
- //normal.normalize();
for(int i=0; i<3; ++i){
ne[i]=e[i].cross(normal); ne[i].normalize();
vl[i]=vertices[i].norm();
=== modified file 'pkg/common/Facet.hpp'
--- pkg/common/Facet.hpp 2013-05-24 16:48:19 +0000
+++ pkg/common/Facet.hpp 2013-09-08 11:12:42 +0000
@@ -37,7 +37,7 @@
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Facet,Shape,"Facet (triangular particle) geometry.",
((vector<Vector3r>,vertices,vector<Vector3r>(3,Vector3r(NaN,NaN,NaN)),(Attr::triggerPostLoad | Attr::noResize),"Vertex positions in local coordinates."))
- ((Vector3r,normal,Vector3r(NaN,NaN,NaN),(Attr::readonly | Attr::noSave),"Facet's normal"))
+ ((Vector3r,normal,Vector3r(NaN,NaN,NaN),(Attr::readonly | Attr::noSave),"Facet's normal (in local coordinate system)"))
((Real,area,NaN,(Attr::readonly | Attr::noSave),"Facet's area"))
#ifdef FACET_TOPO
((vector<Body::id_t>,edgeAdjIds,vector<Body::id_t>(3,Body::ID_NONE),,"Facet id's that are adjacent to respective edges [experimental]"))
=== added file 'pkg/dem/Ip2_ElastMat.cpp'
--- pkg/dem/Ip2_ElastMat.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Ip2_ElastMat.cpp 2013-09-08 11:12:42 +0000
@@ -0,0 +1,64 @@
+#include "Ip2_ElastMat.hpp"
+
+#include <yade/pkg/common/NormShearPhys.hpp>
+#include <yade/pkg/common/NormShearPhys.hpp>
+#include <yade/pkg/dem/DemXDofGeom.hpp>
+
+YADE_PLUGIN((Ip2_ElastMat_ElastMat_NormPhys)(Ip2_ElastMat_ElastMat_NormShearPhys));
+
+
+void Ip2_ElastMat_ElastMat_NormPhys::go( const shared_ptr<Material>& b1
+ , const shared_ptr<Material>& b2
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return;
+ const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
+ const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ interaction->phys = shared_ptr<NormPhys>(new NormPhys());
+ const shared_ptr<NormPhys>& phys = YADE_PTR_CAST<NormPhys>(interaction->phys);
+ Real Kn;
+ const GenericSpheresContact* geom=dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
+ if (geom) {
+ Real Ra,Rb;//Vector3r normal;
+ Ra=geom->refR1>0?geom->refR1:geom->refR2;
+ Rb=geom->refR2>0?geom->refR2:geom->refR1;
+ //harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
+ Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
+ } else {
+ Kn = 2*Ea*Eb/(Ea+Eb);
+ }
+ phys->kn = Kn;
+};
+
+
+void Ip2_ElastMat_ElastMat_NormShearPhys::go( const shared_ptr<Material>& b1
+ , const shared_ptr<Material>& b2
+ , const shared_ptr<Interaction>& interaction)
+{
+ if(interaction->phys) return;
+ const shared_ptr<ElastMat>& mat1 = YADE_PTR_CAST<ElastMat>(b1);
+ const shared_ptr<ElastMat>& mat2 = YADE_PTR_CAST<ElastMat>(b2);
+ Real Ea = mat1->young;
+ Real Eb = mat2->young;
+ Real Va = mat1->poisson;
+ Real Vb = mat2->poisson;
+ interaction->phys = shared_ptr<NormShearPhys>(new NormShearPhys());
+ const shared_ptr<NormShearPhys>& phys = YADE_PTR_CAST<NormShearPhys>(interaction->phys);
+ Real Kn, Ks;
+ GenericSpheresContact* geom=dynamic_cast<GenericSpheresContact*>(interaction->geom.get());
+ if (geom) {
+ Real Ra,Rb;//Vector3r normal;
+ Ra=geom->refR1>0?geom->refR1:geom->refR2;
+ Rb=geom->refR2>0?geom->refR2:geom->refR1;
+ //harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
+ Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
+ Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);
+ } else {
+ Kn = 2*Ea*Eb/(Ea+Eb);
+ Kn = 2*Ea*Va*Eb*Vb/(Ea*Va+Eb*Vb);
+ }
+ phys->kn = Kn;
+ phys->ks = Ks;
+};
=== added file 'pkg/dem/Ip2_ElastMat.hpp'
--- pkg/dem/Ip2_ElastMat.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Ip2_ElastMat.hpp 2013-09-08 11:12:42 +0000
@@ -0,0 +1,24 @@
+
+#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/common/MatchMaker.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+
+class Ip2_ElastMat_ElastMat_NormPhys: public IPhysFunctor{
+ public:
+ virtual void go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(ElastMat,ElastMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_ElastMat_ElastMat_NormPhys,IPhysFunctor,"Create a :yref:`NormPhys` from two :yref:`ElastMats<ElastMat>`. TODO. EXPERIMENTAL",
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_ElastMat_ElastMat_NormPhys);
+
+
+class Ip2_ElastMat_ElastMat_NormShearPhys: public IPhysFunctor{
+ public:
+ virtual void go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction);
+ FUNCTOR2D(ElastMat,ElastMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_ElastMat_ElastMat_NormShearPhys,IPhysFunctor,"Create a :yref:`NormShearPhys` from two :yref:`ElastMats<ElastMat>`. TODO. EXPERIMENTAL",
+ );
+};
+REGISTER_SERIALIZABLE(Ip2_ElastMat_ElastMat_NormShearPhys);
+
=== modified file 'pkg/dem/Shop.cpp'
--- pkg/dem/Shop.cpp 2013-06-20 14:56:35 +0000
+++ pkg/dem/Shop.cpp 2013-09-08 11:12:42 +0000
@@ -206,7 +206,8 @@
// the tensorial expression http://en.wikipedia.org/wiki/Moment_of_inertia#Moment_of_inertia_tensor
// inertia tensor rotation from http://www.kwon3d.com/theory/moi/triten.html
Matrix3r mI; mI<<state->inertia[0],0,0, 0,state->inertia[1],0, 0,0,state->inertia[2];
- E+=.5*state->angVel.transpose().dot((T.transpose()*mI*T)*state->angVel);
+ //E+=.5*state->angVel.transpose().dot((T.transpose()*mI*T)*state->angVel);
+ E+=.5*state->angVel.transpose().dot((T*mI*T.transpose())*state->angVel);
}
else { E+=0.5*state->angVel.dot(state->inertia.cwiseProduct(state->angVel));}
if(maxId && E>maxE) { *maxId=b->getId(); maxE=E; }
@@ -215,6 +216,27 @@
return ret;
}
+Vector3r Shop::momentum() {
+ Vector3r ret = Vector3r::Zero();
+ Scene* scene=Omega::instance().getScene().get();
+ FOREACH(const shared_ptr<Body> b, *scene->bodies){
+ ret += b->state->mass * b->state->vel;
+ }
+ return ret;
+}
+
+Vector3r Shop::angularMomentum(Vector3r origin) {
+ Vector3r ret = Vector3r::Zero();
+ Scene* scene=Omega::instance().getScene().get();
+ Matrix3r T,Iloc;
+ FOREACH(const shared_ptr<Body> b, *scene->bodies){
+ ret += (b->state->pos-origin).cross(b->state->mass * b->state->vel);
+ ret += b->state->angMom;
+ }
+ return ret;
+}
+
+
shared_ptr<FrictMat> Shop::defaultGranularMat(){
shared_ptr<FrictMat> mat(new FrictMat);
mat->density=2e3;
=== modified file 'pkg/dem/Shop.hpp'
--- pkg/dem/Shop.hpp 2013-05-30 20:19:43 +0000
+++ pkg/dem/Shop.hpp 2013-09-08 11:12:42 +0000
@@ -84,6 +84,10 @@
//! Get unbalanced force of the whole simulation
static Real unbalancedForce(bool useMaxForce=false, Scene* _rb=NULL);
static Real kineticEnergy(Scene* _rb=NULL, Body::id_t* maxId=NULL);
+ //! get total momentum of current simulation
+ static Vector3r momentum();
+ //! get total angular momentum of current simulation
+ static Vector3r angularMomentum(Vector3r origin = Vector3r::Zero());
static Vector3r totalForceInVolume(Real& avgIsoStiffness, Scene *_rb=NULL);
=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp 2012-02-16 16:05:15 +0000
+++ pkg/dem/Tetra.cpp 2013-09-08 11:12:42 +0000
@@ -1,16 +1,23 @@
// © 2007 Václav Šmilauer <eudoxos@xxxxxxxx>
+// © 2013 Jan Stránský <jan.stransky@xxxxxxxxxxx>
#include"Tetra.hpp"
-YADE_PLUGIN(/* self-contained in hpp: */ (Tetra) (TTetraGeom) (Bo1_Tetra_Aabb)
- /* some code in cpp (this file): */ (TetraVolumetricLaw) (Ig2_Tetra_Tetra_TTetraGeom)
+YADE_PLUGIN(/* self-contained in hpp: */ (Tetra) (TTetraGeom) (TTetraSimpleGeom) (Bo1_Tetra_Aabb)
+ /* some code in cpp (this file): */ (TetraVolumetricLaw)
+ (Ig2_Tetra_Tetra_TTetraGeom)
+ #ifdef YADE_CGAL
+ (Ig2_Tetra_Tetra_TTetraSimpleGeom)
+ #endif
#ifdef YADE_OPENGL
(Gl1_Tetra)
#endif
+ (Law2_TTetraSimpleGeom_NormPhys_Simple)
);
Tetra::~Tetra(){}
TTetraGeom::~TTetraGeom(){}
+TTetraSimpleGeom::~TTetraSimpleGeom(){}
#include<yade/core/Interaction.hpp>
#include<yade/core/Omega.hpp>
@@ -20,6 +27,11 @@
#include<yade/pkg/common/Aabb.hpp>
+#ifdef YADE_CGAL
+ #include <CGAL/intersections.h>
+#endif
+
+
void Bo1_Tetra_Aabb::go(const shared_ptr<Shape>& ig, shared_ptr<Bound>& bv, const Se3r& se3, const Body*){
Tetra* t=static_cast<Tetra*>(ig.get());
if(!bv){ bv=shared_ptr<Bound>(new Aabb); }
@@ -32,6 +44,569 @@
#undef __VOP
}
+
+
+
+#ifdef YADE_CGAL
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::psMap[4][3] = { // segments of point
+ {0,2,3},
+ {0,1,4},
+ {1,2,5},
+ {3,4,5}
+};
+
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::ptMap[4][3] = { // triangles of point
+ {0,1,3},
+ {0,1,2},
+ {1,2,3},
+ {0,2,3}
+};
+
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::stMap[6][2] = { // triangles of segments
+ {0,1},
+ {1,2},
+ {1,3},
+ {0,3},
+ {0,2},
+ {2,3}
+};
+
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::ppsMap[4][4] = { // point-point pair to segment
+ {-1, 0, 2, 3},
+ { 0,-1, 1, 4},
+ { 2, 1,-1, 5},
+ { 3, 4, 5,-1}
+};
+
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::tsMap[4][3] = { // segmnts of triangle
+ {0,3,4},
+ {0,1,2},
+ {1,4,5},
+ {2,3,5}
+};
+
+const int Ig2_Tetra_Tetra_TTetraSimpleGeom::sstMap[6][6] = { // segment-segment pair to triangle
+ {-1, 1, 1, 0, 0,-1},
+ { 1,-1, 1,-1, 2, 2},
+ { 1, 1,-1, 3,-1, 3},
+ { 0,-1, 3,-1, 0, 3},
+ { 0, 2,-1, 0,-1, 2},
+ {-1, 2, 3, 3, 2,-1}
+};
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkVertexToTriangleCase(
+ const Triangle tA[4],
+ const Point pB[4],
+ const Segment sB[6],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+ for (int i=0; i<4; i++) { // loop over triangles 1
+ const Triangle& t = tA[i]; // choose triangle 1
+ for (int j=0; j<4; j++) { // loop over vertices 2
+ const Point& p = pB[j]; // choose vertex 2
+ // choose edges posessing p
+ const Segment& sa = sB[psMap[j][0]];
+ const Segment& sb = sB[psMap[j][1]];
+ const Segment& sc = sB[psMap[j][2]];
+ if ( !(do_intersect(t,sa) && do_intersect(t,sb) && do_intersect(t,sc)) ) { continue; }// if all edges intersect with t
+ // evaluate the points
+ CGAL::Object oa = intersection(t,sa);
+ const Point* pa = CGAL::object_cast<Point>(&oa);
+ CGAL::Object ob = intersection(t,sb);
+ const Point* pb = CGAL::object_cast<Point>(&ob);
+ CGAL::Object oc = intersection(t,sc);
+ const Point* pc = CGAL::object_cast<Point>(&oc);
+ if ( !(pa && pb && pc) ) { continue; } // check that the intrsection really exists
+ Vector_3 n = CGAL::normal(t[0],t[1],t[2]); // normal of triangle
+ for (int k=0; k<3; k++) {
+ normal[k] = n[k]; // sets normal of contact = nornal of triangle
+ // contact point is center of mass of overlaping tetrahedron
+ contactPoint[k] = .25*(p[k]+pa->operator[](k)+pb->operator[](k)+pc->operator[](k));
+ }
+ normal.normalize();
+ const Point* v[4] = {&p,pa,pb,pc};
+ penetrationVolume = TetrahedronVolume(v);
+ Real vol = TetrahedronVolume(pB);
+ if (penetrationVolume > .5*vol) { penetrationVolume = vol-penetrationVolume; }
+ return true;
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkEdgeToEdgeCase(
+ const Segment sA[6],
+ const Segment sB[6],
+ const Triangle tA[4],
+ const Triangle tB[4],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+ for (int i=0; i<6; i++) {
+ const Segment& sa = sA[i];
+ bool b[4];
+ const Triangle& ta0 = tA[stMap[i][0]];
+ const Triangle& ta1 = tA[stMap[i][1]];
+ for (int j=0; j<6; j++) {
+ const Segment sb = sB[j];
+ if ( !(do_intersect(sb,ta0) && do_intersect(sb,ta1) ) ) { continue; }
+ const Triangle tb0 = tB[stMap[j][0]];
+ const Triangle tb1 = tB[stMap[j][1]];
+ if ( !(do_intersect(sa,tb0) && do_intersect(sa,tb1) ) ) { continue; }
+ CGAL::Object osb1 = intersection(sb,ta0);
+ CGAL::Object osb2 = intersection(sb,ta1);
+ CGAL::Object osa1 = intersection(sa,tb0);
+ CGAL::Object osa2 = intersection(sa,tb1);
+ const Point* psb1 = CGAL::object_cast<Point>(&osb1);
+ const Point* psb2 = CGAL::object_cast<Point>(&osb2);
+ const Point* psa1 = CGAL::object_cast<Point>(&osa1);
+ const Point* psa2 = CGAL::object_cast<Point>(&osa2);
+ if ( !(psb1 && psb2 && psa1 && psa2) ) { continue; }
+ Vector_3 n = CGAL::cross_product(sa.to_vector(),sb.to_vector());
+ Vector3r nApprox;
+ for (int k=0; k<3; k++) {
+ normal[k] = n[k];
+ #define OP(p) p->operator[](k)
+ nApprox[k] = .5*(OP(psa1)+OP(psa2)) - .5*(OP(psb1)+OP(psb2));
+ contactPoint[k] = .25*(OP(psa1)+OP(psa2)+OP(psb1)+OP(psb2));
+ #undef OP
+ }
+ if ( nApprox.dot(normal) < 0 ) { normal *= (Real)-1.; }
+ normal.normalize();
+ const Point* p[4] = {psb1,psb2,psa1,psa2};
+ penetrationVolume = TetrahedronVolume(p);
+ return true;
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkEdgeToTriangleCase1( // edge smaller than triangle
+ const Triangle tA[4],
+ const Segment sB[6],
+ const Point pB[6],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+ for (int i=0; i<4; i++) {
+ const Triangle& ta = tA[i];
+ int ni = 0;
+ for (int j=0; j<6; j++) {
+ const Segment& s = sB[j];
+ if ( do_intersect(ta,s) ) { ni++; }
+ }
+ if ( ni != 4 ) { continue; }
+ Vector_3 n = CGAL::normal(ta[0],ta[1],ta[2]);
+ for (int j=0; j<3; j++) {
+ const Point& p1 = pB[j];
+ if ( Vector_3(ta[0],p1)*n > 0.) { continue; }
+ const Segment& s10 = sB[psMap[j][0]];
+ const Segment& s11 = sB[psMap[j][1]];
+ const Segment& s12 = sB[psMap[j][2]];
+ bool b10 = do_intersect(ta,s10);
+ bool b11 = do_intersect(ta,s11);
+ bool b12 = do_intersect(ta,s12);
+ if ( !((b10 && b11) || (b11 && b12) || (b12 && b10) ) ) { continue; }
+ for (int k=j+1; k<4; k++) {
+ const Point& p2 = pB[k];
+ if ( Vector_3(ta[0],p2)*n > 0. ) { continue; }
+ const Segment& s20 = sB[psMap[k][0]];
+ const Segment& s21 = sB[psMap[k][1]];
+ const Segment& s22 = sB[psMap[k][2]];
+ bool b20 = do_intersect(ta,s20);
+ bool b21 = do_intersect(ta,s21);
+ bool b22 = do_intersect(ta,s22);
+ if ( !((b20 && b21) || (b21 && b22) || (b22 && b20) ) ) { continue; }
+ int si = ppsMap[j][k];
+ const Segment& sb = sB[si];
+ int l,m;
+ for (l=0; l<3; l++) {
+ if (l!=j && l!=k) { break; }
+ }
+ for (m=l+1; m<4; m++) {
+ if (m!=j && m!=k) { break; }
+ }
+ const Point& p3 = pB[l];
+ const Point& p4 = pB[m];
+ const Segment& s13 = sB[ppsMap[j][l]];
+ const Segment& s14 = sB[ppsMap[j][m]];
+ const Segment& s23 = sB[ppsMap[k][l]];
+ const Segment& s24 = sB[ppsMap[k][m]];
+ CGAL::Object o13 = intersection(ta,s13);
+ CGAL::Object o14 = intersection(ta,s14);
+ CGAL::Object o23 = intersection(ta,s23);
+ CGAL::Object o24 = intersection(ta,s24);
+ const Point* ps13 = CGAL::object_cast<Point>(&o13);
+ const Point* ps14 = CGAL::object_cast<Point>(&o14);
+ const Point* ps23 = CGAL::object_cast<Point>(&o23);
+ const Point* ps24 = CGAL::object_cast<Point>(&o24);
+ if ( !(ps13 && ps14 && ps23 && ps24) ) { continue; }
+ const Point* pp1[4] = {&p1,&p2,ps13,ps14};
+ const Point* pp2[4] = {&p2,ps23,ps24,ps14};
+ const Point* pp3[4] = {&p2,ps23,ps13,ps14};
+ Real v1 = TetrahedronVolume(pp1);
+ Real v2 = TetrahedronVolume(pp2);
+ Real v3 = TetrahedronVolume(pp3);
+ Vector3r cg1,cg2,cg3;
+ for (l=0; l<3; l++) {
+ normal[l] = n[l];
+ #define OP(p) p->operator[](l)
+ cg1[l] = .25*(p1[l]+p2[l]+OP(ps13)+OP(ps14));
+ cg2[l] = .25*(p2[l]+OP(ps23)+OP(ps24)+OP(ps14));
+ cg3[l] = .25*(p2[l]+OP(ps23)+OP(ps13)+OP(ps14));
+ #undef OP
+ }
+ penetrationVolume = v1 + v2 + v3;
+ contactPoint = (v1*cg1 + v2*cg2 + v3*cg3) / penetrationVolume;
+ normal.normalize();
+ return true;
+ }
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkEdgeToTriangleCase2( // edge larger than triangle
+ const Triangle tA[4],
+ const Triangle tB[4],
+ const Segment sA[6],
+ const Segment sB[6],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+
+ for (int i=0; i<6; i++) {
+ const Segment& sb = sB[i];
+ int ni = 0;
+ for (int j=0; j<4; j++) {
+ const Triangle& t = tA[j];
+ if ( do_intersect(t,sb) ) { ni++; }
+ }
+ if ( ni != 2 ) { continue; }
+ for (int j=0; j<3; j++) {
+ const Triangle& ta1 = tA[j];
+ if ( !(do_intersect(ta1,sb) ) ) { continue; }
+ for (int k=j+1; k<4; k++) {
+ const Triangle& ta2 = tA[k];
+ if ( !(do_intersect(ta2,sb) ) ) { continue; }
+ const Triangle& tb1 = tB[stMap[i][0]];
+ const Triangle& tb2 = tB[stMap[i][1]];
+ const Segment& sa1a = sA[tsMap[j][0]];
+ const Segment& sa1b = sA[tsMap[j][1]];
+ const Segment& sa1c = sA[tsMap[j][2]];
+ bool b1a = do_intersect(sa1a,tb1) && do_intersect(sa1a,tb2);
+ bool b1b = do_intersect(sa1b,tb1) && do_intersect(sa1b,tb2);
+ bool b1c = do_intersect(sa1c,tb1) && do_intersect(sa1c,tb2);
+ if ( !(b1a || b1b || b1c) ) { continue; }
+ const Segment& sa2a = sA[tsMap[k][0]];
+ const Segment& sa2b = sA[tsMap[k][1]];
+ const Segment& sa2c = sA[tsMap[k][2]];
+ bool b2a = do_intersect(sa2a,tb1) && do_intersect(sa2a,tb2);
+ bool b2b = do_intersect(sa2b,tb1) && do_intersect(sa2b,tb2);
+ bool b2c = do_intersect(sa2c,tb1) && do_intersect(sa2c,tb2);
+ if ( !(b2a || b2b || b2c) ) { continue; }
+ int l = b1a? tsMap[j][0] : b1b? tsMap[j][1] : tsMap[j][2];
+ int m = b2a? tsMap[k][0] : b2b? tsMap[k][1] : tsMap[k][2];
+ const Segment& sa1 = sA[l];
+ const Segment& sa2 = sA[m];
+ const Triangle& taN = tA[sstMap[l][m]];
+ CGAL::Object o1 = intersection(sb,ta1);
+ CGAL::Object o2 = intersection(sb,ta2);
+ CGAL::Object o11 = intersection(sa1,tb1);
+ CGAL::Object o12 = intersection(sa1,tb2);
+ CGAL::Object o21 = intersection(sa2,tb1);
+ CGAL::Object o22 = intersection(sa2,tb2);
+ const Point* p1 = CGAL::object_cast<Point>(&o1);
+ const Point* p2 = CGAL::object_cast<Point>(&o2);
+ const Point* p11 = CGAL::object_cast<Point>(&o11);
+ const Point* p12 = CGAL::object_cast<Point>(&o12);
+ const Point* p21 = CGAL::object_cast<Point>(&o21);
+ const Point* p22 = CGAL::object_cast<Point>(&o22);
+ if ( !(p1 && p2 && p11 && p12 && p21 && p22) ) { continue; }
+ const Point* pp1[4] = {p1,p2,p11,p12};
+ const Point* pp2[4] = {p2,p21,p22,p12};
+ const Point* pp3[4] = {p2,p21,p11,p12};
+ Real v1 = TetrahedronVolume(pp1);
+ Real v2 = TetrahedronVolume(pp2);
+ Real v3 = TetrahedronVolume(pp3);
+ Vector3r cg1,cg2,cg3;
+ Vector_3 n = CGAL::normal(taN[0],taN[1],taN[2]);
+ for (int l=0; l<3; l++) {
+ normal[l] = n[l];
+ #define OP(p) p->operator[](l)
+ cg1[l] = .25*(OP(p1)+OP(p2)+OP(p11)+OP(p12));
+ cg2[l] = .25*(OP(p2)+OP(p21)+OP(p22)+OP(p12));
+ cg3[l] = .25*(OP(p2)+OP(p21)+OP(p11)+OP(p12));
+ #undef OP
+ }
+ penetrationVolume = v1 + v2 + v3;
+ contactPoint = (v1*cg1 + v2*cg2 + v3*cg3) / penetrationVolume;
+ normal.normalize();
+ return true;
+ }
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkVertexToEdgeCase(
+ const Point pA[4],
+ const Segment sA[6],
+ const Triangle tA[4],
+ const Segment sB[6],
+ const Triangle tB[4],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+ for (int i=0; i<4; i++) {
+ const Point& pa = pA[i];
+ if ( Vector_3(tB[0][0],pa)*CGAL::normal(tB[0][0],tB[0][1],tB[0][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[1][0],pa)*CGAL::normal(tB[1][0],tB[1][1],tB[1][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[2][0],pa)*CGAL::normal(tB[2][0],tB[2][1],tB[2][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[3][0],pa)*CGAL::normal(tB[3][0],tB[3][1],tB[3][2]) > 0. ) { continue; }
+ const Segment& sa1 = sA[psMap[i][0]];
+ const Segment& sa2 = sA[psMap[i][1]];
+ const Segment& sa3 = sA[psMap[i][2]];
+ for (int j=0; j<6; j++) {
+ const Segment& sb = sB[j];
+ const Triangle& tb1 = tB[stMap[j][0]];
+ const Triangle& tb2 = tB[stMap[j][1]];
+ const Triangle& ta1 = tA[ptMap[i][0]];
+ const Triangle& ta2 = tA[ptMap[i][1]];
+ const Triangle& ta3 = tA[ptMap[i][2]];
+ bool bsa1tb1 = do_intersect(sa1,tb1);
+ bool bsa1tb2 = do_intersect(sa1,tb2);
+ bool bsa2tb1 = do_intersect(sa2,tb1);
+ bool bsa2tb2 = do_intersect(sa2,tb2);
+ bool bsa3tb1 = do_intersect(sa3,tb1);
+ bool bsa3tb2 = do_intersect(sa3,tb2);
+ bool bsbta1 = do_intersect(sb,ta1);
+ bool bsbta2 = do_intersect(sb,ta2);
+ bool bsbta3 = do_intersect(sb,ta3);
+ if ( !( (bsa1tb1 || bsa1tb2) && (bsa2tb1 || bsa2tb2) && (bsa3tb1 || bsa3tb2) && ((bsbta1 && bsbta2) || (bsbta2 && bsbta3) || (bsbta3 && bsbta1)) ) ) { continue; }
+ CGAL::Object oa1 = intersection(sa1,bsa1tb1? tb1 : tb2);
+ CGAL::Object oa2 = intersection(sa2,bsa2tb1? tb1 : tb2);
+ CGAL::Object oa3 = intersection(sa3,bsa3tb1? tb1 : tb2);
+ CGAL::Object ob1 = intersection(sb, (bsbta1 && bsbta2)? ta1 : (bsbta2 && bsbta3)? ta2 : ta3);
+ CGAL::Object ob2 = intersection(sb, (bsbta1 && bsbta2)? ta2 : (bsbta2 && bsbta3)? ta3 : ta1);
+ const Point* pa1 = CGAL::object_cast<Point>(&oa1);
+ const Point* pa2 = CGAL::object_cast<Point>(&oa2);
+ const Point* pa3 = CGAL::object_cast<Point>(&oa3);
+ const Point* pb1 = CGAL::object_cast<Point>(&ob1);
+ const Point* pb2 = CGAL::object_cast<Point>(&ob2);
+ if ( !(pa1 && pa2 && pa3 && pb1 && pb2) ) { continue; }
+ Segment sa(*pa1,*pa2);
+ Real d1 = sqrt(CGAL::squared_distance(sa,*pb1));
+ Real d2 = sqrt(CGAL::squared_distance(sa,*pb2));
+ const Point* ppb1 = d1<d2? pb1 : pb2;
+ const Point* ppb2 = d1<d2? pb2 : pb1;
+ const Point* pp1[4] = {&pa,pa1,pa2,pa3};
+ const Point* pp2[4] = {pa1,pa2,pa3,ppb2};
+ const Point* pp3[4] = {pa1,pa2,ppb1,ppb2};
+ Real v1 = TetrahedronVolume(pp1);
+ Real v2 = TetrahedronVolume(pp2);
+ Real v3 = TetrahedronVolume(pp3);
+ Vector3r cg1,cg2,cg3;
+ Vector_3 n(pa,sb.supporting_line().projection(pa));
+ for (int l=0; l<3; l++) {
+ normal[l] = n[l];
+ #define OP(p) p->operator[](l)
+ cg1[l] = .25*(pa[l]+OP(pa1)+OP(pa2)+OP(pa3));
+ cg2[l] = .25*(OP(pa1)+OP(pa2)+OP(pa3)+OP(ppb2));
+ cg3[l] = .25*(OP(pa1)+OP(pa2)+OP(ppb1)+OP(ppb2));
+ #undef OP
+ }
+ penetrationVolume = v1 + v2 + v3;
+ contactPoint = (v1*cg1 + v2*cg2 + v3*cg3) / penetrationVolume;
+ normal.normalize();
+ return true;
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::checkVertexToVertexCase(
+ const Point pA[4],
+ const Point pB[4],
+ const Segment sA[6],
+ const Triangle tA[4],
+ const Triangle tB[4],
+ Vector3r& normal,
+ Vector3r& contactPoint,
+ Real& penetrationVolume)
+{
+ for (int i=0; i<4; i++) {
+ const Point& pa = pA[i];
+ if ( Vector_3(tB[0][0],pa)*CGAL::normal(tB[0][0],tB[0][1],tB[0][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[1][0],pa)*CGAL::normal(tB[1][0],tB[1][1],tB[1][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[2][0],pa)*CGAL::normal(tB[2][0],tB[2][1],tB[2][2]) > 0. ) { continue; }
+ if ( Vector_3(tB[3][0],pa)*CGAL::normal(tB[3][0],tB[3][1],tB[3][2]) > 0. ) { continue; }
+ const Segment& sa1 = sA[psMap[i][0]];
+ const Segment& sa2 = sA[psMap[i][1]];
+ const Segment& sa3 = sA[psMap[i][2]];
+ for (int j=0; j<4; j++) {
+ const Point& pb = pB[j];
+ if ( Vector_3(tA[0][0],pb)*CGAL::normal(tA[0][0],tA[0][1],tA[0][2]) > 0. ) { continue; }
+ if ( Vector_3(tA[1][0],pb)*CGAL::normal(tA[1][0],tA[1][1],tA[1][2]) > 0. ) { continue; }
+ if ( Vector_3(tA[2][0],pb)*CGAL::normal(tA[2][0],tA[2][1],tA[2][2]) > 0. ) { continue; }
+ if ( Vector_3(tA[3][0],pb)*CGAL::normal(tA[3][0],tB[3][1],tB[3][2]) > 0. ) { continue; }
+ const Triangle& tb1 = tB[ptMap[j][0]];
+ const Triangle& tb2 = tB[ptMap[j][1]];
+ const Triangle& tb3 = tB[ptMap[j][2]];
+ bool b11 = do_intersect(sa1,tb1);
+ bool b12 = do_intersect(sa1,tb2);
+ bool b13 = do_intersect(sa1,tb3);
+ bool b21 = do_intersect(sa2,tb1);
+ bool b22 = do_intersect(sa2,tb2);
+ bool b23 = do_intersect(sa2,tb3);
+ bool b31 = do_intersect(sa3,tb1);
+ bool b32 = do_intersect(sa3,tb2);
+ bool b33 = do_intersect(sa3,tb3);
+ if ( !(b11 || b12 || b13) && (b21 || b22 || b23) && (b31 || b32 || b33) ) { continue; }
+ CGAL::Object o1 = intersection(sa1, b11? tb1: b12? tb2 : tb3);
+ CGAL::Object o2 = intersection(sa2, b21? tb1: b22? tb2 : tb3);
+ CGAL::Object o3 = intersection(sa3, b31? tb1: b32? tb2 : tb3);
+ const Point* p1 = CGAL::object_cast<Point>(&o1);
+ const Point* p2 = CGAL::object_cast<Point>(&o2);
+ const Point* p3 = CGAL::object_cast<Point>(&o3);
+ if ( !(p1 && p2 && p3) ) { continue; }
+ const Point* pp1[4] = {&pa,p1,p2,p3};
+ const Point* pp2[4] = {&pb,p2,p3,p3};
+ Real v1 = TetrahedronVolume(pp1);
+ Real v2 = TetrahedronVolume(pp2);
+ Vector3r cg1,cg2;
+ Vector_3 n(pa,pb);
+ for (int l=0; l<3; l++) {
+ normal[l] = n[l];
+ #define OP(p) p->operator[](l)
+ cg1[l] = .25*(pa[l]+OP(p1)+OP(p2)+OP(p3));
+ cg2[l] = .25*(pb[l]+OP(p1)+OP(p2)+OP(p3));
+ #undef OP
+ }
+ penetrationVolume = v1 + v2;
+ contactPoint = (v1*cg1 + v2*cg2) / penetrationVolume;
+ normal.normalize();
+ return true;
+ }
+ }
+ return false;
+}
+
+bool Ig2_Tetra_Tetra_TTetraSimpleGeom::go(
+ const shared_ptr<Shape>& cm1,
+ const shared_ptr<Shape>& cm2,
+ const State& state1,
+ const State& state2,
+ const Vector3r& shift2,
+ const bool& force,
+ const shared_ptr<Interaction>& interaction)
+{
+ const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
+ Tetra* shape1 = static_cast<Tetra*>(cm1.get());
+ Tetra* shape2 = static_cast<Tetra*>(cm2.get());
+
+ Point p1[4], p2[4];
+ Vector3r vTemp;
+ // vertices in global coordinates
+ for (int i=0; i<4; i++) {
+ vTemp = se31.position + se31.orientation*shape1->v[i];
+ p1[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
+ vTemp = se32.position + se32.orientation*shape2->v[i];
+ p2[i] = Point(vTemp[0],vTemp[1],vTemp[2]);
+ }
+
+ // Faces (CGAL triangles) of each tetra
+ #define T(p,i,j,k) Triangle(p[i],p[j],p[k])
+ const Triangle t1[4] = {T(p1,0,1,3),T(p1,0,2,1),T(p1,1,2,3),T(p1,0,3,2)};
+ const Triangle t2[4] = {T(p2,0,1,3),T(p2,0,2,1),T(p2,1,2,3),T(p2,0,3,2)};
+ #undef T
+ // Edges (CGAL segments) of each tetra
+ #define S(p,i,j) Segment(p[i],p[j])
+ const Segment s1[6] = {S(p1,0,1),S(p1,1,2),S(p1,0,2),S(p1,0,3),S(p1,1,3),S(p1,2,3)};
+ const Segment s2[6] = {S(p2,0,1),S(p2,1,2),S(p2,0,2),S(p2,0,3),S(p2,1,3),S(p2,2,3)};
+ #undef S
+
+ Vector3r n;
+ Vector3r cp;
+ Real V;
+ int flag;
+
+ # define SET_GEOM_AND_RETURN_TRUE \
+ shared_ptr<TTetraSimpleGeom> geom; \
+ if (!interaction->geom) geom=shared_ptr<TTetraSimpleGeom>(new TTetraSimpleGeom()); \
+ else geom=YADE_PTR_CAST<TTetraSimpleGeom>(interaction->geom); \
+ interaction->geom=geom; \
+ geom->normal = n; \
+ geom->contactPoint = cp; \
+ geom->penetrationVolume = V; \
+ geom->flag = flag; \
+ return true;
+
+ if (checkVertexToTriangleCase(t1,p2,s2,n,cp,V)) {
+ flag = 1;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkVertexToTriangleCase(t2,p1,s1,n,cp,V)) {
+ n *= -1.;
+ flag = 2;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkEdgeToEdgeCase(s1,s2,t1,t2,n,cp,V)) {
+ flag = 3;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkEdgeToTriangleCase1(t1,s2,p2,n,cp,V)) {
+ flag = 4;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkEdgeToTriangleCase1(t2,s1,p1,n,cp,V)) {
+ n *= -1.;
+ flag = 5;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkEdgeToTriangleCase2(t1,t2,s1,s2,n,cp,V)) {
+ flag = 6;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkEdgeToTriangleCase2(t2,t1,s2,s1,n,cp,V)) {
+ n *= -1.;
+ flag = 7;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkVertexToEdgeCase(p1,s1,t1,s2,t2,n,cp,V)) {
+ n *= -1.;
+ flag = 8;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+ if (checkVertexToEdgeCase(p2,s2,t2,s1,t1,n,cp,V)) {
+ flag = 9;
+ SET_GEOM_AND_RETURN_TRUE
+ }
+
+ #undef SET_GEOM_AND_RETURN_TRUE
+
+ if (interaction->geom) {
+ TTetraSimpleGeom* geom = static_cast<TTetraSimpleGeom*>(interaction->geom.get());
+ geom->penetrationVolume = (Real)-1.;
+ geom->flag = 0;
+ return true;
+ }
+ return false;
+}
+#endif
+
+
+
+
+
+
+
CREATE_LOGGER(Ig2_Tetra_Tetra_TTetraGeom);
/*! Calculate configuration of Tetra - Tetra intersection.
@@ -107,7 +682,9 @@
tA=Tetra(Vector3r(0,0,0),Vector3r(1,0,0),Vector3r(0,1,0),Vector3r(0,0,1));
#endif
list<Tetra> tAB=Tetra2TetraIntersection(tA,tB);
- if(tAB.size()==0) { /* LOG_DEBUG("No intersection."); */ return false;} //no intersecting volume
+ if (!interaction->isReal() && !force) {
+ if(tAB.size()==0) { /* LOG_DEBUG("No intersection."); */ return false;} //no intersecting volume
+ }
Real V(0); // volume of intersection (cummulative)
Vector3r Sg(0,0,0); // static moment of intersection
@@ -186,6 +763,8 @@
Real maxPenetrationDepthB=sqrt(6*(IB(ix,ix)+IB(ixx,ixx)-IB(ixxx,ixxx))/V);
TRVAR2(maxPenetrationDepthA,maxPenetrationDepthB);
+ //normal = se32.position - se31.position; normal.normalize();
+
/* store calculated stuff in bang; some is redundant */
bang->normal=normal;
bang->equivalentCrossSection=equivalentCrossSection;
@@ -367,6 +946,7 @@
return(ret); // prevent warning
}
+
CREATE_LOGGER(TetraVolumetricLaw);
/*! Apply forces on tetrahedra in collision based on geometric configuration provided by Ig2_Tetra_Tetra_TTetraGeom.
@@ -413,12 +993,13 @@
#ifdef YADE_OPENGL
#include<yade/lib/opengl/OpenGLWrapper.hpp>
- void Gl1_Tetra::go(const shared_ptr<Shape>& cm, const shared_ptr<State>&,bool,const GLViewInfo&)
+ bool Gl1_Tetra::wire;
+ void Gl1_Tetra::go(const shared_ptr<Shape>& cm, const shared_ptr<State>&,bool wire2,const GLViewInfo&)
{
glMaterialv(GL_FRONT,GL_AMBIENT_AND_DIFFUSE,Vector3r(cm->color[0],cm->color[1],cm->color[2]));
glColor3v(cm->color);
Tetra* t=static_cast<Tetra*>(cm.get());
- if (0) { // wireframe, as for Tetrahedron
+ if (wire && wire2) { // wireframe, as for Tetrahedron
glDisable(GL_LIGHTING);
glBegin(GL_LINES);
#define __ONEWIRE(a,b) glVertex3v(t->v[a]);glVertex3v(t->v[b])
@@ -432,10 +1013,10 @@
glDisable(GL_CULL_FACE); glEnable(GL_LIGHTING);
glBegin(GL_TRIANGLES);
#define __ONEFACE(a,b,c) n=(t->v[b]-t->v[a]).cross(t->v[c]-t->v[a]); n.normalize(); faceCenter=(t->v[a]+t->v[b]+t->v[c])/3.; if((faceCenter-center).dot(n)<0)n=-n; glNormal3v(n); glVertex3v(t->v[a]); glVertex3v(t->v[b]); glVertex3v(t->v[c]);
- __ONEFACE(3,0,1);
- __ONEFACE(0,1,2);
+ __ONEFACE(0,2,1);
+ __ONEFACE(0,1,3);
__ONEFACE(1,2,3);
- __ONEFACE(2,3,0);
+ __ONEFACE(0,3,2);
#undef __ONEFACE
glEnd();
}
@@ -599,6 +1180,47 @@
}
-Real TetrahedronVolume(const Vector3r v[4]) { return fabs((Vector3r(v[3])-Vector3r(v[0])).dot((Vector3r(v[3])-Vector3r(v[1])).cross(Vector3r(v[3])-Vector3r(v[2]))))/6.; }
-Real TetrahedronVolume(const vector<Vector3r>& v) { return fabs(Vector3r(v[1]-v[0]).dot(Vector3r(v[2]-v[0]).cross(v[3]-v[0])))/6.; }
+
+
+Real TetrahedronSignedVolume(const Vector3r v[4]) { return (Vector3r(v[3])-Vector3r(v[0])).dot((Vector3r(v[3])-Vector3r(v[1])).cross(Vector3r(v[3])-Vector3r(v[2])))/6.; }
+Real TetrahedronVolume(const Vector3r v[4]) { return fabs(TetrahedronSignedVolume(v)); }
+Real TetrahedronSignedVolume(const vector<Vector3r>& v) { return Vector3r(v[1]-v[0]).dot(Vector3r(v[2]-v[0]).cross(v[3]-v[0]))/6.; }
+Real TetrahedronVolume(const vector<Vector3r>& v) { return fabs(TetrahedronSignedVolume(v)); }
+Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> >* v[4]) {
+ Vector3r vv[4];
+ for (int i=0; i<4; i++) {
+ for (int j=0; j<3; j++) {
+ vv[i][j] = v[i]->operator[](j);
+ }
+ }
+ return TetrahedronVolume(vv);
+}
+Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> > v[4]) {
+ Vector3r vv[4];
+ for (int i=0; i<4; i++) {
+ for (int j=0; j<3; j++) {
+ vv[i][j] = v[i][j];
+ }
+ }
+ return TetrahedronVolume(vv);
+}
+
+
+
+void Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+ int id1 = contact->getId1(), id2 = contact->getId2();
+ TTetraSimpleGeom* geom= static_cast<TTetraSimpleGeom*>(ig.get());
+ NormPhys* phys = static_cast<NormPhys*>(ip.get());
+ if ( geom->flag == 0 || geom->penetrationVolume <= 0. ) {
+ scene->interactions->requestErase(contact);
+ return;
+ }
+ Real& un=geom->penetrationVolume;
+ phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
+
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
+ applyForceAtContactPoint(-phys->normalForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
+ // TODO periodic
+}
=== modified file 'pkg/dem/Tetra.hpp'
--- pkg/dem/Tetra.hpp 2011-01-09 16:34:50 +0000
+++ pkg/dem/Tetra.hpp 2013-09-08 11:12:42 +0000
@@ -1,4 +1,5 @@
// © 2007 Václav Šmilauer <eudoxos@xxxxxxxx>
+// © 2013 Jan Stránský <jan.stransky@xxxxxxxxxxx>
#pragma once
@@ -10,6 +11,12 @@
#include<yade/pkg/common/Aabb.hpp>
#include<yade/pkg/common/Dispatching.hpp>
+#include<yade/pkg/common/NormShearPhys.hpp>
+
+#ifdef YADE_CGAL
+ #include <CGAL/Cartesian.h>
+ //#include <CGAL/intersections.h>
+#endif
/* Our mold of tetrahedron: just 4 vertices.
@@ -19,16 +26,16 @@
public:
Tetra(Vector3r v0, Vector3r v1, Vector3r v2, Vector3r v3) { createIndex(); v.resize(4); v[0]=v0; v[1]=v1; v[2]=v2; v[3]=v3; }
virtual ~Tetra();
- protected:
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Tetra,Shape,"Tetrahedron geometry.",
- ((std::vector<Vector3r>,v,std::vector<Vector3r>(4),,"Tetrahedron vertices in global coordinate system.")),
- /*ctor*/createIndex();
- );
- REGISTER_CLASS_INDEX(Tetra,Shape);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Tetra,Shape,"Tetrahedron geometry.",
+ ((std::vector<Vector3r>,v,std::vector<Vector3r>(4),,"Tetrahedron vertices (in local coordinate system).")),
+ /*ctor*/createIndex();
+ );
+ REGISTER_CLASS_INDEX(Tetra,Shape);
};
REGISTER_SERIALIZABLE(Tetra);
+
/*! Collision configuration for Tetra and something.
* This is expressed as penetration volume properties: centroid, volume, orientation of principal axes, inertia.
*
@@ -36,22 +43,37 @@
class TTetraGeom: public IGeom{
public:
virtual ~TTetraGeom();
- protected:
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(TTetraGeom,IGeom,"Geometry of interaction between 2 :yref:`tetrahedra<Tetra>`, including volumetric characteristics",
- ((Real,penetrationVolume,NaN,,"Volume of overlap [m³]"))
- ((Real,equivalentCrossSection,NaN,,"Cross-section of the overlap (perpendicular to the axis of least inertia"))
- ((Real,maxPenetrationDepthA,NaN,,"??"))
- ((Real,maxPenetrationDepthB,NaN,,"??"))
- ((Real,equivalentPenetrationDepth,NaN,,"??"))
- ((Vector3r,contactPoint,,,"Contact point (global coords)"))
- ((Vector3r,normal,,,"Normal of the interaction, directed in the sense of least inertia of the overlap volume")),
- createIndex();
- );
- //FUNCTOR2D(Tetra,Tetra);
- REGISTER_CLASS_INDEX(TTetraGeom,IGeom);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(TTetraGeom,IGeom,"Geometry of interaction between 2 :yref:`tetrahedra<Tetra>`, including volumetric characteristics",
+ ((Real,penetrationVolume,NaN,,"Volume of overlap [m³]"))
+ ((Real,equivalentCrossSection,NaN,,"Cross-section of the overlap (perpendicular to the axis of least inertia"))
+ ((Real,maxPenetrationDepthA,NaN,,"??"))
+ ((Real,maxPenetrationDepthB,NaN,,"??"))
+ ((Real,equivalentPenetrationDepth,NaN,,"??"))
+ ((Vector3r,contactPoint,,,"Contact point (global coords)"))
+ ((Vector3r,normal,,,"Normal of the interaction, directed in the sense of least inertia of the overlap volume")),
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(TTetraGeom,IGeom);
};
REGISTER_SERIALIZABLE(TTetraGeom);
+
+class TTetraSimpleGeom: public IGeom{
+ public:
+ virtual ~TTetraSimpleGeom();
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(TTetraSimpleGeom,IGeom,"EXPERIMENTAL. Geometry of interaction between 2 :yref:`tetrahedra<Tetra>`",
+ ((Real,penetrationVolume,NaN,,"Volume of overlap [m³]"))
+ ((Vector3r,contactPoint,,,"Contact point (global coords)"))
+ ((Vector3r,normal,,,"Normal of the interaction TODO"))
+ ((int,flag,0,,"TODO")),
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(TTetraSimpleGeom,IGeom);
+};
+REGISTER_SERIALIZABLE(TTetraSimpleGeom);
+
+
+
/*! Creates Aabb from Tetra.
*
* Self-contained. */
@@ -69,7 +91,9 @@
class Gl1_Tetra: public GlShapeFunctor{
public:
virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
- YADE_CLASS_BASE_DOC(Gl1_Tetra,GlShapeFunctor,"Renders :yref:`Tetra` object");
+ YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Tetra,GlShapeFunctor,"Renders :yref:`Tetra` object",
+ ((bool,wire,true,,"TODO"))
+ );
RENDERS(Tetra);
};
REGISTER_SERIALIZABLE(Gl1_Tetra);
@@ -108,11 +132,67 @@
REGISTER_SERIALIZABLE(Ig2_Tetra_Tetra_TTetraGeom);
+
+
+#ifdef YADE_CGAL
+class Ig2_Tetra_Tetra_TTetraSimpleGeom: public IGeomFunctor
+{
+ protected:
+ typedef CGAL::Cartesian<Real> K;
+ typedef K::Point_3 Point;
+ typedef CGAL::Tetrahedron_3<K> CGALTetra;
+ typedef CGAL::Segment_3<K> Segment;
+ typedef CGAL::Triangle_3<K> Triangle;
+ typedef CGAL::Vector_3<K> Vector_3;
+ bool checkVertexToTriangleCase(const Triangle tA[4], const Point pB[4], const Segment sB[6], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ bool checkEdgeToEdgeCase(const Segment sA[6], const Segment sB[6], const Triangle tA[4], const Triangle tB[4], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ bool checkEdgeToTriangleCase1(const Triangle tA[4], const Segment sB[6], const Point pB[4], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ bool checkEdgeToTriangleCase2(const Triangle tA[4], const Triangle tB[4], const Segment sA[6], const Segment sB[6], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ bool checkVertexToEdgeCase(const Point pA[4], const Segment sA[6], const Triangle tA[4], const Segment sB[6], const Triangle tB[4], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ bool checkVertexToVertexCase(const Point pA[4], const Point pB[4], const Segment sA[6], const Triangle tA[4], const Triangle tB[4], Vector3r& normal, Vector3r& contactPoint, Real& penetrationVolume);
+ static const int psMap[4][3];
+ static const int ptMap[4][3];
+ static const int stMap[6][2];
+ static const int tsMap[4][3];
+ static const int ppsMap[4][4];
+ static const int sstMap[6][6];
+ public:
+ virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+ virtual bool goReverse( const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw std::logic_error("Ig2_Tetra_Tetra_TTetraSimpleGeom::goReverse called, but the functor is symmetric."); }
+ FUNCTOR2D(Tetra,Tetra);
+ DEFINE_FUNCTOR_ORDER_2D(Tetra,Tetra);
+ YADE_CLASS_BASE_DOC(Ig2_Tetra_Tetra_TTetraSimpleGeom,IGeomFunctor,"EXPERIMANTAL. Create/update geometry of collision between 2 :yref:`tetrahedra<Tetra>` (:yref:`TTetraSimpleGeom` instance)");
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Ig2_Tetra_Tetra_TTetraSimpleGeom);
+#endif
+
+
+
+
+class Law2_TTetraSimpleGeom_NormPhys_Simple: public LawFunctor{
+ public:
+ virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ YADE_CLASS_BASE_DOC(Law2_TTetraSimpleGeom_NormPhys_Simple,LawFunctor,"EXPERIMENTAL. TODO");
+ FUNCTOR2D(TTetraSimpleGeom,NormPhys);
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_TTetraSimpleGeom_NormPhys_Simple);
+
+
+
+
// Miscillaneous functions
//! Tetrahedron's volume.
/// http://en.wikipedia.org/wiki/Tetrahedron#Surface_area_and_volume
+Real TetrahedronSignedVolume(const Vector3r v[4]);
Real TetrahedronVolume(const Vector3r v[4]);
+Real TetrahedronSignedVolume(const vector<Vector3r>& v);
Real TetrahedronVolume(const vector<Vector3r>& v);
+#ifdef YADE_CGAL
+Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> >* v[4]);
+Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> > v[4]);
+#endif
Matrix3r TetrahedronInertiaTensor(const vector<Vector3r>& v);
//Matrix3r TetrahedronInertiaTensor(const Vector3r v[4]);
Matrix3r TetrahedronCentralInertiaTensor(const vector<Vector3r>& v);
=== modified file 'py/_utils.cpp'
--- py/_utils.cpp 2013-05-31 13:27:13 +0000
+++ py/_utils.cpp 2013-09-08 11:12:42 +0000
@@ -7,6 +7,7 @@
#include<yade/pkg/dem/ScGeom.hpp>
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/pkg/common/Facet.hpp>
+#include<yade/pkg/dem/Tetra.hpp>
#include<yade/pkg/common/Sphere.hpp>
#include<yade/pkg/common/NormShearPhys.hpp>
#include<yade/lib/computational-geometry/Hull2d.hpp>
@@ -543,4 +544,11 @@
py::def("growParticles",Shop::growParticles,(py::args("multiplier"), py::args("updateMass")=true, py::args("dynamicOnly")=true), "Change the size of spheres and sphere clumps by the multiplier. If updateMass=True, then the mass is updated. dynamicOnly=True is mandatory in many cases since in current implementation the function would crash on the non-spherical and non-dynamic bodies (e.g. facets, boxes, etc.)");
py::def("intrsOfEachBody",intrsOfEachBody,"returns list of lists of interactions of each body");
py::def("numIntrsOfEachBody",numIntrsOfEachBody,"returns list of number of interactions of each body");
+ py::def("TetrahedronSignedVolume",static_cast<Real (*)(const vector<Vector3r>&)>(&TetrahedronSignedVolume),"TODO");
+ py::def("TetrahedronVolume",static_cast<Real (*)(const vector<Vector3r>&)>(&TetrahedronVolume),"TODO");
+ py::def("TetrahedronInertiaTensor",TetrahedronInertiaTensor,"TODO");
+ py::def("TetrahedronCentralInertiaTensor",TetrahedronCentralInertiaTensor,"TODO");
+ py::def("TetrahedronWithLocalAxesPrincipal",TetrahedronWithLocalAxesPrincipal,"TODO");
+ py::def("momentum",Shop::momentum,"TODO");
+ py::def("angularMomentum",Shop::angularMomentum,(py::args("origin")=Vector3r(Vector3r::Zero())),"TODO");
}
=== modified file 'py/utils.py'
--- py/utils.py 2013-08-14 09:03:28 +0000
+++ py/utils.py 2013-09-08 11:12:42 +0000
@@ -348,6 +348,37 @@
b.chain=chain
return b
+def tetra(vertices,strictCheck=True,dynamic=True,fixed=False,wire=True,color=None,highlight=False,noBound=False,material=-1,mask=1,chain=-1):
+ """Create tetrahedron with given parameters.
+
+ :param [Vector3,Vector3,Vector3,Vector3] vertices: coordinates of vertices in the global coordinate system.
+ :param bool strictCheck: checks vertices order, raise RuntimeError for negative volume
+
+ See :yref:`yade.utils.sphere`'s documentation for meaning of other parameters."""
+ b=Body()
+ center = .25*sum(vertices,Vector3.Zero)
+ volume = TetrahedronSignedVolume(vertices)
+ if volume < 0:
+ if strictCheck:
+ raise RuntimeError, "tetra: wrong order of vertices"
+ temp = vertices[3]
+ vertices[3] = vertices[2]
+ vertices[2] = temp
+ volume = TetrahedronSignedVolume(vertices)
+ assert(volume>0)
+ b.shape = Tetra(v=vertices,color=color if color else randomColor(),wire=wire,highlight=highlight)
+ # modifies pos, ori and inertia
+ ori = TetrahedronWithLocalAxesPrincipal(b)
+ _commonBodySetup(b,volume,b.state.inertia,material,noBound=noBound,pos=center,fixed=fixed)
+ b.state.ori = b.state.refOri = ori
+ b.aspherical = True
+ b.mask = mask
+ b.chain = chain
+ return b
+
+
+
+
#def setNewVerticesOfFacet(b,vertices):
# center = inscribedCircleCenter(vertices[0],vertices[1],vertices[2])
# vertices = Vector3(vertices[0])-center,Vector3(vertices[1])-center,Vector3(vertices[2])-center
