diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index c9b8719..594157f 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -34,6 +34,7 @@ export @debug, set_debug_on!, set_debug_off!
 
 using ForwardDiff
 autodiffcache = ForwardDiffCache()
+export derivative, jacobian, hessian
 
 """ Simple linspace extension to arrays.
 
@@ -60,6 +61,7 @@ include("types.jl")  # type definitions
 include("elements.jl")
 include("lagrange.jl") # Lagrange elements
 #include("hierarchical.jl") # P-elements
+include("mortar_elements.jl") # Mortar elements
 
 ### EQUATIONS ###
 include("integrate.jl")  # default integration points for elements
@@ -69,6 +71,7 @@ include("problems.jl")
 
 ### FORMULATIION ###
 include("dirichlet.jl")
+include("mortar.jl")  # mortar projection
 include("heat.jl")
 include("elasticity.jl")
 
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index f80c8b8..4a220b7 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -15,7 +15,7 @@ type DirichletProblem <: BoundaryProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
     equations :: Vector{DirichletEquation}
-    element_mapping :: Dict{Element, Equation}
+#   element_mapping :: Dict{DataType, DataType}
     field_value :: Function
 end
 
@@ -36,14 +36,15 @@ Create u(X) = 0.0 boundary condition for three-dimensional elasticity problem:
 >>> u(X) = [0.0, 0.0, 0.0]
 >>> bc = DirichletProblem(3, u)
 """
-function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0])
-    element_mapping = nothing
-    if dimension == 1
-        element_mapping = Dict(
-            Seg2 => DBC2D2
-            )
-    end
-    DirichletProblem("reaction force", dimension, [], element_mapping, field_value)
+function DirichletProblem(dimension::Int=1, field_value::Function=(X)->[0.0,0.0,0.0], equations=[])
+#    element_mapping = nothing
+#    if dimension == 1
+#        element_mapping = Dict(
+#            Seg2 => DBC2D2
+#            )
+#    end
+    DirichletProblem("reaction force", dimension, equations, field_value)
+#    DirichletProblem("reaction force", dimension, [], element_mapping, field_value)
 end
 
 """ Dirichlet boundary condition element for 2 node line segment """
@@ -51,17 +52,22 @@ type DBC2D2 <: DirichletEquation
     element :: Seg2
     integration_points :: Vector{IntegrationPoint}
 end
-function DBC2D2(element::Seg2)
-    integration_points = default_integration_points(element)
-    if !haskey(element, "reaction force")
-        element["reaction force"] = zeros(1, 2)
-    end
+
+function Base.size(equation::DBC2D2)
+    return (1, 2)
+end
+
+#function DBC2D2(element::Seg2)
+function Base.convert(::Type{DirichletEquation}, element::Seg2)
+    integration_points = line3()
+    haskey(element, "reaction force") || (element["reaction force"] = zeros(1, 2))
     DBC2D2(element, integration_points)
 end
-Base.size(equation::DBC2D2) = (1, 2)
+
 
 function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number=0.0, problem=nothing)
     gdofs = get_gdofs(equation)
+#   info("gdofs = $gdofs")
     element = get_element(equation)
     basis = get_basis(element)
     detJ = det(basis)
