From 6d3e33c3ff8766ff08671d9d5e9bb1fe732780ad Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Wed, 22 Jun 2016 01:39:28 +0300
Subject: [PATCH] modal solver + tie contact works now. dirichlet boundary and
 mpcs are eliminated properly before solution to get reduced system

---
 src/JuliaFEM.jl                               |   5 +-
 src/integrate.jl                              |  55 ++-
 src/modal.jl                                  |  58 ++-
 src/mortar.jl                                 | 394 +++++++++++++++++-
 src/postprocess_utils.jl                      |  19 +
 src/preprocess_aster_reader.jl                |   1 +
 src/solvers.jl                                | 147 ++++---
 src/sparse.jl                                 |  11 +-
 ..._elasticity_2d_linear_with_surface_load.jl |   2 +-
 test/test_modal_analysis.jl                   |  76 ++++
 test/test_mortar_2d_mesh_tie.jl               | 119 ++++++
 test/test_mortar_3d_mesh_tie.jl               |  47 +++
 test/test_mortar_3d_polygon_clip.jl           |  41 ++
 test/test_projection.jl                       |  32 ++
 14 files changed, 887 insertions(+), 120 deletions(-)
 create mode 100644 test/test_mortar_2d_mesh_tie.jl
 create mode 100644 test/test_mortar_3d_mesh_tie.jl
 create mode 100644 test/test_mortar_3d_polygon_clip.jl
 create mode 100644 test/test_projection.jl

diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 68d021c..d489f8b 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -66,7 +66,8 @@ include("solvers.jl")
 export AbstractSolver, Solver, Nonlinear, 
        get_unknown_field_name, get_formulation_type,
        get_field_problems, get_boundary_problems,
-       get_field_assembly, get_boundary_assembly
+       get_field_assembly, get_boundary_assembly,
+       initialize!, create_projection
 include("modal.jl")
 export Modal
 
@@ -107,7 +108,7 @@ end
 
 module Postprocess
 include("postprocess_utils.jl")
-export calc_nodal_values!
+export calc_nodal_values!, get_nodal_vector
 include("postprocess_xdmf.jl")
 export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
 end
diff --git a/src/integrate.jl b/src/integrate.jl
index 2041925..e012a65 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -102,7 +102,11 @@ function get_integration_points(element::TriangularElement, ::Type{Val{2}})
 end
 
 function get_integration_points(element::TriangularElement, ::Type{Val{3}})
-    weights = 0.5*[-0.5625, 0.5208333333333333, 0.5208333333333333, 0.5208333333333333]
+    weights = 0.5*[
+        -0.5625,
+         0.5208333333333333,
+         0.5208333333333333,
+         0.5208333333333333]
     points = Vector{Float64}[
         [1.0/3.0, 1.0/3.0],
         [0.2, 0.2],
@@ -112,26 +116,41 @@ function get_integration_points(element::TriangularElement, ::Type{Val{3}})
 end
 
 function get_integration_points(element::TriangularElement, ::Type{Val{4}})
-    [
-        IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801),
-        IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801),
-        IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801),
-        IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532),
-        IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532),
-        IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532)
-    ]
+    weights = 0.5*[
+        0.22338158967801,
+        0.22338158967801,
+        0.22338158967801,
+        0.10995174365532,
+        0.10995174365532,
+        0.10995174365532]
+    points = Vector{Float64}[
+        [0.44594849091597, 0.44594849091597],        
+        [0.44594849091597, 0.10810301816807],
+        [0.10810301816807, 0.44594849091597],
+        [0.09157621350977, 0.09157621350977],
+        [0.09157621350977, 0.81684757298046],
+        [0.81684757298046, 0.09157621350977]]
+    return zip(weights, points)
 end
 
 function get_integration_points(element::TriangularElement, ::Type{Val{5}})
-    [
-        IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000),
-        IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851),
-        IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851),
-        IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851),
-        IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483),
-        IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483),
-        IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483)
-    ]
+    weights = 0.5*[
+        0.22500000000000,
+        0.13239415278851,
+        0.13239415278851,
+        0.13239415278851,
+        0.12593918054483,
+        0.12593918054483,
+        0.12593918054483]
+    points = Vector{Float64}[
+        [0.33333333333333, 0.33333333333333], 
+        [0.47014206410511, 0.47014206410511], 
+        [0.47014206410511, 0.05971587178977],
+        [0.05971587178977, 0.47014206410511],
+        [0.10128650732346, 0.10128650732346],
+        [0.10128650732346, 0.79742698535309],
+        [0.79742698535309, 0.10128650732346]]
+    return zip(weights, points)
 end
 
 ### 3d elements
diff --git a/src/modal.jl b/src/modal.jl
index 2464b5a..7e0eaf2 100644
--- a/src/modal.jl
+++ b/src/modal.jl
@@ -1,6 +1,16 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+"""
+Examples
+--------
+
+julia> problems = get_problems()
+julia> solver = Solver(Modal)
+julia> push!(solver, problems...)
+julia> call(solver)
+
+"""
 type Modal <: AbstractSolver
     geometric_stiffness :: Bool
     eigvals :: Vector
@@ -21,40 +31,50 @@ function call(solver::Solver{Modal}; debug=false)
         assemble!(problem, solver.time)
         assemble!(problem, solver.time, Val{:mass_matrix})
     end
+    for problem in get_boundary_problems(solver)
+        assemble!(problem, solver.time)
+    end
     t1 = round(toq(), 2)
     info("Assembled in $t1 seconds.")
     M, K, Kg, f = get_field_assembly(solver; with_mass_matrix=true)
+    Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
+    K = K + Kb
+    f = f + fb
     if solver.properties.geometric_stiffness
         K += Kg
     end
