notebook.community

Edit and run



In [3]:

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









    Out[3]:





HPFEM



In [4]:

    
nel = 10
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 [2]:

    
fun(x) = 1
#resp(x) = sin(2*pi*x)
λ(x) = x









    Out[2]:





λ (generic function with 1 method)



In [5]:

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









    Out[5]:





11-element Array{Float64,1}:
 -1.0
 -0.8
 -0.6
 -0.4
 -0.2
  0.0
  0.2
  0.4
  0.6
  0.8
  1.0



In [6]:

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



In [7]:

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



In [8]:

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

    
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] = airyai(a)
Fe[2,nel]= airyai(b)









    



LoadError: MethodError: `add_rhs!` has no method matching add_rhs!(::HPFEM.Basis1d{Float64,HPFEM.ModalC01d}, ::HPFEM.Element1d{Float64}, ::Int64, ::SubArray{Float64,1,Array{Float64,2},Tuple{Colon,Int64},2})
Closest candidates are:
  add_rhs!{T<:Number}(::HPFEM.GenBasis1d, ::HPFEM.Element1d{T<:Number}, !Matched::AbstractArray{T<:Number,1}, ::AbstractArray{T<:Number,1})
  add_rhs!{T<:Number}(!Matched::HPFEM.SEM1d{T<:Number}, ::HPFEM.Element1d{T<:Number}, !Matched::AbstractArray{T<:Number,1}, ::AbstractArray{T<:Number,1})
while loading In[11], in expression starting on line 2

 [inlined code] from In[11]:4
 in anonymous at no file:0



In [88]:

    
HPFEM.solve!(solver, Fe)









    Out[88]:





15x2 Array{Float64,2}:
  0.0  -0.0
 -0.0   0.0
 -0.0   0.0
 -0.0   0.0
 -0.0   0.0
 -0.0   0.0
 -0.0   0.0
  0.0   0.0
  0.0   0.0
  0.0   0.0
 -0.0   0.0
  0.0   0.0
  0.0   0.0
  0.0   0.0
  0.0   0.0



In [89]:

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

Ue = ϕ * Fe;



In [90]:

    
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



In [ ]:



In [ ]: