diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 7bd7075..2c5484f 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -9,8 +9,10 @@ module JuliaFEM
 importall Base
 using ForwardDiff
 using JLD
-autodiffcache = ForwardDiffCache()
-# export derivative, jacobian, hessian
+
+#Grad = Val{:Grad}
+#detJ = Val{:detJ}
+#export Grad, detJ
 
 include("common.jl")
 
@@ -169,4 +171,5 @@ module Interfaces
 include("interfaces.jl")
 end
 
+
 end # module
diff --git a/src/assembly.jl b/src/assembly.jl
index 6e8de19..79a6654 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -44,17 +44,30 @@ end
 function assemble_posthook!
 end
 
-function assemble!(problem::Problem, time::Real)
+function assemble!(problem::Problem, time=0.0; auto_initialize=true)
     if !isempty(problem.assembly)
-        warn("problem.assembly is not empty and assembling, are you sure you know what are you doing?")
+        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
     end
-    if method_exists(assemble_prehook!, Tuple{typeof(problem), Real})
+    if isempty(problem.elements)
+        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
+    else
+        first_element = first(problem.elements)
+        unknown_field_name = get_unknown_field_name(problem)
+        if !haskey(first_element, unknown_field_name)
+            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
+            if auto_initialize
+                info("Initializing problem $(problem.name) at time $time automatically.")
+                initialize!(problem, time)
+            end
+        end
+    end
+    if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64})
         assemble_prehook!(problem, time)
     end
     for element in get_elements(problem)
         assemble!(problem.assembly, problem, element, time)
     end
-    if method_exists(assemble_posthook!, Tuple{typeof(problem), Real})
+    if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
         assemble_posthook!(problem, time)
     end
 end
diff --git a/src/elements.jl b/src/elements.jl
index 8191ee5..9cd92f6 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -29,10 +29,18 @@ function setindex!(element::Element, data, field_name::ASCIIString)
     element.fields[field_name] = Field(data)
 end
 
+function call(element::Element, field_name::ASCIIString)
+    return element[field_name]
+end
+
 function call(element::Element, field_name::ASCIIString, time)
     return element[field_name](time)
 end
 
+function last(element::Element, field_name::ASCIIString)
+    return last(element[field_name])
+end
+
 function call(element::Element, ip, time)
     return get_basis(element, ip, time)
 end
@@ -80,6 +88,10 @@ function call(element::Element, field::DCTI, ip, time::Float64)
     return field.data
 end
 
+function call(element::Element, field::DCTV, ip, time::Float64)
+    return field(time).data
+end
+
 function call(element::Element, field::CVTV, ip, time::Float64)
     return field(ip, time)
 end
@@ -89,7 +101,9 @@ function call(element::Element, field::Field, ip, time::Float64)
     basis = element(ip, time)
     n = length(element)
     m = length(field_)
-    @assert n == m
+    if n != m
+        error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+    end
     return sum([field_[i]*basis[i] for i=1:n])
 end
 
@@ -118,7 +132,7 @@ function update!{K,V}(element::Element, field_name::ASCIIString, data::Pair{Floa
     update!(element, field_name, time => element_data)
 end
 
-function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair{Float64, Union{Real, Vector{Any}}}}...)
+function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
     for data in datas
         if haskey(element, field_name)
             update!(element[field_name], data)
@@ -132,6 +146,12 @@ function update!(element::Element, field_name::ASCIIString, datas::Union{Real, V
     end
 end
 
+function update!(element::Element, field_name::ASCIIString, datas::Pair...)
+    for data in datas
+        update!(element, field_name, data)
+    end
+end
+
 function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Vector{Any}})
     if haskey(element, field_name)
         update!(element[field_name], data)
@@ -200,11 +220,12 @@ function update!(elements::Vector, field_name::ASCIIString, data)
     end
 end
 
+dbasis_cache = ForwardDiff.jacobian
 """ Evaluate partial derivatives of basis functions using ForwardDiff. """
 function get_dbasis(element::Element, ip, time)
     xi = isa(ip, IP) ? ip.coords : ip
     basis(xi) = vec(get_basis(element, xi, time))
-    return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)'
+    return ForwardDiff.jacobian(basis, xi)'
 end
 
 """ Check existence of field. """
diff --git a/src/heat.jl b/src/heat.jl
index 82bccbb..647cd4f 100644
--- a/src/heat.jl
+++ b/src/heat.jl
@@ -22,23 +22,59 @@ where
     f = temperature load                    defined on volume elements
     g = temperature flux                    defined on boundary elements
 
+Parameters
+----------
+temperature thermal conductivity
+temperature load
+temperature flux
+
+thermal conductivity
+heat source
+heat flux
+heat transfer coefficient
+external temperature
+
+
+Formulations
+------------
+1D, 2D, 3D
+
 References
 ----------
 https://en.wikipedia.org/wiki/Heat_equation
-
+https://en.wikipedia.org/wiki/Heat_capacity
+https://en.wikipedia.org/wiki/Thermal_conduction
+https://en.wikipedia.org/wiki/Thermal_conductivity
+https://en.wikipedia.org/wiki/Thermal_diffusivity
+https://en.wikipedia.org/wiki/Volumetric_heat_capacity
 """
 type Heat <: FieldProblem
+    formulation :: ASCIIString
+end
+
+function Heat()
+    return Heat("3D")
 end
 
 function get_unknown_field_name(problem::Problem{Heat})
     return "temperature"
 end
 
-function get_unknown_field_type(problem::Problem{Heat})
-    return Float64
+function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
+    formulation = Val{Symbol(problem.properties.formulation)}
+    assemble!(assembly, problem, element, time, formulation)
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
+# 3d heat problems
+
+function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
+    info("Unknown element type $E for 3d heat problem!")
+end
+
+typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8}
+typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4}
+
+function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
     gdofs = get_gdofs(problem, element)
     field_name = get_unknown_field_name(problem)
     nnodes = length(element)
@@ -57,6 +93,85 @@ function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element,
             f = element("$field_name load", ip, time)
             fq += w*N'*f
         end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        if haskey(element, "$field_name flux")
+            q = element("$field_name flux", ip, time)
+            fq += w*N'*q
+        end
+        if haskey(element, "$field_name heat transfer coefficient")
+            h = element("$field_name heat transfer coefficient", ip, time)
+            Tu = element("$field_name external temperature", ip, time)
+            K += w*h*N'*N
+            fq += w*N'*h*Tu
+        end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+# 2d heat problems
+
+function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    info("Unknown element type $E for 2d heat problem!")
+end
+
+typealias Heat2DVolumeElements Union{Tri3, Tri6, Quad4}
+typealias Heat2DSurfaceElements Union{Seg2, Seg3}
+
+function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        if haskey(element, "$field_name thermal conductivity")
+            dN = element(ip, time, Val{:Grad})
+            k = element("$field_name thermal conductivity", ip, time)
+            K += w*k*dN'*dN
+        end
+        if haskey(element, "$field_name load")
+            f = element("$field_name load", ip, time)
+            fq += w*N'*f
+        end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
         if haskey(element, "$field_name flux")
             g = element("$field_name flux", ip, time)
             fq += w*N'*g
diff --git a/src/lagrange_macro.jl b/src/lagrange_macro.jl
index 043f0fb..e9b7ab2 100644
--- a/src/lagrange_macro.jl
+++ b/src/lagrange_macro.jl
@@ -95,17 +95,36 @@ end
     (xi) -> [1.0, xi[1], xi[2], xi[1]^2, xi[2]^2, xi[1]*xi[2]])
 
 @create_lagrange_element(Quad4, "4 node bilinear quadrangle element",
-    [-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]])
 
+@create_lagrange_element(Quad9, "9 node bilinear quadrangle element",
+    [-1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0
+     -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0],
+    (xi) -> [1.0,     xi[1],   xi[2],         xi[1]*xi[2],
+             xi[1]^2, xi[2]^2, xi[1]^2*xi[2], xi[1]*xi[2]^2])
+
 # 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]])
+    (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(Hex20, "20 node hexahedra",
+    [
+        -1.0  1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0
+        -1.0 -1.0  1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.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  0.0  0.0  0.0  0.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]
+             x[1]^2, 
+    ])
+=#
 
 @create_lagrange_element(Tet4, "4 node tetrahedron",
     [0.0 1.0 0.0 0.0
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index 4555607..429d980 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -309,21 +309,17 @@ function get_element_sets(med::MEDFile, mesh_name)
     return es
 end
 
-# hex8 nodes rotating cw first in yz plane then x+1
-
 global const med_elmap = Dict{Symbol, Vector{Int}}(
+    :PO1 => [1],
     :SE2 => [1, 2],
     :SE3 => [1, 2, 3],
     :TR3 => [1, 2, 3],
-    :TR6 => [1, 2, 3, 4, 5, 6],
     :QU4 => [1, 2, 3, 4],
-    :HE8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
     :TE4 => [3, 2, 1, 4],
-    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8],
-    :PO1 => [1]
-#   :T10 => [3, 4, 1, 2, 10, 8, 7, 6, 9, 5]
-#   :T10 => [5, 9, 6, 7, 8, 10, 2, 1, 4, 3]
-)
+    :TR6 => [1, 2, 3, 4, 5, 6],
+    :QU8 => [1, 2, 3, 4, 5, 6, 7, 8],
+    :HE8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
+    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8])
 
 function get_connectivity(med::MEDFile, elsets, mesh_name)
     elsets[0] = :OTHER
diff --git a/src/problems.jl b/src/problems.jl
index f1269d3..efedf4f 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -145,7 +145,7 @@ end
     take last known value and set it as a initial quess for next
     time increment.
 """
