In [1]:

    

include("../src/HPFEM.jl")


    

    
Out[1]:





HPFEM

In [2]:

    

nel = 2
nnodes = nel + 1
idir = [1,nnodes]
M = 15
Q = M+2
bas = HPFEM.Basis1d(M,Q)
lmap = HPFEM.locmap(bas)
dof = HPFEM.DofMap1d(lmap, nnodes, idir);


    

In [3]:

    

fun(x) = (1 + 4*pi^2)*sin(2*pi*x)
resp(x) = sin(2*pi*x)
λ(x) = x


    

    
Out[3]:





λ (generic function with 1 method)

In [ ]:

    




    

In [4]:

    

a = -1
b = 1
nodes = collect(linspace(a, b, nnodes))


    

    
Out[4]:





3-element Array{Float64,1}:
 -1.0
  0.0
  1.0

In [5]:

    

elems = [HPFEM.Element1d(e, nodes[e], nodes[e+1], bas) for e = 1:nel];


    

In [6]:

    

solver = HPFEM.CholeskySC(dof, HPFEM.BBMatrix);


    

In [7]:

    

for e = 1:nel
    x  = elems[e].x
    lambda = λ(x)
    Ae = HPFEM.mass_matrix(bas, elems[e],lambda)
    Se = HPFEM.stiff_matrix(bas,elems[e],lambda)
    Ae = Ae + Se
    HPFEM.add_local_matrix(solver, e, Ae)
end


    

In [8]:

    

Fe = zeros(HPFEM.nmodes(lmap), nel)
for e = 1:nel
    fe = fun(elems[e].x)
    HPFEM.add_rhs!(bas, elems[e], fe, sub(Fe, :, e))
end

# Apply Dirichilet BCs:
Fe[1,1] = resp(a)
Fe[2,nel]= resp(b)
Fe = Fe


    

    
Out[8]:





15x2 Array{Float64,2}:
  2.44929e-16   6.44234    
 -6.44234      -2.44929e-16
 -2.81676e-16   2.94903e-16
 -3.91647      -3.91647    
  1.58294e-16  -2.31586e-16
  1.6677        1.6677     
  1.81279e-16   1.06685e-16
 -0.196584     -0.196584   
  1.94289e-16   1.249e-16  
  0.0109091     0.0109091  
  1.82146e-17  -5.20417e-17
 -0.000353398  -0.000353398
 -2.10769e-16  -2.56739e-16
  7.51952e-6    7.51952e-6 
 -3.70363e-16   9.06393e-17

In [9]:

    

HPFEM.solve!(solver, Fe)


    

    
Out[9]:





15x2 Array{Float64,2}:
   2.44929e-16    4.64442    
   4.64442       -2.44929e-16
  10.9548        11.5856     
 -22.0226       -21.8509     
  -4.49119       13.1692     
   0.626251      -7.32845    
   0.167009       6.17237    
  -0.0364096     -5.27026    
  -0.0127435      4.27419    
  -0.00219567    -3.42455    
  -0.000946069    2.69396    
  -0.000461664   -2.0458     
  -0.000197058    1.46283    
  -8.17722e-5    -0.933174   
  -2.93503e-5     0.44793    

In [10]:

    

nξ = 101
ξ = collect(linspace(-1,1,nξ));
ϕ = zeros(nξ, M)
for i = 1:M
    ϕ[:,i] = bas(ξ, i)
end

Ue = ϕ * Fe;


    

In [11]:

    

using PyPlot
x = [(1-ξ)*el.a/2 + (1+ξ)*el.b/2 for el in elems]
maxerr = 0.0

for e = 1:nel
    uu = resp(x[e])
    err = maxabs(uu-Ue[:,e])
    if err > maxerr maxerr = err end
    plot(x[e], Ue[:,e])
    plot(x[e], uu, "b")
end
maxerr


    

    













    
Out[11]:





9.727360206464393

In [ ]:

    




    

In [ ]: