mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-29 07:32:54 +00:00
Mortar autodiff tests...
This commit is contained in:
Jukka Aho
+3
-421
@@ -82,416 +82,9 @@ function project_from_slave_to_master{E<:MortarElements2D}(
|
||||
end
|
||||
|
||||
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
|
||||
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real,
|
||||
::Type{Val{:forwarddiff_old}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
nnodes = size(slave_element, 2)
|
||||
X1 = slave_element("geometry", time)
|
||||
#u1 = slave_element("displacement", time)
|
||||
#x1 = X1 + u1
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
|
||||
X = spzeros(10000, 1)
|
||||
for element in get_elements(problem)
|
||||
conn = get_connectivity(element)
|
||||
geom = element("geometry", time)
|
||||
for (c, g) in zip(conn, geom)
|
||||
dofs = [field_dim*(c-1)+1, field_dim*(c-1)+2]
|
||||
X[dofs] = g
|
||||
end
|
||||
end
|
||||
|
||||
# here x does not mean deformed configuration
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
if length(x) == 0
|
||||
info("mortar_2d_autodiff: length(x) == 0")
|
||||
# resize solution vectors according to initial configuration of this problem
|
||||
X = vec(full(sparse(findnz(X)...)))
|
||||
x = zeros(length(X)*2)
|
||||
else
|
||||
# resize initial configuration to match real dimension
|
||||
I, J, V = findnz(X)
|
||||
X = vec(full(sparse(I, J, V, length(problem.assembly.u), 1)))
|
||||
end
|
||||
ndofs = round(Int, length(x)/2)
|
||||
# at the end we should have
|
||||
# info("mortar_2d_autodiff: size of x = $(size(x))")
|
||||
# info("mortar_2d_autodiff: size of X = $(size(X))")
|
||||
# info("mortar_2d_autodiff: ndofs = $ndofs")
|
||||
|
||||
""" Calculate normal vector for slave element nodes in current configuration. """
|
||||
function calculate_normals(u::Matrix)
|
||||
normals = zeros(u)
|
||||
# 1. update nodal normals
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
slave_normals = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
return slave_normals
|
||||
end
|
||||
|
||||
function calculate_mortar_projection(u::Matrix, n1::DVTI)
|
||||
B = SparseMatrixCOO{Real}([], [], [])
|
||||
|
||||
u1 = Field([u[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
u2 = Field([u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
#info("master element coordinate 1 = $(ForwardDiff.get_value(x2[1]))")
|
||||
#info("master element coordinate 2 = $(ForwardDiff.get_value(x2[2]))")
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution
|
||||
|
||||
# integrate slave side
|
||||
D = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
# jacobian of slave element in deformed state
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N = get_basis(slave_element, xi_s)
|
||||
D += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
Ae = D*inv(Me)
|
||||
|
||||
# integrate master side
|
||||
M = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
# jacobian of slave element in deformed state
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_g = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_g)
|
||||
x_g = vec(N1)*x1
|
||||
n_g = vec(N1)*n1
|
||||
#info("slave gauss point coordinates $(ForwardDiff.get_value(x_g))")
|
||||
#info("slave gauss point normal direction $(ForwardDiff.get_value(n_g))")
|
||||
xi_m = project_from_slave_to_master(master_element, x_g, n_g, x2)
|
||||
N2 = get_basis(master_element, xi_m)
|
||||
M += w*kron(Ae*N1', N2)
|
||||
end
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
for i=1:field_dim
|
||||
add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D)
|
||||
add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
return B
|
||||
end
|
||||
|
||||
function calculate_contact_rhs(x::Vector)
|
||||
ndofs = round(Int, length(x)/2)
|
||||
u = x[1:ndofs]
|
||||
la = x[ndofs+1:end]
|
||||
# info("calculate_contact_rhs: size of u = $(size(u))")
|
||||
# info("calculate_contact_rhs: size of X = $(size(X))")
|
||||
# info("calculate_contact_rhs: size of la = $(size(la))")
|
||||
# info("calculate_contact_rhs: ndofs = $ndofs")
|
||||
u2 = reshape(u, field_dim, round(Int, length(u)/field_dim))
|
||||
|
||||
normals = calculate_normals(u2)
|
||||
B = calculate_mortar_projection(u2, normals)
|
||||
B = sparse(B, ndofs, ndofs)
|
||||
fc = B' * la # contact force residual for r = fint + fc - fext = 0
|
||||
|
||||
N = SparseMatrixCOO{Real}([], [], [])
|
||||
T = SparseMatrixCOO{Real}([], [], [])
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
add!(N, dofs, [dofs[1]], reshape(normals[i], 2, 1))
|
||||
add!(T, dofs, [dofs[2]], reshape(Q'*normals[i], 2, 1))
|
||||
end
|
||||
N = sparse(N, ndofs, ndofs)
|
||||
T = sparse(T, ndofs, ndofs)
|
||||
gn = -N*B*(X+u)
|
||||
gt = -T*B*(X+u)
|
||||
# gn = -N*B*u
|
||||
lan = N*la
|
||||
lat = T*la
|
||||
cn = 1.0e3
|
||||
C = lan - max(0, lan - cn*gn) + lat
|
||||
# C = lan - gn + lat - gt <-- ihan viturallensa
|
||||
# C = gn+gt <-- not working
|
||||
# C = B*(X+u)
|
||||
# C = B*u <- pitää kiinni, "tie".
|
||||
# C = N*B*u + T*la <- palikat menee väärään suuntaan
|
||||
# C = -N*B*(X+u) + T*la <- toimii suht hyvin mut kääntyy väärään suuntaan (t-suunnassa)
|
||||
# C = -N*B*(X+u) - T*la <- sama
|
||||
# C = -N*B*(X+u) - T*B*la <- sama
|
||||
# C = N*la + T*la - max(0, N*la + N*B*(X+u))
|
||||
cond = lan[1:field_dim:end] - cn*gn[1:field_dim:end]
|
||||
all_nodes = slave_element_nodes
|
||||
inactive_nodes = find(cond .<= 0)
|
||||
active_nodes = find(cond .> 0)
|
||||
inactive_nodes = setdiff(all_nodes, inactive_nodes)
|
||||
active_nodes = setdiff(all_nodes, active_nodes)
|
||||
info("S = $all_nodes, I = $inactive_nodes, A = $active_nodes")
|
||||
info("lambda = $(ForwardDiff.get_value(lan[slave_dofs]))")
|
||||
info("gn = $(ForwardDiff.get_value(gn[slave_dofs]))")
|
||||
#for j in active_nodes
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# C[dofs] = 0
|
||||
#end
|
||||
|
||||
# C = -(N+T)*B*(X+u)
|
||||
# C = -B*(X+u)
|
||||
|
||||
return [fc; C]
|
||||
|
||||
end
|
||||
|
||||
A, allresults = ForwardDiff.jacobian(calculate_contact_rhs, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = A[1:ndofs,ndofs+1:end]'
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = b[1:ndofs]
|
||||
g = b[ndofs+1:end]
|
||||
C2[2:field_dim:end] = 0
|
||||
g[2:field_dim:end] = 0
|
||||
|
||||
# C2[2:field_dim:end] = 0
|
||||
# g[2:field_dim:end] = 0
|
||||
# inactives = find(g[1:field_dim:end] .<= 0)
|
||||
# actives = find(g[1:field_dim:end] .> 0)
|
||||
# inactives = setdiff(slave_element_nodes, inactives)
|
||||
# actives = setdiff(slave_element_nodes, actives)
|
||||
# info("all nodes = $slave_element_nodes, inactives = $inactives, actives = $actives")
|
||||
# info("g = $g")
|
||||
|
||||
# for j in inactives
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# K[dofs,:] = 0
|
||||
# C1[dofs,:] = 0
|
||||
# C2[dofs,:] = 0
|
||||
# D[dofs,:] = 0
|
||||
# g[dofs,:] = 0
|
||||
#end
|
||||
#for j in actives
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# C2[dofs[1],:] = 0
|
||||
#end
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff_old2}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(x::Vector)
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc = zeros(u)
|
||||
C = zeros(la)
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
# 1. update nodal normals for this element. average nodes from adjacent elements
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
|
||||
# project gauss point from slave element to master element
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s))
|
||||
x_s = N1*x1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = Q'*n_s # tangent direction in gauss point
|
||||
R_s = [n_s t_s]
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = vec(get_basis(master_element, xi_m))
|
||||
x_m = N2*x2
|
||||
Phi = Ae*N1
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
|
||||
la_s = Phi*la1 # traction force in gauss point
|
||||
#lan = dot(n_s, la_s) # normal component
|
||||
#lat = dot(t_s, la_s) # tangential component
|
||||
la_nt = R_s*la_s
|
||||
g = x_s-x_m
|
||||
gn = props.gap_sign*dot(n_s, g) # normal gap
|
||||
#gu = dot(n_s, u_s - u_m) # normal displacement gap
|
||||
#gt = dot(t_s, x_s - x_m) # tangential gap
|
||||
|
||||
fc[:,slave_element_nodes] += w*la_s*N1'
|
||||
fc[:,master_element_nodes] -= w*la_s*N2'
|
||||
C[1,slave_element_nodes] += w*gn*Phi'
|
||||
#C[2,slave_element_nodes] += w*la_nt[2]*Phi'
|
||||
#C[2,slave_element_nodes] += w*la_nt[2,:]*Phi'
|
||||
#C[2,slave_element_nodes] += w*dot(t_s, u_s - u_m)*Phi' # <-- for tie
|
||||
|
||||
#R += w*R_s
|
||||
|
||||
end
|
||||
|
||||
end # master elements done
|
||||
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
n = n1[i]
|
||||
t = Q'*n
|
||||
R = [n t]
|
||||
la_nt = R*la[:,j]
|
||||
C[2,j] += la_nt[2]
|
||||
end
|
||||
|
||||
return vec([fc C])
|
||||
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
#if ndofs == 0
|
||||
# info("INITIALIZING THINGS")
|
||||
# problem.assembly.u = zeros(16)
|
||||
# problem.assembly.la = zeros(16)
|
||||
# x = [problem.assembly.u; problem.assembly.la]
|
||||
# ndofs = round(Int, length(x)/2)
|
||||
#end
|
||||
|
||||
A, allresults = ForwardDiff.jacobian(calculate_interface, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
#b = -calculate_interface(x)
|
||||
#info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
#println(A)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(A[1:ndofs,ndofs+1:end])
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = b[1:ndofs]
|
||||
g = b[ndofs+1:end]
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff}})
|
||||
@@ -508,7 +101,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc = zeros(u)
|
||||
C = zeros(la)
|
||||
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
@@ -601,7 +194,7 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
fc[:,slave_element_nodes] += w*la_s*N1'
|
||||
fc[:,master_element_nodes] -= w*la_s*N2'
|
||||
#C[1,slave_element_nodes] += w*gn*Phi'
|
||||
|
||||
|
||||
lan_tot += w*la_nt[1]*Phi
|
||||
gap_tot += w*gn*Phi
|
||||
|
||||
@@ -609,15 +202,6 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
|
||||
end # master elements done
|
||||
|
||||
# ncf = lan_tot - max(0, lan_tot - gap_tot)
|
||||
# info("pressure in nodes: $(ForwardDiff.get_value(lan_tot))")
|
||||
# info("weighted gap in nodes: $(ForwardDiff.get_value(gap_tot))")
|
||||
# info("ncf: $(ForwardDiff.get_value(ncf))")
|
||||
# cond = +lan_tot + gap_tot
|
||||
# cond = +lan_tot - gap_tot # singular
|
||||
# cond = -lan_tot + gap_tot
|
||||
# cond = -lan_tot - gap_tot
|
||||
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
n = n1[i]
|
||||
t = Q'*n
|
||||
@@ -627,9 +211,8 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
|
||||
# if -lan_tot[i] + gap_tot[i] < 0
|
||||
if -la_nt[1] + gap_tot[i] < 0
|
||||
#if j in [31, 32, 33, 34, 35, 36]
|
||||
info("set node $j active")
|
||||
C[1,j] -= gap_tot[i]
|
||||
C[1,j] += gap_tot[i]
|
||||
# C[1,j] += la_nt[1] - max(0, la_nt[1] - gap_tot[i])
|
||||
C[2,j] += la_nt[2]
|
||||
else
|
||||
@@ -671,4 +254,3 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
|
||||
+1
-1
@@ -177,7 +177,7 @@ function update_assembly!(problem, u, la)
|
||||
elseif get_formulation_type(problem) == :forwarddiff
|
||||
info("$(problem.name): forwarddiff formulation, adding increment to solution vector")
|
||||
assembly.u += u
|
||||
assembly.la += la
|
||||
assembly.la = la
|
||||
else
|
||||
info("$(problem.name): unknown formulation type, don't know what to do with results")
|
||||
error("serious failure with problem formulation: $(get_formulation_type(problem))")
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM.Core: Node, Seg2, update!, Problem, Mortar, assemble!
|
||||
|
||||
function get_testproblems(u, la)
|
||||
nodes = Dict{Int64, Node}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [2.0, 0.0],
|
||||
3 => [1.0, 2.0],
|
||||
4 => [0.0, 2.0],
|
||||
5 => [2.5, 0.0],
|
||||
6 => [4.5, 0.0],
|
||||
7 => [4.5, 1.0],
|
||||
8 => [2.5, 1.0])
|
||||
displacement = Dict{Int64, Vector{Float64}}()
|
||||
reaction_force = Dict{Int64, Vector{Float64}}()
|
||||
for i=1:8
|
||||
displacement[i] = u[:,i]
|
||||
reaction_force[i] = la[:,i]
|
||||
end
|
||||
bc5 = Seg2([3, 2])
|
||||
bc6 = Seg2([8, 5])
|
||||
update!([bc5, bc6], "geometry", nodes)
|
||||
update!([bc5, bc6], "displacement", displacement)
|
||||
update!([bc5, bc6], "reaction force", reaction_force)
|
||||
bc5["master elements"] = [bc6]
|
||||
contact1 = Problem(Mortar, "contact between bodies", 2, "displacement")
|
||||
contact2 = Problem(Mortar, "contact between bodies", 2, "displacement")
|
||||
contact2.properties.formulation = :forwarddiff
|
||||
contact2.assembly.u = vec(u)
|
||||
contact2.assembly.la = vec(la)
|
||||
push!(contact1, bc5, bc6)
|
||||
push!(contact2, bc5, bc6)
|
||||
return contact1, contact2
|
||||
end
|
||||
|
||||
@testset "test linearization of contact force in undeformed state" begin
|
||||
u = zeros(2, 8)
|
||||
la = zeros(2, 8)
|
||||
contact1, contact2 = get_testproblems(u, la)
|
||||
assemble!(contact1, 0.0)
|
||||
assemble!(contact2, 0.0)
|
||||
@test isapprox(full(contact1.assembly.C1), full(contact2.assembly.C1))
|
||||
@test isapprox(full(contact1.assembly.K), full(contact2.assembly.K))
|
||||
end
|
||||
Reference in New Issue
Block a user