Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-14 16:14:17 +00:00
Code Issues Packages Projects Releases Wiki Activity

all tests pass now

Browse Source
This commit is contained in:
Jukka Aho
2016-07-03 21:16:03 +03:00
parent eec6db657d
commit 757ba2b3d9
52 changed files with 607 additions and 853 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/JuliaFEM.jl
+13 -19
View File
@@ -7,17 +7,9 @@ This is JuliaFEM -- Finite Element Package
module JuliaFEM
importall Base
using ForwardDiff
using JLD
#Grad = Val{:Grad}
#detJ = Val{:detJ}
#export Grad, detJ
include("common.jl")
include("fields.jl")
export Field, DCTI, DVTI, DCTV, DVTV
export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV
include("types.jl") # data types: Point, IntegrationPoint, ...
export AbstractPoint, Point, IntegrationPoint, IP, Node
#include("basis.jl") # interpolation of discrete fields
@@ -26,15 +18,15 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
### ELEMENTS ###
include("elements.jl") # common element routines
export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis
include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
include("lagrange.jl") # Continuous Galerkin (Lagrange) elements
include("elements_lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
export get_reference_coordinates
export Poi1,
Seg2, Seg3,
Tri3, Tri6, Quad4, Quad8, Quad9,
Tet4, Tet10, Hex8, Hex20, Hex27
include("nurbs.jl")
include("elements_nurbs.jl")
export NSeg, NSurf, NSolid, is_nurbs
#include("hierarchical.jl") # P-elements
@@ -50,13 +42,13 @@ export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
get_unknown_field_dimension, get_gdofs, Assembly,
get_parent_field_name, get_elements
include("elasticity.jl")
include("problems_elasticity.jl")
export Elasticity
include("dirichlet.jl")
include("problems_dirichlet.jl")
export Dirichlet
include("heat.jl")
include("problems_heat.jl")
export Heat
export assemble!, postprocess!
@@ -74,22 +66,24 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
get_field_problems, get_boundary_problems,
get_field_assembly, get_boundary_assembly,
initialize!, create_projection, eliminate_interior_dofs
include("modal.jl")
include("solvers_modal.jl")
export Modal
include("optics.jl")
export find_intersection, calc_reflection, calc_normal
### Mortar methods ###
include("mortar.jl")
include("problems_mortar.jl")
include("problems_mortar_2d_autodiff.jl")
export calculate_normals,
calculate_normals!,
project_from_slave_to_master,
project_from_master_to_slave,
Mortar, get_slave_elements
Mortar, get_slave_elements,
get_polygon_clip
### Mortar methods, contact mechanics extension ###
include("contact.jl")
include("problems_contact.jl")
export Contact
# rest of things
src/assembly.jl
+25
View File
@@ -26,6 +26,31 @@ function append!(assembly::Assembly, sub_assembly::Assembly)
append!(assembly.c, sub_assembly.c)
end
""" Calculate norm of assembly, i.e., norm of each block of matrix. """
function norm(assembly::Assembly, p=2)
N1 = norm(assembly.M, p)
N2 = norm(assembly.K, p)
N3 = norm(assembly.Kg, p)
N4 = norm(assembly.f, p)
N5 = norm(assembly.fg, p)
N6 = norm(assembly.C1, p)
N7 = norm(assembly.C2, p)
N8 = norm(assembly.D, p)
N9 = norm(assembly.g, p)
N10 = norm(assembly.c, p)
return [N1, N2, N3, N4, N5, N6, N7, N8, N9, N10]
end
function isapprox(a1::Assembly, a2::Assembly)
T = isapprox(a1.K, a2.K)
T &= isapprox(a1.C1, a2.C1)
T &= isapprox(a1.C2, a2.C2)
T &= isapprox(a1.D, a2.D)
T &= isapprox(a1.f, a2.f)
T &= isapprox(a1.g, a2.g)
return T
end
function assemble_prehook!
end
src/common.jl
-80
View File
@@ -1,83 +1,3 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
"""
A very simple debugging macro. It executes commands if environment variable DEBUG is set.
Usage
-----
Instead of starting session `julia file.jl`, do `DEBUG=1 julia file.jl`.
Or set `export DEBUG=1` for your `.bashrc`.
Running inside code
-------------------
julia> @debug info("moimoi heihei")
will get executed iff environment variable DEBUG is set.
Examples
--------
julia> @debug info("moimoi")
(empty)
julia> ENV["DEBUG"] = 1
julia> @debug info("moimoi")
INFO: moimoi
julia> @debug begin
... info("matrix is")
... dump([1 2; 3 4])
... end
INFO: matrix is
Array(Int64(2,2)) 2x2 Array{Int64,2}:
1 2
3 4
"""
macro debug(msg)
haskey(ENV, "DEBUG") || return
return msg
end
function set_debug_on!()
ENV["DEBUG"] = 1;
end
function set_debug_off!()
pop!(ENV, "DEBUG");
end
#=
""" Simple linspace extension to arrays.
Examples
--------
>>> linspace([0.0], [1.0], 3)
3-element Array{Array{Float64,1},1}:
[0.0]
[0.5]
[1.0]
"""
function linspace{T<:Array}(X1::T, X2::T, n)
[1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
end
=#
function resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
(n == A.n) && (m == A.m) && return
@assert n >= A.n
@assert m >= A.m
append!(A.colptr, A.colptr[end]*ones(Int, m-A.m))
A.n = n
A.m = m
end
function ForwardDiff.derivative{T}(f::Function, S::Matrix{T}, args...)
shape = size(S)
wrapper(S::Vector) = f(reshape(S, shape))
deriv = ForwardDiff.gradient(wrapper, vec(S), args...)
return reshape(deriv, shape)
end
export @debug, set_debug_on!, set_debug_off!
src/elements_lagrange.jl
+8
View File
@@ -26,6 +26,14 @@ function call(element::Element{Poi1}, ip, time, ::Type{Val{:detJ}})
return 1.0
end
function get_integration_order(element::Poi1)
return 1
end
function get_integration_points(element::Poi1, order::Int64)
return [ (1.0, [] ) ]
end
### 1d elements
type Seg2 <: AbstractElement
src/elements_nurbs.jl
+11
View File
@@ -1,6 +1,9 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using ForwardDiff
# TODO: evaluate partial derivatives of basis functions without forwarddiff
""" NURBS segment. """
type NSeg <: AbstractElement
order :: Int
@@ -98,6 +101,14 @@ function get_basis(element::Element{NSolid}, xi::Vector, time)
return N / sum(N)
end
# TODO: evaluate partial derivatives of basis functions without forwarddiff
""" Evaluate partial derivatives of basis functions using ForwardDiff. """
function get_dbasis{E<:Union{NSeg, NSurf, NSolid}}(element::Element{E}, ip, time)
xi = isa(ip, IP) ? ip.coords : ip
basis(xi) = vec(get_basis(element, xi, time))
return ForwardDiff.jacobian(basis, xi)'
end
function length(element::Element{NSeg})
nu = length(element.properties.knots) - element.properties.order - 1
return nu
src/preprocess_aster_reader.jl
-12
View File
@@ -128,28 +128,22 @@ end
function aster_renumber_nodes!(mesh1, mesh2)
reserved_node_ids = Set(collect(keys(mesh1["nodes"])))
@debug info("already reserved node ids: $reserved_node_ids")
mesh2_node_numbering = Dict{Int64, Int64}()
# find new node ids assigned for mesh 2
k = 1
for node_id in sort(collect(keys(mesh2["nodes"])))
@debug info("mesh2: processing node $node_id")
# if node id is reserved in mesh 1, find new number
if node_id in reserved_node_ids
@debug info("node id conflict, $node_id already defined in mesh 1, renumbering")
while k in reserved_node_ids
k += 1
end
@debug info("mesh2: node $node_id -> $k")
mesh2_node_numbering[node_id] = k
push!(reserved_node_ids, k)
else
mesh2_node_numbering[node_id] = node_id
end
end
@debug info("new node numering:")
@debug println(mesh2_node_numbering)
aster_renumber_nodes_!(mesh2, mesh2_node_numbering)
#=
@@ -175,28 +169,22 @@ end
function aster_renumber_elements!(mesh1, mesh2)
reserved_element_ids = Set(collect(keys(mesh1["connectivity"])))
@debug info("already reserved element ids: $reserved_element_ids")
mesh2_element_numbering = Dict{Int64, Int64}()
# find new element ids assigned for mesh 2
k = 1
for element_id in sort(collect(keys(mesh2["connectivity"])))
@debug info("mesh2: processing element $element_id")
# if node id is reserved in mesh 1, find new number
if element_id in reserved_element_ids
@debug info("element id conflict, $element_id already defined in mesh 1, renumbering")
while k in reserved_element_ids
k += 1
end
@debug info("mesh2: element $element_id -> $k")
mesh2_element_numbering[element_id] = k
push!(reserved_element_ids, k)
else
mesh2_element_numbering[element_id] = element_id
end
end
@debug info("element numbering for mesh 2:")
@debug info(mesh2_element_numbering)
# create new elements
mesh2_old_elements = mesh2["connectivity"]
src/problems.jl
+7 -1
View File
@@ -152,7 +152,13 @@ function initialize!(problem::Problem, time=0.0)
gdofs = get_gdofs(problem, element)
if haskey(element, field_name)
# if field is found, copy last known solution to new time as initial guess
if !isapprox(last(element[field_name]).time, time)
field = last(element[field_name])
if !isa(field, TimeVariantField)
info("Unable to initialize field $field_name for problem, is not time variant?")
continue
end
if !isapprox(field.time, time)
last_data = copy(last(element[field_name]).data)
push!(element[field_name], time => last_data)
end
src/dirichlet.jl → src/problems_dirichlet.jl
View File
src/problems_elasticity.jl
+1 -1
View File
@@ -181,7 +181,7 @@ function assemble{El<:Elasticity2DVolumeElements}(problem::Problem{Elasticity},
if haskey(element, "displacement load")
b = element("displacement load", ip, time)
f += w*vec(N'*b)
f += w*vec(b*N)
end
for i=1:dim
src/problems_mortar.jl
+13 -262
View File
@@ -165,14 +165,14 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
for slave_element in slave_elements
nsl = length(slave_element)
X1 = slave_element["geometry"](time)
n1 = slave_element["normal"](time)
X1 = slave_element("geometry", time)
n1 = slave_element("normal", time)
# 3. loop all master elements
for master_element in slave_element["master elements"](time)
for master_element in slave_element("master elements", time)
nm = length(master_element)
X2 = master_element["geometry"](time)
X2 = master_element("geometry", time)
# 3.1 calculate segmentation
xi1a = project_from_master_to_slave(slave_element, X2[1], time)
@@ -225,8 +225,8 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
haskey(master_element, "displacement") || continue
norm(mean(X1) - X2[1]) / norm(X1[2] - X1[1]) < props.distval || continue
norm(mean(X1) - X2[2]) / norm(X1[2] - X1[1]) < props.distval || continue
u1 = slave_element["displacement"](time)
u2 = master_element["displacement"](time)
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)*Phi')
@@ -253,256 +253,7 @@ function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Ty
end
# mesh tie 2d end
# mesh tie 2d forwarddiff start
function project_from_master_to_slave{E<:MortarElements2D}(
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector, time::Float64;
tol=1.0e-10, max_iterations=20)
x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
xi1 = 0.0
dxi1 = 0.0
for i=1:max_iterations
dxi1 = -R(xi1)/dR(xi1)
xi1 += dxi1
if norm(dxi1) < tol
return xi1
end
end
info("x1 = $(ForwardDiff.get_value(x1_.data))")
info("n1 = $(ForwardDiff.get_value(n1_.data))")
info("x2 = $(ForwardDiff.get_value(x2))")
info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
error("find projection from master to slave: did not converge")
end
function project_from_slave_to_master{E<:MortarElements2D}(
master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI, time::Float64;
tol=1.0e-10, max_iterations=20)
x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi2) = cross2(x2(xi2)-x1, n1)
dR(xi2) = cross2(dx2(xi2), n1)
xi2 = 0.0
dxi2 = 0.0
for i=1:max_iterations
dxi2 = -R(xi2) / dR(xi2)
xi2 += dxi2
if norm(dxi2) < tol
return xi2
end
end
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
end
""" 2d mesh tie using ForwardDiff.
Construct .. + fc*la and C(d,la)=0
"""
function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{true}})
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
slave_elements = get_slave_elements(problem)
if field_name != "displacement"
error("mortar forwarddiff assembly: only displacement field with adjust=yes supported")
end
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)
gap = zeros(u)
C = zeros(la)
S = Set{Int64}()
# 1. update nodal normals for slave elements
tangents = zeros(u)
for element in slave_elements
conn = get_connectivity(element)
push!(S, conn...)
X1 = element("geometry", time)
u1 = Field([u[:,i] for i in conn])
x1 = X1 + u1
dN = get_dbasis(element, [0.0], time)
tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
for nid in conn
tangents[:,nid] += tangent[:]
end
end
Q = [0.0 -1.0; 1.0 0.0]
normals = zeros(u)
for j in S
tangents[:,j] /= norm(tangents[:,j])
normals[:,j] = Q*tangents[:,j]
end
if props.rotate_normals
for j in S
normals[:,j] = -normals[:,j]
end
end
normals2 = Dict()
tangents2 = Dict()
for j in S
normals2[j] = normals[:,j]
tangents2[j] = tangents[:,j]
end
update!(slave_elements, "normal", time => normals2)
update!(slave_elements, "tangent", time => tangents2)
# 2. loop all slave elements
for slave_element in slave_elements
nsl = length(slave_element)
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])
# 3. loop all master elements
for master_element in slave_element["master elements"](time)
nm = length(master_element)
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
# 3.1 calculate segmentation
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1], time)
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2], time)
# xi1a = project_from_master_to_slave(slave_element, X2[1], time)
# xi1b = project_from_master_to_slave(slave_element, X2[2], time)
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
# 3.2. bi-orthogonal basis
De = zeros(nsl, nsl)
Me = zeros(nsl, nsl)
Ae = zeros(nsl, nsl)
if props.dual_basis
for ip in get_integration_points(slave_element, 3)
detJ = slave_element(ip, time, Val{:detJ})
w = ip.weight*detJ*l
xi = ip.coords[1]
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
N1 = vec(get_basis(slave_element, xi_s, time))
De += w*diagm(N1)
Me += w*N1*N1'
end
Ae = De*inv(Me)
else
Ae = eye(nsl)
end
# 3.3. loop integration points of one integration segment and calculate
# local mortar matrices
for ip in get_integration_points(slave_element, 3)
detJ = slave_element(ip, time, Val{:detJ})
w = ip.weight*detJ*l
#dN = get_dbasis(slave_element, ip, time)
#j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
#w = ip.weight*norm(j)*l
xi = ip.coords[1]
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
N1 = vec(get_basis(slave_element, xi_s, time))
Phi = Ae*N1
# project gauss point from slave element to master element in direction n_s
x_s = N1*x1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
#xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2, time)
N2 = vec(get_basis(master_element, xi_m, time))
x_m = N2*x2
la_s = Phi*la1
gn = dot(n_s, x_s-x_m)
u_s = N1*u1
u_m = N2*u2
X_s = N1*X1
X_m = N2*X2
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[:,slave_element_nodes] += w*(u_s-u_m)*Phi'
if props.adjust
G = ForwardDiff.get_value(w*(X_s-X_m)*Phi')
gap[:,slave_element_nodes] += G
end
end
end # master elements done
end # slave elements done, contact virtual work ready
C = gap
info("interface residual ready")
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]
empty!(problem.assembly)
problem.assembly.K = K
problem.assembly.C1 = C1
problem.assembly.C2 = C2
problem.assembly.D = D
problem.assembly.f = f
problem.assembly.g = g
end
## Mesh tie 2d end
## 3d Mortar mesh tie
@@ -730,7 +481,7 @@ function check_orientation!(P, n; debug=false)
end)
end
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=true)
function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}, ::Type{Val{false}}; debug=true)
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
@@ -748,7 +499,7 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=t
slave_element_nodes = get_connectivity(slave_element)
nsl = length(slave_element)
X1 = slave_element["geometry"](time)
X1 = slave_element("geometry", time)
n1 = Field([normals[j] for j in slave_element_nodes])
# project slave nodes to auxiliary plane (x0, Q)
@@ -760,11 +511,11 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=t
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)
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)
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]
@@ -815,8 +566,8 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}}; debug=t
De += w*N1*N1'
Me += w*N1*N2'
if props.adjust
u1 = slave_element["displacement"](time)
u2 = master_element["displacement"](time)
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')
src/problems_mortar_2d_autodiff.jl
+261
View File
@@ -1,6 +1,264 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using ForwardDiff
# forwarddiff version of mesh tying in 2d
function project_from_master_to_slave{E<:MortarElements2D}(
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector, time::Float64;
tol=1.0e-10, max_iterations=20)
x1(xi1) = vec(get_basis(slave_element, [xi1], time))*x1_
dx1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*x1_
n1(xi1) = vec(get_basis(slave_element, [xi1], time))*n1_
dn1(xi1) = vec(get_dbasis(slave_element, [xi1], time))*n1_
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi1) = cross2(x1(xi1)-x2, n1(xi1))
dR(xi1) = cross2(dx1(xi1), n1(xi1)) + cross2(x1(xi1)-x2, dn1(xi1))
xi1 = 0.0
dxi1 = 0.0
for i=1:max_iterations
dxi1 = -R(xi1)/dR(xi1)
xi1 += dxi1
if norm(dxi1) < tol
return xi1
end
end
info("x1 = $(ForwardDiff.get_value(x1_.data))")
info("n1 = $(ForwardDiff.get_value(n1_.data))")
info("x2 = $(ForwardDiff.get_value(x2))")
info("xi1 = $(ForwardDiff.get_value(xi1)), dxi1 = $(ForwardDiff.get_value(dxi1))")
info("-R(xi1) = $(ForwardDiff.get_value(-R(xi1)))")
info("dR(xi1) = $(ForwardDiff.get_value(dR(xi1)))")
error("find projection from master to slave: did not converge")
end
function project_from_slave_to_master{E<:MortarElements2D}(
master_element::Element{E}, x1::Vector, n1::Vector, x2_::DVTI, time::Float64;
tol=1.0e-10, max_iterations=20)
x2(xi2) = vec(get_basis(master_element, [xi2], time))*x2_
dx2(xi2) = vec(get_dbasis(master_element, [xi2], time))*x2_
cross2(a, b) = cross([a; 0], [b; 0])[3]
R(xi2) = cross2(x2(xi2)-x1, n1)
dR(xi2) = cross2(dx2(xi2), n1)
xi2 = 0.0
dxi2 = 0.0
for i=1:max_iterations
dxi2 = -R(xi2) / dR(xi2)
xi2 += dxi2
if norm(dxi2) < tol
return xi2
end
end
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
end
""" 2d mesh tie using ForwardDiff.
Construct .. + fc*la and C(d,la)=0
"""
function assemble!(problem::Problem{Mortar}, time::Float64, ::Type{Val{1}}, ::Type{Val{true}})
props = problem.properties
field_dim = get_unknown_field_dimension(problem)
field_name = get_parent_field_name(problem)
slave_elements = get_slave_elements(problem)
if field_name != "displacement"
error("mortar forwarddiff assembly: only displacement field with adjust=yes supported")
end
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)
gap = zeros(u)
C = zeros(la)
S = Set{Int64}()
# 1. update nodal normals for slave elements
tangents = zeros(u)
for element in slave_elements
conn = get_connectivity(element)
push!(S, conn...)
X1 = element("geometry", time)
u1 = Field([u[:,i] for i in conn])
x1 = X1 + u1
dN = get_dbasis(element, [0.0], time)
tangent = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
for nid in conn
tangents[:,nid] += tangent[:]
end
end
Q = [0.0 -1.0; 1.0 0.0]
normals = zeros(u)
for j in S
tangents[:,j] /= norm(tangents[:,j])
normals[:,j] = Q*tangents[:,j]
end
if props.rotate_normals
for j in S
normals[:,j] = -normals[:,j]
end
end
normals2 = Dict()
tangents2 = Dict()
for j in S
normals2[j] = normals[:,j]
tangents2[j] = tangents[:,j]
end
update!(slave_elements, "normal", time => normals2)
update!(slave_elements, "tangent", time => tangents2)
# 2. loop all slave elements
for slave_element in slave_elements
nsl = length(slave_element)
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])
# 3. loop all master elements
for master_element in slave_element("master elements", time)
nm = length(master_element)
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
# 3.1 calculate segmentation
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1], time)
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[2], time)
# xi1a = project_from_master_to_slave(slave_element, X2[1], time)
# xi1b = project_from_master_to_slave(slave_element, X2[2], time)
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
# 3.2. bi-orthogonal basis
De = zeros(nsl, nsl)
Me = zeros(nsl, nsl)
Ae = zeros(nsl, nsl)
if props.dual_basis
for ip in get_integration_points(slave_element, 3)
detJ = slave_element(ip, time, Val{:detJ})
w = ip.weight*detJ*l
xi = ip.coords[1]
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
N1 = vec(get_basis(slave_element, xi_s, time))
De += w*diagm(N1)
Me += w*N1*N1'
end
Ae = De*inv(Me)
else
Ae = eye(nsl)
end
# 3.3. loop integration points of one integration segment and calculate
# local mortar matrices
for ip in get_integration_points(slave_element, 3)
detJ = slave_element(ip, time, Val{:detJ})
w = ip.weight*detJ*l
#dN = get_dbasis(slave_element, ip, time)
#j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
#w = ip.weight*norm(j)*l
xi = ip.coords[1]
xi_s = dot([1/2*(1-xi); 1/2*(1+xi)], xi1)
N1 = vec(get_basis(slave_element, xi_s, time))
Phi = Ae*N1
# project gauss point from slave element to master element in direction n_s
x_s = N1*x1 # coordinate in gauss point
n_s = N1*n1 # normal direction in gauss point
#xi_m = project_from_slave_to_master(master_element, X_s, n_s, time)
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2, time)
N2 = vec(get_basis(master_element, xi_m, time))
x_m = N2*x2
la_s = Phi*la1
gn = dot(n_s, x_s-x_m)
u_s = N1*u1
u_m = N2*u2
X_s = N1*X1
X_m = N2*X2
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[:,slave_element_nodes] += w*(u_s-u_m)*Phi'
if props.adjust
G = w*(X_s-X_m)*Phi'
gap[:,slave_element_nodes] += G
end
end
end # master elements done
end # slave elements done, contact virtual work ready
C = gap
info("interface residual ready")
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)
#out = ForwardDiff.JacobianResult(x)
#ForwardDiff.jacobian!(out, calculate_interface)
A = ForwardDiff.jacobian(calculate_interface, x)
#b = -ForwardDiff.value(calculate_interface, x)
b = -calculate_interface(x)
# 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]
empty!(problem.assembly)
problem.assembly.K = K
problem.assembly.C1 = C1
problem.assembly.C2 = C2
problem.assembly.D = D
problem.assembly.f = f
problem.assembly.g = g
end
#=
""" Find segment from slave element corresponding to master element nodes.
Parameters
@@ -270,3 +528,6 @@ function assemble!(problem::Problem{Mortar}, time::Real, ::Type{Val{2}})
return problem.assembly
end
=#
src/solvers.jl
+16
View File
@@ -103,6 +103,9 @@ function get_field_assembly(solver::Solver; show_info=true)
M = sparse(M, solver.ndofs, solver.ndofs)
K = sparse(K, solver.ndofs, solver.ndofs)
if nnz(K) == 0
warn("Field assembly seems to be empty. Check that elements are pushed to problem and formulation is correct.")
end
Kg = sparse(Kg, solver.ndofs, solver.ndofs)
f = sparse(f, solver.ndofs, 1)
fg = sparse(fg, solver.ndofs, 1)
@@ -171,6 +174,14 @@ function get_boundary_assembly(solver::Solver)
return K, C1, C2, D, f, g
end
function resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
(n == A.n) && (m == A.m) && return
@assert n >= A.n
@assert m >= A.m
append!(A.colptr, A.colptr[end]*ones(Int, m-A.m))
A.n = n
A.m = m
end
"""
Given C and g, construct new basis such that v = P*u + g
@@ -450,6 +461,11 @@ function NonlinearSolver(problems...)
end
return solver
end
function NonlinearSolver(name::ASCIIString, problems::Problem...)
solver = NonlinearSolver(problems...)
solver.name = name
return solver
end
### Linear quasistatic solver
src/solvers_modal.jl
+7 -2
View File
@@ -23,7 +23,12 @@ function Modal(nev=10, which=:SM)
solver = Modal(false, Vector(), Matrix(), nev, which)
end
function call(solver::Solver{Modal}; debug=false)
function call(solver::Solver{Modal}; show_info=true, debug=false)
show_info && info(repeat("-", 80))
show_info && info("Starting natural frequency solver")
show_info && info("Increment time t=$(round(solver.time, 3))")
show_info && info(repeat("-", 80))
initialize!(solver)
# assemble all field problems
info("Assembling problems ...")
tic()
@@ -36,7 +41,7 @@ function call(solver::Solver{Modal}; debug=false)
end
t1 = round(toq(), 2)
info("Assembled in $t1 seconds.")
M, K, Kg, f = get_field_assembly(solver; with_mass_matrix=true)
M, K, Kg, f = get_field_assembly(solver)
Kb, C1, C2, D, fb, g = get_boundary_assembly(solver)
K = K + Kb
f = f + fb
src/sparse.jl
+20
View File
@@ -172,3 +172,23 @@ function size(A::SparseMatrixCOO, idx::Int)
return size(A)[idx]
end
""" Matrix norm. Automatically convert to dense when asking for 2-norm for small matrices. """
function Base.norm(A::SparseMatrixCOO, p=Inf; maxdim=1000)
dim = size(A, 1)
if p == 2 && dim > maxdim
info("Assembly norm: dim = $dim > $maxdim and p=$p, not making dense matrices for operation.")
return 0.0
end
if p == 2
return norm(full(A), p)
else
return norm(sparse(A), p)
end
end
function isapprox(A::SparseMatrixCOO, B::SparseMatrixCOO)
A2 = sparse(A)
B2 = sparse(B, size(A2)...)
return isapprox(A2, B2)
end