#5511: arbitrary mesh functions
-------------------------+--------------------------------------------------
Reporter: jason | Owner: was
Type: enhancement | Status: new
Priority: major | Milestone: sage-3.4.1
Component: graphics | Keywords:
-------------------------+--------------------------------------------------
Old description:
> The following worksheet implements arbitrary mesh functions for 3d
> graphics objects that are triangulated. This could probably ought to go
> into sage somewhere:
>
> Arbitrary Mesh Functions
> system:sage
>
> {{{id=15|
> from sage.ext.fast_eval import fast_float
> ///
> }}}
>
> <p>Simple bisection root finder for a function of 3 variables.</p>
>
> {{{id=16|
> %hide
> %auto
> %cython
>
> cpdef find_root(f, double target, point_1, point_2, double epsilon=1e-4):
> """
> Given two 3d points point_1 and point_2, find a point (x,y,z) on the
> segment between them where abs(f(x,y,z)-target)<=epsilon.
>
> Returns (x,y,z) if the point lies on the segment between point_1 and
> point_2 and satisfies the above inequality, otherwise returns None.
>
> The code assumes the function f is continuous.
> """
> cdef double s0,s1,s2
> cdef double e0,e1,e2
> cdef double new_0,new_1,new_2
> cdef double target_0 = target-epsilon
> cdef double target_1 = target+epsilon
> cdef double val
> cdef int i
>
> cdef double min = f(*point_1)
> cdef double max = f(*point_2)
> s0,s1,s2 = point_1
> e0,e1,e2 = point_2
> if min>max:
> min,max=max,min
> s0,s1,s2,e0,e1,e2 = e0,e1,e2,s0,s1,s2
> # Check to see if one of the endpoints satisfies it
> if target_0<=min and min<=target_1:
> return (s0,s1,s2)
> if target_0<=max and max<=target_0:
> return (e0,e1,e2)
> if min>target_1 or max<target_0:
> return None
> i=0
> while True:
> if i>100:
> return None
> else:
> i+=1
> # Get half-way point
> new_0 = s0+(e0-s0)/2.0
> new_1 = s1+(e1-s1)/2.0
> new_2 = s2+(e2-s2)/2.0
> val = f(new_0, new_1, new_2)
> if val<target_0:
> s0, s1, s2 = new_0, new_1, new_2
> min = val
> elif target_1<val:
> e0, e1, e2 = new_0, new_1, new_2
> max=val
> else:
> return (new_0,new_1,new_2)
> ///
> }}}
>
> {{{id=24|
>
> ///
> }}}
>
> {{{id=41|
> def calculate_crossing(f,target,v0,v1,vertices, cache_dict):
> """
> Calculate, for an edge (v0,v1), where f is "close" to target. Use
> cache_dict to cache the values.
> """
> # Make a canonical ordering of the vertices since (v0,v1) is the same
> edge as (v1,v0)
> if v0>v1:
> v0,v1=v1,v0
> if (v0,v1) in cache_dict:
> return cache_dict[(v0,v1)]
> else:
> pt = find_root(f, target, vertices[v0], vertices[v1])
> cache_dict[(v0,v1)]=pt
> return pt
> ///
> }}}
>
> {{{id=34|
> %time
> var('x,y,z')
> p=parametric_plot((x,y,9-x^2-y^2), (x,-3,3), (y,-3,3), mesh=True)
> f=x^2-sin(x*y^2)+cos(z)
> ff=fast_float(f, 'x', 'y','z')
> p.triangulate()
> vertices=p.vertex_list()
>
> # I am still calculating the function on each vertex multiple times; that
> could be optimized
> # vertex_values=[ff(*v) for v in vertices]
>
> mesh=[]
> for target in [0,2,..,20]:
> cache={}
> for face in [f+[f[0]] for f in p.index_faces()]:
> # Adding the [0] takes care of the edge from vertices[0] to
> vertices[-1]
> pts = [calculate_crossing(ff,target,face[i],
> face[i+1],vertices=vertices, cache_dict=cache) for i in
> range(len(face)-1)]
> pts = [i for i in pts if i is not None]
> mesh+=[line([pt1,pt2],thickness=3, color='black') for pt1,pt2 in
> Subsets(pts, 2)]
> ///
>
> CPU time: 0.91 s, Wall time: 1.26 s
> }}}
>
> {{{id=39|
> p+sum(mesh)
> ///
> }}}
>
> {{{id=42|
> (p+sum(mesh)).show(viewer='tachyon')
> ///
> }}}
>
> {{{id=44|
>
> ///
> }}}
>
> {{{id=45|
>
> ///
> }}}
New description:
The attached worksheet implements and illustrates arbitrary mesh functions
for 3d graphics objects that are triangulated. This could probably ought
to go into sage somewhere.
--
Comment(by jason):
Boy, that was all mes
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/5511#comment:1>
Sage <http://sagemath.org/>
Sage - Open Source Mathematical Software: Building the Car Instead of
Reinventing the Wheel
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sage-trac" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/sage-trac?hl=en
-~----------~----~----~----~------~----~------~--~---
