yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3968
committer: Jan Stránský <jan.stransky@xxxxxxxxxxx>
timestamp: Thu 2016-11-10 15:38:54 +0100
message:
  added pack.inHalfSpace and pack.inConvexPolygon + example script
added:
  examples/test/pack-inConvexPolyhedron.py
modified:
  py/pack/pack.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

=== added file 'examples/test/pack-inConvexPolyhedron.py'
--- examples/test/pack-inConvexPolyhedron.py	1970-01-01 00:00:00 +0000
+++ examples/test/pack-inConvexPolyhedron.py	2016-11-10 14:38:54 +0000
@@ -0,0 +1,62 @@
+from yade import pack
+import random
+random.seed(1)
+
+# a cube
+p1 = Vector3(-9,-2,-2)
+p2 = p1 + (1,1,1)
+pred1 = pack.inConvexPolyhedron((
+	(p1, (+1, 0, 0)),
+	(p1, ( 0,+1, 0)),
+	(p1, ( 0, 0,+1)),
+	(p2, (-1, 0, 0)),
+	(p2, ( 0,-1, 0)),
+	(p2, ( 0, 0,-1)),
+))
+
+# a tetrahedron
+p1 = Vector3(-7,-2,-2)
+p2 = p1 + (1,1,1)
+pred2 = pack.inConvexPolyhedron((
+	(p1, (+1, 0, 0)),
+	(p1, ( 0,+1, 0)),
+	(p1, ( 0, 0,+1)),
+	(p2, (-1,-1,-1)),
+))
+
+# random polyhedron
+center = Vector3(-2,-2,-2)
+dMin = 1
+dMax = 1.2
+n = 20
+planes = []
+for i in xrange(n):
+	d = random.random()*(dMax-dMin) + dMin # distance of plane from center, random number in range (dMin,dMax)
+	# http://mathworld.wolfram.com/HyperspherePointPicking.html
+	normal = Vector3([random.gauss(0,10) for _ in xrange(3)]) # random normal vector
+	normal.normalize()
+	planes.append((center-d*normal,normal))
+pred3 = pack.inConvexPolyhedron(planes)
+
+try: # should be ValueError, since the 3 planes does not form closed polyhedron and finding its bounds should fail
+	pred4 = pack.inConvexPolyhedron((
+		((0,0,0), (+1, 0, 0)),
+		((0,0,0), ( 0,+1, 0)),
+		((0,0,0), ( 0, 0,+1)),
+	))
+except ValueError:
+	print 'ValueError successfully detected'
+else:
+	raise RuntimeError, "ValueError should have been detected..."
+
+r = .05
+for p in (pred1,pred2,pred3):
+	O.bodies.append(pack.regularHexa(p,r,0,color=randomColor()))
+
+
+try:
+	from yade import qt
+	v = qt.View()
+	v.axes = True
+except:
+	pass

=== modified file 'py/pack/pack.py'
--- py/pack/pack.py	2015-12-09 00:13:58 +0000
+++ py/pack/pack.py	2016-11-10 14:38:54 +0000
@@ -174,6 +174,78 @@
       inf=float('inf'); return Vector3(inf,inf,inf)
     def __call__(self,pt,pad): return True
 
+  class inHalfSpace(Predicate):
+    """Predicate returning True  any points, with infinite bounding box."""
+    def __init__(self, _center=Vector3().Zero, _dir=Vector3(1,0,0)):
+      self._center = Vector3(_center)
+      self._dir = Vector3(_dir)
+      assert self._dir.norm() > 0., "Direction has to be nonzero vector"
+      self._dir.normalize()
+      self._d = -self._dir.dot(self._center)
+    def aabb(self):
+      d,c = self._dir, self._center
+      inf=float('inf')
+      min = Vector3(-inf,-inf,-inf)
+      max = Vector3(+inf,+inf,+inf)
+      for i in xrange(3):
+        j = (i+1)%3
+        k = (i+2)%3
+        if d[i]==0 and d[j]==0:
+          if d[k] > 0: min[k] = c[k]
+          else: max[k] = c[k]
+      return min,max
+    def center(self):
+      return self._center
+    def dim(self):
+      inf=float('inf')
+      return Vector3(inf,inf,inf)
+    def __call__(self,pt,pad):
+      v = self._dir.dot(pt) + self._d
+      return v > pad
+
+  class inConvexPolyhedron(Predicate):
+    def __init__(self,planes):
+      self._inHalfSpaces = [inHalfSpace(c,d) for c,d in planes]
+      self._min, self._max = self._computeAabb()
+    def _computeAabb(self):
+      try:
+        import scipy.optimize
+      except ImportError:
+        raise ImportError, "scipy (package python-scipy) needed for pack.inConvexPolyhedron"
+      min,max = Vector3.Zero, Vector3.Zero
+      A,b = [],[]
+      for h in self._inHalfSpaces:
+        A.append(tuple(-h._dir))
+        b.append(h._d)
+      inf = float('inf')
+      for i in xrange(3):
+        c = Vector3.Zero
+        #
+        c[i] = 1
+        opt = scipy.optimize.linprog(c,A_ub=A,b_ub=b,bounds=(-inf,inf))
+        errmsg = "Something wrong in pack.inConvexPolyhedron defining planes.\nThe scipy.optimize.linprog output:\n{}\n"
+        if not opt.success:
+          raise ValueError, errmsg.format(opt)
+        min[i] = opt.x[i]
+        #
+        c[i] = -1
+        opt = scipy.optimize.linprog(c,A_ub=A,b_ub=b,bounds=(-inf,inf))
+        if not opt.success:
+          raise ValueError, errmsg.format(opt)
+        max[i] = opt.x[i]
+      return min,max
+    def aabb(self):
+      return Vector3(self._min), Vector3(self._max)
+    def center(self):
+      return .5*(self._min + self._max)
+    def dim(self):
+      return self._max - self._min
+    def __call__(self,pt,pad):
+      for p in self._inHalfSpaces:
+        if not p(pt,pad):
+          return False
+      return True
+
 #####
 ## surface construction and manipulation
 #####