yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2604
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Fri 2010-12-10 11:13:47 +0100
message:
  1. Add support for ScGeom6D to LawTester (not tested...)
modified:
  pkg/dem/DomainLimiter.cpp
  pkg/dem/DomainLimiter.hpp


--
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/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp	2010-12-09 12:38:42 +0000
+++ pkg/dem/DomainLimiter.cpp	2010-12-10 10:13:47 +0000
@@ -77,6 +77,8 @@
 	Dem3DofGeom* d3dGeom=dynamic_cast<Dem3DofGeom*>(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 && !d3dGeom && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom (or L6Geom).");
@@ -96,45 +98,41 @@
 		return;
 	}
 	/* initialize or update local axes and trsf */
-	//DEL rotGeom=Vector3r(NaN,NaN,NaN); // this is meaningful only with l6geom
 	uGeom.end<3>()=Vector3r(NaN,NaN,NaN);
-	if(!l3Geom){
+	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
+		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(axX[1]>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.end<3>()=Vector3r::Zero();
 		} else { // udpate of an existing interaction
 			if(scGeom){
 				scGeom->rotate(axY); scGeom->rotate(axZ);
 				scGeom->rotate(shearTot);
 				shearTot+=scGeom->shearIncrement();
-				//DEL ptGeom=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ)); 
 				uGeom.start<3>()=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
+				if(scGeom6d) uGeom.end<3>()=Vector3r(scGeom6d->getTwist(),scGeom6d->getBending().dot(axY),scGeom6d->getBending().dot(axZ));
 			}
 			else{ // d3dGeom
 				throw runtime_error("LawTester: Dem3DofGeom not yet supported.");
 				// essentially copies code from ScGeom, which is not very nice indeed; oh well…
-				Vector3r vRel=(state2->vel+state2->angVel.cross(-gsc->refR2*gsc->normal))-(state1->vel+state1->angVel.cross(gsc->refR1*gsc->normal));
+				// Vector3r vRel=(state2->vel+state2->angVel.cross(-gsc->refR2*gsc->normal))-(state1->vel+state1->angVel.cross(gsc->refR1*gsc->normal));
 			}
 		}
 		// update the transformation
-		// the matrix is orthonormal, since axX, axY are normalized and and axZ is their corss-product
+		// 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=l3Geom->trsf;
 		axX=trsf.row(0); axY=trsf.row(1); axZ=trsf.row(2);
-		//DEL ptGeom=l3Geom->u;
 		uGeom.start<3>()=l3Geom->u;
-		if(l6Geom){
-			//rotGeom=l6Geom->phi;
-			uGeom.end<3>()=l6Geom->phi;
-			// perform allshearing by translation, as it does not induce bending
-			if(rotWeight!=0){ LOG_INFO("LawTester.rotWeight set to 0 (was"<<rotWeight<<"), since rotational DoFs are in use."); rotWeight=0; }
-		}
+		if(l6Geom) uGeom.end<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;
@@ -150,7 +148,6 @@
 		dU*=refLength;
 	}
 	LOG_DEBUG("Absolute diff is: displacement "<<dU<<", rotation "<<dPhi);
-	//DEL ptOurs+=dU; rotOurs+=dPhi;
 	uTest=uTestNext; // the value that was next in the previous step is the current one now
 	uTestNext.start<3>()+=dU; uTestNext.end<3>()+=dPhi;
 

=== modified file 'pkg/dem/DomainLimiter.hpp'
--- pkg/dem/DomainLimiter.hpp	2010-12-09 12:38:42 +0000
+++ pkg/dem/DomainLimiter.hpp	2010-12-10 10:13:47 +0000
@@ -33,7 +33,7 @@
 		((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."))
+		//((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."))
@@ -46,7 +46,7 @@
 		((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,0,,"Step number in which this engine is active; determines position in path, using pathSteps."))
+		((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."))