diff --git a/docs/links.md b/docs/links.md
index 06b175e..2063fb9 100644
--- a/docs/links.md
+++ b/docs/links.md
@@ -1,12 +1,13 @@
-Discretization
-==============
+Links
+=====
 
+Discretization
+--------------
 http://code.activestate.com/recipes/579021-delaunay-triangulation/
 
 
 Interpolation
-=============
-
+-------------
 - http://www.cs.rpi.edu/~flaherje/pdf/fea4.pdf
 - http://www.sd.ruhr-uni-bochum.de/downloads/Shape_funct.pdf
 - http://what-when-how.com/the-finite-element-method/fem-for-3d-solids-finite-element-method-part-1/
@@ -18,10 +19,13 @@ Interpolation
 - http://www.researchgate.net/publication/267082822_Unified_isoparametric_3D_Lagrange_finite_elements
 
 Integration
-===========
+-----------
 - http://arxiv.org/pdf/1411.1341.pdf
 
 Hierarchial shape functions
+---------------------------
 - https://www.math.vt.edu/people/adjerids/research/papers/basis.pdf
 
-- edited by Ari
\ No newline at end of file
+Solvers
+-------
+- https://github.com/JuliaSparse/MultiFrontalCholesky.jl
diff --git a/src/assembly.jl b/src/assembly.jl
index 33e2b3f..57838e4 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -61,7 +61,9 @@ function assemble(problem::AllProblems, time::Real, nchunks=10)
     for (j, elrange) in enumerate(slices)
         sub_assembly = assemble(problem, elrange, time)
         append!(assembly, sub_assembly)
-        info("Assembly: ", round(j/nchunks*100,1), " % done. ")
+        if ne > 100
+            info("Assembly: ", round(j/nchunks*100,1), " % done. ")
+        end
     end
 #   optimize!(assembly)
 #   dim = length(assembly.stiffness_matrix.I)
diff --git a/src/elements.jl b/src/elements.jl
index 2cff190..2ea6884 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -210,3 +210,26 @@ function Base.haskey(element::Element, what)
     haskey(element.fields, what)
 end
 
+""" Calculate local normal-tangential coordinates for element. """
+function calculate_normal_tangential_coordinates!{E}(element::Element{E}, time::Real)
+    proj(u, v) = dot(v, u) / dot(u, u) * u
+    ntcoords = Matrix[]
+    refcoords = get_reference_element_coordinates(E)
+    x = element("geometry", time)
+    for xi in refcoords
+        dN = get_dbasis(E, xi)*x
+        normal = cross(dN[:,1], dN[:,2])
+        normal /= norm(normal)
+        u1 = normal
+        j = indmax(abs(u1))
+        v2 = zeros(3)
+        v2[mod(j,3)+1] = 1.0
+        u2 = v2 - proj(u1, v2)
+        u3 = cross(u1, u2)
+        tangent1 = u2/norm(u2)
+        tangent2 = u3/norm(u3)
+        push!(ntcoords, [normal tangent1 tangent2])
+    end
+    element["normal-tangential coordinates"] = ntcoords
+end
+
diff --git a/src/lagrange.jl b/src/lagrange.jl
index fcac227..02dd9e7 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -27,6 +27,10 @@ function calculate_lagrange_basis_coefficients(P, X)
     return inv(A)'
 end
 