-    for problem in get_boundary_problems(solver)
-        assemble!(problem, solver.time)
-        # FIXME: Check for tie contacts and rhs. Here we just
-        # remove all fixed dofs giving funny results if problem
-        # is having MPCs or non-homogeneous Dirichlet conditions
-#       eliminate!(M, K, Kg, f, problem)
-        fixed_dofs = get_nonzero_rows(problem.assembly.C2)
-        K[fixed_dofs, :] = 0
-        K[:, fixed_dofs] = 0
-        M[fixed_dofs, :] = 0
-        M[:, fixed_dofs] = 0
-        f[fixed_dofs, :] = 0
-    end
-    fd = get_nonzero_rows(K)
+    
+    @assert nnz(D) == 0
+    @assert C1 == C2
+
+    tic()
+    P, h = create_projection(C1, g)
+    K_red = P'*K*P
+    M_red = P'*M*P
+    t1 = round(toq(), 2)
+    info("Eliminated dirichlet boundaries in $t1 seconds.")
+
+    nz = get_nonzero_rows(K_red)
     ndofs = solver.ndofs
     props = solver.properties
     info("Calculate $(props.nev) eigenvalues...")
-    if debug
+    if debug && length(nz) < 100
         info("Stiffness matrix:")
-        dump(round(full(K[fd, fd])))
+        dump(round(full(K[nz, nz])))
         info("Mass matrix:")
-        dump(round(full(M[fd, fd])))
+        dump(round(full(M[nz, nz])))
     end
+
     tic()
-    om2, X = eigs(K[fd, fd], M[fd, fd]; nev=props.nev, which=props.which)
+    om2, X = eigs(K_red[nz,nz], M_red[nz,nz]; nev=props.nev, which=props.which)
     props.eigvals = om2
     props.eigvecs = zeros(ndofs, length(om2))
-    props.eigvecs[fd, :] = X
+    v = zeros(ndofs)
+    for i=1:length(om2)
+        fill!(v, 0.0)
+        v[nz] = X[:,i]
+        props.eigvecs[:,i] = P*v + g
+    end
     t1 = round(toq(), 2)
     info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
     return true
diff --git a/src/mortar.jl b/src/mortar.jl
index ffb79f0..ddf97de 100644
--- a/src/mortar.jl
+++ b/src/mortar.jl
@@ -2,13 +2,14 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 type Mortar <: BoundaryProblem
+    dimension :: Int
     rotate_normals :: Bool
     adjust :: Bool
     tolerance :: Float64
 end
 
 function Mortar()
-    return Mortar(false, false, 0.0)
+    return Mortar(-1, false, false, 0.0)
 end
 
 function get_unknown_field_name(::Type{Mortar})
@@ -37,6 +38,10 @@ function cross2(a, b)
     cross([a; 0], [b; 0])[3]
 end
 
+function get_slave_elements(problem::Problem{Mortar})
+    filter(el -> haskey(el, "master elements"), get_elements(problem))
+end
+
 function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
     x1_ = slave_element["geometry"](time)
     n1_ = slave_element["normal"](time)
@@ -61,7 +66,7 @@ function project_from_slave_to_master{E<:MortarElements2D}(master_element::Eleme
     return xi2
 end
 
-function calculate_normals(elements, time, rotate_normals=false)
+function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false)
     tangents = Dict{Int64, Vector{Float64}}()
     for element in elements
         conn = get_connectivity(element)
@@ -79,7 +84,7 @@ function calculate_normals(elements, time, rotate_normals=false)
 
     Q = [0.0 -1.0; 1.0 0.0]
     normals = Dict{Int64, Vector{Float64}}()
-    S = sort(collect(keys(tangents)))
+    S = collect(keys(tangents))
     for j in S
         tangents[j] /= norm(tangents[j])
         normals[j] = Q*tangents[j]
@@ -94,8 +99,8 @@ function calculate_normals(elements, time, rotate_normals=false)
     return normals, tangents
 end
 
-function calculate_normals!(elements, time, rotate_normals=false)
-    normals, tangents = calculate_normals(elements, time, rotate_normals)
+function calculate_normals!(elements, time, ::Type{Val{1}}; rotate_normals=false)
+    normals, tangents = calculate_normals(elements, time, Val{1}; rotate_normals=rotate_normals)
     for element in elements
         conn = get_connectivity(element)
         update!(element, "normal", time => [normals[j] for j in conn])
@@ -104,32 +109,34 @@ function calculate_normals!(elements, time, rotate_normals=false)
 end
 
 function assemble!(problem::Problem{Mortar}, time::Real)
+    if problem.dimension == -1
+        error("set interface dimension: problem.properties.dimension = 1 or 2")
+    end
+    assemble!(problem, time, Val{problem.properties.dimension})
+end
+
+function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
 
     props = problem.properties
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
-    slave_elements = filter(el -> haskey(el, "master elements"), get_elements(problem))
+    slave_elements = get_slave_elements(problem)
 
     # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
-    normals, tangents = calculate_normals(slave_elements, time, props.rotate_normals)
+    normals, tangents = calculate_normals(slave_elements, time, Val{1};
+                                          rotate_normals=props.rotate_normals)
     update!(slave_elements, "normal", normals)
-    update!(slave_elements, "tangent", tangents)
-    S = sort(collect(keys(normals)))
 
     # 2. loop all slave elements
     for slave_element in slave_elements
         haskey(slave_element, "master elements") || continue
 
-        slave_element_nodes = get_connectivity(slave_element)
-        nsl = length(slave_element)
         X1 = slave_element["geometry"](time)
-        n1 = Field([normals[j] for j in slave_element_nodes])
+        n1 = slave_element["normal"](time)
 
         # 3. loop all master elements
         for master_element in slave_element["master elements"](time)
 
-            master_element_nodes = get_connectivity(master_element)
-            nm = length(master_element)
             X2 = master_element["geometry"](time)
 
             # 3.1 calculate segmentation
