2821.sagews

from sage.ext.fast_eval import fast_float

Simple bisection root finder for a function of 3 variables.

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

%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: 5.74 s, Wall time: 5.74 s

p+sum(mesh)

(p+sum(mesh)).show(viewer='tachyon')