From 6b32d364e5bc8d66d5591eae76cb3af14ab99509 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Mon, 30 Nov 2015 18:49:45 +0200 Subject: [PATCH] updated JuliaFEM.jl --- src/JuliaFEM.jl | 11 +++++++- src/elements.jl | 5 ++++ src/lagrange.jl | 35 +++++++++++++++++------- test/test_directsolver.jl | 56 ++++++++++++++++++++++++++++++++++++--- 4 files changed, 93 insertions(+), 14 deletions(-) diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 28bcf70..0a27ee5 100644 --- a/src/JuliaFEM.jl +++ b/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 diff --git a/src/elements.jl b/src/elements.jl index 55436b0..e9693e2 100644 --- a/src/elements.jl +++ b/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 diff --git a/src/lagrange.jl b/src/lagrange.jl index dbdbf01..fcac227 100644 --- a/src/lagrange.jl +++ b/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 diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl index 6ec6958..6245b1e 100644 --- a/test/test_directsolver.jl +++ b/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