@@ -141,9 +148,11 @@ function assemble!(problem::Problem{Mortar}, time::Real)
 
             # 3.3. loop integration points of one integration segment and calculate
             # local mortar matrices
+            nsl = length(slave_element)
+            nm = length(master_element)
             De = zeros(nsl, nsl)
             Me = zeros(nsl, nm)
-            ge = zeros(nsl)
+            ge = zeros(field_dim*nsl)
             for ip in get_integration_points(slave_element, 2)
                 detJ = slave_element(ip, time, Val{:detJ})
                 w = ip.weight*detJ*l
@@ -159,10 +168,11 @@ function assemble!(problem::Problem{Mortar}, time::Real)
                 De += w*N1*N1'
                 Me += w*N1*N2'
                 if props.adjust
-                    g = X_s-X_m
-                    if g < props.tol
-                        ge += w*g
-                    end
+                    u1 = slave_element["displacement"](time)
+                    u2 = master_element["displacement"](time)
+                    x_s = X_s + N1*u1
+                    x_m = X_m + N2*u2
+                    ge += w*vec((x_m-x_s)*N1')
                 end
             end
 
@@ -177,12 +187,354 @@ function assemble!(problem::Problem{Mortar}, time::Real)
                 add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
                 add!(problem.assembly.C2, lsdofs, lsdofs, De)
                 add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
-                add!(problem.assembly.g, lsdofs, ge)
             end
-
+            add!(problem.assembly.g, sdofs, ge)
 
         end # master elements done
 
     end # slave elements done, contact virtual work ready
 
 end
+
+function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
+    return p - dot(p-x0, n0)*n0
+end
+
+function inv3(P::Matrix)
+    n, m = size(P)
+    @assert n == m == 3
+    a, b, c, d, e, f, g, h, i = P
+    A = e*i - f*h
+    B = -d*i + f*g
+    C = d*h - e*g
+    D = -b*i + c*h
+    E = a*i - c*g
+    F = -a*h + b*g
+    G = b*f - c*e
+    H = -a*f + c*d
+    I = a*e - b*d
+    return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I]
+end
+
+function vertex_inside_polygon(q, P; atol=1.0e-6)
+    N = length(P)
+    angle = 0.0
+    for i=1:N
+        A = P[i] - q
+        B = P[mod(i,N)+1] - q
+        c = norm(A)*norm(B)
+        isapprox(c, 0.0; atol=atol) && return true
+        cosa = dot(A,B)/c
+        isapprox(cosa, 1.0; atol=atol) && return false
+        isapprox(cosa, -1.0; atol=atol) && return true
+        try
+            angle += acos(cosa)
+        catch
+            info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
+            info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
+            info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
+            info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
+            rethrow()
+        end
+    end
+    return isapprox(angle, 2*pi; atol=atol)
+end
+
+function calculate_centroid(P)
+    N = length(P)
+    P0 = P[1]
+    areas = [norm(1/2*cross(P[i]-P0, P[mod(i,N)+1]-P0)) for i=2:N]
+    centroids = [1/3*(P0+P[i]+P[mod(i,N)+1]) for i=2:N]
+    C = 1/sum(areas)*sum(areas.*centroids)
+    return C
+end
+
+function get_cells(P, C)
+    N = length(P)
+    cells = Vector[]
+    # shared edge etc.
+    N < 3 && return cells
+    # trivial case, polygon already triangle / quadrangle
+    #N == 3 && return Vector[P]
+    #N == 4 && return Vector[P]
+    #V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N])
+    #A = 1/2*abs(dot(n, V))
+    #info("A = $A")
+    cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N]
+    return cells
+
+    maxa = 0.0
+    maxj = 0
+    for i=1:N
+        A = P[i] - C
+        B = P[mod(i,N)+1] - C
+        theta = acos(dot(A,B)/(norm(A)*norm(B)))
+        if theta > maxa
+            maxa = theta
+            maxj = i
+        end
+    end
+    info("max angle $(maxa/pi*180) at index $maxj, N=$N")
+    indices = mod(collect(maxj:maxj+N), N)
+    info("indices = $indices")
+end
+
+function get_polygon_clip(xs, xm, n; debug=false)
+    # objective: search does line xm1 - xm2 clip xs
+    nm = length(xm)
+    ns = length(xs)
+    P = Vector{Float64}[]
+
+    # 1. test is master point inside slave, if yes, add to clip
+    for i=1:nm
+        if vertex_inside_polygon(xm[i], xs)
+            debug && info("1. $(xm[i]) inside S -> push")
+            push!(P, xm[i])
+        end
+    end
+
+    # 2. test is slave point inside master, if yes, add to clip
+    for i=1:ns
+        if vertex_inside_polygon(xs[i], xm)
+            xs[i] in P && continue
+            debug && info("2. $(xs[i]) inside M -> push")
+            push!(P, xs[i])
+        end
+    end
+
+    for i=1:nm
+        # 2. find possible intersection
+        xm1 = xm[i]
+        xm2 = xm[mod(i,nm)+1]
+        #info("intersecting line $xm1 -> $xm2")
+        for j=1:ns
+            xs1 = xs[j]
+            xs2 = xs[mod(j,ns)+1]
+            #info("clipping polygon edge $xs1 -> $xs2")
+            tnom = dot(cross(xm1-xs1, xm2-xm1), n)
+            tdenom = dot(cross(xs2-xs1, xm2-xm1), n)
+            isapprox(tdenom, 0) && continue
+            t = tnom/tdenom
+            (0 <= t <= 1) || continue
+            q = xs1 + t*(xs2 - xs1)
+            #info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
+            if vertex_inside_polygon(q, xm)
+                q in P && continue
+                debug && info("3. $q inside M -> push")
+                push!(P, q)
+            end
+        end
+    end
+
+    return P
+end
+
+function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
+    element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
+    basis(xi) = get_basis(element, xi, time)
+    dbasis(xi) = get_dbasis(element, xi, time)
+    f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p
+    L(theta) = inv3([dbasis(theta[1:2])*x -n0])
+#   L2(theta) = inv(ForwardDiff.get_value([dbasis(theta[2:3])*x -n0]))
+    # FIXME: for some reason forwarddiff gives NaN's here.
+    theta = zeros(3)
+    dtheta = zeros(3)
+    for i=1:max_iterations
+        dtheta = L(theta) * f(theta)
+        theta -= dtheta
+        if norm(dtheta) < iter_tol
+            return theta[1:2], theta[3]
+        end
+    end
+
+    info("failed to project vertex from auxiliary plane back to surface")
+    info("element type: $E")
+    info("element connectivity: $(get_connectivity(element))")
+    info("auxiliary plane: x0 = $x0, n0 = $n0")
+    info("element geometry: $(x.data)")
+    info("vertex to project: $p")
+    info("parameter vector before giving up: $theta")
+    info("increment in parameter vector before giving up: $dtheta")
+    info("norm(dtheta) before giving up: $(norm(dtheta))")
+    info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
+    info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")
+
+    info("iterations:")
+    theta = zeros(3)
+    dtheta = zeros(3)
+    for i=1:max_iterations
+        info("iter $i, theta = $theta")
+        info("f = $(f(theta))")
+        info("L = $(L(theta))")
+        dtheta = L(theta) * f(theta)
+        info("dtheta = $(dtheta)")
+        theta -= dtheta
+    end
+
+    error("project_point_to_surface: did not converge in $max_iterations iterations!")
+end
+
+function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false)
+    normals = Dict{Int64, Vector{Float64}}()
+    for element in elements
+        conn = get_connectivity(element)
+        J = transpose(element([0.0, 0.0], time, Val{:Jacobian}))
+        normal = cross(J[:,1], J[:,2])
+        for nid in conn
+            if haskey(normals, nid)
+                normals[nid] += normal
+            else
+                normals[nid] = normal
+            end
+        end
+    end
+    # normalize to unit normal
+    S = collect(keys(normals))
+    for j in S
+        normals[j] /= norm(normals[j])
+    end
+    if rotate_normals
+        for j in S
+            normals[j] = -normals[j]
+        end
+    end
+    return normals
+end
+
+function check_orientation!(P, n; debug=false)
+    C = mean(P)
+    np = length(P)
+    s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np]
+    all(s .< 0) && return
+    debug && info("polygon not in ccw order, fixing")
+    # project points to new orthogonal basis Q and sort there
+    t1 = (P[1]-C)/norm(P[1]-C)
+    t2 = cross(n, t1)
+    Q = [n t1 t2]
+    sort!(P, lt=(A, B) -> begin
+        A_proj = Q'*(A-C)
+        B_proj = Q'*(B-C)
+        a = atan2(A_proj[3], A_proj[2])
+        b = atan2(B_proj[3], B_proj[2])
+        return a > b
+    end)
+end
+
+function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=true)
+
+    props = problem.properties
+    field_dim = get_unknown_field_dimension(problem)
+    field_name = get_parent_field_name(problem)
+    slave_elements = get_slave_elements(problem)
+    area = 0.0
+
+    # 1. calculate nodal normals and tangents for slave element nodes j ∈ S
+    normals = calculate_normals(slave_elements, time, Val{2};
+                                rotate_normals=props.rotate_normals)
+    update!(slave_elements, "normal", normals)
+
+    # 2. loop all slave elements
+    for slave_element in slave_elements
+
+        slave_element_nodes = get_connectivity(slave_element)
+        nsl = length(slave_element)
+        X1 = slave_element["geometry"](time)
+        n1 = Field([normals[j] for j in slave_element_nodes])
+
+        # project slave nodes to auxiliary plane (x0, Q)
+        #xi = get_reference_element_midpoint(slave_element)
+        xi = [1/3, 1/3]
+        N = vec(get_basis(slave_element, xi, time))
+        x0 = N*X1
+        n0 = N*n1
+        S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
+
+        # 3. loop all master elements
+        for master_element in slave_element["master elements"](time)
+
+            master_element_nodes = get_connectivity(master_element)
+            nm = length(master_element)
+            X2 = master_element["geometry"](time)
+
+            # 3.1 project master nodes to auxiliary plane and create polygon clipping
+            M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
+            P = get_polygon_clip(S, M, n0)
+            length(P) < 3 && continue # no clipping or shared edge (no volume)
+            check_orientation!(P, n0)
+            C0 = calculate_centroid(P)
+
+            De = zeros(nsl, nsl)
+            Me = zeros(nsl, nm)
+            ge = zeros(field_dim*nsl)
+
+            # 4. loop integration cells
+            for cell in get_cells(P, C0)
+                virtual_element = Element(Tri3)
+                update!(virtual_element, "geometry", cell)
+                #x_cell = Field(cell)
+
+                # 5. loop integration point of integration cell
+                for ip in get_integration_points(virtual_element, 3)
+                    N = vec(get_basis(virtual_element, ip, time))
+                    #dN = vec(get_dbasis(virtual_element, ip, time))
+                    #JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
+                    #wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
+                    detJ = virtual_element(ip, time, Val{:detJ})
+                    w = ip.weight*detJ
+
+                    # project gauss point from auxiliary plane to master and slave element
+                    #x_gauss = N*x_cell
+                    x_gauss = virtual_element("geometry", ip, time)
+                    if isnan(x_gauss[1])
+                        info("is nan")
+                        info("x_gauss = $x_gauss")
+                        info("cell = $cell")
+                        info("C0 = $C0")
+                        info("P = $P")
+                        info("S = $S")
+                        info("M = $M")
+                        info("n0 = $n0")
+                        error("nan, unable to continue")
+                    end
+                    xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
+                    xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
+
+                    # add contributions
+                    N1 = vec(get_basis(slave_element, xi_s, time))
+                    N2 = vec(get_basis(master_element, xi_m, time))
+                    De += w*N1*N1'
+                    Me += w*N1*N2'
+                    if props.adjust
+                        u1 = slave_element["displacement"](time)
+                        u2 = master_element["displacement"](time)
+                        x_s = N1*(X1+u1)
+                        x_m = N2*(X2+u2)
+                        ge += w*vec((x_m-x_s)*N1')
+                    end
+                    area += w
+                end # integration points done
+
+            end # integration cells done
+
+            # 6. add contribution to contact virtual work
+            sdofs = get_gdofs(problem, slave_element)
+            mdofs = get_gdofs(problem, master_element)
+            
+            for i=1:field_dim
+                lsdofs = sdofs[i:field_dim:end]
+                lmdofs = mdofs[i:field_dim:end]
+                add!(problem.assembly.C1, lsdofs, lsdofs, De)
+                add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
+                add!(problem.assembly.C2, lsdofs, lsdofs, De)
+                add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
+            end
+            add!(problem.assembly.g, sdofs, ge)
+
+        end # master elements done
+
+    end # slave elements done, contact virtual work ready
+
+    debug && info("area of interface: $area")
+
+end
+
diff --git a/src/postprocess_utils.jl b/src/postprocess_utils.jl
index 46a7ed4..26de9cd 100644
--- a/src/postprocess_utils.jl
+++ b/src/postprocess_utils.jl
@@ -31,3 +31,22 @@ function calc_nodal_values!(elements, field_name, field_dim, time)
     end
     update!(elements, field_name, nodal_values)
 end
+
+"""
+Return node ids + vector of values 
+"""
+function get_nodal_vector(elements, field_name, time)
+    f = Dict{Int64, Vector{Float64}}()
+    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
+    node_ids = sort(collect(keys(f)))
+    field = [f[nid] for nid in node_ids]
+    return node_ids, field
+end
+
diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl
index f93db48..31cce45 100644
--- a/src/preprocess_aster_reader.jl
+++ b/src/preprocess_aster_reader.jl
@@ -380,6 +380,7 @@ function aster_read_mesh(fn::ASCIIString, mesh_name=nothing)
         add_node!(mesh, nid, ncoords)
     end
     mapping = Dict(
+        :PO1 => :Poi1,
         :SE2 => :Seg2,
         :TR3 => :Tri3,
         :TR6 => :Tri6,
diff --git a/src/solvers.jl b/src/solvers.jl
index 2191a0a..6bac12a 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -216,80 +216,113 @@ function get_boundary_assembly(solver::Solver)
 end
 
 
-""" Solve linear system using LU factorization (UMFPACK).
 """
-function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver_UMFPACK}})
-    info("solving linear system of $(length(solver.problems)) problems.")
-    t0 = time()
+Construct new basis such that u = P*uh + g
 
