diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 37d29fe..b7d25ea 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -15,7 +15,7 @@ autodiffcache = ForwardDiffCache()
 include("common.jl")
 
 include("fields.jl")
-export DCTI
+export DCTI, Field
 #include("basis.jl")  # interpolation of discrete fields
 #include("symbolic.jl") # a thin symbolic layer for fields
 #include("types.jl")  # type definitions
@@ -36,7 +36,8 @@ export get_integration_points
 include("sparse.jl")
 
 include("problems.jl") # common problem routines
-export Problem
+export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
+       get_unknown_field_dimension
 
 include("elasticity.jl") # elasticity equations
 export Elasticity
@@ -79,7 +80,7 @@ module Postprocess
 include("xdmf.jl")
 export xdmf_new_temporal_collection, xdmf_new_grid,
        xdmf_new_mesh!, xdmf_new_nodal_field!,
-       xdmf_save_model, xdmf_new_model
+       xdmf_save_model, xdmf_new_model, xdmf_dump
 end
 
 """ JuliaFEM testing routines. """
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index 6902244..956d56e 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -22,16 +22,28 @@ function get_formulation_type(problem::Problem{Dirichlet})
     return problem.properties.formulation
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time::Real)
-
-    @assert problem.properties.dual_basis
+function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Element, time)
 
     # get dimension and name of PARENT field
+    nnodes = length(element)
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
     gdofs = get_gdofs(element, field_dim)
 
-    De, Me, Ae = get_dualbasis(element, time)
+  
+#   if problem.properties.formulation == :dual_basis
+        De, Me, Ae = get_dualbasis(element, time)
+#   else
+#       Ae = eye(nnodes)
+#       De = zeros(nnodes, nnodes)
+#       for (w, xi) in get_integration_points(element, Val{3})
+#           N = element(xi, time)
+#           detJ = element(xi, time, Val{:detJ})
+#           De += w*N'*N*detJ
+#       end
+#   end
+  
+#   De = Ae = eye(nnodes)
 
     # left hand side
     for i=1:field_dim
