← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2625: 1. Fix NewtonIntegrator for clumps, blockedDOFs &c; msg will be sent to yade-dev

 

------------------------------------------------------------
revno: 2625
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Tue 2010-12-21 23:50:34 +0100
message:
  1. Fix NewtonIntegrator for clumps, blockedDOFs &c; msg will be sent to yade-dev
  2. Make static attributes properly hyperlinked and documented in the ui
  3. Fix yref in L3Geom
modified:
  gui/qt4/Inspector.py
  gui/qt4/SerializableEditor.py
  pkg/dem/L3Geom.hpp
  pkg/dem/NewtonIntegrator.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 'gui/qt4/Inspector.py'
--- gui/qt4/Inspector.py	2010-12-21 15:31:40 +0000
+++ gui/qt4/Inspector.py	2010-12-21 22:50:34 +0000
@@ -67,15 +67,17 @@
 		return -1
 
 class BodyInspector(QWidget):
-	def __init__(self,bodyId=0,parent=None,bodyLinkCallback=None,intrLinkCallback=None):
+	def __init__(self,bodyId=-1,parent=None,bodyLinkCallback=None,intrLinkCallback=None):
 		QWidget.__init__(self,parent)
+		v=yade.qt.views()
 		self.bodyId=bodyId
-		self.idGlSync=-1
+		if bodyId<0 and len(v)>0 and v[0].selection>0: self.bodyId=v[0].selection
+		self.idGlSync=self.bodyId
 		self.bodyLinkCallback,self.intrLinkCallback=bodyLinkCallback,intrLinkCallback
 		self.bodyIdBox=QSpinBox(self)
 		self.bodyIdBox.setMinimum(0)
 		self.bodyIdBox.setMaximum(100000000)
-		self.bodyIdBox.setValue(0)
+		self.bodyIdBox.setValue(self.bodyId)
 		self.intrWithCombo=QComboBox(self);
 		self.gotoBodyButton=QPushButton(u'→ #',self)
 		self.gotoIntrButton=QPushButton(u'→ #+#',self)
@@ -166,6 +168,7 @@
 		# no body shown yet, try to get the first one...
 		if self.bodyId<0 and len(O.bodies)>0:
 			try:
+				print 'SET ZERO'
 				b=O.bodies[0]; self.bodyIdBox.setValue(0)
 			except IndexError: pass
 		v=yade.qt.views()

=== modified file 'gui/qt4/SerializableEditor.py'
--- gui/qt4/SerializableEditor.py	2010-12-21 15:31:40 +0000
+++ gui/qt4/SerializableEditor.py	2010-12-21 22:50:34 +0000
@@ -15,14 +15,15 @@
 seqSerializableShowType=True # show type headings in serializable sequences (takes vertical space, but makes the type hyperlinked)
 
 
-def makeWrapperHref(text,className,attr=None):
+def makeWrapperHref(text,className,attr=None,static=False):
 	"""Create clickable HTML hyperlink to a Yade class or its attribute.
 	
 	:param className: name of the class to link to.
 	:param attr: attribute to link to. If given, must exist directly in given *className*; if not given or empty, link to the class itself is created and *attr* is ignored.
 	:return: HTML with the hyperref.
 	"""
-	return '<a href="%s#yade.wrapper.%s%s">%s</a>'%(yade.qt.sphinxDocWrapperPage,className,(('.'+attr) if attr else ''),text)
+	if not static: return '<a href="%s#yade.wrapper.%s%s">%s</a>'%(yade.qt.sphinxDocWrapperPage,className,(('.'+attr) if attr else ''),text)
+	else:          return '<a href="%s#ystaticattr-%s.%s">%s</a>'%(yade.qt.sphinxDocWrapperPage,className,attr,text)
 
 def serializableHref(ser,attr=None,text=None):
 	"""Return HTML href to a *ser* optionally to the attribute *attr*.
@@ -48,7 +49,7 @@
 	else:
 		klass=ser.__class__
 		if not text: text=klass.__name__
-	return makeWrapperHref(text,klass.__name__,attr)
+	return makeWrapperHref(text,klass.__name__,attr,static=(attr and getattr(klass,attr)==getattr(ser,attr)))
 
 class AttrEditor():
 	"""Abstract base class handing some aspects common to all attribute editors.
@@ -346,8 +347,7 @@
 			t=None
 			doc=getattr(self.ser.__class__,attr).__doc__;
 			if '|yhidden|' in doc: continue
