diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 2c5484f..6c0d502 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -27,33 +27,17 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis
 include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
+include("lagrange.jl") # Continuous Galerkin (Lagrange) elements
+export get_reference_coordinates
+export Poi1,
+       Seg2, Seg3,
+       Tri3, Tri6, Quad4, Quad8, Quad9,
+       Tet4, Tet10, Hex8, Hex20, Hex27
 
-type Poi1 <: AbstractElement
-end
-
-function size(element::Element{Poi1})
-    return (0, 1)
-end
-
-function length(element::Element{Poi1})
-    return 1
-end
-
-function get_basis(element::Element{Poi1}, ip, time)
-    return [1]
-end
-
-function call(element::Element{Poi1}, ip, time, ::Type{Val{:detJ}})
-    return 1.0
-end
-
-export Poi1, Seg2, Seg3, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
 include("nurbs.jl")
 export NSeg, NSurf, NSolid, is_nurbs
 
 #include("hierarchical.jl") # P-elements
-#include("mortar_elements.jl") # Mortar elements
-#include("equations.jl")
 
 include("integrate.jl")  # default integration points for elements
 export get_integration_points
@@ -66,7 +50,7 @@ export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
        get_unknown_field_dimension, get_gdofs, Assembly,
        get_parent_field_name, get_elements
 
-include("elasticity.jl") # elasticity equations
+include("elasticity.jl")
 export Elasticity
 
 include("dirichlet.jl")
@@ -75,7 +59,7 @@ export Dirichlet
 include("heat.jl")
 export Heat
 
-export assemble, assemble!
+export assemble!, postprocess!
 
 function assemble!(problem::Problem, element::Element, time=0.0)
     assemble!(problem.assembly, problem, element, time)
@@ -89,7 +73,7 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
        get_unknown_field_name, get_formulation_type,
        get_field_problems, get_boundary_problems,
        get_field_assembly, get_boundary_assembly,
-       initialize!, create_projection
+       initialize!, create_projection, eliminate_interior_dofs
 include("modal.jl")
 export Modal
 
@@ -124,9 +108,11 @@ export create_elements, Mesh,
        add_element!, add_elements!,
        add_element_to_element_set!,
        add_node_to_node_set!,
-       find_nearest_nodes
+       find_nearest_nodes,
+       reorder_element_connectivity!
 include("preprocess_abaqus_reader.jl")
 include("preprocess_abaqus_reader_old.jl")
+export parse_abaqus, parse_section, parse_element_section
 include("preprocess_aster_reader.jl")
 export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
        parse_aster_header, aster_parse_nodes, aster_renumber_nodes!,
@@ -146,10 +132,14 @@ export get_mesh, get_model
 
 module Postprocess
 include("postprocess_utils.jl")
-export calc_nodal_values!, get_nodal_vector, copy_field!
+export calc_nodal_values!, 
+       get_nodal_vector,
+       get_nodal_dict,
+       copy_field!
 include("postprocess_xdmf.jl")
 export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
 end
+export Postprocessor
 
 """ JuliaFEM testing routines. """
 module Test
@@ -171,5 +161,4 @@ module Interfaces
 include("interfaces.jl")
 end
 
-
 end # module
diff --git a/src/assembly.jl b/src/assembly.jl
index 79a6654..afc827a 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -1,18 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# Functions to handle global assembly of problem
-
-type CAssembly
-    interior_dofs :: Vector{Int}
-    boundary_dofs :: Vector{Int}
-    F :: Union{Factorization, Matrix}
-    Kc :: SparseMatrixCSC
-    fc :: SparseMatrixCSC
-    Kib :: SparseMatrixCSC
-    fi :: SparseMatrixCSC
-end
-
 function optimize!(assembly::Assembly)
     optimize!(assembly.K)
     optimize!(assembly.Kg)
@@ -99,54 +87,28 @@ function assemble!(problem::Problem, time::Real, ::Type{Val{:mass_matrix}}; dens
     end
 end
 
-""" Calculate reduced stiffness matrix.
-mindofs: if dofs < mindofs, do not reduce
-"""
-function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}, mindofs=100000)
-    all_dofs = unique(assembly.stiffness_matrix.I)
-    boundary_dofs = intersect(all_dofs, boundary_dofs_)
-    interior_dofs = setdiff(all_dofs, boundary_dofs_)
+# Static condensation routines
 
-    K = sparse(assembly.stiffness_matrix)
-    f = sparse(assembly.force_vector)
+function eliminate_interior_dofs(K::SparseMatrixCSC, f::SparseMatrixCSC, B::Vector{Int64}, I::Vector{Int64}; F=nothing, chunk_size=100000)
 
     dim = size(K, 1)
+    Kib = K[I,B]
 
