yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #02353
[Branch ~yade-dev/yade/trunk] Rev 1817: merge with clump branch:
Merge authors:
Sergei D. <sega@think>
------------------------------------------------------------
revno: 1817 [merge]
committer: Sergei D. <sega@think>
branch nick: trunk
timestamp: Fri 2009-11-27 12:31:55 +0300
message:
merge with clump branch:
1. handle DOFs for new rotation alg.;
2. add new clump-hopper-viscoelastic script;
3. small optimization of acc.rot.alg; fix velocity and add angVel in VTKRecorder.
added:
scripts/test/clump-hopper-viscoelastic.py
modified:
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp
pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp
pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription.
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-11-26 12:57:26 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.cpp 2009-11-27 09:27:08 +0000
@@ -20,6 +20,20 @@
acceleration [i]*= 1 - damping*Mathr::Sign ( f[i]*(velocity [i] + (Real) 0.5 *dt*acceleration [i]) );
}
}
+void NewtonsDampedLaw::blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v) {
+ if(blockedDOFs==State::DOF_NONE) return;
+ if(blockedDOFs==State::DOF_ALL) v = Vector3r::ZERO;
+ if((blockedDOFs & State::DOF_X)!=0) v[0]=0;
+ if((blockedDOFs & State::DOF_Y)!=0) v[1]=0;
+ if((blockedDOFs & State::DOF_Z)!=0) v[2]=0;
+}
+void NewtonsDampedLaw::blockRotateDOFs(unsigned blockedDOFs, Vector3r& v) {
+ if(blockedDOFs==State::DOF_NONE) return;
+ if(blockedDOFs==State::DOF_ALL) v = Vector3r::ZERO;
+ if((blockedDOFs & State::DOF_RX)!=0) v[0]=0;
+ if((blockedDOFs & State::DOF_RY)!=0) v[1]=0;
+ if((blockedDOFs & State::DOF_RZ)!=0) v[2]=0;
+}
void NewtonsDampedLaw::handleClumpMember(MetaBody* ncb, const body_id_t memberId, State* clumpState){
const shared_ptr<Body>& b=Body::byId(memberId,ncb);
assert(b->isClumpMember());
@@ -98,20 +112,26 @@
}
// blocking DOFs
- if(state->blockedDOFs==State::DOF_NONE){
- state->angVel+=dt*state->angAccel;
- state->vel+=+dt*state->accel;
- } else if(state->blockedDOFs==State::DOF_ALL){
- /* do nothing */
- } else {
- // handle more complicated cases here
- if((state->blockedDOFs & State::DOF_X)==0) state->vel[0]+=dt*state->accel[0];
- if((state->blockedDOFs & State::DOF_Y)==0) state->vel[1]+=dt*state->accel[1];
- if((state->blockedDOFs & State::DOF_Z)==0) state->vel[2]+=dt*state->accel[2];
- if((state->blockedDOFs & State::DOF_RX)==0) state->angVel[0]+=dt*state->angAccel[0];
- if((state->blockedDOFs & State::DOF_RY)==0) state->angVel[1]+=dt*state->angAccel[1];
- if((state->blockedDOFs & State::DOF_RZ)==0) state->angVel[2]+=dt*state->angAccel[2];
- }
+ blockTranslateDOFs( state->blockedDOFs, state->accel );
+ state->vel+=+dt*state->accel;
+
+ blockRotateDOFs( state->blockedDOFs, state->angAccel );
+ state->angVel+=dt*state->angAccel;
+
+ //if(state->blockedDOFs==State::DOF_NONE){
+ //state->angVel+=dt*state->angAccel;
+ //state->vel+=+dt*state->accel;
+ //} else if(state->blockedDOFs==State::DOF_ALL){
+ //[> do nothing <]
+ //} else {
+ //// handle more complicated cases here
+ //if((state->blockedDOFs & State::DOF_X)==0) state->vel[0]+=dt*state->accel[0];
+ //if((state->blockedDOFs & State::DOF_Y)==0) state->vel[1]+=dt*state->accel[1];
+ //if((state->blockedDOFs & State::DOF_Z)==0) state->vel[2]+=dt*state->accel[2];
+ //if((state->blockedDOFs & State::DOF_RX)==0) state->angVel[0]+=dt*state->angAccel[0];
+ //if((state->blockedDOFs & State::DOF_RY)==0) state->angVel[1]+=dt*state->angAccel[1];
+ //if((state->blockedDOFs & State::DOF_RZ)==0) state->angVel[2]+=dt*state->angAccel[2];
+ //}
// velocities are ready now, save maxima
if(haveBins) {velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(state));}
@@ -141,12 +161,10 @@
}
void NewtonsDampedLaw::accurateRigidBodyRotationIntegrator(MetaBody* ncb, const shared_ptr<Body>& rb){
- Real dt=Omega::instance().getTimeStep();
+ const Real dt=Omega::instance().getTimeStep();
State* state=rb->state.get();
const body_id_t id=rb->getId();
-
state->accel=state->angAccel=Vector3r::ZERO; // to make sure; should be reset in Clump::moveMembers
-
// sum of forces and torques
const Vector3r& f=ncb->bex.getForce(id); const Vector3r& m=ncb->bex.getTorque(id);
Vector3r F(f), M(m);
@@ -159,21 +177,23 @@
F+=f; M+=(bs->pos-state->pos).Cross(f)+m;
}
// translation equation
- state->accel=F/state->mass;
- state->vel+=dt*state->accel;
- state->pos+=state->vel*dt;//+ncb->bex.getMove(id);
+ state->accel=F/state->mass; blockTranslateDOFs( state->blockedDOFs, state->accel );
+ state->vel+=dt*state->accel; state->pos+=state->vel*dt+ncb->bex.getMove(id);
// rotation equation
- Vector3r l_n = state->prevAngMom + dt/2 * M; // global angular momentum at time n
- Vector3r l_b_n = state->ori.Conjugate().Rotate(l_n); // local angular momentum at time n
- Vector3r angVel_b_n = diagDiv(l_b_n,state->inertia); // local angular velocity at time n
- Quaternionr dotQ_n=DotQ(angVel_b_n,state->ori); // dQ/dt at time n
- Quaternionr Q_half = state->ori + dt/2 * dotQ_n; // Q at time n+1/2
+ Matrix3r A; state->ori.Conjugate().ToRotationMatrix(A); // rotation matrix from global to local r.f.
+ const Vector3r l_n = state->prevAngMom + dt/2 * M; // global angular momentum at time n
+ //Vector3r l_b_n = state->ori.Conjugate().Rotate(l_n); // local angular momentum at time n
+ const Vector3r l_b_n = A*l_n; // local angular momentum at time n
+ const Vector3r angVel_b_n = diagDiv(l_b_n,state->inertia); // local angular velocity at time n
+ const Quaternionr dotQ_n=DotQ(angVel_b_n,state->ori); // dQ/dt at time n
+ const Quaternionr Q_half = state->ori + dt/2 * dotQ_n; // Q at time n+1/2
state->prevAngMom+=dt*M; // global angular momentum at time n+1/2
- Vector3r l_b_half = state->ori.Conjugate().Rotate(state->prevAngMom); // local angular momentum at time n+1/2
+ //Vector3r l_b_half = state->ori.Conjugate().Rotate(state->prevAngMom); // local angular momentum at time n+1/2
+ const Vector3r l_b_half = A*state->prevAngMom; // local angular momentum at time n+1/2
Vector3r angVel_b_half = diagDiv(l_b_half,state->inertia); // local angular velocity at time n+1/2
- Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
- state->ori+=dt*dotQ_half; // Q at time n+1
- state->ori.Normalize();
+ blockRotateDOFs( state->blockedDOFs, angVel_b_half );
+ const Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
+ state->ori+=dt*dotQ_half; state->ori.Normalize(); // Q at time n+1
state->angVel=state->ori.Rotate(angVel_b_half); // global angular velocity at time n+1/2
//if(rb->isClump()) static_cast<Clump*>(rb.get())->moveMembers();
static_cast<Clump*>(rb.get())->moveMembers();
=== modified file 'pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp'
--- pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-11-26 12:57:26 +0000
+++ pkg/dem/Engine/StandAloneEngine/NewtonsDampedLaw.hpp 2009-11-27 09:27:08 +0000
@@ -43,6 +43,9 @@
bool haveBins;
void accurateRigidBodyRotationIntegrator(MetaBody* ncb, const shared_ptr<Body>& rb);
Quaternionr DotQ(const Vector3r& angVel, const Quaternionr& Q);
+ inline void blockTranslateDOFs(unsigned blockedDOFs, Vector3r& v);
+ inline void blockRotateDOFs(unsigned blockedDOFs, Vector3r& v);
+
public:
///damping coefficient for Cundall's non viscous damping
Real damping;
=== modified file 'pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp'
--- pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-11-26 12:57:26 +0000
+++ pkg/dem/Engine/StandAloneEngine/VTKRecorder.cpp 2009-11-27 07:01:33 +0000
@@ -60,13 +60,16 @@
vtkSmartPointer<vtkCellArray> spheresCells = vtkSmartPointer<vtkCellArray>::New();
vtkSmartPointer<vtkFloatArray> radii = vtkSmartPointer<vtkFloatArray>::New();
radii->SetNumberOfComponents(1);
- radii->SetName("Radii");
+ radii->SetName("radii");
vtkSmartPointer<vtkFloatArray> spheresColors = vtkSmartPointer<vtkFloatArray>::New();
spheresColors->SetNumberOfComponents(3);
- spheresColors->SetName("Colors");
+ spheresColors->SetName("colors");
vtkSmartPointer<vtkFloatArray> spheresVelocity = vtkSmartPointer<vtkFloatArray>::New();
spheresVelocity->SetNumberOfComponents(3);
spheresVelocity->SetName("velocity");
+ vtkSmartPointer<vtkFloatArray> sphAngVel = vtkSmartPointer<vtkFloatArray>::New();
+ sphAngVel->SetNumberOfComponents(3);
+ sphAngVel->SetName("angVel");
// facets
vtkSmartPointer<vtkPoints> facetsPos = vtkSmartPointer<vtkPoints>::New();
@@ -147,7 +150,16 @@
spheresColors->InsertNextTupleValue(c);
}
if(recActive[REC_VELOCITY]){
- spheresVelocity->InsertNextTupleValue((float*)(&b->state->vel));
+ const Vector3r& vel = b->state->vel;
+ float v[3] = { vel[0],vel[1],vel[2] };
+ spheresVelocity->InsertNextTupleValue(v);
+
+ const Vector3r& angVel = b->state->angVel;
+ float av[3] = { angVel[0],angVel[1],angVel[2] };
+ sphAngVel->InsertNextTupleValue(av);
+
+ //spheresVelocity->InsertNextTupleValue((float*)(&b->state->vel));
+ //sphAngVel->InsertNextTupleValue((float*)(&b->state->angVel));
}
if (recActive[REC_CPM]) {
cpmDamage->InsertNextValue(YADE_PTR_CAST<CpmState>(b->state)->normDmg);
@@ -199,7 +211,10 @@
spheresUg->SetCells(VTK_VERTEX, spheresCells);
spheresUg->GetPointData()->AddArray(radii);
if (recActive[REC_COLORS]) spheresUg->GetPointData()->AddArray(spheresColors);
- if (recActive[REC_VELOCITY]) spheresUg->GetPointData()->AddArray(spheresVelocity);
+ if (recActive[REC_VELOCITY]) {
+ spheresUg->GetPointData()->AddArray(spheresVelocity);
+ spheresUg->GetPointData()->AddArray(sphAngVel);
+ }
if (recActive[REC_CPM]) {
spheresUg->GetPointData()->AddArray(cpmDamage);
spheresUg->GetPointData()->AddArray(cpmSigma);
=== added file 'scripts/test/clump-hopper-viscoelastic.py'
--- scripts/test/clump-hopper-viscoelastic.py 1970-01-01 00:00:00 +0000
+++ scripts/test/clump-hopper-viscoelastic.py 2009-11-27 09:27:08 +0000
@@ -0,0 +1,87 @@
+# -*- coding: utf-8
+
+from yade import utils,pack,export,qt
+import gts,os,random,itertools
+from numpy import *
+import yade.log
+
+#yade.log.setLevel('NewtonsDampedLaw',yade.log.TRACE)
+
+# Parameters
+tc=0.001# collision time
+en=.3 # normal restitution coefficient
+es=.3 # tangential restitution coefficient
+frictionAngle=radians(35)#
+density=2700
+# facets material
+params=utils.getViscoelasticFromSpheresInteraction(10e3,tc,en,es)
+facetMat=O.materials.append(SimpleViscoelasticMat(frictionAngle=frictionAngle,**params)) # **params sets kn, cn, ks, cs
+# default spheres material
+dfltSpheresMat=O.materials.append(SimpleViscoelasticMat(density=density,frictionAngle=frictionAngle))
+
+O.dt=.1*tc # time step
+
+Rs=0.05 # particle radius
+
+# Create geometry
+
+x0=0.; y0=0.; z0=0.; ab=.7; at=2.; h=1.; hl=h; al=at*3
+
+zb=z0; x0b=x0-ab/2.; y0b=y0-ab/2.; x1b=x0+ab/2.; y1b=y0+ab/2.
+zt=z0+h; x0t=x0-at/2.; y0t=y0-at/2.; x1t=x0+at/2.; y1t=y0+at/2.
+zl=z0-hl;x0l=x0-al/2.; y0l=y0-al/2.; x1l=x0+al/2.; y1l=y0+al/2.
+
+left = pack.sweptPolylines2gtsSurface([[Vector3(x0b,y0b,zb),Vector3(x0t,y0t,zt),Vector3(x0t,y1t,zt),Vector3(x0b,y1b,zb)]],capStart=True,capEnd=True)
+lftIds=O.bodies.append(pack.gtsSurface2Facets(left.faces(),material=facetMat,color=(0,1,0)))
+
+right = pack.sweptPolylines2gtsSurface([[Vector3(x1b,y0b,zb),Vector3(x1t,y0t,zt),Vector3(x1t,y1t,zt),Vector3(x1b,y1b,zb)]],capStart=True,capEnd=True)
+rgtIds=O.bodies.append(pack.gtsSurface2Facets(right.faces(),material=facetMat,color=(0,1,0)))
+
+near = pack.sweptPolylines2gtsSurface([[Vector3(x0b,y0b,zb),Vector3(x0t,y0t,zt),Vector3(x1t,y0t,zt),Vector3(x1b,y0b,zb)]],capStart=True,capEnd=True)
+nearIds=O.bodies.append(pack.gtsSurface2Facets(near.faces(),material=facetMat,color=(0,1,0)))
+
+far = pack.sweptPolylines2gtsSurface([[Vector3(x0b,y1b,zb),Vector3(x0t,y1t,zt),Vector3(x1t,y1t,zt),Vector3(x1b,y1b,zb)]],capStart=True,capEnd=True)
+farIds=O.bodies.append(pack.gtsSurface2Facets(far.faces(),material=facetMat,color=(0,1,0)))
+
+table = pack.sweptPolylines2gtsSurface([[Vector3(x0l,y0l,zl),Vector3(x0l,y1l,zl),Vector3(x1l,y1l,zl),Vector3(x1l,y0l,zl)]],capStart=True,capEnd=True)
+tblIds=O.bodies.append(pack.gtsSurface2Facets(table.faces(),material=facetMat,color=(0,1,0)))
+
+# Create clumps...
+for k,l in itertools.product(arange(0,10),arange(0,10)):
+ clpId,sphId=O.bodies.appendClumped([utils.sphere(Vector3(x0t+Rs*(k*4+2),y0t+Rs*(l*4+2),i*Rs*2+zt),Rs,color=(0,0,1),material=dfltSpheresMat) for i in xrange(4)])
+ for id in sphId:
+ s=O.bodies[id]
+ p=utils.getViscoelasticFromSpheresInteraction(s.state['mass'],tc,en,es)
+ s.mat['kn'],s.mat['cn'],s.mat['ks'],s.mat['cs']=p['kn'],p['cn'],p['ks'],p['cs']
+ #O.bodies[clpId].state.blockedDOFs=['rx','ry','rz']
+ #O.bodies[clpId].state.blockedDOFs=['x','y']
+
+# ... and spheres
+for k,l in itertools.product(arange(0,9),arange(0,9)):
+ sphAloneId=O.bodies.append( [utils.sphere( Vector3(x0t+Rs*(k*4+4),y0t+Rs*(l*4+4),i*Rs*2.3+zt),Rs,color=(0,1,0),material=dfltSpheresMat) for i in xrange(4) ] )
+ for id in sphAloneId:
+ s=O.bodies[id]
+ p=utils.getViscoelasticFromSpheresInteraction(s.state['mass'],tc,en,es)
+ s.mat['kn'],s.mat['cn'],s.mat['ks'],s.mat['cs']=p['kn'],p['cn'],p['ks'],p['cs']
+ #s.state.blockedDOFs=['rx','ry','rz']
+ #s.state.blockedDOFs=['x','y']
+
+# Create engines
+O.engines=[
+ BexResetter(),
+ BoundingVolumeMetaEngine([InteractingSphere2AABB(),InteractingFacet2AABB()]),
+ InsertionSortCollider(),
+ InteractionDispatchers(
+ [InteractingSphere2InteractingSphere4SpheresContactGeometry(), InteractingFacet2InteractingSphere4SpheresContactGeometry()],
+ [Ip2_SimleViscoelasticMat_SimpleViscoelasticMat_SimpleViscoelasticPhys()],
+ [Law2_Spheres_Viscoelastic_SimpleViscoelastic()],
+ ),
+ GravityEngine(gravity=[0,0,-9.81]),
+ NewtonsDampedLaw(damping=0,accRigidBodyRot=True),
+ #VTKRecorder(virtPeriod=0.01,fileName='/tmp/',recorders=['spheres','velocity','facets'])
+]
+
+renderer = qt.Renderer()
+qt.View()
+O.saveTmp()
+
