Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-15 18:14:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

added wedge15 element + integration rules

Browse Source
This commit is contained in:
Jukka Aho
2016-11-30 16:18:57 +02:00
parent 8ed0f20d26
commit be5a1bc8de
3 changed files with 124 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elements_lagrange.jl
+52
View File
@@ -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
src/integrate.jl
+71 -10
View File
@@ -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
src/problems_elasticity.jl
+1 -1
View File
@@ -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)