From 3241dbb98bc506e2fdd81a9a9d63dfe0f9282df7 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Thu, 17 May 2018 19:16:29 +0300 Subject: [PATCH] Fix coverage (#2) Code coverage is back to 100 % now. --- src/contact2d.jl | 47 ++++--- src/mortar2d.jl | 143 +++++++++----------- test/runtests.jl | 24 +++- test/test_contact2d.jl | 134 +++++++++++++++++++ test/test_contact_segmentation.jl | 82 ++++++++++++ test/test_mortar2d_calculate_projection.jl | 145 ++++++++++----------- test/test_mortar2d_ex1.jl | 12 +- test/test_mortar2d_ex2.jl | 55 ++++++++ 8 files changed, 454 insertions(+), 188 deletions(-) create mode 100644 test/test_contact2d.jl create mode 100644 test/test_contact_segmentation.jl create mode 100644 test/test_mortar2d_ex2.jl diff --git a/src/contact2d.jl b/src/contact2d.jl index 8d7b650..8889dcb 100644 --- a/src/contact2d.jl +++ b/src/contact2d.jl @@ -81,27 +81,33 @@ function create_contact_segmentation(problem::Problem{Contact2D}, slave_element: x1 = slave_element("geometry", time) if deformed && haskey(slave_element, "displacement") - x1 += slave_element("displacement", time) + x1 = map(+, x1, slave_element("displacement", time)) end + slave_midpoint = 1/2*(x1[1] + x1[2]) + slave_length = norm(x1[2]-x1[1]) + for master_element in master_elements x2 = master_element("geometry", time) if deformed && haskey(master_element, "displacement") - x2 += master_element("displacement", time) + x2 = map(+, x2, master_element("displacement", time)) end - if norm(mean(x1) - x2[1]) / norm(x1[2] - x1[1]) > max_distance - continue - end - if norm(mean(x1) - x2[2]) / norm(x1[2] - x1[1]) > max_distance + master_midpoint = 1/2*(x2[1] + x2[2]) + master_length = norm(x2[2]-x2[1]) + # charasteristic length + dist = norm(slave_midpoint - master_midpoint) + cl = dist/max(slave_length, master_length) + if cl > max_distance continue end # 3.1 calculate segmentation - xi1a = project_from_master_to_slave(slave_element, x2[1], time) - xi1b = project_from_master_to_slave(slave_element, x2[2], time) + x21, x22 = x2 + xi1a = project_from_master_to_slave(slave_element, x21, time) + xi1b = project_from_master_to_slave(slave_element, x22, time) xi1 = clamp.([xi1a; xi1b], -1.0, 1.0) l = 1/2*abs(xi1[2]-xi1[1]) if isapprox(l, 0.0) @@ -123,9 +129,9 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assem field_dim = get_unknown_field_dimension(problem) field_name = get_parent_field_name(problem) - # 1. calculate nodal normals and tangents for slave element nodes j ∈ S - normals, tangents = calculate_normals(slave_elements, time, Val{1}; - rotate_normals=props.rotate_normals) + # 1. calculate nodal normals for slave element nodes j ∈ S + normals = calculate_normals([slave_elements...], time; rotate_normals=props.rotate_normals) + tangents = Dict(j => [t2, -t1] for (j, (t1,t2)) in normals) update!(slave_elements, "normal", time => normals) update!(slave_elements, "tangent", time => tangents) @@ -136,15 +142,14 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assem nsl = length(slave_element) X1 = slave_element("geometry", time) + x1 = slave_element("geometry", time) if haskey(slave_element, "displacement") u1 = slave_element("displacement", time) - else - u1 = (zeros(2), zeros(2)) + x1 = map(+, x1, u1) end - x1 = map(+, X1, u1) la1 = slave_element("lambda", time) n1 = slave_element("normal", time) - t1 = slave_element("tangent", time) + t1 = [n1[2], -n1[1]] contact_area = 0.0 contact_error = 0.0 @@ -179,12 +184,11 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assem nm = length(master_element) X2 = master_element("geometry", time) + x2 = master_element("geometry", time) if haskey(master_element, "displacement") u2 = master_element("displacement", time) - else - u2 = (zeros(2), zeros(2)) + x2 = map(+, x2, u2) end - x2 = map(+, X2, u2) # 3.3. loop integration points of one integration segment and calculate # local mortar matrices @@ -277,7 +281,8 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assem la = problem.assembly.la # FIXME: for matrix operations, we need to know the dimensions of the # final matrices - ndofs = 0 + ndofs = 2*length(S) + ndofs = max(ndofs, length(la)) ndofs = max(ndofs, size(problem.assembly.K, 2)) ndofs = max(ndofs, size(problem.assembly.C1, 2)) ndofs = max(ndofs, size(problem.assembly.C2, 2)) @@ -380,9 +385,9 @@ function FEMBase.assemble_elements!(problem::Problem{Contact2D}, assembly::Assem update!(slave_elements, "slip nodes", time => is_slip) end - info("# | active | inactive | stick | slip | gap | pres | comp") + info("# | A | I | St | Sl | gap | pres | comp") for j in S - str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | " + str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | " str2 = "$(round(weighted_gap[j][1], 3)) | $(round(contact_pressure[j][1], 3)) | $(round(complementarity_condition[j][1], 3))" info(str1 * str2) end diff --git a/src/mortar2d.jl b/src/mortar2d.jl index 9e291fb..c3ec350 100644 --- a/src/mortar2d.jl +++ b/src/mortar2d.jl @@ -52,104 +52,89 @@ function get_master_dofs(problem::Problem{Mortar2D}) return sort(unique(dofs)) end -function newton(f, df, x; tol=1.0e-6, max_iterations=10) - for i=1:max_iterations - dx = -f(x)/df(x) - x += dx - if norm(dx) < tol - return x - end - end - error("Newton iteration did not converge in $max_iterations iterations") -end - -function cross2(a, b) - cross([a; 0], [b; 0])[3] -end +function project_from_master_to_slave(slave_element::Element{E}, x2, time; + tol=1.0e-9, max_iterations=5) where {E<:MortarElements2D} -function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time) x1_ = slave_element("geometry", time) n1_ = slave_element("normal", time) + cross2(a, b) = cross([a; 0.0], [b; 0.0])[3] x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_) dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_) n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_) dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_) R(xi1) = cross2(x1(xi1)-x2, n1(xi1)) dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1)) - xi1 = nothing - try - xi1 = newton(R, dR, 0.0) - catch - warn("projection from master to slave failed with following arguments:") - warn("slave element x1: $x1_") - warn("slave element n1: $n1_") - warn("master element x2: $x2") - warn("time: $time") - len = norm(x1_[2] - x1_[1]) - midpnt = mean(x1_) - dist = norm(midpnt - x2) - distval = dist/len - warn("midpoint of slave element: $midpnt") - warn("length of slave element: $len") - warn("distance between midpoint of slave element and x2: $dist") - warn("charasteristic measure: $distval") - rethrow() + + xi1 = 0.0 + dxi1 = 0.0 + for i=1:max_iterations + dxi1 = -R(xi1)/dR(xi1) + xi1 += dxi1 + if norm(dxi1) < tol + return xi1 + end end - return xi1 + + len = norm(x1_[2] - x1_[1]) + midpnt = mean(x1_) + dist = norm(midpnt - x2) + distval = dist/len + info("Projection from the master element to the slave failed.") + info("Arguments:") + info("Time: $time") + info("Point to project to the slave element: (x2): $x2") + info("Slave element geometry (x1): $x1_") + info("Slave element normal direction (n1): $n1_") + info("Midpoint of slave element: $midpnt") + info("Length of slave element: $len") + info("Distance between midpoint of slave element and x2: $dist") + info("Distance/Lenght-ratio: $distval") + info("Number of iterations: $max_iterations") + info("Wanted convergence tolerance: $tol") + info("Achieved convergence tolerance: $(norm(dxi1))") + error("Failed to project point to slave surface in $max_iterations.") end -function project_from_slave_to_master{E<:MortarElements2D}(master_element::Element{E}, x1, n1, time) +function project_from_slave_to_master(master_element::Element{E}, x1, n1, time; + max_iterations=5, tol=1.0e-9) where {E<:MortarElements2D} x2_ = master_element("geometry", time) x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_) dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_) cross2(a, b) = cross([a; 0], [b; 0])[3] R(xi2) = cross2(x2(xi2)-x1, n1) dR(xi2) = cross2(dx2(xi2), n1) - xi2 = newton(R, dR, 0.0) - return xi2 -end - -function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false) - tangents = Dict{Int64, Vector{Float64}}() - for element in elements - conn = get_connectivity(element) - #X1 = element("geometry", time) - #dN = get_dbasis(element, [0.0], time) - #tangent = vec(sum([kron(dN[:,i], X1[i]') for i=1:length(X1)])) - tangent = vec(element([0.0], time, Val{:Jacobian})) - for nid in conn - if haskey(tangents, nid) - tangents[nid] += tangent - else - tangents[nid] = tangent - end + xi2 = 0.0 + dxi2 = 0.0 + for i=1:max_iterations + dxi2 = -R(xi2)/dR(xi2) + xi2 += dxi2 + if norm(dxi2) < tol + return xi2 end end + error("Failed to project point to master surface in $max_iterations iterations.") +end + + +function calculate_normals(elements::Vector{Element{Seg2}}, time; rotate_normals=false) - Q = [0.0 -1.0; 1.0 0.0] normals = Dict{Int64, Vector{Float64}}() - S = collect(keys(tangents)) - for j in S - tangents[j] /= norm(tangents[j]) - normals[j] = Q*tangents[j] - end - if rotate_normals - for j in S - normals[j] = -normals[j] + for element in elements + X1, X2 = element("geometry", time) + t1, t2 = (X2-X1)/norm(X2-X1) + normal = [-t2, t1] + for j in get_connectivity(element) + normals[j] = get(normals, j, zeros(2)) + normal end end - return normals, tangents -end - -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]) - update!(element, "tangent", time => [tangents[j] for j in conn]) + s = rotate_normals ? -1.0 : 1.0 + for j in keys(normals) + normals[j] = s*normals[j]/norm(normals[j]) end + + return normals end @@ -162,11 +147,9 @@ function FEMBase.assemble_elements!(problem::Problem{Mortar2D}, assembly::Assemb field_dim = get_unknown_field_dimension(problem) field_name = get_parent_field_name(problem) - # 1. calculate nodal normals and tangents for slave element nodes j ∈ S - normals, tangents = calculate_normals(slave_elements, time, Val{1}; - rotate_normals=false) + # 1. calculate nodal normals for slave element nodes j ∈ S + normals = calculate_normals([slave_elements...], time) update!(slave_elements, "normal", time => normals) - update!(slave_elements, "tangent", time => tangents) # 2. loop all slave elements for slave_element in slave_elements @@ -256,10 +239,7 @@ function get_mortar_matrix_P(problem::Problem{Mortar2D}) end """ Eliminate mesh tie constraints from matrices K, M, f. """ -function FEMBase.eliminate_boundary_conditions!(problem::Problem{Mortar2D}, - K::SparseMatrixCSC, - M::SparseMatrixCSC, - f::SparseVector) +function FEMBase.eliminate_boundary_conditions!(problem::Problem{Mortar2D}, K, M, f) m = get_master_dofs(problem) s = get_slave_dofs(problem) @@ -269,8 +249,9 @@ function FEMBase.eliminate_boundary_conditions!(problem::Problem{Mortar2D}, Id[s] = 0.0 Q = spdiagm(Id) Q[m,s] += P' - K[:,:] = Q*K_red*Q' - M[:,:] = Q*M_red*Q' + K[:,:] = Q*K*Q' + M[:,:] = Q*M*Q' + f[:] = Q*f return true end diff --git a/test/runtests.jl b/test/runtests.jl index 501c564..40232ab 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,4 +4,26 @@ using MortarContact2D using Base.Test -@testset "Mortar coupling, example 1" begin include("test_mortar2d_ex1.jl") end +@testset "test MortarContact2D.jl" begin + + @testset "Projecting vertices between surfaces" begin + include("test_mortar2d_calculate_projection.jl") + end + + @testset "Mortar coupling, example 1" begin + include("test_mortar2d_ex1.jl") + end + + @testset "Mortar coupling, example 2" begin + include("test_mortar2d_ex2.jl") + end + + @testset "Contact basic test" begin + include("test_contact2d.jl") + end + + @testset "Testing contact segmentation" begin + include("test_contact_segmentation.jl") + end + +end diff --git a/test/test_contact2d.jl b/test/test_contact2d.jl new file mode 100644 index 0000000..ca2b1ef --- /dev/null +++ b/test/test_contact2d.jl @@ -0,0 +1,134 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE + +using MortarContact2D +using MortarContact2D: get_slave_dofs, get_master_dofs +using Base.Test + +# Original problem: +# +# using JuliaFEM +# using JuliaFEM.Preprocess +# using JuliaFEM.Testing +# +# mesh = Mesh() +# add_node!(mesh, 1, [0.0, 0.0]) +# add_node!(mesh, 2, [1.0, 0.0]) +# add_node!(mesh, 3, [1.0, 0.5]) +# add_node!(mesh, 4, [0.0, 0.5]) +# gap = 0.1 +# add_node!(mesh, 5, [0.0, 0.5+gap]) +# add_node!(mesh, 6, [1.0, 0.5+gap]) +# add_node!(mesh, 7, [1.0, 1.0+gap]) +# add_node!(mesh, 8, [0.0, 1.0+gap]) +# add_element!(mesh, 1, :Quad4, [1, 2, 3, 4]) +# add_element!(mesh, 2, :Quad4, [5, 6, 7, 8]) +# add_element!(mesh, 3, :Seg2, [1, 2]) +# add_element!(mesh, 4, :Seg2, [7, 8]) +# add_element!(mesh, 5, :Seg2, [4, 3]) +# add_element!(mesh, 6, :Seg2, [6, 5]) +# add_element_to_element_set!(mesh, :LOWER, 1) +# add_element_to_element_set!(mesh, :UPPER, 2) +# add_element_to_element_set!(mesh, :LOWER_BOTTOM, 3) +# add_element_to_element_set!(mesh, :UPPER_TOP, 4) +# add_element_to_element_set!(mesh, :LOWER_TOP, 5) +# add_element_to_element_set!(mesh, :UPPER_BOTTOM, 6) +# upper = Problem(Elasticity, "UPPER", 2) +# upper.properties.formulation = :plane_stress +# upper.elements = create_elements(mesh, "UPPER") +# update!(upper.elements, "youngs modulus", 288.0) +# update!(upper.elements, "poissons ratio", 1/3) +# lower = Problem(Elasticity, "LOWER", 2) +# lower.properties.formulation = :plane_stress +# lower.elements = create_elements(mesh, "LOWER") +# update!(lower.elements, "youngs modulus", 288.0) +# update!(lower.elements, "poissons ratio", 1/3) +# bc_upper = Problem(Dirichlet, "UPPER_TOP", 2, "displacement") +# bc_upper.elements = create_elements(mesh, "UPPER_TOP") +# #update!(bc_upper.elements, "displacement 1", -17/90) +# #update!(bc_upper.elements, "displacement 1", -17/90) +# update!(bc_upper.elements, "displacement 1", -0.2) +# update!(bc_upper.elements, "displacement 2", -0.2) +# bc_lower = Problem(Dirichlet, "LOWER_BOTTOM", 2, "displacement") +# bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM") +# update!(bc_lower.elements, "displacement 1", 0.0) +# update!(bc_lower.elements, "displacement 2", 0.0) +# contact = Problem(Contact2D, "LOWER_TO_UPPER", 2, "displacement") +# contact_slave_elements = create_elements(mesh, "LOWER_TOP") +# contact_master_elements = create_elements(mesh, "UPPER_BOTTOM") +# add_slave_elements!(contact, contact_slave_elements) +# add_master_elements!(contact, contact_master_elements) +# solver = Solver(Nonlinear) +# push!(solver, upper, lower, bc_upper, bc_lower, contact) +# solver() +# master = first(contact_master_elements) +# slave = first(contact_slave_elements) +# um = master("displacement", (0.0,), 0.0) +# us = slave("displacement", (0.0,), 0.0) +# la = slave("lambda", (0.0,), 0.0) +# info("um = $um, us = $us, la = $la") +# @test isapprox(um, [-0.20, -0.15]) +# @test isapprox(us, [0.0, -0.05]) +# @test isapprox(la, [0.0, 30.375]) + +# Equivalent test problem +slave = Element(Seg2, [1, 2]) +update!(slave, "geometry", ([0.0, 0.5], [1.0, 0.5])) +update!(slave, "displacement", 0.0 => (zeros(2), zeros(2))) +update!(slave, "lambda", 0.0 => (zeros(2), zeros(2))) +update!(slave, "displacement", 1.0 => ([-0.01875, -0.05], [0.01875, -0.05])) +update!(slave, "lambda", 1.0 => ([0.0, 30.375], [0.0, 30.375])) +master = Element(Seg2, [3, 4]) +update!(master, "geometry", ([1.0, 0.6], [0.0, 0.6])) +update!(master, "displacement", 0.0 => (zeros(2), zeros(2))) +update!(master, "displacement", 1.0 => ([-0.18125, -0.15], [-0.21875, -0.15])) +problem = Problem(Contact2D, "LOWER_TO_UPPER", 2, "displacement") +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +um = master("displacement", (0.0,), 1.0) +us = slave("displacement", (0.0,), 1.0) +la = slave("lambda", (0.0,), 1.0) +@test isapprox(um, [-0.20, -0.15]) +@test isapprox(us, [0.0, -0.05]) +@test isapprox(la, [0.0, 30.375]) + +assemble!(problem, 0.0) +@test isempty(full(problem.assembly.K)) +@test isempty(full(problem.assembly.C1)) +@test isempty(full(problem.assembly.C2)) +@test isempty(full(problem.assembly.D)) +@test isempty(full(problem.assembly.f)) +@test isempty(full(problem.assembly.g)) + +empty!(problem.assembly) +assemble!(problem, 1.0) +@test isempty(full(problem.assembly.K)) +@test isempty(full(problem.assembly.f)) +C1 = full(problem.assembly.C1, 4, 8) +C2 = full(problem.assembly.C2, 4, 8) +D = full(problem.assembly.D, 4, 8) +g = full(problem.assembly.g, 4, 1) +C1_expected = [ + 0.5 0.0 0.0 0.0 0.0 0.0 -0.5 0.0 + 0.0 0.5 0.0 0.0 0.0 0.0 0.0 -0.5 + 0.0 0.0 0.5 0.0 -0.5 0.0 0.0 0.0 + 0.0 0.0 0.0 0.5 0.0 -0.5 0.0 0.0] +C2_expected = [ + 0.0 0.5 0.0 0.0 0.0 0.0 0.0 -0.5 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 + 0.0 0.0 0.0 0.5 0.0 -0.5 0.0 0.0 + 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] +D_expected = [ + 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 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] +g_expected = zeros(4) +@test isapprox(C1, C1_expected) +@test isapprox(C2, C2_expected) +@test isapprox(D, D_expected) +@test isapprox(g, g_expected; atol=1.0e-12) + +@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] diff --git a/test/test_contact_segmentation.jl b/test/test_contact_segmentation.jl new file mode 100644 index 0000000..d208950 --- /dev/null +++ b/test/test_contact_segmentation.jl @@ -0,0 +1,82 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE + +using MortarContact2D +using MortarContact2D: create_contact_segmentation, calculate_normals + +using FEMBase.Test + +X = Dict(1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [0.0, 1.0], + 4 => [1.0, 1.0]) + +slave = Element(Seg2, [1, 2]) +master = Element(Seg2, [3, 4]) +elements = [slave, master] +update!(elements, "geometry", X) +update!(slave, "displacement", 0.0 => ([0.0,0.0], [0.0,0.0])) +update!(master, "displacement", 0.0 => ([0.0,0.0], [0.0,0.0])) +update!(slave, "displacement", 1.0 => ([0.0,0.0], [0.0,0.0])) +update!(master, "displacement", 1.0 => ([2.0,0.0], [2.0,0.0])) +problem = Problem(Contact2D, "test problem", 2, "displacement") +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +normals = calculate_normals([slave], 0.0) +update!([slave], "normal", 0.0 => normals) +seg1 = create_contact_segmentation(problem, slave, [master], 0.0; deformed=true) +seg2 = create_contact_segmentation(problem, slave, [master], 1.0; deformed=true) +problem.properties.max_distance = 0.0 +seg3 = create_contact_segmentation(problem, slave, [master], 0.0; deformed=false) +@test length(seg1) == 1 +@test length(seg2) == 0 +@test length(seg3) == 0 + +slave = Element(Seg2, [1, 2]) +master = Element(Seg2, [3, 4]) +elements = [slave, master] +update!(elements, "geometry", X) +problem = Problem(Contact2D, "test problem", 2, "displacement") +problem.assembly.la = zeros(8) +problem.properties.store_fields = [ + "contact area", + "contact error", + "weighted gap", + "contact pressure", + "complementarity condition", + "active nodes", + "inactive nodes", + "stick nodes", + "slip nodes"] +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +assemble!(problem, 0.0) + +X[3] += [2.0, 0.0] +X[4] += [2.0, 0.0] +slave = Element(Seg2, [1, 2]) +master = Element(Seg2, [3, 4]) +elements = [slave, master] +update!(elements, "geometry", X) +problem = Problem(Contact2D, "test problem", 2, "displacement") +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +assemble!(problem, 0.0) + +X = Dict(1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [0.0, 0.0], + 4 => [1.0, 0.0]) + +slave = Element(Seg2, [1, 2]) +master = Element(Seg2, [3, 4]) +elements = [slave, master] +update!(elements, "geometry", X) +problem = Problem(Contact2D, "test problem", 2, "displacement") +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +assemble!(problem, 0.0) +empty!(problem.assembly) +problem.properties.iteration = 0 +problem.properties.contact_state_in_first_iteration = :INACTIVE +assemble!(problem, 0.0) diff --git a/test/test_mortar2d_calculate_projection.jl b/test/test_mortar2d_calculate_projection.jl index 606ef1f..3029604 100644 --- a/test/test_mortar2d_calculate_projection.jl +++ b/test/test_mortar2d_calculate_projection.jl @@ -1,94 +1,81 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE -using JuliaFEM -using JuliaFEM.Testing -using JuliaFEM: calculate_normals +using MortarContact2D +using MortarContact2D: calculate_normals, + project_from_master_to_slave, + project_from_slave_to_master -function get_test_2d_model() - X = Dict( - 7 => [0.0, 1.0], - 8 => [5/4, 1.0], - 9 => [2.0, 1.0], - 10 => [0.0, 1.0], - 11 => [3/4, 1.0], - 12 => [2.0, 1.0]) - mel1 = Element(Seg2, [7, 8]) - mel2 = Element(Seg2, [8, 9]) - sel1 = Element(Seg2, [10, 11]) - sel2 = Element(Seg2, [11, 12]) - update!([mel1, mel2, sel1, sel2], "geometry", X) - update!([sel1, sel2], "master elements", [sel1, sel2]) - slave_elements = [sel1, sel2] - time = 0.0 - normals, tangents = calculate_normals(slave_elements, time, Val{1}) - update!(slave_elements, "normal", time => normals) - return [sel1, sel2], [mel1, mel2] -end +using Base.Test -@testset "calculate flat 2d projection from slave to master" begin - (sel1, sel2), (mel1, mel2) = get_test_2d_model() +X = Dict( + 1 => [0.0, 1.0], + 2 => [5/4, 1.0], + 3 => [2.0, 1.0], + 4 => [0.0, 1.0], + 5 => [3/4, 1.0], + 6 => [2.0, 1.0]) +mel1 = Element(Seg2, [1, 2]) +mel2 = Element(Seg2, [2, 3]) +sel1 = Element(Seg2, [4, 5]) +sel2 = Element(Seg2, [5, 6]) +update!([mel1, mel2, sel1, sel2], "geometry", X) +master_elements = [mel1, mel2] +slave_elements = [sel1, sel2] +normals = calculate_normals(slave_elements, 0.0) +update!(slave_elements, "normal", 0.0 => normals) - time = 0.0 - X1 = sel1("geometry", [-1.0], time) - n1 = sel1("normal", [-1.0], time) - println("X1 = ", X1) - println("n1 = ", n1) - xi2 = project_from_slave_to_master(mel1, X1, n1, time) - @test isapprox(xi2, -1.0) +X1 = [0.0, 1.0] +n1 = [0.0, 1.0] +xi2 = project_from_slave_to_master(mel1, X1, n1, 0.0) +@test isapprox(xi2, -1.0) - X1 = sel1("geometry", [1.0], time) - n1 = sel1("normal", [1.0], time) - xi2 = project_from_slave_to_master(mel1, X1, n1, time) - @test isapprox(xi2, 0.2) +X1 = [0.75, 1.0] +n1 = [0.00, 1.0] +xi2 = project_from_slave_to_master(mel1, X1, n1, 0.0) +@test isapprox(xi2, 0.2) - X2 = mel1("geometry", xi2, time) - @test isapprox(X2, [3/4, 1.0]) -end +X2 = mel1("geometry", xi2, 0.0) +@test isapprox(X2, [3/4, 1.0]) -@testset "calculate flat 2d projection from master to slave" begin - (sel1, sel2), (mel1, mel2) = get_test_2d_model() - time = 0.0 - x2 = mel1("geometry", [-1.0], time) - xi1 = project_from_master_to_slave(sel1, x2, time) - @test isapprox(xi1, -1.0) - x2 = mel1("geometry", [1.0], time) - xi1 = project_from_master_to_slave(sel1, x2, time) - X1 = sel1("geometry", xi1, time) - @test isapprox(X1, [5/4, 1.0]) -end +X2 = [0.0, 1.0] +xi1 = project_from_master_to_slave(sel1, X2, 0.0) +@test isapprox(xi1, -1.0) +X2 = [1.25, 1.00] +xi1 = project_from_master_to_slave(sel1, X2, 0.0) +X1 = sel1("geometry", xi1, 0.0) +@test isapprox(X1, [5/4, 1.0]) -@testset "calculate flat 2d projection rotated 90 degrees" begin - X = Dict{Int64, Vector{Float64}}( - 1 => [0.0, 0.0], - 2 => [0.0, 1.0], - 3 => [0.0, 1.0], - 4 => [0.0, 0.0]) - sel1 = Element(Seg2, [1, 2]) - mel1 = Element(Seg2, [3, 4]) - update!([sel1, mel1], "geometry", X) - time = 0.0 - slave_elements = [sel1] - time = 0.0 - normals, tangents = calculate_normals(slave_elements, time, Val{1}) - update!(slave_elements, "normal", time => normals) +X = Dict( + 1 => [0.0, 0.0], + 2 => [0.0, 1.0], + 3 => [0.0, 1.0], + 4 => [0.0, 0.0]) - X2 = mel1("geometry", [-1.0], time) - xi = project_from_master_to_slave(sel1, X2, time) - @test isapprox(xi, 1.0) +sel1 = Element(Seg2, [1, 2]) +mel1 = Element(Seg2, [3, 4]) +update!([sel1, mel1], "geometry", X) +slave_elements = [sel1] +normals = calculate_normals(slave_elements, 0.0) +update!(slave_elements, "normal", 0.0 => normals) - X2 = mel1("geometry", [1.0], time) - xi = project_from_master_to_slave(sel1, X2, time) - @test isapprox(xi, -1.0) +X2 = mel1("geometry", (-1.0, ), 0.0) +xi = project_from_master_to_slave(sel1, X2, 0.0) +@test isapprox(xi, 1.0) - X1 = sel1("geometry", [-1.0], time) - n1 = sel1("normal", [-1.0], time) - xi = project_from_slave_to_master(mel1, X1, n1, time) - @test isapprox(xi, 1.0) +X2 = mel1("geometry", (1.0, ), 0.0) +xi = project_from_master_to_slave(sel1, X2, 0.0) +@test isapprox(xi, -1.0) - X1 = sel1("geometry", [1.0], time) - n1 = sel1("normal", [1.0], time) - xi = project_from_slave_to_master(mel1, X1, n1, time) - @test isapprox(xi, -1.0) -end +X1 = sel1("geometry", (-1.0, ), 0.0) +n1 = sel1("normal", (-1.0, ), 0.0) +xi = project_from_slave_to_master(mel1, X1, n1, 0.0) +@test isapprox(xi, 1.0) +X1 = sel1("geometry", (1.0, ), 0.0) +n1 = sel1("normal", (1.0, ), 0.0) +xi = project_from_slave_to_master(mel1, X1, n1, 0.0) +@test isapprox(xi, -1.0) + +@test_throws Exception project_from_master_to_slave(sel1, [0.0, 0.0], 0.0; tol=0.0) +@test_throws Exception project_from_slave_to_master(mel1, [0.0, 1.0], [-1.0, 0.0], 0.0; tol=0.0) diff --git a/test/test_mortar2d_ex1.jl b/test/test_mortar2d_ex1.jl index 89c8c7e..f4bf7f1 100644 --- a/test/test_mortar2d_ex1.jl +++ b/test/test_mortar2d_ex1.jl @@ -9,19 +9,19 @@ using MortarContact2D: get_slave_dofs, get_master_dofs, get_mortar_matrix_P using Base.Test -#' Slave side coordinates: +# Slave side coordinates: Xs = Dict( 1 => [0.0, 1.0], 2 => [1.0, 1.0], 3 => [2.0, 1.0]) -#' Master side coordinates: +# Master side coordinates: Xm = Dict( 4 => [0.0, 1.0], 5 => [5/4, 1.0], 6 => [2.0, 1.0]) -#' Define elements, update geometry: +# Define elements, update geometry: X = merge(Xm , Xs) es1 = Element(Seg2, [1, 2]) es2 = Element(Seg2, [2, 3]) @@ -30,13 +30,13 @@ em2 = Element(Seg2, [5, 6]) elements = [es1, es2, em1, em2] update!(elements, "geometry", X) -#' Define new problem `Mortar2D`, coupling elements using mortar methods +# Define new problem `Mortar2D`, coupling elements using mortar methods problem = Problem(Mortar2D, "test interface", 1, "u") @test_throws ErrorException add_elements!(problem, [es1, es2]) add_slave_elements!(problem, [es1, es2]) add_master_elements!(problem, [em1, em2]) -#' Assemble problem +# Assemble problem assemble!(problem, 0.0) s = get_slave_dofs(problem) m = get_master_dofs(problem) @@ -52,7 +52,7 @@ D_expected = [1/3 1/6 0; 1/6 2/3 1/6; 0 1/6 1/3] M_expected = [11/30 2/15 0; 41/160 13/20 3/32; 1/480 13/60 9/32] P_expected = 1/320*[329 -24 15; 46 304 -30; -21 56 285] -#' Results: +# Results: @test isapprox(full(D), D_expected) @test isapprox(full(M), M_expected) diff --git a/test/test_mortar2d_ex2.jl b/test/test_mortar2d_ex2.jl new file mode 100644 index 0000000..6be14d8 --- /dev/null +++ b/test/test_mortar2d_ex2.jl @@ -0,0 +1,55 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/MortarContact2D.jl/blob/master/LICENSE + +# Construct discrete coupling operator between two elements and eliminate +# boundary conditions from global matrix assembly. Thesis, p. 73. + +using MortarContact2D +using Base.Test + +X = Dict( + 1 => [1.0, 2.0], + 2 => [3.0, 2.0], + 3 => [2.0, 2.0], + 4 => [0.0, 2.0], + 5 => [3.0, 4.0], + 6 => [1.0, 4.0], + 7 => [0.0, 0.0], + 8 => [2.0, 0.0]) + +slave = Element(Seg2, [1, 2]) +master = Element(Seg2, [3, 4]) +elements = [slave, master] +update!(elements, "geometry", X) + +problem = Problem(Mortar2D, "test interface", 2, "displacement") +add_slave_elements!(problem, [slave]) +add_master_elements!(problem, [master]) +assemble!(problem, 0.0) + +K_SS = K_MM = [ + 180.0 81.0 -126.0 27.0 + 81.0 180.0 -27.0 36.0 +-126.0 -27.0 180.0 -81.0 + 27.0 36.0 -81.0 180.0] +M_SS = M_MM = [ + 40.0 0.0 20.0 0.0 + 0.0 40.0 0.0 20.0 + 20.0 0.0 40.0 0.0 + 0.0 20.0 0.0 40.0] +f_S = [-27.0, -81.0, 81.0, -135.0] +f_M = [ 0.0, 0.0, 0.0, 0.0] +K = [K_SS zeros(4,4); zeros(4,4) K_MM] +M = [M_SS zeros(4,4); zeros(4,4) M_MM] +f = [f_S; f_M] + +eliminate_boundary_conditions!(problem, K, M, f) +K_expected = [ + 441.0 -81.0 -279.0 81.0 + -81.0 684.0 81.0 -36.0 + -279.0 81.0 333.0 -81.0 + 81.0 -36.0 -81.0 252.0] +f_expected = [108.0, -243.0, -54.0, 27.0] + +@test isapprox(K[5:8,5:8], K_expected) +@test isapprox(f[5:8], f_expected)