-function initialize!(problem::Problem, time::Real)
+function initialize!(problem::Problem, time=0.0)
     field_name = get_unknown_field_name(problem)
     field_dim = get_unknown_field_dimension(problem)
     for element in get_elements(problem)
@@ -289,8 +289,12 @@ function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P})
     return problem.parent_field_name
 end
 
-function push!(problem::Problem, element)
-    push!(problem.elements, element)
+function push!(problem::Problem, elements...)
+    push!(problem.elements, elements...)
+end
+
+function push!(problem::Problem, elements::Vector)
+    push!(problem.elements, elements...)
 end
 
 function get_gdofs(element::Element, dim::Int)
diff --git a/src/solvers.jl b/src/solvers.jl
index 5acd55f..4d412be 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -511,13 +511,18 @@ function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factoriz
 end
 
 """ Convenience function to call linear solver. """
-function LinearSolver(problems...)
+function LinearSolver(problems::Problem...)
     solver = Solver(Linear, "default linear solver")
     if length(problems) != 0
         push!(solver, problems...)
     end
     return solver
 end
+function LinearSolver(name::ASCIIString, problems::Problem...)
+    solver = LinearSolver(problems...)
+    solver.name = name
+    return solver
+end
 
 ### End of linear quasistatic solver
 
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 170597a..8c018c4 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -100,3 +100,11 @@ end
     @test isapprox(fb, 1.0)
 end
 
+@testset "add two time dependent fields to element at once" begin
+    el = Element(Seg2, [1, 2])
+    update!(el, "foo1", 1.0 => 1.0)
+    update!(el, "foo1", 2.0 => 2.0)
+    update!(el, "foo2", 1.0 => 1.0, 2.0 => 2.0)
+    @test isapprox(el("foo1", 1.5), el("foo2", 1.5))
+end
+
diff --git a/test/test_heat.jl b/test/test_heat.jl
index bbfd036..2f2967b 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -3,7 +3,7 @@
 
 using JuliaFEM
 using JuliaFEM.Test
-
+using JuliaFEM.Preprocess
 
 @testset "test one element heat problem" begin
 
@@ -19,7 +19,6 @@ using JuliaFEM.Test
     update!(el1, "geometry", X)
     update!(el1, "temperature thermal conductivity", 6.0)
     update!(el1, "temperature load", 12.0)
-    update!(el1, "density", 36.0)
 
     # define boundary element for flux
     el2 = Element(Seg2, [1, 2])
@@ -29,6 +28,7 @@ using JuliaFEM.Test
 
     # define heat problem and push elements to problem
     problem = Problem(Heat, "one element heat problem", 1)
+    problem.properties.formulation = "2D"
     push!(problem, el1, el2)
 
     # define boundary element for dirichlet boundary condition
@@ -53,8 +53,7 @@ using JuliaFEM.Test
     @test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])
 
     # using Solver
-    solver = Solver("solve heat problem")
-    solver.is_linear_system = true
+    solver = LinearSolver("solve heat problem")
     push!(solver, problem, boundary_condition)
 
     # Set constant source f=12 with k=6. Accurate solution is
@@ -75,3 +74,87 @@ using JuliaFEM.Test
     @test isapprox(T[1], 2.0)
 end
 