+function refcoords(X::Matrix)
+    return Vector{Float64}[X[:,i] for i=1:size(X,2)]
+end
+
 """
 Create new Lagrange element
 
@@ -37,7 +41,10 @@ Examples
 macro create_lagrange_element(element_name, element_description, X, P)
     eltype = esc(element_name)
     quote
-        global get_basis, get_dbasis
+        global get_basis, get_dbasis,
+               get_reference_element_coordinates,
+               get_reference_element_midpoint
+
         #basis, dbasis = calculate_lagrange_basis($P, $X)
         C = calculate_lagrange_basis_coefficients($P, $X)
         basis(xi) = C*$P(xi)
@@ -49,6 +56,15 @@ macro create_lagrange_element(element_name, element_description, X, P)
             return basis(xi)'
         end
 
+        XX = refcoords($X)
+        function get_reference_element_coordinates(::Type{$eltype})
+            return XX
+        end
+
+        XXX = vec(mean($X, 2))
+        function get_reference_element_midpoint(::Type{$eltype})
+            return XXX
+        end
 #=
         function get_dbasis(::Type{$eltype}, xi::Vector{Float64})
             return dbasis(xi)'
@@ -81,6 +97,10 @@ end
      0.0 0.0 1.0],
     (xi) -> [1.0, xi[1], xi[2]])
 
+#function get_reference_element_midpoint(::Type{Tri3})
+#    return [1.0/3.0, 1.0/3.0]
+#end
+
 @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],
diff --git a/src/mortar.jl b/src/mortar.jl
index 1713d88..934cb60 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -11,7 +11,7 @@ function project_from_slave_to_master{S,M}(slave::Element{S}, master::Element{M}
 
     # slave side geometry and normal direction at xi1
     X1 = slave("geometry", xi1, time)
-    N1 = slave("nodal ntsys", xi1, time)[:,1]
+    N1 = slave("normal-tangential coordinates", xi1, time)[:,1]
 
     # master side geometry at xi2
     #master_basis = master.basis.data.basis
@@ -59,7 +59,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
     # slave side geometry and normal direction at xi1
 
     slave_geometry = slave("geometry")(time)
-    slave_normals = slave("nodal ntsys")(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])
@@ -86,7 +86,7 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
     end
 
     #X1(xi1) = slave_basis("geometry", [xi1], time)
-    #N1(xi1) = slave_basis("nodal ntsys", [xi1], time)[:,1]
+    #N1(xi1) = slave_basis("normal-tangential coordinates", [xi1], time)[:,1]
 
     #master_basis = get_basis(master)
 
@@ -187,11 +187,16 @@ Notes
 
 """
 # function create_auxiliary_plane(x, ximp, normals, basis)
-function create_auxiliary_plane(element::Element{Tri3}, time::Real)
-    proj(u, v) = dot(v, u) / dot(u, u) * u
-    xi = [1.0/3.0, 1.0/3.0]
+function create_auxiliary_plane{E}(element::Element{E}, time::Real)
+#   proj(u, v) = dot(v, u) / dot(u, u) * u
+#    xi = [1.0/3.0, 1.0/3.0]
+
+    xi = get_reference_element_midpoint(E)
     x0 = element("geometry", xi, time)
-    n = element("nodal ntsys", xi, time)[:, 1]
+    ntbasis = element("normal-tangential coordinates", xi, time)
+    return x0, ntbasis
+#=
+    n = element("normal-tangential coordinates", xi, time)[:, 1]
     n /= norm(n)
     # gram-schmidt
     u1 = n
@@ -204,6 +209,7 @@ function create_auxiliary_plane(element::Element{Tri3}, time::Real)
     t2 = u3/norm(u3)
     new_basis = [n t1 t2]
     return x0, new_basis
+=#
 end
 
 """
@@ -243,7 +249,14 @@ function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix)
     n = Q[:,1]
     ph = p - dot(p-x0, n)*n
     qproj = Q'*(ph-x0)
-    @assert isapprox(qproj[1], 0.0)
+    if !isapprox(qproj[1], 0.0; atol=1.0e-12)
+        info("project_point_to_auxiliary_plane(): point not projected correctly.")
+        info("p: $p")
+        info("x0: $x0")
+        info("Q: \n$Q")
+        info("qproj: $qproj")
+        error("Failed to project point to auxiliary plane.")
+    end
     return qproj[2:3]
 end
 
@@ -626,6 +639,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::BoundaryPro
             N2 = master_element(xi_projected, time)
             S = w*N1'*N1
             M = w*(N1'*N2)'
+#           M = w*N1'*N2
             # FIXME: why this needs now to be transpose?
             # assembly / repeat
             for i=1:field_dim
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 2c8d2f9..6930d39 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -6,7 +6,8 @@ module ElementTests
 using JuliaFEM.Test
 
 using JuliaFEM.Core: AbstractElement, Element, Field, FieldSet, test_element
-import JuliaFEM.Core: get_basis, get_dbasis
+using JuliaFEM.Core: Tri3
+import JuliaFEM.Core: get_basis, get_dbasis, calculate_normal_tangential_coordinates!
 import Base: size
 
 """ Prototype element
@@ -75,4 +76,19 @@ function test_interpolate()
 #   @test isapprox(gradT, 1/2*gradT_expected)
 end
 