-			if attr in self.ignoredAttrs or attr=='blockedDOFs': # HACK here
-				continue
+			if attr in self.ignoredAttrs: continue
 			if isinstance(val,list):
 				t=self.getListTypeFromDocstring(attr)
 				if not t and len(val)==0: t=(val[0].__class__,) # 1-tuple is list of the contained type
@@ -364,10 +364,28 @@
 		doc=(getattr(self.ser.__class__,attr).__doc__ if attr else self.ser.__class__.__doc__)
 		if not doc: return ''
 		doc=re.sub(':y(attrtype|default|attrflags):`[^`]*`','',doc)
+		statAttr=re.compile('^.. ystaticattr::.*$',re.MULTILINE|re.DOTALL)
+		doc=re.sub(statAttr,'',doc) # static classes have their proper docs at the beginning, discard static memeber docs
+		# static: attribute of the type is the same object as attribute of the instance
+		# in that case, get docstring from the class documentation by parsing it
+		if attr and getattr(self.ser.__class__,attr)==getattr(self.ser,attr): doc=self.getStaticAttrDocstring(attr)
 		doc=re.sub(':yref:`([^`]*)`','\\1',doc)
 		import textwrap
 		wrapper=textwrap.TextWrapper(replace_whitespace=False)
 		return wrapper.fill(textwrap.dedent(doc))
+	def getStaticAttrDocstring(self,attr):
+		ret=''; c=self.ser.__class__
+		while hasattr(c,attr) and hasattr(c.__base__,attr): c=c.__base__
+		start='.. ystaticattr:: %s.%s('%(c.__name__,attr)
+		if start in c.__doc__:
+			ll=c.__doc__.split('\n')
+			for i in range(len(ll)):
+				if ll[i].startswith(start): break
+			for i in range(i+1,len(ll)):
+				if len(ll[i])>0 and ll[i][0] not in ' \t': break
+				ret+=ll[i]
+			return ret
+		else: return '[no documentation found]'
 	def mkWidget(self,entry):
 		if not entry.T: return None
 		# single fundamental object
@@ -407,8 +425,6 @@
 		grid.setVerticalSpacing(0)
 		grid.setLabelAlignment(Qt.AlignRight)
 		if self.showType:
-			#lab=QLabel(makeSerializableLabel(self.ser,addr=True,href=True))
-			#lab.setToolTip(('<b>'+self.path+'</b><br>' if self.path else '')+self.getDocstring())
 			lab=SerQLabel(self,makeSerializableLabel(self.ser,addr=True,href=True),tooltip=self.getDocstring(),path=self.path)
 			lab.setFrameShape(QFrame.Box); lab.setFrameShadow(QFrame.Sunken); lab.setLineWidth(2); lab.setAlignment(Qt.AlignHCenter); lab.linkActivated.connect(yade.qt.openUrl)
 			grid.setWidget(0,QFormLayout.SpanningRole,lab)
@@ -416,8 +432,6 @@
 			entry.widget=self.mkWidget(entry)
 			objPath=(self.path+'.'+entry.name) if self.path else None
 			label=SerQLabel(self,serializableHref(self.ser,entry.name),tooltip=self.getDocstring(entry.name),path=objPath)
-			#label=SerAttrQLabel(self,serializableHref(self.ser,entry.name),tooltip=('<b>'+self.path+'.'+entry.name+'</b><br>' if self.path else '')+self.getDocstring(entry.name),
-			#label=QLabel(self); label.setText(); label.setToolTip(; label.linkActivated.connect(yade.qt.openUrl)
 			grid.addRow(label,entry.widget if entry.widget else QLabel('<i>unhandled type</i>'))
 		self.setLayout(grid)
 		self.refreshEvent()

=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp	2010-12-13 12:11:43 +0000
+++ pkg/dem/L3Geom.hpp	2010-12-21 22:50:34 +0000
@@ -74,7 +74,7 @@
 	FUNCTOR1D(L6Geom);
 	void go(const shared_ptr<IGeom>&, const shared_ptr<Interaction>&, const shared_ptr<Body>&, const shared_ptr<Body>&, bool);
 	YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_L6Geom,Gl1_L3Geom,"Render :yref:`L6Geom` geometry.",
-		((Real,phiScale,1.,,"Scale local rotations (:yref:`phi<L3Geom.phi>` - :yref:`phi0<L3Geom.phi0>`). The default scale is to draw $\\pi$ rotation with length equal to minimum radius."))
+		((Real,phiScale,1.,,"Scale local rotations (:yref:`phi<L6Geom.phi>` - :yref:`phi0<L6Geom.phi0>`). The default scale is to draw $\\pi$ rotation with length equal to minimum radius."))
 	);
 };
 REGISTER_SERIALIZABLE(Gl1_L6Geom);

