Merge pull request #56 from KristofferC/kc/tetras

add function space and gauss rule for tetraheder

src/utilities/VTK.jl

@@ -1,3 +1,28 @@
+function get_cell_type(nen, ndim)
+    # TODO: This is a bad way of figuring out the eltype
+    if nen == 3 && ndim == 2
+        cell =  VTKCellType.VTK_TRIANGLE
+    elseif nen == 4 && ndim == 2
+        cell = VTKCellType.VTK_QUAD
+    elseif nen == 4 && ndim == 2
+        cell = VTKCellType.VTK_HEXAHEDRON
+    elseif nen == 4 && ndim == 3
+        cell = VTKCellType.VTK_TETRA
+    end
+    return cell
+end
+
+function pad_zeros(points, ndim, nnodes)
+    if ndim == 3
+        points = points
+    elseif ndim == 2
+        points = [points; zeros(nnodes)']
+    elseif ndim == 1
+        points = [points; zeros(nnodes)'; zeros(nnodes)']
+    end
+    return points
+end
+
 """
     vtk_grid(Edof, Coord, Dof, nen, filename::AbstractString) -> vtkgrid
 Creates an unstructured VTK grid. `nen` is the number of nodes per element
@@ -12,23 +37,11 @@ function vtk_grid(Edof, Coord, Dof, nen, filename::AbstractString)
     nnodes = size(Coord, 2)
     ndim = size(Coord, 1)
 
-    # TODO: This is a bad way of figuring out the eltype
-    if nen == 3 && ndim == 2
-        cell =  VTKCellType.VTK_TRIANGLE
-    elseif nen == 4 && ndim == 2
-        cell = VTKCellType.VTK_QUAD
-    elseif nen == 4 && ndim == 2
-        cell = VTKCellType.VTK_HEXAHEDRON
-    end
+    cell = get_cell_type(nen, ndim)
 
     points = Coord
 
-    if ndim == 2
-        points = [points; zeros(nnodes)']
-    elseif ndim == 1
-        points = [points; zeros(nnodes)'; zeros(nnodes)']
-    end
-
+    points = pad_zeros(points, ndim, nnodes)
     cells = MeshCell[MeshCell(cell, top[:,i]) for i = 1:nele]
 
     vtk = vtk_grid(filename, points, cells)
@@ -48,22 +61,10 @@ function vtk_grid(topology::Matrix{Int}, Coord, filename::AbstractString)
     nnodes = size(Coord, 2)
     ndim = size(Coord, 1)
 
-    # TODO: This is a bad way of figuring out the eltype
-    if nen == 3 && ndim == 2
-        cell =  VTKCellType.VTK_TRIANGLE
-    elseif nen == 4 && ndim == 2
-        cell = VTKCellType.VTK_QUAD
-    elseif nen == 4 && ndim == 2
-        cell = VTKCellType.VTK_HEXAHEDRON
-    end
+    cell = get_cell_type(nen, ndim)
 
     points = Coord
-
-    if ndim == 2
-        points = [points; zeros(nnodes)']
-    elseif ndim == 1
-        points = [points; zeros(nnodes)'; zeros(nnodes)']
-    end
+    points = pad_zeros(points, ndim, nnodes)
 
     cells = MeshCell[MeshCell(cell, topology[:,i]) for i = 1:nele]
 

src/utilities/function_spaces.jl

@@ -177,9 +177,11 @@ function value!(fs::Lagrange{2, Triangle, 1}, N::Vector, ξ::Vec{2})
         ξ_x = ξ[1]
         ξ_y = ξ[2]
 
+        γ = 1 - ξ_x - ξ_y
+
         N[1] = ξ_x
         N[2] = ξ_y
-        N[3] = 1.0 - ξ_x - ξ_y
+        N[3] = γ
     end
 
     return N
@@ -261,6 +263,53 @@ function reference_coordinates(fs::Lagrange{2, Triangle, 2})
 end
 
 
+###################################
+# Lagrange dim 3 Triangle order 1 #
+###################################
+
+n_basefunctions(::Lagrange{3, Triangle, 1}) = 4
+
+function value!(fs::Lagrange{3, Triangle, 1}, N::Vector, ξ::Vec{3})
+    checkdim_value(fs, N, ξ)
+
+    @inbounds begin
+        ξ_x = ξ[1]
+        ξ_y = ξ[2]
+        ξ_z = ξ[3]
+
+        γ = 1. - ξ_x - ξ_y - ξ_z
+
+        N[1] = ξ_x
+        N[2] = ξ_y
+        N[3] = ξ_z
+        N[4] = γ
+    end
+
+    return N
+end
+
+function derivative!{T}(fs::Lagrange{3, Triangle, 1}, dN::Vector{Vec{3, T}}, ξ::Vec{3, T})
+    checkdim_derivative(fs, dN, ξ)
+
+    @inbounds begin
+        dN[1] = Vec{3, T}(( 1.0,  0.0, 0.0))
+        dN[2] = Vec{3, T}(( 0.0,  1.0, 0.0))
+        dN[3] = Vec{3, T}(( 0.0,  0.0, 1.0))
+        dN[4] = Vec{3, T}((-1.0, -1.0, -1.0))
+    end
+
+    return dN
+end
+
+function reference_coordinates(fs::Lagrange{3, Triangle, 1})
+    return (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.0, 0.0, 0.0)))
+end
+
+VTK_type(fs::Lagrange{3, Triangle, 1}) = VTKCellType.VTK_TETRA
+
 ###################################
 # Lagrange dim 3 Square order 1 #
 ###################################
@@ -376,4 +425,4 @@ function reference_coordinates(fs::Serendipity{2, Square, 2})
             Vec{2, Float64}(( 1.0,  0.0)),
             Vec{2, Float64}(( 0.0,  1.0)),
             Vec{2, Float64}((-1.0,  0.0)))
-end
+end

src/utilities/gaussquad_tet_table.jl

@@ -0,0 +1,30 @@
+# Patrick Keast, MODERATE-DEGREE TETRAHEDRAL QUADRATURE FORMULAS
+function _get_gauss_tetdata(n::Int)
+    if n == 1
+        a = 1. / 4.
+        w = 1. / 6.
+        xw = [a a a w]
+    elseif n == 2
+        a = ( 5. + 3. * √(5.) ) / 20.
+        b = ( 5. - √(5.) ) / 20.
+        w = 1. / 24.
+        xw = [a b b w
+              b a b w
+              b b a w
+              b b b w]
+    elseif n == 3
+        a1 = 1. / 4.
+        a2 = 1. / 2.
+        b2 = 1. / 6.
+        w1 = -2. / 15.
+        w2 = 3. / 40.
+        xw = [a1 a1 a1 w1
+              a2 b2 b2 w2
+              b2 a2 b2 w2
+              b2 b2 a2 w2
+              b2 b2 b2 w2]
+    else
+        throw(ArgumentError("Unsupported tetraheder gauss integration"))
+    end
+    return xw
+end

src/utilities/quadrature.jl

@@ -1,3 +1,6 @@
+include("gaussquad_tri_table.jl")
+include("gaussquad_tet_table.jl")
+
 """
 A type that perform quadrature integration
 """
@@ -33,6 +36,14 @@ function get_gaussrule(::Type{Dim{2}}, ::Square, order::Int)
     end
 end
 
+function get_gaussrule(::Type{Dim{3}}, ::Triangle, order::Int)
+    if order <= 3
+        return tetrules[order]
+    else
+        return make_tetrule(order)
+    end
+end
+
 function get_gaussrule(::Type{Dim{3}}, ::Square, order::Int)
     if order <= 5
         return cuberules[order]
@@ -60,13 +71,7 @@ function make_cuberule(order::Int)
 end
 
 const cuberules = [make_cuberule(i) for i = 1:5]
-function get_cuberule(order::Int)
-    if order <= 5
-        return cuberules[order]
-    else
-        return make_cuberule(order)
-    end
-end
+
 
 """
 Creates a `QuadratureRule` that integrates
@@ -86,13 +91,7 @@ function make_quadrule(order::Int)
 end
 
 const quadrules = [make_quadrule(i) for i = 1:5]
-function get_quadrule(order::Int)
-    if order <= 5
-        return quadrules[order]
-    else
-        return make_quadrule(order)
-    end
-end
+
 
 """
 Creates a `QuadratureRule` that integrates
@@ -108,15 +107,6 @@ function make_linerule(order::Int)
 end
 
 const linerules = [make_linerule(i) for i = 1:5]
-function get_linerule(order::Int)
-    if order <= 5
-        return linerules[order]
-    else
-        return make_linerule(order)
-    end
-end
-
-include("gaussquad_tri_table.jl")
 
 function make_trirule(order::Int)
     data = _get_gauss_tridata(order)
@@ -136,10 +126,22 @@ end
 
 const trirules = [make_trirule(i) for i = 1:5]
 
-function get_trirule(order::Int)
-    if order <= 5
-        return trirules[order]
-    else
-        return make_trirule(order)
+
+function make_tetrule(order::Int)
+    data = _get_gauss_tetdata(order)
+    n_points = size(data,1)
+    weights = Array(Float64, n_points)
+    points = Array(Vec{3, Float64}, n_points)
+
+    for p in 1:size(data, 1)
+        points[p] = Vec{3, Float64}((data[p, 1], data[p, 2], data[p, 3]))
     end
+
+    weights = data[:, 4]
+
+    QuadratureRule(weights, points)
 end
+
+const tetrules = [make_tetrule(i) for i = 1:3]
+
+

test/test_fevalues.jl

@@ -86,7 +86,8 @@ end
                                          (Lagrange{2, Triangle, 1}(), get_gaussrule(Dim{2}, Triangle(), 2)),
                                          (Lagrange{2, Triangle, 2}(), get_gaussrule(Dim{2}, Triangle(), 2)),
                                          (Lagrange{3, Square, 1}(), get_gaussrule(Dim{3}, Square(), 2)),
-                                         (Serendipity{2, Square, 2}(), get_gaussrule(Dim{2}, Square(), 2)))
+                                         (Serendipity{2, Square, 2}(), get_gaussrule(Dim{2}, Square(), 2)),
+                                         (Lagrange{3, Triangle, 1}(), get_gaussrule(Dim{3}, Triangle(), 2)))
 
         fev = FEValues(quad_rule, function_space)
         ndim = n_dim(function_space)

test/test_utilities.jl

@@ -150,7 +150,8 @@ end
                           Lagrange{2, Triangle, 1}(),
                           Lagrange{2, Triangle, 2}(),
                           Lagrange{3, Square, 1}(),
-                          Serendipity{2, Square, 2}())
+                          Serendipity{2, Square, 2}(),
+                          Lagrange{3, Triangle, 1}())
         ndim = n_dim(functionspace)
         n_basefuncs = n_basefunctions(functionspace)
         x = rand(Tensor{1, ndim})