add support for quadratic tetrahedron (#129)

Ferrite-FEM · Apr 4, 2017 · 2e3c7c3 · 2e3c7c3
1 parent 31d38b1
commit 2e3c7c3
Show file tree

Hide file tree

Showing 8 changed files with 92 additions and 11 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -1,15 +1,16 @@
 @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)),
                                     (Lagrange{2, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)),
                                     (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
@@ -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
@@ -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
@@ -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