# -*- coding: utf-8 -*- # 3.0 # # Finite Elements (1D) # import sympy as sp import numpy as np import matplotlib.pyplot as plt sp.init_printing() # # Building a basis Φ of test and ansatz functions # # Template for basis functions f0 # x = sp.Symbol('x') f0 = sp.Piecewise((0,x<-1),(x+1,x<0),(1-x,x<=1),(0,True)) print "f0:" sp.pprint(f0) lam_f0 = sp.lambdify((x), f0) plt.figure() xa = np.linspace(-2,2,200) f0a = [lam_f0(xae) for xae in xa] plt.plot(xa,f0a,color='brown',linewidth=3) plt.show() # # Now, specify mesh width h and position x0 # n = 3 h = sp.Matrix(1,n,[sp.Rational(1,4) for i in range(1,n+1)]) x0 = sp.Matrix(1,n,[sp.Rational(i,(n+1)) for i in range(1,n+1)]) print 'x0 =', sp.pretty(x0) print 'h =', sp.pretty(h) # # Define the basis functions # phi = [f0.subs(x,(x-x0[i])/h[i]) for i in range(0,n)] sp.pprint(phi[0]) # # Plot the basis functions # lam_phi = [sp.lambdify((x), phie) for phie in phi] xa = np.linspace(0,1,200) plt.figure() plt.hold(True) for lam_phie in lam_phi: phia = [lam_phie(xae) for xae in xa] plt.plot(xa,phia,linewidth=3) plt.hold(False) plt.show() # # Stiffness matrix and mass matrix # # Stiffness matrix A # # Stiffness matrix: $A_{ij} = \int_0^1 \frac{\partial \phi_i}{\partial x} \frac{\partial \phi_j}{\partial x} \mathrm{d}x$. # from scipy import integrate #A = [[sp.integrate(diff(phiI)*diff(phiJ),x) for phiI in phi] for phiJ in phi] # sympy 'sometimes' incorrectly integrates piecewise functions! A = [[integrate.quad(sp.lambdify(x,sp.diff(phiI)*sp.diff(phiJ)),0.,1.0)[0] for phiI in phi] for phiJ in phi] A = sp.Matrix(A) sp.pprint(A) # # Mass matrix B # # Mass matrix: $B_{ij} = \int_0^1 \phi_i \phi_j \mathrm{d}x$. # B = [[integrate.quad(sp.lambdify(x,phiI*phiJ),0.,1.0)[0] for phiI in phi] for phiJ in phi] B = sp.Matrix(B) sp.pprint(B) # # A very simple example # # We solve the (embarrassingly simple) equation $u = x(1-x)$ using the bilinear form $b(u,v) = \int_0^1 u(x) v(x) \mathrm{d}x $ # The result uh in Vh minimises the L2 norm of uh Kx 1 Kx . # bb = x*(1-x) b = sp.Matrix(n,1,[integrate.quad(sp.lambdify(x,phiI*bb),0.,1.)[0] for phiI in phi]) sp.pprint(b) # # Solution # u = B.LUsolve(b) sp.pprint(u) # # $B$ is the system matrix, $b$ the right hand side. # # The respective piecewise linear function $uu$ is obtained as a linear combination of the basis functions $bs$ # uu = u.dot(phi) lam_uu = sp.lambdify(x, uu) lam_bb = sp.lambdify(x, bb) xa = np.linspace(0,1,200) plt.figure() plt.hold(True) lam_uua = [lam_uu(xae) for xae in xa] lam_bba = [lam_bb(xae) for xae in xa] plt.plot(xa,lam_bba,linewidth=3) plt.plot(xa,lam_uua,linewidth=3) plt.hold(False) plt.show() # # The difference between $bb$ and the approximation $uu$ # dbbuu = [lam_bba[i]-lam_uua[i] for i in range(0,len(lam_bba))] print len(dbbuu), len(xa) plt.figure() plt.plot(xa,dbbuu,linewidth=3) plt.show() # # Poisson's equation in 1D # # For Poisson's equation, $A$ (the stiffness matrix) is the system matrix. # # We can retain $x(1-x)$ as the right hand side (thus solving the equation $-u^{''} = x(1-x)$ ), so $b$ is still our right hand side of the Linear System. # # The equation can still be solved exactly: # usol = -x**3/6 + x**4/12 + x/12 sp.pprint(-sp.diff(sp.diff(usol))) sp.pprint(usol.subs(x,0.0).evalf()) sp.pprint(usol.subs(x,1.0).evalf()) # # Use parabolic function as right-hand side: # bb = x*(1-x) b = sp.Matrix(n,1,[integrate.quad(sp.lambdify(x,phiI*bb),0.,1.)[0] for phiI in phi]) sp.pprint(b) # # The FE solution, however, is: # u = A.LUsolve(b) sp.pprint(u) # uu = u.dot(phi) lam_uu = sp.lambdify(x, uu) lam_usol = sp.lambdify(x, usol) xa = np.linspace(0,1,200) plt.figure() plt.hold(True) lam_uua = [lam_uu(xae) for xae in xa] lam_usol = [lam_usol(xae) for xae in xa] plt.plot(xa,lam_usol,linewidth=3) plt.plot(xa,lam_uua,linewidth=3) plt.hold(False) plt.show()