+function test_calculate_normal_tangential_coordinates()
+    el = Tri3([1, 2, 3])
+    el["geometry"] = Vector{Float64}[
+        [0.0, 0.0, 0.0],
+        [1.0, 0.0, 0.0],
+        [0.0, 1.0, 0.0]]
+    calculate_normal_tangential_coordinates!(el, 0.0)
+    n = [0.0 0.0 1.0]'
+    t1 = [1.0 0.0 0.0]'
+    t2 = [0.0 1.0 0.0]'
+    R = [n t1 t2]
+    @test isapprox(el("normal-tangential coordinates", [0.0, 0.0], 0.0), R)
+end
+#test_calculate_normal_tangential_coordinates()
+
 end
diff --git a/test/test_mortar.jl b/test/test_mortar.jl
index 8a2eb40..4c9f493 100644
--- a/test/test_mortar.jl
+++ b/test/test_mortar.jl
@@ -15,7 +15,8 @@ using JuliaFEM.Core: project_from_slave_to_master, project_from_master_to_slave
 using JuliaFEM.Core: create_auxiliary_plane, project_point_to_auxiliary_plane,
                      get_edge_intersections, get_points_inside_triangle,
                      clip_polygon, calculate_polygon_centerpoint,
-                     project_point_from_plane_to_surface, assemble
+                     project_point_from_plane_to_surface, assemble,
+                     calculate_normal_tangential_coordinates!
 
 
 function get_test_2d_model()
@@ -42,13 +43,13 @@ function get_test_2d_model()
     slave1 = Seg2([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)]
+    slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
     slave1["master elements"] = Element[master1, master2]
 
     slave2 = Seg2([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)]
+    slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
     slave2["master elements"] = Element[master1, master2]
 
     return [slave1, slave2], [master1, master2]
@@ -88,7 +89,7 @@ function test_calc_flat_2d_projection_rotated()
     master1["geometry"] = Vector{Float64}[[0.0, 1.0], [0.0, 0.0]]
     slave1 = Seg2([1, 2])
     slave1["geometry"] = Vector{Float64}[[0.0, 0.0], [0.0, 1.0]]
-    slave1["nodal ntsys"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]]
+    slave1["normal-tangential coordinates"] = Matrix{Float64}[[1.0 0.0; 0.0 1.0], [1.0 0.0; 0.0 1.0]]
     xi = project_from_master_to_slave(slave1, master1, [-1.0])
     info("xi = $xi")
     @test xi == [ 1.0]
@@ -209,7 +210,7 @@ function test_2d_mortar_multiple_bodies_multiple_dirichlet_bc()
 
     slave1 = Seg2([5, 6])
     slave1["geometry"] = Vector[N[5], N[6]]
-    slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
     slave1["master elements"] = Element[master1]
 
     boundary3 = MortarProblem("displacement", 2)
@@ -237,13 +238,19 @@ end
 
 
 function test_2d_mortar_three_bodies_shared_nodes()
