Commits

Nathan Collier committed bd58d49

lecture 2, interpolation, plotting, p-fem, quadrature

  • Participants
  • Parent commits 30c3875

Comments (0)

Files changed (4)

File current/element.py

         self.x1 = x1
         self.LtoG = LtoG
         self.dx = x1-x0
-
+        self.p = len(LtoG)-1
+        
     def basis(self,x):
         if x < self.x0 or x > self.x1:
-            return [0.0, 0.0]
-        return [ (self.x1-x)/self.dx, (x-self.x0)/self.dx ]
+            return [0.0]*(self.p+1)
+        u = (x-self.x0)/self.dx
+        if self.p == 1:
+            return [ 1-u, u ]
+        if self.p == 2:
+            return [ (1-u)**2, 2*u*(1-u), u**2 ]
 
     def dbasis(self,x):
         if x < self.x0 or x > self.x1:
-            return [0.0, 0.0]
-        return [ -1.0/self.dx, 1.0/self.dx ]
+            return [0.0]*(self.p+1)
+        u = (x-self.x0)/self.dx
+        if self.p == 1:
+            return [ -1./self.dx, 1./self.dx ]
+        if self.p == 2:
+            return [ -2.*(1.-u)/self.dx, (2.-4.*u)/self.dx, 2.*u/self.dx ]
 
     def interp(self,x,d):
         N = self.basis(x)
         dN = self.dbasis(x)
         u = 0.
         du = 0.
-        for a in range(2):
+        for a in range(self.p+1):
             u += N[a]*d[self.LtoG[a]]
             du += dN[a]*d[self.LtoG[a]]
         return u,du
 
-    def plot(self,d,res=100):
+    def plot(self,d,res=10):
         xs = np.linspace(self.x0,self.x1,res)
         u = []
         du = []
             u.append(sol)
             du.append(der)
         plt.subplot(211)
-        plt.plot(xs,u,'-b')
+        plt.plot(xs,u,'-b',lw=2)
+        plt.xlabel('x')
+        plt.ylabel('u(x)')
         plt.subplot(212)
-        plt.plot(xs,du,'-g')
+        plt.plot(xs,du,'-g',lw=2)
+        plt.xlabel('x')
+        plt.ylabel("u'(x)")

File current/gauss.py

