diff --git a/test/runtests.jl b/test/runtests.jl
index d6e177d..ca75445 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -293,6 +293,54 @@ end
         C = edgecrossproduct(F,V) 
         @test C == [Vec3{Float64}(0.0,0.0,1.0),Vec3{Float64}(0.0,0.0,-1.0)]
     end
+    @testset "Mesh" begin
+        M = cube(1)
+        @test facenormal(M) == Vec3{Float64}[[0.0, 0.0, -1.0], [0.0, 0.0, 1.0], [-1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0], [0.0, -1.0, 0.0]]
+    end
+end
+
+@testset "facenormal" verbose = true begin
+    @testset "Single triangle" begin
+        F = [TriangleFace{Int64}(1, 2, 3)]
+        V = Vector{GeometryBasics.Point{3, Float64}}(undef,3)
+        V[1] = GeometryBasics.Point{3, Float64}(0.0, 0.0, 0.0)
+        V[2] = GeometryBasics.Point{3, Float64}(1.0, 0.0, 0.0)
+        V[3] = GeometryBasics.Point{3, Float64}(1.0, 1.0, 0.0)
+        C = facenormal(F,V) 
+        @test C == [Vec3{Float64}(0.0,0.0,1.0)]
+    end
+    @testset "Single quad" begin
+        F = [QuadFace{Int64}(1, 2, 3, 4)]
+        V = Vector{GeometryBasics.Point{3, Float64}}(undef,4)
+        V[1] = GeometryBasics.Point{3, Float64}(0.0, 0.0, 0.0)
+        V[2] = GeometryBasics.Point{3, Float64}(1.0, 0.0, 0.0)
+        V[3] = GeometryBasics.Point{3, Float64}(1.0, 1.0, 0.0)
+        V[4] = GeometryBasics.Point{3, Float64}(0.0, 1.0, 0.0)
+        C = facenormal(F,V) 
+        @test C == [Vec3{Float64}(0.0,0.0,1.0)]
+    end
+    @testset "Triangles" begin
+        F = [TriangleFace{Int64}(1, 2, 3),TriangleFace{Int64}(1, 4, 3)]
+        V = Vector{GeometryBasics.Point{3, Float64}}(undef,4)
+        V[1] = GeometryBasics.Point{3, Float64}(0.0, 0.0, 0.0)
+        V[2] = GeometryBasics.Point{3, Float64}(1.0, 0.0, 0.0)
+        V[3] = GeometryBasics.Point{3, Float64}(1.0, 1.0, 0.0)
+        V[4] = GeometryBasics.Point{3, Float64}(0.0, 1.0, 0.0)
+        C = facenormal(F,V) 
+        @test C == [Vec3{Float64}(0.0,0.0,1.0),Vec3{Float64}(0.0,0.0,-1.0)]
+    end
+    @testset "Quads" begin
+        F = [QuadFace{Int64}(1, 2, 3, 4),QuadFace{Int64}(6, 5, 4, 3)]
+        V = Vector{GeometryBasics.Point{3, Float64}}(undef,6)
+        V[1] = GeometryBasics.Point{3, Float64}(0.0, 0.0, 0.0)
+        V[2] = GeometryBasics.Point{3, Float64}(1.0, 0.0, 0.0)
+        V[3] = GeometryBasics.Point{3, Float64}(1.0, 1.0, 0.0)
+        V[4] = GeometryBasics.Point{3, Float64}(0.0, 1.0, 0.0)
+        V[5] = GeometryBasics.Point{3, Float64}(2.0, 0.0, 0.0)
+        V[6] = GeometryBasics.Point{3, Float64}(2.0, 1.0, 0.0)
+        C = facenormal(F,V) 
+        @test C == [Vec3{Float64}(0.0,0.0,1.0),Vec3{Float64}(0.0,0.0,-1.0)]
+    end
 end
 
 @testset "elements2indices" verbose = true begin