mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-10-03 14:47:55 +00:00
Contact algorithms testing & develpoment (#95)
* contact 3d patch test, standard lagrange, small sliding, linear tet4 elements * tet4 dual basis contact patch test pass * contact 3d patch test, standard lagrange, small sliding, linear tet4 elements * tet4 dual basis contact patch test pass * Patch test for linear elements standard lagrange / dual lagrange pass now * Patch test for quadratic contact surfaces for standard + dual basis pass * refactoring * renamed files * Improvements to preprocess scripts * convert several elements to node sets in one command * possibility to find particular node from mesh filtered by node set * 2d small sliding contact patch test, linear elements * Added backward compatibility * 2d contact algorithms pass patch tests * test data for 2d contacts * no common models in different tests. testing generalized alpha stabilization * Preprocess tests * moved tests from test_preprocess_aster_reader.jl to test_preprocess.jl * generalized-alpha time integration, alpha=0.0 by default * Improvements to logging * JuliaFEM.jl: can set environment variable to one of logging levels: OFF, CRITICAL, ERROR, WARNING, INFO, DEBUG * problems_contact_2d_autodiff.jl: do not loop over nodes if logging level != DEBUG
This commit is contained in:
Jukka Aho
committed by
Tero Frondelius
Tero Frondelius
parent
28e09f0305
commit
a0d18568ce
+1
-12
@@ -17,8 +17,7 @@ using Logging
|
||||
Logging.configure(level=INFO)
|
||||
|
||||
if haskey(ENV, "JULIAFEM_LOGLEVEL")
|
||||
ENV["JULIAFEM_LOGLEVEL"] == "INFO" && Logging.configure(level=INFO)
|
||||
ENV["JULIAFEM_LOGLEVEL"] == "DEBUG" && Logging.configure(level=DEBUG)
|
||||
Logging.configure(level=LogLevel(ENV["JULIAFEM_LOGLEVEL"]))
|
||||
end
|
||||
|
||||
export info, debug
|
||||
@@ -141,16 +140,6 @@ export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
|
||||
aster_read_mesh_names, aster_read_node_sets, aster_read_nodes, RMEDFile
|
||||
end
|
||||
|
||||
function get_mesh(mesh_name::AbstractString, args...; kwargs...)
|
||||
return get_mesh(Val{Symbol(mesh_name)}, args...; kwargs...)
|
||||
end
|
||||
|
||||
function get_model(model_name::AbstractString, args...; kwargs...)
|
||||
return get_model(Val{Symbol(model_name)}, args...; kwargs...)
|
||||
end
|
||||
|
||||
export get_mesh, get_model
|
||||
|
||||
module Postprocess
|
||||
|
||||
include("postprocess_utils.jl")
|
||||
|
||||
+15
-2
@@ -25,6 +25,10 @@ type Contact <: BoundaryProblem
|
||||
remove_nodes :: Vector{Int}
|
||||
always_in_contact :: Bool
|
||||
update_contact_pairing :: Bool
|
||||
iteration :: Int
|
||||
contact_state_in_first_iteration :: Symbol
|
||||
alpha :: Float64
|
||||
drop_tolerance :: Float64
|
||||
store_fields :: Vector{AbstractString}
|
||||
end
|
||||
|
||||
@@ -47,6 +51,10 @@ function Contact()
|
||||
[], # remove these nodes always from set
|
||||
false, # mainly for debugging, do not remove inactive nodes
|
||||
true, # update contact pairing on each loop
|
||||
1, # iteration counter
|
||||
:AUTO, # contact state in first iteration, AUTO, INACTIVE, ACTIVE
|
||||
0.0, # alpha basis transform parameter
|
||||
1.0e-9, # drop tolerance
|
||||
default_fields)
|
||||
end
|
||||
|
||||
@@ -55,24 +63,29 @@ function get_unknown_field_name(problem::Problem{Contact})
|
||||
end
|
||||
|
||||
function get_formulation_type(problem::Problem{Contact})
|
||||
#=
|
||||
if problem.properties.use_forwarddiff
|
||||
return :forwarddiff
|
||||
else
|
||||
return :incremental
|
||||
end
|
||||
=#
|
||||
return :incremental
|
||||
#return :forwarddiff
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Contact}, time::Real)
|
||||
if problem.properties.dimension == -1
|
||||
problem.properties.dimension = dim = size(first(problem.elements), 1)
|
||||
info("assuming dimension of mesh tie surface is $dim")
|
||||
info("if this is wrong set is manually using problem.properties.dimension")
|
||||
debug("assuming dimension of mesh tie surface is $dim")
|
||||
debug("if this is wrong set is manually using problem.properties.dimension")
|
||||
end
|
||||
dimension = Val{problem.properties.dimension}
|
||||
finite_sliding = Val{problem.properties.finite_sliding}
|
||||
friction = Val{problem.properties.friction}
|
||||
use_forwarddiff = Val{problem.properties.use_forwarddiff}
|
||||
assemble!(problem, time, dimension, finite_sliding, friction, use_forwarddiff)
|
||||
problem.properties.iteration += 1
|
||||
end
|
||||
|
||||
typealias ContactElements2D Union{Seg2}
|
||||
|
||||
+123
-66
@@ -1,14 +1,62 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
function create_rotation_matrix(element::Element{Seg2}, time::Float64)
|
||||
n = element("normal", time)
|
||||
R = [0.0 -1.0; 1.0 0.0]
|
||||
t1 = R'*n[1]
|
||||
t2 = R'*n[2]
|
||||
Q1 = [n[1] t1]
|
||||
Q2 = [n[2] t2]
|
||||
Z = zeros(2, 2)
|
||||
Q = [Q1 Z; Z Q2]
|
||||
return Q
|
||||
end
|
||||
|
||||
function create_contact_segmentation(problem::Problem{Contact}, slave_element::Element{Seg2}, master_elements::Vector, time::Float64; deformed=false)
|
||||
|
||||
result = []
|
||||
|
||||
x1 = slave_element("geometry", time)
|
||||
|
||||
if deformed
|
||||
x1 += slave_element("displacement", time)
|
||||
end
|
||||
|
||||
for master_element in master_elements
|
||||
|
||||
x2 = master_element("geometry", time)
|
||||
|
||||
if deformed
|
||||
x2 += master_element("displacement", time)
|
||||
end
|
||||
|
||||
if norm(mean(x1) - x2[1]) / norm(x1[2] - x1[1]) > problem.properties.distval
|
||||
continue
|
||||
end
|
||||
if norm(mean(x1) - x2[2]) / norm(x1[2] - x1[1]) > problem.properties.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])
|
||||
if isapprox(l, 0.0)
|
||||
continue # no contribution in this master element
|
||||
end
|
||||
push!(result, (master_element, xi1, l))
|
||||
end
|
||||
return result
|
||||
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)
|
||||
function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{1}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
@@ -33,51 +81,22 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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
|
||||
Q2 = create_rotation_matrix(slave_element, time)
|
||||
|
||||
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)
|
||||
master_elements = slave_element("master elements", time)
|
||||
segmentation = create_contact_segmentation(problem, slave_element, master_elements, time)
|
||||
if length(segmentation) == 0 # no overlapping in master and slave surfaces with this slave element
|
||||
continue
|
||||
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
|
||||
Ae = eye(nsl)
|
||||
if props.dual_basis
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
Ae = zeros(nsl, nsl)
|
||||
if props.dual_basis
|
||||
for (master_element, xi1, l) in segmentation
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
@@ -88,14 +107,21 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
Me += w*N1*N1'
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
else
|
||||
Ae = eye(nsl)
|
||||
end
|
||||
end
|
||||
|
||||
# loop all segments
|
||||
for (master_element, xi1, l) in segmentation
|
||||
|
||||
nm = length(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
# 3.3. loop integration points of one integration segment and calculate
|
||||
# local mortar matrices
|
||||
fill!(De, 0.0)
|
||||
fill!(Me, 0.0)
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
Ne = zeros(nsl, 2*nsl)
|
||||
Te = zeros(nsl, 2*nsl)
|
||||
He = zeros(nsl, 2*nsl)
|
||||
@@ -117,7 +143,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
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
|
||||
X_m = N2*X2
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
@@ -182,14 +208,34 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
|
||||
la = problem.assembly.la
|
||||
ndofs = length(la)
|
||||
# info("contact ndofs: $ndofs")
|
||||
# info("Rn = $Rn")
|
||||
|
||||
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)
|
||||
|
||||
for j in S
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
weighted_gap[j] = g[dofs]
|
||||
end
|
||||
|
||||
state = problem.properties.contact_state_in_first_iteration
|
||||
if problem.properties.iteration == 1
|
||||
info("First contact iteration, initial contact state = $state")
|
||||
|
||||
if state == :AUTO
|
||||
avg_gap = mean([weighted_gap[j][1] for j in S])
|
||||
std_gap = std([weighted_gap[j][1] for j in S])
|
||||
if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
|
||||
state = :ACTIVE
|
||||
else
|
||||
state = :UNKNOWN
|
||||
end
|
||||
info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
@@ -204,7 +250,6 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
contact_pressure[j] = [0.0, 0.0]
|
||||
end
|
||||
|
||||
# contact_pressure[j] = c[dofs]
|
||||
complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
if complementarity_condition[j][1] < 0
|
||||
is_inactive[j] = 1
|
||||
@@ -216,11 +261,24 @@ 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")
|
||||
end
|
||||
end
|
||||
|
||||
if (problem.properties.iteration == 1) && (state == :ACTIVE)
|
||||
for j in S
|
||||
is_inactive[j] = 0
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
if (problem.properties.iteration == 1) && (state == :INACTIVE)
|
||||
for j in S
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
@@ -246,34 +304,33 @@ 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
|
||||
debug("# | 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][1], 3)) | $(round(contact_pressure[j][1], 3)) | $(round(complementarity_condition[j][1], 3))"
|
||||
debug(str1 * str2)
|
||||
end
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
debug("normals: ", normals)
|
||||
|
||||
# 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])")
|
||||
debug("$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")
|
||||
debug("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
|
||||
@@ -163,6 +163,38 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
|
||||
nnodes = size(slave_element, 2)
|
||||
|
||||
# construct dual basis
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip, time)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s, time)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
end
|
||||
|
||||
Ae = De*inv(Me)
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
@@ -183,21 +215,6 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip, time)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi = ip.coords[1]
|
||||
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s, time)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
@@ -237,27 +254,67 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
# at this point we have calculated contact force fc and gap for all slave elements.
|
||||
# next task is to find out are they in contact or not and remove inactive nodes
|
||||
|
||||
#nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
|
||||
#info("gap: $nzgap")
|
||||
state = problem.properties.contact_state_in_first_iteration
|
||||
if problem.properties.iteration == 1
|
||||
info("First contact iteration, initial contact state = $state")
|
||||
if state == :AUTO
|
||||
avg_gap = ForwardDiff.value(mean([gap[1, j] for j in S]))
|
||||
std_gap = ForwardDiff.value(std([gap[1, j] for j in S]))
|
||||
if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
|
||||
state = :ACTIVE
|
||||
else
|
||||
state = :UNKNOWN
|
||||
end
|
||||
info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
|
||||
end
|
||||
end
|
||||
|
||||
for (i, j) in enumerate(sort(collect(S)))
|
||||
# if j in props.always_inactive
|
||||
# info("special node $j always inactive")
|
||||
# C[:,j] = la[:,j]
|
||||
# continue
|
||||
# end
|
||||
n = normals[:,j]
|
||||
t = Q'*n
|
||||
lan = dot(n, la[:,j])
|
||||
lat = dot(t, la[:,j])
|
||||
is_active = Dict{Int, Bool}()
|
||||
condition = Dict()
|
||||
|
||||
if lan - gap[1, j] > 0
|
||||
# info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))")
|
||||
for j in S
|
||||
if j in props.always_in_contact
|
||||
is_active[j] = true
|
||||
continue
|
||||
end
|
||||
lan = dot(normals[:,j], la[:,j])
|
||||
condition[j] = ForwardDiff.value(lan - gap[1, j])
|
||||
is_active[j] = condition[j] > 0
|
||||
end
|
||||
|
||||
if problem.properties.iteration == 1 && state == :ACTIVE
|
||||
for j in S
|
||||
is_active[j] = true
|
||||
end
|
||||
end
|
||||
|
||||
if problem.properties.iteration == 1 && state == :INACTIVE
|
||||
for j in S
|
||||
is_active[j] = false
|
||||
end
|
||||
end
|
||||
|
||||
if Logging._root.level == DEBUG
|
||||
debug("Summary of nodes")
|
||||
for j in sort(collect(keys(is_active)))
|
||||
n = map(ForwardDiff.value, normals[:,j])
|
||||
debug("$j, c=$(condition[j]), s=$(is_active[j]), n=$n")
|
||||
end
|
||||
end
|
||||
|
||||
for j in S
|
||||
|
||||
if is_active[j]
|
||||
n = normals[:,j]
|
||||
t = Q'*n
|
||||
lan = dot(n, la[:,j])
|
||||
lat = dot(t, la[:,j])
|
||||
C[1,j] += gap[1, j]
|
||||
C[2,j] += lat
|
||||
else
|
||||
C[:,j] = la[:,j]
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
return vec([fc C])
|
||||
@@ -269,6 +326,7 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
if length(x) == 0
|
||||
error("2d autodiff contact problem: initialize problem.assembly.u & la before solution")
|
||||
end
|
||||
|
||||
A = ForwardDiff.jacobian(calculate_interface, x)
|
||||
b = calculate_interface(x)
|
||||
A = sparse(A)
|
||||
@@ -278,19 +336,40 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(A[1:ndofs,ndofs+1:end])
|
||||
C1 = A[1:ndofs,ndofs+1:end]
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = -b[1:ndofs]
|
||||
g = -b[ndofs+1:end]
|
||||
|
||||
f += C1*problem.assembly.la
|
||||
g += D*problem.assembly.la
|
||||
|
||||
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)
|
||||
#=
|
||||
if !haskey(problem, "contact force")
|
||||
problem.fields["contact force"] = Field(time => f)
|
||||
else
|
||||
update!(problem.fields["contact force"], time => f)
|
||||
end
|
||||
|
||||
fc = problem.fields["contact force"]
|
||||
|
||||
if length(fc) > 1
|
||||
# kick in generalized alpha rule for time integration
|
||||
alpha = 0.5
|
||||
info("Applying Generalized alpha time integration")
|
||||
K = (1-alpha)*K
|
||||
C1 = (1-alpha)*C1
|
||||
f = alpha*fc[end-1].data
|
||||
end
|
||||
=#
|
||||
|
||||
problem.assembly.K = K
|
||||
problem.assembly.C1 = transpose(C1)
|
||||
problem.assembly.C2 = C2
|
||||
problem.assembly.D = D
|
||||
problem.assembly.f = sparse(f)
|
||||
problem.assembly.g = sparse(g)
|
||||
|
||||
end
|
||||
|
||||
|
||||
+565
-193
@@ -11,6 +11,459 @@ function create_orthogonal_basis(n)
|
||||
return t1, t2
|
||||
end
|
||||
|
||||
""" Create rotation matrix Q for element nodes rotating quantities to nt coordinaet system. """
|
||||
function create_rotation_matrix(element::Element{Tri3}, time::Float64)
|
||||
n = element("normal", time)
|
||||
t11, t21 = create_orthogonal_basis(n[1])
|
||||
t12, t22 = create_orthogonal_basis(n[2])
|
||||
t13, t23 = create_orthogonal_basis(n[3])
|
||||
Q1_ = [n[1] t11 t21]
|
||||
Q2_ = [n[2] t12 t22]
|
||||
Q3_ = [n[3] t13 t23]
|
||||
Z = zeros(3, 3)
|
||||
Q = [
|
||||
Q1_ Z Z
|
||||
Z Q2_ Z
|
||||
Z Z Q3_]
|
||||
return Q
|
||||
end
|
||||
|
||||
function create_rotation_matrix(element::Element{Quad4}, time::Float64)
|
||||
n = element("normal", time)
|
||||
t11, t21 = create_orthogonal_basis(n[1])
|
||||
t12, t22 = create_orthogonal_basis(n[2])
|
||||
t13, t23 = create_orthogonal_basis(n[3])
|
||||
t14, t24 = create_orthogonal_basis(n[4])
|
||||
Q1_ = [n[1] t11 t21]
|
||||
Q2_ = [n[2] t12 t22]
|
||||
Q3_ = [n[3] t13 t23]
|
||||
Q4_ = [n[4] t14 t24]
|
||||
Z = zeros(3, 3)
|
||||
Q = [
|
||||
Q1_ Z Z Z
|
||||
Z Q2_ Z Z
|
||||
Z Z Q3_ Z
|
||||
Z Z Z Q4_]
|
||||
return Q
|
||||
end
|
||||
|
||||
function create_rotation_matrix(element::Element{Tri6}, time::Float64)
|
||||
n = element("normal", time)
|
||||
t11, t21 = create_orthogonal_basis(n[1])
|
||||
t12, t22 = create_orthogonal_basis(n[2])
|
||||
t13, t23 = create_orthogonal_basis(n[3])
|
||||
t14, t24 = create_orthogonal_basis(n[4])
|
||||
t15, t25 = create_orthogonal_basis(n[5])
|
||||
t16, t26 = create_orthogonal_basis(n[6])
|
||||
Q1_ = [n[1] t11 t21]
|
||||
Q2_ = [n[2] t12 t22]
|
||||
Q3_ = [n[3] t13 t23]
|
||||
Q4_ = [n[4] t14 t24]
|
||||
Q5_ = [n[5] t15 t25]
|
||||
Q6_ = [n[6] t16 t26]
|
||||
Z = zeros(3, 3)
|
||||
Q = [
|
||||
Q1_ Z Z Z Z Z
|
||||
Z Q2_ Z Z Z Z
|
||||
Z Z Q3_ Z Z Z
|
||||
Z Z Z Q4_ Z Z
|
||||
Z Z Z Z Q5_ Z
|
||||
Z Z Z Z Z Q6_]
|
||||
return Q
|
||||
end
|
||||
|
||||
""" Create a contact segmentation between one slave element and list of master elements.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
||||
Vector with tuples: (master_element, polygon_clip_vertices, polygon_clip_centroid, polygon_clip_area)
|
||||
"""
|
||||
function create_contact_segmentation(slave_element, master_elements, x0, n0, time::Float64; deformed=false)
|
||||
result = []
|
||||
x1 = slave_element("geometry", time)
|
||||
if deformed
|
||||
x1 += slave_element("displacement", time)
|
||||
end
|
||||
S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x1]
|
||||
for master_element in master_elements
|
||||
x2 = master_element("geometry", time)
|
||||
if deformed
|
||||
x2 += master_element("displacement", time)
|
||||
end
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in x2]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
length(P) < 3 && continue # no clipping or shared edge (no volume)
|
||||
check_orientation!(P, n0)
|
||||
N_P = length(P)
|
||||
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
|
||||
if isapprox(P_area, 0.0)
|
||||
error("Polygon P has zero area")
|
||||
end
|
||||
C0 = calculate_centroid(P)
|
||||
push!(result, (master_element, P, C0, P_area))
|
||||
end
|
||||
return result
|
||||
end
|
||||
|
||||
"Assemble linear surface element to contact problem. """
|
||||
function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time::Float64)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = slave_element("displacement", time)
|
||||
x1 = X1 + u1
|
||||
n1 = slave_element("normal", time)
|
||||
la = slave_element("reaction force", time)
|
||||
|
||||
Q3 = create_rotation_matrix(slave_element, time)
|
||||
|
||||
# project slave nodes to auxiliary plane (x0, Q)
|
||||
xi = mean(get_reference_coordinates(slave_element))
|
||||
N = vec(get_basis(slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
|
||||
# create contact segmentation
|
||||
segmentation = create_contact_segmentation(slave_element, slave_element("master elements", time), x0, n0, time)
|
||||
|
||||
if length(segmentation) == 0 # no overlapping surface in slave and maters
|
||||
return
|
||||
end
|
||||
|
||||
Ae = eye(nsl)
|
||||
|
||||
if problem.properties.dual_basis # construct dual basis
|
||||
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
|
||||
# loop all polygons
|
||||
for (master_element, P, C0, P_area) in segmentation
|
||||
|
||||
# loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3, Int[])
|
||||
update!(virtual_element, "geometry", cell)
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
N1 = slave_element(xi_s, time)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
end # master elements done
|
||||
|
||||
Ae = De*inv(Me)
|
||||
|
||||
debug("Dual basis coeffients = $Ae")
|
||||
end
|
||||
|
||||
# loop all polygons
|
||||
for (master_element, P, C0, P_area) in segmentation
|
||||
|
||||
nm = length(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = master_element("displacement", time)
|
||||
x2 = X2 + u2
|
||||
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
ce = zeros(field_dim*nsl)
|
||||
ge = zeros(field_dim*nsl)
|
||||
|
||||
# loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3, Int[])
|
||||
update!(virtual_element, "geometry", cell)
|
||||
# loop integration point of integration cell
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
|
||||
# project gauss point from auxiliary plane to master and slave element
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
|
||||
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
|
||||
# add contributions
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
Phi = Ae*N1
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
x_s = N1*(X1+u1)
|
||||
x_m = N2*(X2+u2)
|
||||
ge += w*vec((x_m-x_s)*Phi')
|
||||
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
# add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D3 = zeros(nsldofs, nsldofs)
|
||||
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
|
||||
end
|
||||
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D3)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M3)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
|
||||
add!(problem.assembly.g, sdofs, Q3'*ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
end
|
||||
|
||||
|
||||
""" Assemble quadratic surface element to contact problem. """
|
||||
function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time::Float64)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
|
||||
alp = props.alpha
|
||||
|
||||
if alp != 0.0
|
||||
T = [
|
||||
1.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 1.0 0.0 0.0 0.0
|
||||
alp alp 0.0 1.0-2*alp 0.0 0.0
|
||||
0.0 alp alp 0.0 1.0-2*alp 0.0
|
||||
alp 0.0 alp 0.0 0.0 1.0-2*alp
|
||||
]
|
||||
else
|
||||
T = eye(6)
|
||||
end
|
||||
|
||||
nsl = length(slave_element)
|
||||
Xs = slave_element("geometry", time)
|
||||
n1 = slave_element("normal", time)
|
||||
|
||||
Q3 = create_rotation_matrix(slave_element, time)
|
||||
|
||||
Ae = eye(nsl)
|
||||
|
||||
if problem.properties.dual_basis # construct dual basis
|
||||
|
||||
nsl = length(slave_element)
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nsl)
|
||||
|
||||
for sub_slave_element in split_quadratic_element(slave_element, time)
|
||||
|
||||
slave_element_nodes = get_connectivity(sub_slave_element)
|
||||
nsl = length(sub_slave_element)
|
||||
|
||||
X1 = sub_slave_element("geometry", time)
|
||||
#u1 = sub_slave_element("displacement", time)
|
||||
#x1 = X1 + u1
|
||||
n1 = sub_slave_element("normal", time)
|
||||
#la = sub_slave_element("reaction force", time)
|
||||
|
||||
# create auxiliary plane
|
||||
xi = mean(get_reference_coordinates(sub_slave_element))
|
||||
N = vec(get_basis(sub_slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
|
||||
# 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", time)
|
||||
|
||||
Xm = master_element("geometry", time)
|
||||
|
||||
if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
|
||||
continue
|
||||
end
|
||||
|
||||
# split master element to linear sub-elements and loop
|
||||
for sub_master_element in split_quadratic_element(master_element, time)
|
||||
|
||||
master_element_nodes = get_connectivity(sub_master_element)
|
||||
nm = length(sub_master_element)
|
||||
X2 = sub_master_element("geometry", time)
|
||||
#u2 = sub_master_element("displacement", time)
|
||||
#x2 = X2 + u2
|
||||
|
||||
# 3.1 project master nodes to auxiliary plane and create polygon clipping
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
length(P) < 3 && continue # no clipping or shared edge (no volume)
|
||||
check_orientation!(P, n0)
|
||||
|
||||
N_P = length(P)
|
||||
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
|
||||
if isapprox(P_area, 0.0)
|
||||
error("Polygon P has zero area")
|
||||
end
|
||||
|
||||
C0 = calculate_centroid(P)
|
||||
|
||||
# 4. loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3, Int[])
|
||||
update!(virtual_element, "geometry", cell)
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
|
||||
N1 = vec(slave_element(xi_s, time)*T)
|
||||
De += w*diagm(N1)
|
||||
Me += w*N1*N1'
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
end # sub master elements done
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # sub slave elements done
|
||||
|
||||
Ae = De*inv(Me)
|
||||
|
||||
debug("Dual basis coeffients = $Ae")
|
||||
end
|
||||
|
||||
# split slave element to linear sub-elements and loop
|
||||
for sub_slave_element in split_quadratic_element(slave_element, time)
|
||||
|
||||
slave_element_nodes = get_connectivity(sub_slave_element)
|
||||
nsl = length(sub_slave_element)
|
||||
X1 = sub_slave_element("geometry", time)
|
||||
n1 = sub_slave_element("normal", time)
|
||||
|
||||
# create auxiliary plane
|
||||
xi = mean(get_reference_coordinates(sub_slave_element))
|
||||
N = vec(get_basis(sub_slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
|
||||
# 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", time)
|
||||
|
||||
Xm = master_element("geometry", time)
|
||||
|
||||
if norm(mean(Xs) - mean(Xm)) > problem.properties.distval
|
||||
continue
|
||||
end
|
||||
|
||||
# split master element to linear sub-elements and loop
|
||||
for sub_master_element in split_quadratic_element(master_element, time)
|
||||
|
||||
master_element_nodes = get_connectivity(sub_master_element)
|
||||
nm = length(master_element)
|
||||
X2 = sub_master_element("geometry", time)
|
||||
#u2 = master_element("displacement", time)
|
||||
#x2 = X2 + u2
|
||||
|
||||
# 3.1 project master nodes to auxiliary plane and create polygon clipping
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
length(P) < 3 && continue # no clipping or shared edge (no volume)
|
||||
check_orientation!(P, n0)
|
||||
|
||||
N_P = length(P)
|
||||
P_area = sum([norm(1/2*cross(P[i]-P[1], P[mod(i,N_P)+1]-P[1])) for i=2:N_P])
|
||||
if isapprox(P_area, 0.0)
|
||||
error("Polygon P has zero area")
|
||||
end
|
||||
|
||||
C0 = calculate_centroid(P)
|
||||
|
||||
# integration is done in quadratic elements
|
||||
nsl = length(slave_element)
|
||||
nm = length(master_element)
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
ge = zeros(field_dim*nsl)
|
||||
|
||||
# 4. loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3, Int[])
|
||||
update!(virtual_element, "geometry", cell)
|
||||
|
||||
# 5. loop integration point of integration cell
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
|
||||
# project gauss point from auxiliary plane to master and slave element
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, Xs, time)
|
||||
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, Xm, time)
|
||||
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
|
||||
# add contributions
|
||||
N1 = vec(get_basis(slave_element, xi_s, time)*T)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
Phi = Ae*N1
|
||||
|
||||
De += w*Phi*N1'
|
||||
Me += w*Phi*N2'
|
||||
|
||||
us = slave_element("displacement", time)
|
||||
um = master_element("displacement", time)
|
||||
xs = N1*(Xs+us)
|
||||
xm = N2*(Xs+um)
|
||||
ge += w*vec((xm-xs)*Phi')
|
||||
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
# 6. add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D3 = zeros(nsldofs, nsldofs)
|
||||
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
|
||||
end
|
||||
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D3)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M3)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
|
||||
add!(problem.assembly.g, sdofs, Q3'*ge)
|
||||
|
||||
end # sub master elements done
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # sub slave elements done
|
||||
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Frictionless 3d small sliding contact.
|
||||
|
||||
@@ -21,9 +474,7 @@ finite_sliding
|
||||
friction
|
||||
use_forwarddiff
|
||||
"""
|
||||
function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
::Type{Val{2}}, ::Type{Val{false}},
|
||||
::Type{Val{false}}, ::Type{Val{false}}; debug=true)
|
||||
function assemble!(problem::Problem{Contact}, time::Float64, ::Type{Val{2}}, ::Type{Val{false}}, ::Type{Val{false}}, ::Type{Val{false}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
@@ -36,153 +487,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
update!(slave_elements, "normal", time => normals)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for (slave_num, slave_element) in enumerate(slave_elements)
|
||||
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = slave_element("displacement", time)
|
||||
la = slave_element("reaction force", time)
|
||||
n1 = slave_element("normal", time)
|
||||
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
|
||||
|
||||
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
|
||||
if "element area" in props.store_fields
|
||||
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
|
||||
|
||||
# project slave nodes to auxiliary plane (x0, Q)
|
||||
xi = mean(get_reference_coordinates(slave_element))
|
||||
N = vec(get_basis(slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
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", 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
|
||||
norm(mean(X1) - X2[3]) / norm(X1[2] - X1[1]) < props.distval || continue
|
||||
=#
|
||||
|
||||
# 3.1 project master nodes to auxiliary plane and create polygon clipping
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
|
||||
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)
|
||||
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
ge = zeros(field_dim*nsl)
|
||||
|
||||
# 4. loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3, Int[])
|
||||
update!(virtual_element, "geometry", cell)
|
||||
|
||||
# 5. loop integration point of integration cell
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
|
||||
# project gauss point from auxiliary plane to master and slave element
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
if isnan(x_gauss[1])
|
||||
info("is nan")
|
||||
info("x_gauss = $x_gauss")
|
||||
info("cell = $cell")
|
||||
info("C0 = $C0")
|
||||
info("P = $P")
|
||||
info("S = $S")
|
||||
info("M = $M")
|
||||
info("n0 = $n0")
|
||||
error("nan, unable to continue")
|
||||
end
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
|
||||
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
# add contributions
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
De += w*N1*N1'
|
||||
Me += w*N1*N2'
|
||||
|
||||
x_s = N1*(X1+u1)
|
||||
x_m = N2*(X2+u2)
|
||||
ge += w*vec((x_m-x_s)*N1')
|
||||
contact_area += w
|
||||
n_s = N1*n1
|
||||
contact_error += 1/2*w*dot(n_s, x_s-x_m)^2
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
# 6. add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
nsldofs = length(sdofs)
|
||||
nmdofs = length(mdofs)
|
||||
D3 = zeros(nsldofs, nsldofs)
|
||||
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
|
||||
end
|
||||
|
||||
add!(problem.assembly.C1, sdofs, sdofs, D3)
|
||||
add!(problem.assembly.C1, sdofs, mdofs, -M3)
|
||||
add!(problem.assembly.C2, sdofs, sdofs, Q3'*D3)
|
||||
add!(problem.assembly.C2, sdofs, mdofs, -Q3'*M3)
|
||||
add!(problem.assembly.g, sdofs, Q3'*ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
if "contact area" in props.store_fields
|
||||
update!(slave_element, "contact area", time => contact_area)
|
||||
end
|
||||
|
||||
for slave_element in slave_elements
|
||||
assemble!(problem, slave_element, time)
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
S = sort(collect(keys(normals))) # slave element nodes
|
||||
@@ -196,13 +502,87 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
|
||||
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)
|
||||
|
||||
maxdim = maximum(size(C1))
|
||||
if problem.properties.alpha != 0.0
|
||||
debug("mortar_3d: size C1 = ", size(C1), " max dim = $maxdim")
|
||||
debug("alpha != 0.0, applying transformation D = Dh*T^-1")
|
||||
alp = problem.properties.alpha
|
||||
Te = [
|
||||
1.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 1.0 0.0 0.0 0.0
|
||||
alp alp 0.0 1.0-2*alp 0.0 0.0
|
||||
0.0 alp alp 0.0 1.0-2*alp 0.0
|
||||
alp 0.0 alp 0.0 0.0 1.0-2*alp
|
||||
]
|
||||
invTe = [
|
||||
1.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 1.0 0.0 0.0 0.0
|
||||
-alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0 0.0
|
||||
0.0 -alp/(1-2*alp) -alp/(1-2*alp) 0.0 1/(1-2*alp) 0.0
|
||||
-alp/(1-2*alp) 0.0 -alp/(1-2*alp) 0.0 0.0 1/(1-2*alp)
|
||||
]
|
||||
# construct global transformation matrices T and invT
|
||||
T = SparseMatrixCOO()
|
||||
invT = SparseMatrixCOO()
|
||||
for element in slave_elements
|
||||
dofs = get_gdofs(problem, element)
|
||||
for i=1:field_dim
|
||||
ldofs = dofs[i:field_dim:end]
|
||||
add!(T, ldofs, ldofs, Te)
|
||||
add!(invT, ldofs, ldofs, invTe)
|
||||
end
|
||||
end
|
||||
T = sparse(T, maxdim, maxdim, (a, b) -> b)
|
||||
invT = sparse(invT, maxdim, maxdim, (a, b) -> b)
|
||||
# fill diagonal
|
||||
d = ones(size(T, 1))
|
||||
d[get_nonzero_rows(T)] = 0.0
|
||||
T += spdiagm(d)
|
||||
invT += spdiagm(d)
|
||||
#invT2 = sparse(inv(full(T)))
|
||||
#info("invT == invT2? ", invT == invT2)
|
||||
#maxabsdiff = maximum(abs(invT - invT2))
|
||||
#info("max diff = $maxabsdiff")
|
||||
C1 = C1*invT
|
||||
C2 = C2*invT
|
||||
end
|
||||
|
||||
tol = problem.properties.drop_tolerance
|
||||
debug("Dropping small values from C1 & C2, tolerace = $tol")
|
||||
SparseArrays.droptol!(C1, tol)
|
||||
SparseArrays.droptol!(C2, tol)
|
||||
|
||||
for j in S
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
weighted_gap[j] = g[dofs]
|
||||
end
|
||||
|
||||
state = problem.properties.contact_state_in_first_iteration
|
||||
if problem.properties.iteration == 1
|
||||
info("First contact iteration, initial contact state = $state")
|
||||
|
||||
if state == :AUTO
|
||||
avg_gap = mean([weighted_gap[j][1] for j in S])
|
||||
std_gap = std([weighted_gap[j][1] for j in S])
|
||||
if (avg_gap < 1.0e-12) && (std_gap < 1.0e-12)
|
||||
state = :ACTIVE
|
||||
else
|
||||
state = :UNKNOWN
|
||||
end
|
||||
info("Average weighted gap = $avg_gap, std gap = $std_gap, automatically determined contact state = $state")
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
# active / inactive node detection
|
||||
for j in S
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
@@ -218,7 +598,8 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
contact_pressure[j] = [0.0, 0.0, 0.0]
|
||||
end
|
||||
complementarity_condition[j] = contact_pressure[j] - weighted_gap[j]
|
||||
if complementarity_condition[j][1] < 0
|
||||
|
||||
if complementarity_condition[j][1] < 0.0
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
@@ -231,45 +612,49 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
end
|
||||
end
|
||||
|
||||
if "weighted gap" in props.store_fields
|
||||
update!(slave_elements, "weighted gap", time => weighted_gap)
|
||||
if (problem.properties.iteration == 1) && (state == :ACTIVE)
|
||||
for j in S
|
||||
is_inactive[j] = 0
|
||||
is_active[j] = 1
|
||||
is_slip[j] = 1
|
||||
is_stick[j] = 0
|
||||
end
|
||||
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)
|
||||
|
||||
if (problem.properties.iteration == 1) && (state == :INACTIVE)
|
||||
for j in S
|
||||
is_inactive[j] = 1
|
||||
is_active[j] = 0
|
||||
is_slip[j] = 0
|
||||
is_stick[j] = 0
|
||||
end
|
||||
end
|
||||
|
||||
#=
|
||||
info("# | active | inactive | stick | slip | gap | pres | comp")
|
||||
info("# | active | 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))"
|
||||
str1 = "$j | $(is_active[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
|
||||
=#
|
||||
|
||||
# solve variational inequality
|
||||
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
for j in S
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
if is_inactive[j] == 1
|
||||
debug("$j is inactive, removing dofs $dofs")
|
||||
C1[dofs,:] = 0.0
|
||||
C2[dofs,:] = 0.0
|
||||
D[dofs,:] = 0.0
|
||||
g[dofs,:] = 0.0
|
||||
end
|
||||
end
|
||||
|
||||
# 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")
|
||||
debug("$j is in active/slip, removing tangential constraints $tdofs")
|
||||
C2[tdofs,:] = 0.0
|
||||
g[tdofs] = 0.0
|
||||
normal = normals[j]
|
||||
@@ -279,22 +664,9 @@ function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
end
|
||||
end
|
||||
|
||||
# remove inactive nodes from assembly
|
||||
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")
|
||||
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
|
||||
|
||||
|
||||
+26
-3
@@ -13,13 +13,15 @@ type Solver{S<:AbstractSolver}
|
||||
initialized :: Bool
|
||||
u :: Vector{Float64}
|
||||
la :: Vector{Float64}
|
||||
alpha :: Float64 # generalized alpha time integration coefficient
|
||||
fields :: Dict{AbstractString, Field}
|
||||
properties :: S
|
||||
end
|
||||
|
||||
|
||||
function Solver{S<:AbstractSolver}(::Type{S}, name="solver", properties...)
|
||||
variant = S(properties...)
|
||||
solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], variant)
|
||||
solver = Solver{S}(name, 0.0, [], [], 0, nothing, false, [], [], 0.0, Dict(), variant)
|
||||
return solver
|
||||
end
|
||||
|
||||
@@ -33,11 +35,11 @@ function get_problems(solver::Solver)
|
||||
return solver.problems
|
||||
end
|
||||
|
||||
function push!(solver::Solver, problem)
|
||||
function push!(solver::Solver, problem::Problem)
|
||||
push!(solver.problems, problem)
|
||||
end
|
||||
|
||||
function getindex(solver::Solver, problem_name)
|
||||
function getindex(solver::Solver, problem_name::String)
|
||||
for problem in get_problems(solver)
|
||||
if problem.name == problem_name
|
||||
return problem
|
||||
@@ -46,6 +48,10 @@ function getindex(solver::Solver, problem_name)
|
||||
throw(KeyError(problem_name))
|
||||
end
|
||||
|
||||
function haskey(solver::Solver, field_name::String)
|
||||
return haskey(solver.fields, field_name)
|
||||
end
|
||||
|
||||
# one-liner helpers to identify problem types
|
||||
|
||||
is_field_problem(problem) = false
|
||||
@@ -312,6 +318,23 @@ function solve!(solver::Solver; empty_assemblies_before_solution=true, symmetric
|
||||
end
|
||||
gc()
|
||||
end
|
||||
|
||||
if !haskey(solver, "fint")
|
||||
solver.fields["fint"] = Field(time => f)
|
||||
else
|
||||
update!(solver.fields["fint"], time => f)
|
||||
end
|
||||
|
||||
fint = solver.fields["fint"]
|
||||
|
||||
if length(fint) > 1
|
||||
# kick in generalized alpha rule for time integration
|
||||
alpha = solver.alpha
|
||||
debug("Using generalized-α time integration, α=$alpha")
|
||||
K = (1-alpha)*K
|
||||
C1 = (1-alpha)*C1
|
||||
f = (1-alpha)*f + alpha*fint[end-1].data
|
||||
end
|
||||
|
||||
ndofs = solver.ndofs
|
||||
u = zeros(ndofs)
|
||||
|
||||
Reference in New Issue
Block a user