Add Pyr5 elements, including med support and tests

src/JuliaFEM.jl

@@ -58,6 +58,7 @@ export Poi1,
        Tri3, Tri6, Tri7,
        Quad4, Quad8, Quad9,
        Tet4, Tet10,
+       Pyr5,
        Wedge6,
        Hex8, Hex20, Hex27
 

src/elements_lagrange.jl

@@ -13,6 +13,7 @@ global const ELEMENT_DESCRIPTIONS = Dict(
     "Quad9" => "9 node quadratic quadrangle element",
     "Tet4" => "4 node linear tetrahedral element",
     "Tet10" => "10 node quadratic tetrahedral element",
+    "Pyr5" => "5 node linear pyramid element",
     "Wedge6" => "6 node linear prismatic element (wedge)",
     "Wedge15" => "15 node quadratic prismatic element (wedge)",
     "Hex8" => "8 node linear hexahedral element",
@@ -31,6 +32,7 @@ global const ELEMENT_SIZES = Dict(
     "Quad9" => (2, 9),
     "Tet4" => (3, 4),
     "Tet10" => (3, 10),
+    "Pyr5" => (3,5),
     "Wedge6" => (3, 6),
     "Wedge15" => (3, 15),
     "Hex8" => (3, 8),
@@ -355,6 +357,39 @@ end
 
 #
 
+type Pyr5 <: AbstractElement
+end
+
+
+function get_reference_coordinates(::Type{Pyr5})
+    Vector{Float64}[
+        [-1.0,-1.0,-1.0], # N1
+        [ 1.0,-1.0,-1.0], # N2
+        [ 1.0, 1.0,-1.0], # N3
+        [-1.0, 1.0,-1.0], # N4
+        [ 0.0, 0.0, 1.0]] # N5
+end
+
+function get_interpolation_polynomial(::Type{Pyr5}, xi)
+    [
+		1.0/8.0*(1.0-1.0*xi[1])*(1.0-1.0*xi[2])*(1.0-1.0*xi[3])
+		1.0/8.0*(1.0+1.0*xi[1])*(1.0-1.0*xi[2])*(1.0-1.0*xi[3])
+		1.0/8.0*(1.0+1.0*xi[1])*(1.0+1.0*xi[2])*(1.0-1.0*xi[3])
+		1.0/8.0*(1.0-1.0*xi[1])*(1.0+1.0*xi[2])*(1.0-1.0*xi[3])
+		1.0/2.0*(1.0+xi[3])
+    ]'
+end
+
+function get_interpolation_polynomial(::Type{Pyr5}, xi, ::Type{Val{:partial_derivatives}})
+    [
+		-0.125*(1.0-xi[2])*(1.0-xi[3])   0.125*(1.0-xi[2])*(1.0-xi[3])   0.125*(1.0+xi[2])*(1.0-xi[3])  -0.125*(1.0+xi[2])*(1.0-xi[3])  0.0
+		-0.125*(1.0-xi[1])*(1.0-xi[3])  -0.125*(1.0+xi[1])*(1.0-xi[3])   0.125*(1.0+xi[1])*(1.0-xi[3])   0.125*(1.0-xi[1])*(1.0-xi[3])  0.0
+		-0.125*(1.0-xi[1])*(1.0-xi[2])  -0.125*(1.0+xi[1])*(1.0-xi[2])  -0.125*(1.0+xi[1])*(1.0+xi[2])  -0.125*(1.0-xi[1])*(1.0+xi[2])  0.5
+    ]
+end
+
+#
+
 type Wedge6 <: AbstractElement
 end
 
@@ -576,6 +611,7 @@ end
 @create_basis Quad9
 @create_basis Tet4
 @create_basis Tet10
+@create_basis Pyr5
 @create_basis Wedge6
 @create_basis Wedge15
 @create_basis Hex8
@@ -583,7 +619,8 @@ end
 @create_basis Hex27
 
 function inside(::Union{Type{Seg2}, Type{Seg3}, Type{Quad4}, Type{Quad8},
-                        Type{Quad9}, Type{Hex8}, Type{Hex20}, Type{Hex27}}, xi)
+                        Type{Quad9}, Type{Pyr5}, Type{Hex8}, Type{Hex20},
+                        Type{Hex27}}, xi)
     return all(-1.0 .<= xi .<= 1.0)
 end
 

src/integrate.jl

@@ -271,6 +271,27 @@ function get_integration_points(element::TetrahedralElement, ::Type{Val{4}})
     return zip(weights, points)
 end
 
+# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch12.d/AFEM.Ch12.pdf
+
+typealias PyramidalElement Union{Pyr5,}
+
+function get_integration_points(element::PyramidalElement, ::Type{Val{2}})
+    g1 =  0.5842373946721771876874344
+    g2 = -2.0/3.0
+    g3 =  2.0/5.0
+    w1 =  81.0/100.0
+    w2 =  125.0/27.0
+    weights = [w1, w1, w1, w1, w2]
+    points = Vector{Float64}[
+        [-g1, -g1, g2],
+        [ g1, -g1, g2],
+        [ g1,  g1, g2],
+        [-g1,  g1, g2],
+        [0.0, 0.0, g3],
+        ]
+    return zip(weights, points)
+end
+
 typealias PrismaticElement Union{Wedge6, Wedge15}
 
 function get_integration_points(element::PrismaticElement, ::Type{Val{2}})
@@ -408,14 +429,14 @@ function get_integration_points(element::PrismaticElement, ::Type{Val{3}})
 end
 
 function get_integration_points(element::Union{TriangularElement,
-            TetrahedralElement, PrismaticElement}, order::Int64)
+            TetrahedralElement, PyramidalElement, PrismaticElement}, 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
 
-typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Wedge6, Hex8}
+typealias LinearElement Union{Seg2, Tri3, Quad4, Tet4, Pyr5, Wedge6, Hex8}
 
 typealias QuadraticElement Union{Seg3, Tri6, Tri7, Tet10, Quad8, Quad9, Wedge15, Hex20, Hex27}
 

src/preprocess_aster_reader.jl

@@ -176,6 +176,7 @@ end
 global const med_connectivity = Dict{Symbol, Vector{Int}}(
     :Tet4   => [4,3,1,2],
     :Tet10  => [4,3,1,2,10,7,8,9,6,5],
+    :Pyr5   => [1,4,3,2,5],
     :Wedge6 => [4,5,6,1,2,3],
     :Hex8   => [4,8,7,3,1,5,6,2],
     :Hex20  => [4,8,7,3,1,5,6,2,20,15,19,11,12,16,14,10,17,13,18,9],
@@ -209,8 +210,8 @@ global const mapping = Dict(
     :H20 => :Hex20,
     :H27 => :Hex27,
 
-    :PY5 => :Pyramid5,
-    :P13 => :Pyramid13,
+    :PY5 => :Pyr5,
+    :P13 => :Pyr13,
 
     )
 

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, Wedge15, Hex8, Tet10, Hex20, Hex27}
+typealias Elasticity3DVolumeElements Union{Tet4, Pyr5, Wedge6, Wedge15, Hex8, Tet10, Hex20, Hex27}
 
 function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}})
     param_keys = keys(params)

src/problems_heat.jl

@@ -67,7 +67,7 @@ function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Eleme
     info("Unknown element type $E for 3d heat problem!")
 end
 
-typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8, Hex20, Hex27}
+typealias Heat3DVolumeElements Union{Tet4, Tet10, Pyr5, Hex8, Hex20, Hex27}
 typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4, Quad8, Quad9}
 
 function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})