Added checks for test, missing element and support force checks

src/FEMTruss.jl

@@ -12,6 +12,9 @@ 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
@@ -31,19 +34,16 @@ function get_formulation_type(::Problem{Truss})
 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")
@@ -66,22 +66,37 @@ 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)

test/test_problems.jl

@@ -35,8 +35,38 @@ using FEMTruss
     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))
+    # 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