mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-28 20:46:58 +00:00
mortar contact tests
This commit is contained in:
@@ -10,52 +10,7 @@ using JuliaFEM.Core: Seg2, Quad4, Hex8, LinearElasticityProblem, get_connectivit
|
||||
assemble, PlaneStressLinearElasticityProblem, DirichletProblem,
|
||||
LinearSolver
|
||||
using JuliaFEM.Preprocess: aster_parse_nodes
|
||||
using JuliaFEM.Core: PlaneStressLinearElasticPlasticProblem
|
||||
|
||||
function test_plane_stress_linear_elasticity_with_surface_load()
|
||||
nodes = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 1.0],
|
||||
4 => [0.0, 1.0])
|
||||
|
||||
function set_geometry!(element, nodes)
|
||||
element["geometry"] = Vector{Float64}[nodes[i] for i in get_connectivity(element)]
|
||||
end
|
||||
element1 = Quad4([1, 2, 3, 4])
|
||||
set_geometry!(element1, nodes)
|
||||
element1["youngs modulus"] = 9000.0
|
||||
element1["poissons ratio"] = 0.25
|
||||
|
||||
element2 = Seg2([3, 4])
|
||||
set_geometry!(element2, nodes)
|
||||
element2["displacement traction force"] = Vector{Float64}[[0.0, -100.0] for i=1:2]
|
||||
|
||||
problem = PlaneStressLinearElasticityProblem()
|
||||
push!(problem, element1)
|
||||
push!(problem, element2)
|
||||
|
||||
free_dofs = Int64[3, 5, 6, 8]
|
||||
|
||||
ass = assemble(problem, 0.0)
|
||||
f = full(ass.force_vector)
|
||||
K = full(ass.stiffness_matrix)
|
||||
# info("initial force vector")
|
||||
# dump(reshape(f, 2, 4))
|
||||
# info("initial stiffness matrix")
|
||||
# dump(round(Int, K)[free_dofs, free_dofs])
|
||||
|
||||
u = zeros(2, 4)
|
||||
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
|
||||
|
||||
info("result vector")
|
||||
dump(u)
|
||||
# verified using Code Aster.
|
||||
# 2015-10-22-plane-stress/cplan_linear_traction_force.*
|
||||
@test isapprox(u[:,3], [2.77777777777778E-03, -1.11111111111111E-02])
|
||||
end
|
||||
|
||||
#test_plane_stress_linear_elasticity_with_surface_load()
|
||||
using JuliaFEM.Core: PlaneStressLinearElasticPlasticProblem
|
||||
|
||||
function test_plane_stress_linear_elasticplastic_with_surface_load()
|
||||
nodes = Dict{Int64, Vector{Float64}}(
|
||||
|
||||
@@ -88,8 +88,6 @@ function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:primal
|
||||
dim = 12
|
||||
C1 = sparse(assembly.C1, dim, dim)
|
||||
C2 = sparse(assembly.C2, dim, dim)
|
||||
# C1[[7,8], :] = 0
|
||||
# C2[[7,8], :] = 0
|
||||
|
||||
info("PDASS: constraint matrix C1")
|
||||
# dump(round(full(C1[1:2:end,1:2:end]), 3))
|
||||
@@ -109,7 +107,8 @@ function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:primal
|
||||
resize!(X, 12)
|
||||
resize!(u, 12)
|
||||
x = X+u
|
||||
info("x = \n", reshape(round(x, 2), 2, 6))
|
||||
info("x")
|
||||
dump(reshape(round(x, 2), 2, 6))
|
||||
P = sparse(eye(dim))
|
||||
C1 = P*C1
|
||||
C2 = P*C2
|
||||
@@ -142,19 +141,24 @@ function assemble_postprocess!(assembly, problem, time::Real, ::Type{Val{:primal
|
||||
continue
|
||||
end
|
||||
info("dof $i in active set")
|
||||
#g[i] = -gn[i]*4/6
|
||||
g[i] = -gn[i]
|
||||
# g[i] = -gn[i]
|
||||
else
|
||||
info("dof $i not in active set")
|
||||
C1[i,:] = 0
|
||||
C2[i,:] = 0
|
||||
# C1[i,:] = 0
|
||||
# C2[i,:] = 0
|
||||
end
|
||||
end
|
||||
|
||||
g[1] = 2.0
|
||||
g[3] = 2.0
|
||||
#=
|
||||
for i=2:2:dim
|
||||
C1[i,:] = 0
|
||||
C2[i,:] = 0
|
||||
end
|
||||
=#
|
||||
d = [7, 8]
|
||||
C2[d, :] = 0
|
||||
|
||||
assembly.g = sparse(g)
|
||||
assembly.C1 = C1
|
||||
@@ -228,7 +232,7 @@ end
|
||||
# remove tangential direction constraints
|
||||
# add_postprocessor!(cont, :remove_tangential_constraints)
|
||||
# apply PDASS
|
||||
# add_postprocessor!(cont, :primal_dual_active_set_strategy)
|
||||
add_postprocessor!(cont, :primal_dual_active_set_strategy)
|
||||
|
||||
@debug begin
|
||||
info("fel1.fields = $(fel1.fields)")
|
||||
|
||||
Reference in New Issue
Block a user