yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #10134
[Branch ~yade-pkg/yade/git-trunk] Rev 3714: Recover DomainLimiter (LawTester should be moved somewhere).
------------------------------------------------------------
revno: 3714
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Tue 2013-10-15 07:01:48 +0200
message:
Recover DomainLimiter (LawTester should be moved somewhere).
added:
pkg/dem/DomainLimiter.cpp
pkg/dem/DomainLimiter.hpp
--
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 file 'pkg/dem/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/DomainLimiter.cpp 2013-10-15 05:01:48 +0000
@@ -0,0 +1,343 @@
+#include<yade/pkg/dem/DomainLimiter.hpp>
+#include<yade/pkg/dem/Shop.hpp>
+
+YADE_PLUGIN((DomainLimiter)(LawTester)
+ #ifdef YADE_OPENGL
+ (GlExtra_LawTester)(GlExtra_OctreeCubes)
+ #endif
+);
+
+void DomainLimiter::action(){
+ std::list<Body::id_t> out;
+ FOREACH(const shared_ptr<Body>& b, *scene->bodies){
+ if((!b) or ((mask>0) and ((b->groupMask & mask)==0))) continue;
+ const Sphere* sphere = dynamic_cast<Sphere*>(b->shape.get());
+ if (sphere){ //Delete only spheres
+ const Vector3r& p(b->state->pos);
+ if(p[0]<lo[0] || p[0]>hi[0] || p[1]<lo[1] || p[1]>hi[1] || p[2]<lo[2] || p[2]>hi[2]) {
+ out.push_back(b->id);
+ nDeleted++;
+ mDeleted+=b->state->mass;
+ Real r = sphere->radius; vDeleted+=(4/3.)*Mathr::PI*pow(r,3);
+ }
+ }
+ }
+ FOREACH(Body::id_t id, out){
+ scene->bodies->erase(id);
+ }
+}
+
+#include<yade/pkg/dem/DemXDofGeom.hpp>
+#include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/pkg/dem/L3Geom.hpp>
+#include<yade/pkg/common/NormShearPhys.hpp>
+#include<yade/lib/smoothing/LinearInterpolate.hpp>
+#include<yade/lib/pyutil/gil.hpp>
+
+CREATE_LOGGER(LawTester);
+
+void LawTester::postLoad(LawTester&){
+ if(ids.size()==0) return; // uninitialized object, don't do nothing at all
+ if(ids.size()!=2) throw std::invalid_argument("LawTester.ids: exactly two values must be given.");
+ if(disPath.empty() && rotPath.empty()) throw invalid_argument("LawTester.{disPath,rotPath}: at least one point must be given.");
+ if(pathSteps.empty()) throw invalid_argument("LawTester.pathSteps: at least one value must be given.");
+ size_t pathSize=max(disPath.size(),rotPath.size());
+ // update path points
+ _path.clear(); _path.push_back(Vector6r::Zero());
+ for(size_t i=0; i<pathSize; i++) {
+ Vector6r pt;
+ pt.head<3>()=Vector3r(i<disPath.size()?disPath[i]:(disPath.empty()?Vector3r::Zero():*(disPath.rbegin())));
+ pt.tail<3>()=Vector3r(i<rotPath.size()?rotPath[i]:(rotPath.empty()?Vector3r::Zero():*(rotPath.rbegin())));
+ _path.push_back(pt);
+
+ }
+ // update time points from distances, repeat last distance if shorter than path
+ _pathT.clear(); _pathT.push_back(0);
+ for(size_t i=0; i<pathSteps.size(); i++) _pathT.push_back(_pathT[i]+pathSteps[i]);
+ int lastDist=pathSteps[pathSteps.size()-1];
+ for(size_t i=pathSteps.size(); i<pathSize; i++) _pathT.push_back(*(_pathT.rbegin())+lastDist);
+}
+
+void LawTester::action(){
+ Vector2r foo; // avoid undefined ~Vector2r with clang?
+ if(ids.size()!=2) throw std::invalid_argument("LawTester.ids: exactly two values must be given.");
+ LOG_DEBUG("=================== LawTester step "<<step<<" ========================");
+ const shared_ptr<Interaction> Inew=scene->interactions->find(ids[0],ids[1]);
+ string strIds("##"+lexical_cast<string>(ids[0])+"+"+lexical_cast<string>(ids[1]));
+ // interaction not found at initialization
+ if(!I && (!Inew || !Inew->isReal())){
+ LOG_WARN("Interaction "<<strIds<<" does not exist (yet?), no-op."); return;
+ //throw std::runtime_error("LawTester: interaction "+strIds+" does not exist"+(Inew?" (to be honest, it does exist, but it is not real).":"."));
+ }
+ // interaction was deleted meanwhile
+ if(I && (!Inew || !Inew->isReal())) throw std::runtime_error("LawTester: interaction "+strIds+" was deleted"+(Inew?" (is not real anymore).":"."));
+ // different interaction object
+ if(I && Inew && I!=Inew) throw std::logic_error("LawTester: interacion "+strIds+" is a different object now?!");
+ assert(Inew);
+ bool doInit=(!I);
+ if(doInit) I=Inew;
+
+ id1=I->getId1(); id2=I->getId2();
+ // test object types
+ GenericSpheresContact* gsc=dynamic_cast<GenericSpheresContact*>(I->geom.get());
+ ScGeom* scGeom=dynamic_cast<ScGeom*>(I->geom.get());
+ L3Geom* l3Geom=dynamic_cast<L3Geom*>(I->geom.get());
+ L6Geom* l6Geom=dynamic_cast<L6Geom*>(I->geom.get());
+ ScGeom6D* scGeom6d=dynamic_cast<ScGeom6D*>(I->geom.get());
+ bool hasRot=(l6Geom || scGeom6d);
+ //NormShearPhys* phys=dynamic_cast<NormShearPhys*>(I->phys.get()); //Disabled because of warning
+ if(!gsc) throw std::invalid_argument("LawTester: IGeom of "+strIds+" not a GenericSpheresContact.");
+ if(!scGeom && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom (or L6Geom).");
+ assert(!((bool)scGeom && (bool)l3Geom)); // nonsense
+ // get body objects
+ State *state1=Body::byId(id1,scene)->state.get(), *state2=Body::byId(id2,scene)->state.get();
+ scene->forces.sync();
+ if(state1->blockedDOFs!=State::DOF_ALL) { LOG_INFO("Blocking all DOFs for #"<<id1); state1->blockedDOFs=State::DOF_ALL;}
+ if(state2->blockedDOFs!=State::DOF_ALL) { LOG_INFO("Blocking all DOFs for #"<<id2); state2->blockedDOFs=State::DOF_ALL;}
+
+
+ if(step-1>*(_pathT.rbegin())){
+ LOG_INFO("Last step done, setting zero velocities on #"<<id1<<", #"<<id2<<".");
+ state1->vel=state1->angVel=state2->vel=state2->angVel=Vector3r::Zero();
+ uTest=uTestNext;
+ if(doneHook.empty()){ LOG_INFO("No doneHook set, dying."); dead=true; }
+ else{ LOG_INFO("Running doneHook: "<<doneHook); pyRunString(doneHook);}
+ return;
+ }
+ /* initialize or update local axes and trsf */
+ uGeom.tail<3>()=Vector3r(NaN,NaN,NaN);
+ if(!l3Geom){ // IGeom's that don't have local axes
+ axX=gsc->normal; /* just in case */ axX.normalize();
+ if(doInit){ // initialization of the new interaction -- define local axes
+ // take vector in the y or z direction, depending on its length; arbitrary, but one of them is sure to be non-zero
+ axY=axX.cross(abs(axX[1])<abs(axX[2])?Vector3r::UnitY():Vector3r::UnitZ());
+ axY.normalize();
+ axZ=axX.cross(axY);
+ LOG_DEBUG("Initial axes x="<<axX<<", y="<<axY<<", z="<<axZ);
+ if(scGeom6d) uGeom.tail<3>()=Vector3r::Zero();
+ } else { // udpate of an existing interaction
+ if(scGeom){
+ scGeom->rotate(axY); scGeom->rotate(axZ);
+ scGeom->rotate(shearTot);
+ shearTot+=scGeom->shearIncrement();
+ uGeom.head<3>()=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
+ if(scGeom6d) uGeom.tail<3>()=-1.*Vector3r(scGeom6d->getTwist(),scGeom6d->getBending().dot(axY),scGeom6d->getBending().dot(axZ));
+ }
+ else{ // d3dGeom
+ throw runtime_error("Geom type not yet supported.");
+ }
+ }
+ // update the transformation
+ // the matrix is orthonormal, since axX, axY are normalized and and axZ is their cross-product
+ trsf.row(0)=axX; trsf.row(1)=axY; trsf.row(2)=axZ;
+ } else {
+ trsf=Matrix3r(l3Geom->trsf);
+ axX=trsf.row(0); axY=trsf.row(1); axZ=trsf.row(2);
+ uGeom.head<3>()=l3Geom->u;
+ if(l6Geom) uGeom.tail<3>()=l6Geom->phi;
+ }
+ // perform all shearing by translation, as it does not induce bending
+ if(hasRot && rotWeight!=0){ LOG_INFO("LawTester.rotWeight set to 0 (was "<<rotWeight<<"), since rotational DoFs are in use."); rotWeight=0; }
+ contPt=gsc->contactPoint;
+ refLength=gsc->refR1+gsc->refR2;
+ renderLength=.5*refLength;
+
+ // here we go ahead, finally
+ Vector6r uu=linearInterpolate<Vector6r,int>(step,_pathT,_path,_interpPos);
+ Vector6r dUU=uu-uuPrev; uuPrev=uu;
+ Vector3r dU(dUU.head<3>()), dPhi(dUU.tail<3>());
+ //Vector3r dU=u-uPrev.head<3>(); uPrev.head<3>()=u;
+ //Vector3r dPhi=phi-uPrev.tail<3>(); uPrev.tail<3>()=phi;
+ if(displIsRel){
+ LOG_DEBUG("Relative displacement diff is "<<dU<<" (will be normalized by "<<gsc->refR1+gsc->refR2<<")");
+ dU*=refLength;
+ }
+ LOG_DEBUG("Absolute diff is: displacement "<<dU<<", rotation "<<dPhi);
+ uTest=uTestNext; // the value that was next in the previous step is the current one now
+ uTestNext.head<3>()+=dU; uTestNext.tail<3>()+=dPhi;
+
+ // reset velocities where displacement is controlled
+ //for(int i=0; i<3; i++){ if(forceControl[i]==0){ state1.vel[i]=0; state2.vel[i]=0; }
+
+ // shear is applied as rotation of id2: dε=r₁dθ → dθ=dε/r₁;
+ Vector3r vel[2],angVel[2];
+ //State* states[]={state1,state2};
+ for(int i=0; i<2; i++){
+ int sign=(i==0?-1:1);
+ Real weight=(i==0?1-idWeight:idWeight);
+ // FIXME: this should not use refR1, but real CP-particle distance perhaps?
+ Real radius=(i==0?gsc->refR1:gsc->refR2);
+ Real relRad=radius/refLength;
+ // signed and weighted displacement/rotation to be applied on this sphere (reversed for #0)
+ // some rotations must cancel the sign, by multiplying by sign again
+ Vector3r ddU=sign*dU*weight;
+
+ // twist can be still distributed with idWeight (!)
+ Vector3r ddPhi=sign*dPhi*(1-relRad); /* shear angles must distribute to both, otherwise it would induce shear */
+ ddPhi[0]=sign*dPhi[0]*weight; // twist can be still distributed with idWeight
+ vel[i]=angVel[i]=Vector3r::Zero();
+
+ // normal displacement
+
+ vel[i]+=axX*ddU[0]/scene->dt;
+
+ // shear rotation
+
+ // multiplication by sign cancels sign in ddU, since rotation is non-symmetric (to increase shear, both spheres have the same rotation)
+ // (unlike shear displacement, which is symmetric)
+ // rotation around Z (which gives y-shear) must be inverted: +ry gives +εzm while -rz gives +εy
+ Real rotZ=-sign*rotWeight*ddU[1]/radius, rotY=sign*rotWeight*ddU[2]/radius;
+ angVel[i]+=(rotY*axY+rotZ*axZ)/scene->dt;
+
+ // shear displacement
+
+ // angle that is traversed by a sphere in order to give desired ddU when displaced on the branch of r1+r2
+ // FIXME: is the branch value correct here?!
+ Real arcAngleY=atan((1-rotWeight)*ddU[1]/radius), arcAngleZ=atan((1-rotWeight)*ddU[2]/radius);
+ // same, but without the atan, which can be disregarded for small increments:
+ // Real arcAngleY=(1-rotWeight)*ddU[1]/radius, arcAngleZ=(1-rotWeight)*ddU[2]/radius;
+ vel[i]+=axY*radius*sin(arcAngleY)/scene->dt; vel[i]+=axZ*radius*sin(arcAngleZ)/scene->dt;
+
+ // compensate distance increase caused by motion along the perpendicular axis
+ // cos(argAngle*) is always positive, regardless of the orientation
+ // and the compensation is always in the -εx sense (-sign → +1 for #0, -1 for #1)
+ vel[i]+=-sign*axX*radius*((1-cos(arcAngleY))+(1-cos(arcAngleZ)))/scene->dt;
+
+ // rotation, convert from local to global
+ angVel[i]+=trsf.transpose()*ddPhi;
+
+ LOG_DEBUG("vel="<<vel[i]<<", angVel="<<angVel[i]<<", rotY,rotZ="<<rotY<<","<<rotZ<<", arcAngle="<<arcAngleY<<","<<arcAngleZ<<", sign="<<sign<<", weight="<<weight);
+
+ }
+ state1->vel=vel[0]; state1->angVel=angVel[0];
+ state2->vel=vel[1]; state2->angVel=angVel[1];
+ LOG_DEBUG("Body #"<<id1<<", setting vel="<<vel[0]<<", angVel="<<angVel[0]);
+ LOG_DEBUG("Body #"<<id2<<", setting vel="<<vel[1]<<", angVel="<<angVel[1]);
+
+ /* find out where are we at in the path, run hooks if approriate */
+ // _pathT has the first (zero) value added by us, so we skip it
+ int nPathT=_pathT.size();
+ for(int i=1; i<nPathT; i++){
+ // i-th point on _pathT is (i-1)th on path; run corresponding hook, if it exists
+ if(step==_pathT[i] && ((int) hooks.size())>(i-1) && !hooks[i-1].empty()) pyRunString(hooks[i-1]);
+ }
+ step++;
+}
+
+#ifdef YADE_OPENGL
+#include<yade/lib/opengl/OpenGLWrapper.hpp>
+#include<yade/lib/opengl/GLUtils.hpp>
+#include<yade/pkg/common/GLDrawFunctors.hpp>
+#include<yade/pkg/common/OpenGLRenderer.hpp>
+#include<GL/glu.h>
+
+CREATE_LOGGER(GlExtra_LawTester);
+
+void GlExtra_LawTester::render(){
+ // scene object changed (after reload, for instance), for re-initialization
+ if(tester && tester->scene!=scene) tester=shared_ptr<LawTester>();
+
+ if(!tester){ FOREACH(shared_ptr<Engine> e, scene->engines){ tester=dynamic_pointer_cast<LawTester>(e); if(tester) break; } }
+ if(!tester){ LOG_ERROR("No LawTester in O.engines, killing myself."); dead=true; return; }
+
+ //if(tester->renderLength<=0) return;
+ glColor3v(Vector3r(1,0,1));
+
+ // switch to local coordinates
+ glTranslatev(tester->contPt);
+ //glMultMatrixd(Eigen::Affine3d(tester->trsf.transpose()).data());
+ glMultMatrix(Eigen::Transform<Real,3,Eigen::Affine>(tester->trsf.transpose()).data());
+
+
+ glDisable(GL_LIGHTING);
+ //glColor3v(Vector3r(1,0,1));
+ //glBegin(GL_LINES); glVertex3v(Vector3r::Zero()); glVertex3v(.1*Vector3r::Ones()); glEnd();
+ //GLUtils::GLDrawText(string("This is the contact point!"),Vector3r::Zero(),Vector3r(1,0,1));
+
+ // local axes
+ glLineWidth(2.);
+ for(int i=0; i<3; i++){
+ Vector3r pt=Vector3r::Zero(); pt[i]=.5*tester->renderLength; Vector3r color=.3*Vector3r::Ones(); color[i]=1;
+ GLUtils::GLDrawLine(Vector3r::Zero(),pt,color);
+ GLUtils::GLDrawText(string(i==0?"x":(i==1?"y":"z")),pt,color);
+ }
+
+ // put the origin to the initial (no-shear) point, so that the current point appears at the contact point
+ glTranslatev(Vector3r(0,tester->uTestNext[1],tester->uTestNext[2]));
+
+
+ const int t(tester->step); const vector<int>& TT(tester->_pathT); const vector<Vector6r>& VV(tester->_path);
+ size_t numSegments=TT.size();
+ const Vector3r colorBefore=Vector3r(.7,1,.7), colorAfter=Vector3r(1,.7,.7);
+
+ // scale displacement, if they have the strain meaning
+ Real scale=1;
+ if(tester->displIsRel) scale=tester->refLength;
+
+ // find maximum displacement, draw axes in the shear plane
+ Real displMax=0;
+ FOREACH(const Vector6r& v, VV) displMax=max(v.head<3>().squaredNorm(),displMax);
+ displMax=1.2*scale*sqrt(displMax);
+
+ glLineWidth(1.);
+ GLUtils::GLDrawLine(Vector3r(0,-displMax,0),Vector3r(0,displMax,0),Vector3r(.5,0,0));
+ GLUtils::GLDrawLine(Vector3r(0,0,-displMax),Vector3r(0,0,displMax),Vector3r(.5,0,0));
+
+ // draw displacement path
+ glLineWidth(4.);
+ for(size_t segment=0; segment<numSegments-1; segment++){
+ // different colors before and after the current point
+ Real t0=TT[segment],t1=TT[segment+1];
+ const Vector3r &from=-VV[segment].head<3>()*scale, &to=-VV[segment+1].head<3>()*scale;
+ // current segment
+ if(t>t0 && t<t1){
+ Real norm=(t-t0)/(t1-t0);
+ GLUtils::GLDrawLine(from,from+(to-from)*norm,colorBefore);
+ GLUtils::GLDrawLine(from+(to-from)*norm,to,colorAfter);
+ } else { // other segment
+ GLUtils::GLDrawLine(from,to,t<t0?colorAfter:colorBefore);
+ }
+ }
+
+ glLineWidth(1.);
+}
+
+
+void GlExtra_OctreeCubes::postLoad(GlExtra_OctreeCubes&){
+ if(boxesFile.empty()) return;
+ boxes.clear();
+ ifstream txt(boxesFile.c_str());
+ while(!txt.eof()){
+ Real data[8];
+ for(int i=0; i<8; i++){ if(i<7 && txt.eof()) goto done; txt>>data[i]; }
+ OctreeBox ob; Vector3r mn(data[0],data[1],data[2]), mx(data[3],data[4],data[5]);
+ ob.center=.5*(mn+mx); ob.extents=(.5*(mx-mn)); ob.level=(int)data[6]; ob.fill=(int)data[7];
+ // for(int i=0; i<=ob.level; i++) cerr<<"\t"; cerr<<ob.level<<": "<<mn<<"; "<<mx<<"; "<<ob.center<<"; "<<ob.extents<<"; "<<ob.fill<<endl;
+ boxes.push_back(ob);
+ }
+ done:
+ std::cerr<<"GlExtra_OctreeCubes::postLoad: loaded "<<boxes.size()<<" boxes."<<std::endl;
+}
+
+void GlExtra_OctreeCubes::render(){
+ FOREACH(const OctreeBox& ob, boxes){
+ if(ob.fill<fillRangeDraw[0] || ob.fill>fillRangeDraw[1]) continue;
+ if(ob.level<levelRangeDraw[0] || ob.level>levelRangeDraw[1]) continue;
+ bool doFill=(ob.fill>=fillRangeFill[0] && ob.fill<=fillRangeFill[1] && (ob.fill!=0 || !noFillZero));
+ // -2: empty
+ // -1: recursion limit, empty
+ // 0: subdivided
+ // 1: recursion limit, full
+ // 2: full
+ Vector3r color=(ob.fill==-2?Vector3r(1,0,0):(ob.fill==-1?Vector3r(1,1,0):(ob.fill==0?Vector3r(0,0,1):(ob.fill==1)?Vector3r(0,1,0):(ob.fill==2)?Vector3r(0,1,1):Vector3r(1,1,1))));
+ glColor3v(color);
+ glPushMatrix();
+ glTranslatev(ob.center);
+ glScalef(2*ob.extents[0],2*ob.extents[1],2*ob.extents[2]);
+ if (doFill) glutSolidCube(1);
+ else glutWireCube(1);
+ glPopMatrix();
+ }
+}
+
+#endif /* YADE_OPENGL */
=== added file 'pkg/dem/DomainLimiter.hpp'
--- pkg/dem/DomainLimiter.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/DomainLimiter.hpp 2013-10-15 05:01:48 +0000
@@ -0,0 +1,98 @@
+
+#include<yade/pkg/common/PeriodicEngines.hpp>
+#include<yade/core/PartialEngine.hpp>
+#include<yade/pkg/common/Sphere.hpp>
+
+class DomainLimiter: public PeriodicEngine{
+ public:
+ virtual void action();
+ YADE_CLASS_BASE_DOC_ATTRS(DomainLimiter,PeriodicEngine,"Delete particles that are out of axis-aligned box given by *lo* and *hi*.",
+ ((Vector3r,lo,Vector3r(0,0,0),,"Lower corner of the domain."))
+ ((Vector3r,hi,Vector3r(0,0,0),,"Upper corner of the domain."))
+ ((long,nDeleted,0,Attr::readonly,"Cummulative number of particles deleted."))
+ ((Real,mDeleted,0,,"Mass of deleted particles."))
+ ((Real,vDeleted,0,,"Volume of deleted particles."))
+ ((int,mask,-1,,"If mask is defined, only particles with corresponding groupMask will be deleted."))
+ );
+};
+REGISTER_SERIALIZABLE(DomainLimiter);
+
+class LawTester: public PartialEngine{
+ Body::id_t id1,id2; // shorthands for local use
+ public:
+ void init();
+ virtual void action();
+ void postLoad(LawTester&);
+ void warnDeprec(const string& s1, const string& s2){ if(!warnedDeprecPtRot){ warnedDeprecPtRot=true; LOG_WARN("LawTester."<<s1<<" is deprecated, use LawTester."<<s2<<" instead.");} }
+ Vector3r get_ptOurs(){ warnDeprec("ptOurs","uTest.head()"); return uTest.head<3>(); } Vector3r get_ptGeom(){ warnDeprec("ptGeom","uGeom.head()"); return uGeom.head<3>(); }
+ Vector3r get_rotOurs(){ warnDeprec("rotOurs","uTest.tail()"); return uTest.tail<3>(); } Vector3r get_rotGeom(){ warnDeprec("rotGeom","uGeom.tail()"); return uGeom.tail<3>(); }
+ DECLARE_LOGGER;
+ YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(LawTester,PartialEngine,"Prescribe and apply deformations of an interaction in terms of local normal and shear displacements and rotations (using either :yref:`disPpath<LawTester.disPath>` and :yref:`rotPath<LawTester.rotPath>` [or :yref:`path<LawTester.path>` in the future]). Supported :yref:`IGeom` types are :yref:`ScGeom`, :yref:`L3Geom` and :yref:`L6Geom`. \n\nSee :ysrc:`scripts/test/law-test.py` for an example.",
+ ((vector<Vector3r>,disPath,,Attr::triggerPostLoad,"Loading path, where each Vector3 contains desired normal displacement and two components of the shear displacement (in local coordinate system, which is being tracked automatically. If shorter than :yref:`rotPath<LawTester.rotPath>`, the last value is repeated."))
+ ((vector<Vector3r>,rotPath,,Attr::triggerPostLoad,"Rotational components of the loading path, where each item contains torsion and two bending rotations in local coordinates. If shorter than :yref:`path<LawTester.path>`, the last value is repeated."))
+ ((vector<string>,hooks,,,"Python commands to be run when the corresponding point in path is reached, before doing other things in that particular step. See also :yref:`doneHook<LawTester.doneHook>`. "))
+ ((Vector6r,uGeom,Vector6r::Zero(),,"Current generalized displacements (3 displacements, 3 rotations), as stored in the interation itself. They should corredpond to :yref:`uTest<LawTester.uTest>`, otherwise a bug is indicated."))
+ ((Vector6r,uTest,Vector6r::Zero(),,"Current generalized displacements (3 displacements, 3 rotations), as they should be according to this :yref:`LawTester`. Should correspond to :yref:`uGeom<LawTester.uGeom>`."))
+ ((Vector6r,uTestNext,Vector6r::Zero(),Attr::hidden,"The value of uTest in the next step; since uTest is computed before uGeom is updated (in the next time step), uTest and uGeom would be always shifted by one timestep."))
+ ((bool,warnedDeprecPtRot,false,Attr::hidden,"Flag to say that the user was already warned about using deprecated ptOurg/ptGeom/rotOurs/rotGeom."))
+ ((Vector3r,shearTot,Vector3r::Zero(),Attr::hidden,"Current displacement in global coordinates; only used internally with ScGeom."))
+ ((bool,displIsRel,true,,"Whether displacement values in *disPath* are normalized by reference contact length (r1+r2 for 2 spheres)."))
+ //((Vector3i,forceControl,Vector3i::Zero(),,"Select which components of path (non-zero value) have force (stress) rather than displacement (strain) meaning."))
+ ((vector<int>,pathSteps,((void)"(constant step)",vector<int>(1,1)),Attr::triggerPostLoad,"Step number for corresponding values in :yref:`path<LawTester.path>`; if shorter than path, distance between last 2 values is used for the rest."))
+ ((vector<int>,_pathT,,(Attr::readonly|Attr::noSave),"Time value corresponding to points on path, computed from *pathSteps*. Length is always the same as path."))
+ ((vector<Vector6r>,_path,,(Attr::readonly|Attr::noSave),"Generalized displacement path values, computed from *disPath* and *rotPath* by appending zero initial value, and possibly repeating the last value to make lengths of *disPath* and *rotPath* match."))
+ ((shared_ptr<Interaction>,I,,(Attr::hidden),"Interaction object being tracked."))
+
+ // axX, axY, axZ could be replaced by references to rows in trsf; perhaps do that in the future (in such a case, trsf would not be noSave)
+ ((Vector3r,axX,,Attr::hidden,"Local x-axis in global coordinates (normal of the contact) |yupdate|"))
+ ((Vector3r,axY,,Attr::hidden,"Local y-axis in global coordinates; perpendicular to axX; initialized arbitrarily, but tracked to be consistent. |yupdate|"))
+ ((Vector3r,axZ,,(Attr::hidden|Attr::noSave),"Local z-axis in global coordinates; computed from axX and axY. |yupdate|"))
+ ((Matrix3r,trsf,,(Attr::noSave|Attr::readonly),"Transformation matrix for the local coordinate system. |yupdate|"))
+ ((size_t,_interpPos,0,(Attr::readonly|Attr::hidden),"State cookie for the interpolation routine."))
+ ((Vector6r,uuPrev,Vector6r::Zero(),(Attr::readonly),"Generalized displacement values reached in the previous step, for knowing which increment to apply in the current step."))
+ ((int,step,1,,"Step number in which this engine is active; determines position in path, using pathSteps."))
+ ((string,doneHook,,,"Python command (as string) to run when end of the path is achieved. If empty, the engine will be set :yref:`dead<Engine.dead>`."))
+ ((Real,renderLength,0,,"Characteristic length for the purposes of rendering, set equal to the smaller radius."))
+ ((Real,refLength,0,(Attr::readonly),"Reference contact length, for rendering only."))
+ ((Vector3r,contPt,Vector3r::Zero(),Attr::hidden,"Contact point (for rendering only)"))
+ ((Real,idWeight,1,,"Float, usually ∈〈0,1〉, determining on how are displacements distributed between particles (0 for id1, 1 for id2); intermediate values will apply respective part to each of them. This parameter is ignored with 6-DoFs :yref:`IGeom`."))
+ ((Real,rotWeight,1,,"Float ∈〈0,1〉 determining whether shear displacement is applied as rotation or displacement on arc (0 is displacement-only, 1 is rotation-only). Not effective when mutual rotation is specified."))
+ // reset force components along individual axes, instead of blocking DOFs which have no specific direction (for the force control)
+ , /* deprec */ ((path,disPath,"LawTester.path will be used for generalized displacement (6-component) loading path in the future."))
+ , /* init */
+ , /* ctor */
+ , /* py */ .add_property("ptOurs",&LawTester::get_ptOurs,"first 3 components of uTest |ydeprecated|") .add_property("ptGeom",&LawTester::get_ptGeom,"first 3 components of uGeom |ydeprecated|") .add_property("rotOurs",&LawTester::get_rotOurs,"last 3 components of uTest |ydeprecated|") .add_property("rotGeom",&LawTester::get_rotGeom,"last 3 components of uGeom |ydeprecated|")
+ );
+};
+REGISTER_SERIALIZABLE(LawTester);
+
+#ifdef YADE_OPENGL
+#include<yade/pkg/common/OpenGLRenderer.hpp>
+
+class GlExtra_LawTester: public GlExtraDrawer{
+ public:
+ DECLARE_LOGGER;
+ virtual void render();
+ YADE_CLASS_BASE_DOC_ATTRS(GlExtra_LawTester,GlExtraDrawer,"Find an instance of :yref:`LawTester` and show visually its data.",
+ ((shared_ptr<LawTester>,tester,,,"Associated :yref:`LawTester` object."))
+ );
+};
+REGISTER_SERIALIZABLE(GlExtra_LawTester);
+
+class GlExtra_OctreeCubes: public GlExtraDrawer{
+ public:
+ struct OctreeBox{ Vector3r center, extents; int fill; int level; };
+ std::vector<OctreeBox> boxes;
+ void postLoad(GlExtra_OctreeCubes&);
+ virtual void render();
+ YADE_CLASS_BASE_DOC_ATTRS(GlExtra_OctreeCubes,GlExtraDrawer,"Render boxed read from file",
+ ((string,boxesFile,,Attr::triggerPostLoad,"File to read boxes from; ascii files with ``x0 y0 z0 x1 y1 z1 c`` records, where ``c`` is an integer specifying fill (0 for wire, 1 for filled)."))
+ ((Vector2i,fillRangeFill,Vector2i(2,2),,"Range of fill indices that will be filled."))
+ ((Vector2i,fillRangeDraw,Vector2i(-2,2),,"Range of fill indices that will be rendered."))
+ ((Vector2i,levelRangeDraw,Vector2i(-2,2),,"Range of levels that will be rendered."))
+ ((bool,noFillZero,true,,"Do not fill 0-fill boxed (those that are further subdivided)"))
+ );
+};
+REGISTER_SERIALIZABLE(GlExtra_OctreeCubes);
+#endif
+
