In [ ]:

    

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


    

In [ ]:

    

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


    

In [ ]:

    

mp = [5 1 3;
      1 2 4]'


    

In [ ]:

    

a = 0.0
b = π
nodes = collect(linspace(a, b, nnodes))


    

In [ ]:

    

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


    

In [ ]:

    

Me = [HPFEM.mass_matrix(bas, el) for el in elems]


    

In [ ]:

    

A = zeros(5,5)
for e = 1:nel
    for i = 1:M
        ig = mp[i,e]
        for k = 1:M
            kg = mp[k,e]
            A[kg,ig] += Me[e][k,i]
        end
    end
end
A


    

In [ ]:

    

lft = HPFEM.DirichiletLift(A, [5])


    

In [ ]:

    

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

for e = 1:nel
    fe = cos(elems[e].x)
    HPFEM.add_rhs!(bas, elems[e], fe, view(Fe, :, e))
end

# Apply Dirichilet BCs:
Fe[1,1] = 1.0


    

In [ ]:

    

F = zeros(5)
for e = 1:nel
    for i = 1:M
        ig = mp[i,e]
        F[ig] += Fe[i,e]
    end
end
F


    

In [ ]:

    

Ahh = A[1:4,1:4]
Ahd = A[1:4,5:5]


    

In [ ]:

    

F5 = 1.0
Fh = F[1:4] - Ahd*F[5]


    

In [ ]:

    

F2 = copy(F)
F2[5] = 1.0
HPFEM.lift!(lft, F2)
F2[1:4]


    

In [ ]:

    

Uh = Ahh\Fh


    

In [ ]:

    

U = [Uh; F5]


    

In [ ]:

    

# Scatter
Ue = zeros(M, nel)
for e = 1:nel
    for i = 1:M
        ig = mp[i,e]
        Ue[i,e] = U[ig]
    end
end
Ue


    

In [ ]:

    

nξ = 101
ξ = collect(linspace(-1,1, nξ))
x = [(1-ξ)*el.a/2 + (1+ξ)*el.b/2 for el in elems]
ϕ = zeros(nξ, M)
for i = 1:M
    ϕ[:,i] = HPFEM.basis1d(bas, ξ, i)
end
Uh = ϕ*Ue;


    

In [ ]:

    

using PyPlot
for e = 1:nel
    plot(x[e], Uh[:,e], "r", x[e], cos(x[e]), "b")
end


    

In [ ]:

    

# Using static condensation.


    

In [ ]:

    




    

In [ ]: