diff --git a/src/mortar.jl b/src/mortar.jl
index 353ea36..5af1ed0 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -1,137 +1,52 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# Mortar projection calculation for 2d
+"""
+Currently two strategies exists:
 
+a) Remove inactive inequality constraints in element level. This is done in
+   assemble! if normal_condition is set to :Contact. For some reason this
+   leads to convergence issues.
+b) Remove inactive inequality constraints in assembly level. This is done in
+   posthook algorithm if inequality_constraints is set to true. This gives
+   more robust behavior.
+
+   Either use inequality_constraints=True OR :Contact + :Slip, but do not mix.
+
+   minimum_distance can be used to roughly skip integration of mortar
+   projections for elements that are "far enough" from each other. Increases
+   performance.
+
+"""
+type Mortar <: BoundaryProblem
+    formulation :: Symbol # Dual or Standard
+    inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
+    normal_condition :: Symbol # Tie or Contact
+    tangential_condition :: Symbol # Stick or Slip
+    minimum_distance :: Float64 # don't check for a contact if elements are far enough
+    store_debug_info :: Bool # for making debugging easier
+    always_in_contact :: Vector{Int64} # nodes in this list always in contact
+    always_in_stick :: Vector{Int64} # nodes in this list always in stick
+    always_in_slip :: Vector{Int64} # nodes in this list always in slip
+    contact :: Bool
+    friction :: Bool
+end
+
+function Mortar()
+    Mortar(:Dual, false, :Tie, :Stick, Inf, false, [], [], [], false, false)
+end
+
+function get_unknown_field_name(::Type{Mortar})
+    return "reaction force"
+end
 
 macro debug(msg)
     haskey(ENV, "DEBUG") || return
     return msg
 end
 
+include("mortar_2d.jl")
 
-""" 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{S,M}(slave::Element{S}, master::Element{M}, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
-
-    # slave side geometry and normal direction at xi1
-    X1 = slave("geometry", xi1, time)
-    N1 = slave("normal-tangential coordinates", xi1, time)[:,1]
-
-    # master side geometry at xi2
-    master_basis(xi2) = get_basis(M, [xi2])
-    master_dbasis(xi2) = get_dbasis(M, [xi2])
-    master_geometry = master("geometry")(time)
-
-    function X2(xi2)
-        N = master_basis(xi2)
-        return sum([N[i]*master_geometry[i] for i=1:length(N)])
-    end
-
-    function dX2(xi2)
-        dN = master_dbasis(xi2)
-        return sum([dN[i]*master_geometry[i] for i=1:length(dN)])
-    end
-
-    # equation to solve
-    R(xi2) = det([X2(xi2)-X1  N1]')
-    dR(xi2) = det([dX2(xi2) N1]')
-
-    # solve using Newton iterations
-    xi2 = 0.0
-    for i=1:max_iterations
-        dxi2 = -R(xi2) / dR(xi2)
-        xi2 += dxi2
-        if norm(dxi2) < tol
-            return Float64[xi2]
-        end
-    end
-
-    println("slave element geometry")
-    dump(slave("geometry", time).data)
-    println("master element geometry")
-    dump(master("geometry", time).data)
-    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{S,M}(slave::Element{S}, master::Element{M}, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
-#   slave_basis = get_basis(slave)
-
-    # slave side geometry and normal direction at xi1
-
-    slave_geometry = slave("geometry")(time)
-    slave_normals = slave("normal-tangential coordinates")(time)
-    #slave_basis = slave.basis.data.basis
-    #slave_dbasis = slave.basis.data.dbasis
-    slave_basis(xi) = get_basis(S, [xi])
-    slave_dbasis(xi) = get_dbasis(S, [xi])
-
-    function X1(xi1)
-        N = slave_basis(xi1)
-        return sum([N[i]*slave_geometry[i] for i=1:length(N)])
-    end
-
-    function dX1(xi1)
-        dN = slave_dbasis(xi1)
-        return sum([dN[i]*slave_geometry[i] for i=1:length(dN)])
-    end
-
-    function N1(xi1)
-        N = slave_basis(xi1)
-        return sum([N[i]*slave_normals[i] for i=1:length(N)])[:,1]
-    end
-
-    function dN1(xi1)
-        dN = slave_dbasis(xi1)
-        return sum([dN[i]*slave_normals[i] for i=1:length(dN)])[:,1]
-    end
-
-    #X1(xi1) = slave_basis("geometry", [xi1], time)
-    #N1(xi1) = slave_basis("normal-tangential coordinates", [xi1], time)[:,1]
-
-    #master_basis = get_basis(master)
-
-    # master side geometry at xi2
-    #X2 = master_basis("geometry", xi2, time)
-    X2 = master("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)]')
-
-    #=
-    info("R(-1.0) = $(R(-1.0))")
-    info("R( 0.0) = $(R(0.0))")
-    info("R( 1.0) = $(R(1.0))")
-    info("R( 1.5) = $(R(1.5))")
-    info("dR(-1.0) = $(dR(-1.0))")
-    info("dR( 0.0) = $(dR(0.0))")
-    info("dR( 1.0) = $(dR(1.0))")
-    info("dR( 1.5) = $(dR(1.5))")
-    =#
-
-    #dR = ForwardDiff.derivative(R)
-
-    # go!
-    xi1 = 0.0
-    for i=1:max_iterations
-        dxi1 = -R(xi1) / dR(xi1)
-        xi1 += dxi1
-        #info("dxi1 = $dxi1, xi1 = $xi1, norm(dxi1) = $(norm(dxi1))")
-        if norm(dxi1) < tol
-            return Float64[xi1]
-        end
-    end
-
-    println("slave element geometry")
-    dump(slave("geometry", time).data)
-    println("master element geometry")
-    dump(master("geometry", time).data)
-    error("find projection from master to slave: did not converge")
-end
 
 ### Mortar projection calculation for 3d cases
 
@@ -640,274 +555,6 @@ function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix
     error("project_point_to_surface: did not converge in $max_iterations iterations!")
 end
 
-### Mortar problem
-# abstract MortarProblem{T} <: AbstractProblem
-# Mortar assembly 2d
-
-"""
-Currently two strategies exists:
-
-a) Remove inactive inequality constraints in element level. This is done in
-   assemble! if normal_condition is set to :Contact. For some reason this
-   leads to convergence issues.
-b) Remove inactive inequality constraints in assembly level. This is done in
-   posthook algorithm if inequality_constraints is set to true. This gives
-   more robust behavior.
-
-   Either use inequality_constraints=True OR :Contact + :Slip, but do not mix.
-
-   minimum_distance can be used to roughly skip integration of mortar
-   projections for elements that are "far enough" from each other. Increases
-   performance.
-
-"""
-type Mortar <: BoundaryProblem
-    formulation :: Symbol # Dual or Standard
-    inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
-    normal_condition :: Symbol # Tie or Contact
-    tangential_condition :: Symbol # Stick or Slip
-    minimum_distance :: Float64 # don't check for a contact if elements are far enough
-    store_debug_info :: Bool # for making debugging easier
-    always_in_contact :: Vector{Int64} # nodes in this list always in contact
-    always_in_stick :: Vector{Int64} # nodes in this list always in stick
-    always_in_slip :: Vector{Int64} # nodes in this list always in slip
-    contact :: Bool
-    friction :: Bool
-end
-
-function Mortar()
-    Mortar(:Dual, false, :Tie, :Stick, Inf, false, [], [], [], false, false)
-end
-
-function get_unknown_field_name(::Type{Mortar})
-    return "reaction force"
-end
-        
-# quadratic not tested yet
-typealias MortarElements2D Union{Seg2}
-
-function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar},
-                                        slave_element::Element{E}, time::Real)
-
-    # slave element must have a set of master elements
-    haskey(slave_element, "master elements") || return
-
-    # standard formulation for contact is not working at the moment
-    props = problem.properties
-    if props.formulation == :Standard && props.normal_condition == :Contact
-        error("for contact choose Dual formulation.""")
-    end
-
-    # get dimension and name of PARENT field
-    field_dim = problem.dimension
-    field_name = problem.parent_field_name
-
-    slave_dofs = get_gdofs(slave_element, field_dim)
-    nnodes = size(slave_element, 2)
-
-    # slave side quantities: rotation matrix, geometry, displacement, reaction force
-    Q = slave_element("normal-tangential coordinates", time)
-    Z = zeros(nnodes, nnodes)
-    if nnodes == 2
-        Q2 = [Q[1] Z; Z Q[2]]
-    elseif nnodes == 3
-        Q2 = [Q[1] Z Z; Z Q[2] Z; Z Z Q[3]]
-    end
-    X1 = vec(slave_element("geometry", time))
-    u1 = zeros(2*nnodes)
-    if haskey(slave_element, "displacement")
-        u1 = vec(slave_element("displacement", time))
-    end
-    x1 = X1 + u1
-    la = zeros(2*nnodes)
-    if haskey(slave_element, "reaction force")
-        la = vec(slave_element("reaction force", time))
-    end
-    la = Q2'*la
-
-    G = zeros(2*nnodes)
-    u = zeros(2*nnodes)
-    c = zeros(2*nnodes)
-
-    local_assembly = Assembly()
-
-    for master_element in slave_element["master elements"]
-
-        X2 = vec(master_element("geometry", time))
-        u2 = zeros(2*nnodes)
-        if haskey(master_element, "displacement")
-            u2 = vec(master_element("displacement", time))
-        end
-        x2 = X2 + u2
-
-        # if distance between elements is "far enough" cannot expect contact
-        if props.contact && (props.minimum_distance < Inf)
-            slave_midpoint = Float64[mean(x1[1:field_dim:2]), mean(x1[2:field_dim:2])]
-            master_midpoint = Float64[mean(x2[1:field_dim:2]), mean(x2[2:field_dim:2])]
-            if norm(slave_midpoint - master_midpoint) > props.minimum_distance
-                continue
-            end
-        end
-
-        master_dofs = get_gdofs(master_element, field_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*abs(xi1[2]-xi1[1])
-        isapprox(l, 0.0) && continue # no contribution
-
-        # Calculate slave side projection matrix D
-        Ae = zeros(nnodes, nnodes)
-        De = zeros(nnodes, nnodes)
-        Me = zeros(nnodes, nnodes)
-        if problem.properties.formulation == :Dual # Construct dual basis
-            for ip in get_integration_points(slave_element, Val{5})
-                J = get_jacobian(slave_element, ip, time)
-                w = ip.weight*norm(J)*l
-                xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
-                N = slave_element(xi, time)
-                De += w*diagm(vec(N))
-                Me += w*N'*N
-            end
-            Ae = De*inv(Me)
-        else # Standard Lagrange basis
-            for ip in get_integration_points(slave_element, Val{5})
-                J = get_jacobian(slave_element, ip, time)
-                w = ip.weight*norm(J)*l
-                xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
-                N = slave_element(xi, time)
-                De += w*N'*N
-            end
-            Ae = eye(nnodes)
-        end
-
-        C1S2 = zeros(2*nnodes, 2*nnodes)
-        C1M2 = zeros(2*nnodes, 2*nnodes)
-
-        # Slave side already done; it's De
-        for i=1:field_dim
-            C1S2[i:field_dim:end,i:field_dim:end] += De
-        end
-
-        # Calculate master side projection matrix M
-        for ip in get_integration_points(slave_element, Val{5})
-            J = get_jacobian(slave_element, ip, time)
-            w = ip.weight*norm(J)*l
-            # integration point on slave side segment
-            xi_slave = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
-            # projected integration point to master side element
-            xi_master = project_from_slave_to_master(slave_element, master_element, xi_slave)
-            N1 = slave_element(xi_slave, time)
-            N2 = master_element(xi_master, time)
-            M = w*kron(Ae*N1', N2)
-            for i=1:field_dim
-                C1M2[i:field_dim:end,i:field_dim:end] += M
-            end
-        end
-
-        # Calculate normal-tangential constraints
-        C2S2 = Q2'*C1S2
-        C2M2 = Q2'*C1M2
-
-        # initial weighted gap (capital G for "undeformed")
-        G += -(C2S2*X1 - C2M2*X2)
-        # change in weighted gap caused by deformation
-        u += -(C2S2*u1 - C2M2*u2)
-
-        # Add contributions
-        add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2)
-        add!(local_assembly.C1, slave_dofs, master_dofs, -C1M2)
-        add!(local_assembly.C2, slave_dofs, slave_dofs, C2S2)
-        add!(local_assembly.C2, slave_dofs, master_dofs, -C2M2)
-
-    end # all master elements are done
-
-    add!(local_assembly.g, slave_dofs, G)
-
-    # if only equality constraints, i.e., mesh tying problem, we're done for this element.
-    if !props.contact
-        append!(assembly, local_assembly)
-        return
-    end
-
-    C1 = sparse(local_assembly.C1)
-    C2 = sparse(local_assembly.C2)
-    D = spzeros(size(C2)...)
-    g = sparse(local_assembly.g)
-
-    # complementarity condition
-    lan = la[1:field_dim:end]
-    lat = la[2:field_dim:end]
-    Gn = G[1:field_dim:end]
-    un = u[1:field_dim:end]
-    cn = lan - (Gn + un)
-    #Cn = lan - max(0, lan - (Gn+un))
-
-    # normal condition
-    inactive_nodes = find(cn .<= 0)
-    active_nodes = find(cn .> 0)
-    #inactive_nodes = find(Cn .>= 0)
-    #active_nodes = find(Cn .== 0)
-
-    # inactive element
-    if length(active_nodes) == 0
-        return
-    end
-
-    node_ids = get_connectivity(slave_element)
-
-    # normal constraint: remove inactive nodes
-    for j in node_ids[inactive_nodes]
-        if length(props.always_in_contact) != 0
-            j in props.always_in_contact && continue
-        end
-        gdofs = [2*(j-1)+1, 2*(j-1)+2]
-        C1[gdofs,:] = 0
-        C2[gdofs,:] = 0
-        D[gdofs,:] = 0
-        g[gdofs] = 0
-    end
-
-    # frictional contact, see Gitterle2010
-    mu = 0.3
-    ct = lat + c[2:field_dim:end]
-
-    C = max(mu*cn, abs(ct)).*lat - mu*max(0, cn).*ct
-    stick_nodes = find(abs(ct) - mu*cn .< 0)
-    slip_nodes = find(abs(ct) - mu*cn .>= 0)
-    stick_nodes = setdiff(stick_nodes, inactive_nodes)
-    slip_nodes = setdiff(slip_nodes, inactive_nodes)
-
-    for (i, j) in enumerate(node_ids[active_nodes])
-        gdofs = [2*(j-1)+1, 2*(j-1)+2]
-        #D[gdofs[2],gdofs] = C2[gdofs[2],gdofs]
-        D[gdofs[2],gdofs] = Q[i][:,2]'
-        C2[gdofs[2],:] = 0
-        if props.friction
-            g[gdofs[2]] = C[i]
-        else
-            g[gdofs[2]] = 0.0
-        end
-    end
-
-    local_assembly.C1 = C1
-    local_assembly.C2 = C2
-    local_assembly.D = D
-    local_assembly.g = g
-    append!(assembly, local_assembly)
-
-    if props.store_debug_info
-        slave_element["G"] = G
-        slave_element["g"] = g
-        slave_element["c"] = c
-        slave_element["C1"] = C1
-        slave_element["C2"] = C2
-        slave_element["D"] = D2
-        slave_element["active nodes"] = active_nodes
-    end
-
-end
-
 typealias MortarElements3D Union{Tri3, Quad4}
 
 function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar},
