yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3695
committer: Jan Stransky <jan.stransky@xxxxxxxxxxx>
timestamp: Wed 2013-09-11 23:43:09 +0200
message:
  added utils.UnstructuredGrid class for manipulating (external forces evaluation and positin changes) of FEM-like triangular and tetrahedral meshes
  hopefully fixed CGAL compilation problem in pkg/dem/Tetra.*pp reported by Klaus
modified:
  examples/polyhedra/oneTetra.py
  pkg/dem/Tetra.cpp
  pkg/dem/Tetra.hpp
  py/utils.py


--
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

=== modified file 'examples/polyhedra/oneTetra.py'
--- examples/polyhedra/oneTetra.py	2013-09-08 11:12:42 +0000
+++ examples/polyhedra/oneTetra.py	2013-09-11 21:43:09 +0000
@@ -1,6 +1,6 @@
 ################################################################################
 # 
-# Script to test tetra gl functions and prescribe dmotion
+# Script to test tetra gl functions and prescribed motion
 #
 ################################################################################
 v1 = (0,0,0)
@@ -29,6 +29,7 @@
 	NewtonIntegrator(),
 	PyRunner(iterPeriod=1,command="changeVelocity()"),
 ]
+O.step()
 
 
 try:

=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp	2013-09-08 11:12:42 +0000
+++ pkg/dem/Tetra.cpp	2013-09-11 21:43:09 +0000
@@ -3,33 +3,34 @@
 
 #include"Tetra.hpp"
 
+#include<yade/core/Interaction.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/core/State.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
+
+#include<yade/pkg/common/Aabb.hpp>
+
+#ifdef YADE_CGAL
+	#include <CGAL/intersections.h>
+#endif
+
 YADE_PLUGIN(/* self-contained in hpp: */ (Tetra) (TTetraGeom) (TTetraSimpleGeom) (Bo1_Tetra_Aabb) 
 	/* some code in cpp (this file): */ (TetraVolumetricLaw) 
 	(Ig2_Tetra_Tetra_TTetraGeom)
 	#ifdef YADE_CGAL
 		(Ig2_Tetra_Tetra_TTetraSimpleGeom)
+		(Law2_TTetraSimpleGeom_NormPhys_Simple)
 	#endif
 	#ifdef YADE_OPENGL
 		(Gl1_Tetra)
 	#endif	
-	(Law2_TTetraSimpleGeom_NormPhys_Simple)
 	);
 
 Tetra::~Tetra(){}
 TTetraGeom::~TTetraGeom(){}
 TTetraSimpleGeom::~TTetraSimpleGeom(){}
 
-#include<yade/core/Interaction.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-#include<yade/core/State.hpp>
-#include<yade/pkg/common/ElastMat.hpp>
-
-#include<yade/pkg/common/Aabb.hpp>
-
-#ifdef YADE_CGAL
-	#include <CGAL/intersections.h>
-#endif
 
 
 void Bo1_Tetra_Aabb::go(const shared_ptr<Shape>& ig, shared_ptr<Bound>& bv, const Se3r& se3, const Body*){
@@ -1186,6 +1187,7 @@
 Real TetrahedronVolume(const Vector3r v[4]) { return fabs(TetrahedronSignedVolume(v)); }
 Real TetrahedronSignedVolume(const vector<Vector3r>& v) { return Vector3r(v[1]-v[0]).dot(Vector3r(v[2]-v[0]).cross(v[3]-v[0]))/6.; }
 Real TetrahedronVolume(const vector<Vector3r>& v) { return fabs(TetrahedronSignedVolume(v)); }
+#ifdef YADE_CGAL
 Real TetrahedronVolume(const CGAL::Point_3<CGAL::Cartesian<Real> >* v[4]) {
 	Vector3r vv[4];
 	for (int i=0; i<4; i++) {
@@ -1204,9 +1206,11 @@
 	}
 	return TetrahedronVolume(vv);
 }
-
-
-
+#endif
+
+
+
+#ifdef YADE_CGAL
 void Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
 	int id1 = contact->getId1(), id2 = contact->getId2();
 	TTetraSimpleGeom* geom= static_cast<TTetraSimpleGeom*>(ig.get());
@@ -1223,4 +1227,5 @@
 	applyForceAtContactPoint(-phys->normalForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
 	// TODO periodic
 }
+#endif
 

=== modified file 'pkg/dem/Tetra.hpp'
--- pkg/dem/Tetra.hpp	2013-09-08 11:12:42 +0000
+++ pkg/dem/Tetra.hpp	2013-09-11 21:43:09 +0000
@@ -165,7 +165,6 @@
 	DECLARE_LOGGER;
 };
 REGISTER_SERIALIZABLE(Ig2_Tetra_Tetra_TTetraSimpleGeom);