-    # empty assembly to release memory for factorization
-    empty!(assembly.stiffness_matrix)
-    empty!(assembly.force_vector)
-    gc()
-
-    if dim < mindofs
-        # no need to do any reduction of matrix size at all, just \ it.
-        return CAssembly([], all_dofs, Matrix{Float64}(), K, f, spzeros(0, 0), spzeros(0,1))
+    if F == nothing
+        F = cholfact(1/2*(K + K')[I,I])
     end
 
-    # check that matrix is symmetric
-    s = maximum(abs(1/2*(K + K') - K))
-    @assert s < 1.0e-6
-    K = 1/2*(K + K')
-
-    Kib = K[interior_dofs, boundary_dofs]
-    Kbb = K[boundary_dofs, boundary_dofs]
-    fi = f[interior_dofs]
-    fb = f[boundary_dofs]
-
-    F = cholfact(K[interior_dofs, interior_dofs])
-    K = 0
-    gc()
-
-
-    if dim < 100000
+    if dim < chunk_size
         # for small problems we don't need to care about memory usage
         Kd = Kib' * (F \ Kib)
     else
         # for larger problems calculate schur complement in pieces
-        nb = length(boundary_dofs)
+        nb = length(B)
         p = nb > 10 ? round(Int, nb/10) : nb
         Kd = zeros(nb, nb)
         for bi in 1:nb
-            mod(bi, p) == 0 && info("Reduction: ", round(Int, bi/nb*100), " % done")
+            done = round(Int, bi/nb*100)
+            mod(bi, p) == 0 && info("Static condensation: $done % done")
             C = full(F \ Kib[:, bi])
             for bj in bi:nb
                 d = Kib[:, bj]
@@ -155,38 +117,17 @@ function reduce(assembly::Assembly, boundary_dofs_::Vector{Int}, mindofs=100000)
         end
         Kd += tril(Kd, -1)'
     end
-    Kc = spzeros(dim, dim)
-    Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
-
-#=  # this is slightly faster but uses more memory
-    chunks = round(Int, dim/3000)
-    info("Reduction is done in $chunks chunks.")
-    nb = length(boundary_dofs)
-    kk = round(Int, collect(linspace(0, nb, chunks+1)))
-    sl = [kk[j]+1:kk[j+1] for j=1:length(kk)-1]
-    Kd = zeros(Float64, nb, nb)
-    #Kd = SharedArray(Float64, nb, nb)
-    for (k,sli) in enumerate(sl)
-        b1 = boundary_dofs[sli]
-        Sc = F \ Kib[:,sli]
-        for slj in sl
-            b2 = boundary_dofs[slj]
-            #Kc[b2,b1] = Kbb[slj,sli] - Kib[:,slj]'*Sc
-            Kd[slj, sli] = Kib[:,slj]'*Sc
-        end
-        info("Reduction: ", round(k/chunks*100, 0), " % done")
-    end
 
     Kc = spzeros(dim, dim)
-    Kc[boundary_dofs, boundary_dofs] = Kbb - Kd
-=#
+    Kc[B,B] = K[B,B] - Kd
 
     fc = spzeros(dim, 1)
-    fc[boundary_dofs] = fb - Kib' * (F \ fi)
+    fc[B] = f[B] - Kib' * (F \ f[I])
 
-    return CAssembly(interior_dofs, boundary_dofs, F, Kc, fc, Kib, fi)
+    return Kc, fc
 end
 
+#=
 function reconstruct!(ca::CAssembly, x::SparseMatrixCSC)
     if isa(ca.F, Factorization)
         x[ca.interior_dofs] = ca.F \ (ca.fi - ca.Kib*x[ca.boundary_dofs])
@@ -194,10 +135,4 @@ function reconstruct!(ca::CAssembly, x::SparseMatrixCSC)
         x[ca.interior_dofs] = ca.F * (ca.fi - ca.Kib*x[ca.boundary_dofs])
     end
 end
-
-function Base.(:+)(ass1::Assembly, ass2::Assembly)
-    mass_matrix = ass1.mass_matrix + ass2.mass_matrix
-    stiffness_matrix = ass1.stiffness_matrix + ass2.stiffness_matrix
-    force_vector = ass1.force_vector + ass2.force_vector
-    return Assembly(mass_matrix, stiffness_matrix, force_vector)
-end
+=#
diff --git a/src/contact.jl b/src/contact.jl
index 10b0185..b6fd476 100644
--- a/src/contact.jl
+++ b/src/contact.jl
@@ -96,11 +96,11 @@ function assemble!(problem::Problem{Contact}, time::Float64,
         end
 
         # 3. loop all master elements
-        for master_element in slave_element["master elements"](time)
+        for master_element in slave_element("master elements", time)
 
             nm = length(master_element)
-            X2 = master_element["geometry"](time)
-            u2 = master_element["displacement"](time)
+            X2 = master_element("geometry", time)
+            u2 = master_element("displacement", time)
             x2 = X2 + u2
 
             norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index cd79b25..83cafe5 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -36,7 +36,7 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     else
         Ae = eye(nnodes)
         De = zeros(nnodes, nnodes)
-        for ip in get_integration_points(element)
+        for ip in get_integration_points(element, 1)
             N = element(ip, time)
             detJ = element(ip, time, Val{:detJ})
             De += ip.weight*N'*N*detJ
@@ -53,7 +53,7 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
     end
 
     # right hand side
-    for ip in get_integration_points(element)
+    for ip in get_integration_points(element, 1)
         detJ = element(ip, time, Val{:detJ})
         w = ip.weight*detJ
         N = element(ip, time)
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 81f10ba..bc7a1cb 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -1,8 +1,36 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-""" Elasticity problem
+""" Elasticity equations.
 
+Field equation is:
+    
+    m∂²u/∂t² = ∇⋅σ - b
+
+Weak form is: find u∈U such that ∀v in V
+
+    δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0
+
+where
+
+    ρ₀ = density
+    b₀ = displacement load
+    t₀ = displacement traction
+
+Formulations
+------------
+plane stress, plane strain, 3D
+
+References
+----------
+
+https://en.wikipedia.org/wiki/Linear_elasticity
+https://en.wikipedia.org/wiki/Finite_strain_theory
+https://en.wikipedia.org/wiki/Stress_measures
+https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
+https://en.wikipedia.org/wiki/Strain_energy_density_function
+https://en.wikipedia.org/wiki/Plane_stress
+https://en.wikipedia.org/wiki/Hooke's_law
 
 """
 type Elasticity <: FieldProblem
@@ -38,8 +66,14 @@ function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::El
     add!(assembly.f, gdofs, f)
 end
 
+typealias Elasticity2DSurfaceElements Union{Poi1, Seg2, Seg3}
+typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9}
+typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9}
+typealias Elasticity3DVolumeElements Union{Tet4, Tet10, Hex8, Hex20, Hex27}
+
+
 """ Elasticity equations for 2d cases. """
-function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
+function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time, ::Type{Val{:plane}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -162,7 +196,7 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
     return Km, Kg, f
 end
 
-function assemble{El<:Union{Poi1,Seg2,Seg3}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
+function assemble{El<:Elasticity2DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:plane}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -203,7 +237,7 @@ function assemble{El<:Union{Poi1,Seg2,Seg3}}(problem::Problem{Elasticity}, eleme
 end
 
 """ Elasticity equations, 3d, linear. """
-function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
+function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -270,7 +304,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
 end
 
 """ Material and geometric stiffness for linear buckling analysis. """
-function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_buckling}})
+function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_buckling}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -348,7 +382,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
 end
 
 """ Elasticity equations, 3d nonlinear. """
-function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+function assemble{El<:Elasticity3DVolumeElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -487,7 +521,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
 end
 
 """ Elasticity equations, surface traction for continuum formulation. """
-function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
+function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum}})
 
     props = problem.properties
     dim = get_unknown_field_dimension(problem)
@@ -521,39 +555,11 @@ function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, el
     return Km, Kg, f
 end
 
-function assemble{El<:Union{Tri3, Tri6, Quad4}}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
+function assemble{El<:Elasticity3DSurfaceElements}(problem::Problem{Elasticity}, element::Element{El}, time::Real, ::Type{Val{:continuum_linear}})
     return assemble(problem, element, time, Val{:continuum})
 end
 
 """ Elasticity equations using ForwardDiff
-
-Formulation
------------
-
-Field equation is:
-∂u/∂t = ∇⋅f - b
-
-Weak form is: find u∈U such that ∀v in V
-
-    δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0
-
-where
-
-    ρ₀ = density
-    b₀ = displacement load
-    t₀ = displacement traction
-
-References
-----------
-
-https://en.wikipedia.org/wiki/Linear_elasticity
-https://en.wikipedia.org/wiki/Finite_strain_theory
-https://en.wikipedia.org/wiki/Stress_measures
-https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
-https://en.wikipedia.org/wiki/Strain_energy_density_function
-https://en.wikipedia.org/wiki/Plane_stress
-https://en.wikipedia.org/wiki/Hooke's_law
-
 """
 function assemble(problem::Problem{Elasticity}, element::Element, time::Real, ::Type{Val{:forwarddiff}})
 
diff --git a/src/elements.jl b/src/elements.jl
index 9cd92f6..1538777 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -41,7 +41,7 @@ function last(element::Element, field_name::ASCIIString)
     return last(element[field_name])
 end
 
-function call(element::Element, ip, time)
+function call(element::Element, ip, time=0.0)
     return get_basis(element, ip, time)
 end
 
@@ -75,8 +75,17 @@ function call(element::Element, field_name::ASCIIString, ip, time, ::Type{Val{:G
     return element(ip, time, Val{:Grad})*element[field_name](time)
 end
 
+function call(element::Element, field::Field, time)
+    return field(time)
+end
+
+function call(element::Element, field::DCTI, time)
+    return field.data
+end
+
 function call(element::Element, field_name::ASCIIString, time)
-    return element[field_name](time)
+    field = element[field_name]
+    return call(element, field, time)
 end
 
 function call(element::Element, field_name::ASCIIString, ip, time::Float64)
@@ -220,6 +229,7 @@ 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)
@@ -227,6 +237,7 @@ function get_dbasis(element::Element, ip, time)
     basis(xi) = vec(get_basis(element, xi, time))
     return ForwardDiff.jacobian(basis, xi)'
 end
+=#
 
 """ Check existence of field. """
 function haskey(element::Element, field_name)
diff --git a/src/fields.jl b/src/fields.jl
index 6839f26..a73d113 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -142,57 +142,55 @@ end
 
 ### Accessing and manipulating discrete fields
 
-function Base.getindex(field::DVTV, i::Int64)
+function getindex(field::DVTV, i::Int64)
     return field.data[i]
 end
 
-function Base.push!(field::DCTV, data::Pair)
+function 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")
+function push!(field::DVTV, data::Pair)
     push!(field.data, data)
 end
 
-function Base.getindex(field::DVTV, i::Int64)
+function getindex(field::DVTV, i::Int64)
     return field.data[i]
 end
 
-function Base.getindex(field::DVTI, i::Int64)
+function getindex(field::DVTI, i::Int64)
     return field.data[i]
 end
 
-function Base.getindex(field::DCTV, i::Int64)
+function getindex(field::DCTV, i::Int64)
     return field.data[i]
 end
 
-function Base.getindex(field::Field, i::Int64)
+function getindex(field::Field, i::Int64)
     return field.data[i]
 end
 
-function Base.length(field::DVTI)
+function length(field::DVTI)
     return length(field.data)
 end
 
-function Base.length(field::DCTI)
+function length(field::DCTI)
     return 1
 end
 
-function Base.length(field::DVTV)
+function length(field::DVTV)
     return length(field.data)
 end
 
-function Base.length(field::DCTV)
+function length(field::DCTV)
     return length(field.data)
 end
 
-function Base.first(field::Union{DCTV, DVTV})
+function first(field::Union{DCTV, DVTV})
     return field[1]
 end
 
-function Base.isapprox(f1::DCTI, f2::DCTI)
+function isapprox(f1::DCTI, f2::DCTI)
     isapprox(f1.data, f2.data)
 end
 
@@ -236,8 +234,7 @@ function vec(field::DVTI)
 end
 
 function vec(field::DCTV)
-    info("trying to vectorize $field")
-    error("does not make sense")
+    error("trying to vectorize $field does not make sense")
 end
 
 function endof(field::Field)
@@ -321,21 +318,21 @@ end
 ### Interpolation
 
 """ Interpolate time-invariant field in time direction. """
-function Base.call(field::DVTI, time::Float64)
+function call(field::DVTI, time::Float64)
     return field
 end
-function Base.call(field::DCTI, time::Float64)
+function call(field::DCTI, time::Float64)
     return field
 end
-function Base.call(field::CVTI, time::Float64)
+function call(field::CVTI, time::Float64)
     return field.data()
 end
-function Base.call(field::CCTI, time::Float64)
+function call(field::CCTI, time::Float64)
     return field.data()
 end
 
 """ Interpolate constant time-variant field in time direction. """
-function Base.call(field::DCTV, time::Real)
+function call(field::DCTV, time::Real)
     time < first(field).time && return DCTI(first(field).data)
     time > last(field).time && return DCTI(last(field).data)
     for i=reverse(1:length(field))
@@ -355,7 +352,7 @@ function Base.call(field::DCTV, time::Real)
     error("interpolate DCTV: unknown failure when interpolating $(field.data) for time $time")
 end
 
-function Base.call(field::DVTV, time::Float64)
+function call(field::DVTV, time::Float64)
     time < first(field).time && return DVTI(first(field).data)
     time > last(field).time && return DVTI(last(field).data)
     for i=reverse(1:length(field))
@@ -376,17 +373,17 @@ function Base.call(field::DVTV, time::Float64)
 end
 
 """ Interpolate constant field in spatial dimension. """
-function Base.call(basis::CVTI, field::DCTI, xi::Vector)
+function call(basis::CVTI, field::DCTI, xi::Vector)
     return field.data
 end
 
 """ Interpolate variable field in spatial dimension. """
-function Base.call(basis::CVTI, values::DVTI, xi::Vector)
+function 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}})
+function 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)])
     J = sum([kron(dbasis[:,i], geometry[i]') for i=1:length(geometry)])
@@ -395,14 +392,14 @@ function Base.call(basis::CVTI, geometry::DVTI, xi::Vector, ::Type{Val{:grad}})
     return grad
 end
 
-function Base.call(basis::CVTI, geometry::DVTI, values::DVTI, xi::Vector, ::Type{Val{:grad}})
+function 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)])'
     gradf = sum([kron(grad[:,i], values[i]') for i=1:length(values)])'
     return length(gradf) == 1 ? gradf[1] : gradf
 end
 
-function Base.call(basis::CVTI, xi::Vector, time::Number)
+function call(basis::CVTI, xi::Vector, time::Number)
     call(basis, xi)
 end
 
@@ -413,3 +410,4 @@ end
 ### FIELDSET ###
 
 typealias FieldSet Dict{ASCIIString, Field}
+
diff --git a/src/heat.jl b/src/heat.jl
index 647cd4f..fb64321 100644
--- a/src/heat.jl
+++ b/src/heat.jl
@@ -1,13 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# Heat problems
-
 """ Heat equations.
 
-Formulation
------------
-
 Field equation is:
 
     ρc∂u/∂t = ∇⋅(k∇u) + f
@@ -27,7 +22,6 @@ Parameters
 temperature thermal conductivity
 temperature load
 temperature flux
-
 thermal conductivity
 heat source
 heat flux
@@ -43,6 +37,7 @@ References
 ----------
 https://en.wikipedia.org/wiki/Heat_equation
 https://en.wikipedia.org/wiki/Heat_capacity
+https://en.wikipedia.org/wiki/Heat_flux
 https://en.wikipedia.org/wiki/Thermal_conduction
 https://en.wikipedia.org/wiki/Thermal_conductivity
 https://en.wikipedia.org/wiki/Thermal_diffusivity
@@ -50,17 +45,18 @@ https://en.wikipedia.org/wiki/Volumetric_heat_capacity
 """
 type Heat <: FieldProblem
     formulation :: ASCIIString
+    store_fields :: Vector{ASCIIString}
 end
 
 function Heat()
-    return Heat("3D")
+    return Heat("3D", [])
 end
 
 function get_unknown_field_name(problem::Problem{Heat})
     return "temperature"
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
+function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time)
     formulation = Val{Symbol(problem.properties.formulation)}
     assemble!(assembly, problem, element, time, formulation)
 end
@@ -71,8 +67,8 @@ function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Eleme
     info("Unknown element type $E for 3d heat problem!")
 end
 
-typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8}
-typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4}
+typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8, Hex20, Hex27}
+typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4, Quad8, Quad9}
 
 function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
     gdofs = get_gdofs(problem, element)
@@ -100,13 +96,44 @@ function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem
     add!(assembly.f, gdofs, fq)
 end
 
+function postprocess!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time)
+    haskey(element, "temperature thermal conductivity") || return
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    Me = zeros(nnodes, nnodes)
+    De = zeros(nnodes, nnodes)
+    f = zeros(nnodes, 3)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        De += w*diagm(vec(N))
+        Me += w*N'*N
+    end
+    Ae = De*inv(Me)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = vec(element(ip, time))
+        Phi = transpose(Ae*N)
+        k = element("temperature thermal conductivity", ip, time)
+        gradT = element(field_name, ip, time, Val{:Grad})
+        q = -vec(k*gradT)
+        update!(ip, "heat flux", time => q)
+        f += w*Phi'*q'
+    end
+    add!(assembly.M, gdofs, gdofs, De)
+    add!(assembly.f, gdofs, [1, 2, 3], f)
+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)
+    for ip in get_integration_points(element, 1)
         detJ = element(ip, time, Val{:detJ})
         w = ip.weight*detJ
         N = element(ip, time)
diff --git a/src/integrate.jl b/src/integrate.jl
index 2aa8277..63f892e 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -67,8 +67,8 @@ end
 ### 1d elements
 
 typealias CartesianLineElement Union{Seg2, Seg3, NSeg}
-typealias CartesianSurfaceElement Union{Quad4, NSurf}
-typealias CartesianVolumeElement Union{Hex8, NSolid}
+typealias CartesianSurfaceElement Union{Quad4, Quad8, Quad9, NSurf}
+typealias CartesianVolumeElement Union{Hex8, Hex20, Hex27, NSolid}
 
 function get_integration_points(element::CartesianLineElement, order::Int64)
     w, xi = get_integration_points(order)
@@ -238,7 +238,7 @@ end
 
 typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Hex8}
 
-typealias QuadraticElement Union{Seg3, Tri6, Tet10}
+typealias QuadraticElement Union{Seg3, Tri6, Tet10, Quad8, Quad9, Hex20, Hex27}
 
 function get_integration_order(element::LinearElement)
     return 2
diff --git a/src/lagrange.jl b/src/lagrange.jl
new file mode 100644
index 0000000..36026cb
--- /dev/null
+++ b/src/lagrange.jl
@@ -0,0 +1,570 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+### 0d element
+
+type Poi1 <: AbstractElement
+end
+
+function description(::Type{Poi1})
+    "1 node point"
+end
+
+function size(element::Element{Poi1})
+    return (0, 1)
+end
+
+function length(element::Element{Poi1})
+    return 1
+end
+
+function get_basis(element::Element{Poi1}, ip, time)
+    return [1]
+end
+
+function call(element::Element{Poi1}, ip, time, ::Type{Val{:detJ}})
+    return 1.0
+end
+
+### 1d elements
+
+type Seg2 <: AbstractElement
+end
+
+function description(::Type{Seg2})
+    "2 node segment"
+end
+
+function size(element::Element{Seg2})
+    return (1, 2)
+end
+
+function length(element::Element{Seg2})
+    return 2
+end
+
+function get_reference_coordinates(::Type{Seg2})
+    Vector{Float64}[
+        [-1.0], # N1
+        [ 1.0]] # N2
+end
+
+function get_interpolation_polynomial(::Type{Seg2}, xi)
+    [1.0 xi[1]]
+end
+
+function get_interpolation_polynomial(::Type{Seg2}, xi, ::Type{Val{:partial_derivatives}})
+    [0.0 1.0]
+end
+
+#
+
+type Seg3 <: AbstractElement
+end
+
+function description(::Type{Seg3})
+    "3 node segment"
+end
+
+function size(element::Element{Seg3})
+    return (1, 3)
+end
+
+function length(element::Element{Seg3})
+    return 3
+end
+
+function get_reference_coordinates(::Type{Seg3})
+    Vector{Float64}[
+        [-1.0], # N1
+        [ 1.0], # N2
+        [ 0.0]] # N3
+end
+
+function get_interpolation_polynomial(::Type{Seg3}, xi)
+    [1.0 xi[1] xi[1]^2]
+end
+
+function get_interpolation_polynomial(::Type{Seg3}, xi, ::Type{Val{:partial_derivatives}})
+    [0.0 1.0 2.0*xi[1]]
+end
+
+### 2d elements
+
+type Tri3 <: AbstractElement
+end
+
+function description(::Type{Tri3})
+    "3 node triangle"
+end
+
+function size(element::Element{Tri3})
+    return (2, 3)
+end
+
+function length(element::Element{Tri3})
+    return 3
+end
+
+function get_reference_coordinates(::Type{Tri3})
+    Vector{Float64}[
+        [0.0, 0.0], # N1
+        [1.0, 0.0], # N2
+        [0.0, 1.0]] # N3
+end
+
+function get_interpolation_polynomial(::Type{Tri3}, xi)
+    [
+        1 xi[1] xi[2]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Tri3}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0.0 1.0 0.0
+        0.0 0.0 1.0
+    ]
+end
+
+#
+
+type Tri6 <: AbstractElement
+end
+
+function description(::Type{Tri6})
+    "6 node triangle"
+end
+
+function size(element::Element{Tri6})
+    return (2, 6)
+end
+
+function length(element::Element{Tri6})
+    return 6
+end
+
+function get_reference_coordinates(::Type{Tri6})
+    Vector{Float64}[
+        [0.0, 0.0], # N1
+        [1.0, 0.0], # N2
+        [0.0, 1.0], # N3
+        [0.5, 0.0], # N4
+        [0.5, 0.5], # N5
+        [0.0, 0.5]] # N6
+end
+
+function get_interpolation_polynomial(::Type{Tri6}, xi)
+    [
+        1 xi[1] xi[2] xi[1]^2 xi[1]*xi[2] xi[2]^2
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Tri6}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  1  0  2*xi[1]  xi[2]  0
+        0  0  1  0        xi[1]  2*xi[2]
+    ]
+end
+
+#
+
+type Quad4 <: AbstractElement
+end
+
+function description(::Type{Quad4})
+    "4 node quadrangle"
+end
+
+function size(element::Element{Quad4})
+    return (2, 4)
+end
+
+function length(element::Element{Quad4})
+    return 4
+end
+
+function get_reference_coordinates(::Type{Quad4})
+    Vector{Float64}[
+        [-1.0, -1.0], # N1
+        [ 1.0, -1.0], # N2
+        [ 1.0,  1.0], # N3
+        [-1.0,  1.0]] # N4
+end
+
+function get_interpolation_polynomial(::Type{Quad4}, xi)
+    [
+        1.0 xi[1] xi[2] xi[1]*xi[2]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Quad4}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0 1 0 xi[2]
+        0 0 1 xi[1]      
+    ]
+end
+
+#
+
+type Quad8 <: AbstractElement
+end
+
+function description(::Type{Quad8})
+    "8 node Serendip quadrangle"
+end
+
+function size(element::Element{Quad8})
+    return (2, 8)
+end
+
+function length(element::Element{Quad8})
+    return 8
+end
+
+function get_reference_coordinates(::Type{Quad8})
+    Vector{Float64}[
+        [-1.0, -1.0], # N1
+        [ 1.0, -1.0], # N2
+        [ 1.0,  1.0], # N3
+        [-1.0,  1.0], # N4
+        [ 0.0, -1.0], # N5
+        [ 1.0,  0.0], # N6
+        [ 0.0,  1.0], # N7
+        [-1.0,  0.0]] # N8
+end
+
+function get_interpolation_polynomial(::Type{Quad8}, xi)
+    [
+        1 xi[2] xi[1] xi[2]^2 xi[1]*xi[2] xi[1]^2 xi[1]*xi[2]^2 xi[1]^2*xi[2]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Quad8}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  1  0        xi[2]  2*xi[1]  xi[2]^2        2*xi[1]*xi[2]
+        0  1  0  2*xi[2]  xi[1]  0        2*xi[1]*xi[2]  xi[1]^2
+    ]
+end
+
+#
+
+type Quad9 <: AbstractElement
+end
+
+function description(::Type{Quad9})
+    "9 node quadrangle"
+end
+
+function size(element::Element{Quad9})
+    return (2, 9)
+end
+
+function length(element::Element{Quad9})
+    return 9
+end
+
+function get_reference_coordinates(::Type{Quad9})
+    Vector{Float64}[
+        [-1.0, -1.0], # N1
+        [ 1.0, -1.0], # N2
+        [ 1.0,  1.0], # N3
+        [-1.0,  1.0], # N4
+        [ 0.0, -1.0], # N5
+        [ 1.0,  0.0], # N6
+        [ 0.0,  1.0], # N7
+        [-1.0,  0.0], # N8
+        [ 0.0,  0.0]] # N9
+end
+
+function get_interpolation_polynomial(::Type{Quad9}, xi)
+    [
+        1 xi[2] xi[1] xi[2]^2 xi[1]*xi[2] xi[1]^2 xi[1]*xi[2]^2 xi[1]^2*xi[2] xi[1]^2*xi[2]^2
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Quad9}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  1  0        xi[2]  2*xi[1]  xi[2]^2        2*xi[1]*xi[2]  2*xi[1]*xi[2]^2
+        0  1  0  2*xi[2]  xi[1]  0        2*xi[1]*xi[2]  xi[1]^2        2*xi[1]^2*xi[2]
+    ]
+end
+
+### 3d elements
+
+type Tet4 <: AbstractElement
+end
+
+function description(::Type{Tet4})
+    "4 node tetrahedral element"
+end
+
+function size(element::Element{Tet4})
+    return (3, 4)
+end
+
+function length(element::Element{Tet4})
+    return 4
+end
+
+function get_reference_coordinates(::Type{Tet4})
+    Vector{Float64}[
+        [0.0, 0.0, 0.0], # N1
+        [1.0, 0.0, 0.0], # N2
+        [0.0, 1.0, 0.0], # N3
+        [0.0, 0.0, 1.0]] # N4
+end
+
+function get_interpolation_polynomial(::Type{Tet4}, xi)
+    [
+        1.0 xi[1] xi[2] xi[3]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Tet4}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0.0 1.0 0.0 0.0
+        0.0 0.0 1.0 0.0
+        0.0 0.0 0.0 1.0
+    ]
+end
+
+#
+
+type Tet10 <: AbstractElement
+end
+
+function description(::Type{Tet10})
+    "10 node tetrahedral element"
+end
+
+function size(element::Element{Tet10})
+    return (3, 10)
+end
+
+function length(element::Element{Tet10})
+    return 10
+end
+
+function get_reference_coordinates(::Type{Tet10})
+    Vector{Float64}[
+        [0.0, 0.0, 0.0], # N1
+        [1.0, 0.0, 0.0], # N2
+        [0.0, 1.0, 0.0], # N3
+        [0.0, 0.0, 1.0], # N4
+        [0.5, 0.0, 0.0], # N5
+        [0.5, 0.5, 0.0], # N6
+        [0.0, 0.5, 0.0], # N7
+        [0.0, 0.0, 0.5], # N8
+        [0.5, 0.0, 0.5], # N9
+        [0.0, 0.5, 0.5]] # N10
+end
+
+function get_interpolation_polynomial(::Type{Tet10}, xi)
+    [
+        1.0 xi[3] xi[2] xi[1] xi[3]^2 xi[2]*xi[3] xi[2]^2 xi[1]*xi[3] xi[1]*xi[2] xi[1]^2
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Tet10}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  0  1  0        0      0        xi[3]  xi[2]  2*xi[1]
+        0  0  1  0  0        xi[3]  2*xi[2]  0      xi[1]  0      
+        0  1  0  0  2*xi[3]  xi[2]  0        xi[1]  0      0          
+    ]
+end
+
+#
+
+type Hex8 <: AbstractElement
+end
+
+function description(::Type{Hex8})
+    "8 node hexahedral element"
+end
+
+function size(element::Element{Hex8})
+    return (3, 8)
+end
+
+function length(element::Element{Hex8})
+    return 8
+end
+
+function get_reference_coordinates(::Type{Hex8})
+    Vector{Float64}[
+        [-1.0, -1.0, -1.0], # N1
+        [ 1.0, -1.0, -1.0], # N2
+        [ 1.0,  1.0, -1.0], # N3
+        [-1.0,  1.0, -1.0], # N4
+        [-1.0, -1.0,  1.0], # N5
+        [ 1.0, -1.0,  1.0], # N6
+        [ 1.0,  1.0,  1.0], # N7
+        [-1.0,  1.0,  1.0]] # N8
+end
+
+function get_interpolation_polynomial(::Type{Hex8}, xi)
+    [
+        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Hex8}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]
+        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]
+        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]
+    ]
+end
+
+#
+
+type Hex20 <: AbstractElement
+end
+
+function description(::Type{Hex20})
+    "20 node hexahedral element"
+end
+
+function size(element::Element{Hex20})
+    return (3, 20)
+end
+
+function length(element::Element{Hex20})
+    return 20
+end
+
+function get_reference_coordinates(::Type{Hex20})
+    Vector{Float64}[
+        [-1.0, -1.0, -1.0], # N1
+        [ 1.0, -1.0, -1.0], # N2
+        [ 1.0,  1.0, -1.0], # N3
+        [-1.0,  1.0, -1.0], # N4
+        [-1.0, -1.0,  1.0], # N5
+        [ 1.0, -1.0,  1.0], # N6
+        [ 1.0,  1.0,  1.0], # N7
+        [-1.0,  1.0,  1.0], # N8
+        [ 0.0, -1.0, -1.0], # N9
+        [ 1.0,  0.0, -1.0], # N10
+        [ 0.0,  1.0, -1.0], # N11
+        [-1.0,  0.0, -1.0], # N12
+        [-1.0, -1.0,  0.0], # N13
+        [ 1.0, -1.0,  0.0], # N14
+        [ 1.0,  1.0,  0.0], # N15
+        [-1.0,  1.0,  0.0], # N16
+        [ 0.0, -1.0,  1.0], # N17
+        [ 1.0,  0.0,  1.0], # N18
+        [ 0.0,  1.0,  1.0], # N19
+        [-1.0,  0.0,  1.0]] # N20
+end
+
+function get_interpolation_polynomial(::Type{Hex20}, xi)
+    [
+        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3] xi[3]^2 xi[2]^2 xi[1]^2 xi[2]*xi[3]^2 xi[2]^2*xi[3] xi[1]*xi[3]^2 xi[1]*xi[2]^2 xi[1]^2*xi[3] xi[1]^2*xi[2] xi[1]*xi[2]*xi[3]^2 xi[1]*xi[2]^2*xi[3] xi[1]^2*xi[2]*xi[3]
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Hex20}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]  0        0        2*xi[1]  0              0              xi[3]^2        xi[2]^2        2*xi[1]*xi[3]  2*xi[1]*xi[2]  xi[2]*xi[3]^2        xi[2]^2*xi[3]        2*xi[1]*xi[2]*xi[3]
+        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]  0        2*xi[2]  0        xi[3]^2        2*xi[2]*xi[3]  0              2*xi[1]*xi[2]  0              xi[1]^2        xi[1]*xi[3]^2        2*xi[1]*xi[2]*xi[3]  xi[1]^2*xi[3]
+        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]  2*xi[3]  0        0        2*xi[2]*xi[3]  xi[2]^2        2*xi[1]*xi[3]  0              xi[1]^2        0              2*xi[1]*xi[2]*xi[3]  xi[1]*xi[2]^2        xi[1]^2*xi[2]
+    ]
+end
+
+###
+
+type Hex27 <: AbstractElement
+end
+
+function description(::Type{Hex27})
+    "27 node hexahedral element"
+end
+
+function size(element::Element{Hex27})
+    return (3, 27)
+end
+
+function length(element::Element{Hex27})
+    return 27
+end
+
+function get_reference_coordinates(::Type{Hex27})
+    Vector{Float64}[
+        [-1.0, -1.0, -1.0], # N1
+        [ 1.0, -1.0, -1.0], # N2
+        [ 1.0,  1.0, -1.0], # N3
+        [-1.0,  1.0, -1.0], # N4
+        [-1.0, -1.0,  1.0], # N5
+        [ 1.0, -1.0,  1.0], # N6
+        [ 1.0,  1.0,  1.0], # N7
+        [-1.0,  1.0,  1.0], # N8
+        [ 0.0, -1.0, -1.0], # N9
+        [ 1.0,  0.0, -1.0], # N10
+        [ 0.0,  1.0, -1.0], # N11
+        [-1.0,  0.0, -1.0], # N12
+        [-1.0, -1.0,  0.0], # N13
+        [ 1.0, -1.0,  0.0], # N14
+        [ 1.0,  1.0,  0.0], # N15
+        [-1.0,  1.0,  0.0], # N16
+        [ 0.0, -1.0,  1.0], # N17
+        [ 1.0,  0.0,  1.0], # N18
+        [ 0.0,  1.0,  1.0], # N19
+        [-1.0,  0.0,  1.0], # N20
+        [ 0.0,  0.0, -1.0], # N21
+        [ 0.0, -1.0,  0.0], # N22
+        [ 1.0,  0.0,  0.0], # N23
+        [ 0.0,  1.0,  0.0], # N24
+        [-1.0,  0.0,  0.0], # N25
+        [ 0.0,  0.0,  1.0], # N26
+        [ 0.0,  0.0,  0.0]] # N27
+end
+
+function get_interpolation_polynomial(::Type{Hex27}, xi)
+    [
+        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3] xi[3]^2 xi[2]^2 xi[1]^2 xi[2]*xi[3]^2 xi[2]^2*xi[3] xi[1]*xi[3]^2 xi[1]*xi[2]^2 xi[1]^2*xi[3] xi[1]^2*xi[2] xi[2]^2*xi[3]^2 xi[1]*xi[2]*xi[3]^2 xi[1]*xi[2]^2*xi[3] xi[1]^2*xi[3]^2 xi[1]^2*xi[2]*xi[3] xi[1]^2*xi[2]^2 xi[1]*xi[2]^2*xi[3]^2 xi[1]^2*xi[2]*xi[3]^2 xi[1]^2*xi[2]^2*xi[3] xi[1]^2*xi[2]^2*xi[3]^2
+    ]
+end
+
+function get_interpolation_polynomial(::Type{Hex27}, xi, ::Type{Val{:partial_derivatives}})
+    [
+        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]  0        0        2*xi[1]  0              0              xi[3]^2        xi[2]^2        2*xi[1]*xi[3]  2*xi[1]*xi[2]  0                xi[2]*xi[3]^2        xi[2]^2*xi[3]        2*xi[1]*xi[3]^2  2*xi[1]*xi[2]*xi[3]  2*xi[1]*xi[2]^2  xi[2]^2*xi[3]^2        2*xi[1]*xi[2]*xi[3]^2  2*xi[1]*xi[2]^2*xi[3]  2*xi[1]*xi[2]^2*xi[3]^2
+        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]  0        2*xi[2]  0        xi[3]^2        2*xi[2]*xi[3]  0              2*xi[1]*xi[2]  0              xi[1]^2        2*xi[2]*xi[3]^2  xi[1]*xi[3]^2        2*xi[1]*xi[2]*xi[3]  0                xi[1]^2*xi[3]        2*xi[1]^2*xi[2]  2*xi[1]*xi[2]*xi[3]^2  xi[1]^2*xi[3]^2        2*xi[1]^2*xi[2]*xi[3]  2*xi[1]^2*xi[2]*xi[3]^2
+        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]  2*xi[3]  0        0        2*xi[2]*xi[3]  xi[2]^2        2*xi[1]*xi[3]  0              xi[1]^2        0              2*xi[2]^2*xi[3]  2*xi[1]*xi[2]*xi[3]  xi[1]*xi[2]^2        2*xi[1]^2*xi[3]  xi[1]^2*xi[2]        0                2*xi[1]*xi[2]^2*xi[3]  2*xi[1]^2*xi[2]*xi[3]  xi[1]^2*xi[2]^2        2*xi[1]^2*xi[2]^2*xi[3]
+    ]
+end
+
+###
+
+macro create_basis(T)
+    quote
+        T = $T
+        global get_basis, get_dbasis, length, size
+        X = get_reference_coordinates(T)
+        nbasis = length(X)
+        A = zeros(nbasis, nbasis)
+        for i=1:nbasis
+            A[i,:] = get_interpolation_polynomial(T, X[i])
+        end
+        invA = inv(A)
+        function get_basis(element::Element{$T}, ip, time)
+            return get_interpolation_polynomial($T, ip)*invA
+        end
+        function get_dbasis(element::Element{$T}, ip, time)
+            return get_interpolation_polynomial($T, ip, Val{:partial_derivatives})*invA
+        end
+    end
+end
+
+@create_basis Seg2
+@create_basis Seg3
+@create_basis Tri3
+@create_basis Tri6
+@create_basis Quad4
+@create_basis Quad8
+@create_basis Quad9
+@create_basis Tet4
+@create_basis Tet10
+@create_basis Hex8
+@create_basis Hex20
+@create_basis Hex27
+
diff --git a/src/lagrange_macro.jl b/src/lagrange_macro.jl
index e9b7ab2..08545be 100644
--- a/src/lagrange_macro.jl
+++ b/src/lagrange_macro.jl
@@ -30,7 +30,7 @@ Examples
 macro create_lagrange_element(element_name, element_description, X, P)
     eltype = esc(element_name)
     quote
-        global get_basis, length, size
+        global get_basis, length, size, get_reference_coordinates
         #=
                get_reference_element_coordinates,
                get_reference_element_midpoint
@@ -54,91 +54,16 @@ macro create_lagrange_element(element_name, element_description, X, P)
             return size($X, 2)
         end
 
-        #=
         XX = refcoords($X)
-        function get_reference_element_coordinates(::Type{$eltype})
+        function get_reference_coordinates(::Type{$eltype})
             return XX
         end
 
-        XXX = vec(mean($X, 2))
-        function get_reference_element_midpoint(::Type{$eltype})
-            return XXX
-        end
-
-        function $eltype(args...)
-            return Element{$eltype}(args...)
-        end
-
-        =#
-
     end
 end
 
-# 1d Lagrange elements
-
-@create_lagrange_element(Seg2, "2 node linear line element",
-    [-1.0 1.0], (xi) -> [1.0, xi[1]])
-
-@create_lagrange_element(Seg3, "3 node quadratic line element",
-    [-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
-
-# 2d Lagrange elements
-
-@create_lagrange_element(Tri3, "3 node bilinear triangle element",
-    [0.0 1.0 0.0
-     0.0 0.0 1.0],
-    (xi) -> [1.0, xi[1], xi[2]])
-
-@create_lagrange_element(Tri6, "6 node quadratic triangle element",
-    [0.0 1.0 0.0 0.5 0.5 0.0
-     0.0 0.0 1.0 0.0 0.5 0.5],
-    (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],
-    (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]])
-
-#=
-@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
-     0.0 0.0 1.0 0.0
-     0.0 0.0 0.0 1.0],
-    (xi) -> [1.0, xi[1], xi[2], xi[3]])
-
-@create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
-    [0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
-     0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5
-     0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.5 0.5],
-    (xi) -> [    1.0,   xi[1],       xi[2],       xi[3],     xi[1]^2,
-             xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
-
 function get_reference_element_midpoint{E}(element::Element{E})
     get_reference_element_midpoint(E)
 end
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index ca33f40..3f91da1 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -85,6 +85,20 @@ function get_nodal_vector(elements, field_name, time)
     return node_ids, field
 end
 
+""" Return nodal values in Dict format. """
+function get_nodal_dict(T::DataType, elements, field_name, time)
+    f = T()
+    for element in elements
+        for (c, v) in zip(get_connectivity(element), element(field_name, time))
+            if haskey(f, c)
+                @assert isapprox(f[c], v)
+            end
+            f[c] = v
+        end
+    end
+    return f
+end
+
 """ Update nodal field values from set of elements to another. Can be used to
 transform e.g. reaction force from boundary element set to surface of
 volume elements for easier postprocess.
@@ -114,3 +128,17 @@ function copy_field!(src_problem::Problem, dst_problem::Problem, field_name, tim
     copy_field!(src_problem.elements, dst_problem.elements, field_name, time)
 end
 
+""" Return field calculated to nodal points for elements in problem p. """
+function call(problem::Problem, field_name, time=0.0)
+    f = Dict()
+    for element in get_elements(problem)
+        for (c, v) in zip(get_connectivity(element), element(field_name, time))
+            if haskey(f, c)
+                @assert isapprox(f[c], v)
+            end
+            f[c] = v
+        end
+    end
+    return f
+end
+
diff --git a/src/preprocess.jl b/src/preprocess.jl
index aa40855..0590567 100644
--- a/src/preprocess.jl
+++ b/src/preprocess.jl
@@ -1,6 +1,17 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+#=
+- read meshes from different formats
+- reorder connectivity, create element sets, node sets, ...
+- create partitions for parallel runs
+- renumber elements / nodes
+- maybe precheck for bad elements
+- check surface normal direction in boundary elements
+- orientation of 2d elements
+- etc only topology related stuff
+=#
+
 importall Base
 
 using JuliaFEM
@@ -88,8 +99,25 @@ function create_elements(mesh::Mesh)
     return elements
 end
 
-function create_elements(mesh::Mesh, element_set::ASCIIString)
-    return create_elements(filter_by_element_set(mesh, element_set))
+function create_elements(mesh::Mesh, element_sets::ASCIIString...)
+    elements = Element[]
+    for element_set in element_sets
+        new_elements = create_elements(filter_by_element_set(mesh, element_set))
+        push!(elements, new_elements...)
+    end
+    return elements
+end
+
+function create_elements(mesh::Mesh, element_type::Symbol)
+    elements = Element[]
+    for (elid, elcon) in mesh.elements
+        eltype = mesh.element_types[elid]
+        eltype == element_type || continue
+        element = Element(JuliaFEM.(eltype), elcon)
+        update!(element, "geometry", mesh.nodes)
+        push!(elements, element)
+    end
+    return elements
 end
 
 """ find npts nearest nodes form mesh and return id numbers as list. """
@@ -104,3 +132,38 @@ function find_nearest_nodes(mesh::Mesh, coords::Vector, npts=1)
     return node_ids
 end
 
+"""
+Apply new node ordering to elements. In JuliaFEM same node ordering is used
+than in ABAQUS and if mesh is parsed from FEM format with other node ordering
+this can be used to reorder nodes.
+
+Parameters
+----------
+mapping :: Dict{Symbol, Vector{Int}}
+    e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
+
+"""
+function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{Int}})
+    for (elid, eltype) in mesh.element_types
+        haskey(mapping, eltype) || continue
+        new_order = mapping[eltype]
+        element_connectivity = mesh.elements[elid]
+        new_element_connectivity = element_connectivity[new_order]
+        mesh.elements[elid] = new_element_connectivity
+    end
+end
+
+"""
+Swap surface element connectivity s.t. normals point outward
+"""
+function check_orientation!
+    # TODO
+end
+
+"""
+Partition model using METIS
+"""
+function partition_model!
+    # TODO
+end
+
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index 429d980..8dc22b4 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -309,18 +309,6 @@ function get_element_sets(med::MEDFile, mesh_name)
     return es
 end
 
-global const med_elmap = Dict{Symbol, Vector{Int}}(
-    :PO1 => [1],
-    :SE2 => [1, 2],
-    :SE3 => [1, 2, 3],
-    :TR3 => [1, 2, 3],
-    :QU4 => [1, 2, 3, 4],
-    :TE4 => [3, 2, 1, 4],
-    :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
     increments = keys(med.data["ENS_MAA"][mesh_name])
@@ -340,12 +328,14 @@ function get_connectivity(med::MEDFile, elsets, mesh_name)
             eltype = Symbol(eltype)
             elco = element_connectivity[:, i]
             elset = Symbol(elsets[elset_ids[i]])
+#= to more general preprocess
             if haskey(med_elmap, eltype)
                 elco = elco[med_elmap[eltype]]
             else
                 warn("no element mapping info found for element type $eltype")
                 warn("consider this as a warning: element may have french nodal ordering")
             end
+=#
             d[element_ids[i]] = (eltype, elset, elco)
         end
     end
@@ -390,21 +380,56 @@ function parse_aster_med_file(fn::ASCIIString, mesh_name=nothing; debug=false)
     return result
 end
 
+# some glues about ordering, this is still a mystery..
+# http://onelab.info/pipermail/gmsh/2008/003850.html
+# http://caelinux.org/wiki/index.php/Proj:UNVConvert
+
+#global const med_connectivity = Dict{Symbol, Vector{Int}}(
+#    :Tet4 => [3, 2, 1, 4],
+#    :Hex8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
+#    :Tet10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8])
+
+global const med_connectivity = Dict{Symbol, Vector{Int}}(
+    :Tet4  => [4,3,1,2],
+    :Tet10 => [4,3,1,2,10,7,8,9,6,5],
+    :Hex8  => [4,8,7,3,1,5,6,2],
+    :Hex20 => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9],
+    :Hex27 => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9,24,25,26,23,21,22,27])
+
+# element names in CA -> element names in JuliaFEM
 global const mapping = Dict(
+
     :PO1 => :Poi1,
+
     :SE2 => :Seg2,
     :SE3 => :Seg3,
+    :SE4 => :Seg4,
+
     :TR3 => :Tri3,
     :TR6 => :Tri6,
+    :TR7 => :Tru6,
+
     :QU4 => :Quad4,
     :QU8 => :Quad8,
     :QU9 => :Quad9,
+
+    :TE4 => :Tet4,
+    :T10 => :Tet10,
+
+    :PE6 => :Penta6,
+    :P15 => :Penta15,
+    :P18 => :Penta18,
+
     :HE8 => :Hex8,
     :H20 => :Hex20,
-    :TE4 => :Tet4,
-    :T10 => :Tet10)
+    :H27 => :Hex27,
 
-function aster_read_mesh(fn::ASCIIString, mesh_name=nothing)
+    :PY5 => :Pyramid5,
+    :P13 => :Pyramid13,
+
+    )
+
+function aster_read_mesh(fn::ASCIIString, mesh_name=nothing; reorder_element_connectivity=true)
     result = parse_aster_med_file(fn, mesh_name)
     mesh = Mesh()
     for (nid, (nset, ncoords)) in result["nodes"]
@@ -416,6 +441,9 @@ function aster_read_mesh(fn::ASCIIString, mesh_name=nothing)
         add_element!(mesh, elid, mapping[eltype], elcon)
         add_element_to_element_set!(mesh, string(elset), elid)
     end
+    if reorder_element_connectivity
+        reorder_element_connectivity!(mesh, med_connectivity)
+    end
     return mesh
 end
 
diff --git a/src/problems.jl b/src/problems.jl
index efedf4f..0ed6480 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -266,10 +266,14 @@ function update_elements!{P<:BoundaryProblem}(problem::Problem{P}, u, la)
     end
 end
 
-function get_elements(problem)
+function get_elements(problem::Problem)
     return problem.elements
 end
 
+function length(problem::Problem)
+    return length(problem.elements)
+end
+
 function update!(problem::Problem, field_name::ASCIIString, field)
     update!(problem.elements, field_name, field)
 end
@@ -297,6 +301,12 @@ function push!(problem::Problem, elements::Vector)
     push!(problem.elements, elements...)
 end
 
+function push!(problem::Problem, elements_::Vector...)
+    for elements in elements_
+        push!(problem.elements, elements...)
+    end
+end
+
 function get_gdofs(element::Element, dim::Int)
     conn = get_connectivity(element)
     if length(conn) == 0
diff --git a/src/solvers.jl b/src/solvers.jl
index 4d412be..5aa4fc5 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -207,11 +207,11 @@ conditions are first eliminated before solution.
 """
 function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=false)
 
-    nnz(D) == 0 || return false
+    nnz(D) == 0 || return F, false
     nz = get_nonzero_rows(C2)
     B = get_nonzero_rows(C2')
     # C2^-1 exists or this doesn't work
-    length(nz) == length(B) || return false
+    length(nz) == length(B) || return F, false
 
     A = get_nonzero_rows(K)
     I = setdiff(A, B)
@@ -227,8 +227,7 @@ function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; F=nothing, debug=fals
     try
         u[B] = lufact(C2[nz,B]) \ full(g[nz])
     catch
-        info("solver #1 failed to solve boundary dofs (you should not see this message).")
-        return false
+        error("solver #1 failed to solve boundary dofs (you should not see this message).")
     end
 
     # solve interior domain using LDLt factorization
@@ -296,9 +295,7 @@ function solve_linear_system(solver::Solver; F=nothing, empty_assemblies_before_
     i = 0
     for i in [1, 2]
         F, status = solve!(K, C1, C2, D, f, g, u, la, Val{i}; F=F)
-        if status
-            break
-        end
+        status && break
     end
     status || error("Failed to solve linear system!")
 
@@ -467,11 +464,6 @@ Main differences in this solver, compared to nonlinear solver are:
 
 """
 type Linear <: AbstractSolver
-    norms :: Vector{Tuple}
-end
-
-function Linear()
-    solver = Linear([])
 end
 
 function assemble!(solver::Solver{Linear}; show_info=true)
@@ -526,3 +518,66 @@ end
 
 ### End of linear quasistatic solver
 
+### Postprocessor
+
+type Postprocessor <: AbstractSolver
+    assembly :: Assembly
+    F :: Union{Factorization, Void}
+end
+
+function Postprocessor()
+    Postprocessor(Assembly(), nothing)
+end
+
+function assemble!(solver::Solver{Postprocessor}; show_info=true)
+    show_info && info("Assembling problems ...")
+    tic()
+    nproblems = 0
+    ndofs = 0
+    assembly = solver.properties.assembly
+    empty!(assembly)
+    for problem in get_problems(solver)
+        for element in get_elements(problem)
+            postprocess!(assembly, problem, element, solver.time)
+        end
+        nproblems += 1
+        ndofs = max(ndofs, size(problem.assembly.K, 2))
+    end
+    solver.ndofs = ndofs
+    t1 = round(toq(), 2)
+    show_info && info("Assembled $nproblems problems in $t1 seconds. ndofs = $ndofs.")
+end
+
+function call(solver::Solver{Postprocessor}; show_info=true)
+    t0 = Base.time()
+    show_info && info(repeat("-", 80))
+    show_info && info("Starting postprocessor")
+    show_info && info("Increment time t=$(round(solver.time, 3))")
+    show_info && info(repeat("-", 80))
+    initialize!(solver)
+    assemble!(solver)
+    assembly = solver.properties.assembly
+    M = sparse(assembly.M)
+    f = sparse(assembly.f)
+    F = cholfact(M)
+    q = F \ f
+    t1 = round(Base.time()-t0, 2)
+    show_info && info("Postprocess of results ready in $t1 seconds.")
+    return q
+end
+
+""" Convenience function to call postprocessor. """
+function Postprocessor(problems::Problem...)
+    solver = Solver(Postprocessor, "default postprocessor")
+    if length(problems) != 0
+        push!(solver, problems...)
+    end
+    return solver
+end
+
+function Postprocessor(name::ASCIIString, problems::Problem...)
+    solver = Postprocessor(problems...)
+    solver.name = name
+    return solver
+end
+
diff --git a/test/testdata/primitives.hdf b/test/testdata/primitives.hdf
new file mode 100644
index 0000000..233bdc1
Binary files /dev/null and b/test/testdata/primitives.hdf differ
diff --git a/test/testdata/primitives.med b/test/testdata/primitives.med
new file mode 100644
index 0000000..851c4e3
Binary files /dev/null and b/test/testdata/primitives.med differ
diff --git a/test/testdata/rod_short.med b/test/testdata/rod_short.med
new file mode 100644
index 0000000..74e8071
Binary files /dev/null and b/test/testdata/rod_short.med differ