View source

 # Resolution of Poisson 1D using FEM weak form
 
 from numpy import *
 
 # Problem definition
 x0=0.0
 xL=15.0
 Nx=101
 fi0=3.0   # Dirichlet condition
 qL=13.0   # Neumann condition
 Q0=5.0    # Heat load
 km=1.0     # material
 
 # definition of nodes and elements
 long=xL-x0              # length of domain
 nnodes=5                # number of total nodes
 nnod=2                  # number of nodes for each element
 nelems=nnodes-1         # number of elements
 interv=long/(nnodes-1)  # size of each element
 xi=zeros(nnodes,float)
 xi[0]=x0
 for i in range(1,nnodes):
     xi[i]=xi[i-1]+interv
 
 f=zeros((nelems,nnod),float)
 K=zeros((nelems,nnod,nnod),float)
 # obtaining elementary matrix components
 for i in range(0,nelems):
     for j in range(0,nnod):
         f[i,j]=Q0*interv/2
         for k in range(0,nnod):
             K[i,j,k]=pow(-1,j+k)*(km/interv)
         
 # assembly of global matrix components
 fg=zeros(nnodes,float)
 Kg=zeros((nnodes,nnodes),float)
 for i in range(0,nelems):
     for j in range(0,nnod):
         fg[i+j]=fg[i+j]+f[i,j]
         for k in range(0,nnod):
             Kg[i+j,i+k]=Kg[i+j,i+k]+K[i,j,k]
         
 # assembly of boundary conditions
 fg[nnodes-1]=fg[nnodes-1]-qL
 
 # Applying Dirichlet condition (remove the equations related to thr first unknown)
 Kf=Kg[1:,1:]
 cfi=Kg[1:,0]*fi0
 ff=fg[1:]-cfi
 
 fi=inner(linalg.inv(Kf),ff)
 
 fit=hstack((fi0,fi))
 
 # exact solution
 realFi=-(Q0/(2*km))*pow(xi,2)+(-qL+Q0*(xL-x0))*xi/km+fi0;
 
 print "end of the programm"
 







:[[Media:poisson1Dpython.zip |''Download the '''Source code''''']]

:[[Poisson's_Equation#Domain_Discretisation:_.22local.22_Shape_Functions |''Back to the theory'']]


[[Category: Python Code]]