diff --git a/src/abaqus_reader_old.jl b/src/abaqus_reader_old.jl
index fccd7b6..4638b93 100644
--- a/src/abaqus_reader_old.jl
+++ b/src/abaqus_reader_old.jl
@@ -44,7 +44,8 @@ function parse_element_section(model, header, data)
     eldims = Dict(
         "C3D10" => 10,
         "C3D4" => 4,
-        "S3" => 3)
+        "S3" => 3,
+        "STRI65" => 6)
     eltype = header["options"]["TYPE"]
     if !(eltype in keys(eldims))
         throw("Element $eltype dimension information missing")
diff --git a/src/assembly.jl b/src/assembly.jl
index 8222759..703eaf5 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -3,6 +3,16 @@
 
 # Functions to handle global assembly of problem
 
+type CAssembly
+    interior_dofs :: Vector{Int}
+    boundary_dofs :: Vector{Int}
+    F :: Factorization
+    Kc :: SparseMatrixCSC
+    fc :: SparseMatrixCSC
+    Ki :: SparseMatrixCSC
+    fi :: SparseMatrixCSC
+end
+
 function assemble!(assembly::Assembly, problem::AllProblems, time::Float64, empty_assembly::Bool=true)
     if empty_assembly
         empty!(assembly)
@@ -20,6 +30,73 @@ function assemble(problem::AllProblems, time::Float64)
     return assembly
 end
 
+""" Return condensed system. """
+function assemble(problem::FieldProblem, time::Float64, boundary_dofs::Vector{Int})
+    assembly = Assembly()
+    for element in get_elements(problem)
+        assemble!(assembly, problem, element, time)
+    end
+    return condensate(assembly, boundary_dofs)
+end
+
+function condensate(assembly::Assembly, boundary_dofs_::Vector{Int})
+    K = sparse(assembly.stiffness_matrix)
+    all_dofs = unique(assembly.stiffness_matrix.I)
+    boundary_dofs = intersect(all_dofs, boundary_dofs_)
+    interior_dofs = setdiff(all_dofs, boundary_dofs_)
+    dim = size(K, 1)
+    f = sparse(assembly.force_vector, dim, 1)
+
+    # check that matrix is symmetric
+    asdf = maximum(abs(1/2*(K + K') - K))
+    if asdf > 1.0e-6
+        info(full(K))
+        error("asdf $asdf > 1.0e-6")
+    end
+    
+    K = 1/2*(K + K')
+    F::Factorization = cholfact(K[interior_dofs, interior_dofs])
+
+#   info("condensation: all dofs: ", all_dofs)
+#   info("condensation: interior dofs: ", interior_dofs)
+#   info("condensation: boundary dofs: ", boundary_dofs)
+#   info("manually condensated")
+#   Kman = K[boundary_dofs, boundary_dofs] - K[boundary_dofs,interior_dofs] * inv(full(K[interior_dofs, interior_dofs])) * K[interior_dofs, boundary_dofs]
+#   info("\n$(full(Kman))")
+    #info("K = \n$(full(K))")
+    #Ki = K[interior_dofs, boundary_dofs]
+    Ki = K[interior_dofs, boundary_dofs]
+    fi = f[interior_dofs]
+#   info("condensated using factorization")
+#   LL = K[boundary_dofs, boundary_dofs] - K[boundary_dofs, interior_dofs] * (K[interior_dofs, interior_dofs] \ K[interior_dofs, boundary_dofs])
+#   info(LL)
+
+    Ks = F \ Ki
+    Fs = F \ fi
+
+    dim = size(K, 1)
+    Kc = spzeros(dim, dim)
+    fc = spzeros(dim, 1)
+    Kc[boundary_dofs, boundary_dofs] = K[boundary_dofs, boundary_dofs] - Ki' * Ks
+    fc[boundary_dofs] = f[boundary_dofs] - Ki' * Fs
+
+
+    return CAssembly(interior_dofs, boundary_dofs, F, Kc, fc, Ki, fi)
+end
+
+function reconstruct!(ca::CAssembly, x::SparseMatrixCSC)
+#   info("size of la = ", size(la))
+#   info("size of ca.Ki = ", size(ca.Ki))
+#   info("size of ca.fi = ", size(ca.fi))
+#   info("size of la[ca.interior_dofs] = ", size(la[ca.interior_dofs]))
+#   info("interior dofs: $(ca.interior_dofs)")
+#   info("boundary dofs: $(ca.boundary_dofs)")
+#   info("ca.fi = $(ca.fi')")
+#   info("sol1 = ", full(ca.F \ ca.fi)')
+#   info("sol2 = ", full(ca.F \ (ca.Ki*x[ca.boundary_dofs]))')
+    x[ca.interior_dofs] += ca.F \ (ca.fi - ca.Ki*x[ca.boundary_dofs])
+end
+
 function Base.(:+)(ass1::Assembly, ass2::Assembly)
     mass_matrix = ass1.mass_matrix + ass2.mass_matrix
     stiffness_matrix = ass1.stiffness_matrix + ass2.stiffness_matrix
diff --git a/src/directsolver.jl b/src/directsolver.jl
index af2d9a6..10f0334 100644
--- a/src/directsolver.jl
+++ b/src/directsolver.jl
@@ -99,48 +99,93 @@ function call(solver::DirectSolver, time::Number=0.0)
 
         mapper = solver.parallel ? pmap : map
 
+        info("Assembling problems.")
         # assemble boundary problems
         tic(timing, "boundary assembly")
         boundary_assembly = sum(mapper((p)->assemble(p, time), solver.boundary_problems))
         boundary_dofs = unique(boundary_assembly.stiffness_matrix.I)
+#       boundary_dofs = collect(range(1, 12))
+        info("# of interface dofs: $(length(boundary_dofs))")
+        #info("dofs = $boundary_dofs")
         toc(timing, "boundary assembly")
 
-        # assemble field problems
-        # in principle if we want to static condensation we need to pass boundary dofs
-        # to field problems in order to know which dofs are interior dofs and can be
-        # condensated.
-        tic(timing, "field assembly")
-        field_assembly = sum(mapper((p)->assemble(p, time), solver.field_problems))
-        field_dofs = unique(field_assembly.stiffness_matrix.I)
-        info("# of dofs: $(length(field_dofs)), # of interface dofs: $(length(boundary_dofs))")
-        toc(timing, "field assembly")
+        #static_condensation = false
 
+        # assemble field problems
+        dim = 0
+        assemblies = []
+        for (i, problem) in enumerate(solver.field_problems)
+            tic(timing, "field assembly")
+            field_assembly = assemble(problem, time)
+            #info("full assembly body $i")
+            #info(round(full(field_assembly.stiffness_matrix), 3))
+            toc(timing, "field assembly")
+            field_dofs = unique(field_assembly.stiffness_matrix.I)
+            #info("# of dofs in problem $i: $(length(field_dofs))")
+            dim = maximum([dim, maximum(field_dofs)])
+            #tic(timing, "condensate")
+            #cfield_assembly = condensate(field_assembly, boundary_dofs)
+            #toc(timing, "condensate")
+            #push!(assemblies, cfield_assembly)
+            push!(assemblies, field_assembly)
+        end
+
+        #info("dim = $dim")
+
+        #info("assembly done")
         tic(timing, "create sparse matrices")
         # create sparse matrices and saddle point problem
-        K = sparse(field_assembly.stiffness_matrix)
-        dim = size(K, 1)
-        r = sparse(field_assembly.force_vector, dim, 1)
+        #K = sparse(field_assembly.stiffness_matrix)
+        #dim = size(K, 1)
+        #r = sparse(field_assembly.force_vector, dim, 1)
+
+        K = spzeros(dim, dim)
+        r = spzeros(dim, 1)
+
+        for (i, assembly) in enumerate(assemblies)
+            #info("body $i")
+            #info(round(full(assembly.Kc), 3))
+            #resize!(assembly.stiffness_matrix, dim, dim)
+            #resize!(assembly.force_vector, dim, 1)
+            K += sparse(assembly.stiffness_matrix, dim, dim)
+            r += sparse(assembly.force_vector, dim, 1)
+        end
+        #info(round(full(K), 3))
+
         C = sparse(boundary_assembly.stiffness_matrix, dim, dim)
         g = sparse(boundary_assembly.force_vector, dim, 1)
         A = [K C'; C spzeros(dim, dim)]
         b = [r; g]
         toc(timing, "create sparse matrices")
+        #info("problem size = ", size(A))
 
+        info("Solving system")
         tic(timing, "solution of system")
         # solve increment for linearized problem
         nz = unique(rowvals(A))  # take only non-zero rows
         sol = zeros(b)
-        sol[nz] = lufact(A[nz,nz]) \ full(b[nz])
-        info("solved. length of solution vector = $(length(sol))")
-        toc(timing, "solution of system")
-        #info(full(sol[nz]))
+        sol[nz] = A[nz,nz] \ full(b[nz])
+        #info("solution vector before reconstruction")
+        #info(full(sol)')
+        la = sol[dim+1:end]
+        #for assembly in assemblies
+        #    reconstruct!(assembly, sol)
+        #end
+        la = vec(full(la))
+        sol = vec(full(sol))
+        #info("la = ", la')
+        #info("sol = ", sol')
 
+        info("solved. solution norm: $(norm(sol[1:dim]))")
+        toc(timing, "solution of system")
+
+        info("Updating element data")
         tic(timing, "update element data")
         # update elements in field problems
         for field_problem in solver.field_problems
             for element in get_elements(field_problem)
                 gdofs = get_gdofs(element, field_dim)
-                local_sol = vec(full(sol[gdofs]))  # incremental data for element
+                local_sol = sol[gdofs]  # incremental data for element
                 local_sol = reshape(local_sol, field_dim, length(element))
                 local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
                 last(element[field_name]).data += local_sol  # <-- added
@@ -150,8 +195,8 @@ function call(solver::DirectSolver, time::Number=0.0)
         # update elements in boundary problems
         for boundary_problem in solver.boundary_problems
             for element in get_elements(boundary_problem)
-                gdofs = get_gdofs(element, field_dim) + dim
-                local_sol = vec(full(sol[gdofs]))
+                gdofs = get_gdofs(element, field_dim)
+                local_sol = la[gdofs]
                 local_sol = reshape(local_sol, field_dim, length(element))
                 local_sol = Vector{Float64}[local_sol[:,i] for i=1:length(element)]
                 last(element["reaction force"]).data = local_sol  # <-- replaced
@@ -164,6 +209,7 @@ function call(solver::DirectSolver, time::Number=0.0)
             info("timing info for non-linear iteration:")
             info("boundary assembly      : ", time_elapsed(timing, "boundary assembly"))
             info("field assembly         : ", time_elapsed(timing, "field assembly"))
+#           info("condensate             : ", time_elapsed(timing, "condensate"))
             info("create sparse matrices : ", time_elapsed(timing, "create sparse matrices"))
             info("solution of system     : ", time_elapsed(timing, "solution of system"))
             info("update element data    : ", time_elapsed(timing, "update element data"))
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 914cd8c..83a2591 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -76,6 +76,12 @@ function get_residual_vector{P<:ElasticityProblem}(problem::Problem{P}, element:
         end
         E = 1/2*(F'*F - I)  # strain
         S = lambda*trace(E)*I + 2*mu*E
+
+        J = det(element, ip, time)
+        T = J^-1*F*S*F'
+        #ip["cauchy stress"] = T
+        ip["gl strain"] = E
+
         r += F*S*dbasis
     end
 
diff --git a/src/elements.jl b/src/elements.jl
index 00c7b72..55436b0 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -5,7 +5,6 @@ abstract AbstractElement
 
 type Element{E<:AbstractElement}
     connectivity :: Vector{Int}
-#   integration_points :: Vector{IntegrationPoint}
     fields :: Dict{ASCIIString, Field}
 end
 
@@ -15,10 +14,23 @@ function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
 end
 
 function get_integration_points{E}(element::Element{E})
-#   return element.integration_points
     return get_integration_points(E)
 end
 
+function update_gauss_fields!(element::Element, data::Vector{IntegrationPoint}, time::Real)
+    if haskey(element, "integration points")
+        # push or update
+        if !isapprox(last(element["integration points"]).time, time)
+            push!(element["integration points"], time => data)
+        else
+            last(element["integration points"]).data = data
+        end
+    else
+        # create
+        element["integration points"] = Field(time => data)
+    end
+end
+
 """
 Test routine for element. If this passes, element interface is properly
 defined.
@@ -46,7 +58,7 @@ function test_element(element_type)
 
     info("Initializing element")
     try
-        element = element_type(collect(1:n))
+        element = Element{element_type}(collect(1:n))
     catch
         error("""
         Unable to create element with default constructor define function
@@ -55,22 +67,20 @@ function test_element(element_type)
     end
 
     # try to interpolate some scalar field
-    element["field1"] = Field(collect(1:n))
+    element["field1"] = range(1, n)
     # TODO: how to parametrize this?
-    element["geometry"] = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    element["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
 
     # evaluate basis functions at middle point of element
-    basis = get_basis(element)
-    dbasis = grad(basis)
     mid = zeros(dim)
-    val1 = basis(mid, 0.0)
+    val1 = element(mid, 0.0)
     info("basis at $mid: $val1")
-    val2 = basis("field1", mid, 0.0)
+    val2 = element("field1", mid, 0.0)
     info("field val at $mid: $val2")
-    val3 = dbasis(mid, 0.0)
+    val3 = element(mid, 0.0, Val{:grad})
     info("derivative of basis at $mid:\n$val3")
-    val4 = dbasis("field1", mid, 0.0)
-    info("field val at $mid: $val4")
+    #val4 = element("field1", mid, Val{:grad})
+    #info("field val at $mid: $val4")
 
     info("Element $element_type passed tests.")
 end
@@ -122,11 +132,16 @@ function call(element::Element, field_name::ASCIIString, xi::VecOrIP, time::Numb
 end
 
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP)
-    return element.basis(element[field_name], xi)
+    field = element[field_name]
+    basis = get_basis(element)
+    return basis(element[field_name], xi)
 end
 
 function call(element::Element, field_name::ASCIIString, xi::VecOrIP, ::Type{Val{:grad}})
-    return element.basis(element["geometry"], element[field_name], xi, Val{:grad})
+    field = element[field_name]
+    geom = element["geometry"]
+    basis = get_basis(element)
+    return basis(geom, field, xi, Val{:grad})
 end
 
 function call(element::Element, field_name::ASCIIString, time::Number)
@@ -137,6 +152,10 @@ function get_basis{E}(element::Element{E}, ip::IntegrationPoint)
     return get_basis(E, ip.xi)
 end
 
+function get_basis{E}(::Type{Element{E}}, xi::Vector{Float64})
+    return get_basis(E, xi)
+end
+
 function get_basis{E}(element::Element{E}, xi::Vector{Float64})
     return get_basis(E, xi)
 end
@@ -152,9 +171,15 @@ function get_basis{E}(element::Element{E})
     return basis
 end
 
+function call{E}(element::Element{E}, xi::VecOrIP, ::Type{Val{:grad}})
+    basis = get_basis(element)
+    geom = element["geometry"]
+    return basis(geom, xi, Val{:grad})
+end
+
 function call{E}(element::Element{E}, xi::VecOrIP, time::Float64, ::Type{Val{:grad}})
     basis = get_basis(element)
-    return basis(element["geometry"], xi, Val{:grad})
+    return basis(element["geometry"](time), xi, Val{:grad})
 end
 
 function call(element::Element, field_name::ASCIIString)
diff --git a/src/equations.jl b/src/equations.jl
index 4672bad..d816cd7 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -95,7 +95,7 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
     end
 
     # 2. energy form -- user has defined potential energy W -> min!
-    if has_potential_energy(problem, element)
+    if has_potential_energy(problem, element) && haskey(element, unknown_field_name)
         field = element[unknown_field_name](time)
 
         """ Wrapper for potential energy for ForwardDiff. """
@@ -122,7 +122,7 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
     end
 
     # 3. virtual work -- user has defined some residual r = p - f = 0
-    if has_residual_vector(problem, element)
+    if has_residual_vector(problem, element) && haskey(element, unknown_field_name)
 
         field = DVTI(last(element[unknown_field_name]).data)
 
@@ -130,17 +130,24 @@ function assemble!(assembly::Assembly, problem::Problem, element::Element, time:
         function calc_R(data::Vector)
             R = zeros(length(data))
             df = similar(field, data)
+            gauss_fields = IntegrationPoint[]
             # integrate residual vector
             for ip in get_integration_points(element)
                 s = ip.weight*det(element, ip, time)
                 dr = get_residual_vector(problem, element, ip, time; variation=df)
                 R += s*dr
+                if ip.changed
+                    push!(gauss_fields, ip)
+                end
             end
             # external loads -- if any nodal loads is defined, decrease from residual
             if haskey(element, "$unknown_field_name nodal load")
                 R -= vec(element["$unknown_field_name nodal load"](time))
             end
             #info("return = $R")
+            if length(gauss_fields) != 0
+                update_gauss_fields!(element, gauss_fields, time)
+            end
             return R
         end
 
diff --git a/src/problems.jl b/src/problems.jl
index 033e7bc..1021c3f 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -15,6 +15,8 @@ type BoundaryProblem{T<:AbstractProblem}
     elements :: Vector{Element}
 end
 
+typealias FieldProblem Problem
+
 typealias AllProblems Union{Problem, BoundaryProblem}
 
 function get_elements(problem::AllProblems)
diff --git a/src/types.jl b/src/types.jl
index 41e45d3..52b47fa 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -21,34 +21,48 @@ fields
 type IntegrationPoint
     xi :: Vector
     weight :: Float64
-    fields :: FieldSet
+    fields :: Dict{ASCIIString, Field}
+    changed :: Bool
 end
 
 function IntegrationPoint(xi, weight)
-    return IntegrationPoint(xi, weight, FieldSet())
+    return IntegrationPoint(xi, weight, FieldSet(), false)
 end
 
-function Base.convert(::Type{Number}, ip::IntegrationPoint)
+function setindex!{T<:ForwardDiff.ForwardDiffNumber}(ip::IntegrationPoint, data::Array{T,2}, field_name::ASCIIString)
+    data = ForwardDiff.get_value(data)
+    setindex!(ip, data, field_name)
+end
+function setindex!(ip::IntegrationPoint, data, field_name)
+    ip.fields[field_name] = Field(data)
+    ip.changed = true
+end
+
+function getindex(ip::IntegrationPoint, field_name::ASCIIString)
+    ip.fields[field_name]
+end
+
+function convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
 
-function Base.call(field::CVTI, ip::IntegrationPoint)
+function call(field::CVTI, ip::IntegrationPoint)
     return call(field, ip.xi)
 end
 
-function Base.call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
+function call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
     call(basis, field, ip.xi)
 end
 
-function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+function call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
     call(basis, field, ip.xi, Val{:grad})
 end
 
-function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
+function call(basis::CVTI, field::DVTI, ip::IntegrationPoint)
     call(basis, field, ip.xi)
 end
 
-function Base.call(basis::CVTI, geometry::DVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+function call(basis::CVTI, geometry::DVTI, field::Union{DCTI, DVTI}, ip::IntegrationPoint, ::Type{Val{:grad}})
     call(basis, geometry, field, ip.xi, Val{:grad})
 end
 
diff --git a/test/test_assembly.jl b/test/test_assembly.jl
new file mode 100644
index 0000000..297dccb
--- /dev/null
+++ b/test/test_assembly.jl
@@ -0,0 +1,57 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module AssemblyTests
+
+using JuliaFEM
+using JuliaFEM.Test
+using JuliaFEM: Seg2, Quad4, HeatProblem, DirichletProblem, assemble
+using JuliaFEM: condensate, reconstruct!
+
+function test_static_condensation()
+    nodes = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+
+    el1 = Quad4([1, 2, 3, 4])
+    el1["geometry"] = Vector[nodes[1], nodes[2], nodes[3], nodes[4]]
+    el1["temperature thermal conductivity"] = 6.0
+    el1["temperature load"] = [12.0, 12.0, 12.0, 12.0]
+    el2 = Seg2([1, 2])
+    el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
+    el2["temperature flux"] = 6.0
+    field_problem = HeatProblem()
+    push!(field_problem, el1)
+    push!(field_problem, el1)
+
+    el3 = Seg2([3, 4])
+    el3["geometry"] = Vector[nodes[3], nodes[4]]
+    el3["temperature"] = 0.0
+    boundary_problem = DirichletProblem("temperature", 1)
+    push!(boundary_problem, el3)
+
+    fass = assemble(field_problem, 0.0)
+    bass = assemble(boundary_problem, 0.0)
+
+#   interior_dofs = [1, 2]
+    boundary_dofs = [3, 4]
+    cass = condensate(fass, boundary_dofs)
+    @test isapprox(full(cass.Kc), [
+        0.0 0.0  0.0  0.0
+        0.0 0.0  0.0  0.0
+        0.0 0.0  4.8 -4.8
+        0.0 0.0 -4.8  4.8])
+   @test isapprox(full(cass.fc)', [0.0 0.0 12.0 12.0])
+   @test cass.interior_dofs == [1, 2]
+
+   x = sparse(zeros(4))'
+   la = sparse(zeros(4))'
+   la[3] = la[4] = 24.0
+   reconstruct!(cass, x)
+   x = full(x)
+   info(la)
+   info(x)
+   @test isapprox(x[1], 1.0)
+   @test isapprox(x[2], 1.0)
+end
+
+end
+
diff --git a/test/test_basis.jl b/test/test_basis.jl
index e0d745f..a5fa952 100644
--- a/test/test_basis.jl
+++ b/test/test_basis.jl
@@ -6,93 +6,94 @@ module BasisTests
 using JuliaFEM.Test
 
 using JuliaFEM
-using JuliaFEM: Basis, Field
-using JuliaFEM: Increment, TimeStep
+using JuliaFEM: AbstractElement, Element
 
-function get_basis()
+import JuliaFEM: get_basis, get_dbasis
 
-    basis(xi) = 1/4*[
+abstract TestElement <: AbstractElement
+
+function get_basis(::Type{TestElement}, xi::Vector{Float64})
+    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])]'
-    
-    dbasis(xi) = 1/4*[
+end
+
+function get_dbasis(::Type{TestElement}, xi::Vector{Float64})
+    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])]
+end
 
-    return Basis(basis, dbasis)
+function get_element()
+    element = Element{TestElement}([1, 2, 3, 4])
+    element["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    element["temperature"] = Float64[1.0, 2.0, 3.0, 4.0]
+    element["displacement1"] = Vector{Float64}[[0.0, 0.0], [0.0,  0.0], [1/4, 0.0], [0.0, 0.0]]
+    element["displacement2"] = Vector{Float64}[[0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0]]
+    return element
 end
 
 ### Test interpolation in spatial domain
 
 function test_basis_interpolation()
-    N = get_basis()
-    @test N([0.0, 0.0]) == 1/4*[1 1 1 1]
-    @test N([0.0, 0.0], 1.0) == 1/4*[1 1 1 1]
+    element = get_element()
+    info(element([0.0, 0.0]))
+    @test element([0.0, 0.0]) == 1/4*[1 1 1 1]
+    @test element([0.0, 0.0], 1.0) == 1/4*[1 1 1 1]
 end
 
 function test_basis_gradient_interpolation()
-    X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
-#   P(X) = [1.0, X[1], X[2], X[1]*X[2]]
-#   basis2, dbasis2 = JuliaFEM.calculate_lagrange_basis(P, X)
-    N = get_basis()
-    gradN = N(X, [0.0, 0.0], Val{:grad})
-    @test gradN == 1/2*[-1 1 1 -1; -1 -1 1 1]
-#   @test dN([0.0, 0.0]) == dbasis2([0.5, 0.5])
+    element = get_element()
+    grad = element([0.0, 0.0], Val{:grad})
+    info("grad = \n$grad")
+    @test grad == 1/2*[-1 1 1 -1; -1 -1 1 1]
 end
 
-function test_interpolation_of_scalar_increment_in_spatial_domain()
+function test_interpolation_of_scalar_field_in_spatial_domain()
     # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+    element = get_element()
     T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
-    T = Increment([1.0, 2.0, 3.0, 4.0])
-    N = get_basis()
-    T_interpolated = N(T, [0.0, 0.0])
+    T_interpolated = element("temperature", [0.0, 0.0])
     @test T_interpolated == T_known([0.5, 0.5])
 end
 
-function test_interpolation_of_gradient_of_scalar_increment_in_spatial_domain()
+function test_interpolation_of_gradient_of_scalar_field_in_spatial_domain()
     # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
-    X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
-    T = Increment([1.0, 2.0, 3.0, 4.0])
-    N = get_basis()
-    gradT = N(X, T, [0.0, 0.0], Val{:grad})
+    element = get_element()
+    gradT = element("temperature", [0.0, 0.0], Val{:grad})
     gradT_expected(X) = [1-2*X[2] 3-2*X[1]]
     @test gradT == gradT_expected([0.5, 0.5])
 end
 
 function test_interpolation_of_vector_field()
     # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
-    geometry = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
-    displacement = Increment(Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]])
-    N = get_basis()
-    X = N(geometry, [0.0, 0.0])
-    u = N(displacement, [0.0, 0.0])
-    x = X+u
+    element = get_element()
+    u = element("displacement1", [0.0, 0.0])
+#   x = X+u
     u_expected(X) = [1/4*X[1]*X[2], 0]
-    @test isapprox(x, [9/16, 1/2])
+#   @test isapprox(x, [9/16, 1/2])
     @test isapprox(u, u_expected([0.5, 0.5]))
 end
 
+
 function test_interpolation_of_gradient_of_vector_field()
     # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
     # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
-
+    element = get_element()
 #    displacement = Field(
 #        (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
 #        (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
-
-    u = Increment([0.0 0.0; 1.0 -1.0; 2.0 3.0; 0.0 0.0]')
-
-    N = get_basis()
-    gradu(xi) = N(X, u, xi, Val{:grad})
+    gradu = element("displacement2", [0.0, 0.0], Val{:grad})
     gradu_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    @test isapprox(gradu([0.0, 0.0]), gradu_expected([0.5, 0.5]))
+    @test isapprox(gradu, gradu_expected([0.5, 0.5]))
 end
 
 ### Test interpolation in time domain
 
+#=
+
 function test_linear_time_extrapolation_of_field()
     #T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
     T = Field(
@@ -188,8 +189,6 @@ function test_derivative_interpolation_in_temporal_basis_in_variable_velocity_ch
     @test isa(velocity, Increment) == true
 end
 
-#=
-
 function test_time_derivative_gradient_interpolation_of_field()
     # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
     # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
diff --git a/test/test_directsolver.jl b/test/test_directsolver.jl
index 5c6b160..6ec6958 100644
--- a/test/test_directsolver.jl
+++ b/test/test_directsolver.jl
@@ -62,6 +62,7 @@ function test_solver_multiple_dirichlet_bc()
     @test isapprox(disp, [3.17431158889468E-02, -1.38591518927826E-01])
 
 end
+#test_solver_multiple_dirichlet_bc()
 
 function test_solver_multiple_bodies_multiple_dirichlet_bc()
     N = Vector[
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index cad3d69..4a6410c 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -19,10 +19,15 @@ 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)
+    ip1 = last(element["integration points"])[1]
+    ip2 = last(element["integration points"])[2]
+    strain = ip1["gl strain"]
     info("displacement at tip: $disp")
-    # verified using Code Aster.
+    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
     @test isapprox(disp[2], -8.77303119819776)
 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]]
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 822a812..a3bbb5b 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -5,38 +5,37 @@ module ElementTests
 
 using JuliaFEM.Test
 
-using JuliaFEM: Element, Field, FieldSet, test_element
+using JuliaFEM: AbstractElement, Element, Field, FieldSet, test_element
+import JuliaFEM: get_basis, get_dbasis
+import Base: size
 
 """ Prototype element
 
 This should always pass test_element if everything is ok.
 """
-type MockElement <: Element
-    connectivity :: Vector{Int}
-    basis :: Field
-    fields :: FieldSet
-end
+abstract TestElement <: AbstractElement
 
-function MockElement(connectivity, fields...)
-
-    h(xi) = 1/4*[
+function get_basis(::Type{TestElement}, xi::Vector{Float64})
+    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])]'
-
-    dh(xi) = 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])]
-
-    MockElement(connectivity, Field(h, dh), FieldSet(fields...))
 end
 
-Base.size(element::Type{MockElement}) = (2, 4)
+function get_dbasis(::Type{TestElement}, xi::Vector{Float64})
+    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])]
+end
+
+function size(::Type{TestElement})
+    return (2, 4)
+end
 
 """ Return test element with some fields. """
 function get_element()
-    el = MockElement([1, 2, 3, 4])
+    el = Element{TestElement}([1, 2, 3, 4])
     el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
     el["temperature"] = (
         0.0 => [0.0, 0.0, 0.0, 0.0],
@@ -48,7 +47,7 @@ function get_element()
 end
 
 function test_mock_element()
-    test_element(MockElement)
+    test_element(TestElement)
 end
 
 function test_add_fields_to_element()
@@ -69,11 +68,11 @@ function test_interpolate()
     info("gradT(expected) = $gradT_expected")
     @test isapprox(gradT, gradT_expected)
 
-    @test isapprox(el("temperature", [0.0, 0.0], 0.5), 1/2*gradT_expected)
+#   @test isapprox(el("temperature", [0.0, 0.0], 0.5), 1/2*gradT_expected)
 
-    gradT = el("temperature", [0.0, 0.0], 0.5, Val{:grad})
-    info("gradT = $gradT")
-    @test isapprox(gradT, 1/2*gradT_expected)
+#   gradT = el("temperature", [0.0, 0.0], 0.5, Val{:grad})
+#   info("gradT = $gradT")
+#   @test isapprox(gradT, 1/2*gradT_expected)
 end
 
 end
diff --git a/test/test_fields.jl b/test/test_fields.jl
index df4c871..f9cca2b 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -5,378 +5,12 @@
 module FieldTests
 
 using JuliaFEM
-using JuliaFEM: Increment, TimeStep, Field, DefaultDiscreteField, FieldSet
-using JuliaFEM: ContinuousField, DiscreteField, DefaultContinuousField
-
 using JuliaFEM.Test
 
-function test_increment_constant_increment()
-    I = Increment(1)
-    @test isa(I, Increment)
-    @test length(I) == 1
-    @test I == 1
-end
+using JuliaFEM: Field
 
-function test_increments_with_vector_data()
-    I1 = Increment([1, 2, 3])
-    I2 = Increment([2, 3, 4])
-    @test length(I1) == 3
-    @test length(I2) == 3
-    @test I1 == [1, 2, 3]
-    @test I2 == [2, 3, 4]
-end
-
-function test_increments_basic_math()
-    I1 = Increment([1, 2, 3])
-    I2 = Increment([2, 3, 4])
-    @test 1/2*(I1+I2) == [1.5, 2.5, 3.5]
-    @test I1 + 1 == [2, 3, 4]
-    @test I1 - 1 == [0, 1, 2]
-    @test I1*3 == [3, 6, 9]
-    @test I1+I2 == [3, 5, 7]
-end
-
-function test_increment_dot_product()
-    I1 = Increment([1, 2, 3])
-    I2 = Increment([2, 3, 4])
-    @test dot(I1, I2) == 20
-    @test dot([1,2,3], I2) == 20
-    @test dot(I1, [2,3,4]) == 20
-    @test dot([1, 2], Increment[I1, I2])
-end
-
-function test_increment_similarity()
-    f = zeros(Increment, Int, 2, 4)
-    @test length(f) == 4
-    g = similar(f, ones(Int, 8))
-    @test typeof(f) == typeof(g)
-    @test length(f) == length(g)
-    @test size(g) == (2, 4)
-end
-
-function test_increment_vec()
-    g = zeros(Increment, Int, 2, 4)
-    @test vec(g) == ones(Int, 8)
-end
-
-function test_increment_promotion()
-    I1 = Increment([1, 2, 3])
-    I2 = Increment([2, 3, 4])
-    @test isa(I1+1, Increment)
-    @test isa(I1-1, Increment)
-    @test isa(3*I1, Increment)
-    @test isa(1/2*I1, Increment)
-    @test isa(I1+I2, Increment)
-    @test isa(I1-I2, Increment)
-end
-
-function test_timestep_empty_timestep()
-    ts = TimeStep()
-    @test length(ts) == 0
-    @test ts.time == 0.0
-end
-
-function test_timestep_with_two_increments()
-    i1 = Increment([1, 2, 3])
-    i2 = Increment([2, 3, 4])
-    increments = Increment[i1, i2]
-    ts = TimeStep(1.0, increments)
-    @test length(ts) == 2
-end
-
-function test_create_timestep_with_scalar_value()
-    ts = TimeStep(1)
-    @test length(ts) == 1
-    @test ts.time == 0.0
-    @test isa(ts[1], Increment)
-    @test ts[1] == [1]
-end
-
-function test_create_timestep_compactly_for_time_t0()
-    ts = TimeStep([1, 2, 3])
-    @test length(ts) == 1
-    @test ts.time == 0.0
-    @test isa(ts[1], Increment)
-    @test ts[1] == [1, 2, 3]
-end
-
-function test_create_timestep_compactly_add_three_increments_compactly_for_time_t0()
-    ts = TimeStep(1, 2, 3)
-    @test length(ts) == 3
-    @test ts.time == 0.0
-    @test isa(ts[1], Increment)
-end
-
-function test_create_timestep_compactly_add_two_increments()
-    ts = TimeStep([1, 2, 3], [2, 3, 4])
-    @test length(ts) == 2
-    @test ts.time == 0.0
-    @test isa(ts[1], Increment)
-    @test isa(ts[2], Increment)
-    @test ts[1] == [1, 2, 3]
-    @test ts[2] == [2, 3, 4]
-end
-
-function test_create_timesteps_for_different_times()
-    @test TimeStep(0.5, [1, 2]).time == 0.5
-    @test TimeStep(0.5, [1, 2]) == [1, 2]
-    @test TimeStep(0.5, 1).time == 0.5
-    @test TimeStep(0.5, 1) == [1]
-end
-
-function test_default_discrete_field_quick_way_vector()
-    f1 = DefaultDiscreteField([1, 2, 3])
-    @debug("f1 = $f1")
-    @test isa(f1[1], TimeStep)
-    @test isa(f1[1][1], Increment)
-    @test f1[1][1] == [1, 2, 3]
-    @test f1[1].time == 0.0
-end
-
-function test_default_discrete_field_quick_way_scalar()
-    f1 = DefaultDiscreteField(1)
-    @test length(f1) == 1
-    @test isa(f1[1], TimeStep)
-    @test isa(f1[1][1], Increment)
-    @test f1[1][1] == [1]
-    @test f1[1].time == 0.0
-end
-
-function test_default_discrete_field_traditional_way()
-    i1 = Increment([1, 2, 3])
-    i2 = Increment([2, 3, 4])
-    t1 = TimeStep(1.0, Increment[i1, i2])
-    i3 = Increment([2, 3, 4])
-    i4 = Increment([3, 4, 5])
-    t2 = TimeStep(2.0, Increment[i3, i4])
-    timesteps = TimeStep[t1, t2]
-    f1 = DefaultDiscreteField(timesteps)
-    @test length(f1) == 2
-    @test isa(f1, Field)
-    @test f1[1][1] == [1, 2, 3]
-    @test f1[1][2] == [2, 3, 4]
-    @test f1[2][1] == [2, 3, 4]
-    @test f1[2][2] == [3, 4, 5]
-    @test f1[1].time == 1.0
-    @test f1[2].time == 2.0
-end
-
-function test_default_discrete_field_quick_way_two_timesteps_with_constant_value()
-    f1 = DefaultDiscreteField(1, 2)
-    @test length(f1) == 2
-    @test isa(f1[1], TimeStep)
-    @test isa(f1[2], TimeStep)
-    @test isa(f1[1][1], Increment)
-    @test isa(f1[2][1], Increment)
-    @test f1[1][1] == [1]
-    @test f1[2][1] == [2]
-    @test f1[1].time == 0.0
-    @test f1[2].time == 1.0
-end
-
-function test_default_discrete_field_quick_way_two_timesteps_with_vector_value()
-    f1 = DefaultDiscreteField([1, 2, 3], [3, 4, 5])
-    @test length(f1) == 2
-    @test isa(f1[1], TimeStep)
-    @test isa(f1[2], TimeStep)
-    @test isa(f1[1][1], Increment)
-    @test isa(f1[2][1], Increment)
-    @test f1[1][1] == [1, 2, 3]
-    @test f1[2][1] == [3, 4, 5]
-    @test f1[1].time == 0.0
-    @test f1[2].time == 1.0
-end
-
-function test_default_discrete_field_quick_way_set_time_vector_also()
-    f1 = DefaultDiscreteField(
-        (0.5, [1, 2, 3]),
-        (1.0, [3, 4, 5]))
-    @test isa(f1[1], TimeStep)
-    @test isa(f1[2], TimeStep)
-    @test isa(f1[1][1], Increment)
-    @test isa(f1[2][1], Increment)
-    @test f1[1][1] == [1, 2, 3]
-    @test f1[2][1] == [3, 4, 5]
-    @test f1[1].time == 0.5
-    @test f1[2].time == 1.0
-end
-
-function test_default_discrete_field_for_loop()
-    field = DefaultDiscreteField(
-        (0.5, [1, 2, 3]),
-        (1.0, [3, 4, 5]),
-        (1.5, [4, 5, 6]))
-    timesteps = [ts for ts in field]
-    @test timesteps[1].time == 0.5
-    @test timesteps[2].time == 1.0
-    @test timesteps[3].time == 1.5
-    @test timesteps[1][end] == [1, 2, 3]
-    @test timesteps[2][end] == [3, 4, 5]
-    @test timesteps[3][end] == [4, 5, 6]
-end
-
-function test_default_continuous_field()
-
-    function myfield(xi::Vector, time::Float64)
-        time/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])]'
-    end
-
-    f = DefaultContinuousField(myfield)
-    @test f([0.0, 0.0], 1.0) == [0.25 0.25 0.25 0.25]
-    
-end
-
-function test_add_discrete_field_to_fieldset()
-    fs = FieldSet()
-    fs["temperature"] = DefaultDiscreteField([1, 2, 3])
-    @test length(fs) == 1
-    @test fs["temperature"] == [1, 2, 3]
-end
-
-function test_adding_discrete_fields_to_fieldset_quickly()
-    fs = FieldSet()
-    fs["temperature"] = [1, 2, 3, 4]
-    @test fs["temperature"][end][end] == [1, 2, 3, 4]
-    @test last(fs["temperature"]) == [1, 2, 3, 4]
-end
-
-function test_adding_all_kind_of_fields_to_fieldset()
-    fs = FieldSet()
-    fs["constant scalar field"] = 1
-    fs["scalar field"] = [1, 2, 3, 4]
-    fs["vector field"] = reshape(collect(1:8), 2, 4)
-    fs["second order tensor field"] = reshape(collect(1:3*3*4), 3, 3, 4)
-    fs["fourth order tensor field"] = reshape(collect(1:3*3*3*3*4), 3, 3, 3, 3, 4)
-    timestep = fs["vector field"][end]
-    @test fs["vector field"][end].time == 0.0
-end
-
-function test_adding_timesteps()
-    fs = FieldSet()
-    fs["temperature"] = [1, 2, 3, 4]
-    T0 = last(fs["temperature"])  # last increment of last field
-    T1 = Increment(T0 + 1)
-    timestep = TimeStep(1.0, Increment[T1])  # new list of increments for timestep
-    push!(fs["temperature"], timestep)
-    T2 = last(fs["temperature"])
-    @test length(fs["temperature"]) == 2
-    @test last(fs["temperature"]) == [2, 3, 4, 5]
-    @test fs["temperature"][end].time == 1.0
-end
-
-function test_adding_timesteps_compactly()
-    fs = FieldSet()
-    fs["temperature"] = [1, 2, 3, 4]
-    T0 = last(fs["temperature"])
-    T1 = Increment(T0 + 1)
-    push!(fs["temperature"], TimeStep(1.0, T1))
-    @test length(fs["temperature"]) == 2
-    @test last(fs["temperature"]) == [2, 3, 4, 5]
-    @test fs["temperature"][end].time == 1.0
-end
-
-function test_add_several_timesteps_without_time_vector()
-    fs = FieldSet()
-    fs["time series"] = [1, 2, 3, 4], [2, 3, 4, 5]
-    @debug("fieldset = $fs")
-    @test fs["time series"][1].time == 0.0
-    @test fs["time series"][2].time == 1.0
-    @test fs["time series"][1][end] == [1, 2, 3, 4]
-    @test fs["time series"][2][end] == [2, 3, 4, 5]
-end
-
-function test_adding_several_timesteps_at_once_with_time_vector()
-    fs = FieldSet()
-    fs["time series"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5])
-    @test fs3["time series"][1].time == 0.0
-    @test fs3["time series"][2].time == 0.5
-    @test fs3["time series"][1][end] == [1, 2, 3, 4]
-    @test fs3["time series"][2][end] == [2, 3, 4, 5]
-end
-
-function test_adding_continuous_field_to_fieldset()
-    fs = FieldSet()
-    fs["continuous field"] = (xi, t) -> xi[1]*xi[2]*t
-    @test fs["continuous field"]([1.0, 2.0], 3.0) == 6.0
-end
-
-type MyContinuousField <: ContinuousField
-    basis :: Function
-    discrete_field :: DiscreteField
-end
-function Base.call(field::MyContinuousField, xi::Vector, time::Number=1.0)
-    data = last(field.discrete_field) # get the last timestep last increment
-    @debug("data = $data, typeof data = $(typeof(data))")
-    basis = time*field.basis(xi) # evaluate basis at point ξ.
-    sum([basis[i]*data[i] for i=1:length(data)]) # sum results
-end
-function test_continuous_field()
-    fs = FieldSet()
-    fs["discrete field"] = [1, 2, 3, 4]
-    basis(xi) = 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])]
-    fs["continuous field"] = MyContinuousField(basis, fs["discrete field"])
-    @test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(1+2+3+4)
-    T0 = last(fs["discrete field"])
-    T1 = Increment(T0 + 1)
-    push!(fs["discrete field"], TimeStep(1.0, T1))
-    @test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(2+3+4+5)
-end
-
-type MyDiscreteField <: DiscreteField
-    discrete_points :: Vector
-    continuous_field :: ContinuousField
-end
-Base.length(field::MyDiscreteField) = length(field.discrete_points)
-Base.endof(field::MyDiscreteField) = endof(field.discrete_points)
-Base.last(field::MyDiscreteField) = Float64[field[i] for i=1:length(field)]
-function Base.getindex(field::MyDiscreteField, idx::Int64)
-    field.continuous_field(field.discrete_points[idx])
-end
-function test_discrete_field()
-    fs = FieldSet()
-    fs["discrete field"] = [1, 2, 3, 4]
-    basis(xi) = 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])]
-    fs["continuous field"] = MyContinuousField(basis, fs["discrete field"])
-    discrete_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]
-    fs["discrete field 2"] = MyDiscreteField(discrete_points, fs["continuous field"])
-    @test last(fs["discrete field 2"]) ≈ [
-        1.7559830641437073,
-        2.0893163974770410,
-        2.9106836025229590,
-        3.2440169358562922]
-end
-
-function test_field_conversion()
-    i1 = Increment([1, 2, 3])
-    i2 = Increment([2, 3, 4])
-    t1 = TimeStep(1.0, Increment[i1, i2])
-    i3 = Increment([2, 3, 4])
-    i4 = Increment([3, 4, 5])
-    t2 = TimeStep(2.0, Increment[i3, i4])
-    timesteps = TimeStep[t1, t2]
-    info("timesteps = $timesteps")
-    f1 = Field(timesteps)
-    info("field = $f1")
-    @test length(f1) == 2
-    @test isa(f1, Field)
-    @test f1[1][1] == [1, 2, 3]
-    @test f1[1][2] == [2, 3, 4]
-    @test f1[2][1] == [2, 3, 4]
-    @test f1[2][2] == [3, 4, 5]
-    @test f1[1].time == 1.0
-    @test f1[2].time == 2.0
+function test_create_field()
+    f = Field(1.0)
 end
 
 end