-#endif
 
 
 
@@ -178,6 +177,7 @@
 	DECLARE_LOGGER;
 };
 REGISTER_SERIALIZABLE(Law2_TTetraSimpleGeom_NormPhys_Simple);
+#endif
 
 
 

=== modified file 'py/utils.py'
--- py/utils.py	2013-09-08 11:12:42 +0000
+++ py/utils.py	2013-09-11 21:43:09 +0000
@@ -954,3 +954,76 @@
 		self.numCM = len(relRadii)
 		self.relRadii = relRadii
 		self.relPositions = relPositions
+
+class UnstructuredGrid:
+	"""EXPERIMENTAL. TODO docs
+	"""
+	def __init__(self,vertices={},connectivityTable={},**kw):
+		self.vertices = vertices
+		self.connectivityTable = connectivityTable
+		self.forces = dict( (i,Vector3(0,0,0)) for i in vertices)
+		self.build(**kw)
+	def setup(self,vertices,connectivityTable,toSimulation=False,**kw):
+		self.vertices = dict( (i,v) for i,v in vertices.iteritems())
+		self.connectivityTable = dict( (i,e) for i,e in connectivityTable.iteritems())
+		self.build(**kw)
+		if toSimulation:
+			self.toSimulation()
+	def build(self,**kw):
+		self.elements = {}
+		for i,c in self.connectivityTable.iteritems():
+			if len(c) == 3:
+				b = facet([self.vertices[j] for j in c],**kw)
+			elif len(c) == 4:
+				b = tetra([self.vertices[j] for j in c],**kw)
+			else:
+				raise RuntimeError, "TODO"
+			self.elements[i] = b
+	def resetForces(self):
+		for i in self.vertices:
+			self.forces[i] = Vector3(0,0,0)
+	def getForcesOfNodes(self):
+		self.resetForces()
+		for i,e in self.elements.iteritems():
+			ie = self.connectivityTable[i]
+			for i in e.intrs():
+				fn = i.phys.normalForce
+				fs = i.phys.shearForce if hasattr(i.phys,"shearForce") else Vector3.Zero
+				f = (fn+fs) * (-1 if e.id == i.id1 else +1 if e.id == i.id2 else 'ERROR')
+				cp = i.geom.contactPoint
+				if isinstance(e.shape,Facet):
+					v0,v1,v2 = [Vector3(self.vertices[j]) for j in ie]
+					w0 = ((cp-v1).cross(v2-v1)).norm()
+					w1 = ((cp-v2).cross(v0-v2)).norm()
+					w2 = ((cp-v0).cross(v1-v0)).norm()
+					ww = w0+w1+w2
+					self.forces[ie[0]] += f*w0/ww
+					self.forces[ie[1]] += f*w1/ww
+					self.forces[ie[2]] += f*w2/ww
+				elif isinstance(e.shape,Tetra):
+					v0,v1,v2,v3 = [Vector3(self.vertices[j]) for j in ie]
+					w0 = TetrahedronVolume((cp,v1,v2,v3))
+					w1 = TetrahedronVolume((cp,v2,v3,v0))
+					w2 = TetrahedronVolume((cp,v3,v0,v1))
+					w3 = TetrahedronVolume((cp,v0,v1,v2))
+					ww = w0+w1+w2+w3
+					self.forces[ie[0]] += f*w0/ww
+					self.forces[ie[1]] += f*w1/ww
+					self.forces[ie[2]] += f*w2/ww
+					self.forces[ie[3]] += f*w3/ww
+		return self.forces
+	def setPositionsOfNodes(self,newPoss):
+		for i in self.vertices:
+			self.vertices[i] = newPoss[i]
+		self.updateElements()
+	def updateElements(self):
+		for i,c in self.connectivityTable.iteritems():
+			e = self.elements[i]
+			e.state.pos = Vector3(0,0,0)
+			e.state.ori = Quaternion((1,0,0),0)
+			if isinstance(e.shape,Facet):
+				e.shape.vertices = [self.vertices[j] for j in c]
+			elif isinstance(e.shape,Tetra):
+				e.shape.v = [self.vertices[j] for j in c]
+	def toSimulation(self):
+		O.bodies.append(self.elements.values())