=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp	2010-12-21 18:31:15 +0000
+++ pkg/dem/NewtonIntegrator.cpp	2010-12-21 22:50:34 +0000
@@ -98,11 +98,8 @@
 			if(unlikely(!b)) continue;
 			State* state=b->state.get();
 			const Body::id_t& id=b->getId();
-			// clump members are non-dynamic; we only get their velocities here
-			if(unlikely(b->isClumpMember())){
-				saveMaximaVelocity(scene,id,state);
-				continue;
-			}
+			// clump members are handled inside clumps
+			if(unlikely(b->isClumpMember())) continue;
 			const Vector3r& f=scene->forces.getForce(id);
 			const Vector3r& m=scene->forces.getTorque(id);
 			#ifdef YADE_DEBUG
@@ -137,30 +134,27 @@
 
 			if (likely(b->isStandalone())){
 				// translate equation
-				if(b->isDynamic()){
-					state->accel=computeAccel(f,state->mass,state->blockedDOFs);
-					cundallDamp(dt,f,fluctVel,state->accel);}
+				state->accel=computeAccel(f,state->mass,state->blockedDOFs);
+				cundallDamp(dt,f,fluctVel,state->accel);
 				leapfrogTranslate(scene,state,id,dt);
 				// rotate equation
 				// exactAsphericalRot is disabled or the body is spherical
-				if (likely(!exactAsphericalRot || !b->isAspherical() || !b->isDynamic())){
-					if(b->isDynamic()) {
-						state->angAccel=computeAngAccel(m,state->inertia,state->blockedDOFs);
-						cundallDamp(dt,m,state->angVel,state->angAccel);}
+				if (likely(!exactAsphericalRot || !b->isAspherical() || b->state->inertia.minCoeff()<=0.)){
+					state->angAccel=computeAngAccel(m,state->inertia,state->blockedDOFs);
+					cundallDamp(dt,m,state->angVel,state->angAccel);
 					leapfrogSphericalRotate(scene,state,id,dt);
 				} else { // exactAsphericalRot enabled & aspherical body
 					// no damping in this case
 					leapfrogAsphericalRotate(scene,state,id,dt,m);
 				}
 			} else if (b->isClump()){
-				if(b->isDynamic()){
-					// reset acceleration of the clump itself; computed from accels on constituents
-					state->accel=state->angAccel=Vector3r::Zero();
-					// clump mass forces
-					Vector3r dLinAccel=computeAccel(f,state->mass,state->blockedDOFs);
-					cundallDamp(dt,f,fluctVel,dLinAccel);
-					state->accel+=dLinAccel;}
-				Vector3r M(m);
+				// reset acceleration of the clump itself; computed from accels on constituents
+				state->accel=state->angAccel=Vector3r::Zero();
+				// clump mass forces
+				Vector3r dLinAccel=computeAccel(f,state->mass,state->blockedDOFs);
+				cundallDamp(dt,f,fluctVel,dLinAccel);
+				state->accel+=dLinAccel;
+				Vector3r M(m); // torque on the clump itself, will accumulate from members
 				// sum force on clump memebrs
 				// exactAsphericalRot enabled and clump is aspherical
 				if (exactAsphericalRot && b->isAspherical() && b->state->inertia.maxCoeff()>0){
@@ -168,38 +162,35 @@
 						const Body::id_t& memberId=mm.first;
 						const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
 						State* memberState=member->state.get();
-						if(b->isDynamic()){
-							handleClumpMemberAccel(scene,memberId,memberState,state);
-							handleClumpMemberTorque(scene,memberId,memberState,state,M);}
-						//FIXME : we are saving max velocity of previous timestep here, and actually it's already saved above
-						saveMaximaVelocity(scene,memberId,memberState);
+						handleClumpMemberAccel(scene,memberId,memberState,state);
+						handleClumpMemberTorque(scene,memberId,memberState,state,M);
 					}
 					// motion
 					leapfrogTranslate(scene,state,id,dt);
 					leapfrogAsphericalRotate(scene,state,id,dt,M);
 				} else { // exactAsphericalRot disabled or clump is spherical
-					if(b->isDynamic()){
-						Vector3r dAngAccel=computeAngAccel(M,state->inertia,state->blockedDOFs);
-						cundallDamp(dt,M,state->angVel,dAngAccel);
-						state->angAccel+=dAngAccel;}
+					Vector3r dAngAccel=computeAngAccel(M,state->inertia,state->blockedDOFs);
+					cundallDamp(dt,M,state->angVel,dAngAccel);
+					state->angAccel+=dAngAccel;
 					FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
 						const Body::id_t& memberId=mm.first;
 						const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
 						State* memberState=member->state.get();
-						if(b->isDynamic()){
-							handleClumpMemberAccel(scene,memberId,memberState,state);
-							handleClumpMemberAngAccel(scene,memberId,memberState,state);}
-						//FIXME : we are saving max velocity of previous timestep here, and actually it's already saved above
-						saveMaximaVelocity(scene,memberId,memberState);
+						handleClumpMemberAccel(scene,memberId,memberState,state);
+						handleClumpMemberAngAccel(scene,memberId,memberState,state);
 					}
 					// motion
 					leapfrogTranslate(scene,state,id,dt);
 					leapfrogSphericalRotate(scene,state,id,dt);
 				}
 				Clump::moveMembers(b,scene);
