------------------------------------------------------------ revno: 3968 committer: Jan Stránský <jan.stran...@fsv.cvut.cz> 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 #####
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
