improved performance

Makefile

@@ -62,7 +62,7 @@ test:
 # usage: make test_file FILE=test/test_heat.jl
 # or even: watch make test_file FILE=test/test_heat.jl
 test_file:
-	julia -e 'using JuliaFEM.Test; run_test("$(FILE)"); print_test_statistics()'
+	timeout 30 julia -e 'using JuliaFEM.Test; run_test("$(FILE)"); print_test_statistics()'
 
 # usage: make test_function FILE=test/test_fields.jl FUNCTION=test_increment
 # or even: watch ...

src/directsolver.jl

@@ -67,9 +67,13 @@ function solve(K, f, C, g, ::Type{Val{:LDLt}})
     t0 = time()
 
     # make sure K is symmetric
+    K = Symmetric(K)
+#=
     s = maximum(abs(1/2*(K + K') - K))
     @assert s < 1.0e-6
     K = 1/2*(K + K')
+=#
+
     dim = size(K, 1)
 
     # make sure C is square

src/elasticity.jl

@@ -55,16 +55,18 @@ https://en.wikipedia.org/wiki/Hooke's_law
 """
 function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element::Element, ip::IntegrationPoint, time::Number; variation=nothing)
 
-    basis = element(ip, time)
 
-    u = element("displacement", ip, time, variation)
+    # u = element("displacement", ip, time, variation)
 
     r = zeros(Float64, problem.dim, length(element))
 
     # internal forces
     if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
-        dbasis = element(ip, time, Val{:grad})
-        gradu = element("displacement", ip, time, Val{:grad}, variation)
+        u = element("displacement", time, variation)
+        grad = element(ip, time, Val{:grad})
+#        gradu = element("displacement", ip, time, Val{:grad}, variation)
+        gradu = grad*u
+
         F = I + gradu # deformation gradient
 
         young = element("youngs modulus", ip, time)
@@ -77,22 +79,24 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
         E = 1/2*(F'*F - I)  # strain
         S = lambda*trace(E)*I + 2*mu*E
 
-        J = det(element, ip, time)
-        T = J^-1*F*S*F'
+        #J = det(element, ip, time)
+        #T = J^-1*F*S*F'
         #ip["cauchy stress"] = T
         #ip["gl strain"] = E
 
-        r += F*S*dbasis
+        r += F*S*grad
     end
 
     # external forces - volume load
     if haskey(element, "displacement load")
+        basis = element(ip, time)
         b = element("displacement load", ip, time)
         r -= b*basis
     end
 
     # external forces - surface traction force
     if haskey(element, "displacement traction force")
+        basis = element(ip, time)
         T = element("displacement traction force", ip, time)
         r -= T*basis
     end

src/elements.jl

@@ -118,14 +118,20 @@ end
 
 typealias VecOrIP Union{Vector, IntegrationPoint}
 
+function call(element::Element, field_name::ASCIIString, time::Real, variation=nothing)
+    return isa(variation, Void) ? element[field_name](time) : variation
+end
+
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, variation=nothing)
-    field = isa(variation, Void) ? element[field_name](time) : variation
+    field = element(field_name, time, variation)
+#   field = isa(variation, Void) ? element[field_name](time) : variation
     basis = get_basis(element)
     return basis(field, xi)
 end
 
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Number, ::Type{Val{:grad}}, variation=nothing)
-    field = isa(variation, Void) ? element[field_name](time) : variation
+#   field = isa(variation, Void) ? element[field_name](time) : variation
+    field = element(field_name, time, variation)
     basis = get_basis(element)
     geom = element["geometry"](time)
     return basis(geom, field, xi, Val{:grad})
@@ -194,7 +200,7 @@ end
 function LinAlg.det{E<:AbstractElement}(element::Element{E}, ip::IntegrationPoint, time::Number=0.0)
     X = element("geometry", time)
     dN = get_dbasis(E, ip.xi)
-    J = sum([dN[:,i]*X[i]' for i=1:length(X)])
+    J = sum([kron(dN[:,i], X[i]') for i=1:length(X)])
     m, n = size(J)
     return m == n ? det(J) : norm(J)
 end

src/equations.jl

@@ -152,8 +152,9 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
         end
 
         jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
+        residual_vector = -ForwardDiff.value(allresults)
         add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
-        add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))
+        add!(assembly.force_vector, gdofs, residual_vector)
     end
 end
 

src/fields.jl

@@ -290,7 +290,8 @@ end
 
 function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
     dbasis = basis(xi, Val{:grad})
-    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+#    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+    J = sum([kron(dbasis[:,i], geometry[i]') for i=1:length(geometry)])
     invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
     grad = invJ * dbasis
     return grad
@@ -298,7 +299,8 @@ end
 
 function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
     grad = call(basis, geometry, xi, Val{:grad})
-    gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
+#    gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
+    gradf = sum([kron(grad[:,i], values[i]') for i=1:length(values)])'
     return length(gradf) == 1 ? gradf[1] : gradf
 end
 
@@ -306,6 +308,10 @@ function Base.call(basis::CVTI, xi::Vector, time::Number)
     call(basis, xi)
 end
 
+function Base.(:*)(grad::Matrix{Float64}, field::DVTI)
+    return sum([kron(grad[:,i], field[i]') for i=1:length(field)])'
+end
+
 ### FIELDSET ###
 
 typealias FieldSet Dict{ASCIIString, Field}

src/sparse.jl

@@ -85,8 +85,14 @@ Example
 """
 function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix{Float64})
     n, m = size(data)
-    append!(A.I, repeat(dofs1, outer=[m]))
-    append!(A.J, repeat(dofs2, inner=[n]))
+    for i=1:n
+        for j=1:m
+            push!(A.I, dofs1[i])
+            push!(A.J, dofs2[j])
+        end
+    end
+#   append!(A.I, repeat(dofs1, outer=[m]))
+#   append!(A.J, repeat(dofs2, inner=[n]))
     append!(A.V, vec(data))
 end
 

src/types.jl

@@ -18,7 +18,7 @@ weight
 fields
     FieldSet what can be used to store internal variables, stress, strain, ...
 """
-type IntegrationPoint
+immutable IntegrationPoint
     xi :: Vector
     weight :: Float64
     fields :: Dict{ASCIIString, Field}

test/test_elasticity.jl

@@ -19,11 +19,11 @@ function test_elasticity_volume_load()
     free_dofs = [3, 4, 5, 6]
     solve!(problem, free_dofs, 0.0; max_iterations=10)
     disp = element("displacement", [1.0, 1.0], 0.0)
-    ip1 = last(element["integration points"])[1]
-    ip2 = last(element["integration points"])[2]
-    strain = ip1["gl strain"]
+    #ip1 = last(element["integration points"])[1]
+    #ip2 = last(element["integration points"])[2]
+    #strain = ip1["gl strain"]
     info("displacement at tip: $disp")
-    info("strain in first ip: $strain. ip coord = $(ip1.xi) and weight = $(ip1.weight)")
+    #info("strain in first ip: $strain. ip coord = $(ip1.xi) and weight = $(ip1.weight)")
     # verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu
     @test isapprox(disp[2], -8.77303119819776)
 end