mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-20 18:18:31 +00:00
Move contact mechanics to separate package (#196)
Development of auto-differentiated mortar contact mechanics in 2D is moved to own separate package, MortarContact2DAD. Other changes are similar to what is done with MortarContact2D: elements are added to problems using `add_slave_elements!` and `add_master_elements!` instead of `add_elements!`, to make interface more explicit. Also, problem name is `Contact2DAD`, so the dimension is now explicitly stated in problem name. (Also have `Mortar2DAD`, compare to the `Mortar2D` and `Contact2D` of `MortarContact2D.jl`.)
This commit is contained in:
Jukka Aho
@@ -13,3 +13,4 @@ FEMQuad
|
||||
Reexport
|
||||
HeatTransfer
|
||||
MortarContact2D
|
||||
MortarContact2DAD
|
||||
|
||||
+2
-2
@@ -47,10 +47,10 @@ export assemble!, postprocess!
|
||||
### Mortar methods ###
|
||||
|
||||
@reexport using MortarContact2D
|
||||
@reexport using MortarContact2DAD
|
||||
|
||||
include("problems_mortar.jl")
|
||||
include("problems_mortar_3d.jl")
|
||||
include("problems_mortar_2d_autodiff.jl")
|
||||
export calculate_normals, calculate_normals!, project_from_slave_to_master,
|
||||
project_from_master_to_slave, Mortar, get_slave_elements,
|
||||
get_polygon_clip
|
||||
@@ -67,7 +67,7 @@ include("solvers_modal.jl")
|
||||
export Modal
|
||||
include("problems_contact.jl")
|
||||
include("problems_contact_3d.jl")
|
||||
include("problems_contact_2d_autodiff.jl")
|
||||
#include("problems_contact_3d_autodiff.jl")
|
||||
export Contact
|
||||
|
||||
# Preprocess module
|
||||
|
||||
@@ -1,380 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using ForwardDiff
|
||||
|
||||
""" Find segment from slave element corresponding to master element nodes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x1_, n1_
|
||||
slave element geometry and normal direction
|
||||
x2
|
||||
master element node to project onto slave
|
||||
|
||||
Returns
|
||||
-------
|
||||
xi
|
||||
dimensionless coordinate on slave corresponding to
|
||||
projected master
|
||||
|
||||
"""
|
||||
function project_from_master_to_slave_ad{E<:MortarElements2D}(
|
||||
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
|
||||
tol=1.0e-10, max_iterations=20, debug=false)
|
||||
|
||||
""" Multiply basis / dbasis at `xi` with field. """
|
||||
function mul(func, xi, field)
|
||||
B = func(slave_element, [xi], time)
|
||||
return sum(B[i]*field[i] for i=1:length(B))
|
||||
end
|
||||
|
||||
x1(xi1) = mul(get_basis, xi1, x1_)
|
||||
dx1(xi1) = mul(get_dbasis, xi1, x1_)
|
||||
n1(xi1) = mul(get_basis, xi1, n1_)
|
||||
dn1(xi1) = mul(get_dbasis, xi1, n1_)
|
||||
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
|
||||
xi1_next = 0.0
|
||||
dxi1 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi1 = -R(xi1)/dR(xi1)
|
||||
dxi1 = clamp.(dxi1, -0.3, 0.3)
|
||||
xi1_next = clamp.(xi1 + dxi1, -1.0, 1.0)
|
||||
if norm(xi1_next - xi1) < tol
|
||||
return xi1_next
|
||||
end
|
||||
if debug
|
||||
info("xi1 = $xi1")
|
||||
info("R(xi1) = $(R(xi1))")
|
||||
info("dR(xi1) = $(dR(xi1))")
|
||||
info("dxi1 = $dxi1")
|
||||
info("norm = $(norm(xi1_next - xi1))")
|
||||
info("xi1_next = $xi1_next")
|
||||
end
|
||||
xi1 = xi1_next
|
||||
end
|
||||
|
||||
info("x1 = $x1_")
|
||||
info("n1 = $n1_")
|
||||
info("x2 = $x2")
|
||||
info("xi1 = $xi1, dxi1 = $dxi1")
|
||||
info("-R(xi1) = $(-R(xi1))")
|
||||
info("dR(xi1) = $(dR(xi1))")
|
||||
error("find projection from master to slave: did not converge")
|
||||
|
||||
end
|
||||
|
||||
function project_from_slave_to_master_ad{E<:MortarElements2D}(
|
||||
master_element::Element{E}, x1, n1, x2_;
|
||||
tol=1.0e-10, max_iterations=20)
|
||||
|
||||
x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
|
||||
dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
|
||||
cross2(a, b) = cross([a; 0], [b; 0])[3]
|
||||
R(xi2) = cross2(x2(xi2)-x1, n1)
|
||||
dR(xi2) = cross2(dx2(xi2), n1)
|
||||
|
||||
xi2 = 0.0
|
||||
dxi2 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < tol
|
||||
return xi2
|
||||
end
|
||||
end
|
||||
|
||||
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
|
||||
|
||||
end
|
||||
|
||||
"""
|
||||
Frictionless 2d finite sliding contact with forwarddiff.
|
||||
|
||||
true/false flags: finite_sliding, friction, use_forwarddiff
|
||||
"""
|
||||
function assemble!(problem::Problem{Contact}, time::Float64,
|
||||
::Type{Val{1}}, ::Type{Val{true}},
|
||||
::Type{Val{false}}, ::Type{Val{true}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_elements = get_slave_elements(problem)
|
||||
|
||||
function calculate_interface(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 = zeros(u)
|
||||
gap = zeros(u)
|
||||
C = zeros(la)
|
||||
S = Set{Int64}()
|
||||
|
||||
# 1. update nodal normals for slave elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in slave_elements
|
||||
conn = get_connectivity(element)
|
||||
push!(S, conn...)
|
||||
gdofs = get_gdofs(problem, element)
|
||||
X_el = element("geometry", time)
|
||||
x_el = tuple( (X_el[i] + u[:,j] for (i,j) in enumerate(conn))... )
|
||||
#=
|
||||
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])
|
||||
end
|
||||
# swap element normals in 2d if they point to inside of body
|
||||
if props.rotate_normals
|
||||
for i=1:size(normals,2)
|
||||
normals[:,i] = -normals[:,i]
|
||||
end
|
||||
end
|
||||
normals2 = Dict()
|
||||
for j in S
|
||||
normals2[j] = normals[:,j]
|
||||
end
|
||||
update!(slave_elements, "normal", time => normals2)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = ((u[:,i] for i in slave_element_nodes)...)
|
||||
x1 = ((Xi+ui for (Xi,ui) in zip(X1,u1))...)
|
||||
la1 = ((la[:,i] for i in slave_element_nodes)...)
|
||||
n1 = ((normals[:,i] for i in slave_element_nodes)...)
|
||||
nnodes = size(slave_element, 2)
|
||||
|
||||
# construct dual basis
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = ((u[:,i] for i in master_element_nodes)...)
|
||||
x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
|
||||
xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
|
||||
xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
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)
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = ((u[:,i] for i in master_element_nodes)...)
|
||||
x2 = ((Xi+ui for (Xi,ui) in zip(X2,u2))...)
|
||||
|
||||
#x1_midpoint = 1/2*(x1[1]+x1[2])
|
||||
#x2_midpoint = 1/2*(x2[1]+x2[2])
|
||||
#distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
|
||||
#distance > props.maximum_distance && continue
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
xi1a = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[1])
|
||||
xi1b = project_from_master_to_slave_ad(slave_element, field(x1), field(n1), x2[2])
|
||||
xi1 = clamp.([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
slave_dofs = get_gdofs(problem, slave_element)
|
||||
master_dofs = get_gdofs(problem, master_element)
|
||||
|
||||
# 4. loop integration points of segment
|
||||
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
|
||||
|
||||
# project gauss point from slave element to master element
|
||||
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))
|
||||
x_s = interpolate(N1, x1) # coordinate in gauss point
|
||||
n_s = interpolate(N1, n1) # normal direction in gauss point
|
||||
t_s = Q'*n_s # tangent direction in gauss point
|
||||
xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
x_m = interpolate(N2, x2)
|
||||
Phi = Ae*N1
|
||||
|
||||
la_s = interpolate(Phi, la1) # traction force in gauss point
|
||||
gn = -dot(n_s, x_s - x_m) # normal gap
|
||||
|
||||
fc[:,slave_element_nodes] += w*la_s*N1'
|
||||
fc[:,master_element_nodes] -= w*la_s*N2'
|
||||
gap[1,slave_element_nodes] += w*gn*Phi
|
||||
#gap[1,slave_element_nodes] += w*gn*N1'
|
||||
|
||||
end # done integrating segment
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # slave elements done
|
||||
|
||||
# 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
|
||||
|
||||
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
|
||||
|
||||
is_active = Dict{Int, Bool}()
|
||||
condition = Dict()
|
||||
|
||||
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])
|
||||
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
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)
|
||||
b = sparse(b)
|
||||
SparseArrays.droptol!(A, 1.0e-9)
|
||||
SparseArrays.droptol!(b, 1.0e-9)
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
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
|
||||
|
||||
#=
|
||||
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
|
||||
|
||||
@@ -1,247 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using ForwardDiff
|
||||
|
||||
const MortarElements2D = Union{Seg2,Seg3}
|
||||
|
||||
# forwarddiff version of mesh tying in 2d
|
||||
|
||||
function project_from_master_to_slave_ad{E<:MortarElements2D}(
|
||||
slave_element::Element{E}, x1_, n1_, x2, time;
|
||||
tol=1.0e-10, max_iterations=20)
|
||||
|
||||
x1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), x1_)
|
||||
dx1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), x1_)
|
||||
n1(xi1) = interpolate(vec(get_basis(slave_element, [xi1], time)), n1_)
|
||||
dn1(xi1) = interpolate(vec(get_dbasis(slave_element, [xi1], time)), n1_)
|
||||
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
|
||||
dxi1 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi1 = -R(xi1)/dR(xi1)
|
||||
xi1 += dxi1
|
||||
if norm(dxi1) < tol
|
||||
return xi1
|
||||
end
|
||||
end
|
||||
|
||||
info("x1 = $(ForwardDiff.get_value(x1_.data))")
|
||||
info("n1 = $(ForwardDiff.get_value(n1_.data))")
|
||||
info("x2 = $(ForwardDiff.get_value(x2))")
|
||||
info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
|
||||
info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
|
||||
info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
|
||||
error("find projection from master to slave: did not converge")
|
||||
|
||||
end
|
||||
|
||||
function project_from_slave_to_master_ad{E<:MortarElements2D}(
|
||||
master_element::Element{E}, x1, n1, x2_, time;
|
||||
tol=1.0e-10, max_iterations=20)
|
||||
|
||||
x2(xi2) = interpolate(vec(get_basis(master_element, [xi2], time)), x2_)
|
||||
dx2(xi2) = interpolate(vec(get_dbasis(master_element, [xi2], time)), x2_)
|
||||
cross2(a, b) = cross([a; 0], [b; 0])[3]
|
||||
R(xi2) = cross2(x2(xi2)-x1, n1)
|
||||
dR(xi2) = cross2(dx2(xi2), n1)
|
||||
|
||||
xi2 = 0.0
|
||||
dxi2 = 0.0
|
||||
for i=1:max_iterations
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < tol
|
||||
return xi2
|
||||
end
|
||||
end
|
||||
|
||||
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
|
||||
|
||||
end
|
||||
|
||||
""" 2d mesh tie using ForwardDiff.
|
||||
|
||||
Construct .. + fc*la and C(d,la)=0
|
||||
|
||||
"""
|
||||
function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{true}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_elements = get_slave_elements(problem)
|
||||
if field_name != "displacement"
|
||||
error("mortar forwarddiff assembly: only displacement field with adjust=yes supported")
|
||||
end
|
||||
|
||||
function calculate_interface(x::Vector)
|
||||
|
||||
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 = zeros(u)
|
||||
gap = zeros(u)
|
||||
C = zeros(la)
|
||||
|
||||
S = Set{Int64}()
|
||||
# 1. update nodal normals for slave elements
|
||||
tangents = zeros(u)
|
||||
for element in slave_elements
|
||||
conn = get_connectivity(element)
|
||||
push!(S, conn...)
|
||||
X1 = element("geometry", time)
|
||||
u1 = ((u[:,i] for i in conn)...)
|
||||
x1 = map(+, X1, u1)
|
||||
dN = get_dbasis(element, [0.0], time)
|
||||
tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
for nid in conn
|
||||
tangents[:,nid] += tangent[:]
|
||||
end
|
||||
end
|
||||
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for j in S
|
||||
tangents[:,j] /= norm(tangents[:,j])
|
||||
normals[:,j] = Q*tangents[:,j]
|
||||
end
|
||||
|
||||
if props.rotate_normals
|
||||
for j in S
|
||||
normals[:,j] = -normals[:,j]
|
||||
end
|
||||
end
|
||||
|
||||
#update!(slave_elements, "normal", time => Dict(j => normals[:,j] for j in S))
|
||||
#update!(slave_elements, "tangent", time => Dict(j => tangents[:,j] for j in S))
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
nsl = length(slave_element)
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = ((u[:,i] for i in slave_element_nodes)...)
|
||||
x1 = map(+, X1, u1)
|
||||
la1 = ((la[:,i] for i in slave_element_nodes)...)
|
||||
n1 = ((normals[:,i] for i in slave_element_nodes)...)
|
||||
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element("master elements", time)
|
||||
|
||||
nm = length(master_element)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = ((u[:,i] for i in master_element_nodes)...)
|
||||
x2 = map(+, X2, u2)
|
||||
|
||||
# 3.1 calculate segmentation
|
||||
xi1a = project_from_master_to_slave_ad(slave_element, x1, n1, x2[1], time)
|
||||
xi1b = project_from_master_to_slave_ad(slave_element, x1, n1, x2[2], time)
|
||||
# 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
|
||||
for ip in get_integration_points(slave_element, 3)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
#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 = 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 = interpolate(N1, x1) # coordinate in gauss point
|
||||
n_s = interpolate(N1, n1) # normal direction in gauss point
|
||||
#xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
|
||||
xi_m = project_from_slave_to_master_ad(master_element, x_s, n_s, x2, time)
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
x_m = interpolate(N2, x2)
|
||||
|
||||
la_s = interpolate(Phi, la1)
|
||||
gn = dot(n_s, x_s-x_m)
|
||||
|
||||
u_s = interpolate(N1, u1)
|
||||
u_m = interpolate(N2, u2)
|
||||
X_s = interpolate(N1, X1)
|
||||
X_m = interpolate(N2, X2)
|
||||
|
||||
fc[:,slave_element_nodes] += w*la_s*N1'
|
||||
fc[:,master_element_nodes] -= w*la_s*N2'
|
||||
#gap[1,slave_element_nodes] += w*gn*Phi'
|
||||
gap[:,slave_element_nodes] += w*(u_s-u_m)*Phi'
|
||||
if props.adjust
|
||||
G = w*(X_s-X_m)*Phi'
|
||||
gap[:,slave_element_nodes] += G
|
||||
end
|
||||
end
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
C = gap
|
||||
|
||||
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 = ForwardDiff.jacobian(calculate_interface, x)
|
||||
b = -calculate_interface(x)
|
||||
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseArrays.droptol!(A, 1.0e-12)
|
||||
SparseArrays.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]
|
||||
|
||||
empty!(problem.assembly)
|
||||
problem.assembly.K = K
|
||||
problem.assembly.C1 = C1
|
||||
problem.assembly.C2 = C2
|
||||
problem.assembly.D = D
|
||||
problem.assembly.f = f
|
||||
problem.assembly.g = g
|
||||
|
||||
end
|
||||
@@ -33,9 +33,9 @@ function Modal(nev=10, which=:SM)
|
||||
end
|
||||
|
||||
""" Eliminate Dirichlet boundary condition from matrices K, M. """
|
||||
function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Problem{Dirichlet}, ndim::Int)
|
||||
function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Problem{Dirichlet}, ndim::Int)
|
||||
K = sparse(problem.assembly.K, ndim, ndim)
|
||||
C1 = sparse(problem.assembly.C1, ndim, ndim)
|
||||
C2 = sparse(problem.assembly.C2, ndim, ndim)
|
||||
@@ -100,10 +100,10 @@ end
|
||||
|
||||
|
||||
""" Eliminate mesh tie constraints from matrices K, M. """
|
||||
function eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Union{Problem{Mortar}, Problem{Mortar2D}},
|
||||
ndim::Int)
|
||||
function FEMBase.eliminate_boundary_conditions!(K_red::SparseMatrixCSC,
|
||||
M_red::SparseMatrixCSC,
|
||||
problem::Union{Problem{Mortar}, Problem{Mortar2D}},
|
||||
ndim::Int)
|
||||
|
||||
C1 = sparse(problem.assembly.C1, ndim, ndim)
|
||||
C2 = sparse(problem.assembly.C2, ndim, ndim)
|
||||
|
||||
@@ -36,15 +36,13 @@ using JuliaFEM.Testing
|
||||
update!(bc_lower, "displacement 1", 0.0)
|
||||
update!(bc_lower, "displacement 2", 0.0)
|
||||
|
||||
contact = Problem(Contact, "contact between upper and lower block", 2, "displacement")
|
||||
contact = Problem(Contact2DAD, "contact between upper and lower block", 2, "displacement")
|
||||
contact.properties.rotate_normals = true
|
||||
contact.properties.finite_sliding = true
|
||||
contact.properties.friction = false
|
||||
contact.properties.use_forwarddiff = true
|
||||
contact_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
contact_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
update!(contact_slave_elements, "master elements", contact_master_elements)
|
||||
contact.elements = [contact_master_elements; contact_slave_elements]
|
||||
add_slave_elements!(contact, contact_slave_elements)
|
||||
add_master_elements!(contact, contact_master_elements)
|
||||
|
||||
nnodes = length(mesh.nodes)
|
||||
contact.assembly.u = zeros(2*nnodes)
|
||||
contact.assembly.la = zeros(2*nnodes)
|
||||
|
||||
@@ -1,178 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Testing
|
||||
|
||||
function get_model(::Type{Val{Symbol("mesh tie with curved 2d block")}};
|
||||
dy=0.0, adjust=false, tolerance=0.0, rotate_normals=false, swap=false,
|
||||
dual_basis=false, use_forwarddiff=true, finite_strain=false,
|
||||
geometric_stiffness=false)
|
||||
|
||||
meshfile = @__DIR__() * "/testdata/block_2d_curved.med"
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Elasticity, "upper", 2)
|
||||
upper.properties.formulation = :plane_stress
|
||||
upper.properties.finite_strain = finite_strain
|
||||
upper.properties.geometric_stiffness = geometric_stiffness
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper.elements, "youngs modulus", 96.0)
|
||||
update!(upper.elements, "poissons ratio", 1/3)
|
||||
|
||||
lower = Problem(Elasticity, "lower", 2)
|
||||
lower.properties.formulation = :plane_stress
|
||||
lower.properties.finite_strain = finite_strain
|
||||
lower.properties.geometric_stiffness = geometric_stiffness
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower.elements, "youngs modulus", 96.0)
|
||||
update!(lower.elements, "poissons ratio", 1/3)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "upper boundary", 2, "displacement")
|
||||
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
|
||||
update!(bc_upper.elements, "displacement 1", 0.0)
|
||||
update!(bc_upper.elements, "displacement 2", dy)
|
||||
|
||||
bc_lower = Problem(Dirichlet, "lower boundary", 2, "displacement")
|
||||
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
|
||||
update!(bc_lower.elements, "displacement 1", 0.0)
|
||||
update!(bc_lower.elements, "displacement 2", 0.0)
|
||||
|
||||
interface = Problem(Mortar, "interface between upper and lower block", 2, "displacement")
|
||||
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
if swap
|
||||
interface_slave_elements, interface_master_elements = interface_master_elements, interface_slave_elements
|
||||
end
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
interface.properties.adjust = adjust
|
||||
interface.properties.distval = tolerance
|
||||
interface.properties.rotate_normals = rotate_normals
|
||||
interface.properties.dual_basis = dual_basis
|
||||
interface.properties.use_forwarddiff = use_forwarddiff
|
||||
interface.assembly.u = zeros(2*length(mesh.nodes))
|
||||
interface.assembly.la = zeros(2*length(mesh.nodes))
|
||||
|
||||
solver = Solver(Linear)
|
||||
push!(solver, upper, lower, bc_upper, bc_lower, interface)
|
||||
|
||||
return solver
|
||||
|
||||
end
|
||||
|
||||
#=
|
||||
|
||||
@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=0.0" begin
|
||||
# TODO: analytical solution now known, verify using other fem software
|
||||
solver = get_model("mesh tie with curved 2d block";
|
||||
adjust=false, tolerance=10, dy=-0.1, rotate_normals=true,
|
||||
dual_basis=true, use_forwarddiff=true, finite_strain=false,
|
||||
geometric_stiffness=false)
|
||||
solver()
|
||||
interface = solver["interface between upper and lower block"]
|
||||
@test isapprox(norm(interface.assembly.u), 0.11339715157447851)
|
||||
end
|
||||
|
||||
|
||||
@testset "curved surface with adjust=true, dual lagrange, slave=lower surface, dy=0.0" begin
|
||||
# TODO: analytical solution now known, verify using other fem software
|
||||
solver = get_model("mesh tie with curved 2d block";
|
||||
adjust=true, tolerance=10, dy=0.0, rotate_normals=true,
|
||||
dual_basis=true, use_forwarddiff=true)
|
||||
solver()
|
||||
interface = solver["interface between upper and lower block"]
|
||||
@test solver.properties.iteration == 2
|
||||
# differs -- why?
|
||||
@test isapprox(norm(interface.assembly.u), 0.11660422877751599)
|
||||
end
|
||||
|
||||
@testset "curved surface with adjust=true, standard lagrange, slave=lower surface, dy=-0.1" begin
|
||||
# TODO: analytical solution now known, verify using other fem software
|
||||
solver = get_model("mesh tie with curved 2d block";
|
||||
adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
|
||||
dual_basis=false, use_forwarddiff=true)
|
||||
solver()
|
||||
interface = solver["interface between upper and lower block"]
|
||||
@test solver.properties.iteration == 2
|
||||
@test isapprox(norm(interface.assembly.u), 0.34230262165505887)
|
||||
end
|
||||
|
||||
@testset "curved surface, adjust=true, dual basis, slave=lower surface, dy=-0.1" begin
|
||||
# TODO: analytical solution now known, verify using other fem software
|
||||
solver = get_model("mesh tie with curved 2d block";
|
||||
adjust=true, tolerance=10, dy=-0.1, rotate_normals=true,
|
||||
dual_basis=true, use_forwarddiff=true)
|
||||
solver()
|
||||
interface = solver["interface between upper and lower block"]
|
||||
@test solver.properties.iteration == 2
|
||||
@test isapprox(norm(interface.assembly.u), 0.34318800698017704)
|
||||
end
|
||||
|
||||
=#
|
||||
|
||||
|
||||
|
||||
@testset "compare forwarddiff solution to normal" begin
|
||||
X = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [0.0, 1.0],
|
||||
4 => [1.0, 1.0])
|
||||
u = Dict(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [0.0, 0.0],
|
||||
3 => [0.0, 0.0],
|
||||
4 => [0.0, 0.0])
|
||||
sel1 = Element(Seg2, [1, 2])
|
||||
mel1 = Element(Seg2, [3, 4])
|
||||
update!([sel1, mel1], "geometry", X)
|
||||
update!([sel1, mel1], "displacement", u)
|
||||
update!(sel1, "master elements", [mel1])
|
||||
|
||||
p1 = Problem(Mortar2D, "test 1", 2, "displacement")
|
||||
add_slave_elements!(p1, [sel1])
|
||||
add_master_elements!(p1, [mel1])
|
||||
assemble!(p1, 0.0)
|
||||
|
||||
p2 = Problem(Mortar, "test 2", 2, "displacement")
|
||||
push!(p2, sel1, mel1)
|
||||
p2.properties.use_forwarddiff = true
|
||||
p2.assembly.u = zeros(8)
|
||||
p2.assembly.la = zeros(8)
|
||||
assemble!(p2, 0.0)
|
||||
|
||||
@test isapprox(p1.assembly, p2.assembly)
|
||||
|
||||
#=
|
||||
empty!(p1.assembly)
|
||||
empty!(p2.assembly)
|
||||
p1.properties.adjust = true
|
||||
p2.properties.adjust = true
|
||||
assemble!(p1, 0.0)
|
||||
assemble!(p2, 0.0)
|
||||
C11 = full(p1.assembly.C1, 4, 8)
|
||||
C12 = full(p2.assembly.C1, 4, 8)
|
||||
C21 = full(p1.assembly.C2, 4, 8)
|
||||
C22 = full(p2.assembly.C2, 4, 8)
|
||||
D1 = full(p1.assembly.D)
|
||||
D2 = full(p2.assembly.D)
|
||||
g1 = full(p1.assembly.g, 4, 1)
|
||||
g2 = full(p2.assembly.g, 4, 1)
|
||||
println("C1")
|
||||
dump(C11)
|
||||
dump(C12)
|
||||
println("C2")
|
||||
dump(C21)
|
||||
dump(C22)
|
||||
println("D")
|
||||
dump(D1)
|
||||
dump(D2)
|
||||
println("g")
|
||||
dump(g1)
|
||||
dump(g2)
|
||||
@test isapprox(p1.assembly, p2.assembly)
|
||||
=#
|
||||
end
|
||||
@@ -6,88 +6,77 @@ using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Testing
|
||||
|
||||
datadir = first(splitext(basename(@__FILE__)))
|
||||
pkg_dir = Pkg.dir("JuliaFEM")
|
||||
datadir = joinpath(pkg_dir, "test", first(splitext(basename(@__FILE__))))
|
||||
|
||||
function get_model()
|
||||
meshfile = joinpath(datadir, "block_2d.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
#error("mesh has $(length(mesh.nodes)) nodes")
|
||||
meshfile = joinpath(datadir, "block_2d.med")
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(mesh, Elasticity, "UPPER", 2)
|
||||
lower = Problem(mesh, Elasticity, "LOWER", 2)
|
||||
upper = Problem(mesh, Elasticity, "UPPER", 2)
|
||||
lower = Problem(mesh, Elasticity, "LOWER", 2)
|
||||
|
||||
for body in [upper, lower]
|
||||
body.properties.formulation = :plane_stress
|
||||
update!(body, "youngs modulus", 288.0)
|
||||
update!(body, "poissons ratio", 1/3)
|
||||
end
|
||||
|
||||
load = Problem(mesh, Elasticity, "UPPER_TOP", 2)
|
||||
load.properties.formulation = :plane_stress
|
||||
update!(load, "displacement traction force 2", 0.0 => 0.0)
|
||||
update!(load, "displacement traction force 2", 1.0 => -28.8)
|
||||
bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement")
|
||||
update!(bc1, "displacement 2", 0.0)
|
||||
bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement")
|
||||
update!(bc2, "displacement 1", 0.0)
|
||||
bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement")
|
||||
update!(bc3, "displacement 1", 0.0)
|
||||
|
||||
interface = Problem(Contact, "interface", 2, "displacement")
|
||||
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
interface.properties.rotate_normals = true
|
||||
|
||||
# in LOWER_LEFT we have node belonging also to contact interface
|
||||
# let's remove it from dirichlet bc
|
||||
create_node_set_from_element_set!(mesh, "LOWER_LEFT")
|
||||
nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT")
|
||||
coords = mesh.nodes[nid]
|
||||
info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords")
|
||||
dofs = [2*(nid-1)+1, 2*(nid-1)+2]
|
||||
info("removing nid $nid, dofs $dofs from LOWER_LEFT")
|
||||
push!(bc2.assembly.removed_dofs, dofs...)
|
||||
|
||||
solver = Solver(Nonlinear)
|
||||
push!(solver, upper, lower, load, bc1, bc2, bc3, interface)
|
||||
return solver
|
||||
for body in [upper, lower]
|
||||
body.properties.formulation = :plane_stress
|
||||
update!(body, "youngs modulus", 288.0)
|
||||
update!(body, "poissons ratio", 1/3)
|
||||
end
|
||||
|
||||
@testset "finite sliding 2d patch test, linear Seg2 elements, standard basis" begin
|
||||
load = Problem(mesh, Elasticity, "UPPER_TOP", 2)
|
||||
load.properties.formulation = :plane_stress
|
||||
update!(load, "displacement traction force 2", 0.0 => 0.0)
|
||||
update!(load, "displacement traction force 2", 1.0 => -28.8)
|
||||
bc1 = Problem(mesh, Dirichlet, "LOWER_BOTTOM", 2, "displacement")
|
||||
update!(bc1, "displacement 2", 0.0)
|
||||
bc2 = Problem(mesh, Dirichlet, "LOWER_LEFT", 2, "displacement")
|
||||
update!(bc2, "displacement 1", 0.0)
|
||||
bc3 = Problem(mesh, Dirichlet, "UPPER_LEFT", 2, "displacement")
|
||||
update!(bc3, "displacement 1", 0.0)
|
||||
|
||||
solver = get_model()
|
||||
interface = solver["interface"]
|
||||
interface.assembly.u = zeros(48)
|
||||
interface.assembly.la = zeros(48)
|
||||
upper = solver["UPPER"]
|
||||
lower = solver["LOWER"]
|
||||
for body in [upper, lower]
|
||||
body.properties.geometric_stiffness = true
|
||||
body.properties.finite_strain = true
|
||||
end
|
||||
interface.properties.finite_sliding = true
|
||||
interface.properties.use_forwarddiff = true
|
||||
interface = Problem(Contact2DAD, "interface", 2, "displacement")
|
||||
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
add_slave_elements!(interface, interface_slave_elements)
|
||||
add_master_elements!(interface, interface_master_elements)
|
||||
interface.properties.rotate_normals = true
|
||||
|
||||
for time in [0.0, 1/3, 2/3, 1.0]
|
||||
interface.properties.iteration = 1
|
||||
solve!(solver, time)
|
||||
end
|
||||
# in LOWER_LEFT we have node belonging also to contact interface
|
||||
# let's remove it from dirichlet bc
|
||||
create_node_set_from_element_set!(mesh, "LOWER_LEFT")
|
||||
nid = find_nearest_node(mesh, [0.0, 0.5]; node_set="LOWER_LEFT")
|
||||
coords = mesh.nodes[nid]
|
||||
info("nearest node to (0.0, 0.5) = $nid, coordinates = $coords")
|
||||
dofs = [2*(nid-1)+1, 2*(nid-1)+2]
|
||||
info("removing nid $nid, dofs $dofs from LOWER_LEFT")
|
||||
push!(bc2.assembly.removed_dofs, dofs...)
|
||||
|
||||
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0)
|
||||
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0)
|
||||
node_ids, la = get_nodal_vector(get_slave_elements(interface), "lambda", 1.0)
|
||||
u2 = [u[2] for u in displacement]
|
||||
f2 = [f[2] for f in la]
|
||||
maxabsu2 = maximum(abs.(u2))
|
||||
stdabsu2 = std(abs.(u2))
|
||||
info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
|
||||
@test isapprox(stdabsu2, 0.0; atol=1.0e-12)
|
||||
maxabsf2 = maximum(abs.(f2))
|
||||
stdabsf2 = std(abs.(f2))
|
||||
info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
|
||||
@test isapprox(stdabsf2, 0.0; atol=1.0e-12)
|
||||
# for linear case pressure 28.8
|
||||
@test isapprox(mean(abs.(f2)), 27.76616800689944; rtol=1.0e-3)
|
||||
solver = Solver(Nonlinear)
|
||||
push!(solver, upper, lower, load, bc1, bc2, bc3, interface)
|
||||
|
||||
interface.assembly.u = zeros(48)
|
||||
interface.assembly.la = zeros(48)
|
||||
|
||||
for body in [upper, lower]
|
||||
body.properties.geometric_stiffness = true
|
||||
body.properties.finite_strain = true
|
||||
end
|
||||
|
||||
for time in [0.0, 1/3, 2/3, 1.0]
|
||||
interface.properties.iteration = 0
|
||||
solve!(solver, time)
|
||||
end
|
||||
|
||||
node_ids, displacement = get_nodal_vector(interface.elements, "displacement", 1.0)
|
||||
node_ids, geometry = get_nodal_vector(interface.elements, "geometry", 1.0)
|
||||
node_ids, la = get_nodal_vector(interface.elements, "lambda", 1.0)
|
||||
u2 = [u[2] for u in displacement]
|
||||
f2 = [f[2] for f in la]
|
||||
maxabsu2 = maximum(abs.(u2))
|
||||
stdabsu2 = std(abs.(u2))
|
||||
info("max(abs(u2)) = $maxabsu2, std(abs(u2)) = $stdabsu2")
|
||||
@test isapprox(stdabsu2, 0.0; atol=1.0e-12)
|
||||
maxabsf2 = maximum(abs.(f2))
|
||||
stdabsf2 = std(abs.(f2))
|
||||
info("max(abs(f2)) = $maxabsf2, std(abs(f2)) = $stdabsf2")
|
||||
@test isapprox(stdabsf2, 0.0; atol=1.0e-12)
|
||||
# for linear case pressure 28.8
|
||||
@test isapprox(mean(abs.(f2)), 27.76616800689944; rtol=1.0e-3)
|
||||
|
||||
Reference in New Issue
Block a user