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)


    

    
Out[3]:





resp (generic function with 1 method)

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
    Ae = HPFEM.mass_matrix(bas, elems[e])
    Se = HPFEM.stiff_matrix(bas,elems[e])
    Ae = Ae + Se
    HPFEM.add_local_matrix(solver, e, Ae)
end


    

    




LoadError: MethodError: `*` has no method matching *(::Float64, ::Type{AbstractArray{Float64,1}})
Closest candidates are:
  *(::Any, ::Any, !Matched::Any, !Matched::Any...)
  *(::Float64, !Matched::Float64)
  *(::Real, !Matched::Complex{Bool})
  ...
while loading In[7], in expression starting on line 1

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)


    

    




LoadError: KeyError: 1 not found
while loading In[9], in expression starting on line 1

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





6.442340250271481

In [ ]: