From be5a1bc8dee8ac37dbdbfb2358bb8902ed97741e Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 30 Nov 2016 16:18:57 +0200 Subject: [PATCH] added wedge15 element + integration rules --- src/elements_lagrange.jl | 52 ++++++++++++++++++++++++ src/integrate.jl | 81 +++++++++++++++++++++++++++++++++----- src/problems_elasticity.jl | 2 +- 3 files changed, 124 insertions(+), 11 deletions(-) diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl index f3707ce..a9266ba 100644 --- a/src/elements_lagrange.jl +++ b/src/elements_lagrange.jl @@ -473,6 +473,56 @@ end # +type Wedge15 <: AbstractElement +end + +function description(::Type{Wedge15}) + "15 node prismatic element (wedge)" +end + +function size(element::Element{Wedge15}) + return (3, 15) +end + +function length(element::Element{Wedge15}) + return 15 +end + +function get_reference_coordinates(::Type{Wedge15}) + Vector{Float64}[ + [0.0, 0.0, -1.0], # N1 + [1.0, 0.0, -1.0], # N2 + [0.0, 1.0, -1.0], # N3 + [0.0, 0.0, 1.0], # N4 + [1.0, 0.0, 1.0], # N5 + [0.0, 1.0, 1.0], # N6 + [0.5, 0.0, -1.0], # N7 + [0.5, 0.5, -1.0], # N8 + [0.0, 0.5, -1.0], # N9 + [0.5, 0.0, 1.0], # N10 + [0.5, 0.5, 1.0], # N11 + [0.0, 0.5, 1.0], # N12 + [0.0, 0.0, 0.0], # N13 + [1.0, 0.0, 0.0], # N14 + [0.0, 1.0, 0.0]] # N15 +end + +function get_interpolation_polynomial(::Type{Wedge15}, x) + [ + 1 x[1] x[1]^2 x[2] x[1]*x[2] x[2]^2 x[3] x[1]*x[3] x[1]^2*x[3] x[2]*x[3] x[1]*x[2]*x[3] x[2]^2*x[3] x[3]^2 x[1]*x[3]^2 x[2]*x[3]^2 + ] +end + +function get_interpolation_polynomial(::Type{Wedge15}, x, ::Type{Val{:partial_derivatives}}) + [ + 0 1 2*x[1] 0 x[2] 0 0 x[3] 2*x[1]*x[3] 0 x[2]*x[3] 0 0 x[3]^2 0 + 0 0 0 1 x[1] 2*x[2] 0 0 0 x[3] x[1]*x[3] 2*x[2]*x[3] 0 0 x[3]^2 + 0 0 0 0 0 0 1 x[1] x[1]^2 x[2] x[1]*x[2] x[2]^2 2*x[3] 2*x[1]*x[3] 2*x[2]*x[3] + ] +end + +# + type Hex8 <: AbstractElement end @@ -664,6 +714,7 @@ end @create_basis Tet4 @create_basis Tet10 @create_basis Wedge6 +@create_basis Wedge15 @create_basis Hex8 @create_basis Hex20 @create_basis Hex27 @@ -680,3 +731,4 @@ end function get_reference_coordinates{E}(element::Element{E}) get_reference_coordinates(E) end + diff --git a/src/integrate.jl b/src/integrate.jl index 921b2f2..2ca62e9 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -103,17 +103,30 @@ function get_integration_points(element::TriangularElement, ::Type{Val{2}}) return zip(weights, points) end +#= +""" Note, this rule is having negative weight. """ function get_integration_points(element::TriangularElement, ::Type{Val{3}}) - weights = 0.5*[ - -0.5625, - 0.5208333333333333, - 0.5208333333333333, - 0.5208333333333333] + weights = [-27.0/96.0, 25.0/96.0, 25.0/96.0, 25.0/96.0] points = Vector{Float64}[ [1.0/3.0, 1.0/3.0], - [0.2, 0.2], - [0.2, 0.6], - [0.6, 0.2]] + [1.0/5.0, 1.0/5.0], + [1.0/5.0, 3.0/5.0], + [3.0/5.0, 1.0/5.0]] + return zip(weights, points) +end +=# + +function get_integration_points(element::TriangularElement, ::Type{Val{3}}) + weights = [ + 1.5902069087198858469718450103758e-01, + 9.0979309128011415302815498962418e-02, + 1.5902069087198858469718450103758e-01, + 9.0979309128011415302815498962418e-02] + points = Vector{Float64}[ + [1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01], + [6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02], + [1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01], + [6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01]] return zip(weights, points) end @@ -225,8 +238,9 @@ function get_integration_points(element::TetrahedralElement, ::Type{Val{4}}) return zip(weights, points) end -typealias PrismaticElement Union{Wedge6} +typealias PrismaticElement Union{Wedge6, Wedge15} +#= function get_integration_points(element::PrismaticElement, ::Type{Val{2}}) weights = 1/6*[1.0, 1.0, 1.0, 1.0, 1.0, 1.0] points = Vector{Float64}[ @@ -238,6 +252,53 @@ function get_integration_points(element::PrismaticElement, ::Type{Val{2}}) [0.5, 0.5, 1.0/sqrt(3)]] return zip(weights, points) end +=# + +# tensor product of triangular element + segment element +function get_integration_points(element::PrismaticElement, ::Type{Val{2}}) + weights = 1.0/6.0*[1.0, 1.0, 1.0, 1.0, 1.0, 1.0] + points = Vector{Float64}[ + [2.0/3.0, 1.0/6.0, -1.0/sqrt(3.0)], + [1.0/6.0, 2.0/3.0, -1.0/sqrt(3.0)], + [1.0/6.0, 1.0/6.0, -1.0/sqrt(3.0)], + [2.0/3.0, 1.0/6.0, +1.0/sqrt(3.0)], + [1.0/6.0, 2.0/3.0, +1.0/sqrt(3.0)], + [1.0/6.0, 1.0/6.0, +1.0/sqrt(3.0)], + ] + return zip(weights, points) +end + +# tensor product of triangular element + segment element +function get_integration_points(element::PrismaticElement, ::Type{Val{3}}) + weights = [ + 5.0/9.0*1.5902069087198858469718450103758e-01, + 5.0/9.0*9.0979309128011415302815498962418e-02, + 5.0/9.0*1.5902069087198858469718450103758e-01, + 5.0/9.0*9.0979309128011415302815498962418e-02, + 8.0/9.0*1.5902069087198858469718450103758e-01, + 8.0/9.0*9.0979309128011415302815498962418e-02, + 8.0/9.0*1.5902069087198858469718450103758e-01, + 8.0/9.0*9.0979309128011415302815498962418e-02, + 5.0/9.0*1.5902069087198858469718450103758e-01, + 5.0/9.0*9.0979309128011415302815498962418e-02, + 5.0/9.0*1.5902069087198858469718450103758e-01, + 5.0/9.0*9.0979309128011415302815498962418e-02] + points = Vector{Float64}[ + [1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01, -sqrt(3.0/5.0)], + [6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02, -sqrt(3.0/5.0)], + [1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01, -sqrt(3.0/5.0)], + [6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01, -sqrt(3.0/5.0)], + [1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01, 0.0], + [6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02, 0.0], + [1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01, 0.0], + [6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01, 0.0], + [1.5505102572168219018027159252941e-01, 1.7855872826361642311703513337422e-01, sqrt(3.0/5.0)], + [6.4494897427831780981972840747059e-01, 7.5031110222608118177475598324603e-02, sqrt(3.0/5.0)], + [1.5505102572168219018027159252941e-01, 6.6639024601470138670269327409637e-01, sqrt(3.0/5.0)], + [6.4494897427831780981972840747059e-01, 2.8001991549907407200279599420481e-01, sqrt(3.0/5.0)], + ] + return zip(weights, points) +end function get_integration_points(element::Union{TriangularElement, TetrahedralElement, PrismaticElement}, order::Int64) @@ -249,7 +310,7 @@ end typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Wedge6, Hex8} -typealias QuadraticElement Union{Seg3, Tri6, Tri7, Tet10, Quad8, Quad9, Hex20, Hex27} +typealias QuadraticElement Union{Seg3, Tri6, Tri7, Tet10, Quad8, Quad9, Wedge15, Hex20, Hex27} function get_integration_order(element::LinearElement) return 2 diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index de407bb..e3cc30d 100644 --- a/src/problems_elasticity.jl +++ b/src/problems_elasticity.jl @@ -69,7 +69,7 @@ end typealias Elasticity2DSurfaceElements Union{Poi1, Seg2, Seg3} typealias Elasticity2DVolumeElements Union{Tri3, Tri6, Quad4, Quad8, Quad9} typealias Elasticity3DSurfaceElements Union{Poi1, Tri3, Tri6, Quad4, Quad8, Quad9} -typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Hex8, Tet10, Hex20, Hex27} +typealias Elasticity3DVolumeElements Union{Tet4, Wedge6, Wedge15, Hex8, Tet10, Hex20, Hex27} function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}}) param_keys = keys(params)