Author: Jukka Aho
Abstract: 2d tie contact.
Each element is modelled as own "body" and they are connected using tie contacts. Segments 5-6 and 9-10 and 6-7 are slave surfaces, so node 6 or 9 is on at least two tie contacts as slave node. Moreover this model has dirichlet boundary $y=0$ at bottom of body 1 and $x=0$ on left. To get the accurate solution one needs to minimize \begin{equation} \frac{15}{2}u_1^4 + 60 u_1^3 + \frac{15}{4}u_1^2 u_2^2 + 15 u_1^2 u_2 + 120 u_1^2 + 15 u_1 u_2^2 + 60 u_1 u_2 + \frac{15}{2}u_2^4 + 60 u_2^3 + 120 u_2^2 + 50 u_2, \end{equation} which gives approximate $u_1 = 0.0634862$ and $u_2 = -0.277183$ for the displacement of upper right corner. Wolfram.
In [1]:
using JuliaFEM.Core: Element, Seg2, Quad4, PlaneStressElasticityProblem,
DirichletProblem, MortarProblem, DirectSolver
In [10]:
nodes = Dict{Int64, Vector{Float64}}(
1 => [0.0, 0.0], 2 => [2.0, 0.0],
3 => [2.0, 1.0], 4 => [0.0, 1.0],
5 => [0.0, 1.0], 6 => [1.0, 1.0],
7 => [1.0, 2.0], 8 => [0.0, 2.0],
9 => [1.0, 1.0], 10 => [2.0, 1.0],
11 => [2.0, 2.0], 12 => [1.0, 2.0]);
In [11]:
connectivity = Dict{Int64, Vector{Int64}}(
1 => [1, 2, 3, 4],
2 => [5, 6, 7, 8],
3 => [9, 10, 11, 12]);
In [13]:
elements = Element[]
for c in values(connectivity)
element = Quad4(c)
element["geometry"] = Vector{Float64}[nodes[i] for i in c]
element["youngs modulus"] = 900.0
element["poissons ratio"] = 0.25
push!(elements, element)
end
nelements = length(elements)
info("number of elements: $nelements")
Create three bodies, each containing one element.
In [14]:
body1 = PlaneStressElasticityProblem("body 1")
body2 = PlaneStressElasticityProblem("body 2")
body3 = PlaneStressElasticityProblem("body 3")
push!(body1, elements[1])
push!(body2, elements[2])
push!(body3, elements[3]);
Surface traction to the top of bodies 2 and 3:
In [15]:
t2 = Seg2([8, 7])
t2["geometry"] = Vector{Float64}[nodes[8], nodes[7]]
t2["displacement traction force"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]
t3 = Seg2([12, 11])
t3["geometry"] = Vector{Float64}[nodes[12], nodes[11]]
t3["displacement traction force"] = Vector{Float64}[[0.0, -100.0], [0.0, -100.0]]
push!(body2, t2)
push!(body3, t3);
Boundary conditions: $x=0$ for left boundary.
In [16]:
dx1 = Seg2([1, 4])
dx1["geometry"] = Vector[nodes[1], nodes[4]]
dx1["displacement 1"] = 0.0
dx2 = Seg2([5, 8])
dx2["geometry"] = Vector[nodes[5], nodes[8]]
dx2["displacement 1"] = 0.0
bc1 = DirichletProblem("displacement", 2)
push!(bc1, dx1)
push!(bc1, dx2);
$y=0$ for bottom of model
In [17]:
dy1 = Seg2([1, 2])
dy1["geometry"] = Vector[nodes[1], nodes[2]]
dy1["displacement 2"] = 0.0
bc2 = DirichletProblem("displacement", 2)
push!(bc2, dy1);
Mortar boundary conditions: tie contact between body 1 and body 2
In [19]:
#rotation_matrix(phi) = [cos(phi) -sin(phi); sin(phi) cos(phi)]
using JuliaFEM.Core: calculate_normal_tangential_coordinates!
master1 = Seg2([4, 3])
master1["geometry"] = Vector[nodes[4], nodes[3]]
slave1 = Seg2([5, 6])
slave1["geometry"] = Vector[nodes[5], nodes[6]]
slave1["master elements"] = Element[master1]
calculate_normal_tangential_coordinates!(slave1, 0.0)
#slave1["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
contact1 = MortarProblem("displacement", 2)
push!(contact1, slave1);
Tie contact between body 1 and body 3
In [20]:
slave2 = Seg2([9, 10])
slave2["geometry"] = Vector[nodes[9], nodes[10]]
#slave2["nodal ntsys"] = Matrix[rotation_matrix(-pi/2), rotation_matrix(-pi/2)]
slave2["master elements"] = Element[master1]
calculate_normal_tangential_coordinates!(slave2, 0.0)
contact2 = MortarProblem("displacement", 2)
push!(contact2, slave2);
Tie contact between body 2 and body 3
In [21]:
master2 = Seg2([6, 7])
master2["geometry"] = Vector[nodes[6], nodes[7]]
slave3 = Seg2([9, 12])
slave3["geometry"] = Vector[nodes[9], nodes[12]]
#slave3["nodal ntsys"] = Matrix[rotation_matrix(0.0), rotation_matrix(0.0)]
calculate_normal_tangential_coordinates!(slave3, 0.0)
slave3["master elements"] = Element[master2]
contact3 = MortarProblem("displacement", 2)
push!(contact3, slave3);
All defined. Solve it.
In [22]:
solver = DirectSolver()
push!(solver, body1)
push!(solver, body2)
push!(solver, body3)
push!(solver, bc1)
push!(solver, bc2)
push!(solver, contact1)
push!(solver, contact2)
push!(solver, contact3);
In [23]:
iterations, converged = call(solver, 0.0)
Out[23]:
In [24]:
using JuliaFEM.Test
@test converged
X = elements[2]("geometry", [1.0, 1.0], 0.0)
u = elements[2]("displacement", [1.0, 1.0], 0.0)
info("displacement at $X = $u")
@test isapprox(u, [0.0634862, -0.277183], atol=1.0e-5)
Out[24]:
In [1]:
using JuliaFEM.Preprocess: parse_aster_med_file
using JuliaFEM.Core: PlaneStressLinearElasticityProblem, DirichletProblem,
get_connectivity, Quad4, Tri3, Seg2, LinearSolver,
update!, get_elements
In [2]:
using JuliaFEM.Core: Element, MortarProblem, calculate_normal_tangential_coordinates!
phi = 0
rmat = [cos(phi) -sin(phi); sin(phi) cos(phi)]
function create_problems()
mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*"/geometry/2d_beam/BEAM.med")
for (k, v) in mesh["nodes"]
mesh["nodes"][k] = rmat*v
end
field_problem = PlaneStressLinearElasticityProblem()
# field problems
joo = Dict(:QU4 => Quad4, :TR3 => Tri3)
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype in keys(joo) || continue
element = joo[eltype](elcon)
update!(element, "geometry", mesh["nodes"])
element["youngs modulus"] = 900.0
element["poissons ratio"] = 0.25
push!(field_problem, element)
end
# neumann boundary condition -1 on y direction
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == :LOAD || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
# FIXME
# element["displacement traction force"] = rmat*[0.0, -0.01]
f = rmat*[0.0, -0.01]
element["displacement traction force"] = Vector{Float64}[f, f]
push!(field_problem, element)
end
# boundary conditions
boundary_problem = DirichletProblem("displacement", 2)
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == :LEFT || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
# FIXME
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
push!(boundary_problem, element)
end
info("created $(length(get_elements(field_problem))) field elements.")
info("created $(length(get_elements(boundary_problem))) boundary elements.")
# Contact definition: contact pair is `LOWER_TO_UPPER <--> UPPER_TO_LOWER`:
slave_surface = :LOWER_TO_UPPER
master_surface = :UPPER_TO_LOWER
contact_problem = MortarProblem("displacement", 2)
master_elements = JuliaFEM.Core.Element[]
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == master_surface || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
push!(master_elements, element)
push!(contact_problem, element)
end
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == slave_surface || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
element["master elements"] = master_elements
calculate_normal_tangential_coordinates!(element, 0.0)
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
push!(contact_problem, element)
end
info("# of master elements: $(length(master_elements))")
info("# of slave elements: $(length(contact_problem.elements))")
return field_problem, boundary_problem, contact_problem
end
Out[2]:
In [3]:
using JuliaFEM.Core: DirectSolver
field_problem, boundary_problem, contact_problem = create_problems()
solver = DirectSolver()
solver.name = "divided_beam"
solver.method = :UMFPACK
solver.max_iterations = 2
solver.dump_matrices = false
push!(solver, field_problem)
push!(solver, boundary_problem)
push!(solver, contact_problem)
call(solver, 0.0)
Out[3]:
In [4]:
import PyPlot
In [5]:
function plot(field_problem, scaling_factor=30, time=0.0)
fig = PyPlot.figure(figsize=(15, 4))
g_tip = rmat*[100.0, 0.0]
u_tip = nothing
for element in field_problem.elements
conn = get_connectivity(element)
X = element("geometry", time)
u = element("displacement", time)
x = X + scaling_factor*u
if isapprox(element("geometry", [1.0, -1.0], time), g_tip)
u_tip = element("displacement", [1.0, -1.0], time)
info("displacement at tip: $u_tip")
end
# undeformed
for i=1:length(X)
px1 = X[i][1]
py1 = X[i][2]
px2 = X[mod(i,length(X))+1][1]
py2 = X[mod(i,length(X))+1][2]
PyPlot.plot([px1, px2], [py1, py2], "k-", alpha=0.5)
end
# deformed
for i=1:length(x)
px1 = x[i][1]
py1 = x[i][2]
px2 = x[mod(i,length(x))+1][1]
py2 = x[mod(i,length(x))+1][2]
PyPlot.plot([px1, px2], [py1, py2], "r--", alpha=0.5)
end
end
PyPlot.axis("equal")
#PyPlot.axis("off")
return u_tip
end
u_tip = plot(field_problem, 30.0)
#PyPlot.xlim(95, 105)
#PyPlot.ylim(-10, 0)
PyPlot.grid()
In [6]:
isapprox(u_tip, [-0.025032650050963334,-0.19477533256579582]) || warn("different result")
Out[6]:
In [7]:
isapprox(norm(u_tip), 0.19637735038616433)
Out[7]:
Modification to the above: allow beams to slide in tangential direction without a friction. Global assembly operator can be modified by writing own functions preprocess_assembly! and postprocess_assembly!. Here we simply remove all kinematic constraints in tangential direction to achieve frictionless sliding between bodies.
In [8]:
using JuliaFEM.Core: SparseMatrixCOO, Element, get_integration_points, get_jacobian, get_connectivity, add!
using JuliaFEM.Core: BoundaryAssembly, BoundaryProblem, get_elements
function calculate_normal_tangential_coordinates(elements::Vector{Element}, time::Real)
P = SparseMatrixCOO()
field_dim = 2
for element in elements
haskey(element, "normal-tangential coordinates") || continue
for ip in get_integration_points(element, Val{2})
J = get_jacobian(element, ip, time)
w = ip.weight*norm(J)
nt = transpose(element("normal-tangential coordinates", ip, time))
normal = nt[1,:]
tangent = nt[2,:]
for nid in get_connectivity(element)
ndofs = [2*(nid-1)+1, 2*(nid-1)+2]
add!(P, [2*(nid-1)+1], ndofs, normal)
add!(P, [2*(nid-1)+2], ndofs, tangent)
end
end
end
P = sparse(P)
for i=1:size(P,1)
n = norm(P[i,:])
if n > 0.0
P[i,:] = P[i,:] / n
end
end
return P
end
import JuliaFEM.Core: postprocess_assembly!
function JuliaFEM.Core.postprocess_assembly!(
assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},
time::Real)
info("postprocess mortar assembly: remove contraints in tangent direction on boundary.")
C1 = sparse(assembly.C1)
C2 = sparse(assembly.C2)
dim = size(C1, 1)
P = calculate_normal_tangential_coordinates(get_elements(problem), time)
C1 = P*C1
C2 = P*C2
for i=2:2:dim
C1[i,:] = 0
C2[i,:] = 0
end
assembly.C1 = C1
assembly.C2 = C2
info("postprocess mortar assembly: done.")
end
field_problem, boundary_problem, contact_problem = create_problems()
using JuliaFEM.Core: DirectSolver
solver = DirectSolver()
solver.name = "divided_beam_small_sliding_contact"
solver.method = :UMFPACK
solver.max_iterations = 10
solver.dump_matrices = false
push!(solver, field_problem)
push!(solver, boundary_problem)
push!(solver, contact_problem)
call(solver, 0.0)
u_tip = plot(field_problem, 30)
isapprox(u_tip, [-0.03914822922685343,-0.592340359574268]) || warn("result changed")
Out[8]:
In [76]:
function create_problems_2()
mesh = parse_aster_med_file(Pkg.dir("JuliaFEM")*"/geometry/2d_beam/BEAM.med")
upper_body_nodes = Set()
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :TR3 || continue
push!(upper_body_nodes, elcon...)
end
for (nid, coords) in mesh["nodes"]
nid in upper_body_nodes || continue
mesh["nodes"][nid] = mesh["nodes"][nid] + [0.0, 1.0]
end
field_problem = PlaneStressLinearElasticityProblem()
# field problems
joo = Dict(:QU4 => Quad4, :TR3 => Tri3)
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype in keys(joo) || continue
element = joo[eltype](elcon)
update!(element, "geometry", mesh["nodes"])
element["youngs modulus"] = 900.0
element["poissons ratio"] = 0.25
push!(field_problem, element)
end
# neumann boundary condition -1 on y direction
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == :LOAD || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
# FIXME
# element["displacement traction force"] = rmat*[0.0, -0.01]
f = rmat*[0.0, -0.0125]*1.5
element["displacement traction force"] = Vector{Float64}[f, f]
push!(field_problem, element)
end
# boundary conditions
boundary_problem = DirichletProblem("displacement", 2)
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == :LEFT || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
push!(boundary_problem, element)
end
info("created $(length(get_elements(field_problem))) field elements.")
info("created $(length(get_elements(boundary_problem))) boundary elements.")
# Contact definition: contact pair is `LOWER_TO_UPPER <--> UPPER_TO_LOWER`:
master_surface = :UPPER_TO_LOWER
slave_surface = :LOWER_TO_UPPER
contact_problem = MortarProblem("displacement", 2)
master_elements = JuliaFEM.Core.Element[]
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == master_surface || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
push!(master_elements, element)
push!(contact_problem, element)
end
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
eltype == :SE2 || continue
elset == slave_surface || continue
element = Seg2(elcon)
update!(element, "geometry", mesh["nodes"])
element["master elements"] = master_elements
element["displacement"] = (0.0 => Vector{Float64}[[0.0, 0.0], [0.0, 0.0]])
calculate_normal_tangential_coordinates!(element, 0.0)
push!(contact_problem, element)
end
info("# of master elements: $(length(master_elements))")
info("# of slave elements: $(length(contact_problem.elements))")
return field_problem, boundary_problem, contact_problem
end
using JuliaFEM.Core: calculate_nodal_vector, Element
function JuliaFEM.Core.postprocess_assembly!(
assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},
time::Real)
info("postprocess mortar assembly: peforming PDASS.")
C1 = sparse(assembly.C1)
C2 = sparse(assembly.C2)
dim = size(C1, 1)
elements = get_elements(problem)
P = -calculate_normal_tangential_coordinates(elements, time)
resize!(C1, 156, 156)
resize!(C2, 156, 156)
resize!(P, 156, 156)
la = calculate_nodal_vector("reaction force", 2, elements, time)
X = calculate_nodal_vector("geometry", 2, elements, time)
u = calculate_nodal_vector("displacement", 2, elements, time)
x = X+u
#info("x = ", x)
#info("lambda = ", la)
gh = -C1*x # weighted gap vector
info("weighted gap = ", round(gh, 2))
info("lambda = ", round(la, 2))
#=
info("size P = ", size(P))
info("size C1 = ", size(C1))
info("size C2 = ", size(C2))
info("size la = ", size(la))
info("size gh = ", size(gh))
info("displacement = ", round(u, 2))
=#
#nz = sort(unique(rowvals(P)))
#dump(round(full(P[nz,nz]), 3))
C1 = P*C1
C2 = P*C2
gh = P*gh
la = P*la
c = 1.0
# complementarity function
C = la - clamp(la - c*gh, 0, Inf)
C[abs(C) .< 1.0e-12] = 0.0
C[49] = -1 # node in support
#info("complementarity function = ", round(C, 2))
X2 = reshape(X, 2, round(Int, length(X)/2))
j = 0
for i=1:2:dim
j += 1
X2[1,j] == 0 && continue
info("dof $i: normal gap = $(round(gh[i], 3)), X=$(round(X2[:,j], 2)), C=$(round(C[i], 5)), la=$(round(la[i], 5))")
#if C[i] > 0
if i == 45
info("contact dof $i active (C = $(C[i]))")
#gh[i] = -gh[i]
else
C1[i, :] = 0
C2[i, :] = 0
gh[i] = 0
end
end
for i=2:2:dim
C1[i, :] = 0
C2[i, :] = 0
gh[i] = 0
end
#assembly.g = g
assembly.C1 = C1
assembly.C2 = C2
assembly.g = sparse(gh)
info("postprocess mortar assembly: done.")
end
field_problem, boundary_problem, contact_problem = create_problems_2()
using JuliaFEM.Core: DirectSolver
solver = DirectSolver()
solver.name = "divided_beam_inequality_constraints"
solver.method = :UMFPACK
solver.max_iterations = 5
solver.dump_matrices = false
push!(solver, field_problem)
push!(solver, boundary_problem)
push!(solver, contact_problem)
call(solver, 0.0)
u_tip = plot(field_problem, 1)
Out[76]:
In [17]:
A = sparse(rand(3, 3))
resize!(A, 4, 4)
full(A)
Out[17]:
In [34]:
contact_problem.elements[5]("displacement", [0.0], 0.0)
Out[34]:
In [38]:
JuliaFEM.Core.get_gdofs(contact_problem.elements[5], 2)
Out[38]:
In [10]:
function JuliaFEM.Core.postprocess_assembly!(
assembly::BoundaryAssembly, problem::BoundaryProblem{MortarProblem},
time::Real)
info("postprocess mortar assembly: remove contraints in tangent direction on boundary.")
C1 = sparse(assembly.C1)
C2 = sparse(assembly.C2)
g = sparse(assembly.g)
n = round(Int, length(g)/2)
dump(reshape(full(g), 2, n)[:,60:end]')
dim = size(C1, 1)
P = calculate_normal_tangential_coordinates(get_elements(problem), time)
C1 = P*C1
C2 = P*C2
for i=2:2:dim
C1[i,:] = 0
C2[i,:] = 0
end
assembly.g = g
assembly.C1 = C1
assembly.C2 = C2
info("postprocess mortar assembly: done.")
end
el1 = Seg2([1, 2])
el1["geometry"] = Vector{Float64}[[0.0, 0.0], [4.0, 0.0]]
el2 = Seg2([3, 4])
el2["geometry"] = Vector{Float64}[[0.0, 1.0], [4.0, 1.0]]
el1["master elements"] = Element[el2]
JuliaFEM.Core.calculate_normal_tangential_coordinates!(el1, 0.0)
p = MortarProblem("displacement", 2)
push!(p, el1)
a = JuliaFEM.Core.assemble(p, 0.0)
full(a.g)
Out[10]:
In [1]:
using JuliaFEM.Core: Seg2, Tri3, Quad4, Node, update!,
calculate_normal_tangential_coordinates!, PlaneStressLinearElasticityProblem,
ContactProblem, SmallSlidingContact, Element
using JuliaFEM.Core: PlaneStressLinearElasticityProblem, DirichletProblem,
get_connectivity, Quad4, Tri3, Seg2, LinearSolver,
update!, get_elements, BiorthogonalBasis, StandardBasis
nodes = Dict{Int64, Node}(
1 => [0.0, 2.0],
2 => [2.0, 6.0],
3 => [2.0, 2.0],
4 => [2.0, 4.0],
5 => [0.0, 0.0],
6 => [2.0, 0.0],
7 => [2.0, 0.0],
8 => [4.0, 0.0])
Out[1]:
In [2]:
basis = StandardBasis
el1 = Quad4([5, 6, 3, 1])
el2 = Tri3([1, 3, 4])
el3 = Tri3([7, 8, 2])
del1 = Seg2([5, 6])
del2 = Seg2([7, 8])
sel1 = Seg2([4, 3])
sel2 = Seg2([3, 6])
mel1 = Seg2([2, 7])
update!([el1, el2, el3, sel1, sel2, mel1, del1, del2], "geometry", nodes)
for el in [el1, el2, el3]
el["youngs modulus"] = 90.0
el["poissons ratio"] = 0.25
end
for el in [del1, del2]
el["displacement"] = 0.0
end
for el in [sel1, sel2]
el["master elements"] = Element[mel1]
calculate_normal_tangential_coordinates!(el, 0.0)
end
prob = PlaneStressLinearElasticityProblem()
push!(prob, el1, el2, el3)
bc = DirichletProblem("displacement", 2; basis=basis)
push!(bc, del1, del2)
con = ContactProblem("cont", "displacement", 2;
contact_type=SmallSlidingContact)
push!(con, sel1, sel2);
In [3]:
ass1 = assemble(prob, 0.0)
dbc = assemble(bc, 0.0)
cbc = assemble(con, 0.0)
Out[3]:
In [4]:
K = full(ass1.stiffness_matrix)
dims = size(K)
K[1:8,1:8]
Out[4]:
In [6]:
fixed = [9, 10, 11, 12, 13, 14, 15, 16]
ENV["COLUMNS"] = 300
C1 = full(cbc.C1, dims...)
C1[abs(C1) .< 1.0e-9] = 0
#C1[fixed, fixed] = 0
C1*18
Out[6]:
In [7]:
D = full(cbc.D, dims...)
D[abs(D) .< 1.0e-9] = 0
D[fixed, fixed] = 0
D
Out[7]:
In [8]:
C2 = full(cbc.C2, dims...)
C2[abs(C2) .< 1.0e-9] = 0
#C2[fixed, :] = 0
for i=2:2:16
C1[i,:] = 0
C2[i,:] = 0
end
C2*18
Out[8]:
In [9]:
D
Out[9]:
In [10]:
DC1 = full(dbc.C1, dims...)
DC2 = full(dbc.C2, dims...)
DD = full(dbc.D, dims...)
DC1*18
Out[10]:
In [11]:
f = zeros(16)
f[1] = 18.0
g = zeros(16)
A = [K (C1+DC1)'; (C2+DC2) (D+DD)]
b = [f; g]
Out[11]:
In [12]:
for i=1:size(A,1)
if sum(abs(A[i,:])) == 0
info("zero line $i")
A[i,i] = 1.0
end
end
In [13]:
u = A \ b
u = reshape(u, 2, 16)
u[abs(u) .< 1.0e-9] = 0
disp = u[:,1:8]
Out[13]:
In [24]:
# tie contact
@assert isapprox(norm(vec(disp)), 0.8223229928825583)
In [19]:
# frictionless sliding
@assert isapprox(norm(vec(disp)), 0.9363129098080032)
In [15]:
X = [0.0 2; 2 6; 2 2; 2 4; 0 0; 2 0; 2 0; 4 0]'
Out[15]:
In [16]:
x = X + disp
Out[16]:
In [17]:
eldofs = Dict()
eldofs[1] = [5, 6, 3, 1]
eldofs[2] = [1, 3, 4]
eldofs[3] = [7, 8, 2]
Out[17]:
In [18]:
using PyPlot
function get_coords(X)
return [X X[:, 1]]
end
function plot(X, args...; kwargs...)
X1 = X[:, eldofs[1]]
X2 = X[:, eldofs[2]]
X3 = X[:, eldofs[3]]
PyPlot.plot(get_coords(X1)[1,:]', get_coords(X1)[2,:]', args...; kwargs...)
PyPlot.plot(get_coords(X2)[1,:]', get_coords(X2)[2,:]', args...; kwargs...)
PyPlot.plot(get_coords(X3)[1,:]', get_coords(X3)[2,:]', args...; kwargs...)
axis("off")
end
figure(figsize=(4, 5))
plot(X, "-k")
plot(x, "--r")
xlim(-0.1, 4.1)
ylim(-0.1, 6.1)
Out[18]:
In [ ]: