mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
frictionless 2d contact
This commit is contained in:
Jukka Aho
Binary file not shown.
|
After Width: | Height: | Size: 6.1 KiB |
Binary file not shown.
Binary file not shown.
File diff suppressed because one or more lines are too long
+77
-6
@@ -1,16 +1,18 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract DirichletProblem <: AbstractProblem
|
||||
abstract DirichletProblem{T} <: AbstractProblem
|
||||
abstract StandardBasis
|
||||
abstract BiorthogonalBasis
|
||||
|
||||
function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
|
||||
return BoundaryProblem{DirichletProblem}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements)
|
||||
function DirichletProblem(parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]; basis=StandardBasis)
|
||||
return BoundaryProblem{DirichletProblem{basis}}("dirichlet boundary", parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[])
|
||||
return BoundaryProblem{DirichletProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
|
||||
function DirichletProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=Element[]; basis=StandardBasis)
|
||||
return BoundaryProblem{DirichletProblem{basis}}(problem_name, parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
|
||||
function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem}, element::Element, time::Real)
|
||||
function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem{StandardBasis}}, element::Element, time::Real)
|
||||
|
||||
# get dimension and name of PARENT field
|
||||
field_dim = problem.parent_field_dim
|
||||
@@ -52,3 +54,72 @@ function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{Dirichle
|
||||
end
|
||||
end
|
||||
|
||||
function assemble!(assembly::BoundaryAssembly, problem::BoundaryProblem{DirichletProblem{BiorthogonalBasis}}, element::Element, time::Real)
|
||||
|
||||
# get dimension and name of PARENT field
|
||||
field_dim = problem.parent_field_dim
|
||||
field_name = problem.parent_field_name
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
|
||||
# calculate bi-orthogonal basis transformation matrix Ae
|
||||
nnodes = size(element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(element, Val{2})
|
||||
w = ip.weight
|
||||
J = get_jacobian(element, ip, time)
|
||||
JT = transpose(J)
|
||||
if size(JT, 2) == 1 # plane problem
|
||||
w *= norm(JT)
|
||||
else
|
||||
w *= norm(cross(JT[:,1], JT[:,2]))
|
||||
end
|
||||
N = element(ip, time)
|
||||
De += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
# do the actual integration
|
||||
for ip in get_integration_points(element, Val{2})
|
||||
w = ip.weight
|
||||
J = get_jacobian(element, ip, time)
|
||||
JT = transpose(J)
|
||||
if size(JT, 2) == 1 # plane problem
|
||||
w *= norm(JT)
|
||||
else
|
||||
w *= norm(cross(JT[:,1], JT[:,2]))
|
||||
end
|
||||
N = element(ip, time)
|
||||
Phi = (Ae*N')'
|
||||
A = w*Phi'*N
|
||||
A[abs(A) .< 1.0e-9] = 0
|
||||
|
||||
# C1 matrix is always the same
|
||||
for i=1:field_dim
|
||||
ldofs = gdofs[i:field_dim:end]
|
||||
add!(assembly.C1, ldofs, ldofs, A)
|
||||
end
|
||||
|
||||
if haskey(element, field_name)
|
||||
# add all dimensions at once if defined element["blaa"] = 0.0
|
||||
for i=1:field_dim
|
||||
g = element(field_name, ip, time)
|
||||
ldofs = gdofs[i:field_dim:end]
|
||||
add!(assembly.C2, ldofs, ldofs, A)
|
||||
add!(assembly.g, ldofs, w*g*Phi')
|
||||
end
|
||||
else
|
||||
for i=1:field_dim
|
||||
ldofs = gdofs[i:field_dim:end]
|
||||
if haskey(element, field_name*" $i")
|
||||
g = element(field_name*" $i", ip, time)
|
||||
add!(assembly.C2, ldofs, ldofs, A)
|
||||
add!(assembly.g, ldofs, w*g*Phi')
|
||||
else
|
||||
add!(assembly.D, ldofs, ldofs, A)
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -12,6 +12,10 @@ function Base.size{E}(::Element{E})
|
||||
return size(E)
|
||||
end
|
||||
|
||||
function Base.size{E}(::Element{E}, i::Int64)
|
||||
return size(E)[i]
|
||||
end
|
||||
|
||||
function convert{E}(::Type{Element{E}}, connectivity::Vector{Int})
|
||||
# return Element{E}(connectivity, get_integration_points(E), Dict())
|
||||
return Element{E}(connectivity, Dict())
|
||||
|
||||
+131
-1
@@ -658,7 +658,21 @@ function MortarProblem(problem_name::ASCIIString, parent_field_name::ASCIIString
|
||||
return BoundaryProblem{MortarProblem}(problem_name, parent_field_name, parent_field_dim, dim, elements)
|
||||
end
|
||||
|
||||
# Mortar assembly
|
||||
abstract ContactProblem{T} <: AbstractProblem
|
||||
|
||||
abstract AbstractContact
|
||||
abstract TieContact <: AbstractContact
|
||||
abstract SmallSlidingContact <: AbstractContact
|
||||
|
||||
function ContactProblem(problem_name::ASCIIString, parent_field_name::ASCIIString, parent_field_dim::Int, dim::Int=1, elements=[]; contact_type=TieContact)
|
||||
return BoundaryProblem{ContactProblem{contact_type}}(
|
||||
problem_name,
|
||||
parent_field_name,
|
||||
parent_field_dim,
|
||||
dim, elements)
|
||||
end
|
||||
|
||||
# Mortar assembly 2d
|
||||
|
||||
typealias MortarElements2D Union{Seg2, Seg3}
|
||||
|
||||
@@ -707,6 +721,122 @@ function assemble!{E<:MortarElements2D}(assembly::BoundaryAssembly, problem::Bou
|
||||
end
|
||||
end
|
||||
|
||||
""" Calculate bi-orthogonal basis transformation matrix Aₑ. """
|
||||
function get_biorthogonal_transformation_matrix(element::Element, time::Real)
|
||||
nnodes = size(element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(element, Val{5})
|
||||
w = ip.weight
|
||||
J = get_jacobian(element, ip, time)
|
||||
JT = transpose(J)
|
||||
if size(JT, 2) == 1 # plane problem
|
||||
w *= norm(JT)
|
||||
else
|
||||
w *= norm(cross(JT[:,1], JT[:,2]))
|
||||
end
|
||||
N = element(ip, time)
|
||||
De += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
return Ae
|
||||
end
|
||||
|
||||
"""
|
||||
Small strain theory, allow frictionless tangential sliding, keep bodies in contact.
|
||||
"""
|
||||
function assemble!{E<:MortarElements2D}(assembly::BoundaryAssembly, problem::BoundaryProblem{ContactProblem{SmallSlidingContact}}, slave_element::Element{E}, time::Real)
|
||||
|
||||
# get dimension and name of PARENT field
|
||||
field_dim = problem.parent_field_dim
|
||||
field_name = problem.parent_field_name
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
info("slave dofs of element: $slave_dofs")
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
xi1a = project_from_master_to_slave(slave_element, master_element, [-1.0])
|
||||
xi1b = project_from_master_to_slave(slave_element, master_element, [ 1.0])
|
||||
xi1 = clamp([xi1a xi1b], -1.0, 1.0)
|
||||
l = 1/2*(xi1[2]-xi1[1])
|
||||
abs(l) > 1.0e-9 || continue
|
||||
|
||||
# Ae = get_biorthogonal_transformation_matrix(slave_element, time)
|
||||
nnodes = size(element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip_ in get_integration_points(slave_element, Val{5})
|
||||
xi_gauss = 1/2*(1-ip_.xi)*xi1[1] + 1/2*(1+ip_.xi)*xi1[2]
|
||||
ip = IntegrationPoint(xi_gauss, ip_.weight)
|
||||
w = ip.weight
|
||||
J = get_jacobian(slave_element, ip, time)
|
||||
JT = transpose(J)
|
||||
if size(JT, 2) == 1 # plane problem
|
||||
w *= norm(JT)
|
||||
else
|
||||
w *= norm(cross(JT[:,1], JT[:,2]))
|
||||
end
|
||||
N = element(ip, time)
|
||||
De += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
J = get_jacobian(slave_element, ip, time)
|
||||
w = ip.weight*norm(J)*l
|
||||
|
||||
# integration point on slave side segment
|
||||
xi_gauss = 1/2*(1-ip.xi)*xi1[1] + 1/2*(1+ip.xi)*xi1[2]
|
||||
# projected integration point
|
||||
xi_projected = project_from_slave_to_master(slave_element, master_element, xi_gauss)
|
||||
|
||||
# add contribution to C1
|
||||
N1 = slave_element(xi_gauss, time)
|
||||
Phi = (Ae*N1')'
|
||||
N2 = master_element(xi_projected, time)
|
||||
S = w*Phi'*N1
|
||||
M = w*Phi'*N2
|
||||
for i=1:field_dim
|
||||
sd = slave_dofs[i:field_dim:end]
|
||||
md = master_dofs[i:field_dim:end]
|
||||
add!(assembly.C1, sd, sd, S)
|
||||
add!(assembly.C1, sd, md, -M)
|
||||
end
|
||||
|
||||
# construct C2 & D
|
||||
nt = slave_element("normal-tangential coordinates", ip, time)
|
||||
nt = transpose(nt)
|
||||
ntS = nt*S
|
||||
ntM = nt*M
|
||||
info("normal dofs: $(slave_dofs[1:2:end])")
|
||||
info("tangent dofs: $(slave_dofs[2:2:end])")
|
||||
# contribution in normal direction
|
||||
for dof in slave_dofs[1:2:end]
|
||||
add!(assembly.C2, [dofs], sd, ntS[1,:])
|
||||
add!(assembly.C2, sd[1:2:end], md, -ntM[1,:])
|
||||
end
|
||||
# contribution in tangent direction
|
||||
add!(assembly.C2, sd[2:2:end], sd, ntS[2,:])
|
||||
add!(assembly.C2, sd[2:2:end], md, -ntM[2,:])
|
||||
# set lagrange multipliers to zero in tangent direction
|
||||
#tangent = nt[2, :]
|
||||
#add!(assembly.D, sd[2:2:end], sd, tangent)
|
||||
end
|
||||
#=
|
||||
nt = transpose(nt)
|
||||
normal = nt[1,:]
|
||||
tangent = nt[2,:]
|
||||
for nid in get_connectivity(slave_element)
|
||||
ndofs = [2*(nid-1)+1, 2*(nid-1)+2]
|
||||
add!(assembly.C2, [2*(nid-1)+1], ndofs, normal)
|
||||
add!(assembly.D, [2*(nid-1)+2], ndofs, tangent)
|
||||
end
|
||||
=#
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
typealias MortarElements3D Union{Tri3, Quad4}
|
||||
|
||||
|
||||
@@ -1,6 +1,8 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
typealias Node Vector{Float64}
|
||||
|
||||
function ForwardDiff.derivative{T}(f::Function, S::Matrix{T}, args...)
|
||||
shape = size(S)
|
||||
wrapper(S::Vector) = f(reshape(S, shape))
|
||||
|
||||
+69
-27
@@ -4,64 +4,106 @@
|
||||
module TestDirichletBoundaryCondition
|
||||
|
||||
using JuliaFEM.Test
|
||||
using JuliaFEM.Core: Tri3, Seg2, DirichletProblem, Assembly, assemble
|
||||
using JuliaFEM.Core: Tri3, Seg2, DirichletProblem, Assembly, assemble, Node,
|
||||
BiorthogonalBasis
|
||||
|
||||
function test_dirichlet_problem_1_dim()
|
||||
@testset "test dirichlet boundary conditions" begin
|
||||
|
||||
@testset "dirichlet problem in 1 dimension" begin
|
||||
element = Seg2([1, 2])
|
||||
element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["temperature"] = 0.0
|
||||
problem = DirichletProblem("temperature", 1)
|
||||
push!(problem, element)
|
||||
assembly = assemble(problem, 0.0)
|
||||
A = full(assembly.stiffness_matrix)
|
||||
b = full(assembly.force_vector)
|
||||
@test isapprox(A, 1/6*[2 1; 1 2])
|
||||
@test isapprox(b, [0.0, 0.0])
|
||||
C1 = full(assembly.C1)
|
||||
C2 = full(assembly.C2)
|
||||
g = full(assembly.g)
|
||||
@test isapprox(C1, C2)
|
||||
@test isapprox(C1, 1/6*[2 1; 1 2])
|
||||
@test isapprox(g, [0.0, 0.0])
|
||||
end
|
||||
|
||||
function test_dirichlet_problem_2_dim()
|
||||
@testset "dirichlet problem in 2 dimensions" begin
|
||||
element = Seg2([1, 2])
|
||||
element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["displacement"] = 0.0
|
||||
problem = DirichletProblem("displacement", 2)
|
||||
push!(problem, element)
|
||||
assembly = assemble(problem, 0.0)
|
||||
A = full(assembly.stiffness_matrix)
|
||||
b = full(assembly.force_vector)
|
||||
A_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
|
||||
@test isapprox(A, A_expected)
|
||||
@test isapprox(b, [0.0, 0.0, 0.0, 0.0])
|
||||
C1 = full(assembly.C1)
|
||||
C2 = full(assembly.C2)
|
||||
g = full(assembly.g)
|
||||
@test isapprox(C1, C2)
|
||||
C1_expected = 1/6*[2 0 1 0; 0 2 0 1; 1 0 2 0; 0 1 0 2]
|
||||
@test isapprox(C1, C1_expected)
|
||||
@test isapprox(g, [0.0, 0.0, 0.0, 0.0])
|
||||
end
|
||||
|
||||
function test_dirichlet_problem_2_dim_single_dof_fixed()
|
||||
@testset "dirichlet problem in 2 dimensions, with 1 dof fixed" begin
|
||||
element = Seg2([1, 2])
|
||||
element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["geometry"] = Node[[1.0, 1.0], [0.0, 1.0]]
|
||||
element["displacement 2"] = 0.0
|
||||
problem = DirichletProblem("displacement", 2)
|
||||
push!(problem, element)
|
||||
assembly = assemble(problem, 0.0)
|
||||
A = full(assembly.stiffness_matrix)
|
||||
b = full(assembly.force_vector)
|
||||
info(b)
|
||||
info("A = \n$A")
|
||||
A_expected = 1/6*[
|
||||
C1 = full(assembly.C1)
|
||||
C2 = full(assembly.C2)
|
||||
g = full(assembly.g)
|
||||
@test isapprox(C1, C2)
|
||||
C1_expected = 1/6*[
|
||||
0 0 0 0
|
||||
0 2 0 1
|
||||
0 0 0 0
|
||||
0 1 0 2]
|
||||
@test isapprox(A, A_expected)
|
||||
@test isapprox(b, [0.0, 0.0, 0.0, 0.0])
|
||||
@test isapprox(C1, C1_expected)
|
||||
@test isapprox(g, [0.0, 0.0, 0.0, 0.0])
|
||||
end
|
||||
|
||||
function test_dirichlet_surface_tri3()
|
||||
@testset "dirichlet problem using tri3 surface element" begin
|
||||
elem = Tri3([1, 2, 3])
|
||||
elem["geometry"] = Vector{Float64}[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]
|
||||
elem["geometry"] = Node[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
|
||||
elem["temperature"] = 0.0
|
||||
prob = DirichletProblem("temperature", 1)
|
||||
push!(prob, elem)
|
||||
ass = assemble(prob, 0.0)
|
||||
k = full(ass.stiffness_matrix)
|
||||
@test isapprox(k, 1/24*[2 1 1; 1 2 1; 1 1 2])
|
||||
C1 = full(ass.C1)
|
||||
C2 = full(ass.C2)
|
||||
@test isapprox(C1, C2)
|
||||
@test isapprox(C1, 1/24*[2 1 1; 1 2 1; 1 1 2])
|
||||
end
|
||||
|
||||
@testset "dirichlet problem using biorthogonal basis" begin
|
||||
elem = Tri3([1, 2, 3])
|
||||
elem["geometry"] = Node[[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
|
||||
elem["displacement 1"] = 1.0
|
||||
prob = DirichletProblem("displacement", 3; basis=BiorthogonalBasis)
|
||||
push!(prob, elem)
|
||||
ass = assemble(prob, 0.0)
|
||||
C1 = full(ass.C1, 9, 9)
|
||||
C2 = full(ass.C2, 9, 9)
|
||||
D = full(ass.D, 9, 9)
|
||||
g = full(ass.g, 9, 1)
|
||||
C1_expected = eye(9)*1.0/6.0
|
||||
g_expected = zeros(9)
|
||||
g_expected[1] = g_expected[4] = g_expected[7] = 1.0/6.0
|
||||
C2_expected = zeros(9, 9)
|
||||
D_expected = zeros(9, 9)
|
||||
C2_expected[1, 1] = 1.0/6.0
|
||||
D_expected[2, 2] = 1.0/6.0
|
||||
D_expected[3, 3] = 1.0/6.0
|
||||
C2_expected[4, 4] = 1.0/6.0
|
||||
D_expected[5, 5] = 1.0/6.0
|
||||
D_expected[6, 6] = 1.0/6.0
|
||||
C2_expected[7, 7] = 1.0/6.0
|
||||
D_expected[8, 8] = 1.0/6.0
|
||||
D_expected[9, 9] = 1.0/6.0
|
||||
@test isapprox(C1, C1_expected)
|
||||
@test isapprox(C2, C2_expected)
|
||||
@test isapprox(D, D_expected)
|
||||
@test isapprox(g, g_expected)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user