Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-09-19 09:54:55 +00:00
Code Issues Packages Projects Releases Wiki Activity

some refactoring

Browse Source
This commit is contained in:
Jukka Aho
2016-02-24 10:35:18 +02:00
parent dc7497c5c9
commit e0558beffa
3 changed files with 138 additions and 103 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/mortar.jl
+4 -2
View File
@@ -24,18 +24,20 @@ type Mortar <: BoundaryProblem
inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
normal_condition :: Symbol # Tie or Contact
tangential_condition :: Symbol # Stick or Slip
minimum_distance :: Float64 # don't check for a contact if elements are far enough
maximum_distance :: Float64 # don't check for a contact if elements are far enough
store_debug_info :: Bool # for making debugging easier
always_inactive :: Vector{Int64}
always_in_contact :: Vector{Int64} # nodes in this list always in contact
always_in_stick :: Vector{Int64} # nodes in this list always in stick
always_in_slip :: Vector{Int64} # nodes in this list always in slip
contact :: Bool
friction :: Bool
gap_sign :: Int # gap sign convention
rotate_normals :: Bool
end
function Mortar()
Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false, -1)
Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], [], false, false, -1, false)
end
function get_unknown_field_name(::Type{Mortar})
src/mortar_2d_autodiff.jl
+131 -100
View File
@@ -73,10 +73,8 @@ function project_from_slave_to_master{E<:MortarElements2D}(
end
function assemble!{E<:MortarElements2D}(assembly::Assembly,
problem::Problem{Mortar}, slave_element::Element{E},
time::Real, ::Type{Val{:forwarddiff}})
haskey(slave_element, "master elements") || return
function assemble!(problem::Problem{Mortar}, time::Real)
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
@@ -88,21 +86,17 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
u = reshape(x[1:ndofs], field_dim, nnodes)
la = reshape(x[ndofs+1:end], field_dim, nnodes)
fc = zeros(u)
gap = zeros(u)
C = zeros(la)
S = Set{Int64}()
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
# 1. update nodal normals for slave elements
Q = [0.0 -1.0; 1.0 0.0]
normals = zeros(u)
for element in adjacent_elements
for element in get_elements(problem)
haskey(element, "master elements") || continue
conn = get_connectivity(element)
push!(S, conn...)
gdofs = get_gdofs(element, field_dim)
X_el = element("geometry", time)
u_el = Field(Vector[u[:, i] for i in conn])
@@ -114,94 +108,129 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
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])
nnodes = size(slave_element, 2)
lan_tot = zeros(nnodes) # normal pressure
gap_tot = zeros(nnodes) # weighted normal gap
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
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
for i in 1:size(normals,2)
normals[:,i] /= norm(normals[:,i])
end
if props.rotate_normals
for i=1:size(normals,2)
normals[:,i] = -normals[:,i]
end
Ae = De*inv(Me)
end
slave_dofs = get_gdofs(slave_element, field_dim)
master_dofs = get_gdofs(master_element, field_dim)
# 2. loop all slave elements
for slave_element in get_elements(problem)
haskey(slave_element, "master elements") || continue
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
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[normals[:,i] for i in slave_element_nodes])
# 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
nnodes = size(slave_element, 2)
u_s = N1*u1
u_m = N2*u2
# 3. loop all master elements
for master_element in slave_element["master elements"]
la_s = Phi*la1 # traction force in gauss point
la_nt = R_s*la_s
gn = -dot(n_s, x_s - x_m) # normal gap
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
fc[:,slave_element_nodes] += w*la_s*N1'
fc[:,master_element_nodes] -= w*la_s*N2'
#C[1,slave_element_nodes] += w*gn*Phi'
x1_midpoint = 1/2*(x1[1]+x1[2])
x2_midpoint = 1/2*(x2[1]+x2[2])
distance = ForwardDiff.get_value(norm(x2_midpoint - x1_midpoint))
distance > props.maximum_distance && continue
lan_tot += w*la_nt[1]*Phi
gap_tot += w*gn*Phi
# calculate segmentation: we care only about endpoints
# note: these are quadratic/cubic functions, analytical solution possible
xi1a = -Inf
xi1b = -Inf
try
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
catch
info("failed to create projection!!!!")
# TODO
continue
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
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)
# 4. loop integration points of segment
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
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
la_s = Phi*la1 # traction force in gauss point
gn = props.gap_sign*dot(n_s, x_s - x_m) # normal gap
fc[:,slave_element_nodes] += w*la_s*N1'
fc[:,master_element_nodes] -= w*la_s*N2'
gap[1,slave_element_nodes] += w*gn*Phi'
#gap[1,slave_element_nodes] += w*gn*N1'
end
end # master elements done
end # slave elements done
# at this point we have calculated contact force fc and gap for all slave elements.
# next task is to find out are they in contact or not and remove inactive nodes
nzgap = sort(nonzeros(sparse(ForwardDiff.get_value(gap))))
info("gap: $nzgap")
for (i, j) in enumerate(sort(collect(S)))
if j in props.always_inactive
info("special node $j always inactive")
C[:,j] = la[:,j]
continue
end
end # master elements done
for (i, j) in enumerate(slave_element_nodes)
n = n1[i]
n = normals[:,j]
t = Q'*n
R = [n t]
la_nt = R'*la[:,j]
# info("node $j, n=$(ForwardDiff.get_value(n)) lan = $(ForwardDiff.get_value(la_nt[1])) gap = $(ForwardDiff.get_value(gap_tot[i]))")
lan = dot(n, la[:,j])
lat = dot(t, la[:,j])
if la_nt[1] - gap_tot[i] > 0
info("set node $j active")
C[1,j] += gap_tot[i]
C[2,j] += la_nt[2]
if lan - gap[1, j] > 0
info("set node $j active, normal direction = $(ForwardDiff.get_value(n)), tangent plane = $(ForwardDiff.get_value(t))")
C[1,j] += gap[1, j]
C[2,j] += lat
else
info("set node $j inactive")
#info("set node $j inactive")
C[:,j] = la[:,j]
end
end
@@ -213,28 +242,30 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly,
# 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)
A, allresults = ForwardDiff.jacobian(calculate_interface, x,
ForwardDiff.AllResults, cache=autodiffcache)
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
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
src/solver_utils.jl
+3 -1
View File
@@ -70,6 +70,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
# INFO: fixed: new setting is
# INFO: dof 1109: 0.0*u₁₅₃ - 0.0*u₁₅₄ - 0.0*u₁₅₅ + 0.15*u₁₅₆ + 0.0*u₁₁₀₉ - 0.15*u₁₁₁₀ = -0.0
=#
#=
if 555 in nodes
info("overconstraint DIRTY HACK")
# It is possible to selectively remove mortar constraints and the associated
@@ -103,6 +104,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
#C1_[1110,:] = 0
return
end
=#
""" Return all other dofs which connects to overconstrained dofs. """
function get_related_dofs(dofs)
@@ -243,7 +245,7 @@ function handle_overconstraint_error!(problem, nodes, all_dofs, C1_, C1, C2_, C2
return dofs_, true
end
actions = [action1, action2]
actions = [action1, action2, action3, action4]
function show_lambda_coefficients(dofs, C1)
for dof in dofs