yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2055
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Mon 2010-03-01 13:58:22 +0100
message:
  1. Add eudoxos.InteractionLocator and eudoxos.IntrSmooth3 for 3d smoothing of arbitrary quantities defined on interactions.
  2. Adjust some docstring formatting
modified:
  py/_eudoxos.cpp
  py/_packObb.cpp
  py/_packPredicates.cpp
  py/_packSpheres.cpp
  py/eudoxos.py
  py/utils.py


--
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 'py/_eudoxos.cpp'
--- py/_eudoxos.cpp	2010-02-09 16:50:30 +0000
+++ py/_eudoxos.cpp	2010-03-01 12:58:22 +0000
@@ -109,14 +109,86 @@
 	ret["vel"]=vel;
 	return ret;
 }
-
+#ifdef YADE_VTK
+
+#include<vtkPointLocator.h>
+#include<vtkIdList.h>
+#include<vtkUnstructuredGrid.h>
+#include<vtkPoints.h>
+
+/* Fastly locate interactions withing given distance from a given point. See python docs for details. */
+class InteractionLocator{
+	// object doing the work, see http://www.vtk.org/doc/release/5.2/html/a01048.html
+	vtkPointLocator *locator;
+	// used by locator
+	vtkUnstructuredGrid* grid;
+	vtkPoints* points;
+	// maps vtk's id to Interaction objects
+	vector<shared_ptr<Interaction> > intrs;
+	// axis-aligned bounds of all interactions
+	Vector3r mn,mx;
+	// count of interactions we hold
+	int cnt;
+	public:
+	InteractionLocator(){
+		// traverse all real interactions in the current simulation
+		// add them to points, create grid with those points
+		// store shared_ptr's to interactions in intrs separately
+		Scene* scene=Omega::instance().getScene().get();
+		locator=vtkPointLocator::New();
+		points=vtkPoints::New();
+		FOREACH(const shared_ptr<Interaction>& i, *scene->interactions){
+			if(!i->isReal()) continue;
+			Dem3DofGeom* ge=dynamic_cast<Dem3DofGeom*>(i->interactionGeometry.get());
+			if(!ge) continue;
+			const Vector3r& cp(ge->contactPoint);
+			int id=points->InsertNextPoint(cp[0],cp[1],cp[2]);
+			if(intrs.size()<=(size_t)id){intrs.resize(id+1000);}
+			intrs[id]=i;
+		}
+		double bb[6];
+		points->ComputeBounds(); points->GetBounds(bb);
+		mn=Vector3r(bb[0],bb[2],bb[4]); mx=Vector3r(bb[1],bb[3],bb[5]);
+		cnt=points->GetNumberOfPoints();
+
+		grid=vtkUnstructuredGrid::New();
+		grid->SetPoints(points);
+		locator->SetDataSet(grid);
+		locator->BuildLocator();
+	}
+	// cleanup
+	~InteractionLocator(){ locator->Delete(); points->Delete(); grid->Delete(); }
+
+	python::list intrsWithinDistance(const Vector3r& pt, Real radius){
+		vtkIdList *ids=vtkIdList::New();
+		locator->FindPointsWithinRadius(radius,pt,ids);
+		int numIds=ids->GetNumberOfIds();
+		python::list ret;
+		for(int i=0; i<numIds; i++){
+			// LOG_TRACE("Add "<<i<<"th id "<<ids->GetId(i));
+			ret.append(intrs[ids->GetId(i)]);
+		}
+		return ret;
+	}
+	python::tuple getBounds(){ return python::make_tuple(mn,mx);}
+	int getCnt(){ return cnt; }
+};
+#endif
 
 BOOST_PYTHON_MODULE(_eudoxos){
 	import_array();
+	YADE_SET_DOCSTRING_OPTS;
 	def("velocityTowardsAxis",velocityTowardsAxis,velocityTowardsAxis_overloads(args("axisPoint","axisDirection","timeToAxis","subtractDist","perturbation")));
 	def("yieldSigmaTMagnitude",yieldSigmaTMagnitude);
 	// def("spiralSphereStresses2d",spiralSphereStresses2d,(python::arg("dH_dTheta"),python::arg("axis")=2));
 	def("particleConfinement",particleConfinement);
 	def("testNumpy",testNumpy);
+#ifdef YADE_VTK
+	class_<InteractionLocator>("InteractionLocator","Locate all (real) interactions in space by their :yref:`contact point<Dem3DofGeom::contactPoint>`. When constructed, all real interactions are spatially indexed (uses vtkPointLocator internally). Use intrsWithinDistance to use those data. \n\n.. note::\n\tData might become inconsistent with real simulation state if simulation is being run between creation of this object and spatial queries.")
+		.def("intrsWithinDistance",&InteractionLocator::intrsWithinDistance,((python::arg("point"),python::arg("maxDist"))),"Return list of real interactions that are not further than *maxDist* from *point*.")
+		.add_property("bounds",&InteractionLocator::getBounds,"Return coordinates of lower and uppoer corner of axis-aligned abounding box of all interactions")
+		.add_property("count",&InteractionLocator::getCnt,"Number of interactions held")
+	;
+#endif
 }
 

