Elementos finitos linear 1D

In [148]:

    

ϕ₁(x, a, b)  = (b-x)/(b-a)
ϕ₂(x, a, b)  = (x-a)/(b-a)
mass_matrix(a, b) = [(b-a)/3 (b-a)/6; (b-a)/6 (b-a)/3]


    

    
Out[148]:





mass_matrix (generic function with 1 method)

In [149]:

    

mass_matrix(0,1)


    

    
Out[149]:





2x2 Array{Float64,2}:
 0.333333  0.166667
 0.166667  0.333333

In [193]:

    

a = -1.0
b = 1.0
xn = [linspace(a,b, 5);]
Ndof = size(xn,1);


    

In [194]:

    

Nel = size(xn,1)-1
dof_map = zeros(Int, 2, Nel)
for i = 1:Nel
    dof_map[1,i] = i
    dof_map[2,i] = i+1
end
dof_map


    

    
Out[194]:





2x4 Array{Int64,2}:
 1  2  3  4
 2  3  4  5

In [195]:

    

M = zeros(Ndof, Ndof)

for e = 1:Nel
    Me = mass_matrix(xn[e], xn[e+1])
    for i = 1:2
        ig = dof_map[i,e]
        for k = 1:2
            kg = dof_map[k,e]
            M[kg,ig] += Me[k,i]
        end
    end
end
M


    

    
Out[195]:





5x5 Array{Float64,2}:
 0.166667   0.0833333  0.0        0.0        0.0      
 0.0833333  0.333333   0.0833333  0.0        0.0      
 0.0        0.0833333  0.333333   0.0833333  0.0      
 0.0        0.0        0.0833333  0.333333   0.0833333
 0.0        0.0        0.0        0.0833333  0.166667 

In [196]:

    

fun(x) = sin(π*x)


    

    
Out[196]:





fun (generic function with 1 method)

In [197]:

    

fe = fun(xn)


    

    
Out[197]:





5-element Array{Float64,1}:
 -1.22465e-16
 -1.0        
  0.0        
  1.0        
  1.22465e-16

In [198]:

    

F = zeros(Ndof)
for e = 1:Nel
    Me = mass_matrix(xn[e], xn[e+1])
    Fe = Me * fun(xn[e:(e+1)])
    for i = 1:2
        ig = dof_map[i,e]
        F[ig] += Fe[i]
    end
end
        

# Calculando a integral de maneira mais exata:
using Jacobi
Q = 10
z = zgj(Q)
w = wgj(z)

F2 = zeros(Ndof)
for e = 1:Nel
    Fe1 = 0.0
    Fe2 = 0.0
    a1 = xn[e]
    b1 = xn[e+1]
    x = ( (1-z)*a1 + (1+z)*b1 ) / 2
    J = (b1-a1) / 2
    for q = 1:Q
        f = fun(x[q])
        Fe1 += J * w[q] * f * ϕ₁(x[q], a1, b1)
        Fe2 += J * w[q] * f * ϕ₂(x[q], a1, b1)
    end
    i1 = dof_map[1,e]
    i2 = dof_map[2,e]
    F2[i1] += Fe1
    F2[i2] += Fe2
end


    

In [199]:

    

u = M\F


    

    
Out[199]:





5-element Array{Float64,1}:
 -8.32667e-17
 -1.0        
 -4.48359e-17
  1.0        
  9.61532e-17

In [200]:

    

hcat(F,F2)


    

    
Out[200]:





5x2 Array{Float64,2}:
 -0.0833333  -0.115668   
 -0.333333   -0.405285   
  0.0        -1.38778e-17
  0.333333    0.405285   
  0.0833333   0.115668   

In [201]:

    

u2 = M\F2


    

    
Out[201]:





5-element Array{Float64,1}:
 -0.0983749  
 -1.19126    
  8.96719e-17
  1.19126    
  0.0983749  

In [202]:

    

using PyPlot
xx = linspace(a, b, 201)
plot(xx, fun(xx))
plot(xn, u, "rs-")
plot(xn, u2, "bo--")


    

    













    
Out[202]:





1-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7f87d47b5bd0>

In [203]:

    

e1 = u - fe
e2 = u2 - fe

plot(xn, e1, "rs-", xn, e2, "bo--")


    

    













    
Out[203]:





2-element Array{Any,1}:
 PyObject <matplotlib.lines.Line2D object at 0x7f87d46eb4d0>
 PyObject <matplotlib.lines.Line2D object at 0x7f87d46eb710>

In [ ]:

    




    

In [ ]:

    




    

In [ ]: