mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
updated contact codes
This commit is contained in:
Jukka Aho
@@ -109,6 +109,7 @@ include("problems_contact.jl")
|
||||
include("problems_contact_2d.jl")
|
||||
include("problems_contact_3d.jl")
|
||||
include("problems_contact_2d_autodiff.jl")
|
||||
#include("problems_contact_3d_autodiff.jl")
|
||||
export Contact
|
||||
|
||||
module Preprocess
|
||||
|
||||
+1
-1
@@ -603,7 +603,7 @@ function process_output_request(model::Model, solver::Solver, output_request::Ab
|
||||
results = join(results, tables[i], on=:ELEMENT, kind=:outer)
|
||||
end
|
||||
end
|
||||
sort!(results, cols=[:ELEMENT, :IP])
|
||||
#sort!(results; cols=[:ELEMENT, :IP])
|
||||
# filter out elements with id -1, they are automatically created boundary elements
|
||||
fel = find(results[:ELEMENT] .!= Symbol("E-1"))
|
||||
results = results[fel, :]
|
||||
|
||||
+42
-17
@@ -11,11 +11,11 @@ type Element{E<:AbstractElement}
|
||||
properties :: E
|
||||
end
|
||||
|
||||
function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity::Vector{Int64})
|
||||
function Element{E<:AbstractElement}(::Type{E}, id::Int64, connectivity=[])
|
||||
return Element{E}(id, connectivity, [], Dict(), E())
|
||||
end
|
||||
|
||||
function Element{E<:AbstractElement}(::Type{E}, connectivity::Vector{Int64})
|
||||
function Element{E<:AbstractElement}(::Type{E}, connectivity=[])
|
||||
return Element{E}(-1, connectivity, [], Dict(), E())
|
||||
end
|
||||
|
||||
@@ -44,11 +44,11 @@ function setindex!(element::Element, data, field_name)
|
||||
element.fields[field_name] = Field(data)
|
||||
end
|
||||
|
||||
function call(element::Element, field_name)
|
||||
function call(element::Element, field_name::AbstractString)
|
||||
return element[field_name]
|
||||
end
|
||||
|
||||
function call(element::Element, field_name, time)
|
||||
function call(element::Element, field_name::AbstractString, time::Float64)
|
||||
return element[field_name](time)
|
||||
end
|
||||
|
||||
@@ -60,6 +60,31 @@ function call(element::Element, ip, time::Float64=0.0)
|
||||
return get_basis(element, ip, time)
|
||||
end
|
||||
|
||||
"""
|
||||
Examples
|
||||
|
||||
julia> el = Element(Quad4, [1, 2, 3, 4]);
|
||||
|
||||
julia> el([0.0, 0.0], 0.0, 1)
|
||||
1x4 Array{Float64,2}:
|
||||
0.25 0.25 0.25 0.25
|
||||
|
||||
julia> el([0.0, 0.0], 0.0, 2)
|
||||
2x8 Array{Float64,2}:
|
||||
0.25 0.0 0.25 0.0 0.25 0.0 0.25 0.0
|
||||
0.0 0.25 0.0 0.25 0.0 0.25 0.0 0.25
|
||||
|
||||
"""
|
||||
function call(element::Element, ip, time::Float64, dim::Int)
|
||||
dim == 1 && return get_basis(element, ip, time)
|
||||
Ni = get_basis(element, ip, time)
|
||||
N = zeros(dim, length(element)*dim)
|
||||
for i=1:dim
|
||||
N[i,i:dim:end] += Ni
|
||||
end
|
||||
return N
|
||||
end
|
||||
|
||||
function call(element::Element, ip, time::Float64, ::Type{Val{:Jacobian}})
|
||||
X = element("geometry", time)
|
||||
dN = get_dbasis(element, ip, time)
|
||||
@@ -96,15 +121,15 @@ function call(element::Element, field_name::AbstractString, ip, time::Float64, :
|
||||
return element(ip, time, Val{:Grad})*element[field_name](time)
|
||||
end
|
||||
|
||||
function call(element::Element, field::Field, time)
|
||||
function call(element::Element, field::Field, time::Float64)
|
||||
return field(time)
|
||||
end
|
||||
|
||||
function call(element::Element, field::DCTI, time)
|
||||
function call(element::Element, field::DCTI, time::Float64)
|
||||
return field.data
|
||||
end
|
||||
|
||||
function call(element::Element, field_name::AbstractString, time)
|
||||
function call(element::Element, field_name::AbstractString, time::Float64)
|
||||
field = element[field_name]
|
||||
return element(field, time)
|
||||
end
|
||||
@@ -157,7 +182,7 @@ julia> update!(element, "geometry", data)
|
||||
As a result element now have time invariant (variable) vector field "geometry" with data ([0.0, 0.0], [1.0, 2.0]).
|
||||
|
||||
"""
|
||||
function update!(element::Element, field_name, data::Dict)
|
||||
function update!(element::Element, field_name::AbstractString, data::Dict)
|
||||
#element[field_name] = Field(data)
|
||||
element[field_name] = [data[i] for i in get_connectivity(element)]
|
||||
end
|
||||
@@ -195,7 +220,7 @@ function update!(element::Element, field_name, data::Pair...)
|
||||
end
|
||||
=#
|
||||
|
||||
function update!(element::Element, field_name, data::Pair{Float64, Vector{Any}})
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Any}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
@@ -203,7 +228,7 @@ function update!(element::Element, field_name, data::Pair{Float64, Vector{Any}})
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, data::Pair{Float64, Vector{Int64}})
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Int64}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
@@ -211,7 +236,7 @@ function update!(element::Element, field_name, data::Pair{Float64, Vector{Int64}
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, data::Pair{Float64, Vector{Float64}})
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Float64}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
@@ -219,7 +244,7 @@ function update!(element::Element, field_name, data::Pair{Float64, Vector{Float6
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, data::Pair{Float64, Vector{Vector{Float64}}})
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Vector{Vector{Float64}}})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
@@ -227,7 +252,7 @@ function update!(element::Element, field_name, data::Pair{Float64, Vector{Vector
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, data::Pair{Float64, Float64})
|
||||
function update!(element::Element, field_name::AbstractString, data::Pair{Float64, Float64})
|
||||
if haskey(element, field_name)
|
||||
update!(element[field_name], data)
|
||||
else
|
||||
@@ -257,22 +282,22 @@ function update!(element::Element, datas::Pair...)
|
||||
end
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, data::Function)
|
||||
function update!(element::Element, field_name::AbstractString, data::Function)
|
||||
element[field_name] = data
|
||||
end
|
||||
|
||||
function update!(element::Element, field_name, field::Field)
|
||||
function update!(element::Element, field_name::AbstractString, field::Field)
|
||||
element[field_name] = field
|
||||
end
|
||||
|
||||
function update!(elements::Vector, field_name, data)
|
||||
function update!(elements::Vector, field_name::AbstractString, data)
|
||||
for element in elements
|
||||
update!(element, field_name, data)
|
||||
end
|
||||
end
|
||||
|
||||
""" Check existence of field. """
|
||||
function haskey(element::Element, field_name)
|
||||
function haskey(element::Element, field_name::AbstractString)
|
||||
haskey(element.fields, field_name)
|
||||
end
|
||||
|
||||
|
||||
@@ -677,3 +677,6 @@ function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet
|
||||
return all(xi .>= 0.0) && (sum(xi) <= 1.0)
|
||||
end
|
||||
|
||||
function get_reference_coordinates{E}(element::Element{E})
|
||||
get_reference_coordinates(E)
|
||||
end
|
||||
|
||||
+2
-2
@@ -166,13 +166,13 @@ function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{
|
||||
end
|
||||
|
||||
function JuliaFEM.Problem{P<:FieldProblem}(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int64)
|
||||
problem = Problem{P}(name, dimension, "none", [], Dict(), Assembly(), P())
|
||||
problem = Problem(P, name, dimension)
|
||||
problem.elements = create_elements(mesh, name)
|
||||
return problem
|
||||
end
|
||||
|
||||
function JuliaFEM.Problem{P<:BoundaryProblem}(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name)
|
||||
problem = Problem{P}(name, dimension, parent_field_name, [], Dict(), Assembly(), P())
|
||||
problem = Problem(P, name, dimension, parent_field_name)
|
||||
problem.elements = create_elements(mesh, name)
|
||||
return problem
|
||||
end
|
||||
|
||||
+10
-10
@@ -284,11 +284,11 @@ function length(problem::Problem)
|
||||
end
|
||||
|
||||
function update!(problem::Problem, field_name::AbstractString, data)
|
||||
if haskey(problem.fields, field_name)
|
||||
update!(problem.fields[field_name], field_name::AbstractString, data)
|
||||
else
|
||||
problem.fields[field_name] = Field(data)
|
||||
end
|
||||
#if haskey(problem.fields, field_name)
|
||||
# update!(problem.fields[field_name], field_name::AbstractString, data)
|
||||
#else
|
||||
# problem.fields[field_name] = Field(data)
|
||||
#end
|
||||
update!(problem.elements, field_name::AbstractString, data)
|
||||
end
|
||||
|
||||
@@ -302,9 +302,9 @@ end
|
||||
|
||||
""" Return field calculated to nodal points for elements in problem p. """
|
||||
function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
|
||||
if haskey(problem, field_name)
|
||||
return problem[field_name](time)
|
||||
end
|
||||
#if haskey(problem, field_name)
|
||||
# return problem[field_name](time)
|
||||
#end
|
||||
f = nothing
|
||||
for element in get_elements(problem)
|
||||
haskey(element, field_name) || continue
|
||||
@@ -323,8 +323,8 @@ function call(problem::Problem, field_name::AbstractString, time::Float64=0.0)
|
||||
end
|
||||
end
|
||||
end
|
||||
f == nothing && return f
|
||||
update!(problem, field_name, time => f)
|
||||
#f == nothing && return f
|
||||
#update!(problem, field_name, time => f)
|
||||
return f
|
||||
end
|
||||
|
||||
|
||||
+285
-15
@@ -182,10 +182,10 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
info("contact ndofs: $ndofs")
|
||||
info("Rn = $Rn")
|
||||
# info("contact ndofs: $ndofs")
|
||||
# info("Rn = $Rn")
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C1 = sparse(problem.assembly.C1, ndofs, ndofs)
|
||||
C2 = sparse(problem.assembly.C2, ndofs, ndofs)
|
||||
D = sparse(problem.assembly.D, ndofs, ndofs)
|
||||
g = full(problem.assembly.g, ndofs, 1)
|
||||
@@ -216,11 +216,11 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
_c1 = complementarity_condition[j][1]
|
||||
_c2 = c[dofs]
|
||||
_c3 = contact_pressure[j][1]
|
||||
_c4 = g[dofs]
|
||||
info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4")
|
||||
# _c1 = complementarity_condition[j][1]
|
||||
# _c2 = c[dofs]
|
||||
# _c3 = contact_pressure[j][1]
|
||||
# _c4 = g[dofs]
|
||||
# info("active $j: c1 = $_c1, c2 = $_c2, c3 = $_c3, c4 = $_c4")
|
||||
end
|
||||
end
|
||||
|
||||
@@ -246,12 +246,12 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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
|
||||
# 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
|
||||
|
||||
@@ -273,7 +273,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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")
|
||||
# info("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
@@ -287,3 +287,273 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
problem.assembly.g = g
|
||||
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Frictionless 2d small sliding contact without forwarddiff.
|
||||
|
||||
true/false flags: finite_sliding, friction, use_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)
|
||||
|
||||
Rn = 0.0
|
||||
|
||||
# 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)
|
||||
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
|
||||
|
||||
if norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) > props.distval
|
||||
continue
|
||||
end
|
||||
|
||||
if norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) > props.distval
|
||||
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)
|
||||
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
|
||||
t_s = N1*t1 # tangent condition in gauss point
|
||||
n_s /= norm(n_s)
|
||||
t_s /= norm(t_s)
|
||||
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
|
||||
|
||||
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')
|
||||
|
||||
# virtual work
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
contact_area += w
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end
|
||||
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
|
||||
# add contribution to contact virtual work
|
||||
D2 = zeros(field_dim*nsl, field_dim*nsl)
|
||||
M2 = zeros(field_dim*nsl, field_dim*nsl)
|
||||
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)
|
||||
ge = -D2*vec(x1)+M2*vec(x2)
|
||||
add!(problem.assembly.g, sdofs, Q2'*ge)
|
||||
ce = vec(la1) + ge
|
||||
add!(problem.assembly.c, sdofs, Q2'*ce)
|
||||
|
||||
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}()
|
||||
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
|
||||
C1 = sparse(problem.assembly.C1)
|
||||
C2 = sparse(problem.assembly.C2, ndofs, ndofs)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
c = full(problem.assembly.c, ndofs, 1)
|
||||
g = full(problem.assembly.g, ndofs, 1)
|
||||
|
||||
# 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]
|
||||
complementarity_condition[j] = c[dofs]
|
||||
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
|
||||
|
||||
# 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
|
||||
|
||||
end
|
||||
|
||||
|
||||
@@ -121,12 +121,21 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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, 3)
|
||||
dN = get_dbasis(element, ip, time)
|
||||
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
|
||||
=#
|
||||
dN = get_dbasis(element, [0.0], time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
n = Q*t'
|
||||
n /= norm(n)
|
||||
for c in conn
|
||||
normals[:,c] += n
|
||||
end
|
||||
end
|
||||
for i in 1:size(normals,2)
|
||||
normals[:,i] /= norm(normals[:,i])
|
||||
|
||||
+47
-25
@@ -12,7 +12,7 @@ function create_orthogonal_basis(n)
|
||||
end
|
||||
|
||||
"""
|
||||
Frictionless 2d small sliding contact.
|
||||
Frictionless 3d small sliding contact.
|
||||
|
||||
problem
|
||||
time
|
||||
@@ -33,7 +33,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
|
||||
normals = calculate_normals(slave_elements, time, Val{2};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", normals)
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for (slave_num, slave_element) in enumerate(slave_elements)
|
||||
@@ -43,14 +43,29 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
u1 = slave_element("displacement", time)
|
||||
la = slave_element("reaction force", time)
|
||||
n1 = slave_element("normal", time)
|
||||
t11, t21 = create_orthogonal_basis(n1[1])
|
||||
t12, t22 = create_orthogonal_basis(n1[2])
|
||||
t13, t23 = create_orthogonal_basis(n1[3])
|
||||
Q1_ = [n1[1] t11 t21]
|
||||
Q2_ = [n1[2] t12 t22]
|
||||
Q3_ = [n1[3] t13 t23]
|
||||
Z = zeros(3, 3)
|
||||
Q3 = [Q1_ Z Z; Z Q2_ Z; Z Z Q3_]
|
||||
if nsl == 3
|
||||
t11, t21 = create_orthogonal_basis(n1[1])
|
||||
t12, t22 = create_orthogonal_basis(n1[2])
|
||||
t13, t23 = create_orthogonal_basis(n1[3])
|
||||
Q1_ = [n1[1] t11 t21]
|
||||
Q2_ = [n1[2] t12 t22]
|
||||
Q3_ = [n1[3] t13 t23]
|
||||
Z = zeros(3, 3)
|
||||
Q3 = [Q1_ Z Z; Z Q2_ Z; Z Z Q3_]
|
||||
elseif nsl == 4
|
||||
t11, t21 = create_orthogonal_basis(n1[1])
|
||||
t12, t22 = create_orthogonal_basis(n1[2])
|
||||
t13, t23 = create_orthogonal_basis(n1[3])
|
||||
t14, t24 = create_orthogonal_basis(n1[4])
|
||||
Q1_ = [n1[1] t11 t21]
|
||||
Q2_ = [n1[2] t12 t22]
|
||||
Q3_ = [n1[3] t13 t23]
|
||||
Q4_ = [n1[4] t14 t24]
|
||||
Z = zeros(3, 3)
|
||||
Q3 = [Q1_ Z Z Z; Z Q2_ Z Z; Z Z Q3_ Z; Z Z Z Q4_]
|
||||
else
|
||||
error("nsl = $nsl")
|
||||
end
|
||||
contact_area = 0.0
|
||||
contact_error = 0.0
|
||||
|
||||
@@ -64,15 +79,19 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
update!(slave_element, "element area", time => element_area)
|
||||
end
|
||||
|
||||
if slave_num == 1
|
||||
info("First slave element area = $element_area")
|
||||
info("NT basis of first slave element")
|
||||
dump(Q3)
|
||||
end
|
||||
# if slave_num == 1
|
||||
# info("First slave element area = $element_area")
|
||||
# info("NT basis of first slave element")
|
||||
# dump(Q3)
|
||||
# end
|
||||
|
||||
# project slave nodes to auxiliary plane (x0, Q)
|
||||
#xi = get_reference_element_midpoint(slave_element)
|
||||
xi = [1/3, 1/3]
|
||||
if nsl == 3
|
||||
xi = [1/3, 1/3]
|
||||
else
|
||||
xi = [1/4, 1/4]
|
||||
end
|
||||
N = vec(get_basis(slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
@@ -151,7 +170,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D3 = zeros(nsldofs, nsldofs)
|
||||
M3 = zeros(nmdofs, nmdofs)
|
||||
M3 = zeros(nsldofs, nmdofs)
|
||||
for i=1:field_dim
|
||||
D3[i:field_dim:end, i:field_dim:end] += De
|
||||
M3[i:field_dim:end, i:field_dim:end] += Me
|
||||
@@ -180,8 +199,14 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
is_slip = Dict{Int64, Int}()
|
||||
is_stick = Dict{Int64, Int}()
|
||||
|
||||
g = full(problem.assembly.g)
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
|
||||
C1 = sparse(problem.assembly.C1, ndofs, ndofs)
|
||||
C2 = sparse(problem.assembly.C2, ndofs, ndofs)
|
||||
D = sparse(problem.assembly.D, ndofs, ndofs)
|
||||
g = full(problem.assembly.g, ndofs, 1)
|
||||
c = full(problem.assembly.c, ndofs, 1)
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
@@ -233,26 +258,23 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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 = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
tdofs = dofs[[2,3]]
|
||||
if (is_active[j] == 1) && (is_slip[j] == 1)
|
||||
info("$j is in active/slip, removing tangential constraints $tdofs")
|
||||
# info("$j is in active/slip, removing tangential constraints $tdofs")
|
||||
C2[tdofs,:] = 0.0
|
||||
g[tdofs] = 0.0
|
||||
normal = normals[j]
|
||||
@@ -266,7 +288,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
for j in S
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
if is_inactive[j] == 1
|
||||
info("$j is inactive, removing dofs $dofs")
|
||||
# info("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
|
||||
@@ -105,7 +105,7 @@ function get_polygon_clip(xs, xm, n; debug=false)
|
||||
# objective: search does line xm1 - xm2 clip xs
|
||||
nm = length(xm)
|
||||
ns = length(xs)
|
||||
P = Vector{Float64}[]
|
||||
P = Vector{Number}[]
|
||||
|
||||
# 1. test is master point inside slave, if yes, add to clip
|
||||
for i=1:nm
|
||||
|
||||
@@ -75,9 +75,11 @@ using JuliaFEM.Testing
|
||||
|
||||
# element details
|
||||
el = first(block.elements)
|
||||
X = el("geometry", 0.0)
|
||||
debug("X = $X")
|
||||
S1 = block(el, [0.0, 0.0], 0.0, Val{:S})
|
||||
S1 = S1[[1,4,2]]
|
||||
E1= block(el, [0.0, 0.0], 0.0, Val{:E})
|
||||
E1 = block(el, [0.0, 0.0], 0.0, Val{:E})
|
||||
E1 = E1[[1,4,2]]
|
||||
C1 = block(el, [0.0, 0.0], 0.0, Val{:COORD})
|
||||
info("strain = $E1, stress = $S1, at $C1")
|
||||
@@ -103,12 +105,12 @@ using JuliaFEM.Testing
|
||||
info("u3 = $u")
|
||||
@test isapprox(u, u3_expected)
|
||||
|
||||
info("calling nonlinear solver")
|
||||
solver2 = NonlinearSolver(block, traction, bc_sym_23, bc_sym_13)
|
||||
solver2()
|
||||
u = solver2("displacement", 0.0)[3]
|
||||
info("nlsolver u3 = $u, expected = $u3_expected")
|
||||
@test isapprox(u, u3_expected; rtol=1.0e-5)
|
||||
# info("calling nonlinear solver")
|
||||
# solver2 = NonlinearSolver(block, traction, bc_sym_23, bc_sym_13)
|
||||
# solver2()
|
||||
# u = solver2("displacement", 0.0)[3]
|
||||
# info("nlsolver u3 = $u, expected = $u3_expected")
|
||||
# @test isapprox(u, u3_expected; rtol=1.0e-5)
|
||||
end
|
||||
|
||||
#= TODO: to other file
|
||||
|
||||
Reference in New Issue
Block a user