@@ -69,11 +75,13 @@ function assemble!(assembly::Assembly, equation::DirichletEquation, time::Number
         w = ip.weight * detJ(ip)
         N = basis(ip, time)
         add!(assembly.stiffness_matrix, gdofs, gdofs, w*N'*N)
-        if !isa(problem, Void)
-            X = basis("geometry", ip, time)
-            u = problem.field_value(X)
-            add!(assembly.force_vector, gdofs, w*N'*u)
-        end
+#       info("added $(w*N'*N)")
+#        if !isa(problem, Void)
+#            X = basis("geometry", ip, time)
+#            u = problem.field_value(X)[1:length(gdofs)]
+#            add!(assembly.force_vector, gdofs, w*N'*u)
+#        end
     end
+#   info("assembly done")
 end
 
diff --git a/src/elasticity.jl b/src/elasticity.jl
index 0d2f5b0..2ba5bfe 100644
--- a/src/elasticity.jl
+++ b/src/elasticity.jl
@@ -19,8 +19,7 @@ Saint Venant-Kirchhoff material model, which is simply
 
     S(E) = λtr(E) + 2μE
 """
-function get_internal_energy(equation::Equation, ip::IntegrationPoint,
-                             time::Number, F::Matrix)
+function get_internal_energy(equation::ElasticityEquation, ip::IntegrationPoint, time::Number, F::Matrix)
     element = get_element(equation)
     basis = get_basis(element)
     dbasis = grad(basis)
@@ -72,8 +71,7 @@ https://en.wikipedia.org/wiki/Plane_stress
 https://en.wikipedia.org/wiki/Hooke's_law
 
 """
-function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint,
-                             time::Number; variation=nothing)
+function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint, time::Number; variation=nothing)
 
     element = get_element(equation)
     basis = get_basis(element)
@@ -82,9 +80,10 @@ function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint,
     u = basis("displacement", ip, time, variation)
     gradu = dbasis("displacement", ip, time, variation)
     F = I + gradu # deformation gradient
-
+    #info("Deformation gradient: $F")
     # residual vector - internal energy
     r = get_internal_energy(equation, ip, time, F)
+    #info("boundary element")
 
     # external forces - volume load
     if haskey(element, "displacement load")
@@ -94,40 +93,71 @@ function get_residual_vector(equation::ElasticityEquation, ip::IntegrationPoint,
 
     return vec(r)
 end
-has_residual_vector(equation::ElasticityEquation) = true
 
-### Problem 1 - plane elasticity ###
+### Plane stress elasticity ###
 
 abstract PlaneElasticityProblem <: ElasticityProblem
+abstract PlaneStressElasticityEquation <: ElasticityEquation
 
 type PlaneStressElasticityProblem <: PlaneElasticityProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
-    equations :: Array{ElasticityEquation, 1}
-    element_mapping :: Dict{DataType, DataType}
+    equations :: Vector{PlaneStressElasticityEquation}
 end
 
 function PlaneStressElasticityProblem(equations=[])
-    element_mapping = Dict(
-        Quad4 => CPS4)
-    return PlaneStressElasticityProblem("displacement", 2, equations, element_mapping)
+    return PlaneStressElasticityProblem("displacement", 2, equations)
 end
 
 ### Equations ###
 
-abstract PlaneElasticityEquation <: ElasticityEquation
-abstract PlaneStressElasticityEquation <: PlaneElasticityEquation
-
+""" 4-node plane stress element. """
 type CPS4 <: PlaneStressElasticityEquation
     element :: Quad4
     integration_points :: Array{IntegrationPoint, 1}
 end
-function CPS4(element::Quad4)
+
+function Base.size(equation::CPS4)
+    return (2, 4)
+end
+
+function Base.convert(::Type{PlaneStressElasticityEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    if !haskey(element, "displacement")
-        element["displacement"] = zeros(2, 4)
-    end
+    haskey(element, "displacement") || (element["displacement"] = zeros(2, 4))
     CPS4(element, integration_points)
 end
-Base.size(equation::CPS4) = (2, 4)
+
+""" Boundary element for plane stress problem for surface loads. """
+type CPS2 <: PlaneStressElasticityEquation
+    element :: Seg2
+    integration_points :: Vector{IntegrationPoint}
+end
+
+function Base.size(equation::CPS2)
+    return (2, 2)
+end
+
+function Base.convert(::Type{PlaneStressElasticityEquation}, element::Seg2)
+    integration_points = get_default_integration_points(element)
+    haskey(element, "displacement") || (element["displacement"] = zeros(2, 2))
+    CPS2(element, integration_points)
+end
+
+function get_residual_vector(equation::CPS2, ip::IntegrationPoint, time::Number; variation=nothing)
+
+    element = get_element(equation)
+    basis = get_basis(element)
+
+    u = basis("displacement", ip, time, variation)
+    r = zeros(size(equation))
+
+    if haskey(element, "displacement traction force")
+        T = basis("displacement traction force", ip, time)
+#       info("traction force = $T")
+#       info("basis = $(basis(ip, time))")
+        r -= T*basis(ip, time)
+    end
+
+    return vec(r)
+end
 
diff --git a/src/equations.jl b/src/equations.jl
index 9c89b80..3190025 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -97,6 +97,12 @@ function get_gdofs(equation::Equation)
     return gdofs
 end
 
+function get_gdofs(element::Element, dim::Int)
+    conn = get_connectivity(element)
+    gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
+    return gdofs
+end
+
 """ Assemble element. """
 function assemble!(assembly::Assembly, equation::Equation, time::Number=0.0, problem=nothing)
 
diff --git a/src/heat.jl b/src/heat.jl
index 6c56410..e1e746d 100644
--- a/src/heat.jl
+++ b/src/heat.jl
@@ -72,46 +72,46 @@ end
 """ Diffusive heat transfer for 4-node bilinear element. """
 type DC2D4 <: HeatEquation
     element :: Quad4
-    integration_points :: Array{IntegrationPoint, 1}
+    integration_points :: Vector{IntegrationPoint}
 end
-function DC2D4(element::Quad4)
-    integration_points = get_default_integration_points(element)
-    if !haskey(element, "temperature")
-        element["temperature"] = zeros(4)
-    end
-    DC2D4(element, integration_points)
+
+function Base.size(equation::DC2D4)
+    return (1, 4)
 end
-Base.size(equation::DC2D4) = (1, 4)
 
 """ Diffusive heat transfer for 2-node linear segment. """
 type DC2D2 <: HeatEquation
     element :: Seg2
     integration_points :: Vector{IntegrationPoint}
 end
-function DC2D2(element::Seg2)
+
+function Base.size(equation::DC2D2)
+    return (1, 2)
+end
+
+# Conversions element -> equation
+
+function Base.convert(::Type{HeatEquation}, element::Quad4)
     integration_points = get_default_integration_points(element)
-    if !haskey(element, "temperature")
-        element["temperature"] = zeros(2)
-    end
+    haskey(element, "temperature") || (element["temperature"] = 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))
     DC2D2(element, integration_points)
 end
-Base.size(equation::DC2D2) = (1, 2)
 
 ### Problems ###
 
 type PlaneHeatProblem <: HeatProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
-    equations :: Vector{Equation}
-    #element_mapping :: Dict{Element, Equation}
-    # FIXME: Why is not working ^
-    element_mapping :: Dict{Any, Any}
+    equations :: Vector{HeatEquation}
 end
 
 """ Default constructor for problem takes no arguments. """
-function PlaneHeatProblem()
-    element_mapping = Dict(
-        Quad4 => DC2D4,
-        Seg2 => DC2D2)
-    return PlaneHeatProblem("temperature", 1, [], element_mapping)
+function PlaneHeatProblem(equations=[])
+    return PlaneHeatProblem("temperature", 1, equations)
 end
diff --git a/src/integrate.jl b/src/integrate.jl
index f5de31e..2cb807c 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -11,12 +11,20 @@ function get_default_integration_points(element::Quad4)
     ]
 end
 
-function get_default_integration_points(element::Seg2)
+
+function line1()
     [
         IntegrationPoint([0.0], 2.0)
     ]
 end
 
+function line2()
+    [
+        IntegrationPoint([-sqrt(1/3)], 1)
+        IntegrationPoint([+sqrt(1/3)], 1)
+    ]
+end
+
 function line3()
     [
         IntegrationPoint([0.0], 8/9),
@@ -25,6 +33,15 @@ function line3()
     ]
 end
 
+function line4()
+    [
+        IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+        IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+        IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+        IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+    ]
+end
+
 function line5()
     [
         IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
@@ -35,29 +52,10 @@ function line5()
     ]
 end
 
-    #integration_points = [
-    #    IntegrationPoint([ 0.0000000000000000], 0.5688888888888889),
-    #    IntegrationPoint([-0.5384693101056831], 0.4786286704993665),
-    #    IntegrationPoint([ 0.5384693101056831], 0.4786286704993665),
-    #    IntegrationPoint([-0.9061798459386640], 0.2369268850561891),
-    #    IntegrationPoint([ 0.9061798459386640], 0.2369268850561891)
-    #]
-    #integration_points = [
-    #    IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
-    #    IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
-    #    IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
-    #    IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
-    #]
-    #integration_points = [
-    #    IntegrationPoint([0.0], 8/9),
-    #    IntegrationPoint([-sqrt(3/5)], 5/9),
-    #    IntegrationPoint([+sqrt(3/5)], 5/9)
-    #]
-    #integration_points = [
-    #    IntegrationPoint([-sqrt(1/3)], 1)
-    #    IntegrationPoint([+sqrt(1/3)], 1)
-    #]
-    #integration_points = [
-    #    IntegrationPoint([0.0], 2)
-    #]
+function get_default_integration_points(element::Seg2)
+    return line1()
+end
 
+function get_default_integration_points(element::MSeg2)
+    return line3()
+end
diff --git a/src/mortar.jl b/src/mortar.jl
index 871080a..900ff30 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -1,78 +1,86 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# Mortar projection integration
+# Mortar equations
 
 abstract MortarEquation <: Equation
 
+function get_unknown_field_name(equation::MortarEquation)
+    return "reaction force"
+end
+
+""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """
+type MBC2D2 <: MortarEquation
+    element :: MSeg2
+    integration_points :: Vector{IntegrationPoint}
+end
+
+function Base.size(equation::MBC2D2)
+    return (1, 2)
+end
+
+function Base.convert(::Type{MortarEquation}, element::MSeg2)
+    return MBC2D2(element, get_default_integration_points(element))
+end
+
+# Mortar problem
+
 """
 Parameters
 ----------
 node_csys
     coordinate system in node, normal + tangent + "binormal"
-element_pairs
-    m x s matrix of boolean values, indicating elements sharing
-    common surface. s is number of slave elements and m is number
-    of master elements.
+    in 3d 3x3 matrix, in 2d 2x2 matrix, respectively
 """
 type MortarProblem <: BoundaryProblem
     unknown_field_name :: ASCIIString
     unknown_field_dimension :: Int
     equations :: Vector{MortarEquation}
-    element_mapping :: Dict{Element, MortarEquation}
-    master_elements :: Vector{Element}  # mortar surface
-    node_csys :: Dict{Int, Matrix{Float64}}
-    element_pairs :: Matrix{Bool}
 end
 
-function MortarProblem(dimension::Int=1, equations=[], master_elements=[])
-    element_mapping = Dict(
-        Seg2 => MBC2D2,
-    )
-    MortarProblem("reaction force", dimension, equations, element_mapping, master_elements, Dict(), zeros(0,0))
+function MortarProblem(dimension::Int=1, equations=[])
+    MortarProblem("reaction force", dimension, equations)
 end
 
-""" Mortar boundary condition element for 2-dimensional problem, 2 node line segment. """
-type MBC2D2 <: MortarEquation
-    element :: Seg2  # == non-mortar surface element
-    integration_points :: Vector{IntegrationPoint}
-end
-function MBC2D2(element::Seg2)
-    integration_points = default_integration_points(element)
-    if !haskey(element, "reaction force")
-        element["reaction force"] = zeros(1, 2)
-    end
-    MBC2D2(element, integration_points)
-end
-Base.size(equation::MBC2D2) = (1, 2)
+# Mortar projection calculation
 
-function find_master_elements(slave_element, problem)
-    # find slave element "position" in element pairs matrix
-    all_elements = map((equation) -> get_element(equation), problem.equations)
-    seid = findfirst(slave_element, all_elements)
-    info("slave element id = $seid")
-    # find master element "positions" in element pairs matrix
-    meids = find(problem.element_pairs[:, seid])
-    info("master element ids = $meids")
-    # master elements
-    master_elements = problem.master_elements[meids]
-    return master_elements
+""" Find master or "mortar" elements for this slave element. """
+function get_master_elements(element::MortarElement)
+    return element.master_elements
 end
 
-function calculate_local_assembly!(assembly::LocalAssembly, equation::MortarEquation, unknown_field_name::ASCIIString, time::Number=0.0, problem=nothing)
-    # slave element = non-mortar element where integration happens
-    # master element = mortar element projected to non-mortar side
-    isa(problem, Void) && error("Cannot create projection without problem")
-    initialize_local_assembly!(assembly, equation)
+function assemble!(assembly::Assembly, equation::MortarEquation, time::Number=0.0, problem=nothing)
     slave_element = get_element(equation)
-    basis = get_basis(slave_element)
-    detJ = det(basis)
-    master_elements = find_master_elements(equation, problem)
+    master_elements = get_master_elements(slave_element)
+    slave_basis = get_basis(slave_element)
+    detJ = det(slave_basis)
+    dim = size(equation, 1) # number of nodes
+    slave_dofs = get_gdofs(slave_element, dim)
     for master_element in master_elements
-        for ip in get_integration_points(slave_element)
-            mortar_basis = 0 # ...
-            assembly.stiffness_matrix += w*basis'*basis
-            assembly.force_vector += w*N'*gn
+        master_dofs = get_gdofs(master_element, dim)
+        xi1a = project_from_master_to_slave(slave_element, master_element, [-1.0])
+        xi1b = project_from_master_to_slave(slave_element, master_element, [ 1.0])
+        xi1 = clamp([xi1a xi1b], -1.0, 1.0)
+        l = 1/2*(xi1[2]-xi1[1])
+        if abs(l) < 1.0e-6
+            warn("No contribution")
+            continue # no contribution
+        end
+        master_basis = get_basis(master_element)
+        for ip in get_integration_points(equation)
+            w = ip.weight*detJ(ip)*l
+
+            # integration point on slave side segment
+            xi_gauss = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
+            # projected integration point
+            xi_projected = project_from_slave_to_master(slave_element, master_element, xi_gauss)
+
+            # add contribution to left hand side 
+            N1 = slave_basis(xi_gauss, time)
+            N2 = master_basis(xi_projected, time)
+            add!(assembly.lhs, slave_dofs, slave_dofs, w*N1'*N1)
+            add!(assembly.lhs, slave_dofs, master_dofs, -w*N1'*N2)
+
         end
     end
 end
diff --git a/src/mortar_elements.jl b/src/mortar_elements.jl
new file mode 100644
index 0000000..52263c2
--- /dev/null
+++ b/src/mortar_elements.jl
@@ -0,0 +1,84 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Mortar elements
+
+# slave element = non-mortar element where integration happens
+# master element = mortar element projected to non-mortar side
+
+abstract MortarElement <: Element
+
+type MSeg2 <: MortarElement
+    connectivity :: Vector{Int}
+    basis :: Basis
+    fields :: FieldSet
+    master_elements :: Vector{MortarElement}
+end
+
+function MSeg2(connectivity, master_elements=[], biorthogonal=false)
+    basis(xi) = [(1-xi[1])/2 (1+xi[1])/2]
+    dbasisdxi(xi) = [-1/2 1/2]
+    return MSeg2(connectivity, Basis(basis, dbasisdxi), FieldSet(), master_elements)
+end
+
+""" Find projection from slave nodes to master element, i.e. find xi2 from
+master element corresponding to the xi1.
+"""
+function project_from_slave_to_master(slave::MortarElement, master::MortarElement, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
+    slave_basis = get_basis(slave)
+    master_basis = get_basis(master)
+
+    # slave side geometry and normal direction at xi1
+    X1 = slave_basis("geometry", xi1, time)
+    N1 = slave_basis("nodal ntsys", xi1, time)[:,1]
+
+    # master side geometry at xi2
+    X2(xi2) = master_basis("geometry", [xi2], time)
+#   dX2(xi2) = dmaster_basis("geometry", xi2, time)
+
+    # equation to solve
+    R(xi2) = det([X2(xi2)-X1  N1]')
+#   dR(xi2) = det([dX2(xi2) N1]')
+    dR = ForwardDiff.derivative(R)
+
+    # go!
+    xi2 = 0.0
+    for i=1:max_iterations
+        dxi2 = -R(xi2) / dR(xi2)
+        xi2 += dxi2
+        if norm(dxi2) < tol
+            return Float64[xi2]
+        end
+    end
+    error("find projection from slave to master: did not converge")
+end
+
+""" Find projection from master surface to slave point, i.e. find xi1 from slave element corresponding to the xi2. """
+function project_from_master_to_slave(slave::MortarElement, master::MortarElement, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
+    slave_basis = get_basis(slave)
+    master_basis = get_basis(master)
+
+    # slave side geometry and normal direction at xi1
+    X1(xi1) = slave_basis("geometry", [xi1], time)
+    N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
+
+    # master side geometry at xi2
+    X2 = master_basis("geometry", xi2, time)
+
+    # equation to solve
+    R(xi1) = det([X1(xi1)-X2  N1(xi1)]')
+#   dR(xi1) = det([dX1(xi1) N1(xi1)]') + det([X1(xi1)-X2 dN1(xi1)]')
+    dR = ForwardDiff.derivative(R)
+
+    # go!
+    xi1 = 0.0
+    for i=1:max_iterations
+        dxi1 = -R(xi1) / dR(xi1)
+        xi1 += dxi1
+        if norm(dxi1) < tol
+            return Float64[xi1]
+        end
+    end
+    error("find projection from master to slave: did not converge")
+end
+
diff --git a/src/problems.jl b/src/problems.jl
index 17a829e..6d2f4a7 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -33,10 +33,11 @@ Notes
 Equation is automatically created during process based on problem
 element -> equation mapping and element type.
 """
-function Base.push!(problem::Problem, element::Element)
-    element_type = typeof(element)
-    equation_type = problem.element_mapping[element_type]
-    push!(problem.equations, equation_type(element))
+function Base.push!(problem::Problem, element::Element, args...)
+#    element_type = typeof(element)
+#    equation_type = problem.element_mapping[element_type]
+#    push!(problem.equations, equation_type(element, args...))
+    push!(problem.equations, element)
 end
 
 """
diff --git a/src/solvers.jl b/src/solvers.jl
index d376360..c2f7d45 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -111,25 +111,26 @@ common situation, i.e., some main field problem and it's Dirichlet boundary.
     Cu       = g
 
 """
-function call(solver::SimpleSolver, time::Number=Inf)
-    p1, p2 = get_problems(solver)
+function call(solver::SimpleSolver, time::Number=0.0)
+    problem1, problem2 = get_problems(solver)
 
-    ga1 = initialize_global_assembly(p1)
-    calculate_global_assembly!(ga1, p1)
-    ga2 = initialize_global_assembly(p2)
-    calculate_global_assembly!(ga2, p2)
+    assembly1 = Assembly()
+    assemble!(assembly1, problem1, time)
+    assembly2 = Assembly()
+    assemble!(assembly2, problem2, time)
 
-    A1 = ga1.stiffness_matrix
-    b1 = ga1.force_vector
-    A2 = ga2.stiffness_matrix
-    b2 = ga2.force_vector
+#   info("Creating sparse matrices")
+    A1 = sparse(assembly1.stiffness_matrix)
+    b1 = sparse(assembly1.force_vector, size(A1, 1), 1)
+    A2 = sparse(assembly2.stiffness_matrix)
+    b2 = sparse(assembly2.force_vector, size(A2, 1), 1)
     
     # create a saddle point problem
     A = [A1 A2; A2' zeros(A2)]
     b = [b1; b2]
 
     # solve problem
-    nz = unique(rowvals(A))  # here we remove any zero rows
+    nz = unique(rowvals(A))  # take only non-zero rows
     x = zeros(b)
     x[nz] = lufact(A[nz,nz]) \ full(b[nz])
 
@@ -138,23 +139,23 @@ function call(solver::SimpleSolver, time::Number=Inf)
     x2 = x[length(b1)+1:end]
 
     # update field for elements in problem 1
-    for equation in get_equations(p1)
+    for equation in get_equations(problem1)
         element = get_element(equation)
-        field_name = get_unknown_field_name(p1)
-        gdofs = get_gdofs(p1, equation)
+        field_name = get_unknown_field_name(problem1)
+        gdofs = get_gdofs(problem1, equation)
         element_solution = full(x1[gdofs])
-        field = Field(time, element_solution)
-        push!(element[field_name], field)
+        field = Increment(element_solution)
+        push!(element[field_name], TimeStep(time, field))
     end
 
     # update field for elements in problem 2 (Dirichlet boundary)
-    for equation in get_equations(p2)
+    for equation in get_equations(problem2)
         element = get_element(equation)
-        field_name = get_unknown_field_name(p2)
-        gdofs = get_gdofs(p2, equation)
+        field_name = get_unknown_field_name(problem2)
+        gdofs = get_gdofs(problem2, equation)
         element_solution = full(x2[gdofs])
-        field = Field(time, element_solution)
-        push!(element[field_name], field)
+        field = Increment(element_solution)
+        push!(element[field_name], TimeStep(time, field))
     end
 end
 
diff --git a/src/sparse.jl b/src/sparse.jl
index e834800..1e622f1 100644
--- a/src/sparse.jl
+++ b/src/sparse.jl
@@ -14,8 +14,8 @@ function SparseMatrixIJV()
     SparseMatrixIJV([], [], [])
 end
 
-function Base.sparse(A::SparseMatrixIJV)
-    return sparse(A.I, A.J, A.V)
+function Base.sparse(A::SparseMatrixIJV, args...)
+    return sparse(A.I, A.J, A.V, args...)
 end
 
 function Base.push!(A::SparseMatrixIJV, I::Int, J::Int, V::Float64)
@@ -37,7 +37,7 @@ function Base.append!(A::SparseMatrixIJV, I::Vector{Int}, J::Vector{Int}, V::Vec
 end
 
 function Base.full(A::SparseMatrixIJV, args...)
-    return full(sparse(A.I, A.J, A.V), args...)
+    return full(sparse(A.I, A.J, A.V, args...))
 end
 
 """ Add local element matrix to sparse matrix. This basically does:
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index 0704c98..528b9aa 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -4,27 +4,47 @@
 module ElasticityTests
 
 using JuliaFEM.Test
-using JuliaFEM: Quad4, Field, FieldSet, CPS4,
+using JuliaFEM: Seg2, Quad4, Field, FieldSet, CPS4,
                 get_basis, solve!,
                 PlaneStressElasticityProblem
 
 
-function test_elasticity_one_element()
+function test_elasticity_volume_load()
     element = Quad4([1, 2, 3, 4])
     element["geometry"] = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
     element["youngs modulus"] = 500.0
     element["poissons ratio"] = 0.3
     element["displacement load"] = Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]
-    equation = CPS4(element)
     free_dofs = [3, 4, 5, 6]
-    problem = PlaneStressElasticityProblem([equation])
+    problem = PlaneStressElasticityProblem()
+    push!(problem, element)
     solve!(problem, free_dofs; max_iterations=10)
-    #solve!(equation, "displacement", free_dofs; max_iterations=10)
     disp = get_basis(element)("displacement", [1.0, 1.0])[2]
     info("displacement at tip: $disp")
     # verified using Code Aster.
     @test isapprox(disp, -8.77303119819776)
 end
 
+function test_elasticity_surface_load()
+    N = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
+
+    element1 = Quad4([1, 2, 3, 4])
+    element1["geometry"] = Vector[N[1], N[2], N[3], N[4]]
+    element1["youngs modulus"] = 500.0
+    element1["poissons ratio"] = 0.3
+    element2 = Seg2([3, 4])
+    element2["geometry"] = Vector[N[3], N[4]]
+    element2["displacement traction force"] = Vector[[0.0, -10.0], [0.0, -10.0]]
+
+    free_dofs = [3, 4, 5, 6]
+    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]
+    info("displacement at tip: $disp")
+    # verified using Code Aster.
+    @test isapprox(disp, -9.33106637611714)
+end
 
 end
diff --git a/test/test_heat.jl b/test/test_heat.jl
index e654cb2..d47b33b 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -6,6 +6,8 @@
 module HeatTests  # always wrap tests to module ending with "Tests"
 
 using JuliaFEM.Test  # always use JuliaFEM.Test, not Base.Test
+
+using JuliaFEM: HeatEquation
 using JuliaFEM: Seg2, Quad4, DC2D4, DC2D2, Assembly, assemble!
 
 function test_one_element()  # always start test function with name test_
@@ -25,7 +27,7 @@ function test_one_element()  # always start test function with name test_
 
     # 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 = DC2D4(element)
+    equation = convert(HeatEquation, element)
     #la = initialize_local_assembly()
     #calculate_local_assembly!(la, equation, "temperature")
     assembly = Assembly()
@@ -37,7 +39,7 @@ function test_one_element()  # always start test function with name test_
 
     # Set constant flux g=6 on boundary. Accurate solution is
     # u(x,y) = x which equals T=1 on boundary.
-    boundary_equation = DC2D2(boundary_element)
+    boundary_equation = convert(HeatEquation, boundary_element)
     empty!(assembly)
 
     time = 1.0
diff --git a/test/test_mortar.jl b/test/test_mortar.jl
index dab6820..9a788ca 100644
--- a/test/test_mortar.jl
+++ b/test/test_mortar.jl
@@ -6,78 +6,113 @@ module MortarTests
 using JuliaFEM
 using JuliaFEM.Test
 
-function test_calc_flat_2d_assembly()
-    # this is hand calculated and given example in my thesis
+using JuliaFEM: MSeg2, Seg2, MortarProblem, MortarEquation, MortarElement, Assembly, assemble!
+using JuliaFEM: get_basis, grad, project_from_slave_to_master, project_from_master_to_slave
+
+function get_test_2d_model()
+    # this is hand calculated and given as an example in my thesis
     N = Vector[
         [0.0, 2.0], [1.0, 2.0], [2.0, 2.0],
         [0.0, 0.0], [1.0, 0.0], [2.0, 0.0],
         [0.0, 1.0], [5/4, 1.0], [2.0, 1.0],
         [0.0, 1.0], [3/4, 1.0], [2.0, 1.0]]
+    rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
+    slave1 = MSeg2([10, 11])
+    slave1["geometry"] = Vector[N[10], N[11]]
+    # should be n = [0 -1]' and t = [1 0]'
+    slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    slave2 = MSeg2([11, 12])
+    slave2["geometry"] = Vector[N[11], N[12]]
+    # should be n = [0 -1]' and t = [1 0]'
+    slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    master1 = MSeg2([7, 8])
+    master1["geometry"] = Vector[N[7], N[8]]
+    master2 = MSeg2([8, 9])
+    master2["geometry"] = Vector[N[8], N[9]]
+    push!(slave1.master_elements, master1)
+    push!(slave1.master_elements, master2)
+    push!(slave2.master_elements, master1)
+    push!(slave2.master_elements, master2)
+    return [slave1, slave2], [master1, master2]
+end
 
-    slave1 = Seg2([10, 11])
-    slave1["geometry"] = Vector[N10, N11]
 
-    slave2 = Seg2([11, 12])
-    slave2["geometry"] = Vector[N11, N12]
+function test_calc_flat_2d_projection()
+    slaves, masters = get_test_2d_model()
+    slave1, slave2 = slaves
+    master1, master2 = masters
 
-    master1 = Seg2([7, 8])
-    master1["geometry"] = Vector[N7, N8]
+    xi2a = project_from_slave_to_master(slave1, master1, [-1.0])
+    @test xi2a == [-1.0]
 
-    master2 = Seg2([8, 9])
-    master2["geometry"] = Vector[N8, N9]
+    xi2b = project_from_slave_to_master(slave1, master1, [1.0])
+    @test xi2b == [ 0.2]
+    X2 = get_basis(master1)("geometry", xi2b)
+    @test X2 == [3/4, 1.0]
 
+    xi1a = project_from_master_to_slave(slave1, master1, [-1.0])
+    @test xi1a == [-1.0]
+    xi1b = project_from_master_to_slave(slave1, master1, [1.0])
+    X1 = get_basis(slave1)("geometry", xi1b)
+    @test X1 == [5/4, 1.0]
+end
+
+function test_create_flat_2d_assembly()
+    slaves, masters = get_test_2d_model()
+    slave1, slave2 = slaves
+    master1, master2 = masters
+
+    info("creating problem")
     problem = MortarProblem()
+    info("pushing slave elements to problem")
     push!(problem, slave1)
     push!(problem, slave2)
-    push!(problem.master_elements, master1)
-    push!(problem.master_elements, master2)
-
-    rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
-    # should be n = [0 -1]' and t = [1 0]'
-    @test isapprox(rotation_matrix(-phi/2), [[0 -1]' [1 0]'])
-
-    problem.node_csys = Dict(
-        10 => rotation_matrix(-phi/2),
-        11 => rotation_matrix(-phi/2),
-        12 => rotation_matrix(-phi/2))
-
-    # first index = master element id
-    # second index = slave element id
-    problem.element_pairs = zeros(2, 2)
-    # first slave element connects to master element 1
-    problem.element_pairs[1, 1] = true
-    # second slave element connects to master element 1
-    problem.element_pairs[1, 2] = true
-    # second slave element connects to master element 2
-    problem.element_pairs[2, 2] = true
-
-    B_expected = zeros(12, 9)
 
+    B_expected = zeros(12, 12)
     S1 = [10, 11]
     M1 = [7, 8]
     B_expected[S1,S1] += [1/4 1/8; 1/8 1/4]
-    B_expected[S1,M1] += [3/10 3/40; 9/40 3/20]
+    B_expected[S1,M1] -= [3/10 3/40; 9/40 3/20]
 
-    la = initialize_local_assembly(problem)
-    calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem)
-    B = full(la.lhs)
+    info("creating assembly")
+    assembly = Assembly()
+    assemble!(assembly, problem.equations[1], 0.0, problem)
+    B = round(full(assembly.lhs, 12, 12), 6)
+    info("size of B = $(size(B))")
+    info("B matrix in first slave element = \n$(B[10:11,:])")
+    info("B matrix expected = \n$(B_expected[10:11,:])")
     @test isapprox(B, B_expected) 
 
     fill!(B_expected, 0.0)
+    empty!(assembly)
 
     S2 = [11, 12]
     M2 = [7, 8]
     B_expected[S2,S2] += [49/150 11/150; 11/150 2/75]
-    B_expected[S2,M2] += [13/150 47/150; 1/75 13/150]
+    B_expected[S2,M2] -= [13/150 47/150; 1/75 13/150]
     S3 = [11, 12]
     M3 = [8, 9]
     B_expected[S3,S3] += [9/100 27/200; 27/200 39/100]
-    B_expected[S3,M3] += [3/20 3/40; 9/40 3/10]
+    B_expected[S3,M3] -= [3/20 3/40; 9/40 3/10]
+    assemble!(assembly, problem.equations[2], 0.0, problem)
+    B = full(assembly.lhs)
+    info("size of B = $(size(B))")
+    info("B matrix in second slave element = \n$(B[11:12,:])")
+    info("B matrix expected = \n$(B_expected[11:12,:])")
 
-    la = initialize_local_assembly(problem)
-    calculate_local_assembly!(la, problem.equations[1], "reaction force", 0.0, problem=problem)
-    B = full(la.lhs)
     @test isapprox(B, B_expected)
 end
 
+function test_patch_test_heat_2d()
+
+
+    
+    slaves, masters = get_test_2d_model()
+    problem2 = MortarProblem()
+    for slave in slaves:
+        push!(problem2, slave)
+    end
+
+end
+
 end
diff --git a/test/test_solver.jl b/test/test_solver.jl
index e8f4156..15364c2 100644
--- a/test/test_solver.jl
+++ b/test/test_solver.jl
@@ -1,9 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# test SimpleSolver
+module SolverTests
+
+using JuliaFEM.Test
+using JuliaFEM
 
-using FactCheck
 using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, get_element, get_basis
 
 """ Define Problem 1:
@@ -13,25 +15,16 @@ using JuliaFEM: DirichletProblem, Seg2, PlaneHeatProblem, Quad4, SimpleSolver, g
 """
 function get_heatproblem()
     el1 = Quad4([1, 2, 3, 4])
-    # these might look like normal values but believe me, they
-    # are fields with temporal and spatial dimension
     el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
     el1["temperature thermal conductivity"] = 6.0
     el1["density"] = 36.0
 
     el2 = Seg2([1, 2])
     el2["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
-    # Boundary load, linear ramp 0 -> 600 at time 0 -> 1
-    # yet another simplification, if field is given as a tuple,
-    # multiple fields are created. there is 1 second time step between
-    # each field. So the following is basically same as
-    # fieldset = FieldSet("temperature flux")
-    # field1 = Field(0.0, 0.0)
-    # field2 = Field(1.0, 600.0)
-    # push!(fieldset, field1)
-    # push!(fieldset, field2)
-    # element["temperature flux"] = fieldset
-    el2["temperature flux"] = (0.0, 600.0)
+    el2["temperature flux"] = (
+        (0.0 => 0.0),
+        (1.0 => 600.0)
+        )
 
     problem1 = PlaneHeatProblem()
     push!(problem1, el1)
@@ -50,13 +43,17 @@ function get_boundaryproblem()
     return problem2
 end
 
-facts("test simplesolver") do
+function test_simplesolver()
+    info("construct heat problem")
     problem1 = get_heatproblem()
+    info("construct boundary problem")
     problem2 = get_boundaryproblem()
     # Create a solver for a set of problems
+    info("create SimpleSolver with problems.")
     solver = SimpleSolver()
     push!(solver, problem1)
     push!(solver, problem2)
+    info("solve!")
     # Solve problem at time t=1.0 and update fields
     call(solver, 1.0)
     # Postprocess.
@@ -66,6 +63,8 @@ facts("test simplesolver") do
     basis = get_basis(el2)
     X = basis("geometry", xi, 1.0)
     T = basis("temperature", xi, 1.0)
-    Logging.info("Temperature at point X = $X is T = $T")
-    @fact T --> roughly(100.0)
+    info("Temperature at point X = $X is T = $T")
+    @test isapprox(mean(T), 100.0)
+end
+
 end