All--
I'm a new user of sympy, and have been really impressed by how easy it is
to get started. :-) However, I recently came across a problem that's been
causing sympy to hang, and am not sure whether (a) it is an intractable
problem that sympy simply can't solve, (b) there is some hint I could give
sympy which would let it proceed, or (c) I just need to let it run for a
week. Specifically, I'm trying to compute the surface area of a quartic
triangle Bezier patch. The way this works is that you have fifteen control
points (X, a 15 x 3 matrix) and a corresponding fifteen basis functions (b,
a vector length fifteen, in terms of parameters s and t). The position for
a given s and t is given by X^T b. I'm trying to solve for a general
patch, X is a matrix of yet unknown constants. The surface area can be
found by evaluating a pretty intense double integral [1], which is what's
giving sympy trouble.
I've attached the script I'm using--I'd be infinitely grateful if someone
with more experience could offer any advice!
Best, and thanks to everyone who has contributed to this invaluable
software!
--Davis
P.S. I've added a guard to the beginning of the script because when I ran
this with 0.7.4 (default for Ubuntu 14.04 when installing with apt-get) it
resulted in unbounded memory use and crashed my machine--happily, updating
to the development branch fixed this.
[1] https://en.wikipedia.org/wiki/Parametric_surface#Surface_area
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/sympy.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/1488db59-35b6-4fe9-a602-501f304e82b6%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
"""
author: Davis Vigneault
"""
import sympy as sp
import copy
import pickle
import os
#x, y = sp.symbols('x y')
r, s, t = sp.symbols('r s t')
r1, r2, r3, r4 = sp.symbols('r1 r2 r3 r4')
s1, s2, s3, s4 = sp.symbols('s1 s2 s3 s4')
t1, t2, t3, t4 = sp.symbols('t1 t2 t3 t4')
def expand_powers(expr):
rep_r = [(r1,r**1),(r2,r**2),(r3,r**3),(r4,r**4)]
rep_s = [(s1,s**1),(s2,s**2),(s3,s**3),(s4,s**4)]
rep_t = [(t1,t**1),(t2,t**2),(t3,t**3),(t4,t**4)]
expr = expr.subs(rep_r)
expr = expr.subs(rep_s)
expr = expr.subs(rep_t)
return expr
def basis_functions():
position_basis = sp.Matrix([s4, # p0
t4, # p1
r4, # p2
4*r1*s3, # p0-
4*t1*s3, # p0+
4*s1*t3, # p1-
4*r1*t3, # p1+
4*t1*r3, # p2-
4*s1*r3, # p2+
6*s2*t2, # p01
6*t2*r2, # p12
6*r2*s2, # p20
12*s2*t1*r1, # p0m
12*s1*t2*r1, # p0m
12*s1*t1*r2]) # p0m
position_basis = expand_powers(position_basis)
return position_basis
def bezier_surfacearea():
"""
https://en.wikipedia.org/wiki/Parametric_surface#Surface_area
"""
# Create a directory where the output will be saved.
if not os.path.exists("data/"):
os.makedirs("data/")
# Vector of basis functions, length 15, in terms of r, s, and t.
bezier_weights = basis_functions()
# The parametric region is the unit triangle:
# 0 <= t <= 1
# 0 <= s <= 1 - t
# 0 <= r <= 1 - s - t
# R can be eliminated
bezier_weights = bezier_weights.subs( { r : 1 - s - t } )
print "Bezier Basis:"
sp.pprint(bezier_weights)
# Differentiate w.r.t. the remaining parameters, s and t.
bezier_basis_s = copy.deepcopy(bezier_weights).diff(s)
print "Partial derivative with respect to s:"
sp.pprint(bezier_basis_s)
bezier_basis_t = copy.deepcopy(bezier_weights).diff(t)
print "Partial derivative with respect to t:"
sp.pprint(bezier_basis_t)
# The position function is defined in terms of 15 control points.
X = map(lambda i: sp.Symbol('X_({0},{1})'.format(i / 3, i % 3)), xrange(45))
X = sp.Matrix(15,3,X)
print "Control point matrix:"
sp.pprint(X)
Rs = sp.Matrix(X.dot(bezier_basis_s))
print "X^T * bs:"
sp.pprint(Rs)
Rt = sp.Matrix(X.dot(bezier_basis_t))
print "X^T * bt:"
sp.pprint(Rt)
print "Calculating the norm of the cross product..."
cp = Rs.cross(Rt).norm()
with open("data/cross_product.txt", "w") as output_file:
pickle.dump(cp, output_file)
print "Computing first integral..."
int_1 = cp.integrate((s, 0, (1 - t)))
with open("data/inner_integral.txt", "w") as output_file:
pickle.dump(int_1, output_file)
print "Computing second integral..."
int_2 = int_1.integrate((t, 0, 1))
with open("data/outer_integral.txt", "w") as output_file:
pickle.dump(int_2, output_file)
def main():
if (int(sp.__version__[0]) < 1):
print 'ERROR:'
print """Running this script using sympy version 0.7.4 (the version
you get from apt-get on Ubuntu 14.04) resulted in unbounded
memory usage and a system crash. Updating to the development
branch (1.0.1-dev) fixed this issue, though I don\'t know
which specific patch has the fix. Please be aware of this
if you decide to uncomment this block."""
return
bezier_surfacearea()
if __name__ == '__main__':
main()
