mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-09-17 17:22:10 +00:00
first attemps to make properly linearized version of mortar projection for finite sliding. not working at the moment, it has convergence issues.
This commit is contained in:
Jukka Aho
@@ -132,6 +132,35 @@ function call{E}(element::Element{E}, xi::VecOrIP, time::Float64=0.0)
|
||||
return get_basis(element, xi)
|
||||
end
|
||||
|
||||
""" Given a list of elementa and nodes, find a subset of elements
|
||||
containing nodes.
|
||||
"""
|
||||
function find_elements(elements, nodes)
|
||||
s = Set{Element}()
|
||||
for element in elements
|
||||
conn = get_connectivity(element)
|
||||
for j in nodes
|
||||
if j in conn
|
||||
push!(s, element)
|
||||
break
|
||||
end
|
||||
end
|
||||
end
|
||||
return collect(s)
|
||||
end
|
||||
|
||||
function get_gdofs(element::Element)
|
||||
return get_gdofs(element, 1)
|
||||
end
|
||||
|
||||
function get_dbasis{E}(element::Element{E}, ip::IntegrationPoint)
|
||||
return get_dbasis(E, ip.xi)
|
||||
end
|
||||
|
||||
function get_basis{E, T<:Real}(element::Element{E}, xi::T)
|
||||
return get_basis(E, xi)
|
||||
end
|
||||
|
||||
""" Return dual basis transformation matrix Ae. """
|
||||
function get_dualbasis(element::Element, time::Real)
|
||||
if length(element.A) == 0
|
||||
|
||||
@@ -216,6 +216,15 @@ function Base.(:*){T<:Real}(c::T, field::DVTI)
|
||||
return DVTI(c*field.data)
|
||||
end
|
||||
|
||||
""" Multiply DVTI field with another vector T. Vector length
|
||||
must match to the field length and this can be used mainly
|
||||
for interpolation purposes, i.e., u = ∑ Nᵢuᵢ
|
||||
"""
|
||||
function Base.(:*)(T::Vector, f::DVTI)
|
||||
@assert length(T) == length(f)
|
||||
return sum([T[i]*f[i] for i=1:length(f)])
|
||||
end
|
||||
|
||||
function Base.vec(field::DVTI)
|
||||
return [field.data...;]
|
||||
end
|
||||
|
||||
@@ -140,3 +140,11 @@ end
|
||||
(xi) -> [ 1.0, xi[1], xi[2], xi[3], xi[1]^2,
|
||||
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
|
||||
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
|
||||
|
||||
+5
-715
@@ -19,7 +19,7 @@ b) Remove inactive inequality constraints in assembly level. This is done in
|
||||
|
||||
"""
|
||||
type Mortar <: BoundaryProblem
|
||||
formulation :: Symbol # :total or :incremental
|
||||
formulation :: Symbol # :total, :incremental, :autodiff
|
||||
dual_basis :: Bool
|
||||
inequality_constraints :: Bool # Launch PDASS to solve inequality constraints
|
||||
normal_condition :: Symbol # Tie or Contact
|
||||
@@ -31,10 +31,11 @@ type Mortar <: BoundaryProblem
|
||||
always_in_slip :: Vector{Int64} # nodes in this list always in slip
|
||||
contact :: Bool
|
||||
friction :: Bool
|
||||
gap_sign :: Int # gap sign convention
|
||||
end
|
||||
|
||||
function Mortar()
|
||||
Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false)
|
||||
Mortar(:total, true, false, :Tie, :Stick, Inf, false, [], [], [], false, false, -1)
|
||||
end
|
||||
|
||||
function get_unknown_field_name(::Type{Mortar})
|
||||
@@ -51,719 +52,8 @@ macro debug(msg)
|
||||
end
|
||||
|
||||
include("mortar_2d.jl")
|
||||
|
||||
|
||||
### Mortar projection calculation for 3d cases
|
||||
|
||||
"""
|
||||
Construct auxiliary plane for surface.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
x::Array{Float64, 2}
|
||||
Node coordinates
|
||||
ximp::Array{Float64, 1}
|
||||
Element mid-point in dimensionless mother element coordinates ξ
|
||||
normals::Array{Float64, 2}
|
||||
Normal directions in nodes
|
||||
|
||||
Returns
|
||||
-------
|
||||
x0, Q
|
||||
x0::Array{Float64, 1} - origo of auxiliary plane
|
||||
Q::Array{Float64, 2} - orthogonal basis, first vector is normal direction
|
||||
and two rest vectors create orthonormal right-handed basis.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Calculate auxiliary plane given nodal coordinates, midpoint of mother element,
|
||||
node normals and suitable function space:
|
||||
|
||||
julia> xquad = [
|
||||
... -2.5 2.5 2.0 -2.0
|
||||
... -2.0 -2.0 2.3 2.0
|
||||
... 1.0 0.7 0.0 1.0]
|
||||
julia> m_midpoint = [0.0, 0.0]
|
||||
julia> normals = [
|
||||
... 0.05989060 0.0590504 0.225612 0.2445800
|
||||
... -0.00748633 0.1670810 0.182034 -0.0305725
|
||||
... 0.99817700 0.9841730 0.957059 0.9691470]
|
||||
julia> basis(xi) = [
|
||||
... (1-xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1+xi[2])/4
|
||||
... (1-xi[1])(1+xi[2])/4]'
|
||||
julia> x0, Q = create_auxiliary_plane(xquad, mmidpoint, normals, basis)
|
||||
julia> x0
|
||||
3-element Array{Float64,1}:
|
||||
0.0
|
||||
0.075
|
||||
0.675
|
||||
julia> Q
|
||||
3x3 Array{Float64,2}:
|
||||
0.148586 0.988899 0.0
|
||||
0.0784519 -0.0117877 0.996848
|
||||
0.985783 -0.148118 -0.0793325
|
||||
|
||||
Notes
|
||||
-----
|
||||
- Midpoint in mother element typically (0, 0) for quadrangles and (1/3, 1/3)
|
||||
for triangles.
|
||||
- Uses Gram-Schmidt process to find orthogonal basis
|
||||
- [1](http://www.math.umn.edu/~olver/aims_/qr.pdf)
|
||||
- [2](http://www.ecs.umass.edu/ece/ece313/Online_help/gram.pdf)
|
||||
- [3](http://www.terathon.com/code/tangent.html)
|
||||
|
||||
"""
|
||||
# function create_auxiliary_plane(x, ximp, normals, basis)
|
||||
function create_auxiliary_plane{E}(element::Element{E}, time::Real)
|
||||
# proj(u, v) = dot(v, u) / dot(u, u) * u
|
||||
# xi = [1.0/3.0, 1.0/3.0]
|
||||
|
||||
xi = get_reference_element_midpoint(E)
|
||||
x0 = element("geometry", xi, time)
|
||||
ntbasis = element("normal-tangential coordinates", xi, time)
|
||||
return x0, ntbasis
|
||||
#=
|
||||
n = element("normal-tangential coordinates", xi, time)[:, 1]
|
||||
n /= norm(n)
|
||||
# gram-schmidt
|
||||
u1 = n
|
||||
j = indmax(abs(u1))
|
||||
v2 = zeros(3)
|
||||
v2[mod(j,3)+1] = 1.0
|
||||
u2 = v2 - proj(u1, v2)
|
||||
u3 = cross(u1, u2)
|
||||
t1 = u2/norm(u2)
|
||||
t2 = u3/norm(u3)
|
||||
new_basis = [n t1 t2]
|
||||
return x0, new_basis
|
||||
=#
|
||||
end
|
||||
|
||||
"""
|
||||
Project point q onto a plane given by a point p and normal n.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
q::Array{Float64, 2}
|
||||
point to project (row vector)
|
||||
x0::Array{Float64, 2}
|
||||
origo of plane
|
||||
n::Array{Float64, 2}
|
||||
normal vector of plane
|
||||
|
||||
Returns
|
||||
-------
|
||||
y::Array{Float64, 2}
|
||||
projected point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> p = [-0.5 -1.0 4.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> n = [0.1485860 0.0784519 0.9857830]'
|
||||
julia> project_node_to_auxiliary_plane(p, x0, n)
|
||||
3-element Array{Float64,1}:
|
||||
0.963455
|
||||
-1.2447
|
||||
0.925247
|
||||
|
||||
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)
|
||||
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.")
|
||||
end
|
||||
return qproj[2:3]
|
||||
end
|
||||
|
||||
"""
|
||||
Find edge intersections of two planar arbitrary shape polygons.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64,2}
|
||||
M::Array{Float64,2}
|
||||
|
||||
Matrices with size (2, n) where n is number of vertices of each polygon.
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64,2}
|
||||
Intersection points of polygons
|
||||
n::Array{Float64,2}
|
||||
Neighbour info matrix with size (ns, mn). This keeps information which
|
||||
edges of polygons are intersecting. See further explanation in example
|
||||
below.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Find intersection points of two triangles:
|
||||
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 1 3/2]'
|
||||
julia> P, n = get_edge_intersections(S, M)
|
||||
julia> P
|
||||
2x4 Array{Float64,2}:
|
||||
1.0 1.75 0.0 0.0
|
||||
0.0 0.0 0.5 1.25
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 1 0
|
||||
0 0 0
|
||||
1 0 1)
|
||||
|
||||
So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25).
|
||||
"Neighbour matrix" can be interpreted as following:
|
||||
|
||||
1 1 0 <--> First edge of S intersects edges 1 and 2 of M
|
||||
0 0 0 <--> Second edge of S doesn't intersect at all
|
||||
1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M
|
||||
|
||||
"""
|
||||
function get_edge_intersections(S::Matrix, M::Matrix)
|
||||
ns = size(S, 2)
|
||||
nm = size(M, 2)
|
||||
P = zeros(2, 0)
|
||||
n = zeros(Int64, ns, nm)
|
||||
k = 0
|
||||
for i=1:ns
|
||||
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
|
||||
r = A\b
|
||||
if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found
|
||||
k += 1
|
||||
f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i])
|
||||
f = f''
|
||||
P = hcat(P, f)
|
||||
n[i, j] = 1
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return P, n
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Find any points laying inside or border of triangle.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
Y::Array{Float64, 2}
|
||||
Triangle coordinates in 2×3 matrix
|
||||
X::Array{Float64, 2}
|
||||
List of points to test in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64, 2}
|
||||
List of points in triangle in 2×m matrix, where m is number of points inside triangle
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points
|
||||
julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests
|
||||
julia> points_in_triangle(S, pts)
|
||||
2x2 Array{Float64,2}:
|
||||
1.0 0.5
|
||||
1.5 1.5
|
||||
"""
|
||||
function get_points_inside_triangle(Y::Matrix, X::Matrix)
|
||||
@assert size(Y, 2) == 3 # "Point in TRIANGLE..."
|
||||
P = zeros(2, 0)
|
||||
v0 = Y[:,2] - Y[:,1]
|
||||
v1 = Y[:,3] - Y[:,1] # find interior points of X in Y
|
||||
d00 = (v0'*v0)[1]
|
||||
d01 = (v0'*v1)[1]
|
||||
d11 = (v1'*v1)[1] # using baricentric coordinates
|
||||
id = 1/(d00*d11 - d01*d01)
|
||||
for i=1:size(X, 2)
|
||||
v2 = X[:,i] - Y[:,1]
|
||||
d02 = (v0'*v2)[1]
|
||||
d12 = (v1'*v2)[1]
|
||||
u = (d11*d02-d01*d12)*id
|
||||
v = (d00*d12-d01*d02)*id
|
||||
if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary
|
||||
P = hcat(P, X[:,i]'')
|
||||
end
|
||||
end
|
||||
return P
|
||||
end
|
||||
|
||||
"""
|
||||
Determine is point P inside or on boudary of polygon X.
|
||||
|
||||
http://paulbourke.net/geometry/polygonmesh/#insidepoly
|
||||
"""
|
||||
function is_point_inside_convex_polygon(P, X)
|
||||
x, y = P
|
||||
for i=1:length(X)
|
||||
x0, y0 = X[i]
|
||||
x1, y1 = X[mod(i, length(X))+1]
|
||||
if (y-y0)*(x1-x0) - (x-x0)*(y1-y0) < 0
|
||||
return false
|
||||
end
|
||||
end
|
||||
return true
|
||||
end
|
||||
|
||||
function get_points_inside_convex_polygon(pts, X)
|
||||
# TODO: Make more readable
|
||||
X2 = [X[:,i] for i=1:size(X,2)]
|
||||
c = filter(P->is_point_inside_convex_polygon(P, X2), [pts[:,i] for i=1:size(pts, 2)])
|
||||
return length(c) == 0 ? zeros(2, 0) : hcat(c...)
|
||||
end
|
||||
|
||||
""" Return unique objects with some given tolerance. This is used in next function
|
||||
because traditional unique() command returns row vectors as non-unique if they
|
||||
differs only a "little".
|
||||
"""
|
||||
function uniquetol(P, dim::Int; args...)
|
||||
@assert dim == 2
|
||||
items = Vector{Float64}[P[:,i] for i=1:size(P,dim)]
|
||||
new_items = Vector{Float64}[]
|
||||
for item in items
|
||||
has_found = false
|
||||
for new_item in new_items
|
||||
if isapprox(item, new_item; args...)
|
||||
has_found = true
|
||||
break
|
||||
end
|
||||
end
|
||||
if !has_found
|
||||
push!(new_items, item)
|
||||
end
|
||||
end
|
||||
return reshape([new_items...;], length(new_items[]), length(new_items))
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Make polygon clipping of shapes S and M.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64, 2}
|
||||
M::Array{Float64, 2}
|
||||
Shapes to clip. Needs to be triangles at the moment.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64, 2}, Array{Float64, 2}
|
||||
- Polygon vertices in 2×n matrix, sorted in counter-clockwise order.
|
||||
- 3×3 "neighbouring" matrix, see example.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 2 2]'
|
||||
julia> P, n = clip_polygon(S, M)
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]')
|
||||
1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S
|
||||
0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S
|
||||
|
||||
"""
|
||||
function clip_polygon(S::Matrix, M::Matrix)
|
||||
P1, neighbours = get_edge_intersections(M, S)
|
||||
#P2 = get_points_inside_triangle(M, S)
|
||||
#P3 = get_points_inside_triangle(S, M)
|
||||
P2 = get_points_inside_convex_polygon(M, S)
|
||||
P3 = get_points_inside_convex_polygon(S, M)
|
||||
# info("polygon clipping: P1 = $P1")
|
||||
# info("polygon clipping: P2 = $P2")
|
||||
# info("polygon clipping: P3 = $P3")
|
||||
# info("hcat P = $P")
|
||||
P = hcat(P1, P2, P3)
|
||||
if length(P) == 0
|
||||
return nothing, nothing
|
||||
end
|
||||
P = uniquetol(P, 2)
|
||||
meanval = mean(P, 2)
|
||||
tmp = P .- meanval
|
||||
angles = atan2(tmp[2,:], tmp[1,:])
|
||||
angles = reshape(angles, length(angles))
|
||||
order = sortperm(angles)
|
||||
return P[:, order], neighbours
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Calculate polygon geometric center point
|
||||
|
||||
Parameters
|
||||
----------
|
||||
P::Array{Float64, 2}
|
||||
Polygon vertices in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float63, 2}
|
||||
Center point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> C = get_polygon_cp(P)
|
||||
2x1 Array{Float64,2}:
|
||||
1.039740
|
||||
0.804701
|
||||
|
||||
"""
|
||||
function calculate_polygon_centerpoint(P::Matrix)
|
||||
n = size(P, 2)
|
||||
A = 0.0
|
||||
for i=1:n
|
||||
A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i])
|
||||
end
|
||||
Cx = 0.0
|
||||
Cy = 0.0
|
||||
for i=1:n
|
||||
inext = mod(i, n)+1
|
||||
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]
|
||||
end
|
||||
|
||||
"""
|
||||
Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p::Array{Float64,1}
|
||||
point in auxiliary plane, in (n,t1,t2) coordinate system
|
||||
x0::Array{Float64,1}
|
||||
origo of auxiliary plane cs
|
||||
Q::Array{Float64,2}
|
||||
basis of auxiliary plane cs
|
||||
x::Array{Float64,2}
|
||||
surface node coords
|
||||
basis::Array{Float64,2}
|
||||
surface basis functions
|
||||
dbasis::Array{Float64,2}
|
||||
partial derivatives of surface basis functions
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64,2}
|
||||
solution vector (d, ξ₁, ξ₂) where d is distance to surface
|
||||
|
||||
Examples
|
||||
--------
|
||||
Define surface with node points, basis + dbasis
|
||||
|
||||
julia> xquad = [
|
||||
... -2.5 -2.0 1.0
|
||||
... 2.5 -2.0 0.7
|
||||
... 2.0 2.3 0.0
|
||||
... -2.0 2.0 1.0]'
|
||||
julia> basis(xi) = [
|
||||
... (1-xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1+xi[2])/4
|
||||
... (1-xi[1])(1+xi[2])/4]
|
||||
julia> dbasis(xi) = [
|
||||
... -(1-xi[2])/4 -(1-xi[1])/4
|
||||
... (1-xi[2])/4 -(1+xi[1])/4
|
||||
... (1+xi[2])/4 (1+xi[1])/4
|
||||
... -(1+xi[2])/4 (1-xi[1])/4]
|
||||
|
||||
We aim to find point p, which we first project to auxiliary plane defined as following
|
||||
julia> p = [-2.5 -2.0 1.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> Q = [
|
||||
... 0.1485860 0.9888990 0.0000000
|
||||
... 0.0784519 -0.0117877 0.9968480
|
||||
... 0.9857830 -0.1481180 -0.0793325]
|
||||
|
||||
Our projected point is therefore
|
||||
julia> n = Q[:,1] # first component is normal direction
|
||||
julia> ph = project_node_to_auxiliary_plane(p, x0, n)
|
||||
julia> ph = Q'(ph-x0)
|
||||
julia> ph
|
||||
3x1 Array{Float64,2}:
|
||||
1.33264e-7
|
||||
-2.49593
|
||||
-2.09424
|
||||
|
||||
Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we
|
||||
project it back to surface defined by xquad*basis
|
||||
|
||||
julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis)
|
||||
julia> theta
|
||||
3x1 Array{Float64,2}:
|
||||
-0.213874
|
||||
-0.999999
|
||||
-1.0
|
||||
|
||||
We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad
|
||||
[-2.5 -2.0 1.0]' correctly.
|
||||
|
||||
julia> xquad*basis(theta[2:3])
|
||||
3-element Array{Float64,1}:
|
||||
-2.5
|
||||
-2.0
|
||||
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)
|
||||
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]
|
||||
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
|
||||
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
|
||||
|
||||
typealias MortarElements3D Union{Tri3, Quad4}
|
||||
|
||||
function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar},
|
||||
slave_element::Element{E}, time::Real)
|
||||
haskey(slave_element, "master elements") || return
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
|
||||
props = problem.properties
|
||||
if props.formulation == :Standard && props.normal_condition == :Contact
|
||||
error("for contact choose Dual formulation.""")
|
||||
end
|
||||
|
||||
# create auxiliary plane and project slave nodes to it
|
||||
# x0 = origo, Q = local basis
|
||||
x0, Q = create_auxiliary_plane(slave_element, time)
|
||||
|
||||
# 1. project slave nodes to auxiliary plane
|
||||
Sl = Vector{Float64}[]
|
||||
for p in slave_element("geometry", time)
|
||||
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
S = hcat(Sl...)
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
|
||||
# if distance between elements is "far enough" cannot expect contact
|
||||
if (props.normal_condition == :Contact) || props.inequality_constraints
|
||||
slave_midpoint = slave_element("geometry", [0.0, 0.0], time)
|
||||
master_midpoint = master_element("geometry", [0.0, 0.0], time)
|
||||
if norm(slave_midpoint - master_midpoint) > props.minimum_distance
|
||||
continue
|
||||
end
|
||||
end
|
||||
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
# 2. project master nodes to auxiliary plane
|
||||
M = Vector{Float64}[]
|
||||
for p in master_element("geometry", time)
|
||||
push!(M, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
M = hcat(M...)
|
||||
|
||||
# 3. create polygon clipping on auxiliary plane
|
||||
P = nothing
|
||||
neighbours = nothing
|
||||
try
|
||||
P, neighbours = clip_polygon(S, M)
|
||||
catch
|
||||
info("polygon clipping failed")
|
||||
info("S = ")
|
||||
dump(S)
|
||||
info("M = ")
|
||||
dump(M)
|
||||
info("original Sl = ")
|
||||
info(Sl)
|
||||
error("cannot continue")
|
||||
end
|
||||
isa(P, Void) && continue # no clipping
|
||||
|
||||
# shared edge but no shared volume. skipping
|
||||
size(P, 2) < 3 && continue
|
||||
|
||||
C = calculate_polygon_centerpoint(P)
|
||||
npts = size(P, 2) # number of vertices in polygon
|
||||
|
||||
# loop vertices and create temporary integrate cells
|
||||
# TODO: basically when npts == 3 or npts == 4 we could integrate without splitting to cells.
|
||||
nnodes = size(slave_element, 2)
|
||||
C1S3 = zeros(3*nnodes, 3*nnodes)
|
||||
C1M3 = zeros(3*nnodes, 3*nnodes)
|
||||
|
||||
for pnt=1:npts # integration of mortar matrices begin
|
||||
cell = Field(Vector{Float64}[C, P[:,pnt], P[:,mod(pnt,npts)+1]])
|
||||
|
||||
# calculate slave side projection matrix D
|
||||
# construct dual basis
|
||||
Ae = zeros(nnodes, nnodes)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
if problem.properties.formulation == :Dual # Construct dual basis
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = get_basis(Tri3, ip.xi)
|
||||
xi = vec(N*cell)
|
||||
theta = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
|
||||
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], cell[j]') for j=1:length(cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
De += wC*diagm(vec(N1))
|
||||
Me += wC*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
end
|
||||
for i=1:field_dim
|
||||
C1S3[i:field_dim:end,i:field_dim:end] += De
|
||||
end
|
||||
|
||||
# Calculate master side projection matrix M
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
# gauss point in auxiliary plane
|
||||
#N = get_basis(E, ip.xi)
|
||||
N = get_basis(Tri3, ip.xi)
|
||||
xi = vec(N*cell) # xi defined in auxilary plane
|
||||
|
||||
# find projection of gauss point to master and slave elements
|
||||
theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
|
||||
theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time)
|
||||
xi_slave = theta1[2:3]
|
||||
xi_master = theta2[2:3]
|
||||
|
||||
# evaluate shape functions values in gauss point and add contribution to matrices
|
||||
N1 = slave_element(xi_slave, time)
|
||||
N2 = master_element(xi_master, time)
|
||||
|
||||
# jacobian determinant on integration cell
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
|
||||
# extend matrices according to the problem dimension (3)
|
||||
@assert length(slave_dofs) == length(master_dofs)
|
||||
Me = wC*Ae*N1'*N2
|
||||
for k=1:field_dim
|
||||
C1M3[k:field_dim:end,k:field_dim:end] += Me
|
||||
end
|
||||
end
|
||||
end # integration of mortar matrices done.
|
||||
|
||||
# constraints in normal-tangential direction and initial weighted gap
|
||||
X1 = vec(slave_element("geometry", time))
|
||||
X2 = vec(master_element("geometry", time))
|
||||
Q_ = slave_element("normal-tangential coordinates", time)
|
||||
Z = zeros(3, 3)
|
||||
if nnodes == 3
|
||||
Q3 = [Q Z Z; Z Q Z; Z Z Q]
|
||||
elseif nnodes == 4
|
||||
Q3 = [Q Z Z Z; Z Q Z Z; Z Z Q Z; Z Z Z Q]
|
||||
end
|
||||
D3 = zeros(3*nnodes, 3*nnodes)
|
||||
C2S3 = Q3'*C1S3
|
||||
C2M3 = Q3'*C1M3
|
||||
G = -(C2S3*X1 - C2M3*X2)
|
||||
|
||||
# complementarity condition
|
||||
if haskey(slave_element, "displacement")
|
||||
u1 = vec(slave_element("displacement", time))
|
||||
else
|
||||
u1 = zeros(3*nnodes)
|
||||
end
|
||||
if haskey(master_element, "displacement")
|
||||
u2 = vec(master_element("displacement", time))
|
||||
else
|
||||
u2 = zeros(3*nnodes)
|
||||
end
|
||||
x1 = X1 + u1
|
||||
x2 = X2 + u2
|
||||
if haskey(slave_element, "reaction force")
|
||||
la = vec(slave_element("reaction force", time))
|
||||
else
|
||||
la = zeros(3*nnodes)
|
||||
end
|
||||
g = -(C2S3*x1 - C2M3*x2)
|
||||
c = Q3'*la - g
|
||||
inactive_nodes = find(c[1:field_dim:end] .<= 0)
|
||||
active_nodes = find(c[1:field_dim:end] .> 0)
|
||||
|
||||
# normal constraint: remove inactive nodes if normal condition is set to contact
|
||||
if problem.properties.normal_condition == :Contact
|
||||
for j in inactive_nodes
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
G[dofs] = 0
|
||||
C1S3[dofs,:] = 0
|
||||
C1M3[dofs,:] = 0
|
||||
C2S3[dofs,:] = 0
|
||||
C2M3[dofs,:] = 0
|
||||
end
|
||||
end
|
||||
|
||||
# tangential constraint: stick or slip
|
||||
if problem.properties.tangential_condition == :Slip
|
||||
D3 = copy(C2S3)
|
||||
D3[1:field_dim:end, :] = 0
|
||||
C2S3[2:field_dim:end, :] = 0
|
||||
C2M3[2:field_dim:end, :] = 0
|
||||
C2S3[3:field_dim:end, :] = 0
|
||||
C2M3[3:field_dim:end, :] = 0
|
||||
end
|
||||
|
||||
# add contributions
|
||||
add!(assembly.C1, slave_dofs, slave_dofs, C1S3)
|
||||
add!(assembly.C1, slave_dofs, master_dofs, -C1M3)
|
||||
add!(assembly.C2, slave_dofs, slave_dofs, C2S3)
|
||||
add!(assembly.C2, slave_dofs, master_dofs, -C2M3)
|
||||
add!(assembly.D, slave_dofs, slave_dofs, D3)
|
||||
add!(assembly.c, slave_dofs, c)
|
||||
add!(assembly.g, slave_dofs, G)
|
||||
end
|
||||
end
|
||||
|
||||
include("mortar_2d_autodiff.jl")
|
||||
include("mortar_3d.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)
|
||||
|
||||
+2
-3
@@ -217,7 +217,6 @@ function project_from_master_to_slave{S,M}(slave::Element{S}, master::Element{M}
|
||||
error("find projection from master to slave: did not converge")
|
||||
end
|
||||
|
||||
|
||||
# Mortar assembly 2d
|
||||
|
||||
# quadratic not tested yet
|
||||
@@ -546,8 +545,8 @@ function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mor
|
||||
# Calculate normal-tangential constraints and weighted gap
|
||||
C2S2 = Q2'*C1S2
|
||||
C2M2 = Q2'*C1M2
|
||||
G += -(C2S2*X1 - C2M2*X2)
|
||||
g += -(C2S2*x1 - C2M2*x2)
|
||||
G += props.gap_sign*(C2S2*X1 - C2M2*X2)
|
||||
g += props.gap_sign*(C2S2*x1 - C2M2*x2)
|
||||
|
||||
# Add contributions
|
||||
add!(local_assembly.C1, slave_dofs, slave_dofs, C1S2)
|
||||
|
||||
@@ -0,0 +1,674 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
""" Find segment from slave element corresponding to master element nodes.
|
||||
x1_, n1_
|
||||
slave element geometry and normal direction
|
||||
|
||||
x2_ master element nodes to project onto slave
|
||||
"""
|
||||
function project_from_master_to_slave{E<:MortarElements2D}(
|
||||
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector;
|
||||
tol=1.0e-10, max_iterations=20)
|
||||
|
||||
function x1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*x1_
|
||||
end
|
||||
|
||||
function dx1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*x1_
|
||||
end
|
||||
|
||||
function n1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*n1_
|
||||
end
|
||||
|
||||
function dn1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*n1_
|
||||
end
|
||||
|
||||
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;
|
||||
tol=1.0e-10, max_iterations=20)
|
||||
|
||||
function x2(xi2)
|
||||
N = get_basis(E, xi2)
|
||||
return vec(N)*x2_
|
||||
end
|
||||
|
||||
function dx2(xi2)
|
||||
dN = get_dbasis(E, xi2)
|
||||
return vec(dN)*x2_
|
||||
end
|
||||
|
||||
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
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly, problem::Problem{Mortar}, slave_element::Element{E}, time::Real,
|
||||
::Type{Val{:forwarddiff_old}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
nnodes = size(slave_element, 2)
|
||||
X1 = slave_element("geometry", time)
|
||||
#u1 = slave_element("displacement", time)
|
||||
#x1 = X1 + u1
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
|
||||
X = spzeros(10000, 1)
|
||||
for element in get_elements(problem)
|
||||
conn = get_connectivity(element)
|
||||
geom = element("geometry", time)
|
||||
for (c, g) in zip(conn, geom)
|
||||
dofs = [field_dim*(c-1)+1, field_dim*(c-1)+2]
|
||||
X[dofs] = g
|
||||
end
|
||||
end
|
||||
|
||||
# here x does not mean deformed configuration
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
if length(x) == 0
|
||||
info("mortar_2d_autodiff: length(x) == 0")
|
||||
# resize solution vectors according to initial configuration of this problem
|
||||
X = vec(full(sparse(findnz(X)...)))
|
||||
x = zeros(length(X)*2)
|
||||
else
|
||||
# resize initial configuration to match real dimension
|
||||
I, J, V = findnz(X)
|
||||
X = vec(full(sparse(I, J, V, length(problem.assembly.u), 1)))
|
||||
end
|
||||
ndofs = round(Int, length(x)/2)
|
||||
# at the end we should have
|
||||
# info("mortar_2d_autodiff: size of x = $(size(x))")
|
||||
# info("mortar_2d_autodiff: size of X = $(size(X))")
|
||||
# info("mortar_2d_autodiff: ndofs = $ndofs")
|
||||
|
||||
""" Calculate normal vector for slave element nodes in current configuration. """
|
||||
function calculate_normals(u::Matrix)
|
||||
normals = zeros(u)
|
||||
# 1. update nodal normals
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
slave_normals = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
return slave_normals
|
||||
end
|
||||
|
||||
function calculate_mortar_projection(u::Matrix, n1::DVTI)
|
||||
B = SparseMatrixCOO{Real}([], [], [])
|
||||
|
||||
u1 = Field([u[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
u2 = Field([u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
#info("master element coordinate 1 = $(ForwardDiff.get_value(x2[1]))")
|
||||
#info("master element coordinate 2 = $(ForwardDiff.get_value(x2[2]))")
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution
|
||||
|
||||
# integrate slave side
|
||||
D = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
# jacobian of slave element in deformed state
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N = get_basis(slave_element, xi_s)
|
||||
D += w*diagm(vec(N))
|
||||
Me += w*N'*N
|
||||
end
|
||||
Ae = D*inv(Me)
|
||||
|
||||
# integrate master side
|
||||
M = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
# jacobian of slave element in deformed state
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_g = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_g)
|
||||
x_g = vec(N1)*x1
|
||||
n_g = vec(N1)*n1
|
||||
#info("slave gauss point coordinates $(ForwardDiff.get_value(x_g))")
|
||||
#info("slave gauss point normal direction $(ForwardDiff.get_value(n_g))")
|
||||
xi_m = project_from_slave_to_master(master_element, x_g, n_g, x2)
|
||||
N2 = get_basis(master_element, xi_m)
|
||||
M += w*kron(Ae*N1', N2)
|
||||
end
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
for i=1:field_dim
|
||||
add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D)
|
||||
add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M)
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
return B
|
||||
end
|
||||
|
||||
function calculate_contact_rhs(x::Vector)
|
||||
ndofs = round(Int, length(x)/2)
|
||||
u = x[1:ndofs]
|
||||
la = x[ndofs+1:end]
|
||||
# info("calculate_contact_rhs: size of u = $(size(u))")
|
||||
# info("calculate_contact_rhs: size of X = $(size(X))")
|
||||
# info("calculate_contact_rhs: size of la = $(size(la))")
|
||||
# info("calculate_contact_rhs: ndofs = $ndofs")
|
||||
u2 = reshape(u, field_dim, round(Int, length(u)/field_dim))
|
||||
|
||||
normals = calculate_normals(u2)
|
||||
B = calculate_mortar_projection(u2, normals)
|
||||
B = sparse(B, ndofs, ndofs)
|
||||
fc = B' * la # contact force residual for r = fint + fc - fext = 0
|
||||
|
||||
N = SparseMatrixCOO{Real}([], [], [])
|
||||
T = SparseMatrixCOO{Real}([], [], [])
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
add!(N, dofs, [dofs[1]], reshape(normals[i], 2, 1))
|
||||
add!(T, dofs, [dofs[2]], reshape(Q'*normals[i], 2, 1))
|
||||
end
|
||||
N = sparse(N, ndofs, ndofs)
|
||||
T = sparse(T, ndofs, ndofs)
|
||||
gn = -N*B*(X+u)
|
||||
gt = -T*B*(X+u)
|
||||
# gn = -N*B*u
|
||||
lan = N*la
|
||||
lat = T*la
|
||||
cn = 1.0e3
|
||||
C = lan - max(0, lan - cn*gn) + lat
|
||||
# C = lan - gn + lat - gt <-- ihan viturallensa
|
||||
# C = gn+gt <-- not working
|
||||
# C = B*(X+u)
|
||||
# C = B*u <- pitää kiinni, "tie".
|
||||
# C = N*B*u + T*la <- palikat menee väärään suuntaan
|
||||
# C = -N*B*(X+u) + T*la <- toimii suht hyvin mut kääntyy väärään suuntaan (t-suunnassa)
|
||||
# C = -N*B*(X+u) - T*la <- sama
|
||||
# C = -N*B*(X+u) - T*B*la <- sama
|
||||
# C = N*la + T*la - max(0, N*la + N*B*(X+u))
|
||||
cond = lan[1:field_dim:end] - cn*gn[1:field_dim:end]
|
||||
all_nodes = slave_element_nodes
|
||||
inactive_nodes = find(cond .<= 0)
|
||||
active_nodes = find(cond .> 0)
|
||||
inactive_nodes = setdiff(all_nodes, inactive_nodes)
|
||||
active_nodes = setdiff(all_nodes, active_nodes)
|
||||
info("S = $all_nodes, I = $inactive_nodes, A = $active_nodes")
|
||||
info("lambda = $(ForwardDiff.get_value(lan[slave_dofs]))")
|
||||
info("gn = $(ForwardDiff.get_value(gn[slave_dofs]))")
|
||||
#for j in active_nodes
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# C[dofs] = 0
|
||||
#end
|
||||
|
||||
# C = -(N+T)*B*(X+u)
|
||||
# C = -B*(X+u)
|
||||
|
||||
return [fc; C]
|
||||
|
||||
end
|
||||
|
||||
A, allresults = ForwardDiff.jacobian(calculate_contact_rhs, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = A[1:ndofs,ndofs+1:end]'
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = b[1:ndofs]
|
||||
g = b[ndofs+1:end]
|
||||
C2[2:field_dim:end] = 0
|
||||
g[2:field_dim:end] = 0
|
||||
|
||||
# C2[2:field_dim:end] = 0
|
||||
# g[2:field_dim:end] = 0
|
||||
# inactives = find(g[1:field_dim:end] .<= 0)
|
||||
# actives = find(g[1:field_dim:end] .> 0)
|
||||
# inactives = setdiff(slave_element_nodes, inactives)
|
||||
# actives = setdiff(slave_element_nodes, actives)
|
||||
# info("all nodes = $slave_element_nodes, inactives = $inactives, actives = $actives")
|
||||
# info("g = $g")
|
||||
|
||||
# for j in inactives
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# K[dofs,:] = 0
|
||||
# C1[dofs,:] = 0
|
||||
# C2[dofs,:] = 0
|
||||
# D[dofs,:] = 0
|
||||
# g[dofs,:] = 0
|
||||
#end
|
||||
#for j in actives
|
||||
# dofs = [field_dim*(j-1)+1, field_dim*(j-1)+2]
|
||||
# C2[dofs[1],:] = 0
|
||||
#end
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff_old2}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(x::Vector)
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc = zeros(u)
|
||||
C = zeros(la)
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
# 1. update nodal normals for this element. average nodes from adjacent elements
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
X2 = master_element("geometry", time)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
xi1a = project_from_master_to_slave(slave_element, x1, n1, x2[1])
|
||||
xi1b = project_from_master_to_slave(slave_element, x1, n1, x2[end])
|
||||
xi1 = clamp([xi1a; xi1b], -1.0, 1.0)
|
||||
l = 1/2*abs(xi1[2]-xi1[1])
|
||||
isapprox(l, 0.0) && continue # no contribution in this master element
|
||||
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
|
||||
# project gauss point from slave element to master element
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s))
|
||||
x_s = N1*x1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = Q'*n_s # tangent direction in gauss point
|
||||
R_s = [n_s t_s]
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = vec(get_basis(master_element, xi_m))
|
||||
x_m = N2*x2
|
||||
Phi = Ae*N1
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
|
||||
la_s = Phi*la1 # traction force in gauss point
|
||||
#lan = dot(n_s, la_s) # normal component
|
||||
#lat = dot(t_s, la_s) # tangential component
|
||||
la_nt = R_s*la_s
|
||||
g = x_s-x_m
|
||||
gn = props.gap_sign*dot(n_s, g) # normal gap
|
||||
#gu = dot(n_s, u_s - u_m) # normal displacement gap
|
||||
#gt = dot(t_s, x_s - x_m) # tangential gap
|
||||
|
||||
fc[:,slave_element_nodes] += w*la_s*N1'
|
||||
fc[:,master_element_nodes] -= w*la_s*N2'
|
||||
C[1,slave_element_nodes] += w*gn*Phi'
|
||||
#C[2,slave_element_nodes] += w*la_nt[2]*Phi'
|
||||
#C[2,slave_element_nodes] += w*la_nt[2,:]*Phi'
|
||||
#C[2,slave_element_nodes] += w*dot(t_s, u_s - u_m)*Phi' # <-- for tie
|
||||
|
||||
#R += w*R_s
|
||||
|
||||
end
|
||||
|
||||
end # master elements done
|
||||
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
n = n1[i]
|
||||
t = Q'*n
|
||||
R = [n t]
|
||||
la_nt = R*la[:,j]
|
||||
C[2,j] += la_nt[2]
|
||||
end
|
||||
|
||||
return vec([fc C])
|
||||
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
#if ndofs == 0
|
||||
# info("INITIALIZING THINGS")
|
||||
# problem.assembly.u = zeros(16)
|
||||
# problem.assembly.la = zeros(16)
|
||||
# x = [problem.assembly.u; problem.assembly.la]
|
||||
# ndofs = round(Int, length(x)/2)
|
||||
#end
|
||||
|
||||
A, allresults = ForwardDiff.jacobian(calculate_interface, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
#b = -calculate_interface(x)
|
||||
#info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
#println(A)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(A[1:ndofs,ndofs+1:end])
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = b[1:ndofs]
|
||||
g = b[ndofs+1:end]
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(x::Vector)
|
||||
|
||||
ndofs = round(Int, length(x)/2)
|
||||
nnodes = round(Int, ndofs/field_dim)
|
||||
u = reshape(x[1:ndofs], field_dim, nnodes)
|
||||
la = reshape(x[ndofs+1:end], field_dim, nnodes)
|
||||
fc = zeros(u)
|
||||
C = zeros(la)
|
||||
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
X1 = slave_element("geometry", time)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
# 1. update nodal normals for this element. average nodes from adjacent elements
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[ForwardDiff.get_value(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
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
|
||||
# project gauss point from slave element to master element
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = vec(get_basis(slave_element, xi_s))
|
||||
x_s = N1*x1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = Q'*n_s # tangent direction in gauss point
|
||||
R_s = [n_s t_s]
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = vec(get_basis(master_element, xi_m))
|
||||
x_m = N2*x2
|
||||
Phi = Ae*N1
|
||||
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
|
||||
la_s = Phi*la1 # traction force in gauss point
|
||||
la_nt = R_s*la_s
|
||||
gn = -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'
|
||||
#C[1,slave_element_nodes] += w*gn*Phi'
|
||||
|
||||
lan_tot += w*la_nt[1]*Phi
|
||||
gap_tot += w*gn*Phi
|
||||
|
||||
end
|
||||
|
||||
end # master elements done
|
||||
|
||||
# ncf = lan_tot - max(0, lan_tot - gap_tot)
|
||||
# info("pressure in nodes: $(ForwardDiff.get_value(lan_tot))")
|
||||
# info("weighted gap in nodes: $(ForwardDiff.get_value(gap_tot))")
|
||||
# info("ncf: $(ForwardDiff.get_value(ncf))")
|
||||
# cond = +lan_tot + gap_tot
|
||||
# cond = +lan_tot - gap_tot # singular
|
||||
# cond = -lan_tot + gap_tot
|
||||
# cond = -lan_tot - gap_tot
|
||||
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
n = n1[i]
|
||||
t = Q'*n
|
||||
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]))")
|
||||
|
||||
# if -lan_tot[i] + gap_tot[i] < 0
|
||||
if -la_nt[1] + gap_tot[i] < 0
|
||||
#if j in [31, 32, 33, 34, 35, 36]
|
||||
info("set node $j active")
|
||||
C[1,j] -= gap_tot[i]
|
||||
# C[1,j] += la_nt[1] - max(0, la_nt[1] - gap_tot[i])
|
||||
C[2,j] += la_nt[2]
|
||||
else
|
||||
info("set node $j inactive")
|
||||
C[1,j] += la1[i][1]
|
||||
C[2,j] += la1[i][2]
|
||||
end
|
||||
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)
|
||||
b = -ForwardDiff.value(allresults)
|
||||
#b = -calculate_interface(x)
|
||||
#info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
#println(A)
|
||||
K = A[1:ndofs,1:ndofs]
|
||||
C1 = transpose(A[1:ndofs,ndofs+1:end])
|
||||
C2 = A[ndofs+1:end,1:ndofs]
|
||||
D = A[ndofs+1:end,ndofs+1:end]
|
||||
f = b[1:ndofs]
|
||||
g = b[ndofs+1:end]
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
@@ -0,0 +1,643 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# Mortar projection calculation for 3d cases
|
||||
|
||||
""" Construct auxiliary plane for surface. """
|
||||
function create_auxiliary_plane{E}(element::Element{E}, time::Real)
|
||||
xi = get_reference_element_midpoint(E)
|
||||
x0 = element("geometry", xi, time)
|
||||
ntbasis = element("normal-tangential coordinates", xi, time)
|
||||
return x0, ntbasis
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Project point q onto a plane given by a point p and normal n.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
q::Array{Float64, 2}
|
||||
point to project (row vector)
|
||||
x0::Array{Float64, 2}
|
||||
origo of plane
|
||||
n::Array{Float64, 2}
|
||||
normal vector of plane
|
||||
|
||||
Returns
|
||||
-------
|
||||
y::Array{Float64, 2}
|
||||
projected point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> p = [-0.5 -1.0 4.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> n = [0.1485860 0.0784519 0.9857830]'
|
||||
julia> project_node_to_auxiliary_plane(p, x0, n)
|
||||
3-element Array{Float64,1}:
|
||||
0.963455
|
||||
-1.2447
|
||||
0.925247
|
||||
|
||||
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)
|
||||
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.")
|
||||
end
|
||||
return qproj[2:3]
|
||||
end
|
||||
|
||||
"""
|
||||
Find edge intersections of two planar arbitrary shape polygons.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64,2}
|
||||
M::Array{Float64,2}
|
||||
|
||||
Matrices with size (2, n) where n is number of vertices of each polygon.
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64,2}
|
||||
Intersection points of polygons
|
||||
n::Array{Float64,2}
|
||||
Neighbour info matrix with size (ns, mn). This keeps information which
|
||||
edges of polygons are intersecting. See further explanation in example
|
||||
below.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Find intersection points of two triangles:
|
||||
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 1 3/2]'
|
||||
julia> P, n = get_edge_intersections(S, M)
|
||||
julia> P
|
||||
2x4 Array{Float64,2}:
|
||||
1.0 1.75 0.0 0.0
|
||||
0.0 0.0 0.5 1.25
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 1 0
|
||||
0 0 0
|
||||
1 0 1)
|
||||
|
||||
So intersection points are: (1.00, 0.00), (1.75, 0.00), (0.00, 0.50), (0.00, 1.25).
|
||||
"Neighbour matrix" can be interpreted as following:
|
||||
|
||||
1 1 0 <--> First edge of S intersects edges 1 and 2 of M
|
||||
0 0 0 <--> Second edge of S doesn't intersect at all
|
||||
1 0 1 <--> Third edge of S intersects with edges 1 and 3 of M
|
||||
|
||||
"""
|
||||
function get_edge_intersections(S::Matrix, M::Matrix)
|
||||
ns = size(S, 2)
|
||||
nm = size(M, 2)
|
||||
P = zeros(2, 0)
|
||||
n = zeros(Int64, ns, nm)
|
||||
k = 0
|
||||
for i=1:ns
|
||||
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
|
||||
r = A\b
|
||||
if (r[1]>=0) & (r[1]<=1) & (r[2]>=0) & (r[2]<=1) # intersection found
|
||||
k += 1
|
||||
f = S[:,i]+r[1]*(S[:,mod(i,ns)+1] - S[:,i])
|
||||
f = f''
|
||||
P = hcat(P, f)
|
||||
n[i, j] = 1
|
||||
end
|
||||
end
|
||||
end
|
||||
end
|
||||
return P, n
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Find any points laying inside or border of triangle.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
Y::Array{Float64, 2}
|
||||
Triangle coordinates in 2×3 matrix
|
||||
X::Array{Float64, 2}
|
||||
List of points to test in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
P::Array{Float64, 2}
|
||||
List of points in triangle in 2×m matrix, where m is number of points inside triangle
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0.0 0.0; 3.0 0.0; 0.0 3.0]' # triangle corner points
|
||||
julia> pts = [-1.0 1.0; 2.0 -0.5; 1.0 1.5; 0.5 1.5]' # points to tests
|
||||
julia> points_in_triangle(S, pts)
|
||||
2x2 Array{Float64,2}:
|
||||
1.0 0.5
|
||||
1.5 1.5
|
||||
"""
|
||||
function get_points_inside_triangle(Y::Matrix, X::Matrix)
|
||||
@assert size(Y, 2) == 3 # "Point in TRIANGLE..."
|
||||
P = zeros(2, 0)
|
||||
v0 = Y[:,2] - Y[:,1]
|
||||
v1 = Y[:,3] - Y[:,1] # find interior points of X in Y
|
||||
d00 = (v0'*v0)[1]
|
||||
d01 = (v0'*v1)[1]
|
||||
d11 = (v1'*v1)[1] # using baricentric coordinates
|
||||
id = 1/(d00*d11 - d01*d01)
|
||||
for i=1:size(X, 2)
|
||||
v2 = X[:,i] - Y[:,1]
|
||||
d02 = (v0'*v2)[1]
|
||||
d12 = (v1'*v2)[1]
|
||||
u = (d11*d02-d01*d12)*id
|
||||
v = (d00*d12-d01*d02)*id
|
||||
if (u>=0) & (v>=0) & (u+v<=1) # also include nodes on the boundary
|
||||
P = hcat(P, X[:,i]'')
|
||||
end
|
||||
end
|
||||
return P
|
||||
end
|
||||
|
||||
"""
|
||||
Determine is point P inside or on boudary of polygon X.
|
||||
|
||||
http://paulbourke.net/geometry/polygonmesh/#insidepoly
|
||||
"""
|
||||
function is_point_inside_convex_polygon(P, X)
|
||||
x, y = P
|
||||
for i=1:length(X)
|
||||
x0, y0 = X[i]
|
||||
x1, y1 = X[mod(i, length(X))+1]
|
||||
if (y-y0)*(x1-x0) - (x-x0)*(y1-y0) < 0
|
||||
return false
|
||||
end
|
||||
end
|
||||
return true
|
||||
end
|
||||
|
||||
function get_points_inside_convex_polygon(pts, X)
|
||||
# TODO: Make more readable
|
||||
X2 = [X[:,i] for i=1:size(X,2)]
|
||||
c = filter(P->is_point_inside_convex_polygon(P, X2), [pts[:,i] for i=1:size(pts, 2)])
|
||||
return length(c) == 0 ? zeros(2, 0) : hcat(c...)
|
||||
end
|
||||
|
||||
""" Return unique objects with some given tolerance. This is used in next function
|
||||
because traditional unique() command returns row vectors as non-unique if they
|
||||
differs only a "little".
|
||||
"""
|
||||
function uniquetol(P, dim::Int; args...)
|
||||
@assert dim == 2
|
||||
items = Vector{Float64}[P[:,i] for i=1:size(P,dim)]
|
||||
new_items = Vector{Float64}[]
|
||||
for item in items
|
||||
has_found = false
|
||||
for new_item in new_items
|
||||
if isapprox(item, new_item; args...)
|
||||
has_found = true
|
||||
break
|
||||
end
|
||||
end
|
||||
if !has_found
|
||||
push!(new_items, item)
|
||||
end
|
||||
end
|
||||
return reshape([new_items...;], length(new_items[]), length(new_items))
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Make polygon clipping of shapes S and M.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
S::Array{Float64, 2}
|
||||
M::Array{Float64, 2}
|
||||
Shapes to clip. Needs to be triangles at the moment.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64, 2}, Array{Float64, 2}
|
||||
- Polygon vertices in 2×n matrix, sorted in counter-clockwise order.
|
||||
- 3×3 "neighbouring" matrix, see example.
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> S = [0 0; 3 0; 0 3]'
|
||||
julia> M = [-1 1; 2 -1/2; 2 2]'
|
||||
julia> P, n = clip_polygon(S, M)
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> n
|
||||
3x3 Array{Int64,2}:
|
||||
1 0 1 <- first edge of M ([-1 1; 2 -1/2]') intersects with edges 1 and 3 of S ([0 0; 3 0]' and [0 3; 0 0]')
|
||||
1 1 0 <- second edge of M ([2 -1/2; 2 2]') intersects with edges 1 and 2 of S
|
||||
0 1 1 <- third edge of M ([2 2; -1 1]') intersects with edgse 2 and 3 of S
|
||||
|
||||
"""
|
||||
function clip_polygon(S::Matrix, M::Matrix)
|
||||
P1, neighbours = get_edge_intersections(M, S)
|
||||
#P2 = get_points_inside_triangle(M, S)
|
||||
#P3 = get_points_inside_triangle(S, M)
|
||||
P2 = get_points_inside_convex_polygon(M, S)
|
||||
P3 = get_points_inside_convex_polygon(S, M)
|
||||
# info("polygon clipping: P1 = $P1")
|
||||
# info("polygon clipping: P2 = $P2")
|
||||
# info("polygon clipping: P3 = $P3")
|
||||
# info("hcat P = $P")
|
||||
P = hcat(P1, P2, P3)
|
||||
if length(P) == 0
|
||||
return nothing, nothing
|
||||
end
|
||||
P = uniquetol(P, 2)
|
||||
meanval = mean(P, 2)
|
||||
tmp = P .- meanval
|
||||
angles = atan2(tmp[2,:], tmp[1,:])
|
||||
angles = reshape(angles, length(angles))
|
||||
order = sortperm(angles)
|
||||
return P[:, order], neighbours
|
||||
end
|
||||
|
||||
|
||||
"""
|
||||
Calculate polygon geometric center point
|
||||
|
||||
Parameters
|
||||
----------
|
||||
P::Array{Float64, 2}
|
||||
Polygon vertices in 2×n matrix
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float63, 2}
|
||||
Center point
|
||||
|
||||
Examples
|
||||
--------
|
||||
julia> P
|
||||
2x6 Array{Float64,2}:
|
||||
0.0 1.0 2.0 2.0 1.25 0.0
|
||||
0.5 0.0 0.0 1.0 1.75 1.33333,
|
||||
julia> C = get_polygon_cp(P)
|
||||
2x1 Array{Float64,2}:
|
||||
1.039740
|
||||
0.804701
|
||||
|
||||
"""
|
||||
function calculate_polygon_centerpoint(P::Matrix)
|
||||
n = size(P, 2)
|
||||
A = 0.0
|
||||
for i=1:n
|
||||
A += 1/2*(P[1,i]*P[2,mod(i,n)+1] - P[1,mod(i,n)+1]*P[2,i])
|
||||
end
|
||||
Cx = 0.0
|
||||
Cy = 0.0
|
||||
for i=1:n
|
||||
inext = mod(i, n)+1
|
||||
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]
|
||||
end
|
||||
|
||||
"""
|
||||
Project point from auxiliary plane to parametric surface given by (ξ₁, ξ₂)
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p::Array{Float64,1}
|
||||
point in auxiliary plane, in (n,t1,t2) coordinate system
|
||||
x0::Array{Float64,1}
|
||||
origo of auxiliary plane cs
|
||||
Q::Array{Float64,2}
|
||||
basis of auxiliary plane cs
|
||||
x::Array{Float64,2}
|
||||
surface node coords
|
||||
basis::Array{Float64,2}
|
||||
surface basis functions
|
||||
dbasis::Array{Float64,2}
|
||||
partial derivatives of surface basis functions
|
||||
|
||||
Returns
|
||||
-------
|
||||
Array{Float64,2}
|
||||
solution vector (d, ξ₁, ξ₂) where d is distance to surface
|
||||
|
||||
Examples
|
||||
--------
|
||||
Define surface with node points, basis + dbasis
|
||||
|
||||
julia> xquad = [
|
||||
... -2.5 -2.0 1.0
|
||||
... 2.5 -2.0 0.7
|
||||
... 2.0 2.3 0.0
|
||||
... -2.0 2.0 1.0]'
|
||||
julia> basis(xi) = [
|
||||
... (1-xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1-xi[2])/4
|
||||
... (1+xi[1])(1+xi[2])/4
|
||||
... (1-xi[1])(1+xi[2])/4]
|
||||
julia> dbasis(xi) = [
|
||||
... -(1-xi[2])/4 -(1-xi[1])/4
|
||||
... (1-xi[2])/4 -(1+xi[1])/4
|
||||
... (1+xi[2])/4 (1+xi[1])/4
|
||||
... -(1+xi[2])/4 (1-xi[1])/4]
|
||||
|
||||
We aim to find point p, which we first project to auxiliary plane defined as following
|
||||
julia> p = [-2.5 -2.0 1.0]'
|
||||
julia> x0 = [0.0 0.075 0.675]'
|
||||
julia> Q = [
|
||||
... 0.1485860 0.9888990 0.0000000
|
||||
... 0.0784519 -0.0117877 0.9968480
|
||||
... 0.9857830 -0.1481180 -0.0793325]
|
||||
|
||||
Our projected point is therefore
|
||||
julia> n = Q[:,1] # first component is normal direction
|
||||
julia> ph = project_node_to_auxiliary_plane(p, x0, n)
|
||||
julia> ph = Q'(ph-x0)
|
||||
julia> ph
|
||||
3x1 Array{Float64,2}:
|
||||
1.33264e-7
|
||||
-2.49593
|
||||
-2.09424
|
||||
|
||||
Our point ph is now in auxiliary plane in n,t1,t2 coordinate system. Next we
|
||||
project it back to surface defined by xquad*basis
|
||||
|
||||
julia> theta = project_point_from_plane_to_surface(ph, x0, Q, xquad, basis, dbasis)
|
||||
julia> theta
|
||||
3x1 Array{Float64,2}:
|
||||
-0.213874
|
||||
-0.999999
|
||||
-1.0
|
||||
|
||||
We see that our ξ₁ = ξ₂ = -1 so we found first point of xquad
|
||||
[-2.5 -2.0 1.0]' correctly.
|
||||
|
||||
julia> xquad*basis(theta[2:3])
|
||||
3-element Array{Float64,1}:
|
||||
-2.5
|
||||
-2.0
|
||||
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)
|
||||
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]
|
||||
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
|
||||
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
|
||||
|
||||
typealias MortarElements3D Union{Tri3, Quad4}
|
||||
|
||||
function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar},
|
||||
slave_element::Element{E}, time::Real, ::Type{Val{:total}})
|
||||
assemble!(assembly, problem, slave_element, time, Val{problem.properties.formulation})
|
||||
end
|
||||
|
||||
function assemble!{E<:MortarElements3D}(assembly::Assembly, problem::Problem{Mortar},
|
||||
slave_element::Element{E}, time::Real, ::Type{Val{:total}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
|
||||
props = problem.properties
|
||||
if props.formulation == :Standard && props.normal_condition == :Contact
|
||||
error("for contact choose Dual formulation.""")
|
||||
end
|
||||
|
||||
# create auxiliary plane and project slave nodes to it
|
||||
# x0 = origo, Q = local basis
|
||||
x0, Q = create_auxiliary_plane(slave_element, time)
|
||||
|
||||
# 1. project slave nodes to auxiliary plane
|
||||
Sl = Vector{Float64}[]
|
||||
for p in slave_element("geometry", time)
|
||||
push!(Sl, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
S = hcat(Sl...)
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
|
||||
# if distance between elements is "far enough" cannot expect contact
|
||||
if (props.normal_condition == :Contact) || props.inequality_constraints
|
||||
slave_midpoint = slave_element("geometry", [0.0, 0.0], time)
|
||||
master_midpoint = master_element("geometry", [0.0, 0.0], time)
|
||||
if norm(slave_midpoint - master_midpoint) > props.minimum_distance
|
||||
continue
|
||||
end
|
||||
end
|
||||
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
# 2. project master nodes to auxiliary plane
|
||||
M = Vector{Float64}[]
|
||||
for p in master_element("geometry", time)
|
||||
push!(M, project_point_to_auxiliary_plane(p, x0, Q))
|
||||
end
|
||||
M = hcat(M...)
|
||||
|
||||
# 3. create polygon clipping on auxiliary plane
|
||||
P = nothing
|
||||
neighbours = nothing
|
||||
try
|
||||
P, neighbours = clip_polygon(S, M)
|
||||
catch
|
||||
info("polygon clipping failed")
|
||||
info("S = ")
|
||||
dump(S)
|
||||
info("M = ")
|
||||
dump(M)
|
||||
info("original Sl = ")
|
||||
info(Sl)
|
||||
error("cannot continue")
|
||||
end
|
||||
isa(P, Void) && continue # no clipping
|
||||
|
||||
# shared edge but no shared volume. skipping
|
||||
size(P, 2) < 3 && continue
|
||||
|
||||
C = calculate_polygon_centerpoint(P)
|
||||
npts = size(P, 2) # number of vertices in polygon
|
||||
|
||||
# loop vertices and create temporary integrate cells
|
||||
# TODO: basically when npts == 3 or npts == 4 we could integrate without splitting to cells.
|
||||
nnodes = size(slave_element, 2)
|
||||
C1S3 = zeros(3*nnodes, 3*nnodes)
|
||||
C1M3 = zeros(3*nnodes, 3*nnodes)
|
||||
|
||||
for pnt=1:npts # integration of mortar matrices begin
|
||||
cell = Field(Vector{Float64}[C, P[:,pnt], P[:,mod(pnt,npts)+1]])
|
||||
|
||||
# calculate slave side projection matrix D
|
||||
# construct dual basis
|
||||
Ae = zeros(nnodes, nnodes)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
if problem.properties.formulation == :Dual # Construct dual basis
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
N = get_basis(Tri3, ip.xi)
|
||||
xi = vec(N*cell)
|
||||
theta = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
|
||||
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], cell[j]') for j=1:length(cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
De += wC*diagm(vec(N1))
|
||||
Me += wC*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
end
|
||||
for i=1:field_dim
|
||||
C1S3[i:field_dim:end,i:field_dim:end] += De
|
||||
end
|
||||
|
||||
# Calculate master side projection matrix M
|
||||
for ip in get_integration_points(Tri3, Val{5})
|
||||
# gauss point in auxiliary plane
|
||||
#N = get_basis(E, ip.xi)
|
||||
N = get_basis(Tri3, ip.xi)
|
||||
xi = vec(N*cell) # xi defined in auxilary plane
|
||||
|
||||
# find projection of gauss point to master and slave elements
|
||||
theta1 = project_point_from_plane_to_surface(xi, x0, Q, slave_element, time)
|
||||
theta2 = project_point_from_plane_to_surface(xi, x0, Q, master_element, time)
|
||||
xi_slave = theta1[2:3]
|
||||
xi_master = theta2[2:3]
|
||||
|
||||
# evaluate shape functions values in gauss point and add contribution to matrices
|
||||
N1 = slave_element(xi_slave, time)
|
||||
N2 = master_element(xi_master, time)
|
||||
|
||||
# jacobian determinant on integration cell
|
||||
dNC = get_dbasis(Tri3, ip.xi)
|
||||
JC = sum([kron(dNC[:,j], cell[j]') for j=1:length(cell)])
|
||||
wC = ip.weight*det(JC)
|
||||
|
||||
# extend matrices according to the problem dimension (3)
|
||||
@assert length(slave_dofs) == length(master_dofs)
|
||||
Me = wC*Ae*N1'*N2
|
||||
for k=1:field_dim
|
||||
C1M3[k:field_dim:end,k:field_dim:end] += Me
|
||||
end
|
||||
end
|
||||
end # integration of mortar matrices done.
|
||||
|
||||
# constraints in normal-tangential direction and initial weighted gap
|
||||
X1 = vec(slave_element("geometry", time))
|
||||
X2 = vec(master_element("geometry", time))
|
||||
Q_ = slave_element("normal-tangential coordinates", time)
|
||||
Z = zeros(3, 3)
|
||||
if nnodes == 3
|
||||
Q3 = [Q Z Z; Z Q Z; Z Z Q]
|
||||
elseif nnodes == 4
|
||||
Q3 = [Q Z Z Z; Z Q Z Z; Z Z Q Z; Z Z Z Q]
|
||||
end
|
||||
D3 = zeros(3*nnodes, 3*nnodes)
|
||||
C2S3 = Q3'*C1S3
|
||||
C2M3 = Q3'*C1M3
|
||||
G = -(C2S3*X1 - C2M3*X2)
|
||||
|
||||
# complementarity condition
|
||||
if haskey(slave_element, "displacement")
|
||||
u1 = vec(slave_element("displacement", time))
|
||||
else
|
||||
u1 = zeros(3*nnodes)
|
||||
end
|
||||
if haskey(master_element, "displacement")
|
||||
u2 = vec(master_element("displacement", time))
|
||||
else
|
||||
u2 = zeros(3*nnodes)
|
||||
end
|
||||
x1 = X1 + u1
|
||||
x2 = X2 + u2
|
||||
if haskey(slave_element, "reaction force")
|
||||
la = vec(slave_element("reaction force", time))
|
||||
else
|
||||
la = zeros(3*nnodes)
|
||||
end
|
||||
g = -(C2S3*x1 - C2M3*x2)
|
||||
c = Q3'*la - g
|
||||
inactive_nodes = find(c[1:field_dim:end] .<= 0)
|
||||
active_nodes = find(c[1:field_dim:end] .> 0)
|
||||
|
||||
# normal constraint: remove inactive nodes if normal condition is set to contact
|
||||
if problem.properties.normal_condition == :Contact
|
||||
for j in inactive_nodes
|
||||
dofs = [3*(j-1)+1, 3*(j-1)+2, 3*(j-1)+3]
|
||||
G[dofs] = 0
|
||||
C1S3[dofs,:] = 0
|
||||
C1M3[dofs,:] = 0
|
||||
C2S3[dofs,:] = 0
|
||||
C2M3[dofs,:] = 0
|
||||
end
|
||||
end
|
||||
|
||||
# tangential constraint: stick or slip
|
||||
if problem.properties.tangential_condition == :Slip
|
||||
D3 = copy(C2S3)
|
||||
D3[1:field_dim:end, :] = 0
|
||||
C2S3[2:field_dim:end, :] = 0
|
||||
C2M3[2:field_dim:end, :] = 0
|
||||
C2S3[3:field_dim:end, :] = 0
|
||||
C2M3[3:field_dim:end, :] = 0
|
||||
end
|
||||
|
||||
# add contributions
|
||||
add!(assembly.C1, slave_dofs, slave_dofs, C1S3)
|
||||
add!(assembly.C1, slave_dofs, master_dofs, -C1M3)
|
||||
add!(assembly.C2, slave_dofs, slave_dofs, C2S3)
|
||||
add!(assembly.C2, slave_dofs, master_dofs, -C2M3)
|
||||
add!(assembly.D, slave_dofs, slave_dofs, D3)
|
||||
add!(assembly.c, slave_dofs, c)
|
||||
add!(assembly.g, slave_dofs, G)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
@@ -0,0 +1,446 @@
|
||||
using JuliaFEM.Core: MortarElements2D, DVTI, Assembly
|
||||
|
||||
import JuliaFEM.Core: project_from_master_to_slave, project_from_slave_to_master, assemble!,
|
||||
get_unknown_field_dimension, get_parent_field_name, get_gdofs, find_elements, get_nodes, Field,
|
||||
get_integration_points, get_basis, get_dbasis, add!
|
||||
|
||||
""" Find segment from slave element corresponding to master element nodes.
|
||||
x1_, n1_
|
||||
slave element geometry and normal direction
|
||||
|
||||
x2_ master element nodes to project onto slave
|
||||
"""
|
||||
function project_from_master_to_slave{E<:MortarElements2D}(
|
||||
slave_element::Element{E}, x1_::DVTI, n1_::DVTI, x2::Vector)
|
||||
|
||||
function x1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*x1_
|
||||
end
|
||||
|
||||
function dx1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*x1_
|
||||
end
|
||||
|
||||
function n1(xi1)
|
||||
N = get_basis(E, xi1)
|
||||
return vec(N)*n1_
|
||||
end
|
||||
|
||||
function dn1(xi1)
|
||||
dN = get_dbasis(E, xi1)
|
||||
return vec(dN)*n1_
|
||||
end
|
||||
|
||||
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
|
||||
for i=1:5
|
||||
dxi1 = -R(xi1)/dR(xi1)
|
||||
xi1 += dxi1
|
||||
if norm(dxi1) < 1.0e-10
|
||||
return xi1
|
||||
end
|
||||
end
|
||||
|
||||
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)
|
||||
|
||||
function x2(xi2)
|
||||
N = get_basis(E, xi2)
|
||||
return vec(N)*x2_
|
||||
end
|
||||
|
||||
function dx2(xi2)
|
||||
dN = get_dbasis(E, xi2)
|
||||
return vec(dN)*x2_
|
||||
end
|
||||
|
||||
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:5
|
||||
dxi2 = -R(xi2) / dR(xi2)
|
||||
xi2 += dxi2
|
||||
if norm(dxi2) < 1.0e-10
|
||||
return xi2
|
||||
end
|
||||
end
|
||||
|
||||
error("find projection from slave to master: did not converge, last val: $xi2 and $dxi2")
|
||||
|
||||
end
|
||||
|
||||
|
||||
function assemble!{E<:MortarElements2D}(assembly::Assembly,
|
||||
problem::Problem{Mortar}, slave_element::Element{E},
|
||||
time::Real, ::Type{Val{:forwarddiff}})
|
||||
haskey(slave_element, "master elements") || return
|
||||
props = problem.properties
|
||||
field_dim = get_unknown_field_dimension(problem)
|
||||
field_name = get_parent_field_name(problem)
|
||||
|
||||
function calculate_interface(u::Matrix, la::Matrix)
|
||||
|
||||
X1 = slave_element("geometry", time)
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
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
|
||||
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
# 1. update nodal normals for this element
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
fc = SparseMatrixCOO{Real}([], [], []) # interface virtual work
|
||||
C = SparseMatrixCOO{Real}([], [], []) # constraints
|
||||
B = SparseMatrixCOO{Real}([], [], [])
|
||||
|
||||
#info("u1.data = ", ForwardDiff.get_value(u1.data))
|
||||
info("normal calculations done. looping master elements.")
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
u2 = Field(Vector[u[:,i] for i in master_element_nodes])
|
||||
x2 = X2 + u2
|
||||
info("master element ready.")
|
||||
|
||||
# calculate segmentation: we care only about endpoints
|
||||
# note: these are quadratic/cubic functions, analytical solution possible
|
||||
info("calculating segmentation.")
|
||||
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
|
||||
info("xi1 = $xi1")
|
||||
|
||||
info("create bi-orthogonal basis")
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
info("bi-orthogonal basis done. integrating fc.")
|
||||
|
||||
slave_dofs = get_gdofs(slave_element, field_dim)
|
||||
master_dofs = get_gdofs(master_element, field_dim)
|
||||
|
||||
info("integrate fc")
|
||||
D = zeros(nnodes, nnodes)
|
||||
M = zeros(nnodes, nnodes)
|
||||
gn = zeros(nnodes)
|
||||
lan = zeros(nnodes)
|
||||
lat = zeros(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)
|
||||
# project gauss point to master element to evaluate shape function there
|
||||
x_s = vec(N1)*x1 # coordinate in gauss point
|
||||
n_s = vec(N1)*n1 # normal direction in gauss point
|
||||
xi_m = project_from_slave_to_master(master_element, x_s, n_s, x2)
|
||||
N2 = get_basis(master_element, xi_m)
|
||||
x_m = vec(N2)*x2
|
||||
Phi = vec(Ae*N1')
|
||||
la_s = Phi*la1 # traction force in gauss point
|
||||
u_s = vec(N1)*u1
|
||||
u_m = vec(N2)*u2
|
||||
#info("la_s = $(ForwardDiff.get_value(la_s))")
|
||||
SM = [slave_dofs; master_dofs]
|
||||
#N1N2 = [N1 -N2]
|
||||
#info("all_dofs = $(SM)")
|
||||
#info("shape functions = $(ForwardDiff.get_value(N1N2))")
|
||||
#for i=1:field_dim
|
||||
# #info("add to slave dofs $(slave_dofs[i:field_dim:end])")
|
||||
# #info("add to master dofs $(master_dofs[i:field_dim:end])")
|
||||
# add!(fc, slave_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*N1)
|
||||
# add!(fc, master_dofs[i:field_dim:end], [1, 1], +w*la_s[i]*N2)
|
||||
#add!(fc, slave_dofs[i:field_dim:end], [1, 1], w*la_s'*u_s[i])
|
||||
#add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s'*u_m[i])
|
||||
# add!(fc, master_dofs[i:field_dim:end], [1, 1], -w*la_s[i]*u_m)
|
||||
#end
|
||||
#add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1, 1, 1, 1, 1], w*la_s*[u_s' -u_m'])
|
||||
D += w*kron(Ae*N1', N1)
|
||||
M += w*kron(Ae*N1', N2)
|
||||
gn += -w*dot(n_s, x_s-x_m)*Phi
|
||||
lan += w*dot(n_s, la_s)*Phi
|
||||
t_s = Q'*n_s
|
||||
lat += w*dot(t_s, la_s)*Phi
|
||||
end
|
||||
|
||||
#D2 = zeros(2*nnodes, 2*nnodes)
|
||||
#M2 = zeros(2*nnodes, 2*nnodes)
|
||||
#for i=1:field_dim
|
||||
# D2[i:field_dim:end, i:field_dim:end] += D
|
||||
# M2[i:field_dim:end, i:field_dim:end] += M
|
||||
#end
|
||||
#info("size of D2 = $(size(D2))")
|
||||
#fco = [D2 -M2]*vec(la1)
|
||||
#fco = [D -M]*la[:,slave_element_nodes]
|
||||
#info("fco = $(ForwardDiff.get_value(fco))")
|
||||
#add!(fc, [slave_dofs; master_dofs], [1, 1, 1, 1], fco)
|
||||
|
||||
for i=1:field_dim
|
||||
add!(B, slave_dofs[i:field_dim:end], slave_dofs[i:field_dim:end], D)
|
||||
add!(B, slave_dofs[i:field_dim:end], master_dofs[i:field_dim:end], -M)
|
||||
end
|
||||
info("gn = $gn")
|
||||
#Cj = lan - max(0, lan - gn) + lat
|
||||
add!(C, slave_dofs[1:field_dim:end], [1, 1], gn')
|
||||
|
||||
end # master elements done
|
||||
|
||||
ndofs = prod(size(la))
|
||||
N = SparseMatrixCOO{Real}([], [], [])
|
||||
T = SparseMatrixCOO{Real}([], [], [])
|
||||
for (i, j) in enumerate(slave_element_nodes)
|
||||
dofs = [2*(j-1)+1, 2*(j-1)+2]
|
||||
add!(N, [dofs[1]], dofs, reshape(n1[i], 1, 2))
|
||||
add!(T, [dofs[2]], dofs, reshape(Q'*n1[i], 1, 2))
|
||||
end
|
||||
N = sparse(N, ndofs, ndofs)
|
||||
T = sparse(T, ndofs, ndofs)
|
||||
B = sparse(B, ndofs, ndofs)
|
||||
fc = B'*vec(la)
|
||||
#println(sparse(fc))
|
||||
#fc = sparse(fc, ndofs, 1)
|
||||
#println(fc)
|
||||
#dump(full(fc))
|
||||
#C = sparse(C, ndofs, 1)
|
||||
C = N*B*vec(u) + T*vec(la)
|
||||
return fc, C
|
||||
|
||||
end
|
||||
|
||||
|
||||
function calculate_interface_PE(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)
|
||||
#fixed_la = ForwardDiff.get_value(la)
|
||||
#fixed_u = ForwardDiff.get_value(u)
|
||||
#u = ForwardDiff.get_value(u)
|
||||
|
||||
X1 = slave_element("geometry", time)
|
||||
slave_element_nodes = get_connectivity(slave_element)
|
||||
u1 = Field(Vector[u[:,i] for i in slave_element_nodes])
|
||||
la1 = Field(Vector[la[:,i] for i in slave_element_nodes])
|
||||
#fixed_la1 = Field(Vector[fixed_la[:,i] for i in slave_element_nodes])
|
||||
x1 = X1 + u1
|
||||
|
||||
# 1. update nodal normals for this element
|
||||
adjacent_elements = find_elements(get_elements(problem), slave_element_nodes)
|
||||
adjacent_nodes = get_nodes(adjacent_elements) # including also nodes from adjacent elements
|
||||
Q = [0.0 -1.0; 1.0 0.0]
|
||||
normals = zeros(u)
|
||||
for element in adjacent_elements
|
||||
conn = get_connectivity(element)
|
||||
gdofs = get_gdofs(element, field_dim)
|
||||
X_el = element("geometry", time)
|
||||
u_el = Field(Vector[u[:, i] for i in conn])
|
||||
x_el = X_el + u_el
|
||||
for ip in get_integration_points(element, Val{3})
|
||||
dN = get_dbasis(element, ip)
|
||||
N = element(ip, time)
|
||||
t = sum([kron(dN[:,i], x_el[i]') for i=1:length(x_el)])
|
||||
normals[:, conn] += ip.weight*Q*t'*N
|
||||
end
|
||||
end
|
||||
# --> slave side normals in deformed state
|
||||
n1 = Field(Vector[normals[:,i]/norm(normals[:,i]) for i in slave_element_nodes])
|
||||
|
||||
Wco = 0.0
|
||||
Wla = 0.0
|
||||
|
||||
for master_element in slave_element["master elements"]
|
||||
X2 = master_element("geometry", time)
|
||||
master_element_nodes = get_connectivity(master_element)
|
||||
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
|
||||
|
||||
nnodes = size(slave_element, 2)
|
||||
De = zeros(nnodes, nnodes)
|
||||
Me = zeros(nnodes, nnodes)
|
||||
for ip in get_integration_points(slave_element, Val{5})
|
||||
# jacobian of slave element in deformed state
|
||||
dN = get_dbasis(slave_element, ip)
|
||||
j = sum([kron(dN[:,i], x1[i]') for i=1:length(x1)])
|
||||
w = ip.weight*norm(j)*l
|
||||
xi_s = dot([1/2*(1-ip.xi); 1/2*(1+ip.xi)], xi1)
|
||||
N1 = get_basis(slave_element, xi_s)
|
||||
De += w*diagm(vec(N1))
|
||||
Me += w*N1'*N1
|
||||
end
|
||||
Ae = De*inv(Me)
|
||||
|
||||
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 = vec(get_basis(slave_element, xi_s))
|
||||
|
||||
# project gauss point to master element to evaluate shape function there
|
||||
x_s = N1*x1 # coordinate in gauss point
|
||||
n_s = N1*n1 # normal direction in gauss point
|
||||
t_s = Q'*n_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
|
||||
|
||||
gn = -dot(n_s, x_s - x_m)
|
||||
gt = dot(t_s, x_s - x_m)
|
||||
lan = dot(n_s, Phi*la1)
|
||||
lat = dot(t_s, Phi*la1)
|
||||
u_s = N1*u1
|
||||
u_m = N2*u2
|
||||
gu = dot(n_s, u_s - u_m)
|
||||
|
||||
Wco += w*dot(Phi*la1, N1*u1 - N2*u2)
|
||||
#Wla += w*(lan*gn + lat*gt)
|
||||
#gn = min(0, gn)
|
||||
Wla += 1/2*w*1e6*gn*gn
|
||||
#info("gn = $(ForwardDiff.get_value(gn))")
|
||||
#Wla += w*dot(dot(n_s, Phi*la1), dot(n_s, N1*u1 - N2*u2))
|
||||
end
|
||||
|
||||
end
|
||||
return Wco, Wla
|
||||
end
|
||||
|
||||
|
||||
function calculate_contact_rhs(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, C = calculate_interface(u, la)
|
||||
info("interface vector calculated.")
|
||||
return vec(full([fc; C]))
|
||||
end
|
||||
|
||||
# x doesn't mean deformed configuration here
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
if ndofs == 0
|
||||
info("INITIALIZING THINGS")
|
||||
problem.assembly.u = zeros(16)
|
||||
problem.assembly.la = zeros(16)
|
||||
x = [problem.assembly.u; problem.assembly.la]
|
||||
ndofs = round(Int, length(x)/2)
|
||||
end
|
||||
|
||||
function add_fco!()
|
||||
get_PI(x::Vector) = calculate_interface_PE(x)[1]
|
||||
A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.gradient(allresults)
|
||||
info("PE = $(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]
|
||||
add!(assembly.K, K)
|
||||
add!(assembly.C1, C1)
|
||||
add!(assembly.C2, C2)
|
||||
add!(assembly.D, D)
|
||||
add!(assembly.f, f)
|
||||
add!(assembly.g, g)
|
||||
end
|
||||
#add_fco!()
|
||||
|
||||
function add_wla!()
|
||||
get_PI(x::Vector) = calculate_interface_PE(x)[2]
|
||||
A, allresults = ForwardDiff.hessian(get_PI, x, ForwardDiff.AllResults)
|
||||
b = -ForwardDiff.gradient(allresults)
|
||||
info("PE = $(ForwardDiff.value(allresults))")
|
||||
A = sparse(A)
|
||||
b = sparse(b)
|
||||
SparseMatrix.droptol!(A, 1.0e-12)
|
||||
SparseMatrix.droptol!(b, 1.0e-12)
|
||||
#info("A")
|
||||
#println(full(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)
|
||||
end
|
||||
add_wla!()
|
||||
|
||||
return
|
||||
|
||||
end
|
||||
|
||||
mesh, body1, body2, bc_top, bc_bottom, contact = divided_block_problem()
|
||||
bc_top.properties.formulation = :incremental
|
||||
bc_bottom.properties.formulation = :incremental
|
||||
contact.properties.formulation = :forwarddiff
|
||||
contact.assembly.u = zeros(16)
|
||||
contact.assembly.la = zeros(16)
|
||||
assemble!(contact.assembly, contact, contact.elements[1], 0.0, Val{:forwarddiff})
|
||||
|
||||
+13
-7
@@ -166,16 +166,22 @@ function update_assembly!(problem, u, la)
|
||||
# copy current solutions to previous ones and add/replace new solution
|
||||
assembly.u_prev = copy(assembly.u)
|
||||
assembly.la_prev = copy(assembly.la)
|
||||
if get_formulation_type(problem) == :incremental
|
||||
info("$(problem.name): incremental formulation, adding increment to solution vector")
|
||||
#info("solution vector:")
|
||||
#dump(round(u, 3)')
|
||||
assembly.u += u
|
||||
else
|
||||
if get_formulation_type(problem) == :total
|
||||
info("$(problem.name): total formulation, replacing solution vector with new values")
|
||||
assembly.u = u
|
||||
assembly.la = la
|
||||
elseif get_formulation_type(problem) == :incremental
|
||||
info("$(problem.name): incremental formulation, adding increment to solution vector")
|
||||
assembly.u += u
|
||||
assembly.la = la
|
||||
elseif get_formulation_type(problem) == :forwarddiff
|
||||
info("$(problem.name): forwarddiff formulation, adding increment to solution vector")
|
||||
assembly.u += u
|
||||
assembly.la += la
|
||||
else
|
||||
info("$(problem.name): unknown formulation type, don't know what to do with results")
|
||||
error("serious failure with problem formulation: $(get_formulation_type(problem))")
|
||||
end
|
||||
assembly.la = la
|
||||
|
||||
# calculate change of norm
|
||||
assembly.u_norm_change = norm(assembly.u - assembly.u_prev)
|
||||
|
||||
+17
-4
@@ -141,9 +141,11 @@ type Solver
|
||||
name :: ASCIIString # some descriptive name for problem
|
||||
time :: Real # current time
|
||||
iteration :: Int # iteration counter
|
||||
norms :: Vector{Tuple} # solution norms for convergence studies
|
||||
ndofs :: Int # total dimension of global stiffness matrix, i.e., dim*nnodes
|
||||
problems :: Vector{Problem}
|
||||
is_linear_system :: Bool # setting this to true makes assumption of one step convergence
|
||||
nonlinear_system_min_iterations :: Int64
|
||||
nonlinear_system_max_iterations :: Int64
|
||||
nonlinear_system_convergence_tolerance :: Float64
|
||||
nonlinear_system_error_if_no_convergence :: Bool
|
||||
@@ -155,9 +157,11 @@ function Solver(name::ASCIIString="default solver", time::Real=0.0)
|
||||
name,
|
||||
time,
|
||||
0, # iteration #
|
||||
[], # solution norms in (norm(u), norm(la)) tuples
|
||||
0, # ndofs
|
||||
[], # array of problems
|
||||
false, # is_linear_system
|
||||
1, # min nonlinear iterations
|
||||
10, # max nonlinear iterations
|
||||
5.0e-5, # nonlinear iteration convergence tolerance
|
||||
true, # throw error if no convergence
|
||||
@@ -277,12 +281,14 @@ crosspoints.
|
||||
function get_boundary_assembly(solver::Solver)
|
||||
ndofs = solver.ndofs
|
||||
@assert ndofs != 0
|
||||
Kc = spzeros(ndofs, ndofs)
|
||||
C1 = spzeros(ndofs, ndofs)
|
||||
C2 = spzeros(ndofs, ndofs)
|
||||
D = spzeros(ndofs, ndofs)
|
||||
g = spzeros(ndofs, 1)
|
||||
for problem in get_boundary_problems(solver)
|
||||
assembly = problem.assembly
|
||||
Kc_ = sparse(assembly.K, ndofs, ndofs)
|
||||
C1_ = sparse(assembly.C1, ndofs, ndofs)
|
||||
C2_ = sparse(assembly.C2, ndofs, ndofs)
|
||||
D_ = sparse(assembly.D, ndofs, ndofs)
|
||||
@@ -297,12 +303,13 @@ function get_boundary_assembly(solver::Solver)
|
||||
handle_overconstraint_error!(problem, overconstrained_nodes,
|
||||
overconstrained_dofs, C1, C1_, C2, C2_, D, D_, g, g_)
|
||||
end
|
||||
Kc += Kc_
|
||||
C1 += C1_
|
||||
C2 += C2_
|
||||
D += D_
|
||||
g += g_
|
||||
end
|
||||
return C1, C2, D, g
|
||||
return Kc, C1, C2, D, g
|
||||
end
|
||||
|
||||
|
||||
@@ -316,10 +323,10 @@ function solve_linear_system(solver::Solver, ::Type{Val{:DirectLinearSolver}})
|
||||
K, f = get_field_assembly(solver)
|
||||
|
||||
# assemble boundary problems
|
||||
C1, C2, D, g = get_boundary_assembly(solver)
|
||||
Kc, C1, C2, D, g = get_boundary_assembly(solver)
|
||||
|
||||
# construct global system Ax=b and solve using lu factorization
|
||||
A = [K C1'; C2 D]
|
||||
A = [K+Kc C1'; C2 D]
|
||||
b = [f; g]
|
||||
|
||||
nz = get_nonzero_rows(A)
|
||||
@@ -379,6 +386,7 @@ end
|
||||
""" Main solver loop.
|
||||
"""
|
||||
function call(solver::Solver)
|
||||
|
||||
# 1. initialize each problem so that we can start nonlinear iterations
|
||||
for problem in solver.problems
|
||||
initialize!(problem, solver.time)
|
||||
@@ -394,6 +402,7 @@ function call(solver::Solver)
|
||||
|
||||
# 2.2 call solver for linearized system (default: direct lu factorization)
|
||||
u, la = solve_linear_system(solver, Val{solver.linear_system_solver})
|
||||
push!(solver.norms, (norm(u), norm(la)))
|
||||
|
||||
# 2.3 update solution back to elements
|
||||
for problem in solver.problems
|
||||
@@ -404,7 +413,11 @@ function call(solver::Solver)
|
||||
# 2.4 check convergence
|
||||
if has_converged(solver)
|
||||
info("Converged in $(solver.iteration) iterations.")
|
||||
return true
|
||||
if solver.iteration < solver.nonlinear_system_min_iterations
|
||||
info("Converged but continuing")
|
||||
else
|
||||
return true
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
+17
-8
@@ -4,18 +4,22 @@
|
||||
# Sparse utils to make assembly of local and global matrices easier.
|
||||
# Unoptimized but should do all necessary stuff for at start.
|
||||
|
||||
type SparseMatrixCOO
|
||||
type SparseMatrixCOO{T<:Real}
|
||||
I :: Vector{Int}
|
||||
J :: Vector{Int}
|
||||
V :: Vector{Float64}
|
||||
V :: Vector{T}
|
||||
end
|
||||
|
||||
typealias SparseMatrixIJV SparseMatrixCOO
|
||||
|
||||
function SparseMatrixCOO()
|
||||
SparseMatrixCOO([], [], [])
|
||||
SparseMatrixCOO{Float64}([], [], [])
|
||||
end
|
||||
|
||||
#function SparseMatrixCOO{T}()
|
||||
# SparseMatrixCOO{T}([], [], [])
|
||||
#end
|
||||
|
||||
function Base.convert(::Type{SparseMatrixCOO}, A::SparseMatrixCSC)
|
||||
return SparseMatrixCOO(findnz(A)...)
|
||||
end
|
||||
@@ -80,7 +84,7 @@ function Base.(:+)(A::SparseMatrixIJV, B::SparseMatrixIJV)
|
||||
return C
|
||||
end
|
||||
|
||||
function Base.full(A::SparseMatrixIJV, args...)
|
||||
function Base.full(A::SparseMatrixCOO, args...)
|
||||
return full(sparse(A.I, A.J, A.V, args...))
|
||||
end
|
||||
|
||||
@@ -94,7 +98,7 @@ Example
|
||||
>>> S = [3, 4]
|
||||
>>> M = [6, 7, 8]
|
||||
>>> data = Float64[5 6 7; 8 9 10]
|
||||
>>> A = SparseMatrixIJV()
|
||||
>>> A = SparseMatrixCOO()
|
||||
>>> add!(A, S, M, data)
|
||||
>>> full(A)
|
||||
4x8 Array{Float64,2}:
|
||||
@@ -104,7 +108,7 @@ Example
|
||||
0.0 0.0 0.0 0.0 0.0 8.0 9.0 10.0
|
||||
|
||||
"""
|
||||
function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix{Float64})
|
||||
function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data::Matrix)
|
||||
n, m = size(data)
|
||||
for j=1:m
|
||||
for i=1:n
|
||||
@@ -112,11 +116,16 @@ function add!(A::SparseMatrixIJV, dofs1::Vector{Int}, dofs2::Vector{Int}, data::
|
||||
push!(A.J, dofs2[j])
|
||||
end
|
||||
end
|
||||
# append!(A.I, repeat(dofs1, outer=[m]))
|
||||
# append!(A.J, repeat(dofs2, inner=[n]))
|
||||
append!(A.V, vec(data))
|
||||
end
|
||||
|
||||
""" Add sparse matrix of CSC to COO. """
|
||||
function add!(A::SparseMatrixCOO, B::SparseMatrixCSC)
|
||||
I, J, V = findnz(B)
|
||||
C = SparseMatrixCOO(I, J, V)
|
||||
append!(A, C)
|
||||
end
|
||||
|
||||
""" Add new data to COO Sparse vector. """
|
||||
function add!(A::SparseMatrixCOO, dofs::Vector{Int}, data::Array{Float64}, dim::Int=1)
|
||||
if length(dofs) != length(data)
|
||||
|
||||
+2
-2
@@ -73,11 +73,11 @@ function xdmf_new_grid(temporal_collection; time=0)
|
||||
return grid
|
||||
end
|
||||
|
||||
function xdmf_new_mesh!(grid, nodes, elements)
|
||||
function xdmf_new_mesh!(grid, nodes, elements; datatype="XYZ")
|
||||
|
||||
# 1. write nodes
|
||||
geometry = new_child(grid, "Geometry")
|
||||
set_attribute(geometry, "Type", "XYZ")
|
||||
set_attribute(geometry, "Type", datatype)
|
||||
dataitem = new_child(geometry, "DataItem")
|
||||
set_attribute(dataitem, "DataType", "Float")
|
||||
ndim = sum([length(node) for node in nodes])
|
||||
|
||||
Reference in New Issue
Block a user