# FEMLectures

committed d1fde81

lecture 3, further abstraction, sparse linear algebra, supg

# Files changed (2)

`+"""`
`+This code solves the advection-dominated diffusion problem:`
`+`
`+   -u'' + adv*u' = 0`
`+            u(0) = 1`
`+            u(1) = 0`
`+`
`+It optionally includes stabilization using SUPG.`
`+"""`
`+import numpy as np`
`+import pylab as plt`
`+from element import *`
`+from scipy.sparse import linalg`
`+`
`+# encode the mesh, a list of elements`
`+n = 11`
`+dx = 1./(n-1.)`
`+`
`+# encode the mesh, a list of elements`
`+E = []`
`+for e in range(n-1):`
`+    E.append( element(e*dx,(e+1)*dx,[e,e+1]) )`
`+`
`+# problem specific data`
`+adv  = 100.`
`+Pe   = adv*dx*0.5`
`+supg = True`
`+ibc  = [0 ,-1]`
`+vbc  = [1.,0.]`
`+`
`+# Define problem physics`
`+def LHS(x,w,w_x,u,u_x):`
`+    rhs = w_x*u_x + adv*w*u_x`
`+    if supg:`
`+        tau = 0.5*dx/adv*(1./np.tanh(Pe)-1./Pe)`
`+        rhs += (adv*w_x) * tau * (adv*u_x)`
`+    return rhs`
`+`
`+def RHS(x,w,w_x): return 0.`
`+`
`+K,F = FormSystem(E,LHS,RHS)`
`+`
`+# Dirichlet condition`
`+for i,v in zip(ibc,vbc):`
`+    F -= K[:,i]*v`
`+    K[:,i] = K[i,:] = 0`
`+    K[i,i] = 1`
`+    F[i] = v`
`+`
`+# Solve`
`+u = linalg.spsolve(K.tocsr(),F)[0,:]`
`+`
`+# Plot`
`+for e in E:`
`+    e.plot(u,res=10)    `
`+plt.show()`
`+`

# current/element.py

`-import numpy as np`
`-import pylab as plt`
`+from numpy import zeros,linspace`
`+from gauss import GetGauss1D`
`+from scipy.sparse import lil_matrix`
` `
` class element:`
`     def __init__(self,x0,x1,LtoG):`
`         return u,du`
` `
`     def plot(self,d,res=10):`
`-        xs = np.linspace(self.x0,self.x1,res)`
`+        import pylab as plt`
`+        xs = linspace(self.x0,self.x1,res)`
`         u = []`
`         du = []`
`         for x in xs:`
`         plt.plot(xs,du,'-g',lw=2)`
`         plt.xlabel('x')`
`         plt.ylabel("u'(x)")`
`+`
`+def FormSystem(E,LHS,RHS,ng=1):`
`+    # number of basis functions`
`+    n = E[-1].LtoG[-1]+1`
`+`
`+    # get quadrature`
`+    Xg,Wg = GetGauss1D(ng)`
`+`
`+    # left- and right-hand side`
`+    K = lil_matrix((n,n))`
`+    F = zeros((n,1))`
`+    for e in E:`
`+        Ke = zeros((e.p+1,e.p+1))`
`+        Fe = zeros(e.p+1)`
`+        for q in range(ng):`
`+            xq = (Xg[q]+1.)*0.5*e.dx+e.x0`
`+            wq = Wg[q]*e.dx*0.5`
`+            N  = e.basis(xq)`
`+            dN = e.dbasis(xq)`
`+            for a in range(e.p+1):`
`+                Fe[a] += RHS(xq,N[a],dN[a])*wq`
`+                for b in range(e.p+1):`
`+                    Ke[a,b] += LHS(xq,N[a],dN[a],N[b],dN[b])*wq`
`+        for a in range(e.p+1):`
`+            A = e.LtoG[a]`
`+            F[A,0]   += Fe[a]`
`+            for b in range(e.p+1):`
`+                B = e.LtoG[b]`
`+                K[A,B] += Ke[a,b]`
`+    return K,F`
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.