In [46]:

    

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


    

    




WARNING: replacing module HPFEM






    
Out[46]:





HPFEM

In [47]:

    

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


    

In [48]:

    

fun(x) = airyai(x)


    

    
Out[48]:





fun (generic function with 1 method)

In [49]:

    

a = -25
b = 15.0
nodes = collect(linspace(a, b, nnodes))


    

    
Out[49]:





11-element Array{Float64,1}:
 -25.0
 -21.0
 -17.0
 -13.0
  -9.0
  -5.0
  -1.0
   3.0
   7.0
  11.0
  15.0

In [ ]:

    




    

In [50]:

    

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


    

In [51]:

    

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


    

In [52]:

    

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


    

In [53]:

    

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] = fun(a)
#Fe[2,2] = cos(b)
Fe


    

    
Out[53]:





100x10 Array{Float64,2}:
  0.163527      0.0332303     0.0580859    …   2.48647e-7    1.15841e-12
 -0.0280055    -0.066778      0.0725187        2.43277e-8    9.12603e-14
 -0.00436763   -0.00159089   -0.000507684      2.00791e-8    7.8086e-14 
 -0.00412942    0.00424068   -0.0073788       -2.75292e-8   -1.14758e-13
 -0.0104873    -0.00839867    0.0162055        2.37517e-8    1.0932e-13 
 -0.0140066    -0.0138158    -0.000446388  …  -1.54568e-8   -8.06286e-14
  0.00337871   -0.0116266     0.0509802        8.08848e-9    4.89568e-14
 -0.0119664    -0.0437976     0.0209714       -3.51851e-9   -2.5259e-14 
  0.0413794     0.0280696    -0.0561988        1.29867e-9    1.12974e-14
  0.0272303     0.0439412    -0.0105161       -4.12287e-10  -4.44276e-15
 -0.0369904    -0.0172684     0.0215409    …   1.13599e-10   1.55268e-15
 -0.0137093    -0.0161907     0.00177091      -2.73094e-11  -4.86308e-16
  0.0143774     0.00548012   -0.00444534       5.73746e-12   1.37419e-16
  ⋮                                        ⋱                            
  2.47783e-17   6.53973e-18   1.14616e-17      7.50933e-24   3.31373e-29
  8.27763e-18  -2.19526e-17   6.75763e-18     -6.88018e-24  -3.38135e-29
 -1.1793e-17   -1.79982e-17   2.10575e-17  …   3.58107e-24   2.28987e-29
 -2.08595e-17  -4.34219e-17   2.51728e-17      1.22259e-24  -2.94042e-30
  1.48512e-17   1.59208e-17  -2.68192e-18     -6.06217e-24  -1.80197e-29
  1.08935e-17   1.33647e-17  -1.07804e-17      9.69137e-24   3.64282e-29
  1.15656e-17  -1.61384e-17  -4.99707e-18     -1.15744e-23  -4.84434e-29
  3.24215e-17   3.17638e-17  -4.20202e-18  …   1.23923e-23   5.59837e-29
 -5.3441e-17   -1.5671e-17    3.46299e-18     -1.27206e-23  -6.32333e-29
 -2.33079e-17   1.54134e-17  -1.55116e-17      1.22594e-23   7.20786e-29
  8.89374e-18   1.48173e-17  -9.3991e-18      -1.02037e-23  -7.7626e-29 
 -3.57841e-17  -2.12734e-17   7.51938e-18      7.19834e-24   7.39721e-29

In [54]:

    

HPFEM.solve!(solver, Fe)


    

    
Out[54]:





100x10 Array{Float64,2}:
  0.163527      0.226358     -0.105262     …   7.49213e-7    4.22626e-12
  0.226358     -0.105262      0.17151          4.22626e-12   2.29492e-18
 -1.15591      -0.41361      -0.00283743      -1.83819e-6   -1.08043e-11
 -0.401481      0.454021     -0.637809         1.03182e-6    6.56383e-12
 -0.692728     -0.461711      0.777423        -5.40114e-7   -3.84659e-12
 -0.711978     -0.545923     -0.169829     …   2.46084e-7    2.02222e-12
  0.181518     -0.284932      1.30745         -9.68228e-8   -9.43949e-13
 -0.465109     -1.22927       0.472043         3.29771e-8    3.91789e-13
  1.14404       0.646607     -1.04267         -9.76065e-9   -1.4523e-13 
  0.665163      0.87818      -0.135337         2.51761e-9    4.83153e-14
 -0.809062     -0.331075      0.310784     …  -5.66103e-10  -1.44874e-14
 -0.247818     -0.249372      0.00896027       1.10579e-10   3.92114e-15
  0.258741      0.0883759    -0.0519506       -1.85778e-11  -9.50319e-16
  ⋮                                        ⋱                            
  4.9169e-14    1.50407e-13   3.76602e-13     -1.21727e-15   5.86593e-18
  5.93831e-14   4.82114e-14   1.49644e-13      1.21757e-15  -5.92828e-18
  3.25434e-14   1.22784e-13   3.1801e-13   …  -1.02511e-15   4.93992e-18
  4.23828e-14   3.28834e-14   1.29258e-13      1.02525e-15  -4.99189e-18
  2.75752e-14   9.99686e-14   2.58377e-13     -8.28665e-16   3.99326e-18
  3.69288e-14   3.43128e-14   1.00572e-13      8.2869e-16   -4.03484e-18
  1.73125e-14   6.37684e-14   1.9889e-13      -6.27934e-16   3.02595e-18
  2.93862e-14   3.34038e-14   7.91311e-14  …   6.27886e-16  -3.05714e-18
 -1.0339e-14    3.58602e-14   1.39024e-13     -4.22914e-16   2.03798e-18
 -8.25728e-16   1.7264e-14    5.52475e-14      4.22838e-16  -2.05877e-18
  2.21351e-15   2.35922e-14   6.70037e-14     -2.13603e-16   1.02933e-18
 -1.41415e-14  -4.50334e-15   3.41394e-14      2.13543e-16  -1.03973e-18

In [55]:

    

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

Ue = ϕ * Fe;


    

In [56]:

    

using PyPlot
x = [(1-ξ)*el.a/2 + (1+ξ)*el.b/2 for el in elems]
maxerr = 0.0
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")
    plot(x[e], uu)
end
maxerr


    

    













    
Out[56]:





7.046585537295869e-13

In [ ]: