mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
separate optomechanical experiments to own package
This commit is contained in:
Jukka Aho
@@ -120,9 +120,6 @@ export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
|
||||
include("solvers_modal.jl")
|
||||
export Modal
|
||||
|
||||
include("optics.jl")
|
||||
export find_intersection, calc_reflection, calc_normal
|
||||
|
||||
### Mortar methods, contact mechanics extension ###
|
||||
include("problems_contact.jl")
|
||||
include("problems_contact_2d.jl")
|
||||
|
||||
-128
@@ -1,128 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
function gen_rtic_grid!{S<:Union{Seg2,Seg3}}(thetas, element::Element{S}; npts=2)
|
||||
xi1 = midpoints(linspace(-1, 1, npts+1))
|
||||
for xi in xi1
|
||||
push!(thetas, [0.0, xi])
|
||||
end
|
||||
end
|
||||
|
||||
function gen_rtic_grid!(thetas, element::Element{NSeg}; npts=2)
|
||||
knots_u = element.properties.knots
|
||||
u = midpoints(linspace(minimum(knots_u), maximum(knots_u), npts+1))
|
||||
for ui in u
|
||||
push!(thetas, [0.0, ui])
|
||||
end
|
||||
end
|
||||
|
||||
function gen_rtic_grid!(thetas, element::Element{NSurf}; npts=2)
|
||||
knots_u = element.properties.knots_u
|
||||
knots_v = element.properties.knots_v
|
||||
u = midpoints(linspace(minimum(knots_u), maximum(knots_u), npts+1))
|
||||
v = midpoints(linspace(minimum(knots_v), maximum(knots_v), npts+1))
|
||||
for ui in u
|
||||
for vi in v
|
||||
push!(thetas, [0.0, ui, vi])
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
"""
|
||||
Find intersection between element surface X(ξ) and ray x(t) = s + t*n by
|
||||
solving equation F(t, ξ) = x(t) - X(ξ) = 0
|
||||
|
||||
Returns
|
||||
-------
|
||||
t, ξ
|
||||
|
||||
References
|
||||
----------
|
||||
[1] https://en.wikipedia.org/wiki/Ray_tracing_%28graphics%29
|
||||
|
||||
"""
|
||||
function find_intersection{S<:Union{Seg2, Seg3, NSeg, NSurf}}(element::Element{S}, s, n, time;
|
||||
max_iterations=10, tolerance=1.0e-12, info_output=false, deformed=true, secondary=false, npts=2)
|
||||
|
||||
x = element["geometry"](time)
|
||||
if deformed && haskey(element, "displacement")
|
||||
x += element["displacement"](time)
|
||||
end
|
||||
|
||||
function calc_intersection!(theta)
|
||||
i = 0
|
||||
dtheta = zeros(theta)
|
||||
for i=1:max_iterations
|
||||
t = theta[1]
|
||||
xi = theta[2:end]
|
||||
N = get_basis(element, xi, time)
|
||||
dN = get_dbasis(element, xi, time)
|
||||
A = [n -dN*x]
|
||||
r = s + t*n - N*x
|
||||
try
|
||||
dtheta = A \ -r
|
||||
catch err
|
||||
info("failed to solve equation.")
|
||||
dump(theta)
|
||||
dump(A)
|
||||
dump(r)
|
||||
throw(err)
|
||||
end
|
||||
theta[:] += dtheta
|
||||
info_output && info("iter $i, theta = $theta, norm(dtheta) = $(norm(dtheta))")
|
||||
norm(dtheta) < tolerance && return theta
|
||||
isnan(theta[1]) && break
|
||||
end
|
||||
theta[1] = NaN
|
||||
return theta
|
||||
# error("didn't converge in $i iterations")
|
||||
end
|
||||
|
||||
thetas = Vector{Float64}[]
|
||||
gen_rtic_grid!(thetas, element; npts=npts)
|
||||
map(calc_intersection!, thetas)
|
||||
filter!(t -> !isnan(t[1]), thetas)
|
||||
info_output && info("thetas: $thetas")
|
||||
secondary && filter!(t -> t[1] > 1.0e-12, thetas)
|
||||
# error("failure finding ray trace, thetas vec is empty")
|
||||
length(thetas) == 0 && return NaN, [NaN, NaN]
|
||||
sort!(thetas, alg=MergeSort, lt=(a,b)->a[1]<b[1])
|
||||
theta = thetas[1]
|
||||
return theta[1], theta[2:end]
|
||||
end
|
||||
|
||||
function calc_normal{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, time; deformed=true)
|
||||
x = element["geometry"](time)
|
||||
if deformed && haskey(element, "displacement")
|
||||
x += element["displacement"](time)
|
||||
end
|
||||
Q = [0.0 1.0; -1.0 0.0]
|
||||
dN = get_dbasis(element, xi, time)
|
||||
n = Q*(dN*x)
|
||||
n /= norm(n)
|
||||
return n
|
||||
end
|
||||
|
||||
function calc_normal{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, time; deformed=true)
|
||||
x = element["geometry"](time)
|
||||
if deformed && haskey(element, "displacement")
|
||||
x += element["displacement"](time)
|
||||
end
|
||||
dN = get_dbasis(element, xi, time)
|
||||
J = transpose(sum([kron(dN[:,i], x[i]') for i=1:length(x)]))
|
||||
n = cross(J[:,1], J[:,2])
|
||||
n /= norm(n)
|
||||
return n
|
||||
end
|
||||
|
||||
"""
|
||||
References
|
||||
----------
|
||||
[1] http://fp.optics.arizona.edu/optomech/Fall13/Notes/6%20Mirror%20matrices.pdf
|
||||
"""
|
||||
function calc_reflection(element::Element, xi, k, time; deformed=true)
|
||||
n = calc_normal(element, xi, time; deformed=deformed)
|
||||
k2 = k - 2*vecdot(k, n)*n
|
||||
return k2
|
||||
end
|
||||
|
||||
@@ -1,79 +0,0 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM.Testing
|
||||
|
||||
@testset "find intersection of Seg3 element" begin
|
||||
el = Element(Seg3, [1, 2, 3])
|
||||
update!(el, "geometry", Vector{Float64}[[0.0, 1.0], [1.0, 0.0], sqrt(2.0)/2.0*[1.0, 1.0]])
|
||||
s = [ 1.0, 0.5]
|
||||
d = [-1.0, 0.0]
|
||||
t, xi = find_intersection(el, s, d, 0.0)
|
||||
x = el("geometry", xi, 0.0)
|
||||
@test isapprox(x, [0.8604093371313943, 0.5])
|
||||
end
|
||||
|
||||
@testset "find intersection of 2. order NSeg element" begin
|
||||
# this is a exact quarter of circle
|
||||
a = 1.0
|
||||
b = 1.0
|
||||
el = Element(NSeg, [1, 2, 3])
|
||||
el.properties.order = 2
|
||||
el.properties.knots = [0.0, 0.0, 0.0, 1.0, 1.0, 1.0]
|
||||
el.properties.weights = [1.0, 1.0, 2.0]
|
||||
update!(el, "geometry", Vector{Float64}[[a, 0], [a, b], [0, b]])
|
||||
s = [ 1.0, 0.5]
|
||||
d = [-1.0, 0.0]
|
||||
t, xi = find_intersection(el, s, d, 0.0)
|
||||
X = el("geometry", xi, 0.0)
|
||||
@test isapprox(X, [sqrt(3)/2, 1/2])
|
||||
end
|
||||
|
||||
@testset "find intersection of 1. order NSurf" begin
|
||||
# node ordering, it's not same as in Quad4
|
||||
el = Element(NSurf, [1, 2, 3, 4])
|
||||
el.properties.order_u = 1
|
||||
el.properties.order_v = 1
|
||||
el.properties.knots_u = [0.0, 0.0, 1.0, 1.0]
|
||||
el.properties.knots_v = [0.0, 0.0, 1.0, 1.0]
|
||||
el.properties.weights = ones(2, 2)
|
||||
# +-- v
|
||||
# |
|
||||
# u
|
||||
nodes = Vector{Float64}[[0.0,0.0,0.0], [1.0,0.0,0.0],
|
||||
[0.0,1.0,0.0], [1.0,1.0,0.0]]
|
||||
update!(el, "geometry", nodes)
|
||||
s = [0.5, 0.5, 0.5]
|
||||
d = [0.0, 0.0, -1.0]
|
||||
t, xi = find_intersection(el, s, d, 0.0)
|
||||
X = el("geometry", xi, 0.0)
|
||||
@test isapprox(X, [0.5, 0.5, 0.0])
|
||||
k = calc_reflection(el, xi, d, 0.0)
|
||||
@test isapprox(k, [0.0, 0.0, 1.0])
|
||||
s = [1.5, 1.5, 0.5]
|
||||
t, xi = find_intersection(el, s, d, 0.0)
|
||||
@test isnan(t)
|
||||
end
|
||||
|
||||
@testset "find intersection of 3. order NSeg element with multiple reflections" begin
|
||||
el = Element(NSeg, [1, 2, 3, 4])
|
||||
el.properties.order = 3
|
||||
el.properties.knots = [0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0]
|
||||
el.properties.weights = [1.0, 1/3, 1/3, 1.0]
|
||||
update!(el, "geometry", Vector{Float64}[[1, 0], [1, 2], [-1, 2], [-1, 0]])
|
||||
s1 = [-sqrt(3)/2.0, 0.0]
|
||||
k1 = [ 0.0, 1.0]
|
||||
t1, xi1 = find_intersection(el, s1, k1, 0.0)
|
||||
s2 = el("geometry", xi1, 0.0)
|
||||
@test isapprox(s2, [-sqrt(3)/2, 0.5])
|
||||
k2 = calc_reflection(el, xi1, k1, 0.0)
|
||||
t2, xi2 = find_intersection(el, s2, k2, 0.0; secondary=true)
|
||||
s3 = el("geometry", xi2, 0.0)
|
||||
@test isapprox(s3, [0.0, 1.0])
|
||||
k3 = calc_reflection(el, xi2, k2, 0.0)
|
||||
t3, xi3 = find_intersection(el, s3, k3, 0.0; secondary=true)
|
||||
s4 = el("geometry", xi3, 0.0)
|
||||
@test isapprox(s4, [sqrt(3)/2, 0.5])
|
||||
end
|
||||
|
||||
Reference in New Issue
Block a user