+				// save max velocity for all members of the clump
+				FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
+					const shared_ptr<Body>& member=Body::byId(mm.first,scene); assert(member->isClumpMember());
+					saveMaximaVelocity(scene,mm.first,member->state.get());
+				}
 			}
 			saveMaximaVelocity(scene,id,state);
-
 			// process callbacks
 			for(size_t i=0; i<callbacksSize; i++){
 				cerr<<"<"<<b->id<<",cb="<<callbacks[i]<<",scene="<<callbacks[i]->scene<<">"; // <<",force="<<callbacks[i]->scene->forces.getForce(b->id)<<">";
@@ -228,13 +219,13 @@
 	} else if (homoDeform==Cell::HOMO_POS){
 		state->pos+=scene->cell->velGrad*state->pos*dt;
 	}
-	if(state->blockedDOFs!=State::DOF_ALL) state->vel+=dt*state->accel;
-	state->pos += state->vel*dt;
+	state->vel+=dt*state->accel;
+	state->pos+=state->vel*dt;
 }
 
 inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const Body::id_t& id, const Real& dt )
 {
-	if(state->blockedDOFs!=State::DOF_ALL) state->angVel+=dt*state->angAccel;
+	state->angVel+=dt*state->angAccel;
 	Vector3r axis = state->angVel;
 
 	if (axis!=Vector3r::Zero()) {							//If we have an angular velocity, we make a rotation
@@ -252,6 +243,14 @@
 }
 
 void NewtonIntegrator::leapfrogAsphericalRotate(Scene* scene, State* state, const Body::id_t& id, const Real& dt, const Vector3r& M){
+	// if all DoFs are blocked then return; warn if there is angular acceleration, since it will be ignored
+	if((state->blockedDOFs & State::DOF_RXRYRZ)==State::DOF_RXRYRZ){
+		/* necessarily angular acceleration is zero, we can call the spherical integrator safely */
+		if(state->angVel!=Vector3r::Zero()) leapfrogSphericalRotate(scene,state,id,dt); /* in most cases, angVel will be zero; no need to integrate anything, just return */
+		return;
+	}
+	// if only some rotations are allowed, warn, since it is not handled (yet?)
+	if((state->blockedDOFs & State::DOF_RXRYRZ)!=0) LOG_WARN("Aspherical body #"<<id<<" with some rotational DoFs blocked is not supported (blocked DoFs ignored).");
 	Matrix3r A=state->ori.conjugate().toRotationMatrix(); // rotation matrix from global to local r.f.
 	const Vector3r l_n = state->angMom + dt/2 * M; // global angular momentum at time n
 	const Vector3r l_b_n = A*l_n; // local angular momentum at time n
@@ -261,7 +260,6 @@
 	state->angMom+=dt*M; // global angular momentum at time n+1/2
 	const Vector3r l_b_half = A*state->angMom; // local angular momentum at time n+1/2
 	Vector3r angVel_b_half = l_b_half.cwise()/state->inertia; // local angular velocity at time n+1/2
-	//blockRotateDOFs( state->blockedDOFs, angVel_b_half ); // FIXME: angVel_b_half is the local velocity, but blockedDOFs need for the _global_ velocity.
 	const Quaternionr dotQ_half=DotQ(angVel_b_half,Q_half); // dQ/dt at time n+1/2
 	state->ori=state->ori+dt*dotQ_half; // Q at time n+1
 	state->angVel=state->ori*angVel_b_half; // global angular velocity at time n+1/2