Calculate bigger assembly (#6)

* calculate bigger assembly with several elements
* be more consistent with types
* throw exception if discriminant is negative
* cheap check to determine overlapping of master and slave element
* added failing test (discriminant of negative for unknown reason)

src/calculate_normals.jl

@@ -12,8 +12,10 @@ in nodes. As a result we get unique normal direction defined to each node.
 
 - Yang2005
 """
-function calculate_normals(elements::Dict{Int, Vector{Int}}, element_types::Dict{Int, Symbol}, X::Dict{Int, Vector{Float64}})
-    normals = Dict{Int64, Vector{Float64}}()
+function calculate_normals{T<:Integer,P<:AbstractFloat}(elements::Dict{T, Vector{T}},
+                                                        element_types::Dict{T, Symbol},
+                                                        X::Dict{T, Vector{P}})
+    normals = Dict{T, Vector{P}}()
     for (elid, elcon) in elements
         @assert element_types[elid] == :Seg2
         d = X[elcon[2]] - X[elcon[1]]

src/calculate_projections.jl

@@ -41,6 +41,11 @@ function project_from_master_to_slave{T<:Number}(::Type{Val{:Seg2}}, xm::Vector{
     b = nom_b / denom_b
     c = nom_c / denom_c
     d = b^2 - 4*a*c
+    if d < 0
+        warn("Mortar2D.calculate_projections(): negative discriminant")
+        warn("xm = $xm, xs1 = $xs1, xs2 = $xs2, ns1 = $ns1, ns2 = $ns2")
+        error("negative discriminant d=$d when calculating projection")
+    end
     sols = [-b + sqrt(d), -b - sqrt(d)]/(2.0*a)
     return sols[indmin(abs.(sols))]
 end

src/calculate_segments.jl

@@ -25,11 +25,16 @@ function calculate_segments(slave_ids, master_ids, elements, element_types, coor
             mcon = elements[mid]
             xm1 = coords[mcon[1]]
             xm2 = coords[mcon[2]]
+            # first project from slave to master, to find out
+            # are we completely outside of domain
+            xi2a = project_from_slave_to_master(Val{:Seg2}, xs1, ns1, xm1, xm2)
+            xi2b = project_from_slave_to_master(Val{:Seg2}, xs2, ns2, xm1, xm2)
+            xi2a >  1.0 && xi2b >  1.0 && continue
+            xi2a < -1.0 && xi2b < -1.0 && continue
             xi1a = project_from_master_to_slave(Val{:Seg2}, xm1, xs1, xs2, ns1, ns2)
             xi1b = project_from_master_to_slave(Val{:Seg2}, xm2, xs1, xs2, ns1, ns2)
             xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
-            l = 1/2*abs(xi1[2]-xi1[1])
-            isapprox(l, 0.0) && continue
+            isapprox(abs(xi1[2]-xi1[1]), 0.0) && continue
             push!(S[sid], (mid, xi1))
         end
     end

test/runtests.jl

@@ -10,6 +10,7 @@ push!(test_files, "test_calculate_normals.jl")
 push!(test_files, "test_calculate_projections.jl")
 push!(test_files, "test_calculate_segments.jl")
 push!(test_files, "test_calculate_mortar_matrices.jl")
+push!(test_files, "test_calculate_assembly.jl")
 
 @testset "Mortar2D.jl" begin
     for fn in test_files

test/test_calculate_assembly.jl

@@ -0,0 +1,57 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE
+
+using Base.Test
+using Mortar2D: calculate_mortar_matrices, calculate_segments
+
+@testset "calculate bigger interface" begin
+    ns = 4 # number of slaves in interface
+    nm = 5 # 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)
+
+    println(Xs)
+    println(Xm)
+    println(Es)
+    println(Em)
+
+    slave_ids = keys(Es)
+    master_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)
+    um = linspace(0, 1, nm)
+    us = P*um
+    println(us)
+    @test isapprox(us, linspace(0, 1, ns))
+end
+

test/test_calculate_projections.jl

@@ -22,3 +22,13 @@ end
     xi1 = project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2)
     @test isapprox(xi1, -0.281575016087237)
 end
+
+@testset "project slave node to master surface" begin
+    xm = [111.013, 283.152]
+    xs1 = [244.82, 366.732]
+    xs2 = [291.474, 359.169]
+    ns1 = [-0.0765221, 0.997068]
+    ns2 = [0.251347, 0.967897]
+    # FIXME: fails for unknown reason
+    # project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2)
+end