diff --git a/test/test_global_assembly.jl b/test/test_global_assembly.jl
index 7318bd3..c8f5349 100644
--- a/test/test_global_assembly.jl
+++ b/test/test_global_assembly.jl
@@ -4,7 +4,7 @@
 module AssemblyTests
 
 using JuliaFEM.Test
-using JuliaFEM: Quad4, Seg2, FieldSet, Field, PlaneHeatProblem
+using JuliaFEM: Quad4, Seg2, FieldSet, Field, HeatProblem
 using JuliaFEM: Assembly, assemble!
 
 """assemble a simple two element problem and solve"""
@@ -21,7 +21,7 @@ function test_assembly()
     el2["temperature flux"] = ((0.0 => 0.0), (1.0 => 600.0))
     info("element created")
 
-    problem = PlaneHeatProblem()
+    problem = HeatProblem()
     info("problem created. pushing elements")
     push!(problem, el1)
     push!(problem, el2)
diff --git a/test/test_heat.jl b/test/test_heat.jl
index 61a60c7..e9a1cab 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -35,6 +35,14 @@ function test_one_element()  # always start test function with name test_
     fdofs = [1, 2]
     A = full(assembly.stiffness_matrix)
     b = full(assembly.force_vector)
+
+    @test isapprox(A, [
+         4.0 -1.0 -2.0 -1.0
+        -1.0  4.0 -1.0 -2.0
+        -2.0 -1.0  4.0 -1.0
+        -1.0 -2.0 -1.0  4.0
+    ])
+
     @test isapprox(A[fdofs, fdofs] \ b[fdofs], [1.0, 1.0])
 
     # Set constant flux g=6 on boundary. Accurate solution is
