diff --git a/src/Export/VTK.jl b/src/Export/VTK.jl index aafa66ecef..c60aaad666 100644 --- a/src/Export/VTK.jl +++ b/src/Export/VTK.jl @@ -90,6 +90,7 @@ getVTKtype(::Serendipity{2, RefCube, 2}) = VTKCellTypes.VTK_QUADRATIC_QUAD getVTKtype(::Lagrange{3, RefCube, 1}) = VTKCellTypes.VTK_HEXAHEDRON getVTKtype(::Lagrange{3, RefTetrahedron, 1}) = VTKCellTypes.VTK_TETRA +getVTKtype(::Lagrange{3, RefTetrahedron, 2}) = VTKCellTypes.VTK_QUADRATIC_TETRA getVTKtype(::Type{Line}) = VTKCellTypes.VTK_LINE getVTKtype(::Type{QuadraticLine}) = VTKCellTypes.VTK_QUADRATIC_EDGE @@ -99,6 +100,7 @@ getVTKtype(::Type{QuadraticTriangle}) = VTKCellTypes.VTK_QUADRATIC_TRIANGLE getVTKtype(::Type{QuadraticQuadrilateral}) = VTKCellTypes.VTK_BIQUADRATIC_QUAD getVTKtype(::Type{Hexahedron}) = VTKCellTypes.VTK_HEXAHEDRON getVTKtype(::Type{Tetrahedron}) = VTKCellTypes.VTK_TETRA +getVTKtype(::Type{QuadraticTetrahedron}) = VTKCellTypes.VTK_QUADRATIC_TETRA getVTKtype(::Type{Cell{2,8,4}}) = VTKCellTypes.VTK_QUADRATIC_QUAD diff --git a/src/Grid/grid.jl b/src/Grid/grid.jl index f92143b7ae..d880bc1a6f 100644 --- a/src/Grid/grid.jl +++ b/src/Grid/grid.jl @@ -45,7 +45,7 @@ nfaces{dim, N, M}(::Type{Cell{dim, N, M}}) = M @compat const QuadraticQuadrilateral = Cell{2, 9, 4} @compat const Tetrahedron = Cell{3, 4, 4} -@compat const QuadraticTetrahedron = Cell{3, 10, 4} # Function interpolation for this doesn't exist in JuAFEM yet +@compat const QuadraticTetrahedron = Cell{3, 10, 4} @compat const Hexahedron = Cell{3, 8, 6} @compat const QuadraticHexahedron = Cell{3, 20, 6} # Function interpolation for this doesn't exist in JuAFEM yet diff --git a/src/interpolations.jl b/src/interpolations.jl index 0bb9f29ae7..8e9f5d5f15 100644 --- a/src/interpolations.jl +++ b/src/interpolations.jl @@ -25,6 +25,7 @@ The following interpolations are implemented: * `Lagrange{3, RefCube, 1}` * `Serendipity{2, RefCube, 2}` * `Lagrange{3, RefTetrahedron, 1}` +* `Lagrange{3, RefTetrahedron, 2}` **Common methods:** @@ -405,6 +406,61 @@ function derivative!{T}(ip::Lagrange{3, RefTetrahedron, 1}, dN::AbstractVector{V return dN end +######################################### +# Lagrange dim 3 RefTetrahedron order 2 # +######################################### +getnbasefunctions(::Lagrange{3, RefTetrahedron, 2}) = 10 +getnfacenodes(::Lagrange{3, RefTetrahedron, 2}) = 6 +getfacelist(::Lagrange{3, RefTetrahedron, 2}) = ((1,2,3,5,6,7),(1,2,4,5,8,9),(2,3,4,6,9,10),(1,3,4,7,8,10)) + +# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch09.d/AFEM.Ch09.pdf +# http://www.colorado.edu/engineering/CAS/courses.d/AFEM.d/AFEM.Ch10.d/AFEM.Ch10.pdf +function value!(ip::Lagrange{3, RefTetrahedron, 2}, N::AbstractVector, ξ::Vec{3}) + checkdim_value(ip, N, ξ) + + @inbounds begin + ξ_x = ξ[1] + ξ_y = ξ[2] + ξ_z = ξ[3] + + N[1] = (-2 * ξ_x - 2 * ξ_y - 2 * ξ_z + 1) * (-ξ_x - ξ_y - ξ_z + 1) + N[2] = ξ_x * (2 * ξ_x - 1) + N[3] = ξ_y * (2 * ξ_y - 1) + N[4] = ξ_z * (2 * ξ_z - 1) + N[5] = ξ_x * (-4 * ξ_x - 4 * ξ_y - 4 * ξ_z + 4) + N[6] = 4 * ξ_x * ξ_y + N[7] = 4 * ξ_y * (-ξ_x - ξ_y - ξ_z + 1) + N[8] = ξ_z * (-4 * ξ_x - 4 * ξ_y - 4 * ξ_z + 4) + N[9] = 4 * ξ_x * ξ_z + N[10] = 4 * ξ_y * ξ_z + end + + return N +end + +function derivative!{T}(ip::Lagrange{3, RefTetrahedron, 2}, dN::AbstractVector{Vec{3, T}}, ξ::Vec{3, T}) + checkdim_derivative(ip, dN, ξ) + + @inbounds begin + ξ_x = ξ[1] + ξ_y = ξ[2] + ξ_z = ξ[3] + + dN[1] = Vec{3, T}((4 * ξ_x + 4 * ξ_y + 4 * ξ_z - 3, 4 * ξ_x + 4 * ξ_y + 4 * ξ_z - 3, 4 * ξ_x + 4 * ξ_y + 4 * ξ_z - 3)) + dN[2] = Vec{3, T}((4 * ξ_x - 1, 0.0, 0.0)) + dN[3] = Vec{3, T}((0.0, 4 * ξ_y - 1, 0.0)) + dN[4] = Vec{3, T}((0.0, 0.0, 4 * ξ_z - 1)) + dN[5] = Vec{3, T}((-8 * ξ_x - 4 * ξ_y - 4 * ξ_z + 4, -4 * ξ_x, -4 * ξ_x)) + dN[6] = Vec{3, T}((4 * ξ_y, 4 * ξ_x, 0.0)) + dN[7] = Vec{3, T}((-4 * ξ_y, -4 * ξ_x - 8 * ξ_y - 4 * ξ_z + 4, -4 * ξ_y)) + dN[8] = Vec{3, T}((-4 * ξ_z, -4 * ξ_z, -4 * ξ_x - 4 * ξ_y - 8 * ξ_z + 4)) + dN[9] = Vec{3, T}((4 * ξ_z, 0.0, 4 * ξ_x)) + dN[10] = Vec{3, T}((0.0, 4 * ξ_z, 4 * ξ_y)) + end + + return dN +end + ################################## # Lagrange dim 3 RefCube order 1 # ################################## diff --git a/test/test_VTK.jl b/test/test_VTK.jl index 6f450cc04a..5c4fa7e39c 100644 --- a/test/test_VTK.jl +++ b/test/test_VTK.jl @@ -73,7 +73,8 @@ using SHA Lagrange{2, RefTetrahedron, 2}(), Lagrange{3, RefCube, 1}(), Serendipity{2, RefCube, 2}(), - Lagrange{3, RefTetrahedron, 1}()) + Lagrange{3, RefTetrahedron, 1}(), + Lagrange{3, RefTetrahedron, 2}()) @test getVTKtype(interpolation).nodes == getnbasefunctions(interpolation) end diff --git a/test/test_cellvalues.jl b/test/test_cellvalues.jl index df80e123ae..f78f978e8e 100644 --- a/test/test_cellvalues.jl +++ b/test/test_cellvalues.jl @@ -1,6 +1,6 @@ @testset "CellValues" begin - -for (func_interpol, quad_rule) in ((Lagrange{1, RefCube, 1}(), QuadratureRule{1, RefCube}(2)), +for (func_interpol, quad_rule) in ( + (Lagrange{1, RefCube, 1}(), QuadratureRule{1, RefCube}(2)), (Lagrange{1, RefCube, 2}(), QuadratureRule{1, RefCube}(2)), (Lagrange{2, RefCube, 1}(), QuadratureRule{2, RefCube}(2)), (Lagrange{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)), @@ -8,8 +8,9 @@ for (func_interpol, quad_rule) in ((Lagrange{1, RefCube, 1}(), QuadratureRule{1 (Lagrange{2, RefTetrahedron, 2}(), QuadratureRule{2, RefTetrahedron}(2)), (Lagrange{3, RefCube, 1}(), QuadratureRule{3, RefCube}(2)), (Serendipity{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)), - (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{3, RefTetrahedron}(2))) - + (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{3, RefTetrahedron}(2)), + (Lagrange{3, RefTetrahedron, 2}(), QuadratureRule{3, RefTetrahedron}(2)) + ) for fe_valtype in (CellScalarValues, CellVectorValues) cv = fe_valtype(quad_rule, func_interpol) diff --git a/test/test_facevalues.jl b/test/test_facevalues.jl index 02cd47bd63..ba67724f4c 100644 --- a/test/test_facevalues.jl +++ b/test/test_facevalues.jl @@ -8,8 +8,9 @@ for (func_interpol, quad_rule) in ( (Lagrange{2, RefTetrahedron, 2}(), QuadratureRule{1, RefTetrahedron}(2)), (Lagrange{3, RefCube, 1}(), QuadratureRule{2, RefCube}(2)), (Serendipity{2, RefCube, 2}(), QuadratureRule{1, RefCube}(2)), - (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)) - ) + (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)), + (Lagrange{3, RefTetrahedron, 2}(), QuadratureRule{2, RefTetrahedron}(2)) + ) for fe_valtype in (FaceScalarValues, FaceVectorValues) fv = fe_valtype(quad_rule, func_interpol) diff --git a/test/test_interpolations.jl b/test/test_interpolations.jl index 683eacf2f8..055c59175f 100644 --- a/test/test_interpolations.jl +++ b/test/test_interpolations.jl @@ -8,7 +8,8 @@ for interpolation in (Lagrange{1, RefCube, 1}(), Lagrange{2, RefTetrahedron, 2}(), Lagrange{3, RefCube, 1}(), Serendipity{2, RefCube, 2}(), - Lagrange{3, RefTetrahedron, 1}()) + Lagrange{3, RefTetrahedron, 1}(), + Lagrange{3, RefTetrahedron, 2}()) # Test of utility functions ndim = getdim(interpolation) diff --git a/test/test_utils.jl b/test/test_utils.jl index 6c7999e95a..a17af05b68 100644 --- a/test/test_utils.jl +++ b/test/test_utils.jl @@ -82,7 +82,7 @@ function reference_coordinates(::Lagrange{2, RefTetrahedron, 2}) end # Lagrange{3, RefTetrahedron} -function reference_normals(::Lagrange{3, RefTetrahedron, 1}) +function reference_normals(::Lagrange{3, RefTetrahedron}) return [Vec{3, Float64}((0.0, 0.0, -1.0)), Vec{3, Float64}((0.0 ,-1.0, 0.0)), Vec{3, Float64}((1/√3, 1/√3, 1/√3)), @@ -96,6 +96,19 @@ function reference_coordinates(::Lagrange{3, RefTetrahedron, 1}) Vec{3, Float64}((0.0, 0.0, 1.0))] end +function reference_coordinates(::Lagrange{3, RefTetrahedron, 2}) + return [Vec{3, Float64}((0.0, 0.0, 0.0)), + Vec{3, Float64}((1.0, 0.0, 0.0)), + Vec{3, Float64}((0.0, 1.0, 0.0)), + Vec{3, Float64}((0.0, 0.0, 1.0)), + Vec{3, Float64}((0.5, 0.0, 0.0)), + Vec{3, Float64}((0.5, 0.5, 0.0)), + Vec{3, Float64}((0.0, 0.5, 0.0)), + Vec{3, Float64}((0.0, 0.0, 0.5)), + Vec{3, Float64}((0.5, 0.0, 0.5)), + Vec{3, Float64}((0.0, 0.5, 0.5))] +end + # Lagrange{3, Cube} function reference_normals(::Lagrange{3, RefCube, 1}) return [Vec{3, Float64}(( 0.0, 0.0, -1.0)), @@ -201,7 +214,7 @@ function calculate_volume{T, dim}(::Lagrange{2, RefTetrahedron, 2}, x::Vector{Ve return vol end -function calculate_volume{T}(::Lagrange{3, RefTetrahedron, 1}, x::Vector{Vec{3, T}}) +function calculate_volume{T, order}(::Lagrange{3, RefTetrahedron, order}, x::Vector{Vec{3, T}}) vol = norm((x[2] - x[1]) ⋅ ((x[3] - x[1]) × (x[4] - x[1]))) / 6.0 return vol end @@ -279,3 +292,9 @@ function topology_test_nodes(::Lagrange{3, RefTetrahedron, 1}) face_nodes = [[3,4,8], [3,4,1], [4,8,1], [3,8,1], [1,4,1337]] return face_nodes, cell_nodes end + +function topology_test_nodes(::Lagrange{3, RefTetrahedron, 2}) + cell_nodes = [3,4,8,1,5,2,6,7,9,10] + face_nodes = [[3,4,8,5,2,6], [3,4,1,5,9,7], [4,8,1,2,10,9], [3,8,1,6,10,7], [1,4,1337,2,3,5]] + return face_nodes, cell_nodes +end