Cleanup code (#7)

Remove tests not used. Calculate mortar matrix P and test that it's
able to transfer constant and linear load between two non-conforming
surfaces. Create function to calculate assembly.

src/Mortar2D.jl

@@ -6,4 +6,5 @@ include("calculate_normals.jl")
 include("calculate_projections.jl")
 include("calculate_segments.jl")
 include("calculate_mortar_matrices.jl")
+include("calculate_assembly.jl")
 end

src/calculate_assembly.jl

@@ -0,0 +1,58 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE
+
+""" Given data, calculate projection matrix `P`. """
+function calculate_mortar_assembly(elements, element_types, coords,
+                                   slave_element_ids, master_element_ids)
+    S = Int[]
+    M = Int[]
+    for sid in slave_element_ids
+        push!(S, elements[sid]...)
+    end
+    for mid in master_element_ids
+        push!(M, elements[mid]...)
+    end
+    S = sort(unique(S))
+    M = sort(unique(M))
+    ns = length(S)
+    nm = length(M)
+    println("$ns slave nodes, $nm master nodes")
+    normals = calculate_normals(elements, element_types, coords)
+    segmentation = calculate_segments(slave_element_ids, master_element_ids,
+                                      elements, element_types, coords, normals)
+    D_I = Int[]
+    D_J = Int[]
+    D_V = Float64[]
+    M_I = Int[]
+    M_J = Int[]
+    M_V = Float64[]
+    for sid in slave_element_ids
+        mids = [mid for (mid, xi) in segmentation[sid]]
+        De, Me = calculate_mortar_matrices(sid, mids, elements, element_types,
+                                           coords, normals, segmentation)
+        for (i,I) in enumerate(elements[sid])
+            for (j,J) in enumerate(elements[sid])
+                push!(D_I, I)
+                push!(D_J, J)
+                push!(D_V, De[i,j])
+            end
+        end
+        for mid in mids
+            for (i,I) in enumerate(elements[sid])
+                for (j,J) in enumerate(elements[mid])
+                   push!(M_I, I)
+                   push!(M_J, J)
+                   push!(M_V, Me[mid][i,j])
+                end
+            end
+        end
+    end
+
+    D = sparse(D_I, D_J, D_V)
+    Df = cholfact(1/2*(D+D')[S,S])
+
+    M_ = sparse(M_I, M_J, M_V)
+    P =  Df \ M_[S,M]
+    SparseArrays.droptol!(P, 1.0e-12)
+    return P
+end

test/test_calculate_assembly.jl

@@ -2,56 +2,35 @@
 # License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE
 
 using Base.Test
-using Mortar2D: calculate_mortar_matrices, calculate_segments
+using Mortar2D: calculate_mortar_assembly
 
-@testset "calculate bigger interface" begin
-    ns = 4 # number of slaves in interface
-    nm = 5 # number of masters in interface
+@testset "calculate mortar assembly" begin
+    ns = 10 # number of slaves in interface
+    nm = 9 # number of masters in interface
 
     Xs = Dict(i => [x, 1/2] for (i, x) in enumerate(linspace(0, 1, ns)))
     Es = Dict(i => [i, i+1] for i=1:(ns-1))
     Xm = Dict(ns+i => [x, 1/2] for (i, x) in enumerate(linspace(0, 1, nm)))
     Em = Dict(i => [i, i+1]+1 for i=ns:ns+nm-2)
+    element_types = Dict(i => :Seg2 for i=1:(ns+nm))
 
-    println(Xs)
-    println(Xm)
-    println(Es)
-    println(Em)
-
-    slave_ids = keys(Es)
-    master_ids = keys(Em)
+    slave_element_ids = keys(Es)
+    master_element_ids = keys(Em)
     elements = merge(Es, Em)
     coords = merge(Xs, Xm)
-    element_types = Dict(i => :Seg2 for i=1:(ns+nm))
-    normals = Dict(i => [0.0, 1.0] for (i, x) in enumerate(linspace(0, 1, ns)))
-    segmentation = calculate_segments(slave_ids, master_ids,
-                               elements, element_types,
-                               coords, normals)
-    B = zeros(ns, ns+nm)
-    for sid in keys(Es)
-        mids = []
-        for (mid, xi) in segmentation[sid]
-            push!(mids, mid)
-        end
-        De, Me = calculate_mortar_matrices(sid, mids, elements, element_types,
-                                           coords, normals, segmentation)
-        B[elements[sid], elements[sid]] += De
-        for mid in mids
-            B[elements[sid], elements[mid]] += Me[mid]
-        end
-    end
 
-    # next, calculate projection P
-    D = B[1:ns,1:ns]
-    M = B[1:ns,ns+1:ns+nm]
-    println(D)
-    println(M)
-    P = inv(D) * M
-    P[abs.(P) .< 1.0e-9] = 0
-    println(P)
+    P = calculate_mortar_assembly(elements, element_types, coords,
+                                  slave_element_ids, master_element_ids)
+
+    # test case 1, constant pressure load
+    um = ones(nm)
+    us = P*um
+    us_expected = ones(ns)
+    @test isapprox(us, us_expected)
+
+    # test case 2, linear load
     um = linspace(0, 1, nm)
     us = P*um
-    println(us)
-    @test isapprox(us, linspace(0, 1, ns))
+    us_expected = linspace(0, 1, ns)
+    @test isapprox(us, us_expected)
 end
-