mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-29 07:32:54 +00:00
3d mortar autodiff not working
This commit is contained in:
Jukka Aho
+1
-1
@@ -124,7 +124,7 @@ function get_basis{E}(::Type{Element{E}}, xi::Vector{Float64})
|
||||
return get_basis(E, xi)
|
||||
end
|
||||
|
||||
function get_basis{E}(element::Element{E}, xi::Vector{Float64})
|
||||
function get_basis{E}(element::Element{E}, xi::Vector)
|
||||
return get_basis(E, xi)
|
||||
end
|
||||
|
||||
|
||||
+6
-1
@@ -141,10 +141,15 @@ end
|
||||
xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
|
||||
|
||||
|
||||
# some helpers to make accessing 1d basis functions more easily
|
||||
# some helpers to make accessing 1d basis functions more easy
|
||||
function get_basis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T)
|
||||
get_basis(E, [xi])
|
||||
end
|
||||
|
||||
function get_dbasis{T<:Real, E<:Union{Seg2,Seg3}}(::Type{E}, xi::T)
|
||||
get_dbasis(E, [xi])
|
||||
end
|
||||
|
||||
function get_reference_element_midpoint{E}(element::Element{E})
|
||||
get_reference_element_midpoint(E)
|
||||
end
|
||||
|
||||
@@ -53,9 +53,21 @@ macro debug(msg)
|
||||
return msg
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
elements = get_elements(problem)
|
||||
if length(elements) == 0
|
||||
info("$(typeof(problem)) : forget to add elements?")
|
||||
return
|
||||
end
|
||||
# returns 3 if eldim 2 (tri3, quad4, ...) for 3d problems etc.
|
||||
eldim = size(elements[1], 1)+1
|
||||
assemble!(problem, time, Val{eldim})
|
||||
end
|
||||
|
||||
include("mortar_2d.jl")
|
||||
include("mortar_2d_autodiff.jl")
|
||||
include("mortar_3d.jl")
|
||||
include("mortar_3d_autodiff.jl")
|
||||
|
||||
""" Remove inactive inequality constraints by using primal-dual active set strategy. """
|
||||
function boundary_assembly_posthook!(solver::Solver, problem::Problem{Mortar}, C1, C2, D, g)
|
||||
|
||||
@@ -73,7 +73,8 @@ function project_from_slave_to_master{E<:MortarElements2D}(
|
||||
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
""" Assemble Mortar problem for two-dimensional problems, i.e. for Seg2 and Seg3 elements. """
|
||||
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
@@ -99,7 +100,7 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
push!(S, conn...)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
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)
|
||||
@@ -111,6 +112,7 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
for i in 1:size(normals,2)
|
||||
normals[:,i] /= norm(normals[:,i])
|
||||
end
|
||||
# swap element normals in 2d if they point to inside of body
|
||||
if props.rotate_normals
|
||||
for i=1:size(normals,2)
|
||||
normals[:,i] = -normals[:,i]
|
||||
@@ -127,8 +129,8 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
x1 = X1 + u1
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
n1 = Field(Vector[normals[:,i] for i in slave_element_nodes])
|
||||
|
||||
nnodes = size(slave_element, 2)
|
||||
update!(slave_element, "normals", time => ForwardDiff.get_value(n1))
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element["master elements"]
|
||||
@@ -202,7 +204,7 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
gap[1,slave_element_nodes] += w*gn*Phi'
|
||||
#gap[1,slave_element_nodes] += w*gn*N1'
|
||||
|
||||
end
|
||||
end # done integrating segment
|
||||
|
||||
end # master elements done
|
||||
|
||||
@@ -230,7 +232,6 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
C[1,j] += gap[1, j]
|
||||
C[2,j] += lat
|
||||
else
|
||||
#info("set node $j inactive")
|
||||
C[:,j] = la[:,j]
|
||||
end
|
||||
end
|
||||
|
||||
+59
-30
@@ -45,21 +45,25 @@ Notes
|
||||
[1](http://stackoverflow.com/questions/8942950/how-do-i-find-the-orthogonal-projection-of-a-point-onto-a-plane)
|
||||
|
||||
"""
|
||||
function project_point_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix)
|
||||
function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, Q::Matrix)
|
||||
n = Q[:,1]
|
||||
ph = p - dot(p-x0, n)*n
|
||||
qproj = Q'*(ph-x0)
|
||||
if !isapprox(qproj[1], 0.0; atol=1.0e-12)
|
||||
info("project_point_to_auxiliary_plane(): point not projected correctly.")
|
||||
info("p: $p")
|
||||
info("x0: $x0")
|
||||
info("Q: \n$Q")
|
||||
info("qproj: $qproj")
|
||||
error("Failed to project point to auxiliary plane.")
|
||||
if abs(qproj[1]) > 1.0e-2
|
||||
# we should have something very little for normal direction if projected
|
||||
# properly
|
||||
info("project_point_to_auxiliary_plane(): vertex not projected correctly.")
|
||||
info("p: $(ForwardDiff.get_value(p))")
|
||||
info("x0: $(ForwardDiff.get_value(x0))")
|
||||
info("Q: \n$(ForwardDiff.get_value(Q))")
|
||||
info("qproj: $(ForwardDiff.get_value(qproj))")
|
||||
error("Failed to project vertex to auxiliary plane.")
|
||||
end
|
||||
return qproj[2:3]
|
||||
end
|
||||
|
||||
project_point_to_auxiliary_plane = project_vertex_to_auxiliary_plane
|
||||
|
||||
"""
|
||||
Find edge intersections of two planar arbitrary shape polygons.
|
||||
|
||||
@@ -114,7 +118,7 @@ function get_edge_intersections(S::Matrix, M::Matrix)
|
||||
for j=1:nm
|
||||
b = M[:,j]-S[:,i]
|
||||
A = [S[:,mod(i,ns)+1]-S[:,i] -M[:,mod(j,nm)+1]+M[:,j]]
|
||||
if rank(A) == 2
|
||||
if rank(ForwardDiff.get_value(A)) == 2
|
||||
r = A\b
|
||||
if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found
|
||||
k += 1
|
||||
@@ -206,12 +210,12 @@ end
|
||||
"""
|
||||
function uniquetol(P, dim::Int; args...)
|
||||
@assert dim == 2
|
||||
items = Vector{Float64}[P[:,i] for i=1:size(P,dim)]
|
||||
new_items = Vector{Float64}[]
|
||||
items = Vector[P[:,i] for i=1:size(P,dim)]
|
||||
new_items = Vector[]
|
||||
for item in items
|
||||
has_found = false
|
||||
for new_item in new_items
|
||||
if isapprox(item, new_item; args...)
|
||||
if isapprox(ForwardDiff.get_value(item), ForwardDiff.get_value(new_item); args...)
|
||||
has_found = true
|
||||
break
|
||||
end
|
||||
@@ -317,7 +321,7 @@ function calculate_polygon_centerpoint(P::Matrix)
|
||||
Cx += 1/(6*A)*(P[1,i] + P[1,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i])
|
||||
Cy += 1/(6*A)*(P[2,i] + P[2,inext])*(P[1,i]*P[2,inext] - P[1,inext]*P[2,i])
|
||||
end
|
||||
return Float64[Cx, Cy]
|
||||
return [Cx, Cy]
|
||||
end
|
||||
|
||||
"""
|
||||
@@ -401,33 +405,58 @@ julia> xquad*basis(theta[2:3])
|
||||
1.0
|
||||
|
||||
"""
|
||||
function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix, element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
|
||||
function project_point_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix,
|
||||
element::Element{E}, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
|
||||
x = element("geometry", time)
|
||||
return project_point_from_plane_to_surface(p, x0, Q, element, x, time;
|
||||
max_iterations=max_iterations, iter_tol=iter_tol)
|
||||
end
|
||||
function project_vertex_from_plane_to_surface{E}(p::Vector, x0::Vector, Q::Matrix,
|
||||
element::Element{E}, x, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
|
||||
basis(xi) = get_basis(E, xi)
|
||||
dbasis(xi) = get_dbasis(E, xi)
|
||||
x = element("geometry", time)
|
||||
ph = Q*[0; p] + x0
|
||||
theta = Float64[0.0, 0.0, 0.0]
|
||||
n = Q[:,1]
|
||||
b(theta) = ph + theta[1]*n - basis(theta[2:3])*x
|
||||
J(theta) = [n -dbasis(theta[2:3])*x]
|
||||
theta = zeros(3)
|
||||
dtheta = zeros(3)
|
||||
for i=1:max_iterations
|
||||
b = ph + theta[1]*n - basis(theta[2:3])*x
|
||||
J = [n -dbasis(theta[2:3])*x]
|
||||
dtheta = J \ -b
|
||||
# FIXME: gives NaN if partials in J
|
||||
dtheta = ForwardDiff.get_value(J(theta)) \ -b(theta)
|
||||
theta += dtheta
|
||||
if norm(ForwardDiff.get_value(dtheta)) < iter_tol
|
||||
return theta
|
||||
end
|
||||
end
|
||||
info("failed to project vertex from auxiliary plane back to surface")
|
||||
info("element type: $E")
|
||||
info("element connectivity: $(get_connectivity(element))")
|
||||
info("auxiliary plane: x0 = $(ForwardDiff.get_value(x0)), Q = $(ForwardDiff.get_value(Q))")
|
||||
info("point coordinates on plane: $(ForwardDiff.get_value(p))")
|
||||
info("element geometry: $(ForwardDiff.get_value(x.data))")
|
||||
info("ph: $(ForwardDiff.get_value(ph))")
|
||||
info("normal direction: $(ForwardDiff.get_value(n))")
|
||||
info("parameter vector before giving up: $(ForwardDiff.get_value(theta))")
|
||||
info("increment in parameter vector before giving up: $(ForwardDiff.get_value(dtheta))")
|
||||
info("b([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(b([0.0, 0.0, 0.0])))")
|
||||
info("J([0.0, 0.0, 0.0]) = $(ForwardDiff.get_value(J([0.0, 0.0, 0.0])))")
|
||||
|
||||
info("iterations were")
|
||||
theta = zeros(3)
|
||||
dtheta = zeros(3)
|
||||
for i=1:max_iterations
|
||||
info("iter $i, theta = $(ForwardDiff.get_value(theta))")
|
||||
info("b = $(ForwardDiff.get_value(b(theta)))")
|
||||
info("J = $(ForwardDiff.get_value(J(theta)))")
|
||||
dtheta = ForwardDiff.get_value(J(theta)) \ -b(theta)
|
||||
info("dtheta = $(ForwardDiff.get_value(dtheta))")
|
||||
theta += dtheta
|
||||
if norm(dtheta) < iter_tol
|
||||
return theta
|
||||
end
|
||||
end
|
||||
begin
|
||||
info("projecting point from auxiliary plane back to surface didn't go very well.")
|
||||
info("element type: $E")
|
||||
info("element connectivity: $(get_connectivity(element))")
|
||||
info("auxiliary plane: x0 = $x0, Q = $Q")
|
||||
info("point coordinates on plane: $p")
|
||||
info("element geometry: $x")
|
||||
info("ph: $ph")
|
||||
info("normal direction: $n")
|
||||
info("parameter vector before giving up: $theta")
|
||||
end
|
||||
|
||||
error("project_point_to_surface: did not converge in $max_iterations iterations!")
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,312 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
""" Assemble Mortar problem for three-dimensional problems, i.e. for Tri3, Tri6, Quad4, Quad8, Quad9 elements. """
|
||||
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{3}})
|
||||
|
||||
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) # x = [u; la]
|
||||
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)
|
||||
gap = zeros(u)
|
||||
gap_added = zeros(size(u)...)
|
||||
C = zeros(la)
|
||||
all_slave_nodes = Set{Int64}()
|
||||
|
||||
# 1. calculate and average node normals for slave element nodes
|
||||
normal = zeros(u)
|
||||
tangent1 = zeros(u)
|
||||
tangent2 = zeros(u)
|
||||
for element in get_elements(problem)
|
||||
haskey(element, "master elements") || continue
|
||||
conn = get_connectivity(element)
|
||||
push!(all_slave_nodes, conn...)
|
||||
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)
|
||||
j = transpose(sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)]))
|
||||
#info("size of j = $(size(j))")
|
||||
n = reshape(cross(j[:,1], j[:,2]), 3, 1)
|
||||
normal[:, conn] += ip.weight*n*N
|
||||
end
|
||||
end
|
||||
# calculate tangents
|
||||
for i in 1:size(normal, 2)
|
||||
i in all_slave_nodes || continue
|
||||
normal[:,i] /= norm(normal[:,i])
|
||||
u1 = normal[:,i]
|
||||
j = indmax(abs(u1))
|
||||
v2 = zeros(3)
|
||||
v2[mod(j,3)+1] = 1.0
|
||||
u2 = v2 - dot(u1,v2)/dot(v2,v2)*v2
|
||||
u3 = cross(u1,u2)
|
||||
tangent1[:,i] = u2/norm(u2)
|
||||
tangent2[:,i] = u3/norm(u3)
|
||||
end
|
||||
if props.rotate_normals
|
||||
for i=1:size(normal, 2)
|
||||
normal[:,i] = -normal[:,i]
|
||||
end
|
||||
end
|
||||
|
||||
# 2. loop slave elements and find contact segments
|
||||
for slave_element in get_elements(problem)
|
||||
haskey(slave_element, "master elements") || continue
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
n1 = Field(Vector[normal[:,i] for i in slave_element_nodes])
|
||||
t1 = Field(Vector[tangent1[:,i] for i in slave_element_nodes])
|
||||
t2 = Field(Vector[tangent2[:,i] for i in slave_element_nodes])
|
||||
nnodes = size(slave_element, 2)
|
||||
update!(slave_element, "normals", time => ForwardDiff.get_value(n1.data))
|
||||
|
||||
# create auxiliary plane (x0, Q)
|
||||
xi = get_reference_element_midpoint(slave_element)
|
||||
N = vec(get_basis(slave_element, xi))
|
||||
x0 = N*x1
|
||||
Q = [N*n1 N*t1 N*t2]
|
||||
|
||||
# project slave nodes to auxiliary plane
|
||||
S = hcat([project_vertex_to_auxiliary_plane(p, x0, Q) for p in x1]...)
|
||||
|
||||
# 3. loop all master elements
|
||||
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
|
||||
|
||||
x1_midpoint = mean(x1)
|
||||
x2_midpoint = mean(x2)
|
||||
distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
|
||||
distance > props.maximum_distance && continue
|
||||
|
||||
# project master nodes to auxiliary plane
|
||||
M = hcat([project_vertex_to_auxiliary_plane(p, x0, Q) for p in x2]...)
|
||||
|
||||
# create polygon clipping on auxiliary plane
|
||||
#=
|
||||
P = nothing
|
||||
neighbours = nothing
|
||||
try
|
||||
catch
|
||||
info("polygon clipping failed")
|
||||
info("S = ")
|
||||
dump(ForwardDiff.get_value(S))
|
||||
info("M = ")
|
||||
dump(ForwardDiff.get_value(M))
|
||||
error("cannot continue")
|
||||
end
|
||||
=#
|
||||
P, neighbours = clip_polygon(S, M)
|
||||
isa(P, Void) && continue # no clipping
|
||||
info("polygon clip found: S = $(ForwardDiff.get_value(S)), M = $(ForwardDiff.get_value(M)), P = $(ForwardDiff.get_value(P))")
|
||||
# special case, shared edge but no shared volume
|
||||
size(P, 2) < 3 && continue
|
||||
gap_added += 1
|
||||
|
||||
# clip polygon centerpoint
|
||||
#C0 = calculate_polygon_centerpoint(P)
|
||||
C0 = vec(mean(P, 2))
|
||||
npts = size(P, 2) # number of vertices in polygon
|
||||
for pnt=1:npts # loop integration cells
|
||||
x_cell = Field(Vector[C0, P[:,pnt], P[:,mod(pnt,npts)+1]])
|
||||
#=
|
||||
try
|
||||
catch
|
||||
info("centerpoint: $(ForwardDiff.get_value(C0))")
|
||||
info("P1 = $(ForwardDiff.get_value(P[:,pnt]))")
|
||||
info("P2 = $(ForwardDiff.get_value(P[:,mod(pnt,npts)+1]))")
|
||||
data = Vector[C0, P[:,pnt], P[:,mod(pnt,npts)+1]]
|
||||
info("data = $(ForwardDiff.get_value(data))")
|
||||
rethrow()
|
||||
end
|
||||
=#
|
||||
# create dual basis
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = vec(get_basis(Tri3, ip.xi))
|
||||
x_g = N*x_cell
|
||||
theta = zeros(3)
|
||||
try
|
||||
theta = project_vertex_from_plane_to_surface(x_g, x0, Q, slave_element, x1, time)
|
||||
catch
|
||||
info("creating dual basis did fail for finding projection back to slave surface.")
|
||||
info("cell coords in auxiliary plane are:")
|
||||
info(ForwardDiff.get_value(x_cell.data))
|
||||
info("interpolated value x_g is $(ForwardDiff.get_value(x_g))")
|
||||
info("clip polygon centerpoint is $(ForwardDiff.get_value(C0))")
|
||||
info("clip polygon is $(ForwardDiff.get_value(P))")
|
||||
info("slave side nodes projected onto a plane are $(ForwardDiff.get_value(S))")
|
||||
info("master side nodes projected onto a plane are $(ForwardDiff.get_value(M))")
|
||||
rethrow()
|
||||
end
|
||||
xi_slave = theta[2:3]
|
||||
N1 = slave_element(xi_slave, time)
|
||||
|
||||
# jacobian determinant on integration cell
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
|
||||
De += wC*diagm(vec(N1))
|
||||
Me += wC*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
# loop integration points of cell
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = vec(get_basis(Tri3, ip.xi))
|
||||
x_g = N*x_cell
|
||||
# project gauss point back to element surfaces
|
||||
theta1 = project_vertex_from_plane_to_surface(x_g, x0, Q, slave_element, x1, time)
|
||||
theta2 = project_vertex_from_plane_to_surface(x_g, x0, Q, master_element, x2, time)
|
||||
xi_slave = theta1[2:3]
|
||||
xi_master = theta2[2:3]
|
||||
|
||||
# evaluate shape functions, calculate contact force and gap
|
||||
N1 = vec(get_basis(slave_element, xi_slave))
|
||||
N2 = vec(get_basis(master_element, xi_master))
|
||||
Phi = Ae*N1
|
||||
|
||||
# jacobian determinant of integration cell
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
|
||||
x_s = N1*x1
|
||||
n_s = N1*n1
|
||||
x_m = N2*x2
|
||||
la_s = Phi*la1
|
||||
gn = props.gap_sign*dot(n_s, x_s - x_m)
|
||||
fc[:,slave_element_nodes] += wC*la_s*N1'
|
||||
fc[:,master_element_nodes] -= wC*la_s*N2'
|
||||
wg = wC*gn*Phi'
|
||||
if any(isnan(wg))
|
||||
info("gap has NaNs!")
|
||||
info("wC = $(ForwardDiff.get_value(wC))")
|
||||
info("gn = $(ForwardDiff.get_value(gn))")
|
||||
info("Phi = $(ForwardDiff.get_value(Phi))")
|
||||
info("xi_slave = $(ForwardDiff.get_value(xi_slave))")
|
||||
info("N1 = $(ForwardDiff.get_value(N1))")
|
||||
info("Ae = $(ForwardDiff.get_value(Ae))")
|
||||
info("De = $(ForwardDiff.get_value(De))")
|
||||
info("Me = $(ForwardDiff.get_value(Me))")
|
||||
info("x_cell(data) = $(ForwardDiff.get_value(x_cell.data))")
|
||||
info("C0 = $(ForwardDiff.get_value(C0))")
|
||||
info("P = $(ForwardDiff.get_value(P))")
|
||||
error("fix this")
|
||||
end
|
||||
gap[1,slave_element_nodes] += wg
|
||||
gap_added[1,slave_element_nodes] += 1
|
||||
end # done integrating cell
|
||||
|
||||
end # done for all cells in this segment
|
||||
|
||||
end # done all master elements for this slave element
|
||||
|
||||
end # done all slave elements
|
||||
|
||||
# like in 2d, check contact in nodes based on a complementarity condition
|
||||
|
||||
nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
|
||||
all_slave_nodes = sort(collect(all_slave_nodes))
|
||||
info("gap: $nzgap")
|
||||
info("size of normal = $(size(normal))")
|
||||
info("size of la = $(size(la))")
|
||||
info("size of C = $(size(C))")
|
||||
info("S = $all_slave_nodes")
|
||||
|
||||
for (i, j) in enumerate(all_slave_nodes)
|
||||
if j in props.always_inactive
|
||||
info("special node $j always inactive")
|
||||
C[:,j] = la[:,j]
|
||||
continue
|
||||
end
|
||||
n = normal[:,j]
|
||||
t1 = tangent1[:,j]
|
||||
t2 = tangent2[:,j]
|
||||
lan = dot(n, la[:,j])
|
||||
|
||||
# if lan - gap[1, j] > 0
|
||||
info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t1)) x $(ForwardDiff.get_value(t2))")
|
||||
C[1,j] += gap[1, j]
|
||||
C[2,j] += dot(t1, la[:,j])
|
||||
C[3,j] += dot(t2, la[:,j])
|
||||
# else
|
||||
# C[:,j] = la[:,j]
|
||||
# end
|
||||
end
|
||||
|
||||
for (i, j) in enumerate(all_slave_nodes)
|
||||
Ci = ForwardDiff.get_value(C[:,j])
|
||||
gapi = ForwardDiff.get_value(gap[:,j])
|
||||
fci = ForwardDiff.get_value(fc[:,j])
|
||||
lai = ForwardDiff.get_value(la[:,j])
|
||||
ui = ForwardDiff.get_value(u[:,j])
|
||||
ni = ForwardDiff.get_value(normal[:,j])
|
||||
gai = gap_added[:,j]
|
||||
info("$i/$j: C = $Ci, f = $fci, gap = $gapi, la = $lai, u = $ui, n = $ni, gai = $gai")
|
||||
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)
|
||||
A, allresults = ForwardDiff.jacobian(calculate_interface, x,
|
||||
ForwardDiff.AllResults, cache=autodiffcache)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
|
||||
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]
|
||||
|
||||
slaves = [101,108,111,112,113,120,123,124,125,126,129,130,149,150,151,152]
|
||||
for j in slaves
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
info("slave node $j, dofs $dofs")
|
||||
info("Stiffness: $(K[dofs,:])")
|
||||
info("force fc: $(C1[dofs,:])")
|
||||
info("constraint: $(C2[dofs,:])")
|
||||
info("lambdas: $(D[dofs,:])")
|
||||
info("f = $(f[dofs]), g = $(g[dofs])")
|
||||
end
|
||||
|
||||
empty!(problem.assembly)
|
||||
add!(problem.assembly.K, K)
|
||||
add!(problem.assembly.C1, C1)
|
||||
add!(problem.assembly.C2, C2)
|
||||
add!(problem.assembly.D, D)
|
||||
add!(problem.assembly.f, f)
|
||||
add!(problem.assembly.g, g)
|
||||
|
||||
return problem.assembly
|
||||
end
|
||||
|
||||
@@ -3,13 +3,13 @@
|
||||
|
||||
using HDF5
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Core: Element, Quad4, Tri3, Seg2, Hex8, update!
|
||||
using JuliaFEM.Core: Element, Quad4, Tri3, Tet4, Seg2, Hex8, update!
|
||||
|
||||
|
||||
# TODO: this should be elsewhere
|
||||
function aster_create_elements(mesh, element_set, element_type=nothing)
|
||||
elements = Element[]
|
||||
mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8)
|
||||
mapping = Dict(:QU4 => Quad4, :TR3 => Tri3, :SE2 => Seg2, :HE8 => Hex8, :TE4 => Tet4)
|
||||
for (elid, (eltype, elset, elcon)) in mesh["connectivity"]
|
||||
if !haskey(mapping, eltype)
|
||||
error("aster_create_elements: unknown element mapping $eltype")
|
||||
|
||||
Reference in New Issue
Block a user