mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-10-03 22:57:57 +00:00
Jukka Aho
48006a144d
* 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
93 lines
3.5 KiB
Julia
93 lines
3.5 KiB
Julia
# 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
|
|
|
|
datadir = first(splitext(basename(@__FILE__)))
|
|
|
|
function get_model()
|
|
meshfile = joinpath(datadir, "block_2d.med")
|
|
mesh = aster_read_mesh(meshfile)
|
|
println(mesh.nodes[1])
|
|
|
|
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
|
|
end
|
|
|
|
@testset "finite sliding 2d patch test, linear Seg2 elements, standard basis" begin
|
|
|
|
solver = get_model()
|
|
interface = solver["interface"]
|
|
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
|
|
|
|
for time in [0.0, 1/3, 2/3, 1.0]
|
|
interface.properties.iteration = 1
|
|
solver.time = time
|
|
solver()
|
|
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, reaction_force = get_nodal_vector(get_slave_elements(interface), "reaction force", 1.0)
|
|
u2 = [u[2] for u in displacement]
|
|
f2 = [f[2] for f in reaction_force]
|
|
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)
|
|
end
|