Merge 66a1dda into e2b4f54

.travis.yml

@@ -3,6 +3,10 @@ os:
     - linux
 julia:
     - 0.6
+addons:
+    apt:
+        packages:
+            - hdf5-tools
 before_script:
     - julia --color=yes -e 'Pkg.clone("https://github.com/JuliaFEM/PkgTestSuite.jl.git")'
     - julia --color=yes -e 'using PkgTestSuite; init()'

src/contact2dad.jl

@@ -4,13 +4,20 @@
 type Contact2DAD <: BoundaryProblem
     master_elements :: Vector{Element}
     rotate_normals :: Bool
+    dual_basis :: Bool
     iteration :: Int
     contact_state_in_first_iteration :: Symbol
-    always_in_contact :: Bool
+    always_in_contact :: Set{Int64}
+    use_scaling :: Bool
+    alpha :: Dict{Int64, Real}
+    beta :: Dict{Int64, Real}
+    nadj_nodes :: Dict{Int64, Real}
+    scaling_factors :: Dict{Int64, Real}
 end
 
 function Contact2DAD()
-    return Contact2DAD([], false, 0, :AUTO, false)
+    return Contact2DAD([], false, true, 0, :AUTO, Set(),
+                        true, Dict(), Dict(), Dict(), Dict())
 end
 
 function FEMBase.add_elements!(::Problem{Contact2DAD}, ::Any)
