diff --git a/src/aster_reader.jl b/src/aster_reader.jl
index 7a15fe2..7121ab7 100644
--- a/src/aster_reader.jl
+++ b/src/aster_reader.jl
@@ -2,6 +2,34 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 
+function aster_parse_nodes(section::ASCIIString; strip_characters=true)
+    nodes = Dict{Any, Vector{Float64}}()
+    has_started = false
+    for line in split(section, '\n')
+        m = matchall(r"[\w.-]+", line)
+        if (length(m) != 1) && (!has_started)
+            continue
+        end
+        if length(m) == 1
+            if (m[1] == "COOR_2D") || (m[1] == "COOR_3D")
+                has_started = true
+                continue
+            end
+            if m[1] == "FINSF"
+                break
+            end
+        end
+        if length(m) == 4
+            nid = m[1]
+            if strip_characters
+                nid = matchall(r"\d", nid)
+                nid = parse(Int, nid[1])
+            end
+            nodes[nid] = float(m[2:end])
+        end
+    end
+    return nodes
+end
 
 function parse(mesh::ASCIIString, ::Type{Val{:CODE_ASTER_MAIL}})
     model = Dict{ASCIIString, Any}()
diff --git a/src/core.jl b/src/core.jl
index 766bbc0..ece2689 100644
--- a/src/core.jl
+++ b/src/core.jl
@@ -80,6 +80,7 @@ include("equations.jl")
 include("dirichlet.jl")
 include("heat.jl")
 include("elasticity.jl")
+include("linear_elasticity.jl")
 
 ### ASSEMBLY + SOLVE ###
 include("assembly.jl")
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 05ddffe..83db89e 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -79,10 +79,10 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
     if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
         u = element("displacement", time, variation)
         grad = element(ip, time, Val{:grad})
-#        gradu = element("displacement", ip, time, Val{:grad}, variation)
         gradu = grad*u
 
         F = I + gradu # deformation gradient
+#       info("gradu = \n$(ForwardDiff.get_value(gradu))")
 
         young = element("youngs modulus", ip, time)
         poisson = element("poissons ratio", ip, time)
diff --git a/src/integrate.jl b/src/integrate.jl
index 3098855..8d6ba29 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -6,20 +6,22 @@
 
 ### 1d elements
 
-function get_integration_points(::Type{Seg2}, ::Type{Val{1}})
+typealias LineElement Union{Type{Seg2}, Type{Seg3}}
+
+function get_integration_points(::LineElement, ::Type{Val{1}})
     [
         IntegrationPoint([0.0], 2.0)
     ]
 end
 
-function get_integration_points(::Type{Seg2}, ::Type{Val{2}})
+function get_integration_points(::LineElement, ::Type{Val{2}})
     [
         IntegrationPoint([-sqrt(1/3)], 1)
         IntegrationPoint([+sqrt(1/3)], 1)
     ]
 end
 
-function get_integration_points(::Type{Seg3}, ::Type{Val{3}})
+function get_integration_points(::LineElement, ::Type{Val{3}})
     [
         IntegrationPoint([0.0], 8/9),
         IntegrationPoint([-sqrt(3/5)], 5/9),
@@ -27,7 +29,7 @@ function get_integration_points(::Type{Seg3}, ::Type{Val{3}})
     ]
 end
 
-function get_integration_points(::Type{Seg3}, ::Type{Val{4}})
+function get_integration_points(::LineElement, ::Type{Val{4}})
     [
         IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
         IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
@@ -36,7 +38,7 @@ function get_integration_points(::Type{Seg3}, ::Type{Val{4}})
     ]
 end
 
