SyFi  0.3
fem_sympy.py
Go to the documentation of this file.
00001 
00002 from sympy import *
00003 
00004 
00005 x,y,z = symbols('xyz')
00006 
00007 class ReferenceSimplex: 
00008     def __init__(self, nsd): 
00009         self.nsd = nsd
00010         coords = []
00011         if nsd <= 3: 
00012             coords = symbols('xyz')[:nsd]
00013         else: 
00014             coords = []
00015             for d in range(0,nsd): 
00016                 coords.append(Symbol("x_%d" % d))
00017         self.coords = coords
00018 
00019     def integrate(self,f): 
00020         coords = self.coords
00021         nsd = self.nsd
00022 
00023         limit = 1 
00024         for p in coords:  
00025             limit -= p
00026 
00027         intf = f 
00028         for d in range(0,nsd): 
00029             p = coords[d]
00030             limit += p
00031             intf = integrate(intf, (p, 0, limit))  
00032         return intf
00033 
00034 def bernstein_space(order, nsd): 
00035     if nsd > 3: 
00036         raise RuntimeError("Bernstein only implemented in 1D, 2D, and 3D")
00037     sum = 0 
00038     basis = []
00039     coeff = []
00040 
00041     if nsd == 1: 
00042         b1, b2 = x, 1-x 
00043         for o1 in range(0,order+1): 
00044             for o2 in range(0,order+1): 
00045                 if o1 + o2  == order: 
00046                     aij = Symbol("a_%d_%d" % (o1,o2)) 
00047                     sum += aij*binomial(order,o1)*pow(b1, o1)*pow(b2, o2)
00048                     basis.append(binomial(order,o1)*pow(b1, o1)*pow(b2, o2))
00049                     coeff.append(aij)
00050 
00051 
00052     if nsd == 2: 
00053         b1, b2, b3 = x, y, 1-x-y 
00054         for o1 in range(0,order+1): 
00055             for o2 in range(0,order+1): 
00056                 for o3 in range(0,order+1): 
00057                     if o1 + o2 + o3 == order: 
00058                         aij = Symbol("a_%d_%d_%d" % (o1,o2,o3)) 
00059                         fac = factorial(order)/(factorial(o1)*factorial(o2)*factorial(o3))
00060                         sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)
00061                         basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3))
00062                         coeff.append(aij)
00063 
00064     if nsd == 3: 
00065         b1, b2, b3, b4 = x, y, z, 1-x-y-z 
00066         for o1 in range(0,order+1): 
00067             for o2 in range(0,order+1): 
00068                 for o3 in range(0,order+1): 
00069                     for o4 in range(0,order+1): 
00070                         if o1 + o2 + o3 + o4 == order: 
00071                             aij = Symbol("a_%d_%d_%d_%d" % (o1,o2,o3,o4)) 
00072                             fac = factorial(order)/(factorial(o1)*factorial(o2)*factorial(o3)*factorial(o4))
00073                             sum += aij*fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4)
00074                             basis.append(fac*pow(b1, o1)*pow(b2, o2)*pow(b3, o3)*pow(b4, o4))
00075                             coeff.append(aij)
00076 
00077 
00078     return sum, coeff, basis 
00079 
00080 def create_point_set(order, nsd): 
00081     h = Rational(1,order)
00082     set = []
00083 
00084     if nsd == 1: 
00085         for i in range(0, order+1):
00086             x = i*h 
00087             if x <= 1: 
00088                 set.append((x,y))
00089 
00090     if nsd == 2: 
00091         for i in range(0, order+1):
00092             x = i*h 
00093             for j in range(0, order+1):
00094                 y = j*h 
00095                 if x + y <= 1: 
00096                     set.append((x,y))
00097 
00098     if nsd == 3: 
00099         for i in range(0, order+1):
00100             x = i*h 
00101             for j in range(0, order+1):
00102                 y = j*h 
00103                 for k in range(0, order+1):
00104                     z = j*h 
00105                     if x + y + z  <= 1: 
00106                         set.append((x,y,z))
00107 
00108     return set
00109 
00110 
00111 
00112 def create_matrix(equations, coeffs): 
00113     A = zeronm(len(equations), len(equations))
00114     i = 0;  j = 0 
00115     for j in range(0, len(coeffs)):  
00116         c = coeffs[j]
00117         for i in range(0, len(equations)):  
00118             e = equations[i]
00119             d, r = div(e, c)
00120             A[i,j] = d 
00121     return A 
00122 
00123 
00124 
00125 class Lagrange: 
00126     def __init__(self,nsd, order): 
00127         self.nsd = nsd
00128         self.order = order
00129         self.compute_basis()
00130 
00131     def nbf(self): 
00132         return len(self.N)
00133 
00134     def compute_basis(self): 
00135         order = self.order 
00136         nsd = self.nsd
00137         N = []
00138         pol, coeffs, basis = bernstein_space(order, nsd)
00139         points = create_point_set(order, nsd)
00140 
00141         equations = []
00142         for p in points: 
00143             ex = pol.subs(x, p[0])
00144             if nsd > 1: 
00145                 ex = ex.subs(y, p[1]) 
00146             if nsd > 2: 
00147                 ex = ex.subs(z, p[2]) 
00148             equations.append(ex )
00149 
00150         A = create_matrix(equations, coeffs)
00151         Ainv = A.inv()
00152 
00153         b = eye(len(equations)) 
00154 
00155         xx = Ainv*b 
00156 
00157         for i in range(0,len(equations)):
00158             Ni = pol 
00159             for j in range(0,len(coeffs)):
00160                 Ni = Ni.subs(coeffs[j], xx[j,i]) 
00161             N.append(Ni)
00162 
00163         self.N = N                             
00164 
00165 
00166 
00167 
00168 
00169 
00170 
00171 if __name__ == '__main__': 
00172     t = ReferenceSimplex(2)
00173     x,y,z = symbols('xyz')
00174     f = x+y
00175 
00176     print "integral of f = x+y ", t.integrate(f)
00177     print ""
00178     print "linear bernstein in 1D ", bernstein_space(1,1)
00179     print "linear bernstein in 2D ", bernstein_space(1,2)
00180     print "linear bernstein in 3D ", bernstein_space(1,3)
00181     print ""
00182 
00183     
00184     print "basis functions of a second order Lagrange element"  
00185     fe = Lagrange(2,2)
00186     for i in range(0, fe.nbf()):
00187         print "N[%d]= %s " % (i,fe.N[i])
00188 
00189 
00190     print "The Jacobian matrix of Fi = u * Ni " 
00191 
00192     u = 0 
00193     #compute u = sum_i u_i N_i 
00194     us = []
00195     for i in range(0, fe.nbf()):
00196         ui = Symbol("u_%d" % i)
00197         us.append(ui)
00198         u += ui*fe.N[i]
00199 
00200 
00201     J = zeronm(fe.nbf(), fe.nbf()) 
00202     for i in range(0, fe.nbf()):
00203         Fi = u*fe.N[i]
00204         for j in range(0, fe.nbf()):
00205             uj = us[j] 
00206             integrands = diff(Fi, uj) 
00207             J[j,i] = t.integrate(integrands) 
00208 
00209 
00210     print J
00211 
00212 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines