diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 8794663..9d4828c 100644 --- a/src/JuliaFEM.jl +++ b/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 diff --git a/src/common.jl b/src/common.jl index d78deb0..25e15fd 100644 --- a/src/common.jl +++ b/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 diff --git a/src/elements.jl b/src/elements.jl index c63c9a8..e433def 100644 --- a/src/elements.jl +++ b/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) diff --git a/src/nurbs.jl b/src/nurbs.jl index f6f2bb0..ea995b9 100644 --- a/src/nurbs.jl +++ b/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 + diff --git a/src/optics.jl b/src/optics.jl new file mode 100644 index 0000000..b15bafe --- /dev/null +++ b/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] [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 + diff --git a/test/test_elements_add_fields.jl b/test/test_elements_add_fields.jl index e67c50e..3e527fa 100644 --- a/test/test_elements_add_fields.jl +++ b/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 diff --git a/test/test_fields.jl b/test/test_fields.jl index ea83b14..9ab554c 100644 --- a/test/test_fields.jl +++ b/test/test_fields.jl @@ -24,3 +24,4 @@ end update!(f, 2.0) @test f.data == 2.0 end + diff --git a/test/test_find_edge_intersections.jl b/test/test_find_edge_intersections.jl new file mode 100644 index 0000000..5317ac2 --- /dev/null +++ b/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 +