diff --git a/src/mortar_2d.jl b/src/mortar_2d.jl
new file mode 100644
index 0000000..caab0cf
--- /dev/null
+++ b/src/mortar_2d.jl
@@ -0,0 +1,329 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+# Mortar projection calculation for 2d
+
+""" 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{S,M}(slave::Element{S}, master::Element{M}, xi1::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
+
+    # slave side geometry and normal direction at xi1
+    X1 = slave("geometry", xi1, time)
+    N1 = slave("normal-tangential coordinates", xi1, time)[:,1]
+
+    # master side geometry at xi2
+    master_basis(xi2) = get_basis(M, [xi2])
+    master_dbasis(xi2) = get_dbasis(M, [xi2])
+    master_geometry = master("geometry")(time)
+
+    function X2(xi2)
+        N = master_basis(xi2)
+        return sum([N[i]*master_geometry[i] for i=1:length(N)])
+    end
+
+    function dX2(xi2)
+        dN = master_dbasis(xi2)
+        return sum([dN[i]*master_geometry[i] for i=1:length(dN)])
+    end
+
+    # equation to solve
+    R(xi2) = det([X2(xi2)-X1  N1]')
+    dR(xi2) = det([dX2(xi2) N1]')
+
+    # solve using Newton iterations
+    xi2 = 0.0
+    for i=1:max_iterations
+        dxi2 = -R(xi2) / dR(xi2)
+        xi2 += dxi2
+        if norm(dxi2) < tol
+            return Float64[xi2]
+        end
+    end
+
+    println("slave element geometry")
+    dump(slave("geometry", time).data)
+    println("master element geometry")
+    dump(master("geometry", time).data)
+    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{S,M}(slave::Element{S}, master::Element{M}, xi2::Vector, time::Float64=0.0; max_iterations=5, tol=1.0e-9)
+
+    # slave side geometry and normal direction at xi1
+
+    slave_geometry = slave("geometry")(time)
+    slave_normals = slave("normal-tangential coordinates")(time)
+    slave_basis(xi) = get_basis(S, [xi])
+    slave_dbasis(xi) = get_dbasis(S, [xi])
+
+    function X1(xi1)
+        N = slave_basis(xi1)
+        return sum([N[i]*slave_geometry[i] for i=1:length(N)])
+    end
+
+    function dX1(xi1)
+        dN = slave_dbasis(xi1)
+        return sum([dN[i]*slave_geometry[i] for i=1:length(dN)])
+    end
+
+    function N1(xi1)
+        N = slave_basis(xi1)
+        return sum([N[i]*slave_normals[i] for i=1:length(N)])[:,1]
+    end
+
+    function dN1(xi1)
+        dN = slave_dbasis(xi1)
+        return sum([dN[i]*slave_normals[i] for i=1:length(dN)])[:,1]
+    end
+
+    # master side geometry at xi2
+    X2 = master("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)]')
+
+    # go!
+    xi1 = 0.0
+    for i=1:max_iterations
+        dxi1 = -R(xi1) / dR(xi1)
+        xi1 += dxi1
+        #info("dxi1 = $dxi1, xi1 = $xi1, norm(dxi1) = $(norm(dxi1))")
+        if norm(dxi1) < tol
+            return Float64[xi1]
+        end
+    end
+
+    println("slave element geometry")
+    dump(slave("geometry", time).data)
+    println("master element geometry")
+    dump(master("geometry", time).data)
+    error("find projection from master to slave: did not converge")
+end
+
+
+# Mortar assembly 2d
+
+# quadratic not tested yet
+typealias MortarElements2D Union{Seg2}
+
+function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar},
+                                        slave_element::Element{E}, time::Real)
+
+    # slave element must have a set of master elements
+    haskey(slave_element, "master elements") || return
+    props = problem.properties
+
+    # get dimension and name of PARENT field
+    field_dim = problem.dimension
+    field_name = problem.parent_field_name
+
+    slave_dofs = get_gdofs(slave_element, field_dim)
+    nnodes = size(slave_element, 2)
+
+    # slave side quantities: rotation matrix, geometry, displacement, reaction force
+    Q = slave_element("normal-tangential coordinates", time)
+    Z = zeros(nnodes, nnodes)
+    if nnodes == 2
+        Q2 = [Q[1] Z; Z Q[2]]
+    elseif nnodes == 3
+        Q2 = [Q[1] Z Z; Z Q[2] Z; Z Z Q[3]]
+    end
+    X1 = vec(slave_element("geometry", time))
+    u1 = zeros(2*nnodes)
+    if haskey(slave_element, "displacement")
+        u1 = vec(slave_element("displacement", time))
+    end
+    x1 = X1 + u1
+    la = zeros(2*nnodes)
+    if haskey(slave_element, "reaction force")
+        la = vec(slave_element("reaction force", time))
+    end
+    la = Q2'*la
+
+    G = zeros(2*nnodes)
+    u = zeros(2*nnodes)
+    c = zeros(2*nnodes)
+
+    local_assembly = Assembly()
+
+    for master_element in slave_element["master elements"]
+
+        X2 = vec(master_element("geometry", time))
+        u2 = zeros(2*nnodes)
+        if haskey(master_element, "displacement")
+            u2 = vec(master_element("displacement", time))
+        end
+        x2 = X2 + u2
+
+        # if distance between elements is "far enough" cannot expect contact
+        if props.contact && (props.minimum_distance < Inf)
+            slave_midpoint = Float64[mean(x1[1:field_dim:2]), mean(x1[2:field_dim:2])]
+            master_midpoint = Float64[mean(x2[1:field_dim:2]), mean(x2[2:field_dim:2])]
+            if norm(slave_midpoint - master_midpoint) > props.minimum_distance
+                continue
+            end
+        end
+
+        master_dofs = get_gdofs(master_element, field_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*abs(xi1[2]-xi1[1])
+        isapprox(l, 0.0) && continue # no contribution
+
+        # Calculate slave side projection matrix D
+        Ae = zeros(nnodes, nnodes)
+        De = zeros(nnodes, nnodes)
+        Me = zeros(nnodes, nnodes)
+        if problem.properties.formulation == :Dual # Construct dual basis
+            for ip in get_integration_points(slave_element, Val{5})
+                J = get_jacobian(slave_element, ip, time)
+                w = ip.weight*norm(J)*l
+                xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
+                N = slave_element(xi, time)
+                De += w*diagm(vec(N))
+                Me += w*N'*N
+            end
+            Ae = De*inv(Me)
+        else # Standard Lagrange basis
+            for ip in get_integration_points(slave_element, Val{5})
+                J = get_jacobian(slave_element, ip, time)
+                w = ip.weight*norm(J)*l
+                xi = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
+                N = slave_element(xi, time)
+                De += w*N'*N
+            end
+            Ae = eye(nnodes)
+        end
+
+        C1S2 = zeros(2*nnodes, 2*nnodes)
+        C1M2 = zeros(2*nnodes, 2*nnodes)
+
+        # Slave side already done; it's De
+        for i=1:field_dim
+            C1S2[i:field_dim:end,i:field_dim:end] += De
+        end
+
+        # Calculate master side projection matrix M
+        for ip in get_integration_points(slave_element, Val{5})
+            J = get_jacobian(slave_element, ip, time)
+            w = ip.weight*norm(J)*l
+            # integration point on slave side segment
+            xi_slave = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
+            # projected integration point to master side element
+            xi_master = project_from_slave_to_master(slave_element, master_element, xi_slave)
+            N1 = slave_element(xi_slave, time)
+            N2 = master_element(xi_master, time)
+            M = w*kron(Ae*N1', N2)
+            for i=1:field_dim
+                C1M2[i:field_dim:end,i:field_dim:end] += M
+            end
+        end
+
+        # Calculate normal-tangential constraints
+        C2S2 = Q2'*C1S2
+        C2M2 = Q2'*C1M2
+
+        # initial weighted gap (capital G for "undeformed")
+        G += -(C2S2*X1 - C2M2*X2)
+        # change in weighted gap caused by deformation
+        u += -(C2S2*u1 - C2M2*u2)
+
+        # Add contributions
+        add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2)
+        add!(local_assembly.C1, slave_dofs, master_dofs, -C1M2)
+        add!(local_assembly.C2, slave_dofs, slave_dofs, C2S2)
+        add!(local_assembly.C2, slave_dofs, master_dofs, -C2M2)
+
+    end # all master elements are done
+
+    add!(local_assembly.g, slave_dofs, G)
+
+    # if only equality constraints, i.e., mesh tying problem, we're done for this element.
+    if !props.contact
+        append!(assembly, local_assembly)
+        return
+    end
+
+    C1 = sparse(local_assembly.C1)
+    C2 = sparse(local_assembly.C2)
+    D = spzeros(size(C2)...)
+    g = sparse(local_assembly.g)
+
+    # complementarity condition
+    lan = la[1:field_dim:end]
+    lat = la[2:field_dim:end]
+    Gn = G[1:field_dim:end]
+    un = u[1:field_dim:end]
+    cn = lan - (Gn + un)
+    #Cn = lan - max(0, lan - (Gn+un))
+
+    # normal condition
+    inactive_nodes = find(cn .<= 0)
+    active_nodes = find(cn .> 0)
+    #inactive_nodes = find(Cn .>= 0)
+    #active_nodes = find(Cn .== 0)
+
+    # inactive element
+    if length(active_nodes) == 0
+        return
+    end
+
+    node_ids = get_connectivity(slave_element)
+
+    # normal constraint: remove inactive nodes
+    for j in node_ids[inactive_nodes]
+        if length(props.always_in_contact) != 0
+            j in props.always_in_contact && continue
+        end
+        gdofs = [2*(j-1)+1, 2*(j-1)+2]
+        C1[gdofs,:] = 0
+        C2[gdofs,:] = 0
+        D[gdofs,:] = 0
+        g[gdofs] = 0
+    end
+
+    # frictional contact, see Gitterle2010
+    mu = 0.3
+    ct = lat + c[2:field_dim:end]
+
+    C = max(mu*cn, abs(ct)).*lat - mu*max(0, cn).*ct
+    stick_nodes = find(abs(ct) - mu*cn .< 0)
+    slip_nodes = find(abs(ct) - mu*cn .>= 0)
+    stick_nodes = setdiff(stick_nodes, inactive_nodes)
+    slip_nodes = setdiff(slip_nodes, inactive_nodes)
+
+    for (i, j) in enumerate(node_ids[active_nodes])
+        gdofs = [2*(j-1)+1, 2*(j-1)+2]
+        #D[gdofs[2],gdofs] = C2[gdofs[2],gdofs]
+        D[gdofs[2],gdofs] = Q[i][:,2]'
+        C2[gdofs[2],:] = 0
+        if props.friction
+            g[gdofs[2]] = C[i]
+        else
+            g[gdofs[2]] = 0.0
+        end
+    end
+
+    local_assembly.C1 = C1
+    local_assembly.C2 = C2
+    local_assembly.D = D
+    local_assembly.g = g
+    append!(assembly, local_assembly)
+
+    if props.store_debug_info
+        slave_element["G"] = G
+        slave_element["g"] = g
+        slave_element["c"] = c
+        slave_element["C1"] = C1
+        slave_element["C2"] = C2
+        slave_element["D"] = D2
+        slave_element["active nodes"] = active_nodes
+    end
+
+end
+