-    N = Vector[
-        [0.0, 0.0], [2.0, 0.0],
-        [0.0, 1.0], [2.0, 1.0],
-        [0.0, 1.0], [1.0, 1.0],
-        [0.0, 2.0], [1.0, 2.0],
-        [1.0, 1.0], [2.0, 1.0],
-        [1.0, 2.0], [2.0, 2.0]]
+    N = Dict{Int, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [2.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [2.0, 1.0],
+        5 => [0.0, 1.0],
+        6 => [1.3, 1.0],
+        7 => [0.0, 2.0],
+        8 => [1.3, 2.0],
+        9 => [1.3, 1.0],
+        10 => [2.0, 1.0],
+        11 => [1.3, 2.0],
+        12 => [2.0, 2.0])
 
     e1 = Quad4([1, 2, 4, 3])
     e1["geometry"] = Vector[N[1], N[2], N[4], N[3]]
@@ -307,7 +314,7 @@ function test_2d_mortar_three_bodies_shared_nodes()
 
     slave1 = Seg2([5, 6])
     slave1["geometry"] = Vector[N[5], N[6]]
-    slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    slave1["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
     slave1["master elements"] = Element[master1]
     bc3 = MortarProblem("displacement", 2)
     push!(bc3, slave1)
@@ -315,7 +322,7 @@ function test_2d_mortar_three_bodies_shared_nodes()
     # mortar boundary between body 1 and body 3
     slave2 = Seg2([9, 10])
     slave2["geometry"] = Vector[N[9], N[10]]
-    slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    slave2["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
     slave2["master elements"] = Element[master1]
     bc4 = MortarProblem("displacement", 2)
     push!(bc4, slave2)
@@ -326,8 +333,8 @@ function test_2d_mortar_three_bodies_shared_nodes()
 
     slave3 = Seg2([6, 8])
     slave3["geometry"] = Vector[N[6], N[8]]
-    #slave3["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
-    slave3["nodal ntsys"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
+    #slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
+    slave3["normal-tangential coordinates"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
     slave3["master elements"] = Element[master2]
     bc5 = MortarProblem("displacement", 2)
     push!(bc5, slave3)
@@ -346,12 +353,15 @@ function test_2d_mortar_three_bodies_shared_nodes()
 
     # launch solver
     solver.method = :UMFPACK
+    solver.name = "test_2d_mortar_three_bodies_shared_nodes"
+    solver.dump_matrices = true
     call(solver, 0.0)
 
-    disp = e2("displacement", [1.0, 1.0], 0.0)
-    info("displacement at tip: $disp")
+    X = e3("geometry", [1.0, 1.0], 0.0)
+    u = e3("displacement", [1.0, 1.0], 0.0)
+    info("displacement at $X: $u")
     # code aster verification, two_elements.comm
-    @test isapprox(disp, [3.17431158889468E-02, -2.77183037855653E-01])
+    @test isapprox(u, [2*3.17431158889468E-02, -2.77183037855653E-01])
 
 end
 #test_2d_mortar_three_bodies_shared_nodes()
@@ -368,7 +378,7 @@ function test_auxiliary_plane_transforms()
          0.0 0.0 1.0
          1.0 0.0 0.0]
     e1["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]
-    e1["nodal ntsys"] = Matrix{Float64}[R, R, R]
+    e1["normal-tangential coordinates"] = Matrix{Float64}[R, R, R]
     time::Real = 0.0
     x0, Q = create_auxiliary_plane(e1, time)
     info("x0 = $x0")
@@ -466,6 +476,7 @@ end
 #test_calculate_polygon_centerpoint()
 
 
+
 function test_assemble_3d_problem()
     nodes = Vector{Float64}[
         [0.0, 0.0, 0.0],
@@ -478,12 +489,14 @@ function test_assemble_3d_problem()
     mel["geometry"] = Vector{Float64}[nodes[4], nodes[5], nodes[6]]
     sel = Tri3([1, 2, 3])
     sel["geometry"] = Vector{Float64}[nodes[1], nodes[2], nodes[3]]
-    R = [0.0 1.0 0.0
-         0.0 0.0 1.0
-         1.0 0.0 0.0]
-    sel["nodal ntsys"] = Matrix{Float64}[R, R, R]
+#   Rv = [0.0 1.0 0.0
+#         0.0 0.0 1.0
+#         1.0 0.0 0.0]
+#   sel["normal-tangential coordinates"] = Matrix{Float64}[Rv, Rv, Rv]
+    calculate_normal_tangential_coordinates!(sel, 0.0)
     sel["master elements"] = Element[mel]
     prob = MortarProblem("temperature", 1)
+
     push!(prob, sel)
     stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
     info("stiffness matrix for this problem:\n$stiffness_matrix")
@@ -491,8 +504,36 @@ function test_assemble_3d_problem()
     B = [D -M]  # slave dofs are first in this.
     info("expected matrix for this problem:\n$B")
     @test isapprox(stiffness_matrix, B)
+
+    # rotate and translate surface and check that we are still having same results
+    Rx(t) = [
+        1.0 0.0     0.0
+        0.0 cos(t) -sin(t)
+        0.0 sin(t)  cos(t)]
+    Ry(t) = [
+        cos(t)  0.0 sin(t)
+        0.0     1.0 0.0
+        -sin(t) 0.0 cos(t)
+        ]
+    Rz(t) = [
+        cos(t) -sin(t) 0.0
+        sin(t)  cos(t) 0.0
+        0.0     0.0    1.0]
+    T = [1.0, 1.0, 1.0]
+    tx = pi/3.0
+    ty = pi/4.0
+    tz = pi/5.0
+    for node in nodes
+        node[:] = Rz(tz)*Ry(ty)*Rx(tx)*node + T
+    end
+    calculate_normal_tangential_coordinates!(sel, 0.0)
+    stiffness_matrix = full(assemble(prob, 0.0).stiffness_matrix)
+    info("sel midpnt: ", sel("geometry", [1/3, 1/3], 0.0))
+    info("nt basis: ", sel("normal-tangential coordinates", [1/3, 1/3], 0.0))
+    @test isapprox(stiffness_matrix, B)
+
 end
-#test_assemble_3d_problem()
+test_assemble_3d_problem()
 
 
 end