mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-19 01:48:47 +00:00
modal solver + tie contact works now. dirichlet boundary and mpcs are eliminated properly before solution to get reduced system
This commit is contained in:
Jukka Aho
+3
-2
@@ -66,7 +66,8 @@ include("solvers.jl")
|
||||
export AbstractSolver, Solver, Nonlinear,
|
||||
get_unknown_field_name, get_formulation_type,
|
||||
get_field_problems, get_boundary_problems,
|
||||
get_field_assembly, get_boundary_assembly
|
||||
get_field_assembly, get_boundary_assembly,
|
||||
initialize!, create_projection
|
||||
include("modal.jl")
|
||||
export Modal
|
||||
|
||||
@@ -107,7 +108,7 @@ end
|
||||
|
||||
module Postprocess
|
||||
include("postprocess_utils.jl")
|
||||
export calc_nodal_values!
|
||||
export calc_nodal_values!, get_nodal_vector
|
||||
include("postprocess_xdmf.jl")
|
||||
export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!
|
||||
end
|
||||
|
||||
+37
-18
@@ -102,7 +102,11 @@ function get_integration_points(element::TriangularElement, ::Type{Val{2}})
|
||||
end
|
||||
|
||||
function get_integration_points(element::TriangularElement, ::Type{Val{3}})
|
||||
weights = 0.5*[-0.5625, 0.5208333333333333, 0.5208333333333333, 0.5208333333333333]
|
||||
weights = 0.5*[
|
||||
-0.5625,
|
||||
0.5208333333333333,
|
||||
0.5208333333333333,
|
||||
0.5208333333333333]
|
||||
points = Vector{Float64}[
|
||||
[1.0/3.0, 1.0/3.0],
|
||||
[0.2, 0.2],
|
||||
@@ -112,26 +116,41 @@ function get_integration_points(element::TriangularElement, ::Type{Val{3}})
|
||||
end
|
||||
|
||||
function get_integration_points(element::TriangularElement, ::Type{Val{4}})
|
||||
[
|
||||
IntegrationPoint([0.44594849091597, 0.44594849091597], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.44594849091597, 0.10810301816807], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.10810301816807, 0.44594849091597], 0.5*0.22338158967801),
|
||||
IntegrationPoint([0.09157621350977, 0.09157621350977], 0.5*0.10995174365532),
|
||||
IntegrationPoint([0.09157621350977, 0.81684757298046], 0.5*0.10995174365532),
|
||||
IntegrationPoint([0.81684757298046, 0.09157621350977], 0.5*0.10995174365532)
|
||||
]
|
||||
weights = 0.5*[
|
||||
0.22338158967801,
|
||||
0.22338158967801,
|
||||
0.22338158967801,
|
||||
0.10995174365532,
|
||||
0.10995174365532,
|
||||
0.10995174365532]
|
||||
points = Vector{Float64}[
|
||||
[0.44594849091597, 0.44594849091597],
|
||||
[0.44594849091597, 0.10810301816807],
|
||||
[0.10810301816807, 0.44594849091597],
|
||||
[0.09157621350977, 0.09157621350977],
|
||||
[0.09157621350977, 0.81684757298046],
|
||||
[0.81684757298046, 0.09157621350977]]
|
||||
return zip(weights, points)
|
||||
end
|
||||
|
||||
function get_integration_points(element::TriangularElement, ::Type{Val{5}})
|
||||
[
|
||||
IntegrationPoint([0.33333333333333, 0.33333333333333], 0.5*0.22500000000000),
|
||||
IntegrationPoint([0.47014206410511, 0.47014206410511], 0.5*0.13239415278851),
|
||||
IntegrationPoint([0.47014206410511, 0.05971587178977], 0.5*0.13239415278851),
|
||||
IntegrationPoint([0.05971587178977, 0.47014206410511], 0.5*0.13239415278851),
|
||||
IntegrationPoint([0.10128650732346, 0.10128650732346], 0.5*0.12593918054483),
|
||||
IntegrationPoint([0.10128650732346, 0.79742698535309], 0.5*0.12593918054483),
|
||||
IntegrationPoint([0.79742698535309, 0.10128650732346], 0.5*0.12593918054483)
|
||||
]
|
||||
weights = 0.5*[
|
||||
0.22500000000000,
|
||||
0.13239415278851,
|
||||
0.13239415278851,
|
||||
0.13239415278851,
|
||||
0.12593918054483,
|
||||
0.12593918054483,
|
||||
0.12593918054483]
|
||||
points = Vector{Float64}[
|
||||
[0.33333333333333, 0.33333333333333],
|
||||
[0.47014206410511, 0.47014206410511],
|
||||
[0.47014206410511, 0.05971587178977],
|
||||
[0.05971587178977, 0.47014206410511],
|
||||
[0.10128650732346, 0.10128650732346],
|
||||
[0.10128650732346, 0.79742698535309],
|
||||
[0.79742698535309, 0.10128650732346]]
|
||||
return zip(weights, points)
|
||||
end
|
||||
|
||||
### 3d elements
|
||||
|
||||
+39
-19
@@ -1,6 +1,16 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
"""
|
||||
Examples
|
||||
--------
|
||||
|
||||
julia> problems = get_problems()
|
||||
julia> solver = Solver(Modal)
|
||||
julia> push!(solver, problems...)
|
||||
julia> call(solver)
|
||||
|
||||
"""
|
||||
type Modal <: AbstractSolver
|
||||
geometric_stiffness :: Bool
|
||||
eigvals :: Vector
|
||||
@@ -21,40 +31,50 @@ function call(solver::Solver{Modal}; debug=false)
|
||||
assemble!(problem, solver.time)
|
||||
assemble!(problem, solver.time, Val{:mass_matrix})
|
||||
end
|
||||
for problem in get_boundary_problems(solver)
|
||||
assemble!(problem, solver.time)
|
||||
end
|
||||
t1 = round(toq(), 2)
|
||||
info("Assembled in $t1 seconds.")
|
||||
M, K, Kg, f = get_field_assembly(solver; with_mass_matrix=true)
|
||||
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
|
||||
K = K + Kb
|
||||
f = f + fb
|
||||
if solver.properties.geometric_stiffness
|
||||
K += Kg
|
||||
end
|
||||
for problem in get_boundary_problems(solver)
|
||||
assemble!(problem, solver.time)
|
||||
# FIXME: Check for tie contacts and rhs. Here we just
|
||||
# remove all fixed dofs giving funny results if problem
|
||||
# is having MPCs or non-homogeneous Dirichlet conditions
|
||||
# eliminate!(M, K, Kg, f, problem)
|
||||
fixed_dofs = get_nonzero_rows(problem.assembly.C2)
|
||||
K[fixed_dofs, :] = 0
|
||||
K[:, fixed_dofs] = 0
|
||||
M[fixed_dofs, :] = 0
|
||||
M[:, fixed_dofs] = 0
|
||||
f[fixed_dofs, :] = 0
|
||||
end
|
||||
fd = get_nonzero_rows(K)
|
||||
|
||||
@assert nnz(D) == 0
|
||||
@assert C1 == C2
|
||||
|
||||
tic()
|
||||
P, h = create_projection(C1, g)
|
||||
K_red = P'*K*P
|
||||
M_red = P'*M*P
|
||||
t1 = round(toq(), 2)
|
||||
info("Eliminated dirichlet boundaries in $t1 seconds.")
|
||||
|
||||
nz = get_nonzero_rows(K_red)
|
||||
ndofs = solver.ndofs
|
||||
props = solver.properties
|
||||
info("Calculate $(props.nev) eigenvalues...")
|
||||
if debug
|
||||
if debug && length(nz) < 100
|
||||
info("Stiffness matrix:")
|
||||
dump(round(full(K[fd, fd])))
|
||||
dump(round(full(K[nz, nz])))
|
||||
info("Mass matrix:")
|
||||
dump(round(full(M[fd, fd])))
|
||||
dump(round(full(M[nz, nz])))
|
||||
end
|
||||
|
||||
tic()
|
||||
om2, X = eigs(K[fd, fd], M[fd, fd]; nev=props.nev, which=props.which)
|
||||
om2, X = eigs(K_red[nz,nz], M_red[nz,nz]; nev=props.nev, which=props.which)
|
||||
props.eigvals = om2
|
||||
props.eigvecs = zeros(ndofs, length(om2))
|
||||
props.eigvecs[fd, :] = X
|
||||
v = zeros(ndofs)
|
||||
for i=1:length(om2)
|
||||
fill!(v, 0.0)
|
||||
v[nz] = X[:,i]
|
||||
props.eigvecs[:,i] = P*v + g
|
||||
end
|
||||
t1 = round(toq(), 2)
|
||||
info("Eigenvalues computed in $t1 seconds. Eigenvalues: $om2")
|
||||
return true
|
||||
|
||||
+373
-21
@@ -2,13 +2,14 @@
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
type Mortar <: BoundaryProblem
|
||||
dimension :: Int
|
||||
rotate_normals :: Bool
|
||||
adjust :: Bool
|
||||
tolerance :: Float64
|
||||
end
|
||||
|
||||
function Mortar()
|
||||
return Mortar(false, false, 0.0)
|
||||
return Mortar(-1, false, false, 0.0)
|
||||
end
|
||||
|
||||
function get_unknown_field_name(::Type{Mortar})
|
||||
@@ -37,6 +38,10 @@ function cross2(a, b)
|
||||
cross([a; 0], [b; 0])[3]
|
||||
end
|
||||
|
||||
function get_slave_elements(problem::Problem{Mortar})
|
||||
filter(el -> haskey(el, "master elements"), get_elements(problem))
|
||||
end
|
||||
|
||||
function project_from_master_to_slave{E<:MortarElements2D}(slave_element::Element{E}, x2, time)
|
||||
x1_ = slave_element["geometry"](time)
|
||||
n1_ = slave_element["normal"](time)
|
||||
@@ -61,7 +66,7 @@ function project_from_slave_to_master{E<:MortarElements2D}(master_element::Eleme
|
||||
return xi2
|
||||
end
|
||||
|
||||
function calculate_normals(elements, time, rotate_normals=false)
|
||||
function calculate_normals(elements, time, ::Type{Val{1}}; rotate_normals=false)
|
||||
tangents = Dict{Int64, Vector{Float64}}()
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
@@ -79,7 +84,7 @@ function calculate_normals(elements, time, rotate_normals=false)
|
||||
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = Dict{Int64, Vector{Float64}}()
|
||||
S = sort(collect(keys(tangents)))
|
||||
S = collect(keys(tangents))
|
||||
for j in S
|
||||
tangents[j] /= norm(tangents[j])
|
||||
normals[j] = Q*tangents[j]
|
||||
@@ -94,8 +99,8 @@ function calculate_normals(elements, time, rotate_normals=false)
|
||||
return normals, tangents
|
||||
end
|
||||
|
||||
function calculate_normals!(elements, time, rotate_normals=false)
|
||||
normals, tangents = calculate_normals(elements, time, rotate_normals)
|
||||
function calculate_normals!(elements, time, ::Type{Val{1}}; rotate_normals=false)
|
||||
normals, tangents = calculate_normals(elements, time, Val{1}; rotate_normals=rotate_normals)
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
update!(element, "normal", time => [normals[j] for j in conn])
|
||||
@@ -104,32 +109,34 @@ function calculate_normals!(elements, time, rotate_normals=false)
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
if problem.dimension == -1
|
||||
error("set interface dimension: problem.properties.dimension = 1 or 2")
|
||||
end
|
||||
assemble!(problem, time, Val{problem.properties.dimension})
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{1}})
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_elements = filter(el -> haskey(el, "master elements"), get_elements(problem))
|
||||
slave_elements = get_slave_elements(problem)
|
||||
|
||||
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
|
||||
normals, tangents = calculate_normals(slave_elements, time, props.rotate_normals)
|
||||
normals, tangents = calculate_normals(slave_elements, time, Val{1};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", normals)
|
||||
update!(slave_elements, "tangent", tangents)
|
||||
S = sort(collect(keys(normals)))
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
haskey(slave_element, "master elements") || continue
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element["geometry"](time)
|
||||
n1 = Field([normals[j] for j in slave_element_nodes])
|
||||
n1 = slave_element["normal"](time)
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element["master elements"](time)
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
nm = length(master_element)
|
||||
X2 = master_element["geometry"](time)
|
||||
|
||||
# 3.1 calculate segmentation
|
||||
@@ -141,9 +148,11 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
|
||||
# 3.3. loop integration points of one integration segment and calculate
|
||||
# local mortar matrices
|
||||
nsl = length(slave_element)
|
||||
nm = length(master_element)
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
ge = zeros(nsl)
|
||||
ge = zeros(field_dim*nsl)
|
||||
for ip in get_integration_points(slave_element, 2)
|
||||
detJ = slave_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ*l
|
||||
@@ -159,10 +168,11 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
De += w*N1*N1'
|
||||
Me += w*N1*N2'
|
||||
if props.adjust
|
||||
g = X_s-X_m
|
||||
if g < props.tol
|
||||
ge += w*g
|
||||
end
|
||||
u1 = slave_element["displacement"](time)
|
||||
u2 = master_element["displacement"](time)
|
||||
x_s = X_s + N1*u1
|
||||
x_m = X_m + N2*u2
|
||||
ge += w*vec((x_m-x_s)*N1')
|
||||
end
|
||||
end
|
||||
|
||||
@@ -177,12 +187,354 @@ function assemble!(problem::Problem{Mortar}, time::Real)
|
||||
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
|
||||
add!(problem.assembly.C2, lsdofs, lsdofs, De)
|
||||
add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
|
||||
add!(problem.assembly.g, lsdofs, ge)
|
||||
end
|
||||
|
||||
add!(problem.assembly.g, sdofs, ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
end
|
||||
|
||||
function project_vertex_to_auxiliary_plane(p::Vector, x0::Vector, n0::Vector)
|
||||
return p - dot(p-x0, n0)*n0
|
||||
end
|
||||
|
||||
function inv3(P::Matrix)
|
||||
n, m = size(P)
|
||||
@assert n == m == 3
|
||||
a, b, c, d, e, f, g, h, i = P
|
||||
A = e*i - f*h
|
||||
B = -d*i + f*g
|
||||
C = d*h - e*g
|
||||
D = -b*i + c*h
|
||||
E = a*i - c*g
|
||||
F = -a*h + b*g
|
||||
G = b*f - c*e
|
||||
H = -a*f + c*d
|
||||
I = a*e - b*d
|
||||
return 1/(a*A + b*B + c*C)*[A B C; D E F; G H I]
|
||||
end
|
||||
|
||||
function vertex_inside_polygon(q, P; atol=1.0e-6)
|
||||
N = length(P)
|
||||
angle = 0.0
|
||||
for i=1:N
|
||||
A = P[i] - q
|
||||
B = P[mod(i,N)+1] - q
|
||||
c = norm(A)*norm(B)
|
||||
isapprox(c, 0.0; atol=atol) && return true
|
||||
cosa = dot(A,B)/c
|
||||
isapprox(cosa, 1.0; atol=atol) && return false
|
||||
isapprox(cosa, -1.0; atol=atol) && return true
|
||||
try
|
||||
angle += acos(cosa)
|
||||
catch
|
||||
info("Unable to calculate acos($(ForwardDiff.get_value(cosa))) when determining is a vertex inside polygon.")
|
||||
info("Polygon is: $(ForwardDiff.get_value(P)) and vertex under consideration is $(ForwardDiff.get_value(q))")
|
||||
info("Polygon corner point in loop: A=$(ForwardDiff.get_value(A)), B=$(ForwardDiff.get_value(B))")
|
||||
info("c = ||A||*||B|| = $(ForwardDiff.get_value(c))")
|
||||
rethrow()
|
||||
end
|
||||
end
|
||||
return isapprox(angle, 2*pi; atol=atol)
|
||||
end
|
||||
|
||||
function calculate_centroid(P)
|
||||
N = length(P)
|
||||
P0 = P[1]
|
||||
areas = [norm(1/2*cross(P[i]-P0, P[mod(i,N)+1]-P0)) for i=2:N]
|
||||
centroids = [1/3*(P0+P[i]+P[mod(i,N)+1]) for i=2:N]
|
||||
C = 1/sum(areas)*sum(areas.*centroids)
|
||||
return C
|
||||
end
|
||||
|
||||
function get_cells(P, C)
|
||||
N = length(P)
|
||||
cells = Vector[]
|
||||
# shared edge etc.
|
||||
N < 3 && return cells
|
||||
# trivial case, polygon already triangle / quadrangle
|
||||
#N == 3 && return Vector[P]
|
||||
#N == 4 && return Vector[P]
|
||||
#V = sum([cross(P[i], P[mod(i,N)+1]) for i=1:N])
|
||||
#A = 1/2*abs(dot(n, V))
|
||||
#info("A = $A")
|
||||
cells = Vector[Vector[C, P[i], P[mod(i,N)+1]] for i=1:N]
|
||||
return cells
|
||||
|
||||
maxa = 0.0
|
||||
maxj = 0
|
||||
for i=1:N
|
||||
A = P[i] - C
|
||||
B = P[mod(i,N)+1] - C
|
||||
theta = acos(dot(A,B)/(norm(A)*norm(B)))
|
||||
if theta > maxa
|
||||
maxa = theta
|
||||
maxj = i
|
||||
end
|
||||
end
|
||||
info("max angle $(maxa/pi*180) at index $maxj, N=$N")
|
||||
indices = mod(collect(maxj:maxj+N), N)
|
||||
info("indices = $indices")
|
||||
end
|
||||
|
||||
function get_polygon_clip(xs, xm, n; debug=false)
|
||||
# objective: search does line xm1 - xm2 clip xs
|
||||
nm = length(xm)
|
||||
ns = length(xs)
|
||||
P = Vector{Float64}[]
|
||||
|
||||
# 1. test is master point inside slave, if yes, add to clip
|
||||
for i=1:nm
|
||||
if vertex_inside_polygon(xm[i], xs)
|
||||
debug && info("1. $(xm[i]) inside S -> push")
|
||||
push!(P, xm[i])
|
||||
end
|
||||
end
|
||||
|
||||
# 2. test is slave point inside master, if yes, add to clip
|
||||
for i=1:ns
|
||||
if vertex_inside_polygon(xs[i], xm)
|
||||
xs[i] in P && continue
|
||||
debug && info("2. $(xs[i]) inside M -> push")
|
||||
push!(P, xs[i])
|
||||
end
|
||||
end
|
||||
|
||||
for i=1:nm
|
||||
# 2. find possible intersection
|
||||
xm1 = xm[i]
|
||||
xm2 = xm[mod(i,nm)+1]
|
||||
#info("intersecting line $xm1 -> $xm2")
|
||||
for j=1:ns
|
||||
xs1 = xs[j]
|
||||
xs2 = xs[mod(j,ns)+1]
|
||||
#info("clipping polygon edge $xs1 -> $xs2")
|
||||
tnom = dot(cross(xm1-xs1, xm2-xm1), n)
|
||||
tdenom = dot(cross(xs2-xs1, xm2-xm1), n)
|
||||
isapprox(tdenom, 0) && continue
|
||||
t = tnom/tdenom
|
||||
(0 <= t <= 1) || continue
|
||||
q = xs1 + t*(xs2 - xs1)
|
||||
#info("t=$t, q=$q, q ∈ xm ? $(vertex_inside_polygon(q, xm))")
|
||||
if vertex_inside_polygon(q, xm)
|
||||
q in P && continue
|
||||
debug && info("3. $q inside M -> push")
|
||||
push!(P, q)
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
return P
|
||||
end
|
||||
|
||||
function project_vertex_to_surface{E}(p::Vector, x0::Vector, n0::Vector,
|
||||
element::Element{E}, x::DVTI, time::Real; max_iterations::Int=10, iter_tol::Float64=1.0e-9)
|
||||
basis(xi) = get_basis(element, xi, time)
|
||||
dbasis(xi) = get_dbasis(element, xi, time)
|
||||
f(theta) = basis(theta[1:2])*x - theta[3]*n0 - p
|
||||
L(theta) = inv3([dbasis(theta[1:2])*x -n0])
|
||||
# L2(theta) = inv(ForwardDiff.get_value([dbasis(theta[2:3])*x -n0]))
|
||||
# FIXME: for some reason forwarddiff gives NaN's here.
|
||||
theta = zeros(3)
|
||||
dtheta = zeros(3)
|
||||
for i=1:max_iterations
|
||||
dtheta = L(theta) * f(theta)
|
||||
theta -= dtheta
|
||||
if norm(dtheta) < iter_tol
|
||||
return theta[1:2], theta[3]
|
||||
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 = $x0, n0 = $n0")
|
||||
info("element geometry: $(x.data)")
|
||||
info("vertex to project: $p")
|
||||
info("parameter vector before giving up: $theta")
|
||||
info("increment in parameter vector before giving up: $dtheta")
|
||||
info("norm(dtheta) before giving up: $(norm(dtheta))")
|
||||
info("f([0.0, 0.0, 0.0]) = $(f([0.0, 0.0, 0.0]))")
|
||||
info("L([0.0, 0.0, 0.0]) = $(L([0.0, 0.0, 0.0]))")
|
||||
|
||||
info("iterations:")
|
||||
theta = zeros(3)
|
||||
dtheta = zeros(3)
|
||||
for i=1:max_iterations
|
||||
info("iter $i, theta = $theta")
|
||||
info("f = $(f(theta))")
|
||||
info("L = $(L(theta))")
|
||||
dtheta = L(theta) * f(theta)
|
||||
info("dtheta = $(dtheta)")
|
||||
theta -= dtheta
|
||||
end
|
||||
|
||||
error("project_point_to_surface: did not converge in $max_iterations iterations!")
|
||||
end
|
||||
|
||||
function calculate_normals(elements, time, ::Type{Val{2}}; rotate_normals=false)
|
||||
normals = Dict{Int64, Vector{Float64}}()
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
J = transpose(element([0.0, 0.0], time, Val{:Jacobian}))
|
||||
normal = cross(J[:,1], J[:,2])
|
||||
for nid in conn
|
||||
if haskey(normals, nid)
|
||||
normals[nid] += normal
|
||||
else
|
||||
normals[nid] = normal
|
||||
end
|
||||
end
|
||||
end
|
||||
# normalize to unit normal
|
||||
S = collect(keys(normals))
|
||||
for j in S
|
||||
normals[j] /= norm(normals[j])
|
||||
end
|
||||
if rotate_normals
|
||||
for j in S
|
||||
normals[j] = -normals[j]
|
||||
end
|
||||
end
|
||||
return normals
|
||||
end
|
||||
|
||||
function check_orientation!(P, n; debug=false)
|
||||
C = mean(P)
|
||||
np = length(P)
|
||||
s = [dot(n, cross(P[i]-C, P[mod(i+1,np)+1]-C)) for i=1:np]
|
||||
all(s .< 0) && return
|
||||
debug && info("polygon not in ccw order, fixing")
|
||||
# project points to new orthogonal basis Q and sort there
|
||||
t1 = (P[1]-C)/norm(P[1]-C)
|
||||
t2 = cross(n, t1)
|
||||
Q = [n t1 t2]
|
||||
sort!(P, lt=(A, B) -> begin
|
||||
A_proj = Q'*(A-C)
|
||||
B_proj = Q'*(B-C)
|
||||
a = atan2(A_proj[3], A_proj[2])
|
||||
b = atan2(B_proj[3], B_proj[2])
|
||||
return a > b
|
||||
end)
|
||||
end
|
||||
|
||||
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=true)
|
||||
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_elements = get_slave_elements(problem)
|
||||
area = 0.0
|
||||
|
||||
# 1. calculate nodal normals and tangents for slave element nodes j ∈ S
|
||||
normals = calculate_normals(slave_elements, time, Val{2};
|
||||
rotate_normals=props.rotate_normals)
|
||||
update!(slave_elements, "normal", normals)
|
||||
|
||||
# 2. loop all slave elements
|
||||
for slave_element in slave_elements
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
nsl = length(slave_element)
|
||||
X1 = slave_element["geometry"](time)
|
||||
n1 = Field([normals[j] for j in slave_element_nodes])
|
||||
|
||||
# project slave nodes to auxiliary plane (x0, Q)
|
||||
#xi = get_reference_element_midpoint(slave_element)
|
||||
xi = [1/3, 1/3]
|
||||
N = vec(get_basis(slave_element, xi, time))
|
||||
x0 = N*X1
|
||||
n0 = N*n1
|
||||
S = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X1]
|
||||
|
||||
# 3. loop all master elements
|
||||
for master_element in slave_element["master elements"](time)
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
nm = length(master_element)
|
||||
X2 = master_element["geometry"](time)
|
||||
|
||||
# 3.1 project master nodes to auxiliary plane and create polygon clipping
|
||||
M = Vector[project_vertex_to_auxiliary_plane(p, x0, n0) for p in X2]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
length(P) < 3 && continue # no clipping or shared edge (no volume)
|
||||
check_orientation!(P, n0)
|
||||
C0 = calculate_centroid(P)
|
||||
|
||||
De = zeros(nsl, nsl)
|
||||
Me = zeros(nsl, nm)
|
||||
ge = zeros(field_dim*nsl)
|
||||
|
||||
# 4. loop integration cells
|
||||
for cell in get_cells(P, C0)
|
||||
virtual_element = Element(Tri3)
|
||||
update!(virtual_element, "geometry", cell)
|
||||
#x_cell = Field(cell)
|
||||
|
||||
# 5. loop integration point of integration cell
|
||||
for ip in get_integration_points(virtual_element, 3)
|
||||
N = vec(get_basis(virtual_element, ip, time))
|
||||
#dN = vec(get_dbasis(virtual_element, ip, time))
|
||||
#JC = transpose(sum([kron(dNC[:,j], x_cell[j]') for j=1:length(x_cell)]))
|
||||
#wC = ip.weight*norm(cross(JC[:,1], JC[:,2]))
|
||||
detJ = virtual_element(ip, time, Val{:detJ})
|
||||
w = ip.weight*detJ
|
||||
|
||||
# project gauss point from auxiliary plane to master and slave element
|
||||
#x_gauss = N*x_cell
|
||||
x_gauss = virtual_element("geometry", ip, time)
|
||||
if isnan(x_gauss[1])
|
||||
info("is nan")
|
||||
info("x_gauss = $x_gauss")
|
||||
info("cell = $cell")
|
||||
info("C0 = $C0")
|
||||
info("P = $P")
|
||||
info("S = $S")
|
||||
info("M = $M")
|
||||
info("n0 = $n0")
|
||||
error("nan, unable to continue")
|
||||
end
|
||||
xi_s, alpha = project_vertex_to_surface(x_gauss, x0, n0, slave_element, X1, time)
|
||||
xi_m, alpha = project_vertex_to_surface(x_gauss, x0, n0, master_element, X2, time)
|
||||
|
||||
# add contributions
|
||||
N1 = vec(get_basis(slave_element, xi_s, time))
|
||||
N2 = vec(get_basis(master_element, xi_m, time))
|
||||
De += w*N1*N1'
|
||||
Me += w*N1*N2'
|
||||
if props.adjust
|
||||
u1 = slave_element["displacement"](time)
|
||||
u2 = master_element["displacement"](time)
|
||||
x_s = N1*(X1+u1)
|
||||
x_m = N2*(X2+u2)
|
||||
ge += w*vec((x_m-x_s)*N1')
|
||||
end
|
||||
area += w
|
||||
end # integration points done
|
||||
|
||||
end # integration cells done
|
||||
|
||||
# 6. add contribution to contact virtual work
|
||||
sdofs = get_gdofs(problem, slave_element)
|
||||
mdofs = get_gdofs(problem, master_element)
|
||||
|
||||
for i=1:field_dim
|
||||
lsdofs = sdofs[i:field_dim:end]
|
||||
lmdofs = mdofs[i:field_dim:end]
|
||||
add!(problem.assembly.C1, lsdofs, lsdofs, De)
|
||||
add!(problem.assembly.C1, lsdofs, lmdofs, -Me)
|
||||
add!(problem.assembly.C2, lsdofs, lsdofs, De)
|
||||
add!(problem.assembly.C2, lsdofs, lmdofs, -Me)
|
||||
end
|
||||
add!(problem.assembly.g, sdofs, ge)
|
||||
|
||||
end # master elements done
|
||||
|
||||
end # slave elements done, contact virtual work ready
|
||||
|
||||
debug && info("area of interface: $area")
|
||||
|
||||
end
|
||||
|
||||
|
||||
@@ -31,3 +31,22 @@ function calc_nodal_values!(elements, field_name, field_dim, time)
|
||||
end
|
||||
update!(elements, field_name, nodal_values)
|
||||
end
|
||||
|
||||
"""
|
||||
Return node ids + vector of values
|
||||
"""
|
||||
function get_nodal_vector(elements, field_name, time)
|
||||
f = Dict{Int64, Vector{Float64}}()
|
||||
for element in elements
|
||||
for (c, v) in zip(get_connectivity(element), element[field_name](time))
|
||||
if haskey(f, c)
|
||||
@assert isapprox(f[c], v)
|
||||
end
|
||||
f[c] = v
|
||||
end
|
||||
end
|
||||
node_ids = sort(collect(keys(f)))
|
||||
field = [f[nid] for nid in node_ids]
|
||||
return node_ids, field
|
||||
end
|
||||
|
||||
|
||||
@@ -380,6 +380,7 @@ function aster_read_mesh(fn::ASCIIString, mesh_name=nothing)
|
||||
add_node!(mesh, nid, ncoords)
|
||||
end
|
||||
mapping = Dict(
|
||||
:PO1 => :Poi1,
|
||||
:SE2 => :Seg2,
|
||||
:TR3 => :Tri3,
|
||||
:TR6 => :Tri6,
|
||||
|
||||
+92
-55
@@ -216,80 +216,113 @@ function get_boundary_assembly(solver::Solver)
|
||||
end
|
||||
|
||||
|
||||
""" Solve linear system using LU factorization (UMFPACK).
|
||||
"""
|
||||
function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver_UMFPACK}})
|
||||
info("solving linear system of $(length(solver.problems)) problems.")
|
||||
t0 = time()
|
||||
Construct new basis such that u = P*uh + g
|
||||
|
||||
# assemble field problems
|
||||
M, K, Kg, f = get_field_assembly(solver)
|
||||
Parameters
|
||||
----------
|
||||
S set of linearly independent dofs.
|
||||
"""
|
||||
function create_projection(C::SparseMatrixCSC, g; S=nothing, tol=1.0e-12)
|
||||
n, m = size(C)
|
||||
@assert n == m
|
||||
if S == nothing
|
||||
S = get_nonzero_rows(C)
|
||||
end
|
||||
# FIXME: this creates dense matrices
|
||||
# efficiency / memory usage is a question
|
||||
P = sparse(C[S,:] \ full(C[S,:]))
|
||||
h = sparse(C[S,:] \ full(g[S]))
|
||||
resize!(P, n, m)
|
||||
resize!(h, n, 1)
|
||||
P = speye(n) - P
|
||||
SparseMatrix.droptol!(P, tol)
|
||||
return P, h
|
||||
end
|
||||
|
||||
# assemble boundary problems
|
||||
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
|
||||
|
||||
# construct global system Ax=b and solve using lu factorization
|
||||
A = [
|
||||
K+Kg+Kb C1'
|
||||
C2 D]
|
||||
b = [f+fb; g]
|
||||
"""
|
||||
Solve linear system using LDLt factorization (SuiteSparse). This version
|
||||
requires that final system is symmetric and positive definite, so boundary
|
||||
conditions are first eliminated before solution.
|
||||
"""
|
||||
function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{1}}; debug=false)
|
||||
|
||||
nnz(D) == 0 || return false
|
||||
nz = get_nonzero_rows(C2)
|
||||
B = get_nonzero_rows(C2')
|
||||
# C2^-1 exists or this doesn't work
|
||||
length(nz) == length(B) || return false
|
||||
|
||||
A = get_nonzero_rows(K)
|
||||
I = setdiff(A, B)
|
||||
|
||||
if debug
|
||||
info("# nz = $(length(nz))")
|
||||
info("# A = $(length(A))")
|
||||
info("# B = $(length(B))")
|
||||
info("# I = $(length(I))")
|
||||
end
|
||||
|
||||
# solver boundary dofs
|
||||
try
|
||||
u[B] = lufact(C2[nz,B]) \ full(g[nz])
|
||||
catch
|
||||
info("solver #1 failed to solve boundary dofs (you should not see this message).")
|
||||
return false
|
||||
end
|
||||
|
||||
# solve interior domain using LDLt factorization
|
||||
u[I] = ldltfact(K[I,I]) \ (f[I] - K[I,B]*u[B])
|
||||
# solve lambda
|
||||
la[B] = lufact(C1[B,nz]) \ full(f[B] - K[B,I]*u[I] - K[B,B]*u[B])
|
||||
|
||||
return true
|
||||
end
|
||||
|
||||
"""
|
||||
Solve linear system using LU factorization (UMFPACK). This version solves
|
||||
directly the saddle point problem without elimination of boundary conditions.
|
||||
"""
|
||||
function solve!(K, C1, C2, D, f, g, u, la, ::Type{Val{2}})
|
||||
# construct global system Ax = b and solve using lufact (UMFPACK)
|
||||
A = [K C1'; C2 D]
|
||||
b = [f; g]
|
||||
nz = get_nonzero_rows(A)
|
||||
x = zeros(length(b))
|
||||
x[nz] = lufact(A[nz,nz]) \ full(b[nz])
|
||||
|
||||
ndofs = solver.ndofs
|
||||
u = x[1:ndofs]
|
||||
la = x[ndofs+1:end]
|
||||
info("UMFPACK: solved in ", time()-t0, " seconds. norm = ", norm(u))
|
||||
return u, la
|
||||
ndofs = size(K, 1)
|
||||
u[:] = x[1:ndofs]
|
||||
la[:] = x[ndofs+1:end]
|
||||
return true
|
||||
end
|
||||
|
||||
""" Solve linear system using LDLt factorization (SuiteSparse). """
|
||||
function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
|
||||
function solve_linear_system(solver::Solver)
|
||||
info("solving linear system of $(length(solver.problems)) problems.")
|
||||
t0 = time()
|
||||
|
||||
# assemble field problems
|
||||
M, K, Kg, f = get_field_assembly(solver)
|
||||
|
||||
# assemble boundary problems
|
||||
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
|
||||
|
||||
K = K + Kb + Kg
|
||||
f = f + fb
|
||||
K = K + Kg + Kb
|
||||
K = 1/2*(K + K')
|
||||
f = f + fb
|
||||
|
||||
u = zeros(solver.ndofs)
|
||||
la = zeros(solver.ndofs)
|
||||
|
||||
# determine interior and boundary dofs
|
||||
all_dofs = get_nonzero_rows(K)
|
||||
boundary_dofs = get_nonzero_rows(C1)
|
||||
boundary_dofs2 = get_nonzero_rows(C2)
|
||||
interior_dofs = setdiff(all_dofs, boundary_dofs)
|
||||
@assert length(boundary_dofs) == length(boundary_dofs2)
|
||||
@assert setdiff(Set(boundary_dofs), Set(boundary_dofs2)) == Set()
|
||||
|
||||
# solve boundary
|
||||
LUF = lufact(C1[boundary_dofs, boundary_dofs])
|
||||
u[boundary_dofs] = LUF \ full(g[boundary_dofs])
|
||||
normub = norm(u[boundary_dofs])
|
||||
if isapprox(normub, 0.0)
|
||||
info("CHOLMOD: homogeneous dirichlet boundary condition.")
|
||||
status = false
|
||||
for i in [1, 2]
|
||||
status = solve!(K, C1, C2, D, f, g, u, la, Val{i})
|
||||
if status
|
||||
info("succesfully solved Ax = b using solver #$i")
|
||||
break
|
||||
end
|
||||
end
|
||||
status || error("Failed to solve linear system!")
|
||||
|
||||
# solver interior
|
||||
CF = ldltfact(K[interior_dofs, interior_dofs])
|
||||
Kib = K[interior_dofs, boundary_dofs]
|
||||
Kbb = K[boundary_dofs, boundary_dofs]
|
||||
fi = f[interior_dofs]
|
||||
u[interior_dofs] = CF \ (fi - Kib*u[boundary_dofs])
|
||||
|
||||
# solve lambda
|
||||
# LUF2 = lufact(C2[boundary_dofs, boundary_dofs])
|
||||
la[boundary_dofs] = LUF \ full(Kib' * u[interior_dofs] - Kbb*u[boundary_dofs])
|
||||
|
||||
info("CHOLMOD: solved in ", time()-t0, " seconds. norm = ", norm(u))
|
||||
info("linear system solver: solved in ", time()-t0, " seconds. norm = ", norm(u))
|
||||
return u, la
|
||||
end
|
||||
|
||||
@@ -350,15 +383,19 @@ function assemble!(solver::Solver; force_assembly=true)
|
||||
info("Assembled in $t1 seconds.")
|
||||
end
|
||||
|
||||
function initialize!(solver::Solver)
|
||||
for problem in solver.problems
|
||||
initialize!(problem, solver.time)
|
||||
end
|
||||
end
|
||||
|
||||
""" Default solver for quasistatic nonlinear problems. """
|
||||
function call(solver::Solver{Nonlinear})
|
||||
|
||||
properties = solver.properties
|
||||
|
||||
# 1. initialize each problem so that we can start nonlinear iterations
|
||||
for problem in solver.problems
|
||||
initialize!(problem, solver.time)
|
||||
end
|
||||
initialize!(solver)
|
||||
|
||||
# 2. start non-linear iterations
|
||||
for properties.iteration=1:properties.max_iterations
|
||||
@@ -370,7 +407,7 @@ function call(solver::Solver{Nonlinear})
|
||||
# 2.2 call solver for linearized system (default: direct lu factorization)
|
||||
info("Solve linear system ...")
|
||||
tic()
|
||||
u, la = solve_linear_system(solver, Val{properties.linear_system_solver})
|
||||
u, la = solve_linear_system(solver)
|
||||
push!(properties.norms, (norm(u), norm(la)))
|
||||
t1 = round(toq(), 2)
|
||||
info("Solved Ax = b in $t1 seconds.")
|
||||
|
||||
+7
-4
@@ -147,16 +147,19 @@ Returns
|
||||
Ordered list of row indices.
|
||||
"""
|
||||
function get_nonzero_rows(A::SparseMatrixCSC)
|
||||
# FIXME: This is probably a very inefficient way to do this.
|
||||
return sort(unique(rowvals(A)))
|
||||
end
|
||||
|
||||
function get_nonzero_rows(A::SparseMatrixCOO)
|
||||
function get_nonzero_columns(A::SparseMatrixCSC)
|
||||
return get_nonzero_rows(transpose(A))
|
||||
end
|
||||
|
||||
function get_nonzero_rows(A::Union{SparseMatrixCOO, Matrix})
|
||||
return get_nonzero_rows(sparse(A))
|
||||
end
|
||||
|
||||
function get_nonzero_rows(A::Matrix)
|
||||
return get_nonzero_rows(sparse(A))
|
||||
function get_nonzero_columns(A::Union{SparseMatrixCOO, Matrix})
|
||||
return get_nonzero_columns(sparse(A))
|
||||
end
|
||||
|
||||
function size(A::SparseMatrixCOO)
|
||||
|
||||
@@ -18,7 +18,7 @@ using JuliaFEM.Test
|
||||
block.elements = create_elements(mesh, "BLOCK")
|
||||
update!(block.elements, "youngs modulus", 288.0)
|
||||
update!(block.elements, "poissons ratio", 1/3)
|
||||
# update!(block.elements, "displacement load 2", 576.0)
|
||||
update!(block.elements, "displacement load 2", 576.0)
|
||||
|
||||
traction = create_elements(mesh, "TOP")
|
||||
update!(traction, "displacement traction force 2", 288.0)
|
||||
|
||||
@@ -42,3 +42,79 @@ using JuliaFEM.Test
|
||||
call(s1; debug=true)
|
||||
@test isapprox(s1.properties.eigvals, [5/3, 2/3])
|
||||
end
|
||||
|
||||
@testset "test poisson problem modal analysis without tie" begin
|
||||
X = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 3.0],
|
||||
4 => [0.0, 3.0],
|
||||
5 => [0.0, 3.0],
|
||||
6 => [1.0, 3.0],
|
||||
7 => [1.0, 9.0],
|
||||
8 => [0.0, 9.0])
|
||||
T = Dict{Int64, Float64}()
|
||||
for i=1:8
|
||||
T[i] = 0.0
|
||||
end
|
||||
el1 = Element(Quad4, [1, 2, 3, 4])
|
||||
el2 = Element(Quad4, [4, 3, 7, 8])
|
||||
el3 = Element(Seg2, [1, 2])
|
||||
el4 = Element(Seg2, [7, 8])
|
||||
update!([el1, el2, el3, el4], "geometry", X)
|
||||
update!([el1, el2], "density", 6.0)
|
||||
update!([el1, el2], "temperature thermal conductivity", 36.0)
|
||||
update!([el1, el2], "temperature", T)
|
||||
update!([el3, el4], "temperature 1", 0.0)
|
||||
p1 = Problem(Heat, "combined body", 1)
|
||||
p2 = Problem(Dirichlet, "fixed ends", 1, "temperature")
|
||||
push!(p1, el1, el2)
|
||||
push!(p2, el3, el4)
|
||||
|
||||
solver = Solver(Modal)
|
||||
push!(solver, p1, p2)
|
||||
call(solver)
|
||||
@test isapprox(solver.properties.eigvals[1], 1.0)
|
||||
end
|
||||
|
||||
@testset "test poisson modal problem with mesh tie" begin
|
||||
X = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 3.0],
|
||||
4 => [0.0, 3.0],
|
||||
5 => [0.0, 3.0],
|
||||
6 => [1.0, 3.0],
|
||||
7 => [1.0, 9.0],
|
||||
8 => [0.0, 9.0])
|
||||
T = Dict{Int64, Float64}()
|
||||
for i=1:8
|
||||
T[i] = 0.0
|
||||
end
|
||||
el1 = Element(Quad4, [1, 2, 3, 4])
|
||||
el2 = Element(Quad4, [5, 6, 7, 8])
|
||||
el3 = Element(Seg2, [1, 2])
|
||||
el4 = Element(Seg2, [7, 8])
|
||||
el5 = Element(Seg2, [3, 4])
|
||||
el6 = Element(Seg2, [5, 6])
|
||||
update!([el1, el2, el3, el4, el5, el6], "geometry", X)
|
||||
update!([el1, el2], "temperature", T)
|
||||
update!([el1, el2], "density", 6.0)
|
||||
update!([el1, el2], "temperature thermal conductivity", 36.0)
|
||||
update!([el3, el4], "temperature 1", 0.0)
|
||||
update!(el5, "master elements", [el6])
|
||||
p1 = Problem(Heat, "body 1", 1)
|
||||
p2 = Problem(Heat, "body 2", 1)
|
||||
p3 = Problem(Dirichlet, "fixed ends", 1, "temperature")
|
||||
p4 = Problem(Mortar, "interface between bodies", 1, "temperature")
|
||||
p4.properties.dimension = 1
|
||||
push!(p1, el1)
|
||||
push!(p2, el2)
|
||||
push!(p3, el3, el4)
|
||||
push!(p4, el5, el6)
|
||||
solver = Solver(Modal)
|
||||
push!(solver, p1, p2, p3, p4)
|
||||
call(solver)
|
||||
@test isapprox(solver.properties.eigvals[1], 1.0)
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,119 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Test
|
||||
|
||||
function get_test_model()
|
||||
X = Dict{Int64, Vector{Float64}}(
|
||||
1 => [0.0, 0.0],
|
||||
2 => [1.0, 0.0],
|
||||
3 => [1.0, 0.5],
|
||||
4 => [0.0, 0.5],
|
||||
5 => [0.0, 0.6],
|
||||
6 => [1.0, 0.6],
|
||||
7 => [1.0, 1.1],
|
||||
8 => [0.0, 1.1])
|
||||
el1 = Element(Quad4, [1, 2, 3, 4])
|
||||
el2 = Element(Quad4, [5, 6, 7, 8])
|
||||
el3 = Element(Seg2, [1, 2])
|
||||
el4 = Element(Seg2, [7, 8])
|
||||
el5 = Element(Seg2, [4, 3])
|
||||
el6 = Element(Seg2, [5, 6])
|
||||
update!([el1, el2, el3, el4, el5, el6], "geometry", X)
|
||||
update!([el1, el2], "youngs modulus", 96.0)
|
||||
update!([el1, el2], "poissons ratio", 1/3)
|
||||
update!([el3], "displacement 1", 0.0)
|
||||
update!([el3], "displacement 2", 0.0)
|
||||
update!([el4], "displacement 1", 0.0)
|
||||
update!([el4], "displacement 2", 0.0)
|
||||
update!(el6, "master elements", [el5])
|
||||
p1 = Problem(Elasticity, "body1", 2)
|
||||
p2 = Problem(Elasticity, "body2", 2)
|
||||
p3 = Problem(Dirichlet, "fixed", 2, "displacement")
|
||||
p4 = Problem(Mortar, "interface", 2, "displacement")
|
||||
push!(p1, el1)
|
||||
push!(p2, el2)
|
||||
push!(p3, el3, el4)
|
||||
push!(p4, el5, el6)
|
||||
return p1, p2, p3, p4
|
||||
end
|
||||
|
||||
@testset "test adjust setting in 2d tie contact" begin
|
||||
p1, p2, p3, p4 = get_test_model()
|
||||
p1.properties.formulation = :plane_stress
|
||||
p2.properties.formulation = :plane_stress
|
||||
p4.properties.dimension = 1
|
||||
p4.properties.adjust = true
|
||||
p4.properties.rotate_normals = false
|
||||
solver = Solver(Nonlinear)
|
||||
push!(solver, p1, p2, p3, p4)
|
||||
call(solver)
|
||||
el5 = p4.elements[1]
|
||||
u = el5("displacement", [0.0], 0.0)
|
||||
info("u = $u")
|
||||
@test isapprox(u, [0.0, 0.05])
|
||||
end
|
||||
|
||||
@testset "test that interface transfers constant field without error" begin
|
||||
meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_2d.med"
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Heat, "upper", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper.elements, "temperature thermal conductivity", 1.0)
|
||||
|
||||
lower = Problem(Heat, "lower", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower.elements, "temperature thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
|
||||
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
|
||||
update!(bc_upper.elements, "temperature 1", 0.0)
|
||||
|
||||
bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
|
||||
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
|
||||
update!(bc_lower.elements, "temperature 1", 1.0)
|
||||
|
||||
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
|
||||
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
interface.properties.dimension = 1
|
||||
|
||||
solver = Solver()
|
||||
push!(solver, upper, lower, bc_upper, bc_lower, interface)
|
||||
call(solver)
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
T = [t[1] for t in temperature]
|
||||
minT = minimum(T)
|
||||
maxT = maximum(T)
|
||||
info("minT = $minT, maxT = $maxT")
|
||||
@test isapprox(minT, 0.5)
|
||||
@test isapprox(maxT, 0.5)
|
||||
end
|
||||
|
||||
|
||||
#=
|
||||
|
||||
@testset "expect clear error when trying to solve 2d model in 3d setting" begin
|
||||
p1, p2, p3, p4 = get_test_model()
|
||||
# p1.properties.formulation = :plane_stress
|
||||
# p2.properties.formulation = :plane_stress
|
||||
p4.properties.adjust = true
|
||||
p4.properties.rotate_normals = false
|
||||
solver = Solver(Nonlinear)
|
||||
solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK
|
||||
push!(solver, p1, p2, p3, p4)
|
||||
call(solver)
|
||||
el5 = p4.elements[1]
|
||||
u = el5("displacement", [0.0], 0.0)
|
||||
info("u = $u")
|
||||
@test isapprox(u, [0.0, 0.05])
|
||||
end
|
||||
|
||||
=#
|
||||
@@ -0,0 +1,47 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Preprocess
|
||||
using JuliaFEM.Postprocess
|
||||
using JuliaFEM.Test
|
||||
|
||||
@testset "test that interface transfers constant field without error" begin
|
||||
meshfile = Pkg.dir("JuliaFEM") * "/test/testdata/block_3d.med"
|
||||
mesh = aster_read_mesh(meshfile)
|
||||
|
||||
upper = Problem(Heat, "upper", 1)
|
||||
upper.elements = create_elements(mesh, "UPPER")
|
||||
update!(upper.elements, "temperature thermal conductivity", 1.0)
|
||||
lower = Problem(Heat, "lower", 1)
|
||||
lower.elements = create_elements(mesh, "LOWER")
|
||||
update!(lower.elements, "temperature thermal conductivity", 1.0)
|
||||
|
||||
bc_upper = Problem(Dirichlet, "upper boundary", 1, "temperature")
|
||||
bc_upper.elements = create_elements(mesh, "UPPER_TOP")
|
||||
update!(bc_upper.elements, "temperature 1", 0.0)
|
||||
|
||||
bc_lower = Problem(Dirichlet, "lower boundary", 1, "temperature")
|
||||
bc_lower.elements = create_elements(mesh, "LOWER_BOTTOM")
|
||||
update!(bc_lower.elements, "temperature 1", 1.0)
|
||||
|
||||
interface = Problem(Mortar, "interface between upper and lower block", 1, "temperature")
|
||||
interface_slave_elements = create_elements(mesh, "LOWER_TOP")
|
||||
interface_master_elements = create_elements(mesh, "UPPER_BOTTOM")
|
||||
update!(interface_slave_elements, "master elements", interface_master_elements)
|
||||
interface.elements = [interface_master_elements; interface_slave_elements]
|
||||
interface.properties.dimension = 2
|
||||
|
||||
solver = Solver()
|
||||
solver.properties.linear_system_solver = :DirectLinearSolver_UMFPACK
|
||||
push!(solver, upper, lower, bc_upper, bc_lower, interface)
|
||||
call(solver)
|
||||
|
||||
node_ids, temperature = get_nodal_vector(interface.elements, "temperature", 0.0)
|
||||
T = [t[1] for t in temperature]
|
||||
minT = minimum(T)
|
||||
maxT = maximum(T)
|
||||
info("minT = $minT, maxT = $maxT")
|
||||
@test isapprox(minT, 0.5)
|
||||
@test isapprox(maxT, 0.5)
|
||||
end
|
||||
@@ -0,0 +1,41 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
@testset "polygon clip case 1" begin
|
||||
S = Vector[
|
||||
[0.375, 0.0, 0.5],
|
||||
[0.6, 0.0, 0.5],
|
||||
[0.5, 0.25, 0.5]]
|
||||
M = Vector[
|
||||
[0.50, 0.0, 0.5],
|
||||
[0.25, 0.0, 0.5],
|
||||
[0.375, 0.25, 0.5]]
|
||||
n0 = [0.0, 0.0, 1.0]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
P_expected = Vector{Float64}[
|
||||
[0.500, 0.0, 0.5],
|
||||
[0.375, 0.0, 0.5],
|
||||
[0.4375, 0.125, 0.5]]
|
||||
@test length(P) == length(P_expected)
|
||||
for (Pi, Pj) in zip(P, P_expected)
|
||||
@test isapprox(Pi, Pj)
|
||||
end
|
||||
end
|
||||
|
||||
@testset "polygon clip case 2" begin
|
||||
S = Vector[
|
||||
[0.25, 0.0, 0.5],
|
||||
[0.75, 0.0, 0.5],
|
||||
[0.50, 0.25, 0.5]]
|
||||
M = Vector[
|
||||
[0.50, 0.0, 0.5],
|
||||
[0.25, 0.0, 0.5],
|
||||
[0.375, 0.25, 0.5]]
|
||||
n0 = [0.0, 0.0, 1.0]
|
||||
P = get_polygon_clip(S, M, n0)
|
||||
@test length(P) == 3
|
||||
end
|
||||
|
||||
@@ -0,0 +1,32 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Test
|
||||
|
||||
@testset "test projection" begin
|
||||
C = [
|
||||
2.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
1.0 2.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 2.0 1.0 -1.0 -2.0 0.0 0.0
|
||||
0.0 0.0 1.0 2.0 -2.0 -1.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
|
||||
g = [3.0, 3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
P, h = create_projection(sparse(C), g)
|
||||
P_expected = [
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0
|
||||
0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0]
|
||||
h_expected = [1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
|
||||
@test isapprox(full(P), P_expected)
|
||||
@test isapprox(full(h), h_expected)
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user