@@ -143,6 +150,7 @@ function project_from_slave_to_master_ad{E<:MortarElements2D}(
 
 end
 
+
 """
 Frictionless 2d finite sliding contact with forwarddiff.
 
@@ -209,6 +217,10 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass
         end
         update!(slave_elements, "normal", time => normals2)
 
+        alpha = empty!(problem.properties.alpha)
+        beta = empty!(problem.properties.beta)
+        nadj_nodes = empty!(problem.properties.nadj_nodes)
+
         # 2. loop all slave elements
         for slave_element in slave_elements
 
@@ -220,45 +232,69 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass
             n1 = ((normals[:,i] for i in slave_element_nodes)...)
             nnodes = size(slave_element, 2)
 
-            # construct dual basis
-            De = zeros(nnodes, nnodes)
-            Me = zeros(nnodes, nnodes)
+            Ae = eye(nnodes)
 
-            nsegments = 0
+            if problem.properties.dual_basis # construct dual basis
 
-            for master_element in get_master_elements(problem)
+                De = zeros(nnodes, nnodes)
+                Me = zeros(nnodes, nnodes)
+                ae = zeros(nnodes)
+                be = zeros(nnodes)
 
-                master_element_nodes = get_connectivity(master_element)
-                X2 = master_element("geometry", time)
-                u2 = ((u[:,i] for i in master_element_nodes)...)
-                x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
-
-                # calculate segmentation: we care only about endpoints
-                xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
-                xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
-                xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
-                l = 1/2*abs(xi1[2]-xi1[1])
-                isapprox(l, 0.0) && continue # no contribution in this master element
-
-                nsegments += 1
                 for ip in get_integration_points(slave_element, 3)
-                    # jacobian of slave element in deformed state
                     dN = get_dbasis(slave_element, ip, time)
                     j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
-                    w = ip.weight*norm(j)*l
-                    xi = ip.coords[1]
-                    xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
-                    N1 = get_basis(slave_element, xi_s, time)
-                    De += w*diagm(vec(N1))
-                    Me += w*N1'*N1
+                    detj = norm(j)
+                    w = ip.weight*detj
+                    N1 = slave_element(ip, time)
+                    ae += w*vec(N1)/detj
                 end
-            end
 
-            if nsegments == 0
-                continue
-            end
+                nsegments = 0
+
+                for master_element in get_master_elements(problem)
+
+                    master_element_nodes = get_connectivity(master_element)
+                    X2 = master_element("geometry", time)
+                    u2 = ((u[:,i] for i in master_element_nodes)...)
+                    x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
+
+                    # calculate segmentation: we care only about endpoints
+                    xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
+                    xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
+                    xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
+                    l = 1/2*abs(xi1[2]-xi1[1])
+                    isapprox(l, 0.0) && continue # no contribution in this master element
+
+                    nsegments += 1
+                    for ip in get_integration_points(slave_element, 3)
+                        # jacobian of slave element in deformed state
+                        dN = get_dbasis(slave_element, ip, time)
+                        j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
+                        detj = norm(j)
+                        w = ip.weight*norm(j)*l
+                        xi = ip.coords[1]
+                        xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
+                        N1 = get_basis(slave_element, xi_s, time)
+                        De += w*diagm(vec(N1))
+                        Me += w*N1'*N1
+                        be += w*vec(N1)/detj
+                    end
+                end
+
+                for (i, j) in enumerate(get_connectivity(slave_element))
+                    alpha[j] = get(alpha, j, 0.0) + ae[i]
+                    beta[j] = get(beta, j, 0.0) + be[i]
+                    nadj_nodes[j] = get(nadj_nodes, j, 0) + 1
+                end
 
-            Ae = De*inv(Me)
+                if nsegments == 0
+                    continue
+                end
+
+                Ae = De*inv(Me)
+
+            end
 
             # 3. loop all master elements
             for master_element in get_master_elements(problem)
@@ -374,14 +410,30 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass
                 t = Q'*n
                 lan = dot(n, la[:,j])
                 lat = dot(t, la[:,j])
-                C[1,j] += gap[1, j]
+                C[1,j] -= gap[1, j]
                 C[2,j] += lat
             else
                 C[:,j] = la[:,j]
             end
 
         end
 
+        # Apply scaling to contact virtual work and contact constraints
+        if problem.properties.use_scaling
+            scaling = empty!(problem.properties.scaling_factors)
+            for j in S
+                isapprox(beta[j], 0.0) && continue
+                #is_active[j] || continue
+                #haskey(alpha, j) || continue
+                #haskey(beta, j) || continue
+                #haskey(nadj_nodes, j) || continue
+                scaling[j] = alpha[j] / (nadj_nodes[j] * beta[j])
+                C[1,j] *= scaling[j]
+                fc[:,j] *= scaling[j]
+            end
+            update!(slave_elements, "contact scaling", time => scaling)
+        end
+
         return vec([fc C])
 
     end
@@ -395,7 +447,7 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass
     A = ForwardDiff.jacobian(calculate_interface, x)
     b = calculate_interface(x)
     A = sparse(A)
-    b = sparse(b)
+    b = -sparse(b)
     SparseArrays.droptol!(A, 1.0e-9)
     SparseArrays.droptol!(b, 1.0e-9)
 
@@ -404,8 +456,8 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2DAD}, assembly::Ass
     C1 = A[1:ndofs,ndofs+1:end]
     C2 = A[ndofs+1:end,1:ndofs]
     D = A[ndofs+1:end,ndofs+1:end]
-    f = -b[1:ndofs]
-    g = -b[ndofs+1:end]
+    f = b[1:ndofs]
+    g = b[ndofs+1:end]
 
     f += C1*problem.assembly.la
     g += D*problem.assembly.la

src/mortar2dad.jl

@@ -3,11 +3,10 @@
 
 type Mortar2DAD <: BoundaryProblem
     master_elements :: Vector{Element}
-    rotate_normals :: Bool
 end
 
 function Mortar2DAD()
-    return Mortar2DAD([], false)
+    return Mortar2DAD([])
 end
 
 function FEMBase.add_elements!(::Problem{Mortar2DAD}, ::Any)
@@ -51,9 +50,8 @@ function get_master_dofs(problem::Problem{Mortar2DAD})
     return sort(unique(dofs))
 end
 
-function project_from_master_to_slave_ad{E<:MortarElements2D}(
-    slave_element::Element{E}, x1_, n1_, x2, time;
-    tol=1.0e-10, max_iterations=20)
+function project_from_master_to_slave_ad(slave_element::Element{E}, x1_, n1_, x2, time;
+                    tol=1.0e-9, max_iterations=5) where {E<:MortarElements2D}
 
     x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_)
     dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_)
@@ -73,12 +71,12 @@ function project_from_master_to_slave_ad{E<:MortarElements2D}(
         end
     end
 
-    info("x1 = $(ForwardDiff.get_value(x1_.data))")
-    info("n1 = $(ForwardDiff.get_value(n1_.data))")
-    info("x2 = $(ForwardDiff.get_value(x2))")
-    info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
-    info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
-    info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
+    info("x1 = $x1")
+    info("n1 = $n1")
+    info("x2 = $x2")
+    info("xi1 = $xi1, dxi1 = $dxi1")
+    info("-R(xi1) = $(-R(xi1))")
+    info("dR(xi1) = $(dR(xi1))")
     error("find projection from master to slave: did not converge")
 
 end
@@ -119,9 +117,6 @@ function FEMBase.assemble_elements!(problem::Problem{Mortar2DAD}, assembly::Asse
     field_dim = get_unknown_field_dimension(problem)
     field_name = get_parent_field_name(problem)
     slave_elements = get_slave_elements(problem)
-    if field_name != "displacement"
-        error("mortar forwarddiff assembly: only displacement field with adjust=yes supported")
-    end
 
     function calculate_interface(x::Vector)
 
@@ -156,12 +151,6 @@ function FEMBase.assemble_elements!(problem::Problem{Mortar2DAD}, assembly::Asse
             normals[:,j] = Q*tangents[:,j]
         end
 
-        if props.rotate_normals
-            for j in S
-                normals[:,j] = -normals[:,j]
-            end
-        end
-
         #update!(slave_elements, "normal", time => Dict(j => normals[:,j] for j in S))
         #update!(slave_elements, "tangent", time => Dict(j => tangents[:,j] for j in S))
 

test/REQUIRE

@@ -0,0 +1 @@
+JLD

test/runtests.jl

@@ -6,4 +6,7 @@ using Base.Test
 
 @testset "MortarContact2DAD.jl tests" begin
     @testset "test_01.jl" begin include("test_01.jl") end
+    @testset "test_02.jl" begin include("test_02.jl") end
+    @testset "test_contact_1.jl" begin include("test_contact_1.jl") end
+    @testset "test_contact_2.jl" begin include("test_contact_2.jl") end
 end

test/test_01.jl

@@ -2,7 +2,9 @@
 # License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE
 
 using MortarContact2DAD
-using Base.Test
+using MortarContact2DAD: get_slave_dofs, get_master_dofs
+using MortarContact2DAD: project_from_master_to_slave_ad, project_from_slave_to_master_ad
+using FEMBase.Test
 
 X = Dict(
     1 => [0.0, 0.0],
@@ -20,6 +22,7 @@ slave = Element(Seg2, [1, 2])
 master = Element(Seg2, [3, 4])
 update!([slave, master], "geometry", X)
 update!([slave, master], "displacement", u)
+update!(slave, "normal", ([0.0, 1.0], [0.0, 1.0]))
 problem = Problem(Mortar2DAD, "test problem", 2, "displacement")
 add_slave_elements!(problem, [slave])
 add_master_elements!(problem, [master])
@@ -35,3 +38,16 @@ C1 = full(problem.assembly.C1)
 C2 = full(problem.assembly.C2)
 @test isapprox(C1, C2)
 @test isapprox(C1, C1_expected)
+
+@test_throws Exception add_elements!(problem, [slave])
+@test get_slave_dofs(problem) == [1, 2, 3, 4]
+@test get_master_dofs(problem) == [5, 6, 7, 8]
+x1 = slave("geometry", 0.0)
+n1 = slave("normal", 0.0)
+xm = 1/3*(X[3]+X[4])
+xs = 1/3*(X[1]+X[2])
+@test_throws ErrorException project_from_master_to_slave_ad(slave, x1, n1, xm, 0.0; tol=0.0, max_iterations=1)
+x1 = [0.5, 0.0]
+n1 = [0.0, 1.0]
+x2 = master("geometry", 0.0)
+@test_throws ErrorException project_from_slave_to_master_ad(master, x1, n1, x2, 0.0; tol=0.0, max_iterations=1)

test/test_02.jl

@@ -0,0 +1,61 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/MortarContact2DAD.jl/blob/master/LICENSE
+
+using MortarContact2DAD
+using FEMBase.Test
+import FEMBase.Test: @test_resource
+
+X = Dict(1 => [-2.0, 0.0], 2 => [ 0.0, 0.0], 3 => [ 2.0, 0.0],
+         4 => [-2.0, 0.0], 5 => [ 0.0, 0.0], 6 => [-2.0, -2.0],
+         7 => [0.0, -2.0], 8 => [2.0, -2.0], 9 => [-2.0, 2.0],
+         10 => [0.0, 2.0])
+
+ D = [0.5, 0.0]
+
+for j in [4, 5, 9, 10]
+    X[j] += D
+end
+
+contact = Problem(Contact2DAD, "contact", 2, "displacement")
+element7 = Element(Seg2, [1, 2])
+element8 = Element(Seg2, [2, 3])
+element9 = Element(Seg2, [4, 5])
+contact_slave_elements = [element7, element8]
+contact_master_elements = [element9]
+add_slave_elements!(contact, contact_slave_elements)
+add_master_elements!(contact, contact_master_elements)
+all_elements = [element7, element8, element9]
+update!(all_elements, "geometry", X)
+
+using JLD
+i1 = load(@test_resource("iter_1.jld"))
+i2 = load(@test_resource("iter_2.jld"))
+i3 = load(@test_resource("iter_3.jld"))
+
+function to_dict(u, ndofs, nnodes)
+    return Dict(j => [u[ndofs*(j-1)+i] for i=1:ndofs] for j=1:nnodes)
+end
+
+ndofs = 20
+contact.properties.iteration = 1
+for (data, time) in zip([i1, i2, i3], [0.0, 1.0, 2.0])
+    println("Testing data for iteration / time $time ")
+    empty!(contact.assembly)
+    contact.assembly.u = data["u"]
+    contact.assembly.la = data["la"]
+    update!(contact, "displacement", time => to_dict(data["u"], 2, 5))
+    update!(contact, "lambda", time => to_dict(data["la"], 2, 5))
+    assemble!(contact, time)
+    K = full(contact.assembly.K, ndofs, ndofs)
+    C1 = full(contact.assembly.C1, ndofs, ndofs)
+    C2 = full(contact.assembly.C2, ndofs, ndofs)
+    D = full(contact.assembly.D, ndofs, ndofs)
+    f = full(contact.assembly.f, ndofs, 1)
+    g = full(contact.assembly.g, ndofs, 1)
+    @test isapprox(K, data["K"])
+    @test isapprox(C1, data["C1"])
+    @test isapprox(C2, data["C2"])
+    @test isapprox(D, data["D"])
+    @test isapprox(f, data["f"])
+    @test isapprox(g, data["g"])
+end