Skip to content

Commit

Permalink
Merge pull request #3 from jtveiten/master
Browse files Browse the repository at this point in the history
* Better test coverage, use same template for running tests.
* Setting nodal boundary conditions / loads implemented.
  • Loading branch information
ahojukka5 committed Feb 15, 2018
2 parents ce7fffb + 61ce067 commit c337df5
Show file tree
Hide file tree
Showing 4 changed files with 267 additions and 60 deletions.
83 changes: 67 additions & 16 deletions src/FEMTruss.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,15 @@ import FEMBase: get_unknown_field_name,

"""
This truss fromulation is from
Concepts and applications of finite element analysis, forth edition
Chapter 2.4
Cook, Malkus, Plesha, Witt
# Features
- Nodal forces can be set using element `Poi1` with field
`nodal force i`, where `i` is dof number.
- Displacements can be fixed using element `Poi1` with field
`fixed displacement i`, where `i` is dof number.
"""
type Truss <: FieldProblem
end
Expand All @@ -21,23 +30,20 @@ function get_unknown_field_name(::Problem{Truss})
end

function get_formulation_type(::Problem{Truss})
return :incremental
return :total
end

"""
assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time)
Start finite element assembly procedure for Elasticity problem.
Function groups elements to arrays by their type and assembles one element type
at time. This makes it possible to pre-allocate matrices common to same type
of elements.
get_K_and_T(element::Element{Seg2}, nnodes, ndim, time)
This function assembles the local stiffness and transformation matrix
Need this to find element forces later on
Will need to change allocation strategy to pre-allocation later
Can discuss if we really need nnodes, since that is always 2 for trusses
"""
function assemble!(assembly::Assembly, problem::Problem{Truss},
element::Element{Seg2}, time)
#Require that the number of nodes = 2 ?
nnodes = length(element)
ndim = get_unknown_field_dimension(problem)
function get_K_and_T(element::Element{Seg2}, nnodes, ndim, time)
ndofs = nnodes*ndim
K = zeros(ndofs,ndofs)
T = zeros(2, ndofs*2)
node_id1 = element.connectivity[1] #First node in element
node_id2 = element.connectivity[2] #second node in elements
pos = element("geometry")
Expand All @@ -60,27 +66,72 @@ function assemble!(assembly::Assembly, problem::Problem{Truss},
E = loc_elem("youngs modulus", ip, time)
K_loc += ip.weight*E*A*dN'*dN*detJ
end
# Should we keep the local transform

if ndim == 1
K = K_loc
elseif ndim == 2
l_theta = dl[1]/l
m_theta = dl[2]/l
L = [l_theta m_theta 0 0; 0 0 l_theta m_theta]
K = L'*K_loc*L
T = [l_theta m_theta 0 0; 0 0 l_theta m_theta]
K = T'*K_loc*T
else # All three dimensions
l_theta = dl[1]/l
m_theta = dl[2]/l
n_theta = dl[3]/l
L = [l_theta m_theta n_theta 0 0 0; 0 0 0 l_theta m_theta n_theta]
K = L'*K_loc*L
T = [l_theta m_theta n_theta 0 0 0; 0 0 0 l_theta m_theta n_theta]
K = T'*K_loc*T
end
return (K,T)
end

"""
assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time)
Start finite element assembly procedure for Elasticity problem.
Function groups elements to arrays by their type and assembles one element type
at time. This makes it possible to pre-allocate matrices common to same type
of elements.
"""
function assemble!(assembly::Assembly, problem::Problem{Truss},
element::Element{Seg2}, time)
#Require that the number of nodes = 2 ?
nnodes = length(element)
ndim = get_unknown_field_dimension(problem)
K,T = get_K_and_T(element, nnodes, ndim, time)
gdofs = get_gdofs(problem, element)
add!(assembly.K, gdofs, gdofs, K)
#add!(assembly.f, gdofs, f)
end

import FEMBase: assemble_elements!

function assemble_elements!(problem::Problem, assembly::Assembly,
elements::Vector{Element{Poi1}}, time::Float64)
dim = ndofs = get_unknown_field_dimension(problem)
Ce = zeros(ndofs,ndofs)
ge = zeros(ndofs)
fe = zeros(ndofs)
for element in elements
fill!(Ce, 0.0)
fill!(ge, 0.0)
fill!(fe, 0.0)
ip = (0.0,)
for i=1:dim
if haskey(element, "fixed displacement $i")
Ce[i,i] = 1.0
ge[i] = element("fixed displacement $i", ip, time)
end
if haskey(element, "nodal force $i")
fe[i] = element("nodal force $i", ip, time)
end
end
gdofs = get_gdofs(problem, element)
add!(assembly.C1, gdofs, gdofs, Ce)
add!(assembly.C2, gdofs, gdofs, Ce)
add!(assembly.g, gdofs, ge)
add!(assembly.f, gdofs, fe)
end
end

export Truss

end
59 changes: 15 additions & 44 deletions test/runtests.jl
Original file line number Diff line number Diff line change
@@ -1,51 +1,24 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using FEMBase
using FEMTruss

using Base.Test
using TimerOutputs
const to = TimerOutput()

@testset "FEMTruss.jl" begin

elem1 = Element(Seg2, [1, 2]) # connects to nodes 1, 2
elem1.id = 1
X = Dict{Int64, Vector{Float64}}(
1 => [0.0],
2 => [10.0])

update!(elem1, "geometry", X)
update!(elem1, "youngs modulus", 288.0)
update!(elem1, "cross section area", 0.1)
p1 = Problem(Truss, "my truss problem", 1) # 1 dofs/node for now;
empty!(p1.assembly)
assemble!(p1.assembly, p1, elem1, 0.0)
K_truss = full(p1.assembly.K)
@test isapprox(K_truss, [2.88 -2.88;-2.88 2.88])

# now for the 2d case
X = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0],
2 => [10.0, 0.0])