@@ -44,27 +56,19 @@ function assemble!(assembly::Assembly, problem::Problem{Dirichlet}, element::Ele
 
     # right hand side
     for (w, xi) in get_integration_points(element, Val{3})
-        J = element(xi, time, Val{:Jacobian})
-        JT = transpose(J)
-        if size(JT, 2) == 1  # plane problem
-            w *= norm(JT)
-        else
-            w *= norm(cross(JT[:,1], JT[:,2]))
-        end
+        detJ = element(xi, time, Val{:detJ})
         N = element(xi, time)
 
         for i=1:field_dim
             ldofs = gdofs[i:field_dim:end]
             if haskey(element, field_name*" $i")
                 g = element(field_name*" $i", xi, time)
-                if get_formulation_type(problem) == :incremental
-                    # if having incremental formulation need to add previous
-                    # displacement to rhs (solving increment Δu !
+                if true
                     haskey(element, "displacement") || continue
                     g_prev = element(field_name, xi, time)
                     g -= g_prev[i]
                 end
-                add!(assembly.g, ldofs, w*g*Ae*N')
+                add!(assembly.g, ldofs, w*g*Ae*N'*detJ)
             end
         end
 
diff --git a/src/elasticity.jl b/src/elasticity.jl
index c4459ea..63a8052 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -27,6 +27,7 @@ end
 function get_formulation_type(problem::Problem{Elasticity})
     # we are solving residual and add increment to previous solution vector
     return :incremental
+    #return :total
 end
 
 function assemble!(assembly::Assembly, problem::Problem{Elasticity}, element::Element, time::Real)
@@ -129,13 +130,22 @@ function assemble{El<:Union{Tri3,Tri6,Quad4}}(problem::Problem{Elasticity}, elem
         if props.finite_strain # add geometric stiffness
             Kt += w*BNL'*S2*BNL*detJ # geometric stiffness
         end
-        f -= w*BL'*S*detJ # internal force
+
+        if get_formulation_type(problem) == :incremental
+            f -= w*BL'*S*detJ # internal force
+        end
 
         # volume load
         if haskey(element, "displacement load")
             b = element("displacement load", xi, time)
             f += w*vec(N'*b)*detJ
         end
+        for i=1:dim
+            if haskey(element, "displacement load $i")
+                b = element("displacement load $i", xi, time)
+                f[i:dim:end] += w*vec(b*N)*detJ
+            end
+        end
 
     end
 
@@ -213,16 +223,13 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
         E = element("youngs modulus", xi, time)
         nu = element("poissons ratio", xi, time)
 
-        a = 1 - nu
-        b = 1 - 2*nu
-        c = 1 + nu
-        D = E/(b*c) .* [
-            a nu nu 0 0 0
-            nu a nu 0 0 0
-            nu nu a 0 0 0
-            0 0 0 b 0 0
-            0 0 0 0 b 0
-            0 0 0 0 0 b]
+        D = E/((1.0+nu)*(1.0-2.0*nu)) * [
+            1.0-nu nu nu 0.0 0.0 0.0
+            nu 1.0-nu nu 0.0 0.0 0.0
+            nu nu 1.0-nu 0.0 0.0 0.0
+            0.0 0.0 0.0 0.5-nu 0.0 0.0
+            0.0 0.0 0.0 0.0 0.5-nu 0.0
+            0.0 0.0 0.0 0.0 0.0 0.5-nu]
 
         # # PK2 stress tensor in voigt notation
         S = D*[strain[1,1]; strain[2,2]; strain[3,3]; 2*strain[2,3]; 2*strain[1,3]; 2*strain[1,2]]
@@ -249,6 +256,7 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
             BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
             BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
         end
+
         fill!(BNL, 0.0)
         for i=1:size(dN, 2)
             BNL[1, 3*(i-1)+1] = dN[1,i]
@@ -274,13 +282,22 @@ function assemble{El<:Union{Tet4, Tet10, Hex8}}(problem::Problem{Elasticity}, el
         if props.finite_strain
             Kt += w*BNL'*S3*BNL*detJ
         end
-        f -= w*BL'*S*detJ
+
+        if get_formulation_type(problem) == :incremental
+            f -= w*BL'*S*detJ
+        end
 
         # volume load
         if haskey(element, "displacement load")
             T = element("displacement load", ip, time)
             f += w*vec(T*N)*detJ
         end
+        for i=1:dim
+            if haskey(element, "displacement load $i")
+                b = element("displacement load $i", xi, time)
+                f[i:dim:end] += w*vec(b*N)*detJ
+            end
+        end
     end
 
     return Kt, f
diff --git a/src/elements.jl b/src/elements.jl
index 6c13080..676a0ca 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -6,13 +6,14 @@ abstract AbstractElement
 typealias Node Vector{Float64}
 
 type Element{E<:AbstractElement}
+    id :: Int
     connectivity :: Vector{Int}
     fields :: Dict{ASCIIString, Field}
     properties :: E
 end
 
-function Element{E<:AbstractElement}(::Type{E}, connectivity=[], fields=Dict(), properties...)
-    Element{E}(connectivity, fields, E(properties...))
+function Element{E<:AbstractElement}(::Type{E}, connectivity=[], id=-1, fields=Dict(), properties...)
+    Element{E}(id, connectivity, fields, E(properties...))
 end
 
 function getindex(element::Element, field_name::ASCIIString)
@@ -59,12 +60,8 @@ function call(element::Element, xi::Vector, time, ::Type{Val{:detJ}})
     end
 end
 
-function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
-    element(xi, time, Val{:Jacobian})
-end
-
 function call(element::Element, xi::Vector, time, ::Type{Val{:Grad}})
-    J = get_jacobian(element, xi, time)
+    J = element(xi, time, Val{:Jacobian})
     return inv(J)*get_dbasis(element, xi, time)
 end
 
@@ -72,6 +69,9 @@ function call(element::Element, field_name, xi::Vector, time, ::Type{Val{:Grad}}
     element(xi, time, Val{:Grad})*element[field_name](time)
 end
 
+#function get_jacobian{E}(element::Element{E}, xi::Vector, time=0.0)
+#    element(xi, time, Val{:Jacobian})
+#end
 #function get_basis(element::Element, xi::Vector, time=0.0)
 #    get_basis(element.properties, xi, time)
 #end
diff --git a/src/xdmf.jl b/src/xdmf.jl
index 578b9c9..f248900 100644
--- a/src/xdmf.jl
+++ b/src/xdmf.jl
@@ -2,6 +2,7 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 using LightXML
+using JuliaFEM
 
 # element codes: http://www.paraview.org/pipermail/paraview/2013-July/028859.html
 # > from  ./VTK/ThirdParty/xdmf2/vtkxdmf2/libsrc/XdmfTopology.h
@@ -124,3 +125,49 @@ function xdmf_save_model(xdoc, filename)
   save_file(xdoc, filename)
 end
 
+function xdmf_dump(all_elements, eltype, elsym, time=0.0, filename="/tmp/xdmf_result.xmf")
+    info("$(length(all_elements)) elements.")
+    xdoc, xmodel = xdmf_new_model()
+    coll = xdmf_new_temporal_collection(xmodel)
+    grid = xdmf_new_grid(coll; time=time)
+
+    Xg = Dict{Int64, Vector{Float64}}()
+    ug = Dict{Int64, Vector{Float64}}()
+    nids = Dict{Int64, Int64}()
+    for element in all_elements
+        conn = get_connectivity(element)
+        for (i, c) in enumerate(conn)
+            nids[c] = c
+        end
+        X = element("geometry", time)
+        for (i, c) in enumerate(conn)
+            Xg[c] = X[i]
+        end
+        haskey(element, "displacement") || continue
+        u = element("displacement", time)
+        for (i, c) in enumerate(conn)
+            ug[c] = u[i]
+        end
+    end
+    perm = sort(collect(keys(Xg)))
+    nodes = Vector{Float64}[Xg[i] for i in perm]
+    disp = Vector{Float64}[ug[i] for i in perm]
+    nids = Int[nids[i] for i in perm]
+    inids = Dict{Int64, Int64}()
+    for (i, nid) in enumerate(nids)
+        inids[nid] = i
+    end
+    elements = []
+    for element in all_elements
+        isa(element, eltype) || continue
+        conn = get_connectivity(element)
+        nconn = [inids[i] for i in conn]
+        push!(elements, (elsym, nconn))
+    end
+
+    xdmf_new_mesh!(grid, nodes, elements)
+    xdmf_new_nodal_field!(grid, "displacement", disp)
+    xdmf_save_model(xdoc, filename)
+    info("model dumped to $filename")
+end
+