mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-30 05:19:19 +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
@@ -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