=== modified file 'py/_packObb.cpp'
--- py/_packObb.cpp	2010-02-09 20:22:04 +0000
+++ py/_packObb.cpp	2010-03-01 12:58:22 +0000
@@ -7,6 +7,8 @@
 #include<yade/extra/boost_python_len.hpp>
 #include<yade/lib-base/Logging.hpp>
 #include<yade/lib-base/Math.hpp>
+#include<yade/lib-pyutil/doc_opts.hpp>
+
 #include<vector>
 #include<stdexcept>
 using namespace boost;
@@ -70,6 +72,7 @@
 }
 
 BOOST_PYTHON_MODULE(_packObb){
+	YADE_SET_DOCSTRING_OPTS;
 	python::scope().attr("__doc__")="Computation of oriented bounding box for cloud of points.";
 	python::def("cloudBestFitOBB",bestFitOBB_py,"Return (Vector3 center, Vector3 halfSize, Quaternion orientation) of\nbest-fit oriented bounding-box for given tuple of points\n(uses brute-force velome minimization, do not use for very large clouds).");
 };

=== modified file 'py/_packPredicates.cpp'
--- py/_packPredicates.cpp	2010-02-09 20:22:04 +0000
+++ py/_packPredicates.cpp	2010-03-01 12:58:22 +0000
@@ -3,6 +3,8 @@
 #include<yade/extra/boost_python_len.hpp>
 #include<yade/lib-base/Logging.hpp>
 #include<yade/lib-base/Math.hpp>
+#include<yade/lib-pyutil/doc_opts.hpp>
+
 // #include<yade/gui-py/_utils.hpp> // will be: yade/lib-py/_utils.hpp> at some point
 
 using namespace boost;