-    # assemble field problems
-    M, K, Kg, f = get_field_assembly(solver)
+Parameters
+----------
+S  set of linearly independent dofs.
+"""
+function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
+    n, m = size(C)
+    @assert n == m
+    if S == nothing
+        S = get_nonzero_rows(C)
+    end
+    # FIXME: this creates dense matrices
+    # efficiency / memory usage is a question
+    P = sparse(C[S,:] \ full(C[S,:]))
+    h = sparse(C[S,:] \ full(g[S]))
+    resize!(P, n, m)
+    resize!(h, n, 1)
+    P = speye(n) - P
+    SparseMatrix.droptol!(P, tol)
+    return P, h
+end
 
-    # assemble boundary problems
-    Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
 
-    # construct global system Ax=b and solve using lu factorization
-    A = [
-        K+Kg+Kb  C1'
-        C2       D]
-    b = [f+fb; g]
+"""
+Solve linear system using LDLt factorization (SuiteSparse). This version
+requires that final system is symmetric and positive definite, so boundary
+conditions are first eliminated before solution.
+"""
+function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
 
+    nnz(D) == 0 || return 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
+
+    A = get_nonzero_rows(K)
+    I = setdiff(A, B)
+
+    if debug
+        info("# nz = $(length(nz))")
+        info("# A = $(length(A))")
+        info("# B = $(length(B))")
+        info("# I = $(length(I))")
+    end
+
+    # solver boundary dofs
+    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
+    end
+
+    # solve interior domain using LDLt factorization
+    u[I] = ldltfact(K[I,I]) \ (f[I] - K[I,B]*u[B])
+    # solve lambda
+    la[B] = lufact(C1[B,nz]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
+
+    return true
+end
+
+"""
+Solve linear system using LU factorization (UMFPACK). This version solves
+directly the saddle point problem without elimination of boundary conditions.
+"""
+function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
+    # construct global system Ax = b and solve using lufact (UMFPACK)
+    A = [K C1'; C2  D]
+    b = [f; g]
     nz = get_nonzero_rows(A)
     x = zeros(length(b))
     x[nz] = lufact(A[nz,nz]) \ full(b[nz])
-
-    ndofs = solver.ndofs
-    u = x[1:ndofs]
-    la = x[ndofs+1:end]
-    info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u))
-    return u, la
+    ndofs = size(K, 1)
+    u[:] = x[1:ndofs]
+    la[:] = x[ndofs+1:end]
+    return true
 end
 
-""" Solve linear system using LDLt factorization (SuiteSparse). """
-function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
+function solve_linear_system(solver::Solver)
     info("solving linear system of $(length(solver.problems)) problems.")
     t0 = time()
 
     # assemble field problems
     M, K, Kg, f = get_field_assembly(solver)
-
     # assemble boundary problems
     Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
-
-    K = K + Kb + Kg
-    f = f + fb
+    K = K + Kg + Kb
     K = 1/2*(K + K')
+    f = f + fb
+
     u = zeros(solver.ndofs)
     la = zeros(solver.ndofs)
 
-    # determine interior and boundary dofs
-    all_dofs = get_nonzero_rows(K)
-    boundary_dofs = get_nonzero_rows(C1)
-    boundary_dofs2 = get_nonzero_rows(C2)
-    interior_dofs = setdiff(all_dofs, boundary_dofs)
-    @assert length(boundary_dofs) == length(boundary_dofs2)
-    @assert setdiff(Set(boundary_dofs), Set(boundary_dofs2)) == Set()
-
-    # solve boundary
-    LUF = lufact(C1[boundary_dofs, boundary_dofs])
-    u[boundary_dofs] = LUF \ full(g[boundary_dofs])
-    normub = norm(u[boundary_dofs])
-    if isapprox(normub, 0.0)
-        info("CHOLMOD: homogeneous dirichlet boundary condition.")
+    status = false
+    for i in [1, 2]
+        status = solve!(K, C1, C2, D, f, g, u, la, Val{i})
+        if status
+            info("succesfully solved Ax = b using solver #$i")
+            break
+        end
     end
+    status || error("Failed to solve linear system!")
 
-    # solver interior
-    CF = ldltfact(K[interior_dofs, interior_dofs])
-    Kib = K[interior_dofs, boundary_dofs]
-    Kbb = K[boundary_dofs, boundary_dofs]
-    fi = f[interior_dofs]
-    u[interior_dofs] = CF \ (fi - Kib*u[boundary_dofs])
-
-    # solve lambda
-    # LUF2 = lufact(C2[boundary_dofs, boundary_dofs])
-    la[boundary_dofs] = LUF \ full(Kib' * u[interior_dofs] - Kbb*u[boundary_dofs])
-
-    info("CHOLMOD: solved in ", time()-t0, " seconds. norm = ", norm(u))
+    info("linear system solver: solved in ", time()-t0, " seconds. norm = ", norm(u))
     return u, la
 end
 
@@ -350,15 +383,19 @@ function assemble!(solver::Solver; force_assembly=true)
     info("Assembled in $t1 seconds.")
 end
 
+function initialize!(solver::Solver)
+    for problem in solver.problems
+        initialize!(problem, solver.time)
+    end
+end
+
 """ Default solver for quasistatic nonlinear problems. """
 function call(solver::Solver{Nonlinear})
 
     properties = solver.properties
 
     # 1. initialize each problem so that we can start nonlinear iterations
-    for problem in solver.problems
-        initialize!(problem, solver.time)
-    end
+    initialize!(solver)
 
     # 2. start non-linear iterations
     for properties.iteration=1:properties.max_iterations
@@ -370,7 +407,7 @@ function call(solver::Solver{Nonlinear})
         # 2.2 call solver for linearized system (default: direct lu factorization)
         info("Solve linear system ...")
         tic()
-        u, la = solve_linear_system(solver, Val{properties.linear_system_solver})
+        u, la = solve_linear_system(solver)
         push!(properties.norms, (norm(u), norm(la)))
         t1 = round(toq(), 2)
         info("Solved Ax = b in $t1 seconds.")
diff --git a/src/sparse.jl b/src/sparse.jl
index 799ad8c..2d09157 100644
--- a/src/sparse.jl
+++ b/src/sparse.jl
@@ -147,16 +147,19 @@ Returns
 Ordered list of row indices.
 """
 function get_nonzero_rows(A::SparseMatrixCSC)
-    # FIXME: This is probably a very inefficient way to do this.
     return sort(unique(rowvals(A)))
 end
 
-function get_nonzero_rows(A::SparseMatrixCOO)
+function get_nonzero_columns(A::SparseMatrixCSC)
+    return get_nonzero_rows(transpose(A))
+end
+
+function get_nonzero_rows(A::Union{SparseMatrixCOO, Matrix})
     return get_nonzero_rows(sparse(A))
 end
 
-function get_nonzero_rows(A::Matrix)
-    return get_nonzero_rows(sparse(A))
+function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix})
+    return get_nonzero_columns(sparse(A))
 end
 
 function size(A::SparseMatrixCOO)
diff --git a/test/test_elasticity_2d_linear_with_surface_load.jl b/test/test_elasticity_2d_linear_with_surface_load.jl
index 5030530..ddb0298 100644
--- a/test/test_elasticity_2d_linear_with_surface_load.jl
+++ b/test/test_elasticity_2d_linear_with_surface_load.jl
@@ -18,7 +18,7 @@ using JuliaFEM.Test
     block.elements = create_elements(mesh, "BLOCK")
     update!(block.elements, "youngs modulus", 288.0)
     update!(block.elements, "poissons ratio", 1/3)
-#   update!(block.elements, "displacement load 2", 576.0)
+    update!(block.elements, "displacement load 2", 576.0)
 
     traction = create_elements(mesh, "TOP")
     update!(traction, "displacement traction force 2", 288.0)
