In [76]:

    

using Jacobi
using PyPlot
include("../src/HPFEM.jl");


    

    




WARNING: replacing module HPFEM

In [77]:

    

M= 10
Q=M
nel = 1
nnodes = nel + 1
a=-1
b= 1
idir = [1]#,nnodes]
lagr = HPFEM.Lagrange1d(M);
quad = HPFEM.QuadType(Q);
base = HPFEM.Basis1d(lagr, quad);
nodes = collect(linspace(a, b, nnodes));
lmap = HPFEM.locmap(base)
dof = HPFEM.DofMap1d(lmap, nnodes, idir);


    

In [78]:

    

elemento = [HPFEM.Element1d(e, nodes[e], nodes[e+1], base) for e = 1:nel]
solver = HPFEM.CholeskySC(dof, HPFEM.BBSymTri);


    

In [79]:

    

for e = 1:nel
    Ae = HPFEM.mass_matrix(base, elemento[e])
    HPFEM.add_local_matrix(solver, e, Ae)
end


    

In [80]:

    

fun(x) = sin(2*pi*x)
x = linspace(-1,1,101)
#plot(x,fun(x));


    

In [81]:

    

Fe = zeros(HPFEM.nmodes(lmap), nel)

for e = 1:nel
    fe = fun(elemento[e].x)
    HPFEM.add_rhs!(base, elemento[e], fe, sub(Fe, :, e))
end


    

In [82]:

    

HPFEM.solve!(solver, Fe)


    

    
Out[82]:





10x10 Array{Float64,2}:
 5.44287e-19  0.951057   0.587785  …   0.587785  -0.587785  -0.951057   
 0.0505368    0.965458   0.546149      0.546149  -0.627919  -0.934225   
 0.163397     0.988767   0.447694      0.447694  -0.712077  -0.887782   
 0.322178     0.999904   0.295796      0.295796  -0.817092  -0.800787   
 0.500755     0.977966   0.103661      0.103661  -0.9139    -0.668482   
 0.668482     0.9139    -0.103661  …  -0.103661  -0.977966  -0.500755   
 0.800787     0.817092  -0.295796     -0.295796  -0.999904  -0.322178   
 0.887782     0.712077  -0.447694     -0.447694  -0.988767  -0.163397   
 0.934225     0.627919  -0.546149     -0.546149  -0.965458  -0.0505368  
 0.951057     0.587785  -0.587785     -0.587785  -0.951057  -2.44929e-16

In [83]:

    

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

Ue = ϕ * Fe;


    

In [86]:

    

using PyPlot
x = [(1-ξ)*el.a/2 + (1+ξ)*el.b/2 for el in elemento]
maxerr = -1000000
for e = 1:nel
    uu = fun(x[e])
    err = maxabs(uu-Ue[:,e])
    if err > maxerr maxerr = err end
        
    plot(x[e], Ue[:,e], "r", x[e], uu, "b")
end
maxerr


    

    













    
Out[86]:





7.180034344855812e-12