diff --git a/test/test_randomfields.jl b/test/test_randomfields.jl
deleted file mode 100644
index a93ab11..0000000
--- a/test/test_randomfields.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
-
-module RandomFieldTests
-
-using JuliaFEM
-using JuliaFEM: DiscreteField, Field, Increment, Quad4
-using JuliaFEM.Test
-
-type RandomField <: DiscreteField
-    mu :: Float64
-    std :: Float64
-end
-
-Base.first(field::RandomField) = Increment(randn(2, 4).*field.std^2 + field.mu)
-
-
-function test_interpolate_in_time()
-    r = RandomField(10.0, 0.0)
-    f = Increment(ones(2, 4)*10.0)
-    @test r(0.0) == f
-    @test r(-Inf) == f
-    @test r(+Inf) == f
-    @test r(1.0) == f
-end
-
-function test_interpolate_in_spatial_domain()
-    basis = Quad4([1, 2, 3, 4]).basis
-    r = RandomField(10.0, 0.0)
-    feval = basis(r(0.0), [0.0, 0.0])
-    @test feval == [10.0, 10.0]
-end
-
-end
diff --git a/test/test_types.jl b/test/test_types.jl
index 135a7ba..7f620c2 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -1,22 +1,15 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using JuliaFEM: Basis, Field, FieldSet, interpolate
-using FactCheck
+module TypesTests
 
-facts("test fields") do
-    # multiple field with some constant
-    u1 = Field(0.0, [0.0, 1.0])
-    u2 = 3.0*u1
-    @fact u1.time --> 0.0
-    @fact u2.time --> 0.0
-    @fact u2.values --> [0.0, 3.0]
+using JuliaFEM
+using JuliaFEM.Test
 
-    # addition of fields together
-    u1 = Field(0.0, [0.0, 1.0])
-    u2 = Field(0.0, [1.0, 2.0])
-    u3 = u1 + u2
-    @fact u3.values --> [1.0, 3.0]
+using JuliaFEM: Field, FieldSet
+
+function test_foo()
+    @test 1+1 == 2
 end
 
 #= to be fixed
@@ -85,3 +78,4 @@ end
 
 =#
 
+end