update!(elem1, "geometry", X)
p1 = Problem(Truss, "my truss problem", 2) # 1 dofs/node for now;
empty!(p1.assembly)
assemble!(p1.assembly, p1, elem1, 0.0)
K_truss = full(p1.assembly.K)
@test isapprox(K_truss, [2.88 0.0 -2.88 0.0; 0.0 0.0 0.0 0.0; -2.88 0.0 2.88 0.0; 0.0 0.0 0.0 0.0])
test_files = String[]
push!(test_files, "test_elements.jl")
push!(test_files, "test_problems.jl")

# now for the 3d case
X = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0, 0.0],
2 => [10.0, 0.0, 0.0])

update!(elem1, "geometry", X)
p1 = Problem(Truss, "my truss problem", 3) # 1 dofs/node for now;
empty!(p1.assembly)
assemble!(p1.assembly, p1, elem1, 0.0)
K_truss = full(p1.assembly.K)
@test isapprox(K_truss, [2.88 0.0 0.0 -2.88 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; -2.88 0.0 0.0 2.88 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0])
@testset "FEMTruss.jl" begin
for fn in test_files
timeit(to, fn) do
include(fn)
end
end
end
println()
println("Test statistics:")
println(to)

# Need some boundary conditions and forces
# simple truss along x axis
Expand All @@ -56,5 +29,3 @@ K_truss = full(p1.assembly.K)
#solver = Solver(Linear, p1, f1, b1)
#solver()
# u = 10*10/(288*0.1)

end
113 changes: 113 additions & 0 deletions test/test_elements.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using FEMBase
using FEMTruss
using FEMBase: get_formulation_type

using Base.Test
# thsi is the same for all tests testing K and orientation
function make_test_problem(coords::Vector{Tuple{Int64, Vector{Float64}}}, dim::Int64, t::Float64=0.0)
elem_id = map(x->x[1], coords)
coord_dict=Dict(coords)
elem1 = Element(Seg2, elem_id) # connects to nodes 1, 2
elem1.id = 1
update!(elem1, "geometry", coord_dict)
update!(elem1, "youngs modulus", 288.0)
update!(elem1, "cross section area", 0.1)
p1 = Problem(Truss, "my truss problem", dim) # 1 dofs/node for now;
empty!(p1.assembly)
assemble!(p1.assembly, p1, elem1, t)
return p1
end

@testset "FEMTrussElements.jl" begin

