In [87]:

    

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


    

    




WARNING: replacing module HPFEM

In [88]:

    

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 [89]:

    

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


    

In [90]:

    

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


    

In [91]:

    

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


    

In [92]:

    

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 [93]:

    

HPFEM.solve!(solver, Fe)


    

    
Out[93]:





10x1 Array{Float64,2}:
  5.44287e-18
  0.484318   
  0.997513   
 -0.138257   
 -0.861633   
  0.861633   
  0.138257   
 -0.997513   
 -0.484318   
 -2.44929e-16

In [94]:

    

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

Ue = ϕ * Fe;


    

In [95]:

    

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[95]:





0.011170062638266476