optics

src/JuliaFEM.jl

@@ -23,11 +23,11 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 
 ### ELEMENTS ###
 include("elements.jl") # common element routines
-export Node, AbstractElement, Element, update!, get_connectivity
+export Node, AbstractElement, Element, update!, get_connectivity, get_basis, get_dbasis
 include("lagrange_macro.jl") # Continuous Galerkin (Lagrange) elements generated using macro
 export Seg2, Seg3, Tri3, Tri6, Quad4, Hex8, Tet4, Tet10
 include("nurbs.jl")
-export NSeg, NSurf, NSolid
+export NSeg, NSurf, NSolid, is_nurbs
 
 #include("hierarchical.jl") # P-elements
 #include("mortar_elements.jl") # Mortar elements
@@ -64,6 +64,9 @@ include("solver_utils.jl")
 include("solvers.jl")
 export Solver
 
+include("optics.jl")
+export find_intersection, calc_reflection
+
 ### MORTAR STUFF ###
 include("mortar.jl")  # mortar projection
 

src/common.jl

@@ -46,6 +46,7 @@ function set_debug_off!()
     pop!(ENV, "DEBUG");
 end
 
+#=
 """ Simple linspace extension to arrays.
 
 Examples
@@ -60,6 +61,7 @@ Examples
 function linspace{T<:Array}(X1::T, X2::T, n)
     [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
 end
+=#
 
 function resize!(A::SparseMatrixCSC, m::Int64, n::Int64)
     (n == A.n) && (m == A.m) && return

src/elements.jl

@@ -116,7 +116,11 @@ function update!(element::Element, field_name::ASCIIString, data::Dict)
 end
 
 function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair})
-    element[field_name] = data
+    if haskey(element, field_name)
+        update!(element[field_name], data)
+    else
+        element[field_name] = data
+    end
 end
 
 function update!(elements::Vector, field_name::ASCIIString, data)

src/nurbs.jl

@@ -1,3 +1,6 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
 """ NURBS segment. """
 type NSeg <: AbstractElement
     order :: Int
@@ -62,7 +65,7 @@ function get_basis(element::Element{NSeg}, xi::Vector, time)
     w = element.properties.weights
     nu = length(tu)-pu-1
     u = xi[1]
-    N = [w[j]*NURBS(j,pu,u,tu) for j=1:nu]
+    N = vec([w[j]*NURBS(j,pu,u,tu) for j=1:nu])'
     return N/sum(N)
 end
 
@@ -75,7 +78,7 @@ function get_basis(element::Element{NSurf}, xi::Vector, time)
     nu = length(tu)-pu-1
     nv = length(tv)-pv-1
     u, v = xi
-    N = [w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv) for i=1:nu, j=1:nv]
+    N = vec([w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv) for i=1:nu, j=1:nv])'
     return N / sum(N)
 end
 
@@ -125,3 +128,11 @@ function size(element::Element{NSolid})
     return (3, length(element))
 end
 
+function is_nurbs(element::Element)
+    return false
+end
+
+function is_nurbs{E<:Union{NSeg, NSurf, NSolid}}(element::Element{E})
+    return true
+end
+

src/optics.jl

@@ -0,0 +1,124 @@
+# 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
+
+"""
+References
+----------
+[1] http://fp.optics.arizona.edu/optomech/Fall13/Notes/6%20Mirror%20matrices.pdf
+"""
+function calc_reflection{S<:Union{Seg2, Seg3, NSeg}}(element::Element{S}, xi, k, 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)
+    k2 = k - 2*vecdot(k, n)*n
+    return k2
+end
+
+function calc_reflection{S<:Union{Quad4, NSurf}}(element::Element{S}, xi, k, 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)
+    k2 = k - 2*vecdot(k, n)*n
+    return k2
+end
+

test/test_elements.jl

@@ -1,61 +1,10 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-module ElementTests
-
+using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM.Core: AbstractElement, Element, Field, FieldSet, test_element
-using JuliaFEM.Core: Tri3, Quad4
-import JuliaFEM.Core: get_basis, get_dbasis, calculate_normal_tangential_coordinates!
-import Base: size
-
-""" Prototype element
-
-This should always pass test_element if everything is ok.
-"""
-abstract TestElement <: AbstractElement
-
-function get_basis(::Type{TestElement}, xi::Vector{Float64})
-    1/4*[
-        (1-xi[1])*(1-xi[2])
-        (1+xi[1])*(1-xi[2])
-        (1+xi[1])*(1+xi[2])
-        (1-xi[1])*(1+xi[2])]'
-end
-
-function get_dbasis(::Type{TestElement}, xi::Vector{Float64})
-    1/4*[
-        -(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
-        -(1-xi[1])   -(1+xi[1])   (1+xi[1])   (1-xi[1])]
-end
-
-function size(::Type{TestElement})
-    return (2, 4)
-end
-
-""" Return test element with some fields. """
-function get_element()
-    el = Element{TestElement}([1, 2, 3, 4])
-    el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
-    el["temperature"] = (
-        0.0 => [0.0, 0.0, 0.0, 0.0],
-        1.0 => [1.0, 2.0, 3.0, 4.0])
-    el["displacement"] = (
-        0.0 => Vector{Float64}[[0.0,0.0], [0.0, 0.0], [0.0,0.0], [0.0,0.0]],
-        1.0 => Vector{Float64}[[0.0,0.0], [1.0,-1.0], [2.0,3.0], [0.0,0.0]])
-    return el
-end
-
-function test_mock_element()
-    test_element(TestElement)
-end
-
-function test_add_fields_to_element()
-    el = get_element()
-    info(el.fields)
-end
-
+#=
 function test_interpolate()
     el = get_element()
     @test isapprox(el("geometry", [0.0, 0.0]), [0.5, 0.5])
@@ -89,7 +38,6 @@ function test_calculate_normal_tangential_coordinates()
     R = [n t1 t2]
     @test isapprox(el("normal-tangential coordinates", [0.0, 0.0], 0.0), R)
 end
-#test_calculate_normal_tangential_coordinates()
 
 function test_manifold_determinant()
     el = Quad4([1, 2, 3, 4])
@@ -100,6 +48,26 @@ function test_manifold_determinant()
     d_expected = 0.25
     @test d == d_expected
 end
-#test_manifold_determinant()
 
+@testset "add new discrete constant time-variant field and interpolate it" begin
+    element = Element(Quad4, [1, 2, 3, 4])
+    element["my field"] = (0.0 => 0.0, 1.0 => 1.0)
+    @test isapprox(element("my field", [0.0, 0.0], 0.5), 0.5)
+    update!(element, "my field 2", 0.0 => 0.0, 1.0 => 1.0)
+    @test isapprox(element("my field 2", [0.0, 0.0], 0.5), 0.5)
 end
+
+=#
+
+@testset "test add time dependent field to element" begin
+    el = Element(Seg2, [1, 2])
+    u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]
+    u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]]
+    update!(el, "displacement", 0.0 => u1)
+    update!(el, "displacement", 1.0 => u2)
+    @test length(el["displacement"]) == 2
+    @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0])
+    @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5])
+    @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0])
+end
+

test/test_elements_add_fields.jl

@@ -1,14 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using JuliaFEM
 using JuliaFEM.Test
 
-using JuliaFEM.Core: Quad4, update!
 
-@testset "add new discrete constant time-variant field and interpolate it" begin
-    element = Quad4([1, 2, 3, 4])
-    element["my field"] = (0.0 => 0.0, 1.0 => 1.0)
-    @test isapprox(element("my field", [0.0, 0.0], 0.5), 0.5)
-    update!(element, "my field 2", 0.0 => 0.0, 1.0 => 1.0)
-    @test isapprox(element("my field 2", [0.0, 0.0], 0.5), 0.5)
-end

test/test_fields.jl

@@ -24,3 +24,4 @@ end
     update!(f, 2.0)
     @test f.data == 2.0
 end
+

test/test_find_edge_intersections.jl

@@ -0,0 +1,79 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM
+using JuliaFEM.Test
+
+@testset "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
+