-function get_integration_points(::Union{Type{Seg2}, Type{Seg3}}, ::Type{Val{5}})
+function get_integration_points(::LineElement, ::Type{Val{5}})
     [
         IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
         IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
@@ -58,16 +60,16 @@ end
 
 # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
 
-typealias TriangularElements Union{Type{Tri3}, Type{Tri6}}
+typealias TriangularElement Union{Type{Tri3}, Type{Tri6}}
 
-function get_integration_points(::TriangularElements, ::Type{Val{1}})
+function get_integration_points(::TriangularElement, ::Type{Val{1}})
     # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
     [
         IntegrationPoint([1.0/3.0, 1.0/3.0], 0.5)
     ]
 end
 
-function get_integration_points(::TriangularElements, ::Type{Val{2}})
+function get_integration_points(::TriangularElement, ::Type{Val{2}})
     # http://libmesh.github.io/doxygen/quadrature__gauss__2D_8C_source.html
     [
         IntegrationPoint([2.0/3.0, 1.0/6.0], 1.0/6.0),
@@ -76,7 +78,7 @@ function get_integration_points(::TriangularElements, ::Type{Val{2}})
     ]
 end
 
-function get_integration_points(::TriangularElements, ::Type{Val{5}})
+function get_integration_points(::TriangularElement, ::Type{Val{5}})
     # http://math2.uncc.edu/~shaodeng/TEACHING/math5172/Lectures/Lect_15.PDF
     # FIXME: something wrong here with weights ..?
     [
@@ -95,7 +97,7 @@ function get_integration_points(::Type{Tri3})
 end
 
 
-function get_integration_points(::Type{Quad4})
+function get_integration_points(::Type{Quad4}, ::Type{Val{2}})
     [
         IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
         IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
@@ -104,8 +106,23 @@ function get_integration_points(::Type{Quad4})
     ]
 end
 
+function get_integration_points(::Type{Quad4})
+    return get_integration_points(Quad4, Val{2})
+end
+
 ### 3d elements
 
+
+function get_integration_points(::Type{Hex8}, ::Type{Val{2}})
+    p = 1.0/sqrt(3.0)*[-1.0, 1.0]
+    w = 1.0
+    return vec([IntegrationPoint([p[i], p[j], p[k]], w) for i=1:2, j=1:2, k=1:2])
+end
+
+function get_integration_points(::Type{Hex8})
+    return get_integration_points(Hex8, Val{2})
+end
+
 function get_integration_points(::Type{Tet4})
     # http://libmesh.github.io/doxygen/quadrature__gauss__3D_8C_source.html
     [
diff --git a/src/lagrange.jl b/src/lagrange.jl
index 02dd9e7..98b278a 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -113,6 +113,20 @@ end
 
 # 3d Lagrange elements
 
+@create_lagrange_element(Hex8, "8 node hexahedra",
+    [-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]])
+
 @create_lagrange_element(Tet4, "4 node tetrahedron",
     [0.0 1.0 0.0 0.0
      0.0 0.0 1.0 0.0
diff --git a/src/linear_elasticity.jl b/src/linear_elasticity.jl
new file mode 100644
index 0000000..40d2dc4
--- /dev/null
+++ b/src/linear_elasticity.jl
@@ -0,0 +1,105 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Linear elasticity
+
+abstract LinearElasticityProblem <: ElasticityProblem
+
+function LinearElasticityProblem(dim::Int=3, elements=[])
+    return Problem{LinearElasticityProblem}(dim, elements)
+end
+
+""" Elasticity equations, general 3D case. """
+function assemble!{E<:CG, P<:LinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real)
+
+    gdofs = get_gdofs(element, problem.dim)
+    ndim, nnodes = size(E)
+    B = zeros(6, 3*nnodes)
+    for ip in get_integration_points(element)
+        w = ip.weight*det(element, ip, time)
+        N = element(ip, time)
+        if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+            v = element("poissons ratio", ip, time)
+            E_ = element("youngs modulus", ip, time)
+            a = 1 - v
+            b = 1 - 2*v
+            c = 1 + v
+            C = E_/(b*c) .* [
+                a v v 0 0 0
+                v a v 0 0 0
+                v v a 0 0 0
+                0 0 0 b 0 0
+                0 0 0 0 b 0
+                0 0 0 0 0 b]
+            dN = element(ip, time, Val{:grad})
+            fill!(B, 0.0)
+            for i=1:size(dN, 2)
+                B[1, 3*(i-1)+1] = dN[1,i]
+                B[2, 3*(i-1)+2] = dN[2,i]
+                B[3, 3*(i-1)+3] = dN[3,i]
+                B[4, 3*(i-1)+1] = dN[2,i]
+                B[4, 3*(i-1)+2] = dN[1,i]
+                B[5, 3*(i-1)+2] = dN[3,i]
+                B[5, 3*(i-1)+3] = dN[2,i]
+                B[6, 3*(i-1)+1] = dN[3,i]
+                B[6, 3*(i-1)+3] = dN[1,i]
+            end
+            add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B)
+        end
+        if haskey(element, "displacement load")
+            b = element("displacement load", ip, time)
+            add!(assembly.force_vector, gdofs, w*N'*b)
+        end
+        if haskey(element, "displacement traction force")
+            T = element("displacement traction force", ip, time)
+            L = w*T*N
+#           dump(L)
+            add!(assembly.force_vector, gdofs, vec(L))
+        end
+    end
+end
+
+abstract PlaneStressLinearElasticityProblem <: LinearElasticityProblem
+
+function PlaneStressLinearElasticityProblem(dim::Int=2, elements=[])
+    return Problem{PlaneStressLinearElasticityProblem}(dim, elements)
+end
+
+""" Elasticity equations, plane stress. """
+function assemble!{E<:CG, P<:PlaneStressLinearElasticityProblem}(assembly::Assembly, problem::Problem{P}, element::Element{E}, time::Real)
+
+    gdofs = get_gdofs(element, problem.dim)
+    ndim, nnodes = size(E)
+    B = zeros(3, 2*nnodes)
+    for ip in get_integration_points(element)
+        w = ip.weight*det(element, ip, time)
+        N = element(ip, time)
+        if haskey(element, "youngs modulus") && haskey(element, "poissons ratio")
+            nu = element("poissons ratio", ip, time)
+            E_ = element("youngs modulus", ip, time)
+            C = E_/(1.0 - nu^2) .* [
+                1.0  nu 0.0
+                nu  1.0 0.0
+                0.0 0.0 (1.0-nu)/2.0]
+            dN = element(ip, time, Val{:grad})
+            fill!(B, 0.0)
+            for i=1:size(dN, 2)
+                B[1, 2*(i-1)+1] = dN[1,i]
+                B[2, 2*(i-1)+2] = dN[2,i]
+                B[3, 2*(i-1)+1] = dN[2,i]
+                B[3, 2*(i-1)+2] = dN[1,i]
+            end
+            add!(assembly.stiffness_matrix, gdofs, gdofs, w*B'*C*B)
+        end
+        if haskey(element, "displacement load")
+            b = element("displacement load", ip, time)
+            add!(assembly.force_vector, gdofs, w*N'*b)
+        end
+        if haskey(element, "displacement traction force")
+            T = element("displacement traction force", ip, time)
+            L = w*T*N
+#           dump(L)
+            add!(assembly.force_vector, gdofs, vec(L))
+        end
+    end
+end
diff --git a/test/test_aster_reader.jl b/test/test_aster_reader.jl
index ec2830f..a1a0dd8 100644
--- a/test/test_aster_reader.jl
+++ b/test/test_aster_reader.jl
@@ -6,7 +6,8 @@ module AsterReaderTests
 using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM: parse
+#using JuliaFEM: parse
+using JuliaFEM.Preprocess: aster_parse_nodes
 
 function test_read_mesh()
 mesh = """
@@ -44,7 +45,30 @@ mesh = """
     @test m["nsets"]["NALL"] == ["N1", "N2"]
 
 end
+#test_read_mesh()
 
-test_read_mesh()
+function test_parse_nodes()
+    section = """
+    N9  2.0 3.0 4.0
+    COOR_3D
+    N1          0.0 0.0 0.0
+    N2          1.0 0.0 0.0
+    N3          1.0 1.0 0.0
+    N4          0.0 1.0 0.0
+    N5          0.0 0.0 1.0
+    N6          1.0 0.0 1.0
+    N7          1.0 1.0 1.0
+    N8          0.0 1.0 1.0
+    FINSF
+    absdflasdf
+    N12 3.0 4.0 5.0 6.0
+    N13 3.0 4.0 5.0
+    """
+    nodes = aster_parse_nodes(section)
+    @test nodes[1] == Float64[0.0, 0.0, 0.0]
+    @test nodes[8] == Float64[0.0, 1.0, 1.0]
+    @test length(nodes) == 8
+end
+test_parse_nodes()
 
 end
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index 6b9936d..36c9268 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -4,7 +4,11 @@
 module ElasticityTests
 
 using JuliaFEM.Test
-using JuliaFEM.Core: Seg2, Quad4, PlaneStressElasticityProblem, solve!
+using JuliaFEM
+using JuliaFEM.Core: Seg2, Quad4, Hex8,
+                     ElasticityProblem, PlaneStressElasticityProblem,
+                     solve!, get_connectivity, DirichletProblem
+
 
 function test_elasticity_volume_load()
     element = Quad4([1, 2, 3, 4])
@@ -19,9 +23,19 @@ function test_elasticity_volume_load()
     free_dofs = [3, 4, 5, 6]
     solve!(problem, free_dofs, 0.0; max_iterations=10)
     disp = element("displacement", [1.0, 1.0], 0.0)
+    # function get_previous_ip(element::Element, current_ip::IntegrationPoint)
+    # end
+    # ipdata = element("integration points", time) => IntegrationPoint[ip1, ip2, ..., ipN]
+    # for some_ip in ipdata
+    #   if isapprox(some_ip.xi, ip.xi)
+    #       info("found")
+    #       last_value = some_ip("material parameter", time)
+    #       break
+    #   end
+    # end
     #ip1 = last(element["integration points"])[1]
     #ip2 = last(element["integration points"])[2]
-    #strain = ip1["gl strain"]
+    # strain = ip1("gl strain")
     info("displacement at tip: $disp")
     #info("strain in first ip: $strain. ip coord = $(ip1.xi) and weight = $(ip1.weight)")
     # verified using Code Aster, verification/2015-10-22-plane-stress/cplan_grot_gdep_volume_force.resu
@@ -29,6 +43,7 @@ function test_elasticity_volume_load()
 end
 #test_elasticity_volume_load()
 
+
 function test_elasticity_surface_load()
     N = Vector[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]
 
@@ -55,4 +70,93 @@ function test_elasticity_surface_load()
     @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
 end
 
+
+function test_continuum_elasticity_with_surface_load()
+    nodes = JuliaFEM.Preprocess.aster_parse_nodes("""
+    COOR_3D
+    N1          0.0 0.0 0.0
+    N2          1.0 0.0 0.0
+    N3          1.0 1.0 0.0
+    N4          0.0 1.0 0.0
+    N5          0.0 0.0 1.0
+    N6          1.0 0.0 1.0
+    N7          1.0 1.0 1.0
+    N8          0.0 1.0 1.0
+    FINSF
+    """)
+
+    function set_geometry!(element, nodes)
+        element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
+    end
+    element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
+    set_geometry!(element1, nodes)
+#    element1["youngs modulus"] = 900.0
+#    element1["poissons ratio"] = 0.25
+    element1["youngs modulus"] = 9000.0
+    element1["poissons ratio"] = 0.25
+    element1["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:8])
+
+    element2 = Quad4([5, 6, 7, 8])
+    set_geometry!(element2, nodes)
+    element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]
+    element2["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0, 0.0] for i=1:4])
+
+    problem = ElasticityProblem()
+    push!(problem, element1)
+    push!(problem, element2)
+
+    #=
+    free_dofs = zeros(Bool, 8, 3)
+    x = 1
+    y = 2
+    z = 3
+    free_dofs[2, x] = true
+    free_dofs[3, [x, y]] = true
+    free_dofs[4, y] = true
+    free_dofs[5, z] = true
+    free_dofs[6, [x, z]] = true
+    free_dofs[7, [x, y, z]] = true
+    free_dofs[8, [y, z]] = true
+    free_dofs = find(vec(free_dofs'))
+    info("free dofs: $free_dofs")
+
+    info("initial force vector")
+    ass = JuliaFEM.Core.assemble(problem, 0.0)
+    info(reshape(full(ass.force_vector), 3, 8))
+    info("initial stiffness matrix")
+    dump(round(Int, full(ass.stiffness_matrix))[free_dofs, free_dofs])
+    solve!(problem, free_dofs, 0.0; max_iterations=10)
+    =#
+
+    dx = Quad4([1, 4, 8, 5])
+    dx["displacement 1"] = 0.0
+    dy = Quad4([1, 5, 6, 2])
+    dy["displacement 2"] = 0.0
+    dz = Quad4([1, 2, 3, 4])
+    dz["displacement 3"] = 0.0
+    bc = DirichletProblem("displacement", 3)
+    for el in [dx, dy, dz]
+        set_geometry!(el, nodes)
+        push!(bc, el)
+    end
+    solver = JuliaFEM.Core.DirectSolver()
+    push!(solver, problem)
+    push!(solver, bc)
+    solver.dump_matrices = true
+    solver.name = "3d_hex8"
+    solver(0.0)
+
+    disp = element1("displacement", [1.0, 1.0, 1.0], 0.0)
+    info("displacement at tip: $disp")
+    info("displacement on element: ")
+    for (i, d) in enumerate(element1("displacement", 0.0))
+        @printf "%d %f %f %f\n" [i;d]...
+    end
+    # verified using Code Aster.
+    # 2015-12-12-continuum-elasticity/vim c3d_grot_gdep_traction_force.comm
+    # @test isapprox(disp, [3.17431158889468E-02, 3.17431158889468E-02, -1.38591518927826E-01])
+    @test isapprox(disp, [2.80559539222183E-03, 2.80559539222183E-03, -1.13019918093242E-02])
+end
+#test_continuum_elasticity_with_surface_load()
+
 end
diff --git a/test/test_linear_elasticity.jl b/test/test_linear_elasticity.jl
new file mode 100644
index 0000000..b273fe7
--- /dev/null
+++ b/test/test_linear_elasticity.jl
@@ -0,0 +1,123 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module LinearElasticityTests
+
+using JuliaFEM
+using JuliaFEM.Test
+
+using JuliaFEM.Core: Seg2, Quad4, Hex8, LinearElasticityProblem, get_connectivity,
+                     assemble, PlaneStressLinearElasticityProblem
+using JuliaFEM.Preprocess: aster_parse_nodes
+
+
+function test_plane_stress_linear_elasticity_with_surface_load()
+    nodes = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+
+    function set_geometry!(element, nodes)
+        element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
+    end
+    element1 = Quad4([1, 2, 3, 4])
+    set_geometry!(element1, nodes)
+    element1["youngs modulus"] = 9000.0
+    element1["poissons ratio"] = 0.25
+
+    element2 = Seg2([3, 4])
+    set_geometry!(element2, nodes)
+    element2["displacement traction force"] = Vector{Float64}[[0.0, -100.0] for i=1:2]
+
+    problem = PlaneStressLinearElasticityProblem()
+    push!(problem, element1)
+    push!(problem, element2)
+
+    free_dofs = Int64[3, 5, 6, 8]
+
+    info("initial force vector")
+    ass = assemble(problem, 0.0)
+    f = full(ass.force_vector)
+    K = full(ass.stiffness_matrix)
+    dump(reshape(f, 2, 4))
+    info("initial stiffness matrix")
+    dump(round(Int, K)[free_dofs, free_dofs])
+
+    u = zeros(2, 4)
+    u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
+
+    info("result vector")
+    dump(u)
+    # verified using Code Aster.
+    # 2015-10-22-plane-stress/cplan_linear_traction_force.*
+    @test isapprox(u[:,3], [2.77777777777778E-03, -1.11111111111111E-02])
+end
+#test_plane_stress_linear_elasticity_with_surface_load()
+
+
+function test_continuum_elasticity_with_surface_load()
+    nodes = JuliaFEM.Preprocess.aster_parse_nodes("""
+    COOR_3D
+    N1          0.0 0.0 0.0
+    N2          1.0 0.0 0.0
+    N3          1.0 1.0 0.0
+    N4          0.0 1.0 0.0
+    N5          0.0 0.0 1.0
+    N6          1.0 0.0 1.0
+    N7          1.0 1.0 1.0
+    N8          0.0 1.0 1.0
+    FINSF
+    """)
+
+    function set_geometry!(element, nodes)
+        element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
+    end
+    element1 = Hex8([1, 2, 3, 4, 5, 6, 7, 8])
+    set_geometry!(element1, nodes)
+    element1["youngs modulus"] = 9000.0
+    element1["poissons ratio"] = 0.25
+
+    element2 = Quad4([5, 6, 7, 8])
+    set_geometry!(element2, nodes)
+    element2["displacement traction force"] = Vector{Float64}[[0.0, 0.0, -100.0] for i=1:4]
+
+    problem = LinearElasticityProblem()
+    push!(problem, element1)
+    push!(problem, element2)
+
+    free_dofs = zeros(Bool, 8, 3)
+    x = 1
+    y = 2
+    z = 3
+    free_dofs[2, x] = true
+    free_dofs[3, [x, y]] = true
+    free_dofs[4, y] = true
+    free_dofs[5, z] = true
+    free_dofs[6, [x, z]] = true
+    free_dofs[7, [x, y, z]] = true
+    free_dofs[8, [y, z]] = true
+    free_dofs = find(vec(free_dofs'))
+    info("free dofs: $free_dofs")
+
+    info("initial force vector")
+    ass = assemble(problem, 0.0)
+    f = full(ass.force_vector)
+    K = full(ass.stiffness_matrix)
+    dump(reshape(f, 3, 8))
+    info("initial stiffness matrix")
+    dump(round(Int, K)[free_dofs, free_dofs])
+
+    u = zeros(3, 8)
+    u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
+
+    info("result vector")
+    dump(u)
+    # verified using Code Aster.
+    # 2015-12-12-continuum-elasticity/c3d_linear.*
+    @test isapprox(u[:,7], [2.77777777777778E-02, 2.77777777777778E-02, -1.11111111111111E-01])
+end
+#test_continuum_elasticity_with_surface_load()
+
+
+end