@@ -321,6 +323,7 @@
 
 BOOST_PYTHON_MODULE(_packPredicates){
 	python::scope().attr("__doc__")="Spatial predicates for volumes (defined analytically or by triangulation).";
+	YADE_SET_DOCSTRING_OPTS;
 	// base predicate class
 	python::class_<PredicateWrap,/* necessary, as methods are pure virtual*/ boost::noncopyable>("Predicate")
 		.def("__call__",python::pure_virtual(&Predicate::operator()))

=== modified file 'py/_packSpheres.cpp'
--- py/_packSpheres.cpp	2010-01-09 12:09:13 +0000
+++ py/_packSpheres.cpp	2010-03-01 12:58:22 +0000
@@ -1,9 +1,11 @@
 // 2009 © Václav Šmilauer <eudoxos@xxxxxxxx>
 
 #include<yade/pkg-dem/SpherePack.hpp>
+#include<yade/lib-pyutil/doc_opts.hpp>
 
 BOOST_PYTHON_MODULE(_packSpheres){
 	python::scope().attr("__doc__")="Creation, manipulation, IO for generic sphere packings.";
+	YADE_SET_DOCSTRING_OPTS;
 	python::class_<SpherePack>("SpherePack","Set of spheres as centers and radii",python::init<python::optional<python::list> >(python::args("list"),"Empty constructor, optionally taking list [ ((cx,cy,cz),r), … ] for initial data." ))
 		.def("add",&SpherePack::add,"Add single sphere to packing, given center as 3-tuple and radius")
 		.def("toList",&SpherePack::toList,"Return packing data as python list.")

=== modified file 'py/eudoxos.py'
--- py/eudoxos.py	2009-12-18 09:48:16 +0000
+++ py/eudoxos.py	2010-03-01 12:58:22 +0000
@@ -13,6 +13,52 @@
 from math import *
 from yade._eudoxos import * ## c++ implementations
 
+
+class IntrSmooth3d():
+	r"""Return spatially weigted gaussian average of arbitrary quantity defined on interactions.
+
+	At construction time, all real interactions are put inside spatial grid, permitting fast search for
+	points in neighbourhood defined by distance.
+
+	Parameters for the distribution are standard deviation :math:`\sigma` and relative cutoff distance
+	*relThreshold* (3 by default) which will discard points farther than *relThreshold* :math:`\times \sigma`.
+
+	Given central point :math:`p_0`, points are weighted by gaussian function
+
+	.. math::
+
+		\rho(p_0,p)=\frac{1}{\sigma\sqrt{2\pi}}\exp\left(\frac{-||p_0-p||^2}{2\sigma^2}\right)
+
+	To get the averaged value, simply call the instance, passing central point and callable object which received interaction object and returns the desired quantity:
+
+		>>> is3d=IntrSmooth3d(0.003)
+		>>> is3d(Vector3r(0,0,0),lambda i: i.phys.omega)
+	
+	"""
+	def __init__(self,stDev):
+		self.locator=InteractionLocator()
+		self.stDev=stDev
+		self.relThreshold=3.
+		self.sqrt2pi=sqrt(2*pi)
+		import yade.config
+		if not 'vtk' in yade.config.features: raise RuntimeError("IntrSmooth3d is only function with VTK-enabled builds.")
+	def _ptpt2weight(self,pt0,pt1):
+		distSq=(pt0-pt1).SquaredLength()
+		return (1./(self.stDev*self.sqrt2pi))*exp(-distSq/(2*self.stDev*self.stDev))
+	def bounds(self): return self.locator.bounds()
+	def count(self): return self.locator.count()
+	def __call__(self,pt,extr):
+		intrs=self.locator.intrsWithinDistance(pt,self.stDev*self.relThreshold)
+		if len(intrs)==0: return float('nan')
+		weights,val=0.,0.
+		for i in intrs:
+			weight=self._ptpt2weight(pt,i.geom.contactPoint)
+			val+=weight*extr(i)
+			weights+=weight
+			#print i,weight,extr(i)
+		return val/weights
+		
+
 def estimateStress(strain,cutoff=0.):
 	"""Use summed stored energy in contacts to compute macroscopic stress over the same volume, provided known strain."""
 	# E=(1/2)σεAl # global stored energy

=== modified file 'py/utils.py'
--- py/utils.py	2010-02-25 10:20:07 +0000
+++ py/utils.py	2010-03-01 12:58:22 +0000
@@ -262,7 +262,7 @@
 		wallMask=16
 	if (((extents[0]==0) and (extents[1]==0)) or ((extents[0]==0) and (extents[2]==0)) or ((extents[1]==0) and (extents[2]==0))):
 		raise RuntimeError("Please, specify at least 2 none-zero dimensions in extents!");
-	"""___________________________"""
+	# ___________________________
 	
 	mn,mx=[-extents[i] for i in 0,1,2],[extents[i] for i in 0,1,2]
 	def doWall(a,b,c,d):
@@ -291,13 +291,20 @@
 	Return List of facets forming the cylinder;
 	
 	:Parameters:
-			`center`: center of the created cylinder; (X,Y,Z) coordinates;
-			`radius`: cylinder radius;
-			`height`: cylinder height;
-			`orientation`: orientation of the cylinder in quaternion format;
-			`segmentsNumber`: the number of edges on the cylinder surface, the minimum is 5;
-			`wallMask`: bitmask; determines which walls will be created, in the order up (1), down (2), side (4). The numbers are ANDed; the default 7 means to create all walls;
-			`**kw`: passed to utils.facet;
+			`center`: Vector3
+				center of the created cylinder
+			`radius`: float
+				cylinder radius
+			`height`: float
+				cylinder height
+			`orientation`: Quaternion
+				orientation of the cylinder
+			`segmentsNumber`: int
+				number of edges on the cylinder surface (>=5)
+			`wallMask`: bitmask
+				determines which walls will be created, in the order up (1), down (2), side (4). The numbers are ANDed; the default 7 means to create all walls;
+			`**kw`: (unused keyword arguments)
+				passed to utils.facet;
 	"""
 	
 	#Defense from zero dimensions
@@ -310,7 +317,7 @@
 	if (wallMask>7):
 		print "wallMask must be 7 or less"
 		wallMask=7
-	"""___________________________"""
+	# ___________________________
 	import numpy
 	anglesInRad = numpy.linspace(0, 2.0*math.pi, segmentsNumber+1, endpoint=True)
 	P1=[]; P2=[]