updated JuliaFEM.jl

src/JuliaFEM.jl

@@ -10,7 +10,7 @@ module JuliaFEM
 #@Logging.configure(level=DEBUG)
 #using Lexicon
 
-import Base: +, -, /, *, push!, convert, getindex, length, similar, call, vec, endof
+import Base: +, -, /, *, push!, convert, getindex, setindex!, length, similar, call, vec, endof
 
 #importall Base
 
@@ -55,6 +55,15 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n)
     [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
 end
 
+function Base.resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
+    (n == A.n) && (m == A.m) && return
+    @assert n >= A.n
+    @assert m >= A.m
+    append!(A.colptr, A.colptr[end]*ones(Int, m-A.m))
+    A.n = n
+    A.m = m
+end
+
 # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
 include("fields.jl")
 #include("basis.jl")  # interpolation of discrete fields

src/elements.jl

@@ -148,6 +148,11 @@ function call(element::Element, field_name::ASCIIString, time::Number)
     return element[field_name](time)
 end
 
+function get_dbasis{E<:AbstractElement}(::Type{E}, xi::Vector)
+    basis(xi) = vec(get_basis(E, xi))
+    return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)'
+end
+
 function get_basis{E}(element::Element{E}, ip::IntegrationPoint)
     return get_basis(E, ip.xi)
 end

src/lagrange.jl

@@ -9,16 +9,22 @@ abstract CG <: AbstractElement
 Given polynomial P and coordinates of reference element, calculate
 Lagrange basis functions
 """
-function calculate_lagrange_basis(P, X)
+function calculate_lagrange_basis_coefficients(P, X)
     dim, nbasis = size(X)
     A = zeros(nbasis, nbasis)
     for i=1:nbasis
         A[i,:] = P(X[:, i])
     end
-    invA = inv(A)'
-    basis(xi) = (invA*P(xi))'
-    dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))'
-    basis, dbasisdxi
+#   invA = inv(A)'
+#   basis(xi) = (invA*P(xi))'
+#   dbasisdxi(xi) = (ForwardDiff.jacobian((xi) -> invA*P(xi), xi, cache=autodiffcache))'
+#   basis, dbasisdxi
+#   info(inv(A))
+#   info(P([0.0, 0.0]))
+#   invA = inv(A)'
+#   basis(xi) = invA*P(xi)
+#   return basis
+    return inv(A)'
 end
 
 """
@@ -32,20 +38,31 @@ macro create_lagrange_element(element_name, element_description, X, P)
     eltype = esc(element_name)
     quote
         global get_basis, get_dbasis
-        basis, dbasis = calculate_lagrange_basis($P, $X)
+        #basis, dbasis = calculate_lagrange_basis($P, $X)
+        C = calculate_lagrange_basis_coefficients($P, $X)
+        basis(xi) = C*$P(xi)
+#       dbasis = ForwardDiff.jacobian(basis)
+
         abstract $eltype <: CG
-        function get_basis(::Type{$eltype}, xi::Vector{Float64})
-            return basis(xi)
+
+        function get_basis(::Type{$eltype}, xi::Vector)
+            return basis(xi)'
         end
+
+#=
         function get_dbasis(::Type{$eltype}, xi::Vector{Float64})
-            return dbasis(xi)
+            return dbasis(xi)'
         end
+=#
+
         function $eltype(args...)
             return Element{$eltype}(args...)
         end
+
         function Base.size(::Type{$eltype})
             return Base.size($X)
         end
+
     end
 end
 

test/test_directsolver.jl

@@ -11,6 +11,7 @@ using JuliaFEM: PlaneStressElasticityProblem, DirichletProblem
 using JuliaFEM: DirectSolver
 
 function test_solver_multiple_dirichlet_bc()
+
     N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
 
     e1 = Quad4([1, 2, 4, 3])
@@ -19,7 +20,9 @@ function test_solver_multiple_dirichlet_bc()
     e1["poissons ratio"] = 0.25
     b1 = Seg2([3, 4])
     b1["geometry"] = Vector[N[3], N[4]]
-    b1["displacement traction force"] = Vector[[0.0, -100.0], [0.0, -100.0]]
+    b1["displacement traction force"] = (
+        0.0 => Vector[[0.0,    0.0], [0.0,    0.0]],
+        1.0 => Vector[[0.0, -100.0], [0.0, -100.0]])
 
     problem = PlaneStressElasticityProblem()
     push!(problem, e1)
@@ -56,14 +59,60 @@ function test_solver_multiple_dirichlet_bc()
 
     # launch solver
     norm = solver(0.0)
-
-    disp = e1("displacement", [1.0, 1.0], 0.0)
+    norm = solver(1.0)
+    disp = e1("displacement", [1.0, 1.0], 1.0)
     info("displacement at tip: $disp")
     @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
 
 end
 #test_solver_multiple_dirichlet_bc()
 
+
+function test_solver_no_convergence()
+
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
+
+    e1 = Quad4([1, 2, 4, 3])
+    e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
+    e1["youngs modulus"] = 900.0
+    e1["poissons ratio"] = 0.25
+    b1 = Seg2([3, 4])
+    b1["geometry"] = Vector[N[3], N[4]]
+    b1["displacement traction force"] = Vector[[100.0, 100.0], [100.0, 100.0]]
+
+    problem = PlaneStressElasticityProblem()
+    push!(problem, e1)
+    push!(problem, b1)
+
+    # boundary elements for dirichlet dx=0
+    dx = Seg2([1, 3])
+    dx["geometry"] = Vector[N[1], N[3]]
+    dx["displacement 1"] = 0.0
+
+    # boundary elements for dirichlet dy=0
+    dy = Seg2([1, 2])
+    dy["geometry"] = Vector[N[1], N[2]]
+    dy["displacement 2"] = 0.0
+  
+    problem2 = DirichletProblem("displacement", 2)
+    push!(problem2, dx)
+
+    problem3 = DirichletProblem("displacement", 2)
+    push!(problem3, dy)
+
+    solver = DirectSolver()
+    solver.max_iterations = 1
+    push!(solver, problem)
+    push!(solver, problem2)
+    push!(solver, problem3)
+
+    # launch solver
+    iterations, status = solver(0.0)
+
+    @test status == false
+
+end
+
 function test_solver_multiple_bodies_multiple_dirichlet_bc()
     N = Vector[
         [0.0, 0.0], [1.0, 0.0],
@@ -129,5 +178,4 @@ end
 
 # test_solver_multiple_bodies_multiple_dirichlet_bc()
 
-
 end