+def GetGauss1D(n):
+    """ Given the number of points desired, this function returns 2
+    lists. The first are the function evaluation gauss points for the
+    domain [-1,1]. The second list are the accompanying weights.
+    
+    Ex:
+    Xg,Wg = GetGauss1D
+    """
+
+    W=[]
+    X=[]
+
+    if n == 1: # integrates p <= 1 exactly
+        X.append(0.0) 
+        W.append(2.0) 
+    elif n == 2: # integrates p <= 3 exactly
+        X.append(-0.577350269189626)
+        X.append(-X[0])
+        W.append(1.0)
+        W.append(1.0)
+    elif n == 3: # integrates p <= 5 exactly
+        X.append(-0.774596669241483)
+        X.append(0.0)
+        X.append(-X[0]) 
+        W.append(0.555555555555556)
+        W.append(0.888888888888889)
+        W.append(W[0]) 
+    elif n == 4: # integrates p <= 7 exactly
+        X.append(-.86113631159405257524)
+        X.append(-.33998104358485626481)
+        X.append(.33998104358485626481)
+        X.append(.86113631159405257524)
+        W.append(.34785484513745385736)
+        W.append(.65214515486254614264)
+        W.append(.65214515486254614264)
+        W.append(.34785484513745385736)
+    elif n == 5: # integrates p <= 9 exactly
+        X.append(-.90617984593866399282)
+        X.append(-.53846931010568309105)
+        X.append(0.0)
+        X.append(-X[1])
+        X.append(-X[0])
+        W.append(.23692688505618908749)
+        W.append(.47862867049936646808)
+        W.append(.56888888888888888888)
+        W.append(W[1])
+        W.append(W[0])
+
+    return X,W
+
+def GetGauss2D(ng1d):
+    """ Two dimensional analog to GetGauss1D. Also returns the number
+    of points returned, that is ng1d*ng1d
+
+    Ex:
+    ng,Xg,Wg = GetGauss2D(ng1d)
+    """
+    Xt,Wt = GetGauss1D(ng1d)
+    ng = ng1d*ng1d
+    Xg=[]
+    Wg=[]
+    for i in range(len(Xt)):
+        for j in range(len(Xt)):
+            Xg.append([Xt[i],Xt[j]])
+            Wg.append(Wt[i]*Wt[j])
+
+    return ng,Xg,Wg
+
+def GetGauss3D(ng1d):
+    """ Three dimensional analog to GetGauss1D. Also returns the number
+    of points returned, that is ng1d*ng1d*ng1d
+
+    Ex:
+    ng,Xg,Wg = GetGauss3D(ng1d)
+    """
+    Xt,Wt = GetGauss1D(ng1d)
+    ng = ng1d*ng1d*ng1d
+    Xg=[]
+    Wg=[]
+    for i in range(len(Xt)):
+        for j in range(len(Xt)):
+            for k in range(len(Xt)):
+                Xg.append([Xt[i],Xt[j],Xt[k]])
+                Wg.append(Wt[i]*Wt[j]*Wt[k])
+
+    return ng,Xg,Wg
+
+def GetTriangleGauss(ng):
+    """ Given the number of points desired, this function returns 2
+    lists. The first are the function evaluation gauss points for the
+    2D triangle domain. The second list is the accompanying weights.
+    
+    Ex:
+    Xg,Wg = GetGauss1D
+    """
+
+    Wg=[]
+    Xg=[]
+
+    if ng == 3: # degree of precision = 2
+        a = 2.0/3.0
+        b = 1.0/6.0
+        c = 1.0/3.0 
+        Xg.append([a,b]) 
+        Wg.append(c) 
+        Xg.append([b,a]) 
+        Wg.append(c) 
+        Xg.append([b,b]) 
+        Wg.append(c) 
+    elif ng == 4:
+        a = 1.0/3.0
+        b = 0.520833333333333
+        Xg.append([a,a])
+        Wg.append(-0.5625)
+        Xg.append([0.6,0.2])
+        Wg.append(b)
+        Xg.append([0.2,0.6])
+        Wg.append(b)
+        Xg.append([0.2,0.2])
+        Wg.append(b)
+
+    return Xg,Wg

File current/poisson1d.py

 F[0] = g
 
 u = np.linalg.solve(K,F)
-plt.plot(np.linspace(0,1,n),u,'-o')
+
+#plt.plot(np.linspace(0,1,n),u,'-o')
+for e in E:
+    e.plot(u)
+    
 plt.show()
 
 

File current/poisson1d_p.py

+"""
+This code solves the Poisson equation:
+
+   -u'' = force
+   u(0) = g
+  u'(1) = h
+
+using the element point-of-view. With this code we introduced the idea
+of the `element' data structure.
+"""
+import numpy as np
+import pylab as plt
+from element import *
+from gauss import GetGauss1D
+
+def force(x):
+    return 1.
+
+# problem specific data
+g = 0. #  u(0) = g
+h = 0. # u'(1) = h
+
+# encode the mesh, a list of elements
+E = []
+E.append( element( 0.0,0.3, [0, 1] ) )
+E.append( element( 0.3,0.8, [1, 2, 3] ) )
+E.append( element( 0.8,1.0, [3, 4] ) )
+
+# number of basis functions
+n = E[-1].LtoG[-1]+1
+
+# get quadrature
+ng = 2
+Xg,Wg = GetGauss1D(ng)
+
+# left- and right-hand side
+K = np.zeros((n,n))
+F = np.zeros(n)
+for e in E:
+    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):
+            A = e.LtoG[a]
+            F[A] += N[a]*force(xq)*wq
+            for b in range(e.p+1):
+                B = e.LtoG[b]
+                K[A,B] += dN[a]*dN[b]*wq
+
+# Neumann condition
+F[-1] += h
+
+# Dirichlet condition
+F -= K[:,0]*g
+K[:,0] = K[0,:] = 0
+K[0,0] = 1
+F[0] = g
+
+u = np.linalg.solve(K,F)
+
+for e in E:
+    e.plot(u,res=100)    
+plt.show()
+