Merge 98bd054 into ce7fffb

src/FEMTruss.jl

@@ -12,6 +12,12 @@ import FEMBase: get_unknown_field_name,
 
 """
 This truss fromulation is from
+
+# 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
@@ -21,7 +27,7 @@ function get_unknown_field_name(::Problem{Truss})
 end
 
 function get_formulation_type(::Problem{Truss})
-    return :incremental
+    return :total
 end
 
 """
@@ -81,6 +87,36 @@ function assemble!(assembly::Assembly, problem::Problem{Truss},
     #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

test/runtests.jl

@@ -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
@@ -56,5 +29,3 @@ K_truss = full(p1.assembly.K)
   #solver = Solver(Linear, p1, f1, b1)
   #solver()
   # u = 10*10/(288*0.1)
-
-end

test/test_elements.jl

@@ -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

test/test_problems.jl

@@ -0,0 +1,42 @@
+# 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)
+    println("K = ", full(ls.K))
+    println("f = ", full(ls.f))
+    println("u = ", full(ls.u))
+    println("la = ", full(ls.la))
+end