zeros_2d = zeros(1,2)
zeros_3d = zeros(1,3)
K_oracle_base = [2.88 -2.88;-2.88 2.88]
#truss test
X = [(1 , [0.0]),(2, [10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)

@test get_unknown_field_name(p1)=="displacement"

#This one fails
#Test threw an exception of type UndefVarError
# Expression: get_formulation_type(p1)
# UndefVarError: get_formulation_type not defined
@test get_formulation_type(p1) == :total

K_truss = full(p1.assembly.K)
K_oracle = K_oracle_base
@test isapprox(K_truss, K_oracle)

# now for the 2d x case
X = [(1 , [0.0,0.0]),(2, [10.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [1 0]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# lets do y as well
X = [(1 , [0.0,0.0]),(2, [0.0,10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [0 1]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# now for the 2d x case reverted
X = [(2 , [0.0,0.0]),(1, [10.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [-1 0]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# lets do y as well reverted
X = [(2 , [0.0,0.0]),(1, [0.0,10.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_2d = [0 -1]
trans = [trans_2d zeros_2d;zeros_2d trans_2d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

# now for the 3d case
X = [(1 , [0.0,0.0,0.0]),(2, [10.0,0.0,0.0])]
ndofs = length(X[1][2])
p1 = make_test_problem(X, ndofs)
K_truss = full(p1.assembly.K)
trans_3d = [1 0 0]
trans = [trans_3d zeros_3d;zeros_3d trans_3d]
K_oracle = trans'*K_oracle_base*trans
@test isapprox(K_truss, K_oracle)

end

@testset "test nodal elements" begin
element = Element(Poi1, [1])
update!(element, "fixed displacement 1", 1.0)
update!(element, "nodal force 2", 2.0)
problem = Problem(Truss, "test problem", 2)
add_elements!(problem, [element])
assemble!(problem, 0.0)
C1 = full(problem.assembly.C1, 2, 2)
C2 = full(problem.assembly.C2, 2, 2)
f = full(problem.assembly.f, 2, 1)
g = full(problem.assembly.g, 2, 1)
@test isapprox(C1, [1.0 0.0; 0.0 0.0])
@test isapprox(C1, C2)
@test isapprox(f, [0.0, 2.0])
@test isapprox(g, [1.0, 0.0])
end
72 changes: 72 additions & 0 deletions test/test_problems.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,72 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMTruss.jl/blob/master/LICENSE

using FEMBase
using FEMBase.Test
using FEMTruss

@testset "FEMTrussProblems.jl" begin
E = 400.0*sqrt(5)
A = 0.1
Fx = 6.4*sqrt(2)
Fy = 1.6*sqrt(2)
X = Dict(1 => [0.0, 0.0],
2 => [0.0, 1.0],
3 => [1.0, 0.5])
el1 = Element(Seg2, [1, 3])
el2 = Element(Seg2, [2, 3])
bel1 = Element(Poi1, [1])
bel2 = Element(Poi1, [2])
bel3 = Element(Poi1, [3])
elements = [el1, el2, bel1, bel2, bel3]
update!(elements, "geometry", X)

update!([bel1, bel2], "fixed displacement 1", 0.0)
update!([bel1, bel2], "fixed displacement 2", 0.0)
update!(bel3, "nodal force 1", Fx)
update!(bel3, "nodal force 2", Fy)
update!(elements, "youngs modulus", E)
update!(elements, "cross section area", A)

problem = Problem(Truss, "test problem", 2)
add_elements!(problem, elements)

# solution
step = Analysis(Static)
add_problems!(step, [problem])
ls, normu, normla = run!(step)
# Lets verify the system
#Stiffness Matrix assemble by hand
l = sqrt(0.5^2+1^2)
K_loc = E*A/l*[1 -1; -1 1]
l1 = 1/l
m1 = 0.5/l
T1 = [l1 m1 0 0; 0 0 l1 m1]
K1 = T1'*K_loc*T1
T2 = [l1 -m1 0 0; 0 0 l1 -m1]
K2 = T2'*K_loc*T2
K_glob = zeros(6,6)
K_glob[3:6,3:6]+=K2
K_glob[1:2, 1:2] += K1[1:2,1:2]
K_glob[5:6, 1:2] += K1[3:4, 1:2]
K_glob[1:2, 5:6] += K1[1:2, 3:4]
K_glob[5:6, 5:6] += K1[3:4, 3:4]

@test isapprox(full(ls.K), K_glob)

# Forces are ok
f_glob = zeros(6);
f_glob[5:6]=[Fx,Fy]
@test isapprox(full(ls.f), f_glob)

# deflections
u_glob = zeros(6)
u_glob[5:6] = K_glob[5:6,5:6]\f_glob[5:6]
@test isapprox(full(ls.u), u_glob)

# support forces
#println("K = ", full(ls.K))
#println("f = ", full(ls.f))
#println("u = ", full(ls.u))
println("la = ", full(ls.la))
end

0 comments on commit c337df5

Please sign in to comment.