diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 594157f..d747ea4 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -52,8 +52,8 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n)
 end
 
 # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
-include("fields.jl")
-include("basis.jl")  # interpolation of discrete fields
+include("fields2.jl")
+#include("basis.jl")  # interpolation of discrete fields
 include("symbolic.jl") # a thin symbolic layer for fields
 include("types.jl")  # type definitions
 
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 2ba5bfe..0885c94 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -114,7 +114,7 @@ end
 """ 4-node plane stress element. """
 type CPS4 <: PlaneStressElasticityEquation
     element :: Quad4
-    integration_points :: Array{IntegrationPoint, 1}
+    integration_points :: Vector{IntegrationPoint}
 end
 
 function Base.size(equation::CPS4)
@@ -123,7 +123,9 @@ end
 
 function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    haskey(element, "displacement") || (element["displacement"] = zeros(2, 4))
+    if !haskey(element, "displacement")
+        element["displacement"] = 0.0 => [zeros(2) for i=1:4]
+    end
     CPS4(element, integration_points)
 end
 
@@ -139,7 +141,9 @@ end
 
 function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
     integration_points = get_default_integration_points(element)
-    haskey(element, "displacement") || (element["displacement"] = zeros(2, 2))
+    if !haskey(element, "displacement")
+        element["displacement"] = 0.0 => [zeros(2) for i=1:2]
+    end
     CPS2(element, integration_points)
 end
 
diff --git a/src/elements.jl b/src/elements.jl
index 2acd923..527424e 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -75,18 +75,21 @@ Examples
 >>> element["temperature"] = (0.0, [0, 0, 0, 0]), (1.0, [1, 2, 3, 4])
 >>> element["temperature"] = (0.0 => [0, 0, 0, 0], 1.0 => [1, 2, 3, 4])
 """
-function Base.setindex!(element::Element, field_data, field_name)
-    setindex!(element.fields, field_data, field_name)
+function Base.setindex!(element::Element, data, name::ASCIIString)
+    element.fields[name] = Field(data)
+end
+function Base.setindex!(element::Element, data::Tuple, name::ASCIIString)
+    element.fields[name] = Field(data...)
 end
 
-function Base.setindex!(element::Element, field_data::Tuple, field_name)
-    field = Field()
-    for (time, data) in field_data
-        ts = TimeStep(time, Increment[Increment(data)])
-        push!(field, ts)
-    end
-    element[field_name] = field
-end
+#function Base.setindex!(element::Element, field_data::Tuple, field_name)
+#    field = Field()
+#    for (time, data) in field_data
+#        ts = TimeStep(time, Increment[Increment(data)])
+#        push!(field, ts)
+#    end
+#    element[field_name] = field
+#end
 
 function get_connectivity(el::Element)
     return el.connectivity
@@ -95,12 +98,12 @@ end
 abstract AbstractFunctionSpace
 
 type FunctionSpace <: AbstractFunctionSpace
-    basis :: Basis
+    basis :: CVTI
     fields :: FieldSet
 end
 
 type GradientFunctionSpace <: AbstractFunctionSpace
-    basis :: Basis
+    basis :: CVTI
     fields :: FieldSet
 end
 
@@ -130,9 +133,6 @@ end
 """ Evaluate field on element function space. """
 function call(u::FunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
     field = !isa(variation, Void) ? variation : u.fields[field_name](t)
-    if length(field) == 1
-        return field.data[1]
-    end
     u.basis(field, xi)
 end
 
@@ -171,7 +171,7 @@ end
 """ Get a determinant of element in point ξ. """
 function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
     X = u.fields["geometry"](time)
-    dN = u.basis.dbasisdxi(xi)
+    dN = u.basis(xi, Val{:grad})
     J = sum([dN[:,i]*X[i]' for i=1:length(X)])
     m, n = size(J)
     return m == n ? det(J) : norm(J)
diff --git a/src/equations.jl b/src/equations.jl
index 3190025..bd8e206 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -161,7 +161,8 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
 
     # 3. virtual work -- user has defined some residual r = p - f = 0
     if has_residual_vector(equation)
-        field = element[unknown_field_name](time)
+
+        field = DVTI(last(element[unknown_field_name]).data)
 
         """ Wrapper for virtual work for ForwardDiff. """
         function calc_R(data::Vector)
@@ -177,9 +178,12 @@ function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, pro
             if haskey(element, "$unknown_field_name nodal load")
                 R -= vec(element["$unknown_field_name nodal load"](time))
             end
+            #info("return = $R")
             return R
         end
-       
+
+        #info("field = $field")
+        #info("vec(field) = $(vec(field))")
         jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
         add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
         add!(assembly.force_vector, gdofs, -ForwardDiff.value(allresults))
diff --git a/src/fields2.jl b/src/fields2.jl
index 14bbfe5..d36713a 100644
--- a/src/fields2.jl
+++ b/src/fields2.jl
@@ -16,7 +16,48 @@ type Field{A<:Union{Discrete,Continuous}, B<:Union{Constant,Variable}, C<:Union{
     data
 end
 
-# Different field combinations
+### Basic data structure for discrete field
+
+type Increment{T}
+    time :: Float64
+    data :: T
+end
+
+function Base.convert{T}(::Type{Increment{T}}, data::Pair{Float64,T})
+    return Increment{T}(data[1], data[2])
+end
+
+function Base.convert{T}(::Type{Increment{Vector{Vector{T}}}}, data::Pair{Float64, Matrix{T}})
+    time = data[1]
+    content = data[2]
+    return Increment(time, Vector{T}[content[:,i] for i=1:size(content,2)])
+end
+
+function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
+    return increment.data[i]
+end
+
+function Base.(:*)(d, increment::Increment)
+    return d*increment.data
+end
+
+### Basic data structure for continuous field
+
+type Basis
+    basis :: Function
+    dbasis :: Function
+end
+
+function Base.call(basis::Basis, xi::Vector)
+    basis.basis(xi)
+end
+
+function Base.call(basis::Basis, xi::Vector, ::Type{Val{:grad}})
+    basis.dbasis(xi)
+end
+
+### Different field combinations and other typealiases
+
 typealias DCTI Field{Discrete,   Constant, TimeInvariant}
 typealias DVTI Field{Discrete,   Variable, TimeInvariant}
 typealias DCTV Field{Discrete,   Constant, TimeVariant}
@@ -26,71 +67,207 @@ typealias CVTI Field{Continuous, Variable, TimeInvariant} # can be used to inter
 typealias CCTV Field{Continuous, Constant, TimeVariant} # can be used to interpolate in time
 typealias CVTV Field{Continuous, Variable, TimeVariant}
 
-# Basic data structure for discrete field
-type Increment{T}
-    time :: Float64
-    data :: T
-end
+typealias ScalarIncrement{T} Increment{T}
+typealias VectorIncrement{T} Increment{Vector{T}}
+typealias TensorIncrement{T} Increment{Matrix{T}}
 
-typealias VectorIncrement Increment{Vector}
+typealias DiscreteField      Union{DCTI, DVTI, DCTV, DVTV}
+typealias ContinuousField    Union{CCTI, CVTI, CCTV, CVTV}
+typealias ConstantField      Union{DCTI, DCTV, CCTI, CCTV}
+typealias VariableField      Union{DVTI, DVTV, CVTI, CVTV}
+typealias TimeInvariantField Union{DCTI, DVTI, CCTI, CVTI}
+typealias TimeVariantField   Union{DCTV, DVTV, CCTV, CVTV}
 
-function Base.getindex{T}(increment::Increment{Vector{T}}, i::Int64)
-    return increment.data[i]
-end
 
-# Basic data structure for continuous field
-type Basis
-    basis :: Function
-    dbasis :: Function
-end
+### Convenient functions to create fields
 
-# Functions simplifying definition of fields.
-
-"""
-All other data than vectors are considered as constant time invariant fields.
-"""
 function Field(data)
-    DCTI(data)
+    return DCTI(data)
 end
 
-"""
-Vector data is considered as variable field time invariant field.
-"""
 function Field(data::Vector)
-    DVTI(data)
+    return DVTI(data)
 end
 
-"""
-Data given in (time, value) pairs, where value is not vector, is considered as
-constant time variant field.
-"""
 function Field{T}(data::Pair{Float64, T}...)
-    increments = [Increment{T}(d[1], d[2]) for d in data]
-    DCTV(increments)
+    return DCTV([Increment{T}(d[1], d[2]) for d in data])
 end
 
-"""
-Data given in (time, value) pairs, where value is a vector, is considered as
-variable time variant field.
-"""
 function Field{T}(data::Pair{Float64, Vector{T}}...)
-    increments = [Increment{Vector{T}}(d[1], d[2]) for d in data]
-    DVTV(increments)
+    return DVTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
-""" Special case, constant time-variant vector, converted automatically. """
 function Base.convert{T}(::Type{DCTV}, data::Pair{Float64, Vector{T}}...)
-    increments = [Increment(d[1], d[2]) for d in data]
-    DCTV(increments)
+    return DCTV([Increment{Vector{T}}(d[1], d[2]) for d in data])
 end
 
-## Other field related functions
+function CVTI(basis::Function, dbasis::Function)
+    return CVTI(Basis(basis, dbasis))
+end
+
+function Field(basis::Function, dbasis::Function)
+    return CVTI(basis, dbasis)
+end
+
+### Accessing and manipulating discrete fields
 
 function Base.getindex(field::DVTV, i::Int64)
     return field.data[i]
 end
 
+function Base.push!(field::DCTV, data::Pair)
+    push!(field.data, data)
+end
+
+function Base.push!(field::DVTV, data::Pair)
+#    info("field.data = \n$(field.data)")
+#    info("data = \n$data")
+    push!(field.data, data)
+end
+
+function Base.getindex(field::DVTV, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::DVTI, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::DCTV, i::Int64)
+    return field.data[i]
+end
+
+function Base.getindex(field::Field, i::Int64)
+    return field.data[i]
+end
+
+function Base.length(field::DVTI)
+    return length(field.data)
+end
+
+function Base.length(field::DCTI)
+    return 1
+end
+
+function Base.length(field::DVTV)
+    return length(field.data)
+end
+
+function Base.length(field::DCTV)
+    return length(field.data)
+end
+
+for op = (:+, :*, :/, :-)
+    @eval ($op)(increment::Increment, field::DCTI) = ($op)(increment.data, field.data)
+    @eval ($op)(field::DCTI, increment::Increment) = ($op)(increment.data, field.data)
+    @eval ($op)(field1::DCTI, field2::DCTI) = ($op)(field1.data, field2.data)
+    @eval ($op)(field::DCTI, k) = ($op)(field.data, k)
+    @eval ($op)(k, field::DCTI) = ($op)(field.data, k)
+end
+
+function Base.vec(field::DVTI)
+    return [field.data...;]
+end
+
+function Base.vec(field::DCTV)
+    info("trying to vectorize $field")
+    error("does not make sense")
+end
+
+function Base.endof(field::Field)
+    return endof(field.data)
+end
+
+#function Base.similar{T}(field::DVTI, data::Vector{T})
+#    return Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
+#end
+
+function Base.similar{T}(field::DVTI, data::Vector{T})
+    n = length(field.data)
+    data = reshape(data, round(Int, length(data)/n), n)
+    newdata = Vector[data[:,i] for i=1:n]
+    return typeof(field)(newdata)
+end
+
+### Accessing continuous fields
+
+function Base.call(field::CVTI, xi::Vector)
+    field.data(xi)
+end
+
+function Base.call(field::CVTI, xi::Vector, ::Type{Val{:grad}})
+    field.data(xi, Val{:grad})
+end
+
+function Base.convert(::Type{Basis}, field::CVTI)
+    return field.data
+end
+
+### Interpolation 
+
+function Base.call(field::DVTI, time::Float64)
+    # interpolating time-invariant field in time direction -> no effect
+    return field
+end
+
+function Base.call(field::DCTI, time::Float64)
+    # interpolating time-invariant field in time direction -> no effect
+    return field
+end
+
+function Base.call(basis::CVTI, field::DCTI, xi::Vector)
+    # try to interpolate constant value -> no effect
+    return field
+end
+
+#function Base.call(basis::Basis, field::DCTI, xi::Vector)
+    # calling constant field with basis -> no effect
+#    return field
+#end
+
+function Base.call(basis::CVTI, values::DVTI, xi::Vector)
+    N = basis(xi)
+    return sum([N[i]*values[i] for i=1:length(N)])
+end
+
+function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
+    dbasis = basis(xi, Val{:grad})
+    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+    invJ = isa(J, Vector) ? inv(J[1]) : inv(J)
+    grad = invJ * dbasis
+    return grad
+end
+
+function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
+    grad = call(basis, geometry, xi, Val{:grad})
+    gradf = sum([grad[:,i]*values[i]' for i=1:length(geometry)])'
+    return length(gradf) == 1 ? gradf[1] : gradf
+end
+
+function Base.call(field::DCTV, time::Float64)
+    for i in length(field)
+        if isapprox(field[i].time, time)
+            return DCTI(field[i].data)
+        end
+    end
+    error("interpolate DCTV: not implemented yet")
+end
+
+function Base.call(field::DVTV, time::Float64, time_extrapolation::Symbol=:linear)
+#    info("length of field DVTV: $(length(field))")
+    for i=reverse(1:length(field))
+        res = isapprox(field[i].time, time)
+        #info("isapprox $(field[i].time) to $time ? $res")
+        if isapprox(field[i].time, time)
+            return DVTI(field[i].data)
+        end
+    end
+    info(field.data)
+    info(time)
+    error("interpolate DVTV: not implemented yet")
+end
 
 ### FIELDSET ###
 
 typealias FieldSet Dict{ASCIIString, Field}
+
diff --git a/src/heat.jl b/src/heat.jl
index e1e746d..e7d63cf 100644
--- a/src/heat.jl
+++ b/src/heat.jl
@@ -60,6 +60,7 @@ function assemble!(assembly::Assembly, equation::HeatEquation, time::Number=0.0,
             add!(assembly.force_vector, gdofs, w*N'*f)
         end
         if haskey(element, "temperature flux")
+            info("assemble boundary flux")
             g = basis("temperature flux", ip, time)
             add!(assembly.force_vector, gdofs, w*N'*g)
         end
@@ -93,13 +94,13 @@ end
 
 function Base.convert(::Type{HeatEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    haskey(element, "temperature") || (element["temperature"] = zeros(4))
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
     DC2D4(element, integration_points)
 end
 
 function Base.convert(::Type{HeatEquation}, element::Seg2)
     integration_points = get_default_integration_points(element)
-    haskey(element, "temperature") || (element["temperature"] = zeros(2))
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
     DC2D2(element, integration_points)
 end
 
diff --git a/src/lagrange.jl b/src/lagrange.jl
index 2d5a824..c97d630 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -36,11 +36,11 @@ macro create_lagrange_element(element_name, element_description, X, P)
         basis, dbasisdxi = calculate_lagrange_basis($P, $X)
         type $eltype <: CG
             connectivity :: Array{Int, 1}
-            basis :: Basis
+            basis :: CVTI
             fields :: FieldSet
         end
         function $eltype(connectivity, args...)
-            $eltype(connectivity, Basis(basis, dbasisdxi), FieldSet())
+            $eltype(connectivity, CVTI(basis, dbasisdxi), FieldSet())
         end
         get_element_description(el::Type{$eltype}) = $element_description
         Base.size(element::Type{$eltype}) = Base.size($X)
diff --git a/src/solvers.jl b/src/solvers.jl
index f0cda11..0ae9831 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -9,8 +9,7 @@ abstract Solver
 Solve field equations for single element with some dofs fixed. This can be used
 to test nonlinear element formulations.
 """
