From 4b07da194afc88cd5431a00039bcef43a9aca8d5 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Fri, 4 Dec 2015 07:27:40 +0200 Subject: [PATCH] improved performance --- Makefile | 2 +- src/directsolver.jl | 4 ++++ src/elasticity.jl | 18 +++++++++++------- src/elements.jl | 12 +++++++++--- src/equations.jl | 3 ++- src/fields.jl | 10 ++++++++-- src/sparse.jl | 10 ++++++++-- src/types.jl | 2 +- test/test_elasticity.jl | 8 ++++---- 9 files changed, 48 insertions(+), 21 deletions(-) diff --git a/Makefile b/Makefile index 007b946..c6f5a5f 100644 --- a/Makefile +++ b/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 ... diff --git a/src/directsolver.jl b/src/directsolver.jl index 5ca8b07..69ae0b3 100644 --- a/src/directsolver.jl +++ b/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 diff --git a/src/elasticity.jl b/src/elasticity.jl index 0afdc0d..a767b43 100644 --- a/src/elasticity.jl +++ b/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 diff --git a/src/elements.jl b/src/elements.jl index e9693e2..fbdfdf4 100644 --- a/src/elements.jl +++ b/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 diff --git a/src/equations.jl b/src/equations.jl index d816cd7..186a538 100644 --- a/src/equations.jl +++ b/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 diff --git a/src/fields.jl b/src/fields.jl index 58b53f7..c0449dd 100644 --- a/src/fields.jl +++ b/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} diff --git a/src/sparse.jl b/src/sparse.jl index 311304f..325fea9 100644 --- a/src/sparse.jl +++ b/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 diff --git a/src/types.jl b/src/types.jl index 52b47fa..cac08e1 100644 --- a/src/types.jl +++ b/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} diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl index 5c4427d..6b9936d 100644 --- a/test/test_elasticity.jl +++ b/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