From 07812f3bc753b68e0a7c3c309fc7bb6992d73554 Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 1 Jun 2016 20:43:25 +0300 Subject: [PATCH] nurbs elements --- src/JuliaFEM.jl | 2 + src/elements.jl | 2 +- src/integrate.jl | 78 ++++++++++++----------------- src/nurbs.jl | 127 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 163 insertions(+), 46 deletions(-) create mode 100644 src/nurbs.jl diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index d8fa744..8794663 100644 --- a/src/JuliaFEM.jl +++ b/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 diff --git a/src/elements.jl b/src/elements.jl index b5dd2ff..c63c9a8 100644 --- a/src/elements.jl +++ b/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 diff --git a/src/integrate.jl b/src/integrate.jl index 56e3b60..916a586 100644 --- a/src/integrate.jl +++ b/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 diff --git a/src/nurbs.jl b/src/nurbs.jl new file mode 100644 index 0000000..f6f2bb0 --- /dev/null +++ b/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 +