mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-30 16:12:51 +00:00
removed obsolete code
This commit is contained in:
Jukka Aho
@@ -86,6 +86,7 @@ export calculate_normals,
|
||||
|
||||
### Mortar methods, contact mechanics extension ###
|
||||
include("problems_contact.jl")
|
||||
include("problems_contact_2d.jl")
|
||||
export Contact
|
||||
|
||||
module API
|
||||
|
||||
@@ -37,8 +37,6 @@ function get_formulation_type(problem::Problem{Contact})
|
||||
return :incremental
|
||||
end
|
||||
|
||||
typealias ContactElements2D Union{Seg2}
|
||||
|
||||
function assemble!(problem::Problem{Contact}, time::Real)
|
||||
if problem.properties.dimension == -1
|
||||
problem.properties.dimension = dim = size(first(problem.elements), 1)
|
||||
@@ -52,247 +50,3 @@ function assemble!(problem::Problem{Contact}, time::Real)
|
||||
assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff)
|
||||
end
|
||||
|
||||
""" Frictionless 2d small sliding contact without forwarddiff. """
|
||||
function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
::Type{Val{1}}, ::Type{Val{false}},
|
||||
::Type{Val{false}}, ::Type{Val{false}}; debug=false)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(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, Val{1};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
update!(slave_elements, "tangent", time => tangents)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element["geometry"](time)
|
||||
u1 = slave_element["displacement"](time)
|
||||
la1 = slave_element["reaction force"](time)
|
||||
n1 = slave_element["normal"](time)
|
||||
t1 = slave_element["tangent"](time)
|
||||
x1 = X1 + u1
|
||||
Q1_ = [n1[1] t1[1]]
|
||||
Q2_ = [n1[2] t1[2]]
|
||||
Z = zeros(2, 2)
|
||||
Q2 = [Q1_ Z; Z Q2_]
|
||||
contact_area = 0.0
|
||||
contact_error = 0.0
|
||||
|
||||
if "element area" in props.store_fields
|
||||
element_area = 0.0
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
element_area += w
|
||||
end
|
||||
update!(slave_element, "element area", time => element_area)
|
||||
end
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
nm = length(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
|
||||
# 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)
|
||||
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
|
||||
|
||||
# 3.2. bi-orthogonal basis
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
Ae = zeros(nsl, nsl)
|
||||
if props.dual_basis
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
De += w*diagm(N1)
|
||||
Me += w*N1*N1'
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
else
|
||||
Ae = eye(nsl)
|
||||
end
|
||||
|
||||
# 3.3. loop integration points of one integration segment and calculate
|
||||
# local mortar matrices
|
||||
fill!(De, 0.0)
|
||||
fill!(Me, 0.0)
|
||||
ge = zeros(field_dim*nsl)
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
Phi = Ae*N1
|
||||
# project gauss point from slave element to master element in direction n_s
|
||||
X_s = N1*X1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
X_m = N2*X2
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
x_s = X_s + u_s
|
||||
x_m = X_m + u_m
|
||||
la_s = Phi*la1
|
||||
ge += w*vec((x_m-x_s)*Phi')
|
||||
contact_area += w
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end
|
||||
|
||||
# add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D2 = zeros(nsldofs, nsldofs)
|
||||
M2 = zeros(nmdofs, nmdofs)
|
||||
for i=1:field_dim
|
||||
D2[i:field_dim:end, i:field_dim:end] += De
|
||||
M2[i:field_dim:end, i:field_dim:end] += Me
|
||||
end
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D2)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M2)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2)
|
||||
add!(problem.assembly.g, sdofs, Q2'*ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
if "contact area" in props.store_fields
|
||||
update!(slave_element, "contact area", time => contact_area)
|
||||
end
|
||||
|
||||
if "contact error" in props.store_fields
|
||||
update!(slave_element, "contact error", time => contact_error)
|
||||
end
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
S = sort(collect(keys(normals))) # slave element nodes
|
||||
weighted_gap = Dict{Int64, Vector{Float64}}()
|
||||
contact_pressure = Dict{Int64, Vector{Float64}}()
|
||||
complementarity_condition = Dict{Int64, Vector{Float64}}()
|
||||
is_active = Dict{Int64, Int}()
|
||||
is_inactive = Dict{Int64, Int}()
|
||||
is_slip = Dict{Int64, Int}()
|
||||
is_stick = Dict{Int64, Int}()
|
||||
|
||||
g = full(problem.assembly.g)
|
||||
la = problem.assembly.la
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
weighted_gap[j] = g[dofs]
|
||||
if length(la) != 0
|
||||
p = dot(normals[j], la[dofs])
|
||||
t = dot(tangents[j], la[dofs])
|
||||
contact_pressure[j] = [p, t]
|
||||
else
|
||||
contact_pressure[j] = [0.0, 0.0]
|
||||
end
|
||||
complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
if complementarity_condition[j][1] < 0
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
is_stick[j] = 0
|
||||
else
|
||||
is_inactive[j] = 0
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
if "weighted gap" in props.store_fields
|
||||
update!(slave_elements, "weighted gap", time => weighted_gap)
|
||||
end
|
||||
if "contact pressure" in props.store_fields
|
||||
update!(slave_elements, "contact pressure", time => contact_pressure)
|
||||
end
|
||||
if "complementarity condition" in props.store_fields
|
||||
update!(slave_elements, "complementarity condition", time => complementarity_condition)
|
||||
end
|
||||
if "active nodes" in props.store_fields
|
||||
update!(slave_elements, "active nodes", time => is_active)
|
||||
end
|
||||
if "inactive nodes" in props.store_fields
|
||||
update!(slave_elements, "inactive nodes", time => is_inactive)
|
||||
end
|
||||
if "stick nodes" in props.store_fields
|
||||
update!(slave_elements, "stick nodes", time => is_stick)
|
||||
end
|
||||
if "slip nodes" in props.store_fields
|
||||
update!(slave_elements, "slip nodes", time => is_slip)
|
||||
end
|
||||
|
||||
info("# | active | inactive | stick | slip | gap | pres | comp")
|
||||
for j in S
|
||||
str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
|
||||
str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
|
||||
info(str1 * str2)
|
||||
end
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
ndofs = size(C1, 1)
|
||||
C2 = sparse(problem.assembly.C2)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
|
||||
# constitutive modelling in tangent direction, frictionless contact
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if (is_active[j] == 1) && (is_slip[j] == 1)
|
||||
info("$j is in active/slip, removing tangential constraint $(dofs[2])")
|
||||
C2[dofs[2],:] = 0.0
|
||||
g[dofs[2]] = 0.0
|
||||
D[dofs[2], dofs] = tangents[j]
|
||||
end
|
||||
end
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if is_inactive[j] == 1
|
||||
info("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
g[dofs,:] = 0.0
|
||||
end
|
||||
end
|
||||
|
||||
problem.assembly.C1 = C1
|
||||
problem.assembly.C2 = C2
|
||||
problem.assembly.D = D
|
||||
problem.assembly.g = g
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,248 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
typealias ContactElements2D Union{Seg2}
|
||||
|
||||
""" Frictionless 2d small sliding contact without forwarddiff. """
|
||||
function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
::Type{Val{1}}, ::Type{Val{false}},
|
||||
::Type{Val{false}}, ::Type{Val{false}}; debug=false)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(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, Val{1};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
update!(slave_elements, "tangent", time => tangents)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element["geometry"](time)
|
||||
u1 = slave_element["displacement"](time)
|
||||
la1 = slave_element["reaction force"](time)
|
||||
n1 = slave_element["normal"](time)
|
||||
t1 = slave_element["tangent"](time)
|
||||
x1 = X1 + u1
|
||||
Q1_ = [n1[1] t1[1]]
|
||||
Q2_ = [n1[2] t1[2]]
|
||||
Z = zeros(2, 2)
|
||||
Q2 = [Q1_ Z; Z Q2_]
|
||||
contact_area = 0.0
|
||||
contact_error = 0.0
|
||||
|
||||
if "element area" in props.store_fields
|
||||
element_area = 0.0
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
element_area += w
|
||||
end
|
||||
update!(slave_element, "element area", time => element_area)
|
||||
end
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
nm = length(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
|
||||
# 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)
|
||||
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
|
||||
|
||||
# 3.2. bi-orthogonal basis
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
Ae = zeros(nsl, nsl)
|
||||
if props.dual_basis
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
De += w*diagm(N1)
|
||||
Me += w*N1*N1'
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
else
|
||||
Ae = eye(nsl)
|
||||
end
|
||||
|
||||
# 3.3. loop integration points of one integration segment and calculate
|
||||
# local mortar matrices
|
||||
fill!(De, 0.0)
|
||||
fill!(Me, 0.0)
|
||||
ge = zeros(field_dim*nsl)
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
Phi = Ae*N1
|
||||
# project gauss point from slave element to master element in direction n_s
|
||||
X_s = N1*X1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
X_m = N2*X2
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
x_s = X_s + u_s
|
||||
x_m = X_m + u_m
|
||||
la_s = Phi*la1
|
||||
ge += w*vec((x_m-x_s)*Phi')
|
||||
contact_area += w
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end
|
||||
|
||||
# add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D2 = zeros(nsldofs, nsldofs)
|
||||
M2 = zeros(nmdofs, nmdofs)
|
||||
for i=1:field_dim
|
||||
D2[i:field_dim:end, i:field_dim:end] += De
|
||||
M2[i:field_dim:end, i:field_dim:end] += Me
|
||||
end
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D2)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M2)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q2'*D2)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q2'*M2)
|
||||
add!(problem.assembly.g, sdofs, Q2'*ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
if "contact area" in props.store_fields
|
||||
update!(slave_element, "contact area", time => contact_area)
|
||||
end
|
||||
|
||||
if "contact error" in props.store_fields
|
||||
update!(slave_element, "contact error", time => contact_error)
|
||||
end
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
S = sort(collect(keys(normals))) # slave element nodes
|
||||
weighted_gap = Dict{Int64, Vector{Float64}}()
|
||||
contact_pressure = Dict{Int64, Vector{Float64}}()
|
||||
complementarity_condition = Dict{Int64, Vector{Float64}}()
|
||||
is_active = Dict{Int64, Int}()
|
||||
is_inactive = Dict{Int64, Int}()
|
||||
is_slip = Dict{Int64, Int}()
|
||||
is_stick = Dict{Int64, Int}()
|
||||
|
||||
g = full(problem.assembly.g)
|
||||
la = problem.assembly.la
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
weighted_gap[j] = g[dofs]
|
||||
if length(la) != 0
|
||||
p = dot(normals[j], la[dofs])
|
||||
t = dot(tangents[j], la[dofs])
|
||||
contact_pressure[j] = [p, t]
|
||||
else
|
||||
contact_pressure[j] = [0.0, 0.0]
|
||||
end
|
||||
complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
if complementarity_condition[j][1] < 0
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
is_stick[j] = 0
|
||||
else
|
||||
is_inactive[j] = 0
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
if "weighted gap" in props.store_fields
|
||||
update!(slave_elements, "weighted gap", time => weighted_gap)
|
||||
end
|
||||
if "contact pressure" in props.store_fields
|
||||
update!(slave_elements, "contact pressure", time => contact_pressure)
|
||||
end
|
||||
if "complementarity condition" in props.store_fields
|
||||
update!(slave_elements, "complementarity condition", time => complementarity_condition)
|
||||
end
|
||||
if "active nodes" in props.store_fields
|
||||
update!(slave_elements, "active nodes", time => is_active)
|
||||
end
|
||||
if "inactive nodes" in props.store_fields
|
||||
update!(slave_elements, "inactive nodes", time => is_inactive)
|
||||
end
|
||||
if "stick nodes" in props.store_fields
|
||||
update!(slave_elements, "stick nodes", time => is_stick)
|
||||
end
|
||||
if "slip nodes" in props.store_fields
|
||||
update!(slave_elements, "slip nodes", time => is_slip)
|
||||
end
|
||||
|
||||
info("# | active | inactive | stick | slip | gap | pres | comp")
|
||||
for j in S
|
||||
str1 = "$j | $(is_active[j]) | $(is_inactive[j]) | $(is_stick[j]) | $(is_slip[j]) | "
|
||||
str2 = "$(round(weighted_gap[j], 3)) | $(round(contact_pressure[j], 3)) | $(round(complementarity_condition[j], 3))"
|
||||
info(str1 * str2)
|
||||
end
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
ndofs = size(C1, 1)
|
||||
C2 = sparse(problem.assembly.C2)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
|
||||
# constitutive modelling in tangent direction, frictionless contact
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if (is_active[j] == 1) && (is_slip[j] == 1)
|
||||
info("$j is in active/slip, removing tangential constraint $(dofs[2])")
|
||||
C2[dofs[2],:] = 0.0
|
||||
g[dofs[2]] = 0.0
|
||||
D[dofs[2], dofs] = tangents[j]
|
||||
end
|
||||
end
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
if is_inactive[j] == 1
|
||||
info("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
g[dofs,:] = 0.0
|
||||
end
|
||||
end
|
||||
|
||||
problem.assembly.C1 = C1
|
||||
problem.assembly.C2 = C2
|
||||
problem.assembly.D = D
|
||||
problem.assembly.g = g
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
@@ -1,513 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
""" Fast inverse of 3x3 matrix. """
|
||||
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
|
||||
|
||||
""" Project vertex `p` from element surface to auxiliary plane defined
|
||||
with centerpoint `x0` and normal direction `n0`.
|
||||
"""
|
||||
function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
|
||||
return p - dot(p-x0, n0)*n0
|
||||
end
|
||||
|
||||
""" Project vertex `p` from auxiliary plane (x0, n0) back to element surface.
|
||||
|
||||
This requires solving nonlinear system of equations
|
||||
|
||||
f(α,ξ₁,ξ₂) = Nₖξₖ - αn₀ - p = 0
|
||||
|
||||
"""
|
||||
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(E, xi)
|
||||
dbasis(xi) = get_dbasis(E, xi)
|
||||
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
|
||||
norm(ForwardDiff.get_value(dtheta)) < iter_tol && return theta[1:2], theta[3]
|
||||
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 = $(ForwardDiff.get_value(x0)), n0 = $(ForwardDiff.get_value(n0))")
|
||||
info("element geometry: $(ForwardDiff.get_value(x.data))")
|
||||
info("vertex to project: $(ForwardDiff.get_value(p))")
|
||||
info("parameter vector before giving up: $(ForwardDiff.get_value(theta)')")
|
||||
info("increment in parameter vector before giving up: $(ForwardDiff.get_value(dtheta)')")
|
||||
info("norm(dtheta) before giving up: $(ForwardDiff.get_value(norm(dtheta)))")
|
||||
info("f([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(f([0.0, 0.0, 0.0]))')")
|
||||
info("L([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(L([0.0, 0.0, 0.0])))")
|
||||
|
||||
info("iterations:")
|
||||
theta = zeros(3)
|
||||
dtheta = zeros(3)
|
||||
for i=1:max_iterations
|
||||
info("iter $i, theta = $(ForwardDiff.get_value(theta)')")
|
||||
info("f = $(ForwardDiff.get_value(f(theta))')")
|
||||
info("L = $(ForwardDiff.get_value(L(theta)))")
|
||||
# info("L2 = $(ForwardDiff.get_value(L2(theta)))")
|
||||
dtheta = L(theta) * f(theta)
|
||||
info("dtheta = $(ForwardDiff.get_value(dtheta)')")
|
||||
theta -= dtheta
|
||||
end
|
||||
|
||||
error("project_point_to_surface: did not converge in $max_iterations iterations!")
|
||||
end
|
||||
|
||||
""" Test is q inside sm.
|
||||
http://bbs.dartmouth.edu/~fangq/MATH/download/source/Determining%20if%20a%20point%20lies%20on%20the%20interior%20of%20a%20polygon.htm
|
||||
"""
|
||||
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_polygon_clip(xs, xm, n)
|
||||
# objective: search does line xm1 - xm2 clip xs
|
||||
nm = length(xm)
|
||||
ns = length(xs)
|
||||
P = []
|
||||
|
||||
# 1. test is master point inside slave, if yes, add to clip
|
||||
for i=1:nm
|
||||
vertex_inside_polygon(xm[i], xs) && push!(P, xm[i])
|
||||
end
|
||||
|
||||
# 2. test is slave point inside master, if yes, add to clip
|
||||
for i=1:ns
|
||||
vertex_inside_polygon(xs[i], xm) && push!(P, xs[i])
|
||||
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))")
|
||||
vertex_inside_polygon(q, xm) && push!(P, q)
|
||||
end
|
||||
end
|
||||
|
||||
return P
|
||||
end
|
||||
|
||||
""" Divide polygon to cells. """
|
||||
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
|
||||
|
||||
""" Check that polygon P is in CCW order for the direction n. Reorder if not. """
|
||||
function check_orientation!(P, n)
|
||||
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
|
||||
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
|
||||
|
||||
""" Assemble Mortar problem for three-dimensional problems, i.e. for Tri3, Tri6, Quad4, Quad8, Quad9 elements. """
|
||||
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{3}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(x::Vector)
|
||||
|
||||
ndofs = round(Int, length(x)/2) # x = [u; la]
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc = zeros(u)
|
||||
gap = zeros(u)
|
||||
C = zeros(la)
|
||||
all_slave_nodes = Set{Int64}()
|
||||
slave_surface_area = 0.0
|
||||
slave_surface_area_2 = 0.0
|
||||
slave_element_areas = []
|
||||
|
||||
# 1. calculate and average node normals for slave element nodes
|
||||
normal = zeros(u)
|
||||
for element in get_elements(problem)
|
||||
haskey(element, "master elements") || continue
|
||||
conn = get_connectivity(element)
|
||||
push!(all_slave_nodes, conn...)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:,i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
J = transpose(sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]))
|
||||
n = reshape(cross(J[:,1], J[:,2]), 3, 1)
|
||||
normal[:, conn] += ip.weight*n*N
|
||||
end
|
||||
end
|
||||
|
||||
all_slave_nodes = sort(collect(all_slave_nodes))
|
||||
|
||||
# normalize to unit normal
|
||||
for i in all_slave_nodes
|
||||
normal[:,i] /= norm(normal[:,i])
|
||||
end
|
||||
#normal = ForwardDiff.get_value(normal)
|
||||
|
||||
if props.rotate_normals
|
||||
for i=1:size(normal, 2)
|
||||
normal[:,i] = -normal[:,i]
|
||||
end
|
||||
end
|
||||
|
||||
# 2. loop slave elements and find contact segments
|
||||
for slave_element in get_elements(problem)
|
||||
slave_element_area = 0.0
|
||||
haskey(slave_element, "master elements") || continue
|
||||
info("new slave element")
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
if haskey(slave_element, "displacement")
|
||||
u1 -= slave_element("displacement", time)
|
||||
end
|
||||
x1 = X1 + u1
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
#if haskey(slave_element, "reaction force")
|
||||
# la1 -= slave_element("reaction force", time)
|
||||
#end
|
||||
n1 = Field(Vector[normal[:,i] for i in slave_element_nodes])
|
||||
nnodes = size(slave_element, 2)
|
||||
update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data))
|
||||
|
||||
# 2.1. create auxiliary plane (x0, Q)
|
||||
xi = get_reference_element_midpoint(slave_element)
|
||||
N = vec(get_basis(slave_element, xi))
|
||||
x0 = N*x1
|
||||
n0 = N*n1
|
||||
|
||||
# 2.2. project slave nodes to auxiliary plane
|
||||
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"]
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
if haskey(master_element, "displacement")
|
||||
u2 -= master_element("displacement", time)
|
||||
end
|
||||
x2 = X2 + u2
|
||||
|
||||
distance = norm(mean(x2) - mean(x1))
|
||||
distance > props.maximum_distance && continue
|
||||
|
||||
# 3.1. project master nodes to auxiliary plane
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2]
|
||||
|
||||
# 3.2. create polygon clipping on auxiliary plane
|
||||
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)
|
||||
|
||||
# 3.3. loop integration cells one at time
|
||||
for cell in get_cells(P, C0)
|
||||
x_cell = Field(cell)
|
||||
|
||||
# 3.3.1. create dual basis
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = vec(get_basis(Tri3, ip.xi))
|
||||
x_gauss = N*x_cell
|
||||
xi_slave, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, x1, time)
|
||||
N1 = slave_element(xi_slave, time)
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
|
||||
wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
|
||||
De += wC*diagm(vec(N1))
|
||||
Me += wC*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
# 3.3.2 loop integration points of cell and calculate fc and gap
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = vec(get_basis(Tri3, ip.xi))
|
||||
x_gauss = N*x_cell
|
||||
|
||||
# project gauss point back to element surfaces
|
||||
xi_slave, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, x1, time)
|
||||
xi_master, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, x2, time)
|
||||
|
||||
# evaluate shape functions, calculate contact force and gap
|
||||
N1 = vec(get_basis(slave_element, xi_slave))
|
||||
N2 = vec(get_basis(master_element, xi_master))
|
||||
Phi = Ae*N1
|
||||
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
|
||||
wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
|
||||
|
||||
x_s = N1*x1
|
||||
x_m = N2*x2
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
n_s = N1*n1
|
||||
la_s = Phi*la1
|
||||
la_n = dot(n_s, la_s)
|
||||
g_s = x_s-x_m
|
||||
#gn = -dot(n_s, g_s)
|
||||
fc[:,slave_element_nodes] += wC*la_s*N1'
|
||||
fc[:,master_element_nodes] -= wC*la_s*N2'
|
||||
#gap[:,slave_element_nodes] += wC*g_s*N1'
|
||||
#gap[:,master_element_nodes] += wC*g_s*N2'
|
||||
gn = props.gap_sign*dot(n_s, g_s)
|
||||
#gap[1,slave_element_nodes] += wC*gn*Phi'
|
||||
#gap[1,slave_element_nodes] += wC*gn*Phi'
|
||||
#C[:,master_element_nodes] -= wC*u_s*N2'
|
||||
gap[:,slave_element_nodes] = wC*props.gap_sign*g_s*Phi'
|
||||
#gap[:,master_element_nodes] -= wC*(u_s-u_m)*N2'
|
||||
|
||||
slave_surface_area += wC
|
||||
slave_element_area += wC
|
||||
end # done integrating cell
|
||||
|
||||
end # done for all cells in this segment
|
||||
|
||||
end # done all master elements for this slave element
|
||||
push!(slave_element_areas, slave_element_area)
|
||||
|
||||
end # done all slave elements
|
||||
|
||||
# like in 2d, check contact in nodes based on a complementarity condition
|
||||
|
||||
nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
|
||||
info("gap: $nzgap")
|
||||
|
||||
info("size of normal = $(size(normal))")
|
||||
info("size of la = $(size(la))")
|
||||
info("size of C = $(size(C))")
|
||||
info("S = $all_slave_nodes")
|
||||
info("slave surface area: $(ForwardDiff.get_value(slave_surface_area))")
|
||||
info("slave surface area 2: $(ForwardDiff.get_value(slave_surface_area_2))")
|
||||
info("slave element areas:")
|
||||
for (i, a) in enumerate(slave_element_areas)
|
||||
info("element $i, area = $(ForwardDiff.get_value(a))")
|
||||
end
|
||||
for (i, element) in enumerate(get_elements(problem))
|
||||
haskey(element, "master elements") || continue
|
||||
info("element $i geometry: $(element("geometry", time).data)")
|
||||
end
|
||||
|
||||
for (i, j) in enumerate(all_slave_nodes)
|
||||
if j in props.always_inactive
|
||||
info("special node $j always inactive")
|
||||
C[:,j] = la[:,j]
|
||||
continue
|
||||
end
|
||||
n = normal[:,j]
|
||||
I = eye(3)
|
||||
k = indmax([norm(cross(n,I[:,k])) for k in 1:3])
|
||||
t1 = cross(n, I[:,k])/norm(cross(n, I[:,k]))
|
||||
t2 = cross(n, t1)
|
||||
Q = [n t1 t2]
|
||||
la_nt = Q'*la[:,j]
|
||||
gap_nt = Q'*gap[:,j]
|
||||
C[1,j] = gap_nt[1]
|
||||
C[2:3,j] = la_nt[2:3]
|
||||
#C[:,j] -= gap[:,j]
|
||||
#=
|
||||
if lan - gn < 0
|
||||
info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t1)) x $(ForwardDiff.get_value(t2))")
|
||||
C[1,j] = gn
|
||||
C[2,j] = dot(t1, la[:,j])
|
||||
C[3,j] = dot(t2, la[:,j])
|
||||
else
|
||||
C[:,j] = la[:,j]
|
||||
end
|
||||
=#
|
||||
|
||||
end
|
||||
|
||||
#=
|
||||
for (i, j) in enumerate(all_slave_nodes)
|
||||
n = normal[:,j]
|
||||
I = eye(3)
|
||||
k = indmax([norm(cross(n,I[:,k])) for k in 1:3])
|
||||
t1 = cross(n, I[:,k])/norm(cross(n, I[:,k]))
|
||||
t2 = cross(n, t1)
|
||||
Q = [n t1 t2]
|
||||
Ci = ForwardDiff.get_value(Q'*C[:,j])
|
||||
gapi = ForwardDiff.get_value(Q'*gap[:,j])
|
||||
fci = ForwardDiff.get_value(Q'*fc[:,j])
|
||||
lai = ForwardDiff.get_value(Q'*la[:,j])
|
||||
ui = ForwardDiff.get_value(Q'*u[:,j])
|
||||
info("$i/$j: \nC = $Ci, \nf = $fci, \ngap = $gapi, \nla = $lai, \nu = $ui")
|
||||
end
|
||||
=#
|
||||
|
||||
return vec([fc C])
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
A, allresults = ForwardDiff.jacobian(calculate_interface, x,
|
||||
ForwardDiff.AllResults, cache=autodiffcache)
|
||||
b = ForwardDiff.value(allresults)
|
||||
|
||||
# dump(round(A, 3))
|
||||
# dump(round(b, 3)')
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(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]
|
||||
|
||||
function joo(x)
|
||||
nz1 = sort(unique(rowvals(x)))
|
||||
nz2 = sort(unique(rowvals(x')))
|
||||
info("nz1 = $nz1, nz2 = $nz2")
|
||||
dump(round(full(x[nz1,nz2]), 3))
|
||||
end
|
||||
println("K")
|
||||
joo(K)
|
||||
println("C1")
|
||||
joo(C1)
|
||||
println("C2")
|
||||
joo(C2)
|
||||
println("D")
|
||||
joo(D)
|
||||
println("f")
|
||||
joo(f)
|
||||
println("g")
|
||||
joo(g)
|
||||
#=
|
||||
slaves = [101,108,111,112,113,120,123,124,125,126,129,130,149,150,151,152]
|
||||
for j in slaves
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
info("slave node $j, dofs $dofs")
|
||||
info("Stiffness: $(K[dofs,:])")
|
||||
info("force fc: $(C1[dofs,:])")
|
||||
info("constraint: $(C2[dofs,:])")
|
||||
info("lambdas: $(D[dofs,:])")
|
||||
info("f = $(f[dofs]), g = $(g[dofs])")
|
||||
end
|
||||
=#
|
||||
|
||||
empty!(problem.assembly)
|
||||
add!(problem.assembly.K, K)
|
||||
add!(problem.assembly.C1, C1)
|
||||
add!(problem.assembly.C2, C2)
|
||||
add!(problem.assembly.D, D)
|
||||
add!(problem.assembly.f, f)
|
||||
add!(problem.assembly.g, g)
|
||||
|
||||
return problem.assembly
|
||||
end
|
||||
@@ -1,446 +0,0 @@
|
||||
using JuliaFEM.Core: MortarElements2D, DVTI, Assembly
|
||||
|
||||
import JuliaFEM.Core: project_from_master_to_slave, project_from_slave_to_master, assemble!,
|
||||
get_unknown_field_dimension, get_parent_field_name, get_gdofs, find_elements, get_nodes, Field,
|
||||
get_integration_points, get_basis, get_dbasis, add!
|
||||
|
||||
""" Find segment from slave element corresponding to master element nodes.
|
||||
x1_, n1_
|
||||
slave element geometry and normal direction
|
||||
|
||||
x2_ master element nodes to project onto slave
|
||||
"""
|
||||
function project_from_master_to_slave{E<:MortarElements2D}(
|
||||
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector)
|
||||
|
||||
function x1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*x1_
|
||||
end
|
||||
|
||||
function dx1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*x1_
|
||||
end
|
||||
|
||||
function n1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*n1_
|
||||
end
|
||||
|
||||
function dn1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*n1_
|
||||
end
|
||||
|
||||
cross2(a, b) = cross([a; 0], [b; 0])[3]
|
||||
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
|
||||
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
|
||||
xi1 = 0.0
|
||||
for i=1:5
|
||||
dxi1 = -R(xi1)/dR(xi1)
|
||||
xi1 += dxi1
|
||||
if norm(dxi1) < 1.0e-10
|
||||
return xi1
|
||||
end
|
||||
end
|
||||
|
||||
error("find projection from master to slave: did not converge")
|
||||
|
||||
end
|
||||
|
||||
function project_from_slave_to_master{E<:MortarElements2D}(
|
||||
master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI)
|
||||
|
||||
function x2(xi2)
|
||||
N = get_basis(E, xi2)
|
||||
return vec(N)*x2_
|
||||
end
|
||||
|
||||
function dx2(xi2)
|
||||
dN = get_dbasis(E, xi2)
|
||||
return vec(dN)*x2_
|
||||
end
|
||||
|
||||
cross2(a, b) = cross([a; 0], [b; 0])[3]
|
||||
R(xi2) = cross2(x2(xi2)-x1, n1)
|
||||
dR(xi2) = cross2(dx2(xi2), n1)
|
||||
xi2 = 0.0
|
||||
dxi2 = 0.0
|
||||
for i=1:5
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < 1.0e-10
|
||||
return xi2
|
||||
end
|
||||
end
|
||||
|
||||
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
|
||||
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(u::Matrix, la::Matrix)
|
||||
|
||||
X1 = slave_element("geometry", time)
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
# 1. update nodal normals for this element
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
fc = SparseMatrixCOO{Real}([], [], []) # interface virtual work
|
||||
C = SparseMatrixCOO{Real}([], [], []) # constraints
|
||||
B = SparseMatrixCOO{Real}([], [], [])
|
||||
|
||||
#info("u1.data = ", ForwardDiff.get_value(u1.data))
|
||||
info("normal calculations done. looping master elements.")
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
info("master element ready.")
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
info("calculating segmentation.")
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution
|
||||
info("xi1 = $xi1")
|
||||
|
||||
info("create bi-orthogonal basis")
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
info("bi-orthogonal basis done. integrating fc.")
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
info("integrate fc")
|
||||
D = zeros(nnodes, nnodes)
|
||||
M = zeros(nnodes, nnodes)
|
||||
gn = zeros(nnodes)
|
||||
lan = zeros(nnodes)
|
||||
lat = zeros(nnodes)
|
||||
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
# project gauss point to master element to evaluate shape function there
|
||||
x_s = vec(N1)*x1 # coordinate in gauss point
|
||||
n_s = vec(N1)*n1 # normal direction in gauss point
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = get_basis(master_element, xi_m)
|
||||
x_m = vec(N2)*x2
|
||||
Phi = vec(Ae*N1')
|
||||
la_s = Phi*la1 # traction force in gauss point
|
||||
u_s = vec(N1)*u1
|
||||
u_m = vec(N2)*u2
|
||||
#info("la_s = $(ForwardDiff.get_value(la_s))")
|
||||
SM = [slave_dofs; master_dofs]
|
||||
#N1N2 = [N1 -N2]
|
||||
#info("all_dofs = $(SM)")
|
||||
#info("shape functions = $(ForwardDiff.get_value(N1N2))")
|
||||
#for i=1:field_dim
|
||||
# #info("add to slave dofs $(slave_dofs[i:field_dim:end])")
|
||||
# #info("add to master dofs $(master_dofs[i:field_dim:end])")
|
||||
# add!(fc, slave_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*N1)
|
||||
# add!(fc, master_dofs[i:field_dim:end], [1, 1], +w*la_s[i]*N2)
|
||||
#add!(fc, slave_dofs[i:field_dim:end], [1, 1], w*la_s'*u_s[i])
|
||||
#add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s'*u_m[i])
|
||||
# add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*u_m)
|
||||
#end
|
||||
#add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1, 1, 1, 1, 1], w*la_s*[u_s' -u_m'])
|
||||
D += w*kron(Ae*N1', N1)
|
||||
M += w*kron(Ae*N1', N2)
|
||||
gn += -w*dot(n_s, x_s-x_m)*Phi
|
||||
lan += w*dot(n_s, la_s)*Phi
|
||||
t_s = Q'*n_s
|
||||
lat += w*dot(t_s, la_s)*Phi
|
||||
end
|
||||
|
||||
#D2 = zeros(2*nnodes, 2*nnodes)
|
||||
#M2 = zeros(2*nnodes, 2*nnodes)
|
||||
#for i=1:field_dim
|
||||
# D2[i:field_dim:end, i:field_dim:end] += D
|
||||
# M2[i:field_dim:end, i:field_dim:end] += M
|
||||
#end
|
||||
#info("size of D2 = $(size(D2))")
|
||||
#fco = [D2 -M2]*vec(la1)
|
||||
#fco = [D -M]*la[:,slave_element_nodes]
|
||||
#info("fco = $(ForwardDiff.get_value(fco))")
|
||||
#add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1], fco)
|
||||
|
||||
for i=1:field_dim
|
||||
add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D)
|
||||
add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M)
|
||||
end
|
||||
info("gn = $gn")
|
||||
#Cj = lan - max(0, lan - gn) + lat
|
||||
add!(C, slave_dofs[1:field_dim:end], [1, 1], gn')
|
||||
|
||||
end # master elements done
|
||||
|
||||
ndofs = prod(size(la))
|
||||
N = SparseMatrixCOO{Real}([], [], [])
|
||||
T = SparseMatrixCOO{Real}([], [], [])
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
add!(N, [dofs[1]], dofs, reshape(n1[i], 1, 2))
|
||||
add!(T, [dofs[2]], dofs, reshape(Q'*n1[i], 1, 2))
|
||||
end
|
||||
N = sparse(N, ndofs, ndofs)
|
||||
T = sparse(T, ndofs, ndofs)
|
||||
B = sparse(B, ndofs, ndofs)
|
||||
fc = B'*vec(la)
|
||||
#println(sparse(fc))
|
||||
#fc = sparse(fc, ndofs, 1)
|
||||
#println(fc)
|
||||
#dump(full(fc))
|
||||
#C = sparse(C, ndofs, 1)
|
||||
C = N*B*vec(u) + T*vec(la)
|
||||
return fc, C
|
||||
|
||||
end
|
||||
|
||||
|
||||
function calculate_interface_PE(x::Vector)
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
#fixed_la = ForwardDiff.get_value(la)
|
||||
#fixed_u = ForwardDiff.get_value(u)
|
||||
#u = ForwardDiff.get_value(u)
|
||||
|
||||
X1 = slave_element("geometry", time)
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
#fixed_la1 = Field(Vector[fixed_la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
# 1. update nodal normals for this element
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
Wco = 0.0
|
||||
Wla = 0.0
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution
|
||||
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s))
|
||||
|
||||
# project gauss point to master element to evaluate shape function there
|
||||
x_s = N1*x1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = Q'*n_s
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = vec(get_basis(master_element, xi_m))
|
||||
x_m = N2*x2
|
||||
Phi = Ae*N1
|
||||
|
||||
gn = -dot(n_s, x_s - x_m)
|
||||
gt = dot(t_s, x_s - x_m)
|
||||
lan = dot(n_s, Phi*la1)
|
||||
lat = dot(t_s, Phi*la1)
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
gu = dot(n_s, u_s - u_m)
|
||||
|
||||
Wco += w*dot(Phi*la1, N1*u1 - N2*u2)
|
||||
#Wla += w*(lan*gn + lat*gt)
|
||||
#gn = min(0, gn)
|
||||
Wla += 1/2*w*1e6*gn*gn
|
||||
#info("gn = $(ForwardDiff.get_value(gn))")
|
||||
#Wla += w*dot(dot(n_s, Phi*la1), dot(n_s, N1*u1 - N2*u2))
|
||||
end
|
||||
|
||||
end
|
||||
return Wco, Wla
|
||||
end
|
||||
|
||||
|
||||
function calculate_contact_rhs(x::Vector)
|
||||
ndofs = round(Int, length(x)/2)
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc, C = calculate_interface(u, la)
|
||||
info("interface vector calculated.")
|
||||
return vec(full([fc; C]))
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
if ndofs == 0
|
||||
info("INITIALIZING THINGS")
|
||||
problem.assembly.u = zeros(16)
|
||||
problem.assembly.la = zeros(16)
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
end
|
||||
|
||||
function add_fco!()
|
||||
get_PI(x::Vector) = calculate_interface_PE(x)[1]
|
||||
A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.gradient(allresults)
|
||||
info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(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]
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
end
|
||||
#add_fco!()
|
||||
|
||||
function add_wla!()
|
||||
get_PI(x::Vector) = calculate_interface_PE(x)[2]
|
||||
A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.gradient(allresults)
|
||||
info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
#info("A")
|
||||
#println(full(A))
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(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]
|
||||
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
end
|
||||
add_wla!()
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
mesh, body1, body2, bc_top, bc_bottom, contact = divided_block_problem()
|
||||
bc_top.properties.formulation = :incremental
|
||||
bc_bottom.properties.formulation = :incremental
|
||||
contact.properties.formulation = :forwarddiff
|
||||
contact.assembly.u = zeros(16)
|
||||
contact.assembly.la = zeros(16)
|
||||
assemble!(contact.assembly, contact, contact.elements[1], 0.0, Val{:forwarddiff})
|
||||
|
||||
Reference in New Issue
Block a user