nurbs elements

src/JuliaFEM.jl

@@ -26,6 +26,8 @@ include("elements.jl") # common element routines
 export Node, AbstractElement, Element, update!, get_connectivity
 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
 
 #include("hierarchical.jl") # P-elements
 #include("mortar_elements.jl") # Mortar elements

src/elements.jl

@@ -156,7 +156,7 @@ of integration scheme mainly for mass matrix.
 function get_integration_points(element::Element, change_order::Int)
     order = get_integration_order(element.properties)
     order += change_order
-    ips = get_integration_points(element.properties, Val{order})
+    ips = get_integration_points(element.properties, order)
     return [IP(i, w, xi) for (i, (w, xi)) in enumerate(ips)]
 end
 

src/integrate.jl

@@ -19,15 +19,15 @@ end
 
 function get_integration_points(::Type{Val{4}})
     weights = 1.0/36.0*[
-        18.0+sqrt(30.0),
-        18.0+sqrt(30.0),
         18.0-sqrt(30.0),
+        18.0+sqrt(30.0),
+        18.0+sqrt(30.0),
         18.0-sqrt(30.0)]
     points = [
-        sqrt( 3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
-        sqrt(-3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
-        sqrt( 3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)),
-        sqrt(-3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0))]
+        -sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0)),
+        -sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
+         sqrt(3.0/7.0 - 2.0/7.0*sqrt(6.0/5.0)),
+         sqrt(3.0/7.0 + 2.0/7.0*sqrt(6.0/5.0))]
     return weights, points
 end
 
@@ -47,53 +47,37 @@ function get_integration_points(::Type{Val{5}})
     return weights, points
 end
 
+function get_integration_points(order::Int64)
+    if order <= 5
+        return get_integration_points(Val{order})
+    else
+        points, weights = Base.QuadGK.gauss(Float64, order)
+        return weights, points
+    end
+end
+
 ### "cartesian" elements, integration rules comes from tensor product
 
 ### 1d elements
 
-typealias CartesianLineElement Union{Seg2, Seg3}
-typealias CartesianSurfaceElement Union{Quad4}
-typealias CartesianVolumeElement Union{Hex8}
+typealias CartesianLineElement Union{Seg2, Seg3, NSeg}
+typealias CartesianSurfaceElement Union{Quad4, NSurf}
+typealias CartesianVolumeElement Union{Hex8, NSolid}
 
-function get_integration_points(element::CartesianLineElement, ::Type{Val{1}})
-    w, xi = get_integration_points(Val{1})
-    [ (w[i], [xi[i]]) for i=1:1 ]
-end
-function get_integration_points(element::CartesianLineElement, ::Type{Val{2}})
-    w, xi = get_integration_points(Val{2})
-    [ (w[i], [xi[i]]) for i=1:2 ]
-end
-function get_integration_points(element::CartesianLineElement, ::Type{Val{3}})
-    w, xi = get_integration_points(Val{3})
-    [ (w[i], [xi[i]]) for i=1:3 ]
+function get_integration_points(element::CartesianLineElement, order::Int64)
+    w, xi = get_integration_points(order)
+    [ (w[i], [xi[i]]) for i=1:order ]
 end
 
-function get_integration_points(element::CartesianSurfaceElement, ::Type{Val{1}})
-    w, xi = get_integration_points(Val{1})
-    [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:1, j=1:1 ]
-end
-function get_integration_points(element::CartesianSurfaceElement, ::Type{Val{2}})
-    w, xi = get_integration_points(Val{2})
-    [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:2, j=1:2 ]
-end
-function get_integration_points(element::CartesianSurfaceElement, ::Type{Val{3}})
-    w, xi = get_integration_points(Val{3})
-    [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:3, j=1:3 ]
+function get_integration_points(element::CartesianSurfaceElement, order::Int64)
+    w, xi = get_integration_points(order)
+    [ (w[i]*w[j], [xi[i], xi[j]]) for i=1:order, j=1:order ]
 end
 
-function get_integration_points(element::CartesianVolumeElement, ::Type{Val{1}})
-    w, xi = get_integration_points(Val{1})
-    [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:1, j=1:1, k=1:1 ]
+function get_integration_points(element::CartesianVolumeElement, order::Int64)
+    w, xi = get_integration_points(order)
+    [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:order, j=1:order, k=1:order ]
 end
-function get_integration_points(element::CartesianVolumeElement, ::Type{Val{2}})
-    w, xi = get_integration_points(Val{2})
-    [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:2, j=1:2, k=1:2 ]
-end
-function get_integration_points(element::CartesianVolumeElement, ::Type{Val{3}})
-    w, xi = get_integration_points(Val{3})
-    [ (w[i]*w[j]*w[k], [xi[i], xi[j], xi[k]]) for i=1:3, j=1:3, k=1:3 ]
-end
-
 
 ### triangular and tetrahedral elements
 
@@ -188,6 +172,10 @@ function get_integration_points(element::TetrahedralElement, ::Type{Val{3}})
     return zip(weights, points)
 end
 
+function get_integration_points(element::Union{TriangularElement, TetrahedralElement}, order::Int64)
+    return get_integration_points(element, Val{order})
+end
+
 ### default number of integration points for each element
 ### 2 for linear elements, 3 for quadratic
 
@@ -205,11 +193,11 @@ end
 
 function get_integration_points(element::LinearElement)
     order = get_integration_order(element)
-    get_integration_points(element, Val{order})
+    get_integration_points(element, order)
 end
 
 function get_integration_points(element::QuadraticElement)
     order= get_integration_order(element)
-    get_integration_points(element, Val{order})
+    get_integration_points(element, order)
 end
 

src/nurbs.jl

@@ -0,0 +1,127 @@
+""" 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
+
+function NURBS(i, p, u, t)
+    p == 0 && return t[i] <= u <= t[i+1] ? 1.0 : 0.0
+    anom = u-t[i]
+    adenom = t[i+p]-t[i]
+    a = isapprox(adenom, 0.0) ? 0.0 : anom/adenom
+    bnom = t[i+p+1]-u
+    bdenom = t[i+p+1]-t[i+1]
+    b = isapprox(bdenom, 0.0) ? 0.0 : bnom/bdenom
+    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 = [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 = [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
+    w = element.properties.weights
+    nu = length(tu)
+    nv = length(tv)
+    nw = length(tw)
+    u, v, w = xi
+    N = [w[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
+
+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
+