I have made a small script that computes the basis of a Lagragian
finite element in 2D
and would like to integrate its basis functions.
Basically, its like this (this works)
f = 1 - 3*x - 3*y + 2*x**2 + 2*y**2 + 4*x*y
intf = integrate(f, (y, 0, 1))
intff = integrate(intf, (x, 0, 1))
When I compute the Lagrangian basis, the first basis function looks
like f, when I print it out. However, there is obviously some
difference between
f and the first basis function, because when I try to integrate I
get:
Traceback (most recent call last):
File "sympy_test.py", line 90, in <module>
intff = integrate(intf, (x, 0, 1))
File "/usr/lib/python2.5/site-packages/sympy/integrals/
integrals.py", line 139, in integrate
return integral.doit()
File "/usr/lib/python2.5/site-packages/sympy/integrals/
integrals.py", line 90, in doit
antideriv = self._eval_integral(function, x)
File "/usr/lib/python2.5/site-packages/sympy/integrals/
integrals.py", line 120, in _eval_integral
return risch_norman(f, x)
File "/usr/lib/python2.5/site-packages/sympy/integrals/risch.py",
line 182, in risch_norman
if not f.has(x):
File "/usr/lib/python2.5/site-packages/sympy/core/basic.py", line
387, in has
p = Basic.sympify(patterns[0])
File "/usr/lib/python2.5/site-packages/sympy/core/basic.py", line
251, in sympify
raise NotImplementedError('matrix support')
NotImplementedError: matrix support
The complete script is as follows:
from sympy import *
x, y = symbols('xy')
def integral_over_reference_square(f):
intf = integrate(f, (y, 0, 1))
intff = integrate(intf, (x, 0, 1))
return intff
def pol_space(order):
sum = 0
basis = []
coeff = []
for o1 in range(0,order+1):
for o2 in range(0,order+1):
if o1 + o2 <= order:
aij = Symbol("a_%d_%d" % (o1,o2))
sum += aij*pow(x, o1)*pow(y, o2)
basis.append(pow(x, o1)*pow(y, o2))
coeff.append(aij)
return sum, coeff, basis
def create_point_set(order):
h = 1.0/(order)
set = []
for i in range(0, order+1):
x = i*h
for j in range(0, order+1):
y = j*h
if x + y <= 1.0 + h/2.0:
set.append((x,y))
return set
def create_matrix(equations, coeffs):
A = zeronm(len(equations), len(equations))
i = 0; j = 0
for j in range(0, len(coeffs)):
c = coeffs[j]
for i in range(0, len(equations)):
e = equations[i]
d, r = div(e, c)
A[i,j] = d
return A
order = 2
pol, coeffs, basis = pol_space(order)
print "polynom ", pol
print "polynom ", coeffs
print "polynom ", basis
points = create_point_set(order)
print "points ", points
equations = []
for p in points:
ex = pol.subs(x, p[0]).subs(y, p[1])
equations.append(ex )
A = create_matrix(equations, coeffs)
Ainv = A.inv()
b = eye(len(equations))
x = Ainv*b
print x
basis = []
for i in range(0,len(equations)):
N = pol
for j in range(0,len(coeffs)):
N = N.subs(coeffs[j], x[j,i])
basis.append(N)
f = basis[0]
print f
intf = integrate(f, (y, 0, 1))
print intf
intff = integrate(intf, (x, 0, 1))
print intff
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"sympy" 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/sympy?hl=en
-~----------~----~----~----~------~----~------~--~---