-function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number=0.0; 
-                 max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     unknown_field_name = get_unknown_field_name(equation)
     element = get_element(equation)
     x0 = element[unknown_field_name](0.0)
@@ -28,7 +27,9 @@ function solve!(equation::Equation, free_dofs::Vector{Int}, time::Number=0.0;
         end
         dx[free_dofs] = A \ b
         x += dx
-        push!(element[unknown_field_name], reshape(x, size(equation)))
+        eqsize = size(equation)
+        data = eqsize[1] != 1 ? reshape(x, eqsize) : x
+        push!(element[unknown_field_name], time => data)
         norm(dx) < tolerance && return
     end
     error("Did not converge in $max_iterations iterations")
@@ -40,7 +41,7 @@ to test nonlinear element formulations. Dirichlet boundary is assumed to be homo
 and degrees of freedom are eliminated. So if boundary condition is known in nodal
 points and everything is zero this should be quite good.
 """
-function solve!(problem::Problem, free_dofs::Vector{Int}, time::Number=1.0; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+function solve!(problem::Problem, free_dofs::Vector{Int}, time::Float64; max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     info("start solver")
     assembly = Assembly()
 #   x = zeros(ga.ndofs)
@@ -68,8 +69,12 @@ function solve!(problem::Problem, free_dofs::Vector{Int}, time::Number=1.0; max_
         for equation in get_equations(problem)
             element = get_element(equation)
             gdofs = get_gdofs(equation)
-            data = reshape(full(x[gdofs]), size(equation))
-            push!(element[field_name], data)
+            data = full(x[gdofs])
+            eqsize = size(equation)
+            if eqsize[1] != 1
+                data = reshape(data, eqsize)
+            end
+            push!(element[field_name], time => data)
         end
         norm(dx) < tolerance && return
     end
diff --git a/src/symbolic.jl b/src/symbolic.jl
index c03eac0..4822b73 100644
--- a/src/symbolic.jl
+++ b/src/symbolic.jl
@@ -3,7 +3,7 @@
 
 # https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
 
-immutable SymbolicField <: Field
+immutable SymbolicField <: AbstractField
     name :: ASCIIString
 end
 
diff --git a/src/types.jl b/src/types.jl
index 9fa6c78..41e45d3 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -25,25 +25,45 @@ type IntegrationPoint
 end
 
 function IntegrationPoint(xi, weight)
-    IntegrationPoint(xi, weight, Dict())
+    return IntegrationPoint(xi, weight, FieldSet())
 end
 
 function Base.convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
 
-function Base.call(basis::Basis, ip::IntegrationPoint)
-    return basis(ip.xi)
+function Base.call(field::CVTI, ip::IntegrationPoint)
+    return call(field, ip.xi)
 end
 
-function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
-    return call(basis, increment, ip.xi)
+function Base.call(basis::CVTI, field::DCTI, ip::IntegrationPoint)
+    call(basis, field, ip.xi)
 end
 
-function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-    return call(basis, increment, ip.xi, Val{:grad})
+function Base.call(basis::CVTI, field::DVTI, ip::IntegrationPoint, ::Type{Val{:grad}})
+    call(basis, field, ip.xi, Val{:grad})
 end
 
-function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
-    return call(basis, geometry, field, ip.xi, Val{:grad})
+function Base.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}})
+    call(basis, geometry, field, ip.xi, Val{:grad})
+end
+
+#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
+#    return call(basis, increment, ip.xi)
+#end
+#function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, increment, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, field, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+#    return call(basis, geometry, field, ip.xi, Val{:grad})
+#end
+#function Base.call(basis::Basis, field::Field, ip::IntegrationPoint)
+#    return call(basis, field, ip.xi)
+#end
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index cb99e6e..d91b506 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -17,8 +17,8 @@ function test_elasticity_volume_load()
     free_dofs = [3, 4, 5, 6]
     problem = PlaneStressElasticityProblem()
     push!(problem, element)
-    solve!(problem, free_dofs; max_iterations=10)
-    disp = get_basis(element)("displacement", [1.0, 1.0])
+    solve!(problem, free_dofs, 0.0; max_iterations=10)
+    disp = get_basis(element)("displacement", [1.0, 1.0], 0.0)
     info("displacement at tip: $disp")
     # verified using Code Aster.
     @test isapprox(disp[2], -8.77303119819776)
@@ -39,11 +39,13 @@ function test_elasticity_surface_load()
     problem = PlaneStressElasticityProblem()
     push!(problem, element1)
     push!(problem, element2)
-    solve!(problem, free_dofs; max_iterations=10)
-    disp = get_basis(element1)("displacement", [1.0, 1.0])[2]
+    solve!(problem, free_dofs, 1.0; max_iterations=10)
+    disp = get_basis(element1)("displacement", [1.0, 1.0], 1.0)[2]
     info("displacement at tip: $disp")
     # verified using Code Aster.
     @test isapprox(disp, -9.33106637611714)
 end
 
+#test_elasticity_volume_load()
+
 end
diff --git a/test/test_heat.jl b/test/test_heat.jl
index d47b33b..8dd4408 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -14,6 +14,7 @@ function test_one_element()  # always start test function with name test_
 
     # volume element
     element = Quad4([1, 2, 3, 4])
+
     element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
     element["temperature thermal conductivity"] = 6.0
     element["temperature load"] = [12.0, 12.0, 12.0, 12.0]
@@ -22,14 +23,12 @@ function test_one_element()  # always start test function with name test_
     # boundary element 
     boundary_element = Seg2([1, 2])
     boundary_element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
-    # linear ramp from 1 to 6 in time 0 to 1
-    boundary_element["temperature flux"] = (0.0, 0.0), (1.0, 6.0)
+    # linear ramp from 0 to 6 in time 0 to 1
+    boundary_element["temperature flux"] = (0.0 => 0.0, 1.0 => 6.0)
 
     # Set constant source f=12 with k=6. Accurate solution is
     # T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
     equation = convert(HeatEquation, element)
-    #la = initialize_local_assembly()
-    #calculate_local_assembly!(la, equation, "temperature")
     assembly = Assembly()
     assemble!(assembly, equation)
     fdofs = [1, 2]
@@ -44,13 +43,7 @@ function test_one_element()  # always start test function with name test_
 
     time = 1.0
     assemble!(assembly, equation, time)
-    info("after first element: $(length(assembly.force_vector.V))")
-    info(full(assembly.force_vector)')
     assemble!(assembly, boundary_equation, time)
-    info("after second element: $(length(assembly.force_vector.V))")
-    info(full(assembly.force_vector)')
-    #calculate_local_assembly!(la, boundary_equation, "temperature")
-    #b = la.force_vector
     A = full(assembly.stiffness_matrix)
     b = full(assembly.force_vector)
     T = A[fdofs, fdofs] \ b[fdofs]
diff --git a/test/test_potential_energy.jl b/test/test_potential_energy.jl
index f461376..3890767 100644
--- a/test/test_potential_energy.jl
+++ b/test/test_potential_energy.jl
@@ -23,14 +23,6 @@ type DC2D4NL <: MyEquation
     integration_points :: Vector{IntegrationPoint}
 end
 
-function DC2D4NL(element::Quad4)
-    integration_points = get_default_integration_points(element)
-    if !haskey(element, "temperature")
-        element["temperature"] = zeros(4)
-    end
-    DC2D4NL(element, integration_points)
-end
-
 function Base.size(equation::DC2D4NL)
     return (1, 4)
 end
@@ -41,18 +33,23 @@ type DC2D2NL <: MyEquation
     integration_points :: Vector{IntegrationPoint}
 end
 
-function DC2D2NL(element::Seg2)
-    integration_points = JuliaFEM.line5()
-    if !haskey(element, "temperature")
-        element["temperature"] = zeros(2)
-    end
-    DC2D2NL(element, integration_points)
-end
-
 function Base.size(equation::DC2D2NL)
     return (1, 2)
 end
 
+function Base.convert(::Type{MyEquation}, element::Quad4)
+    integration_points = get_default_integration_points(element)
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(4))
+    DC2D4NL(element, integration_points)
+end
+
+function Base.convert(::Type{MyEquation}, element::Seg2)
+    integration_points = JuliaFEM.line5()
+    haskey(element, "temperature") || (element["temperature"] = 0.0 => zeros(2))
+    DC2D2NL(element, integration_points)
+end
+
+
 """ Calculate a potential Π = Wint - Wext of system. """
 function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing)
     element = get_element(equation)
@@ -89,27 +86,13 @@ function test_potential_energy_method()
     element["temperature load"] = [0.0, 0.0, 0.0, 0.0]
     element["temperature nodal load"] = [3.0, 3.0, 0.0, 0.0]
     element["temperature nonlinearity coefficient"] = 6.0
-    equation = DC2D4NL(element)
+    equation = convert(MyEquation, element)
     # create model -- end
 
-    ass = Assembly()
-    info("unknown field name: $(get_unknown_field_name(equation))")
-
-    T = zeros(4) # create workspace for solution vector
-    dT = zeros(4) # 
-    fd = [1, 2] # free dofs
-    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.
-    for i=1:10
-        empty!(ass)
-        assemble!(ass, equation)  # calculate local matrices
-        dT[fd] = full(ass.stiffness_matrix)[fd,fd] \ full(ass.force_vector)[fd]
-        T += dT
-        push!(element["temperature"], T) # add new increment to model
-        @printf("increment %2d, |du| = %8.5f\n", i, norm(dT))
-        err = last(element["temperature"])[1] - 2/3
-        isapprox(err, 0.0) && break
-    end
-    err = last(element["temperature"])[1] - 2/3
+    solve!(equation, [1, 2], 0.0)
+    basis = get_basis(element)
+    temp = basis("temperature", [0.0, -1.0], 0.0)
+    err = temp - 2/3
     info("error: $err")
     @test isapprox(err, 0.0)
 end
@@ -118,15 +101,11 @@ end
 type TestProblem <: Problem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
-    equations :: Vector{Equation}
-    element_mapping :: Dict{DataType, DataType}
+    equations :: Vector{MyEquation}
 end
 
 function TestProblem(equations=[])
-    element_mapping = Dict(
-        Quad4 => DC2D4NL,
-        Seg2 => DC2D2NL)
-    TestProblem("temperature", 1, equations, element_mapping)
+    TestProblem("temperature", 1, equations)
 end
 
 function test_potential_energy_method_2()
@@ -138,41 +117,23 @@ function test_potential_energy_method_2()
     element1["temperature thermal conductivity"] = 6.0
     element1["temperature load"] = [0.0, 0.0, 0.0, 0.0]
     element1["temperature nonlinearity coefficient"] = [0.0, 0.0, 0.0, 0.0]
-    element1["temperature"] = ones(4)
 
     element2 = Seg2([1, 2])
     element2["geometry"] = Vector[N[1], N[2]]
     element2["temperature coefficient"] = 3.0e-8 # ~ 5.7e-8 * 0.5
     element2["temperature external"] = 100.0
-    element2["temperature"] = ones(2)
     # create model -- end
-    
-    equation1 = DC2D4NL(element1)
-    equation2 = DC2D2NL(element2)
- 
-    ass = Assembly()
-    info("unknown field name: $(get_unknown_field_name(equation1))")
 
-    T = zeros(4) # create workspace for solution vector
-    dT = zeros(4) # 
-    fd = [1, 2] # free dofs
-    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.
-    for i=1:10
-        empty!(ass)
-        assemble!(ass, equation1)
-        assemble!(ass, equation2)
-        dT[fd] = full(ass.stiffness_matrix)[fd,fd] \ full(ass.force_vector)[fd]
-        T += dT
-        push!(element1["temperature"], T)
-        push!(element2["temperature"], T[fd])
-        @printf("increment %2d, |du| = %8.5f\n", i, norm(dT))
-        err = last(element1["temperature"])[1] - 0.5
-        isapprox(err, 0.0) && break
-    end
+    problem = TestProblem()
+    push!(problem, element1)
+    push!(problem, element2)
+    solve!(problem, [1, 2], 0.0)
 
-    err = last(element1["temperature"])[1] - 0.5
+    basis = get_basis(element1)
+    temp = basis("temperature", [0.0, -1.0], 0.0)
+    err = temp - 0.5
     info("error: $err")
-    @test isapprox(err, 0.0)
+    @test isapprox(err, 0.0, atol=1.0e-6)
     
     # @test isapprox(temp, 2.93509690572300E+00)  # tested using Code Aster
 end
diff --git a/test/test_virtual_work.jl b/test/test_virtual_work.jl
index 3957687..0713004 100644
--- a/test/test_virtual_work.jl
+++ b/test/test_virtual_work.jl
@@ -5,8 +5,7 @@ module TestAutoDiffWeakForm
 
 using JuliaFEM.Test
 using JuliaFEM
-using JuliaFEM: Quad4, Equation, IntegrationPoint, assemble!,
-                Assembly,
+using JuliaFEM: Quad4, Equation, IntegrationPoint, assemble!, Assembly,
                 solve!, get_field, get_element, get_basis,
                 grad, get_default_integration_points
 
@@ -23,7 +22,7 @@ end
 function CPS4(element::Quad4)
     integration_points = get_default_integration_points(element)
     if !haskey(element, "displacement")
-        element["displacement"] = zeros(2, 4)
+        element["displacement"] = 0.0 => Vector{Float64}[[0.0,0.0], [0.0,0.0], [0.0,0.0], [0.0,0.0]]
     end
     CPS4(element, integration_points)
 end
@@ -71,8 +70,8 @@ function test_residual_form()
     # create model -- end
 
     free_dofs = [3, 4, 5, 6]
-    solve!(equation, free_dofs)  # launch a newton solver for single element
-    disp = get_basis(element)("displacement", [1.0, 1.0])[2]
+    solve!(equation, free_dofs, 0.0)  # launch a newton solver for single element
+    disp = get_basis(element)("displacement", [1.0, 1.0], 0.0)[2]
     println("displacement at tip: $disp")
     # verified using Code Aster.
     @test isapprox(disp, -8.77303119819776E+00)