-
Notifications
You must be signed in to change notification settings - Fork 7
/
calculate_projections.jl
53 lines (44 loc) · 2.46 KB
/
calculate_projections.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE
"""
project_from_slave_to_master(Val{:Seg2}, xs, ns, xm1, xm2)
Find the projection of a slave node `xs`, having normal vector `ns`, onto master
elements with nodes (`xm1`, `xm2`).
"""
function project_from_slave_to_master{T<:Number}(::Type{Val{:Seg2}}, xs::Vector{T},
ns::Vector{T}, xm1::Vector{T},
xm2::Vector{T})
nom = ns[1]*(xm1[2] + xm2[2] - 2*xs[2]) - ns[2]*(xm1[1] + xm2[1] - 2*xs[1])
denom = ns[1]*(xm1[2] - xm2[2]) + ns[2]*(xm2[1] - xm1[1])
return nom/denom
end
"""
project_from_master_to_slave(Val{:Seg2}, xm, xs1, xs2, ns1, ns2)
Find the projection of a master node `xm`, to the slave surface with nodes
(`xs1`, `xs2`), in direction of slave surface normal defined by (`ns1`, `ns2`).
"""
function project_from_master_to_slave{T<:Number}(::Type{Val{:Seg2}}, xm::Vector{T},
xs1::Vector{T}, xs2::Vector{T},
ns1::Vector{T}, ns2::Vector{T})
# Special case when normal is constant. Then we have
# unique solution and linear equation to solve.
if isapprox(ns1, ns2)
nom = -2*ns1[1]*xm[2] + ns1[1]*xs1[2] + ns1[1]*xs2[2] + 2*ns1[2]*xm[1] - ns1[2]*xs1[1] - ns1[2]*xs2[1]
denom = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1]
return nom/denom
end
nom_b = 2*ns1[1]*xm[2] - 2*ns1[1]*xs1[2] - 2*ns1[2]*xm[1] + 2*ns1[2]*xs1[1] - 2*ns2[1]*xm[2] + 2*ns2[1]*xs2[2] + 2*ns2[2]*xm[1] - 2*ns2[2]*xs2[1]
denom_b = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1] - ns2[1]*xs1[2] + ns2[1]*xs2[2] + ns2[2]*xs1[1] - ns2[2]*xs2[1]
nom_c = -2*ns1[1]*xm[2] + ns1[1]*xs1[2] + ns1[1]*xs2[2] + 2*ns1[2]*xm[1] - ns1[2]*xs1[1] - ns1[2]*xs2[1] - 2*ns2[1]*xm[2] + ns2[1]*xs1[2] + ns2[1]*xs2[2] + 2*ns2[2]*xm[1] - ns2[2]*xs1[1] - ns2[2]*xs2[1]
denom_c = ns1[1]*xs1[2] - ns1[1]*xs2[2] - ns1[2]*xs1[1] + ns1[2]*xs2[1] - ns2[1]*xs1[2] + ns2[1]*xs2[2] + ns2[2]*xs1[1] - ns2[2]*xs2[1]
a = 1.0
b = nom_b / denom_b
c = nom_c / denom_c
d = b^2 - 4*a*c
if d < 0
warn("Mortar2D.calculate_projections(): negative discriminant $d")
warn("xm = $xm, xs1 = $xs1, xs2 = $xs2, ns1 = $ns1, ns2 = $ns2")
end
sols = [-b + sqrt(d), -b - sqrt(d)]/(2.0*a)
return sols[indmin(abs.(sols))]
end