Commit b92b8461 authored by Jan Stransky's avatar Jan Stransky

Tetra modification (Ig2_Tetra_Tetra_TTtraSimpleGeom,...

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
parent 3aeba293
################################################################################
#
# 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()
################################################################################
#
# 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()
This diff is collapsed.
......@@ -28,7 +28,6 @@ void Facet::postLoad(Facet&)
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();
......
......@@ -37,7 +37,7 @@ class Facet : public Shape {
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]"))
......
#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;
};
#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);
......@@ -206,7 +206,8 @@ Real Shop::kineticEnergy(Scene* _scene, Body::id_t* maxId){
// 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 @@ Real Shop::kineticEnergy(Scene* _scene, Body::id_t* maxId){
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;
......
......@@ -84,6 +84,10 @@ class Shop{
//! 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);
......
This diff is collapsed.
// © 2007 Václav Šmilauer <[email protected]>
// © 2013 Jan Stránský <[email protected]>
#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 @@ class Tetra: public Shape{
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 @@ REGISTER_SERIALIZABLE(Tetra);
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 @@ REGISTER_SERIALIZABLE(Bo1_Tetra_Aabb);
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 @@ class Ig2_Tetra_Tetra_TTetraGeom: public IGeomFunctor
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);
......
......@@ -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 @@ BOOST_PYTHON_MODULE(_utils){
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");
}
......@@ -348,6 +348,37 @@ def facet(vertices,dynamic=None,fixed=True,wire=True,color=None,highlight=False,
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
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment