diff --git a/src/elasticity.jl b/src/elasticity.jl
index be9feab..17ca040 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -184,17 +184,16 @@ function assemble{El<:Union{Seg2,Seg3}}(problem::Problem{Elasticity}, element::E
     return Kt, f
 end
 
-
-""" Elasticity equations, continuum formulation. """
-function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+""" Elasticity equations, 3d, linear. """
+function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
-    nnodes = size(element, 2)
-    BL = zeros(6, dim*nnodes)
-    BNL = zeros(9, dim*nnodes)
-    Kt = zeros(dim*nnodes, dim*nnodes)
-    f = zeros(dim*nnodes)
+    nnodes = length(element)
+    ndofs = dim*nnodes
+    BL = zeros(6, ndofs)
+    Kt = zeros(ndofs, ndofs)
+    f = zeros(ndofs)
 
     for ip in get_integration_points(element)
         detJ = element(ip, time, Val{:detJ})
@@ -202,17 +201,17 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
         N = element(ip, time)
         dN = element(ip, time, Val{:Grad})
 
-        # kinematics; calculate deformation gradient and strain
-        gradu = zeros(dim, dim)
-        if haskey(element, "displacement")
-            gradu += element("displacement", ip, time, Val{:Grad})
-        end
-        strain = zeros(dim , dim)
-        strain += 1/2*(gradu' + gradu)
-        F = eye(dim)
-        if props.finite_strain
-            F += gradu
-            strain += 1/2*gradu'*gradu
+        fill!(BL, 0.0)
+        for i=1:nnodes
+            BL[1, 3*(i-1)+1] = dN[1,i]
+            BL[2, 3*(i-1)+2] = dN[2,i]
+            BL[3, 3*(i-1)+3] = dN[3,i]
+            BL[4, 3*(i-1)+1] = dN[2,i] + dN[1,i]
+            BL[4, 3*(i-1)+2] = dN[2,i] + dN[1,i]
+            BL[5, 3*(i-1)+2] = dN[3,i] + dN[2,i]
+            BL[5, 3*(i-1)+3] = dN[3,i] + dN[2,i]
+            BL[6, 3*(i-1)+1] = dN[1,i] + dN[3,i]
+            BL[6, 3*(i-1)+3] = dN[1,i] + dN[3,i]
         end
 
         E = element("youngs modulus", ip, time)
@@ -226,35 +225,142 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
             0.0 0.0 0.0 0.0 0.5-nu 0.0
             0.0 0.0 0.0 0.0 0.0 0.5-nu]
 
-        # calculate stress
+        Kt += w*BL'*D*BL
+
+        if haskey(element, "displacement load")
+            T = element("displacement load", ip, time)
+            f += w*vec(T*N)
+        end
+        for i=1:dim
+            if haskey(element, "displacement load $i")
+                b = element("displacement load $i", ip, time)
+                f[i:dim:end] += w*vec(b*N)
+            end
+        end
+
+    end
+
+    if get_formulation_type(problem) == :incremental
+        if haskey(element, "displacement")
+            u = vec(element["displacement"](time))
+            f -= Kt*u
+        end
+    end
+
+    return Kt, f
+end
+
+""" Material and geometric stiffness for linear buckling analysis. """
+function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_buckling}})
+
+    props = problem.properties
+    dim = get_unknown_field_dimension(problem)
+    nnodes = length(element)
+    ndofs = dim*nnodes
+    BL = zeros(6, ndofs)
+    BNL = zeros(9, ndofs)
+    Km = zeros(ndofs, ndofs)
+    Kg = zeros(ndofs, ndofs)
+
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        dN = element(ip, time, Val{:Grad})
+
+        gradu = element("displacement", ip, time, Val{:Grad})
+        strain = 1/2*(gradu' + gradu)
+
+        fill!(BL, 0.0)
+        for i=1:nnodes
+            BL[1, 3*(i-1)+1] = dN[1,i]
+            BL[2, 3*(i-1)+2] = dN[2,i]
+            BL[3, 3*(i-1)+3] = dN[3,i]
+            BL[4, 3*(i-1)+1] = dN[2,i] + dN[1,i]
+            BL[4, 3*(i-1)+2] = dN[2,i] + dN[1,i]
+            BL[5, 3*(i-1)+2] = dN[3,i] + dN[2,i]
+            BL[5, 3*(i-1)+3] = dN[3,i] + dN[2,i]
+            BL[6, 3*(i-1)+1] = dN[1,i] + dN[3,i]
+            BL[6, 3*(i-1)+3] = dN[1,i] + dN[3,i]
+        end
+
+        fill!(BNL, 0.0)
+        for i=1:size(dN, 2)
+            BNL[1, 3*(i-1)+1] = dN[1,i]
+            BNL[2, 3*(i-1)+1] = dN[2,i]
+            BNL[3, 3*(i-1)+1] = dN[3,i]
+            BNL[4, 3*(i-1)+2] = dN[1,i]
+            BNL[5, 3*(i-1)+2] = dN[2,i]
+            BNL[6, 3*(i-1)+2] = dN[3,i]
+            BNL[7, 3*(i-1)+3] = dN[1,i]
+            BNL[8, 3*(i-1)+3] = dN[2,i]
+            BNL[9, 3*(i-1)+3] = dN[3,i]
+        end
+
+        E = element("youngs modulus", ip, time)
+        nu = element("poissons ratio", ip, time)
+        D = E/((1.0+nu)*(1.0-2.0*nu)) * [
+            1.0-nu nu nu 0.0 0.0 0.0
+            nu 1.0-nu nu 0.0 0.0 0.0
+            nu nu 1.0-nu 0.0 0.0 0.0
+            0.0 0.0 0.0 0.5-nu 0.0 0.0
+            0.0 0.0 0.0 0.0 0.5-nu 0.0
+            0.0 0.0 0.0 0.0 0.0 0.5-nu]
+
         strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
         stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
 
-        stress = [
-            stress_vec[1] stress_vec[4] stress_vec[6]
-            stress_vec[4] stress_vec[2] stress_vec[5]
-            stress_vec[6] stress_vec[5] stress_vec[3]]
-        cauchy_stress = F'*stress*F/det(F)
-        cauchy_stress_vec = [
-            cauchy_stress[1,1];
-            cauchy_stress[2,2];
-            cauchy_stress[3,3];
-            cauchy_stress[1,2];
-            cauchy_stress[2,3];
-            cauchy_stress[1,3]]
+        S3 = zeros(3*dim, 3*dim)
+        S3[1,1] = stress_vec[1]
+        S3[2,2] = stress_vec[2]
+        S3[3,3] = stress_vec[3]
+        S3[2,3] = S3[3,2] = stress_vec[4]
+        S3[1,3] = S3[3,1] = stress_vec[5]
+        S3[1,2] = S3[2,1] = stress_vec[6]
+        S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
 
-        s = cauchy_stress - 1.0/3.0*trace(cauchy_stress)*I
-        J2 = 1/2*trace(s*s')
+        Km += w*BL'*D*BL
+        Kg += w*BNL'*S3*BNL
 
-        # update values to integration point
-        update!(ip, "strain", time => strain_vec)
-        update!(ip, "stress", time => stress_vec)
-        update!(ip, "cauchy stress", time => cauchy_stress_vec)
-        update!(ip, "von mises stress", time => sqrt(3.0*J2))
+    end
+
+    return Km, Kg
+end
+
+""" Elasticity equations, 3d nonlinear. """
+function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+
+    props = problem.properties
+    dim = get_unknown_field_dimension(problem)
+    nnodes = length(element)
+    ndofs = dim*nnodes
+    BL = zeros(6, ndofs)
+    BNL = zeros(9, ndofs)
+    Kt = zeros(ndofs, ndofs)
+    f = zeros(ndofs)
+
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        dN = element(ip, time, Val{:Grad})
+
+        # kinematics; calculate deformation gradient and strain
+
+        gradu = zeros(dim, dim)
+        if haskey(element, "displacement")
+            gradu += element("displacement", ip, time, Val{:Grad})
+        end
+        strain = zeros(dim , dim)
+        strain += 1/2*(gradu' + gradu)
+        F = eye(dim)
+        if props.finite_strain
+            F += gradu
+            strain += 1/2*gradu'*gradu
+        end
 
-        # add contributions: material and geometric stiffness + internal forces
         fill!(BL, 0.0)
-        for i=1:size(dN, 2)
+        for i=1:nnodes
             BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
             BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
             BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
@@ -275,6 +381,37 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
             BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
         end
 
+        # material stiffness start
+
+        E = element("youngs modulus", ip, time)
+        nu = element("poissons ratio", ip, time)
+
+        D = E/((1.0+nu)*(1.0-2.0*nu)) * [
+            1.0-nu nu nu 0.0 0.0 0.0
+            nu 1.0-nu nu 0.0 0.0 0.0
+            nu nu 1.0-nu 0.0 0.0 0.0
+            0.0 0.0 0.0 0.5-nu 0.0 0.0
+            0.0 0.0 0.0 0.0 0.5-nu 0.0
+            0.0 0.0 0.0 0.0 0.0 0.5-nu]
+
+        # calculate stress
+        strain_vec = [strain[1,1]; strain[2,2]; strain[3,3]; strain[1,2]; strain[2,3]; strain[1,3]]
+        stress_vec = D * ([1.0, 1.0, 1.0, 2.0, 2.0, 2.0].*strain_vec)
+
+        # update values to integration point
+        update!(ip, "strain", time => strain_vec)
+        update!(ip, "stress", time => stress_vec)
+
+        Kt += w*BL'*D*BL
+
+        if get_formulation_type(problem) == :incremental
+            f -= w*BL'*stress_vec
+        end
+
+        # material stiffness end
+
+        # geometric stiffness start
+
         fill!(BNL, 0.0)
         for i=1:size(dN, 2)
             BNL[1, 3*(i-1)+1] = dN[1,i]
@@ -296,16 +433,14 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
         S3[1,2] = S3[2,1] = stress_vec[6]
         S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]
 
-        Kt += w*BL'*D*BL
         if props.finite_strain
             Kt += w*BNL'*S3*BNL
         end
 
-        if get_formulation_type(problem) == :incremental
-            f -= w*BL'*stress_vec
-        end
+        # geometric stiffness end
+
+        # external load start
 
-        # volume load
         if haskey(element, "displacement load")
             T = element("displacement load", ip, time)
             f += w*vec(T*N)
@@ -316,15 +451,10 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
                 f[i:dim:end] += w*vec(b*N)
             end
         end
-    end
 
-#=
-    if get_formulation_type(problem) == :incremental
-        if haskey(element, "displacement")
-            f -= Kt*vec(element["displacement"](time))
-        end
+        # external load end
+
     end
-=#
 
     return Kt, f
 end
@@ -356,6 +486,10 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el
     return Kt, f
 end
 
+function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
+    return assemble(problem, element, time, Val{:continuum})
+end
+
 """ Elasticity equations using ForwardDiff
 
 Formulation
diff --git a/src/lagrange_macro.jl b/src/lagrange_macro.jl
index da4b4be..043f0fb 100644
--- a/src/lagrange_macro.jl
+++ b/src/lagrange_macro.jl
@@ -105,15 +105,7 @@ end
     [-1.0  1.0  1.0 -1.0 -1.0  1.0 1.0 -1.0
      -1.0 -1.0  1.0  1.0 -1.0 -1.0 1.0  1.0
      -1.0 -1.0 -1.0 -1.0  1.0  1.0 1.0  1.0],
-    (xi) -> [
-    1.0,
-    xi[1],
-    xi[2],
-    xi[1]*xi[2],
-    xi[3],
-    xi[1]*xi[3],
-    xi[2]*xi[3],
-    xi[1]*xi[2]*xi[3]])
+    (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2], xi[3], xi[1]*xi[3], xi[2]*xi[3], xi[1]*xi[2]*xi[3]])
 
 @create_lagrange_element(Tet4, "4 node tetrahedron",
     [0.0 1.0 0.0 0.0
@@ -121,14 +113,6 @@ end
      0.0 0.0 0.0 1.0],
     (xi) -> [1.0, xi[1], xi[2], xi[3]])
 
-#=
-@create_lagrange_element(Tet4, "4 node tetrahedron",
-    [0.0 0.0 0.0 1.0
-     1.0 0.0 0.0 0.0
-     0.0 1.0 0.0 0.0],
-    (xi) -> [1.0, xi[1], xi[2], xi[3]])
-=#
-
 @create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
     [0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
      0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5
@@ -136,24 +120,6 @@ end
     (xi) -> [    1.0,   xi[1],       xi[2],       xi[3],     xi[1]^2,
              xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
 
-#=
-@create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
-    [0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.5 0.5
-     1.0 0.0 0.0 0.0 0.5 0.0 0.5 0.5 0.0 0.0
-     0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0],
-    (xi) -> [    1.0,   xi[1],       xi[2],       xi[3],     xi[1]^2,
-             xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
-=#
-
-# some helpers to make accessing 1d basis functions more easy
-function get_basis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T)
-    get_basis(E, [xi])
-end
-
-function get_dbasis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T)
-    get_dbasis(E, [xi])
-end
-
 function get_reference_element_midpoint{E}(element::Element{E})
     get_reference_element_midpoint(E)
 end
diff --git a/src/problems.jl b/src/problems.jl
index 6e7270f..8ce244f 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -8,14 +8,16 @@ abstract MixedProblem <: AbstractProblem
 
 """
 General linearized problem to solve
-    K*u +   C1.T*la  = f
-    C2*u +  D*la    = g
+    (K₁+K₂)*Δu +   C1.T*λ = f₁+f₂
+         C2*Δu +      D*λ = g
 """
 type Assembly
     # for field assembly
     M :: SparseMatrixCOO  # mass matrix
     K :: SparseMatrixCOO  # stiffness matrix
+    Kg :: SparseMatrixCOO  # geometric stiffness matrix
     f :: SparseMatrixCOO  # force vector
+#   f2 :: SparseMatrixCOO
     # for boundary assembly
     C1 :: SparseMatrixCOO
     C2 :: SparseMatrixCOO
@@ -44,14 +46,16 @@ function Assembly()
         SparseMatrixCOO(),
         SparseMatrixCOO(),
         SparseMatrixCOO(),
+        SparseMatrixCOO(),
         [], [], Inf,
         [], [], Inf,
         true)
 end
 
-function Base.empty!(assembly::Assembly)
+function empty!(assembly::Assembly)
     empty!(assembly.M)
     empty!(assembly.K)
+    empty!(assembly.Kg)
     empty!(assembly.f)
     empty!(assembly.C1)
     empty!(assembly.C2)
@@ -287,7 +291,7 @@ element.element connectivity using formula gdofs = [dim*(nid-1)+j for j=1:dim]
 2. if not found, use element.connectivity to update dofmap and 1.
 """
 function get_gdofs(problem::Problem, element::Element)
-    if !haskey(element, problem.dofmap)
+    if !haskey(problem.dofmap, element)
         dim = get_unknown_field_dimension(problem)
         problem.dofmap[element] = get_gdofs(element, dim)
     end
diff --git a/test/test_define_new_element.jl b/test/test_define_new_element.jl
index 39bde9d..69bc592 100644
--- a/test/test_define_new_element.jl
+++ b/test/test_define_new_element.jl
@@ -9,13 +9,13 @@ import JuliaFEM: get_basis, get_dbasis, get_integration_points
 type MyQuad4 <: AbstractElement
 end
 
-function get_basis(element::Element{MyQuad4}, xi, time)
-    1/4*[(1-xi[1])*(1-xi[2]) (1+xi[1])*(1-xi[2]) (1+xi[1])*(1+xi[2]) (1-xi[1])*(1+xi[2])]
+function get_basis(element::Element{MyQuad4}, ip, time)
+    1/4*[(1-ip[1])*(1-ip[2]) (1+ip[1])*(1-ip[2]) (1+ip[1])*(1+ip[2]) (1-ip[1])*(1+ip[2])]
 end
 
-function get_dbasis(element::Element{MyQuad4}, xi, time)
-    1/4*[-(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
-         -(1-xi[1])   -(1+xi[1])   (1+xi[1])   (1-xi[1])]
+function get_dbasis(element::Element{MyQuad4}, ip, time)
+    1/4*[-(1-ip[2])    (1-ip[2])   (1+ip[2])  -(1+ip[2])
+         -(1-ip[1])   -(1+ip[1])   (1+ip[1])   (1-ip[1])]
 end
 
 function get_integration_points(element::MyQuad4)
@@ -39,7 +39,7 @@ end
     el = Element(MyQuad4)
     el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
     el["displacement"] = Vector{Float64}[[0.0,0.0], [0.0,0.0], [1.0,0.0], [0.0,0.0]]
-    @test isapprox(el("geometry", [0.0, 0.0]), [0.5, 0.5])
+    @test isapprox(el("geometry", [0.0, 0.0], 0.0), [0.5, 0.5])
     @test isapprox(el("displacement", [0.0, 0.0], 0.0), [0.25, 0.0])
     el["temperature thermal conductivity"] = 6.0
     dim = length(el)