diff --git a/src/elasticity.jl b/src/elasticity.jl
index 17ca040..440ae8d 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -482,6 +482,13 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el
                 f[i:dim:end] += w*vec(T*N)
             end
         end
+        if haskey(element, "displacement traction force n")
+            J = element(ip, time, Val{:Jacobian})'
+            n = cross(J[:,1], J[:,2])
+            n /= norm(n)
+            p = element("displacement traction force n", ip, time)
+            f += w*p*vec(n*N)
+        end
     end
     return Kt, f
 end
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
new file mode 100644
index 0000000..46a7ed4
--- /dev/null
+++ b/src/postprocess_utils.jl
@@ -0,0 +1,33 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+"""
+Calculate field values to nodal points from Gauss points using least-squares fitting.
+"""
+function calc_nodal_values!(elements, field_name, field_dim, time)
+    A = SparseMatrixCOO()
+    b = SparseMatrixCOO()
+    for element in elements
+        gdofs = get_connectivity(element)
+        for ip in get_integration_points(element)
+            detJ = element(ip, time, Val{:detJ})
+            w = ip.weight*detJ
+            f = ip(field_name, time)
+            N = element(ip, time)
+            add!(A, gdofs, gdofs, w*kron(N', N))
+            for dim=1:field_dim
+                add!(b, gdofs, w*f[dim]*N, dim)
+            end
+        end
+    end
+    A = sparse(A)
+    b = sparse(b)
+    nz = get_nonzero_rows(A)
+    x = zeros(size(b)...)
+    x[nz, :] = A[nz,nz] \ b[nz, :]
+    nodal_values = Dict()
+    for i=1:size(x,1)
+        nodal_values[i] = vec(x[i,:])
+    end
+    update!(elements, field_name, nodal_values)
+end
diff --git a/test/test_elasticity_tet4_volume_load.jl b/test/test_elasticity_tet4_volume_load.jl
deleted file mode 100644
index cd6bc9a..0000000
--- a/test/test_elasticity_tet4_volume_load.jl
+++ /dev/null
@@ -1,34 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using JuliaFEM.Preprocess
-using JuliaFEM.Test
-
-@testset "test tet4 element under volume load" begin
-    X1 = [2.0, 3.0, 4.0]
-    X2 = [6.0, 3.0, 2.0]
-    X3 = [2.0, 5.0, 1.0]
-    X4 = [4.0, 3.0, 6.0]
-    e1 = Element(Tet4, [1, 2, 3, 4])
-    e1["geometry"] = Node[X1, X2, X3, X4]
-    e1["youngs modulus"] = 96.0
-    e1["poissons ratio"] = 1/3
-    e1["displacement load 3"] = 784.0/110.0
-    e2 = Element(Tri3, [1, 2, 3])
-    e2["geometry"] = Node[X2, X1, X3]
-    e2["displacement 1"] = 0.0
-    e2["displacement 2"] = 0.0
-    e2["displacement 3"] = 0.0
-    p1 = Problem(Elasticity, "tetra", 3)
-    p2 = Problem(Dirichlet, "bc", 3, "displacement")
-    p1.properties.finite_strain = false
-    push!(p1, e1)
-    push!(p2, e2)
-    s = Solver()
-    push!(s, p1, p2)
-    call(s)
-    u_4 = p1.assembly.u[10:end]
-    u_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
-    @test isapprox(u_4, u_expected)
-end
diff --git a/test/test_elasticity_tetra.jl b/test/test_elasticity_tetra.jl
new file mode 100644
index 0000000..c4b92da
--- /dev/null
+++ b/test/test_elasticity_tetra.jl
@@ -0,0 +1,66 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Test
+
+@testset "test tet4 + volume load" begin
+    X1 = [2.0, 3.0, 4.0]
+    X2 = [6.0, 3.0, 2.0]
+    X3 = [2.0, 5.0, 1.0]
+    X4 = [4.0, 3.0, 6.0]
+    e1 = Element(Tet4, [1, 2, 3, 4])
+    e1["geometry"] = Node[X1, X2, X3, X4]
+    e1["youngs modulus"] = 96.0
+    e1["poissons ratio"] = 1/3
+    e1["displacement load 3"] = 784.0/110.0
+    e2 = Element(Tri3, [1, 2, 3])
+    e2["geometry"] = Node[X2, X1, X3]
+    e2["displacement 1"] = 0.0
+    e2["displacement 2"] = 0.0
+    e2["displacement 3"] = 0.0
+    p1 = Problem(Elasticity, "tetra", 3)
+    p2 = Problem(Dirichlet, "bc", 3, "displacement")
+    push!(p1, e1)
+    push!(p2, e2)
+    s = Solver()
+    push!(s, p1, p2)
+    call(s)
+    u_4 = p1.assembly.u[10:end]
+    u_expected = [-3.0/220.0, -9.0/220.0, 1.0/10.0]
+    @test isapprox(u_4, u_expected)
+end
+
+@testset "test tet4 + surface load" begin
+    nodes = Dict{Int, Vector{Float64}}(
+        1 => [2.0, 3.0, 4.0],
+        2 => [6.0, 3.0, 2.0],
+        3 => [2.0, 5.0, 1.0],
+        4 => [4.0, 3.0, 6.0])
+    e1 = Element(Tet4, [1, 2, 3, 4])
+    update!(e1, "geometry", nodes)
+    e1["youngs modulus"] = 96.0
+    e1["poissons ratio"] = 1/3
+    e2 = Element(Tri3, [1, 2, 3])
+    update!(e2, "geometry", nodes)
+    e2["displacement 1"] = 0.0
+    e2["displacement 2"] = 0.0
+    e2["displacement 3"] = 0.0
+    e3 = Element(Tri3, [4, 3, 2])
+    update!(e3, "geometry", nodes)
+    e3["displacement traction force n"] = 96.0
+    p1 = Problem(Elasticity, "tetra", 3)
+    p2 = Problem(Dirichlet, "bc", 3, "displacement")
+    push!(p1, e1, e3)
+    push!(p2, e2)
+    s = Solver()
+    push!(s, p1, p2)
+    call(s)
+    u_4 = p1.assembly.u[10:end]
+    u_expected = [-17.0/14.0, -27.0/14.0, 1.0]
+    info("u_4 = $(u_4)")
+    info("u_expected = $(u_expected)")
+    @test isapprox(u_4, u_expected)
+end
+