#5511: arbitrary mesh functions
-------------------------+--------------------------------------------------
Reporter: jason | Owner: was
Type: enhancement | Status: new
Priority: major | Milestone: sage-3.4.1
Component: graphics | Keywords:
-------------------------+--------------------------------------------------
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|
///
}}}
--
Ticket URL: <http://trac.sagemath.org/sage_trac/ticket/5511>
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
-~----------~----~----~----~------~----~------~--~---
