diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 1c9df7f..4d908c7 100644 --- a/src/JuliaFEM.jl +++ b/src/JuliaFEM.jl @@ -58,6 +58,7 @@ export Poi1, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9, Tet4, Tet10, + Pyr5, Wedge6, Hex8, Hex20, Hex27 diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl index 83a3269..9dfad83 100644 --- a/src/elements_lagrange.jl +++ b/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 diff --git a/src/integrate.jl b/src/integrate.jl index 14219f0..86dd72c 100644 --- a/src/integrate.jl +++ b/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} diff --git a/src/preprocess_aster_reader.jl b/src/preprocess_aster_reader.jl index 07c324f..143a485 100644 --- a/src/preprocess_aster_reader.jl +++ b/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, ) diff --git a/src/problems_elasticity.jl b/src/problems_elasticity.jl index da28467..26eb708 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, 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) diff --git a/src/problems_heat.jl b/src/problems_heat.jl index c7956b0..0874020 100644 --- a/src/problems_heat.jl +++ b/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")}}) diff --git a/test/test_elasticity_med_pyr5_point_load.jl b/test/test_elasticity_med_pyr5_point_load.jl new file mode 100644 index 0000000..8cf3ecc --- /dev/null +++ b/test/test_elasticity_med_pyr5_point_load.jl @@ -0,0 +1,46 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Testing + +@testset "test Pyr5 elasticity with point load" begin + fn = Pkg.dir("JuliaFEM") * "/geometry/3d_pyr/Pyr5.med" + mesh = aster_read_mesh(fn) + element_sets = join(keys(mesh.element_sets), ", ") + info("element sets: $element_sets") + + element1 = create_elements(mesh,"Pyr5") + baseQuad = create_elements(mesh,"baseQuad") + tipPoint = Element(Poi1, collect(mesh.node_sets[:tipPoint])) + + update!([element1,baseQuad,tipPoint], "geometry", mesh.nodes) + update!([element1], "youngs modulus", 288.0) + update!([element1], "poissons ratio", 1/3) + + update!([tipPoint], "displacement traction force 1", 5.0) + update!([tipPoint], "displacement traction force 2", -7.0) + update!([tipPoint], "displacement traction force 3", 3.0) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + elasticity_problem.properties.finite_strain = false + push!(elasticity_problem, element1) + push!(elasticity_problem, tipPoint) + + baseQuad[1]["displacement 1"] = 0.0 + baseQuad[1]["displacement 2"] = 0.0 + baseQuad[1]["displacement 3"] = 0.0 + boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement") + push!(boundary_problem, baseQuad) + + solver = LinearSolver(elasticity_problem, boundary_problem) + solver() + + disp = element1[1]("displacement", [0.0, 0.0, 1.0], 0.0) + info("########################################################") + info("displacement at tip: $disp") + # Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt + u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02] + @test isapprox(disp, u_expected) +end diff --git a/test/test_elasticity_pyr5_point_load.jl b/test/test_elasticity_pyr5_point_load.jl new file mode 100644 index 0000000..3ac86dd --- /dev/null +++ b/test/test_elasticity_pyr5_point_load.jl @@ -0,0 +1,47 @@ +# This file is a part of JuliaFEM. +# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md + +using JuliaFEM +using JuliaFEM.Preprocess +using JuliaFEM.Testing + +@testset "test Pyr5 elasticity with point load" begin + nodes = Dict{Int64, Node}( + 1 => [-1.0,-1.0,-1.0], + 2 => [ 1.0,-1.0,-1.0], + 3 => [ 1.0, 1.0,-1.0], + 4 => [-1.0, 1.0,-1.0], + 5 => [ 0.0, 0.0, 1.0]) + + element1 = Element(Pyr5, [1, 2, 3, 4, 5]) + baseQuad = Element(Quad4, [1, 2, 3, 4]) + tipPoint = Element(Poi1, [5,]) + + update!([element1,baseQuad,tipPoint], "geometry", nodes) + update!([element1], "youngs modulus", 288.0) + update!([element1], "poissons ratio", 1/3) + + update!([tipPoint], "displacement traction force 1", 5.0) + update!([tipPoint], "displacement traction force 2", -7.0) + update!([tipPoint], "displacement traction force 3", 3.0) + + elasticity_problem = Problem(Elasticity, "solve continuum block", 3) + elasticity_problem.properties.finite_strain = false + push!(elasticity_problem, element1, baseQuad, tipPoint) + + baseQuad["displacement 1"] = 0.0 + baseQuad["displacement 2"] = 0.0 + baseQuad["displacement 3"] = 0.0 + boundary_problem = Problem(Dirichlet, "Boundary conditions", 3, "displacement") + push!(boundary_problem, baseQuad) + + solver = LinearSolver(elasticity_problem, boundary_problem) + solver() + + disp = element1("displacement", [0.0, 0.0, 1.0], 0.0) + info("########################################################") + info("displacement at tip: $disp") + # Code_Aster Result in verification/2017-05-27-pyramids/Pyr5_displacement.txt + u_expected = [6.9444444444427100E-02,-9.7222222222197952E-02,1.0416666666679683E-02] + @test isapprox(disp, u_expected) +end