fix nurbs basis

src/FEMBasis.jl

@@ -8,6 +8,7 @@ abstract type AbstractBasis end
 
 include("create_basis.jl")
 export get_reference_element_coordinates, eval_basis!, eval_dbasis!
+
 include("lagrange_segments.jl")
 export Seg2, Seg3
 include("lagrange_quadrangles.jl")
@@ -22,4 +23,13 @@ include("lagrange_wedges.jl")
 export Wedge6, Wedge15
 include("lagrange_pyramids.jl")
 export Pyr5
+
+include("nurbs.jl")
+include("nurbs_segment.jl")
+export NSeg
+include("nurbs_surface.jl")
+export NSurf
+include("nurbs_solid.jl")
+export NSolid
+
 end

src/nurbs.jl

@@ -1,54 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
 
-using ForwardDiff
-# TODO: evaluate partial derivatives of basis functions without forwarddiff
-
-""" NURBS segment. """
-type NSeg <: AbstractElement
-    order :: Int
-    knots :: Vector{Float64}
-    weights :: Vector{Float64}
-end
-
-function NSeg()
-    NSeg(1,
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(4))
-end
-
-type NSurf <: AbstractElement
-    order_u :: Int
-    order_v :: Int
-    knots_u :: Vector{Float64}
-    knots_v :: Vector{Float64}
-    weights :: Matrix{Float64}
-end
-
-function NSurf()
-    NSurf(1, 1,
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(2, 2))
-end
-
-type NSolid <: AbstractElement
-    order_u :: Int
-    order_v :: Int
-    order_w :: Int
-    knots_u :: Vector{Float64}
-    knots_v :: Vector{Float64}
-    knots_w :: Vector{Float64}
-    weights :: Array{Float64, 3}
-end
-
-function NSolid()
-    NSolid(1, 1, 1,
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(2, 2, 2))
-end
+import Base: size, length
 
 function NURBS(i, p, u, t)
     p == 0 && return t[i] <= u <= t[i+1] ? 1.0 : 0.0
@@ -61,89 +14,3 @@ function NURBS(i, p, u, t)
     result = a*NURBS(i,p-1,u,t) + b*NURBS(i+1,p-1,u,t)
     return result
 end
-
-function get_basis(element::Element{NSeg}, xi::Vector, time)
-    pu = element.properties.order
-    tu = element.properties.knots
-    w = element.properties.weights
-    nu = length(tu)-pu-1
-    u = xi[1]
-    N = vec([w[j]*NURBS(j,pu,u,tu) for j=1:nu])'
-    return N/sum(N)
-end
-
-function get_basis(element::Element{NSurf}, xi::Vector, time)
-    pu = element.properties.order_u
-    pv = element.properties.order_v
-    tu = element.properties.knots_u
-    tv = element.properties.knots_v
-    w = element.properties.weights
-    nu = length(tu)-pu-1
-    nv = length(tv)-pv-1
-    u, v = xi
-    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
-
-function get_basis(element::Element{NSolid}, xi::Vector, time)
-    pu = element.properties.order_u
-    pv = element.properties.order_v
-    pw = element.properties.order_w
-    tu = element.properties.knots_u
-    tv = element.properties.knots_v
-    tw = element.properties.knots_w
-    weights = element.properties.weights
-    nu = length(tu)-pu-1
-    nv = length(tv)-pv-1
-    nw = length(tw)-pw-1
-    u, v, w = xi
-    N = vec([weights[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw])'
-    return N / sum(N)
-end
-
-# TODO: evaluate partial derivatives of basis functions without forwarddiff
-""" Evaluate partial derivatives of basis functions using ForwardDiff. """
-function get_dbasis{E<:Union{NSeg, NSurf, NSolid}}(element::Element{E}, ip, time)
-    xi = isa(ip, IP) ? ip.coords : ip
-    basis(xi) = vec(get_basis(element, xi, time))
-    return ForwardDiff.jacobian(basis, xi)'
-end
-
-function length(element::Element{NSeg})
-    nu = length(element.properties.knots) - element.properties.order - 1
-    return nu
-end
-
-function size(element::Element{NSeg})
-    return (1, length(element))
-end
-
-function length(element::Element{NSurf})
-    nu = length(element.properties.knots_u) - element.properties.order_u - 1
-    nv = length(element.properties.knots_v) - element.properties.order_v - 1
-    return nu*nv
-end
-
-function size(element::Element{NSurf})
-    return (2, length(element))
-end
-
-function length(element::Element{NSolid})
-    nu = length(element.properties.knots_u) - element.properties.order_u - 1
-    nv = length(element.properties.knots_v) - element.properties.order_v - 1
-    nw = length(element.properties.knots_w) - element.properties.order_w - 1
-    return nu*nv*nw
-end
-
-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/nurbs_segment.jl

@@ -0,0 +1,37 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+""" NURBS segment. """
+type NSeg <: AbstractBasis
+    order :: Int
+    knots :: Vector{Float64}
+    weights :: Vector{Float64}
+end
+
+function NSeg()
+    NSeg(1,
+        [-1.0, -1.0, 1.0, 1.0],
+        ones(4))
+end
+
+function length(basis::NSeg)
+    nu = length(basis.knots) - basis.order - 1
+    return nu
+end
+
+function size(basis::NSeg)
+    return (1, length(basis))
+end
+
+function eval_basis!(basis::NSeg, N::Matrix, xi)
+    pu = basis.order
+    tu = basis.knots
+    w = basis.weights
+    nu = length(tu)-pu-1
+    u = xi[1]
+    for j=1:nu
+        N[j] = w[j]*NURBS(j,pu,u,tu)
+    end
+    N ./= sum(N)
+    return N
+end

src/nurbs_solid.jl

@@ -0,0 +1,59 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+type NSolid <: AbstractBasis
+    order_u :: Int
+    order_v :: Int
+    order_w :: Int
+    knots_u :: Vector{Float64}
+    knots_v :: Vector{Float64}
+    knots_w :: Vector{Float64}
+    weights :: Array{Float64, 3}
+end
+
+function NSolid()
+    NSolid(1, 1, 1,
+        [-1.0, -1.0, 1.0, 1.0],
+        [-1.0, -1.0, 1.0, 1.0],
+        [-1.0, -1.0, 1.0, 1.0],
+        ones(2, 2, 2))
+end
+
+function length(basis::NSolid)
+    nu = length(basis.knots_u) - basis.order_u - 1
+    nv = length(basis.knots_v) - basis.order_v - 1
+    nw = length(basis.knots_w) - basis.order_w - 1
+    return nu*nv*nw
+end
+
+function size(basis::NSolid)
+    return (3, length(basis))
+end
+
+function eval_basis!(basis::NSolid, N::Matrix, xi)
+    pu = basis.order_u
+    pv = basis.order_v
+    pw = basis.order_w
+    tu = basis.knots_u
+    tv = basis.knots_v
+    tw = basis.knots_w
+    weights = basis.weights
+    nu = length(tu)-pu-1
+    nv = length(tv)-pv-1
+    nw = length(tw)-pw-1
+    u, v, w = xi
+    n = 1
+    for i=1:nu
+        for j=1:nv
+            for k=1:nw
+                A = NURBS(i,pu,u,tu)
+                B = NURBS(j,pv,v,tv)
+                C = NURBS(k,pw,w,tw)
+                N[n] = weights[i,j,k]*A*B*C
+                n += 1
+            end
+        end
+    end
+    N ./= sum(N)
+    return N
+end

src/nurbs_surface.jl

@@ -0,0 +1,47 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+type NSurf <: AbstractBasis
+    order_u :: Int
+    order_v :: Int
+    knots_u :: Vector{Float64}
+    knots_v :: Vector{Float64}
+    weights :: Matrix{Float64}
+end
+
+function NSurf()
+    NSurf(1, 1,
+        [-1.0, -1.0, 1.0, 1.0],
+        [-1.0, -1.0, 1.0, 1.0],
+        ones(2, 2))
+end
+
+function length(basis::NSurf)
+    nu = length(basis.knots_u) - basis.order_u - 1
+    nv = length(basis.knots_v) - basis.order_v - 1
+    return nu*nv
+end
+
+function size(basis::NSurf)
+    return (2, length(basis))
+end
+
+function eval_basis!(basis::NSurf, N::Matrix, xi)
+    pu = basis.order_u
+    pv = basis.order_v
+    tu = basis.knots_u
+    tv = basis.knots_v
+    w = basis.weights
+    nu = length(tu)-pu-1
+    nv = length(tv)-pv-1
+    u, v = xi
+    n = 1
+    for i=1:nu
+        for j=1:nv
+            N[n] = w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)
+            n += 1
+        end
+    end
+    N ./= sum(N)
+    return N
+end

test/runtests.jl

@@ -8,6 +8,7 @@ const to = TimerOutput()
 test_files = String[]
 push!(test_files, "test_create_basis.jl")
 push!(test_files, "test_lagrange_elements.jl")
+push!(test_files, "test_nurbs_elements.jl")
 
 @testset "FEMBasis.jl" begin
     for fn in test_files

test/test_nurbs_elements.jl

@@ -0,0 +1,32 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBasis.jl/blob/master/LICENSE
+
+using Base.Test
+using FEMBasis
+
+@testset "NSeg interpolate" begin
+    basis = NSeg()
+    @test size(basis) == (1, 2)
+    @test length(basis) == 2
+    N = zeros(1, 2)
+    eval_basis!(basis, N, (0.0, ))
+    @test isapprox(N, [0.5 0.5])
+end
+
+@testset "NSurf interpolate" begin
+    basis = NSurf()
+    @test size(basis) == (2, 4)
+    @test length(basis) == 4
+    N = zeros(1, 4)
+    eval_basis!(basis, N, (0.0, 0.0))
+    @test isapprox(N, [0.25 0.25 0.25 0.25])
+end
+
+@testset "NSolid interpolate" begin
+    basis = NSolid()
+    @test size(basis) == (3, 8)
+    @test length(basis) == 8
+    N = zeros(1, 8)
+    eval_basis!(basis, N, (0.0, 0.0, 0.0))
+    @test isapprox(N, [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125])
+end