+function T_acc(x)
+    # accurate solution
+    a = 0.01
+    L = 0.20
+    k = 50.0
+    Tᵤ = 20.0
+    h = 10.0
+    P = 4*a
+    A = a^2
+    α = h
+    β = sqrt((h*P)/(k*A))
+    T̂ = 100.0
+    C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T̂-Tᵤ, 0.0]
+    return dot(C, [exp(β*x), exp(-β*x)]) + Tᵤ
+end
+
+#=
+@testset "test 1d heat problem" begin
+    X = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0, 0.0],
+        2 => [0.1, 0.0, 0.0],
+        3 => [0.2, 0.0, 0.0])
+    e1 = Element(Seg2, [1, 2])
+    e2 = Element(Seg2, [2, 3])
+    e3 = Element(Poi1, [3])
+
+    p1 = Problem(Heat, "1d heat problem", 1)
+    p1.properties.formulation = "1D"
+    push!(p1, e1, e2, e3)
+    update!(p1, "geometry", X)
+    a = 0.010
+    update!(p1, "cross-section area", a^2)
+    update!(p1, "cross-section perimeter", 4*a)
+    update!(p1, "temperature thermal conductivity", 50.0) # k [W/(m∘C)]
+    update!(p1, "temperature heat transfer coefficient", 10.0) # h [W/(m²∘C)]
+    update!(p1, "temperature external temperature", 20.0)
+
+    p2 = Problem(Dirichlet, "left boundary", 1, "temperature")
+    e3 = Element(Poi1, [1])
+    update!(e3, "geometry", X)
+    update!(e3, "temperature 1", 100.0)
+    push!(p2, e3)
+
+    solver = LinearSolver(p1, p2)
+    call(solver)
+    T_min = minimum(p1.assembly.u)
+    @test isapprox(T_max, T_acc(0.2); rtol=4.5e-2)
+end
+=#
+
+@testset "test 3d heat problem" begin
+    fn = Pkg.dir("JuliaFEM") * "/test/testdata/rod_short.med"
+    mesh = aster_read_mesh(fn, "SHORT_ROD_RECTANGLE_HEX8")
+
+    p1 = Problem(Heat, "rod", 1)
+    push!(p1, create_elements(mesh, "ROD"))
+    push!(p1, create_elements(mesh, "SIDES"))
+    push!(p1, create_elements(mesh, "RIGHT"))
+    update!(p1, "temperature thermal conductivity", 50.0)
+    update!(p1, "temperature external temperature", 20.0)
+    update!(p1, "temperature heat transfer coefficient", 10.0)
+
+    p2 = Problem(Dirichlet, "left support T=100", 1, "temperature")
+    push!(p2, create_elements(mesh, "LEFT"))
+    update!(p2, "temperature 1", 100.0)
+
+    solver = LinearSolver(p1, p2)
+    call(solver)
+
+    T_min = minimum(p1.assembly.u)
+
+    # Code Aster solution
+    T_CA_HEX20 = 4.58158267950429E+01
+    T_CA_HEX8 = 3.77215189873436E+01
+    info("T_min = $T_min")
+    info("T_acc = $(T_acc(0.2))")
+    rtol1 = norm(T_min-T_CA_HEX8)/max(T_min,T_CA_HEX8)*100.0
+    rtol2 = norm(T_min-T_acc(0.2))/max(T_min,T_acc(0.2))*100.0
+    info("rel. tol to CA solution: $rtol1 %")
+    info("rel. tol to accurate solution: $rtol2 %")
+
+    @test isapprox(T_min, T_acc(0.2); rtol=18.0e-2)
+    @test isapprox(T_min, T_CA_HEX8; rtol=1.0e-9)
+end
diff --git a/test/test_problem.jl b/test/test_problem.jl
new file mode 100644
index 0000000..5192886
--- /dev/null
+++ b/test/test_problem.jl
@@ -0,0 +1,28 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "test initialize field problem" begin
+    el = Element(Seg2, [1, 2])
+    pr = Problem(Heat, 1)
+    push!(pr, el)
+    initialize!(pr)
+    @test haskey(el, "temperature")
+    # one timestep in field "temperature"
+    @test length(el("temperature")) == 1
+    @test length(el["temperature"]) == 1
+    # length of single increment
+    @test length(el("temperature", 0.0)) == 2
+    @test length(last(el, "temperature").data) == 2
+end
+
+@testset "test initialize boundary problem" begin
+    el = Element(Seg2, [1, 2])
+    pr = Problem(Dirichlet, "bc", 1, "temperature")
+    push!(pr, el)
+    initialize!(pr)
+    @test haskey(el, "reaction force")
+    @test haskey(el, "temperature")
+end