diff --git a/test/test_modal_analysis.jl b/test/test_modal_analysis.jl
index 2d7f4f9..3007045 100644
--- a/test/test_modal_analysis.jl
+++ b/test/test_modal_analysis.jl
@@ -42,3 +42,79 @@ using JuliaFEM.Test
     call(s1; debug=true)
     @test isapprox(s1.properties.eigvals, [5/3, 2/3])
 end
+
+@testset "test poisson problem modal analysis without tie" begin
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 3.0],
+        4 => [0.0, 3.0],
+        5 => [0.0, 3.0],
+        6 => [1.0, 3.0],
+        7 => [1.0, 9.0],
+        8 => [0.0, 9.0])
+    T = Dict{Int64, Float64}()
+    for i=1:8
+        T[i] = 0.0
+    end
+    el1 = Element(Quad4, [1, 2, 3, 4])
+    el2 = Element(Quad4, [4, 3, 7, 8])
+    el3 = Element(Seg2, [1, 2])
+    el4 = Element(Seg2, [7, 8])
+    update!([el1, el2, el3, el4], "geometry", X)
+    update!([el1, el2], "density", 6.0)
+    update!([el1, el2], "temperature thermal conductivity", 36.0)
+    update!([el1, el2], "temperature", T)
+    update!([el3, el4], "temperature 1", 0.0)
+    p1 = Problem(Heat, "combined body", 1)
+    p2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
+    push!(p1, el1, el2)
+    push!(p2, el3, el4)
+
+    solver = Solver(Modal)
+    push!(solver, p1, p2)
+    call(solver)
+    @test isapprox(solver.properties.eigvals[1], 1.0)
+end
+
+@testset "test poisson modal problem with mesh tie" begin
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 3.0],
+        4 => [0.0, 3.0],
+        5 => [0.0, 3.0],
+        6 => [1.0, 3.0],
+        7 => [1.0, 9.0],
+        8 => [0.0, 9.0])
+    T = Dict{Int64, Float64}()
+    for i=1:8
+        T[i] = 0.0
+    end
+    el1 = Element(Quad4, [1, 2, 3, 4])
+    el2 = Element(Quad4, [5, 6, 7, 8])
+    el3 = Element(Seg2, [1, 2])
+    el4 = Element(Seg2, [7, 8])
+    el5 = Element(Seg2, [3, 4])
+    el6 = Element(Seg2, [5, 6])
+    update!([el1, el2, el3, el4, el5, el6], "geometry", X)
+    update!([el1, el2], "temperature", T)
+    update!([el1, el2], "density", 6.0)
+    update!([el1, el2], "temperature thermal conductivity", 36.0)
+    update!([el3, el4], "temperature 1", 0.0)
+    update!(el5, "master elements", [el6])
+    p1 = Problem(Heat, "body 1", 1)
+    p2 = Problem(Heat, "body 2", 1)
+    p3 = Problem(Dirichlet, "fixed ends", 1, "temperature")
+    p4 = Problem(Mortar, "interface between bodies", 1, "temperature")
+    p4.properties.dimension = 1
+    push!(p1, el1)
+    push!(p2, el2)
+    push!(p3, el3, el4)
+    push!(p4, el5, el6)
+    solver = Solver(Modal)
+    push!(solver, p1, p2, p3, p4)
+    call(solver)
+    @test isapprox(solver.properties.eigvals[1], 1.0)
+end
+
diff --git a/test/test_mortar_2d_mesh_tie.jl b/test/test_mortar_2d_mesh_tie.jl
new file mode 100644
index 0000000..e40c786
--- /dev/null
+++ b/test/test_mortar_2d_mesh_tie.jl
@@ -0,0 +1,119 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Test
+
+function get_test_model()
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 0.5],
+        4 => [0.0, 0.5],
+        5 => [0.0, 0.6],
+        6 => [1.0, 0.6],
+        7 => [1.0, 1.1],
+        8 => [0.0, 1.1])
+    el1 = Element(Quad4, [1, 2, 3, 4])
+    el2 = Element(Quad4, [5, 6, 7, 8])
+    el3 = Element(Seg2, [1, 2])
+    el4 = Element(Seg2, [7, 8])
+    el5 = Element(Seg2, [4, 3])
+    el6 = Element(Seg2, [5, 6])
+    update!([el1, el2, el3, el4, el5, el6], "geometry", X)
+    update!([el1, el2], "youngs modulus", 96.0)
+    update!([el1, el2], "poissons ratio", 1/3)
+    update!([el3], "displacement 1", 0.0)
+    update!([el3], "displacement 2", 0.0)
+    update!([el4], "displacement 1", 0.0)
+    update!([el4], "displacement 2", 0.0)
+    update!(el6, "master elements", [el5])
+    p1 = Problem(Elasticity, "body1", 2)
+    p2 = Problem(Elasticity, "body2", 2)
+    p3 = Problem(Dirichlet, "fixed", 2, "displacement")
+    p4 = Problem(Mortar, "interface", 2, "displacement")
+    push!(p1, el1)
+    push!(p2, el2)
+    push!(p3, el3, el4)
+    push!(p4, el5, el6)
+    return p1, p2, p3, p4
+end
+
+@testset "test adjust setting in 2d tie contact" begin
+    p1, p2, p3, p4 = get_test_model()
+    p1.properties.formulation = :plane_stress
+    p2.properties.formulation = :plane_stress
+    p4.properties.dimension = 1
+    p4.properties.adjust = true
+    p4.properties.rotate_normals = false
+    solver = Solver(Nonlinear)
+    push!(solver, p1, p2, p3, p4)
+    call(solver)
+    el5 = p4.elements[1]
+    u = el5("displacement", [0.0], 0.0)
+    info("u = $u")
+    @test isapprox(u, [0.0, 0.05])
+end
+
+@testset "test that interface transfers constant field without error" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
+    mesh = aster_read_mesh(meshfile)
+
+    upper = Problem(Heat, "upper", 1)
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper.elements, "temperature thermal conductivity", 1.0)
+
+    lower = Problem(Heat, "lower", 1)
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower.elements, "temperature thermal conductivity", 1.0)
+
+    bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
+    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
+    update!(bc_upper.elements, "temperature 1", 0.0)
+
+    bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
+    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
+    update!(bc_lower.elements, "temperature 1", 1.0)
+
+    interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
+    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
+    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
+    update!(interface_slave_elements, "master elements", interface_master_elements)
+    interface.elements = [interface_master_elements; interface_slave_elements]
+    interface.properties.dimension = 1
+
+    solver = Solver()
+    push!(solver, upper, lower, bc_upper, bc_lower, interface)
+    call(solver)
+
+    node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
+    T = [t[1] for t in temperature]
+    minT = minimum(T)
+    maxT = maximum(T)
+    info("minT = $minT, maxT = $maxT")
+    @test isapprox(minT, 0.5)
+    @test isapprox(maxT, 0.5)
+end
+
+
+#=
+
+@testset "expect clear error when trying to solve 2d model in 3d setting" begin
+    p1, p2, p3, p4 = get_test_model()
+#   p1.properties.formulation = :plane_stress
+#   p2.properties.formulation = :plane_stress
+    p4.properties.adjust = true
+    p4.properties.rotate_normals = false
+    solver = Solver(Nonlinear)
+    solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK
+    push!(solver, p1, p2, p3, p4)
+    call(solver)
+    el5 = p4.elements[1]
+    u = el5("displacement", [0.0], 0.0)
+    info("u = $u")
+    @test isapprox(u, [0.0, 0.05])
+end
+
+=#
diff --git a/test/test_mortar_3d_mesh_tie.jl b/test/test_mortar_3d_mesh_tie.jl
new file mode 100644
index 0000000..3d2305a
--- /dev/null
+++ b/test/test_mortar_3d_mesh_tie.jl
@@ -0,0 +1,47 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Preprocess
+using JuliaFEM.Postprocess
+using JuliaFEM.Test
+
+@testset "test that interface transfers constant field without error" begin
+    meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med"
+    mesh = aster_read_mesh(meshfile)
+
+    upper = Problem(Heat, "upper", 1)
+    upper.elements = create_elements(mesh, "UPPER")
+    update!(upper.elements, "temperature thermal conductivity", 1.0)
+    lower = Problem(Heat, "lower", 1)
+    lower.elements = create_elements(mesh, "LOWER")
+    update!(lower.elements, "temperature thermal conductivity", 1.0)
+
+    bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
+    bc_upper.elements = create_elements(mesh, "UPPER_TOP")
+    update!(bc_upper.elements, "temperature 1", 0.0)
+
+    bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
+    bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
+    update!(bc_lower.elements, "temperature 1", 1.0)
+
+    interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
+    interface_slave_elements = create_elements(mesh, "LOWER_TOP")
+    interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
+    update!(interface_slave_elements, "master elements", interface_master_elements)
+    interface.elements = [interface_master_elements; interface_slave_elements]
+    interface.properties.dimension = 2
+
+    solver = Solver()
+    solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK
+    push!(solver, upper, lower, bc_upper, bc_lower, interface)
+    call(solver)
+    
+    node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
+    T = [t[1] for t in temperature]
+    minT = minimum(T)
+    maxT = maximum(T)
+    info("minT = $minT, maxT = $maxT")
+    @test isapprox(minT, 0.5)
+    @test isapprox(maxT, 0.5)
+end
diff --git a/test/test_mortar_3d_polygon_clip.jl b/test/test_mortar_3d_polygon_clip.jl
new file mode 100644
index 0000000..c5b9a7b
--- /dev/null
+++ b/test/test_mortar_3d_polygon_clip.jl
@@ -0,0 +1,41 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "polygon clip case 1" begin
+    S = Vector[
+        [0.375, 0.0, 0.5],
+        [0.6, 0.0, 0.5],
+        [0.5, 0.25, 0.5]]
+    M = Vector[
+        [0.50, 0.0, 0.5],
+        [0.25, 0.0, 0.5],
+        [0.375, 0.25, 0.5]]
+    n0 = [0.0, 0.0, 1.0]
+    P = get_polygon_clip(S, M, n0)
+    P_expected = Vector{Float64}[
+        [0.500, 0.0, 0.5],
+        [0.375, 0.0, 0.5],
+        [0.4375, 0.125, 0.5]]
+    @test length(P) == length(P_expected)
+    for (Pi, Pj) in zip(P, P_expected)
+        @test isapprox(Pi, Pj)
+    end
+end
+
+@testset "polygon clip case 2" begin
+    S = Vector[
+        [0.25, 0.0, 0.5],
+        [0.75, 0.0, 0.5],
+        [0.50, 0.25, 0.5]]
+    M = Vector[
+        [0.50, 0.0, 0.5],
+        [0.25, 0.0, 0.5],
+        [0.375, 0.25, 0.5]]
+    n0 = [0.0, 0.0, 1.0]
+    P = get_polygon_clip(S, M, n0)
+    @test length(P) == 3
+end
+
diff --git a/test/test_projection.jl b/test/test_projection.jl
new file mode 100644
index 0000000..ea97365
--- /dev/null
+++ b/test/test_projection.jl
@@ -0,0 +1,32 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "test projection" begin
+    C = [
+        2.0  1.0  0.0  0.0   0.0   0.0  0.0  0.0
+        1.0  2.0  0.0  0.0   0.0   0.0  0.0  0.0
+        0.0  0.0  2.0  1.0  -1.0  -2.0  0.0  0.0
+        0.0  0.0  1.0  2.0  -2.0  -1.0  0.0  0.0
+        0.0  0.0  0.0  0.0   0.0   0.0  0.0  0.0
+        0.0  0.0  0.0  0.0   0.0   0.0  0.0  0.0
+        0.0  0.0  0.0  0.0   0.0   0.0  0.0  0.0
+        0.0  0.0  0.0  0.0   0.0   0.0  0.0  0.0]
+    g = [3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    P, h = create_projection(sparse(C), g)
+    P_expected = [
+        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+        0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
+        0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
+        0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
+        0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0
+        0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0
+        0.0  0.0  0.0  0.0  0.0  0.0  1.0  0.0
+        0.0  0.0  0.0  0.0  0.0  0.0  0.0  1.0]
+    h_expected = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    @test isapprox(full(P), P_expected)
+    @test isapprox(full(h), h_expected)
+end
+