diff --git a/Project.toml b/Project.toml
index 220c5f14..2a68b39e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.4.2"
 EarCut_jll = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
diff --git a/src/GeometryBasics.jl b/src/GeometryBasics.jl
index fb60f3a8..69c13569 100644
--- a/src/GeometryBasics.jl
+++ b/src/GeometryBasics.jl
@@ -1,10 +1,12 @@
 module GeometryBasics
 
-using StaticArrays, Tables, StructArrays, IterTools, LinearAlgebra
+using Tables, StructArrays, IterTools, LinearAlgebra
 using EarCut_jll
+import Base: *
 
 using Base: @propagate_inbounds
 
+include("mat.jl")
 include("fixed_arrays.jl")
 include("offsetintegers.jl")
 include("basic_types.jl")
@@ -24,8 +26,6 @@ include("triangulation.jl")
 include("lines.jl")
 include("boundingboxes.jl")
 
-include("deprecated.jl")
-
 export AbstractGeometry, GeometryPrimitive
 export Mat, Point, Vec
 export LineFace, Polytope, Line, NgonFace, convert_simplex
@@ -43,18 +43,12 @@ export PolygonMeta, MultiPointMeta, MultiLineStringMeta, MeshMeta, LineStringMet
 export decompose, coordinates, faces, normals, decompose_uv, decompose_normals,
        texturecoordinates
 export Tesselation, pointmeta, Normal, UV, UVW
-export GLTriangleFace, GLNormalMesh3D, GLPlainTriangleMesh, GLUVMesh3D, GLUVNormalMesh3D
-export AbstractMesh, Mesh, TriangleMesh
-export GLNormalMesh2D, PlainTriangleMesh
+export AbstractMesh, Mesh
 export MetaT, meta_table
 
 # all the different predefined mesh types
 # Note: meshes can contain arbitrary meta information,
-export AbstractMesh, TriangleMesh, PlainMesh, GLPlainMesh, GLPlainMesh2D, GLPlainMesh3D
-export UVMesh, GLUVMesh, GLUVMesh2D, GLUVMesh3D
-export NormalMesh, GLNormalMesh, GLNormalMesh2D, GLNormalMesh3D
-export NormalUVMesh, GLNormalUVMesh, GLNormalUVMesh2D, GLNormalUVMesh3D
-export NormalUVWMesh, GLNormalUVWMesh, GLNormalUVWMesh2D, GLNormalUVWMesh3D
+export AbstractMesh, TriangleMesh, PlainMesh
 
 # mesh creation functions
 export triangle_mesh, triangle_mesh, uv_mesh
@@ -71,7 +65,7 @@ export min_euclideansq, minmax_dist_dim, minmax_euclidean, minmax_euclideansq
 export self_intersections, split_intersections
 
 if Base.VERSION >= v"1.4.2"
-    include("precompile.jl")
+    include("precompiles.jl")
     _precompile_()
 end
 
diff --git a/src/basic_types.jl b/src/basic_types.jl
index 16aa8047..c6123f72 100644
--- a/src/basic_types.jl
+++ b/src/basic_types.jl
@@ -13,26 +13,19 @@ Note That `Polytope{N} where N == 3` denotes a Triangle both as a Simplex or Ngo
 abstract type Polytope{Dim,T} <: AbstractGeometry{Dim,T} end
 abstract type AbstractPolygon{Dim,T} <: Polytope{Dim,T} end
 
-abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
 abstract type AbstractFace{N,T} <: StaticVector{N,T} end
 abstract type AbstractSimplexFace{N,T} <: AbstractFace{N,T} end
 abstract type AbstractNgonFace{N,T} <: AbstractFace{N,T} end
 
 abstract type AbstractSimplex{Dim,N,T} <: StaticVector{Dim,T} end
 
-"""
-Face index, connecting points to form a simplex
-"""
-
 @fixed_vector SimplexFace AbstractSimplexFace
+
 const TetrahedronFace{T} = SimplexFace{4,T}
 Face(::Type{<:SimplexFace{N}}, ::Type{T}) where {N,T} = SimplexFace{N,T}
 
-"""
-Face index, connecting points to form an Ngon
-"""
-
 @fixed_vector NgonFace AbstractNgonFace
+
 const LineFace{T} = NgonFace{2,T}
 const TriangleFace{T} = NgonFace{3,T}
 const QuadFace{T} = NgonFace{4,T}
@@ -58,14 +51,14 @@ Fixed Size Polygon, e.g.
 """
 struct Ngon{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: AbstractPolygon{Dim,T}
 
-    points::SVector{N,Point}
+    points::Vec{N,Point}
 end
 
 const NNgon{N} = Ngon{Dim,T,N,P} where {Dim,T,P}
 
 function (::Type{<:NNgon{N}})(points::Vararg{P,N}) where {P<:AbstractPoint{Dim,T},
                                                           N} where {Dim,T}
-    return Ngon{Dim,T,N,P}(SVector(points))
+    return Ngon{Dim,T,N,P}(Vec(points))
 end
 Base.show(io::IO, x::NNgon{N}) where {N} = print(io, "Ngon{$N}(", join(x, ", "), ")")
 
@@ -138,7 +131,7 @@ to allow embedding in higher-order spaces by parameterizing on `T`.
 """
 struct Simplex{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: Polytope{Dim,T}
 
-    points::SVector{N,Point}
+    points::Vec{N,Point}
 end
 
 const NSimplex{N} = Simplex{Dim,T,N,P} where {Dim,T,P}
@@ -151,7 +144,7 @@ coordinates(x::Simplex) = x.points
 
 function (::Type{<:NSimplex{N}})(points::Vararg{P,N}) where {P<:AbstractPoint{Dim,T},
                                                              N} where {Dim,T}
-    return Simplex{Dim,T,N,P}(SVector(points))
+    return Simplex{Dim,T,N,P}(Vec(points))
 end
 
 # Base Array interface
diff --git a/src/deprecated.jl b/src/deprecated.jl
deleted file mode 100644
index 54336aa0..00000000
--- a/src/deprecated.jl
+++ /dev/null
@@ -1,30 +0,0 @@
-using Base: @deprecate_binding
-
-# Types ...f0 renamed to ...f
-@deprecate_binding Vecf0 Vecf
-@deprecate_binding Pointf0 Pointf
-for i in 1:4
-    for T in [:Point, :Vec]
-        oldname = Symbol("$T$(i)f0")
-        newname = Symbol("$T$(i)f")
-        @eval begin
-            @deprecate_binding $oldname $newname
-        end
-    end
-    oldname = Symbol("Mat$(i)f0")
-    newname = Symbol("Mat$(i)f")
-    @eval begin
-        @deprecate_binding $oldname $newname
-    end
-end
-
-# Rect types
-@deprecate_binding Rect2D  Rect2
-@deprecate_binding Rect3D  Rect3
-@deprecate_binding FRect   Rectf
-@deprecate_binding FRect2D Rect2f
-@deprecate_binding FRect3D Rect3f
-@deprecate_binding IRect   Recti
-@deprecate_binding IRect2D Rect2i
-@deprecate_binding IRect3D Rect3i
-@deprecate_binding TRect   RectT
diff --git a/src/fixed_arrays.jl b/src/fixed_arrays.jl
index c312b829..db9b9215 100644
--- a/src/fixed_arrays.jl
+++ b/src/fixed_arrays.jl
@@ -1,127 +1,237 @@
 
-function unit(::Type{T}, i::Integer) where {T <: StaticVector}
-    tup = ntuple(Val(length(T))) do j
-        return ifelse(i == j, 1, 0)
-    end
-    return T(tup)
-end
+using LinearAlgebra
+import Random
+import Base: setindex
+
+abstract type StaticVector{N, T} end
 
-macro fixed_vector(name, parent)
+macro fixed_vector(VecT, SuperT)
     expr = quote
-        struct $(name){S,T} <: $(parent){S,T}
-            data::NTuple{S,T}
+        struct $(VecT){N, T} <: $(SuperT){N, T}
+            data::NTuple{N,T}
 
-            function $(name){S,T}(x::NTuple{S,T}) where {S,T}
-                return new{S,T}(x)
+            function $(VecT){N,T}(x::NTuple{N,T}) where {N,T}
+                return new{N,T}(x)
             end
 
-            function $(name){S,T}(x::NTuple{S,Any}) where {S,T}
-                return new{S,T}(StaticArrays.convert_ntuple(T, x))
+            function $(VecT){N,T}(x::NTuple{N,Any}) where {N,T}
+                return new{N,T}(convert(NTuple{N, T}, x))
             end
+            $(VecT){1, T}(x::Tuple{T}) where {T} = new{1, T}(x)
         end
 
-        size_or(::Type{$(name)}, or) = or
-        eltype_or(::Type{$(name)}, or) = or
-        eltype_or(::Type{$(name){S,T} where S}, or) where {T} = T
-        eltype_or(::Type{$(name){S,T} where T}, or) where {S} = or
-        eltype_or(::Type{$(name){S,T}}, or) where {S,T} = T
-
-        size_or(::Type{$(name){S,T} where S}, or) where {T} = or
-        size_or(::Type{$(name){S,T} where T}, or) where {S} = Size{(S,)}()
-        size_or(::Type{$(name){S,T}}, or) where {S,T} = (S,)
-
         # Array constructor
-        function $(name){S}(x::AbstractVector{T}) where {S,T}
-            @assert S <= length(x)
-            return $(name){S,T}(ntuple(i -> x[i], Val(S)))
+        $(VecT)(x::AbstractVector) = error("You need to supply size of vector")
+        $(VecT){N}(x::AbstractVector{T}) where {N,T} = $(VecT){N,T}(x)
+        function $(VecT){N,T1}(x::AbstractVector{T2}) where {N,T1,T2}
+            @assert N <= length(x)
+            return $(VecT){N,T1}(ntuple(i -> convert(T1, x[i]), N))
         end
 
-        function $(name){S,T1}(x::AbstractVector{T2}) where {S,T1,T2}
-            @assert S <= length(x)
-            return $(name){S,T1}(ntuple(i -> convert(T1, x[i]), Val(S)))
+        # StaticVector conversion
+        $(VecT)(x::StaticVector{N, T}) where {N,T} = $(VecT){N, T}(x)
+        $(VecT){N}(x::StaticVector{N2, T}) where {N,N2,T} = $(VecT){N,T}(x)
+        function $(VecT){N1,T1}(x::StaticVector{N2, T2}) where {N1,T1,N2,T2}
+            @assert N1 <= N2
+            return $(VecT){N1,T1}(ntuple(i -> convert(T1, x[i]), N1))
         end
-
-        function $(name){S,T}(x) where {S,T}
-            return $(name){S,T}(ntuple(i -> convert(T, x), Val(S)))
+        function $(VecT){1,T1}(x::StaticVector{N2, T2}) where {T1,N2,T2}
+            @assert 1 <= N2
+            return $(VecT){1,T1}(ntuple(i -> convert(T1, x[i]), 1))
         end
 
-        $(name){S}(x::T) where {S,T} = $(name){S,T}(ntuple(i -> x, Val(S)))
-        $(name){1,T}(x::T) where {T} = $(name){1,T}((x,))
-        $(name)(x::NTuple{S}) where {S} = $(name){S}(x)
-        function $(name)(x::T) where {S,T <: Tuple{Vararg{Any,S}}}
-            return $(name){S,StaticArrays.promote_tuple_eltype(T)}(x)
+        # repeat
+        $(VecT){N}(x::T) where {N,T<:Number} = $(VecT){N, T}(x)
+        function $(VecT){N,T}(x::Number) where {N,T}
+            return $(VecT){N,T}(ntuple(i -> convert(T, x), N))
         end
 
-        function $(name){S}(x::T) where {S,T <: Tuple}
-            return $(name){S,StaticArrays.promote_tuple_eltype(T)}(x)
-        end
-        $(name){S,T}(x::StaticVector) where {S,T} = $(name){S,T}(Tuple(x))
+        $(VecT){1, T}(x) where {T} = $(VecT){1, T}((x,))
+        $(VecT){1, T}(x::Tuple{Any}) where T = $(VecT){1, T}((T(x[1]),))
+        $(VecT)(x::Tuple) = $(VecT)(promote(x...))
+        $(VecT){N}(x::Tuple) where {N} = $(VecT){N}(promote(x...))
+        $(VecT){N, T}(x::Tuple) where {N,T} = $(VecT){N,T}(convert(NTuple{N,T}, x))
 
-        @generated function (::Type{$(name){S,T}})(x::$(name)) where {S,T}
-            idx = [:(x[$i]) for i in 1:S]
-            return quote
-                $($(name)){S,T}($(idx...))
-            end
+        $(VecT)(x::NTuple{N, T}) where {N, T} = $(VecT){N,T}(x)
+        $(VecT){N}(x::NTuple{N, T}) where {N, T} = $(VecT){N,T}(x)
+
+        $(VecT)(x::Vararg{<:Any,N}) where {N} = $(VecT){N}(x)
+        $(VecT)(x::Vararg{T,N}) where {T,N} = $(VecT){N,T}(x)
+
+        $(VecT){N}(x::Vararg{<:Any,N}) where {N} = $(VecT){N}(x)
+        $(VecT){N}(x::Vararg{T,N}) where {T,N} = $(VecT){N,T}(x)
+
+        $(VecT){N, T}(x::Vararg{<:Any,N}) where {T,N} = $(VecT){N,T}(x)
+        $(VecT){N, T1}(x::Vararg{T2,N}) where {T1,T2,N} = $(VecT){N, T1}(x)
+
+        Base.convert(::Type{$(VecT){N,T}}, x) where {N,T} = $(VecT){N,T}(x)
+        Base.convert(::Type{$(VecT){N}}, x) where {N} = $(VecT){N}(x)
+        Base.convert(::Type{$(VecT){N}}, x::$(VecT){N}) where {N} = x
+        Base.convert(::Type{$(VecT){N,T}}, x::$(VecT){N,T}) where {N,T} = x
+
+        function Base.convert(::Type{$(VecT){N,T}}, x::NTuple{N,T}) where {N,T}
+            return $(VecT){N,T}(x)
+        end
+        function Base.convert(::Type{$(VecT){N,T}}, x::Tuple) where {N,T}
+            return $(VecT){N,T}(convert(NTuple{N,T}, x))
         end
+        Base.convert(::Type{$(VecT)}, x::Tuple) = $(VecT)(x)
 
-        @generated function Base.convert(::Type{$(name){S,T}}, x::$(name)) where {S,T}
-            idx = [:(x[$i]) for i in 1:S]
-            return quote
-                $($(name)){S,T}($(idx...))
-            end
+        Base.@propagate_inbounds function Base.getindex(v::$(VecT){N,T}, i::Int) where {N,T}
+            return v.data[i]
+        end
+        Base.setindex(c::$(VecT){N, T}, v, i::Integer) where {N,T} = $(VecT){N,T}(Base.setindex(c.data, v, i))
+        Base.@propagate_inbounds function Base.getindex(a::AbstractArray{T}, idx::$(VecT){N, <:Integer}) where {N,T}
+            return $(VecT){N,T}(map(i-> a[i], idx))
         end
 
-        @generated function (::Type{SV})(x::StaticVector) where {SV <: $(name)}
-            len = size_or(SV, size(x))[1]
-            return if length(x) == len
-                :(SV(Tuple(x)))
-            elseif length(x) > len
-                elems = [:(x[$i]) for i in 1:len]
-                :(SV($(Expr(:tuple, elems...))))
-            else
-                error("Static Vector too short: $x, target type: $SV")
-            end
+        Base.@propagate_inbounds function Base.getindex(a::StaticVector{N1, T}, idx::$(VecT){N, <:Integer}) where {N,N1,T}
+            return $(VecT){N,T}(map(i-> a[i], idx))
         end
 
-        Base.@pure StaticArrays.Size(::Type{$(name){S,Any}}) where {S} = Size(S)
-        Base.@pure StaticArrays.Size(::Type{$(name){S,T}}) where {S,T} = Size(S)
+        Base.Tuple(v::$(VecT)) = v.data
 
-        Base.@propagate_inbounds function Base.getindex(v::$(name){S,T}, i::Int) where {S,T}
-            return v.data[i]
+
+        function Base.broadcasted(f, a::AbstractArray{T}, b::$(VecT)) where {T <: $(VecT)}
+            return Base.broadcasted(f, a, (b,))
         end
 
-        Base.Tuple(v::$(name)) = v.data
-        function Base.convert(::Type{$(name){S,T}}, x::NTuple{S,T}) where {S,T}
-            return $(name){S,T}(x)
+        function Base.broadcasted(f, a::$(VecT), b::$(VecT))
+            return $(VecT)(map(f, a.data, b.data))
         end
-        function Base.convert(::Type{$(name){S,T}}, x::Tuple) where {S,T}
-            return $(name){S,T}(convert(NTuple{S,T}, x))
+
+        Base.broadcasted(f, a::$(VecT)) = $(VecT)(f.(a.data))
+        Base.broadcasted(f, a::$(VecT), b) = $(VecT)(f.(a.data, b))
+        Base.broadcasted(f, a, b::$(VecT)) = $(VecT)(f.(a, b.data))
+
+        Base.map(f, b::$(VecT)) = $(VecT)(map(f, b.data))
+
+        (*)(a::Mat{M, N, T1}, b::$(VecT){O, T2}) where {T1, T2, M, N, O} = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))
+
+        # vector * (row vector)
+        @generated function *(a::$VecT{N, T1}, b::Mat{1, M, T2}) where {T1, T2, N, M}
+            elements = Expr(:tuple, [Expr(:tuple, [:(a[$i] * b[$j]) for i in 1:N]...) for j in 1:M]...)
+            return :($($(VecT))($elements))
         end
 
-        @generated function StaticArrays.similar_type(::Type{SV}, ::Type{T},
-                                                      s::Size{S}) where {SV <: $(name),T,S}
-            return if length(S) === 1
-                $(name){S[1],T}
+        # matrix * vector
+        @generated function *(a::Mat{M, N, T1}, b::$VecT{N, T2}) where {T1, T2, M, N}
+            total_terms = M*N
+            if total_terms <= 64
+                # Full unrolling
+                elements = Expr(:tuple, [Expr(:call, :+, [:(a[$i,$k]*b[$k]) for k = 1:N]...) for i in 1:M]...)
             else
-                StaticArrays.default_similar_type(T, s(), Val{length(S)})
+                # Expand as a bunch of dot products
+                elements = Expr(:tuple, [:(bilindot($($(VecT))(row(a,$i)),b)) for i in 1:M]...)
             end
+            return :($($(VecT))($elements))
         end
 
-        Base.:(*)(a::$name, b::$name) = a .* b
-        function Base.broadcasted(f, a::AbstractArray{T}, b::$name) where {T <: $name}
-            return Base.broadcasted(f, a, (b,))
+        Base.:(*)(a::$VecT, b::$VecT) = a .* b
+        Base.:(*)(a::Number, b::$VecT) = a .* b
+        Base.:(*)(a::$VecT, b::Number) = a .* b
+
+        Base.:(+)(a::$VecT, b::$VecT) = a .+ b
+        Base.:(+)(a::Number, b::$VecT) = a .+ b
+        Base.:(+)(a::$VecT, b::Number) = a .+ b
+
+        Base.:(-)(a::$VecT) = (-).(a)
+        Base.:(-)(a::$VecT, b::$VecT) = a .- b
+        Base.:(-)(a::Number, b::$VecT) = a .- b
+        Base.:(-)(a::$VecT, b::Number) = a .- b
+
+        function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{$(VecT){N,T}}) where {N,T}
+            $(VecT){N,T}(ntuple(i-> rand(rng, T), N))
+        end
+        function Random.randn(rng::Random.AbstractRNG, ::Type{$(VecT){N,T}}) where {N,T}
+            $(VecT){N,T}(ntuple(i-> randn(rng, T), N))
+        end
+        function LinearAlgebra.cross(a::$(VecT){3}, b::$(VecT){3})
+            @inbounds elements = (a[2]*b[3]-a[3]*b[2],
+                                    a[3]*b[1]-a[1]*b[3],
+                                    a[1]*b[2]-a[2]*b[1])
+            return $(VecT)(elements)
         end
     end
     return esc(expr)
 end
 
+LinearAlgebra.cross(a::StaticVector{2}, b::StaticVector{2}) = a[1]*b[2]-a[2]*b[1]
+LinearAlgebra.norm(a::StaticVector) = sqrt(dot(a,a))
+LinearAlgebra.normalize(a::StaticVector) = a ./ norm(a)
+
+Base.eltype(::StaticVector{N, T}) where {N, T} = T
+Base.eltype(::Type{<: StaticVector{N, T}}) where {N, T} = T
+
+Base.size(::StaticVector{N}) where {N} = (N,)
+Base.size(::Type{<: StaticVector{N}}) where {N} = (N,)
+Base.length(::StaticVector{N}) where {N} = N
+Base.length(::Type{<: StaticVector{N}}) where {N} = N
+Base.ndims(::Type{<: StaticVector}) = 1
+
+function Base.iterate(A::StaticVector, i=1)
+    i - 1 < length(A) ? (A[i], i + 1) : nothing
+end
+
+function unit(::Type{T}, i::Integer) where {T <: StaticVector}
+    tup = ntuple(Val(length(T))) do j
+        return ifelse(i == j, 1, 0)
+    end
+    return T(tup)
+end
+
+Base.zero(::Type{V}) where {V <:StaticVector} = V(0)
+
+function Base.:(==)(a::StaticVector{N}, b::StaticVector{N}) where N
+    for i in 1:N
+        a[i] == b[i] || return false
+    end
+    return true
+end
+
+function Base.isapprox(a::StaticVector{N}, b::StaticVector{N}; kw...) where N
+    for i in 1:N
+        isapprox(a[i], b[i]; kw...) || return false
+    end
+    return true
+end
+
+@generated function Base.transpose(b::StaticVector{N,T}) where {N,T}
+    expr = [:(transpose(b[$i])) for i=1:N]
+    return quote
+        Mat{1,N,T}($(expr...))
+    end
+end
+Base.reverse(x::P) where P <: StaticVector = P(reverse(x.data))
+
+# Since we don't inherit from AbstractArray, some extra functions need to be overloaded
+LinearAlgebra.promote_leaf_eltypes(x::StaticVector{N, T}) where {N,T} = T
+
+
+
 abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
+
 @fixed_vector Point AbstractPoint
 @fixed_vector Vec StaticVector
 
-const Mat = SMatrix
-const VecTypes{N,T} = Union{StaticVector{N,T},NTuple{N,T}}
-const Vecf{N} = Vec{N,Float32}
+function Base.getindex(mat::Mat{R, C, T}, r::Vec{NR}, c::Vec{NC}) where {R, C, NR, NC, T}
+    idx = CartesianIndices((NR, NC))
+    data = ntuple(NR * NC) do i
+        ri, ci = Tuple(idx[i])
+        return mat[r[ri], c[ci]]
+    end
+    return Mat{NR, NC, T}(data)
+end
+
+Base.lastindex(::StaticVector{N}) where N = N
+
+Base.broadcasted(f, a::Point, b::GeometryBasics.Vec) = Vec(f.(a.data, b.data))
+Base.broadcasted(f, a::Vec, b::Point) = Vec(f.(a.data, b.data))
+
+Base.:(+)(a::Vec{N}, b::Point{N}) where {N} = Point{N}(a.data .+ b.data)
+
+const VecTypes{N,T} = Union{StaticVector{N,T}, NTuple{N,T}}
+const Vecf{N} = Vec{N, Float32}
 const Pointf{N} = Point{N,Float32}
 Base.isnan(p::Union{AbstractPoint,Vec}) = any(x -> isnan(x), p)
 
@@ -138,8 +248,9 @@ for i in 1:4
     end
     name = Symbol("Mat$i")
     namef = Symbol("Mat$(i)f")
+    namef = Symbol("Mat$(i)f")
     @eval begin
-        const $name{T} = $Mat{$i,$i,T,$(i * i)}
+        const $name{T} = $Mat{$i,$i,T, $(i*i)}
         const $namef = $name{Float32}
         export $name
         export $namef
diff --git a/src/lines.jl b/src/lines.jl
index c89ae8da..ae891fb8 100644
--- a/src/lines.jl
+++ b/src/lines.jl
@@ -9,7 +9,7 @@ function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
     T = promote_type(T1, T2)
     v1, v2 = a
     v3, v4 = b
-    MT = Mat{2,2,T,4}
+    MT = Mat{2,2,T}
     p0 = zero(Point2{T})
 
     verticalA = v1[1] == v2[1]
diff --git a/src/mat.jl b/src/mat.jl
new file mode 100644
index 00000000..0838de8f
--- /dev/null
+++ b/src/mat.jl
@@ -0,0 +1,164 @@
+import LinearAlgebra: inv, det
+import Random
+
+struct Mat{Row, Column, T, L} <: AbstractMatrix{T}
+    values::NTuple{L, T}
+    function Mat{R, C, T}(values::NTuple{L, T}) where {R, C, T, L}
+        @assert L == R * C "$R * $C needs to be $L"
+        return new{R, C, T, L}(values)
+    end
+        function Mat{R, C, T, L}(values::NTuple{L, T}) where {R, C, T, L}
+        @assert L == R * C "$R * $C needs to be $L"
+        return new{R, C, T, L}(values)
+    end
+end
+
+Base.size(::Mat{R, C}) where {R, C} = (R, C)
+Base.size(::Type{<: Mat{R, C}}) where {R, C} = (R, C)
+Base.ndims(::Type{<: Mat}) = 2
+Base.getindex(mat::Mat{R, C}, i) where {R, C} = mat.values[i]
+
+
+Base.IndexStyle(::Mat)= Base.IndexLinear()
+
+function Mat{C, R, T}(::LinearAlgebra.UniformScaling) where {C, R, T}
+    idx = CartesianIndices((R, C))
+    data = ntuple(C * R) do i
+        ci, ri = Tuple(idx[i])
+        return ci === ri ? T(1) : T(0)
+    end
+    return Mat{R, C, T}(data)
+end
+
+function Mat{C, R, T, L}(::LinearAlgebra.UniformScaling) where {C, R, T, L}
+    idx = CartesianIndices((R, C))
+    data = ntuple(C * R) do i
+        ci, ri = Tuple(idx[i])
+        return ci === ri ? T(1) : T(0)
+    end
+    return Mat{R, C, T, L}(data)
+end
+
+Mat{C, R, T}(args...) where {C, R, T} = Mat{C, R, T}(args)
+Mat{C, R, T, L}(args...) where {C, R, T, L} = Mat{C, R, T}(args)
+Mat{C}(args...) where {C} = Mat{C, C}(args)
+Mat{C}(arg) where {C} = Mat{C, C}(arg)
+Mat{C, R}(x::Tuple) where {C, R} = Mat{C, R}(promote(x...))
+Mat{C, R, T}(x::Tuple) where {C, R, T} = Mat{C, R, T}(convert(NTuple{length(x), T}, x))
+Mat{C, R}(x::NTuple{L, T}) where {C, R, L, T} = Mat{C, R, T}(x)
+Mat{C, R, T1}(x::NTuple{L, T2}) where {C, R, L, T1, T2} = Mat{C, R, T1}(convert(NTuple{L, T1}, x))
+
+Mat{C, R}(x::AbstractMatrix{T}) where {C, R, T} = Mat{C, R, T}(x)
+Mat{C, R, T}(x::AbstractMatrix) where {C, R, T} = Mat{C, R, T}(ntuple(i-> convert(T, x[i]), C*R))
+
+Base.convert(::Type{Mat{C, R, T, L}}, from::Mat{C, R}) where {C, R, T, L} = Mat{C, R, T}(from.values)
+
+# Matrix products
+# General shape mismatched versions are errors
+(*)(a::Mat{M, N, T1}, b::Mat{O, K, T2}) where {T1, T2, M, N, O, K} = throw(DimensionMismatch("$N != $O in $(typeof(a)) and $(typeof(b))"))
+
+# matrix * matrix
+@generated function *(a::Mat{M, N, T1}, b::Mat{N, P, T2}) where {T1, T2, M, N, P}
+    elements = Expr(:tuple)
+    for j in 1:P
+        for i in 1:M
+            plus_expr = Expr(:call, :+, (:(a[$i,$k]*b[$k,$j]) for k = 1:N)...)
+            push!(elements.args, plus_expr)
+        end
+    end
+    :(Mat{$M, $P}($elements))
+end
+
+@inline det(A::Mat{1, 1, T}) where {T} = @inbounds return ( A[1] )
+@inline det(A::Mat{2, 2, T}) where {T} = @inbounds return ( A[1,1]*A[2,2] - A[1,2]*A[2,1])
+@inline det(A::Mat{3, 3, T}) where {T} = @inbounds return (
+    A[1,1]*(A[2,2]*A[3,3]-A[2,3]*A[3,2]) -
+    A[1,2]*(A[2,1]*A[3,3]-A[2,3]*A[3,1]) +
+    A[1,3]*(A[2,1]*A[3,2]-A[2,2]*A[3,1])
+)
+
+det(A::Mat{4, 4, T}) where {T} = @inbounds return (
+    A[13] * A[10]  * A[7]  * A[4]  - A[9] * A[14] * A[7]  * A[4]   -
+    A[13] * A[6]   * A[11] * A[4]  + A[5] * A[14] * A[11] * A[4]   +
+    A[9]  * A[6]   * A[15] * A[4]  - A[5] * A[10] * A[15] * A[4]   -
+    A[13] * A[10]  * A[3]  * A[8]  + A[9] * A[14] * A[3]  * A[8]   +
+    A[13] * A[2]   * A[11] * A[8]  - A[1] * A[14] * A[11] * A[8]   -
+    A[9]  * A[2]   * A[15] * A[8]  + A[1] * A[10] * A[15] * A[8]   +
+    A[13] * A[6]   * A[3]  * A[12] - A[5] * A[14] * A[3]  * A[12]  -
+    A[13] * A[2]   * A[7]  * A[12] + A[1] * A[14] * A[7]  * A[12]  +
+    A[5]  * A[2]   * A[15] * A[12] - A[1] * A[6]  * A[15] * A[12]  -
+    A[9]  * A[6]   * A[3]  * A[16] + A[5] * A[10] * A[3]  * A[16]  +
+    A[9]  * A[2]   * A[7]  * A[16] - A[1] * A[10] * A[7]  * A[16]  -
+    A[5]  * A[2]   * A[11] * A[16] + A[1] * A[6]  * A[11] * A[16]
+)
+
+det(A::Mat) = det(Matrix(A))
+
+inv(A::Mat{1, 1, T, 1}) where T = @inbounds return Mat{1, 1, T, 1}(inv(A[1]))
+function inv(A::Mat{2, 2, T, L}) where {T, L}
+    determinant = det(A)
+    @inbounds return Mat{2, 2, T}(
+        (A[2,2] /determinant, -A[2,1]/determinant),
+        (-A[1,2]/determinant, A[1,1] /determinant)
+    )
+end
+function inv(A::Mat{3, 3, T, L}) where {T, L}
+    determinant = det(A)
+    @inbounds return  Mat{3, 3, T}(
+        (A[2,2]*A[3,3]-A[2,3]*A[3,2]) /determinant,
+        -(A[2,1]*A[3,3]-A[2,3]*A[3,1])/determinant,
+        (A[2,1]*A[3,2]-A[2,2]*A[3,1]) /determinant,
+
+        -(A[1,2]*A[3,3]-A[1,3]*A[3,2])/determinant,
+        (A[1,1]*A[3,3]-A[1,3]*A[3,1]) /determinant,
+        -(A[1,1]*A[3,2]-A[1,2]*A[3,1])/determinant,
+
+        (A[1,2]*A[2,3]-A[1,3]*A[2,2]) /determinant,
+        -(A[1,1]*A[2,3]-A[1,3]*A[2,1])/determinant,
+        (A[1,1]*A[2,2]-A[1,2]*A[2,1]) /determinant
+    )
+end
+
+function inv(A::Mat{4, 4, T, L}) where {T, L}
+    determinant = det(A)
+    @inbounds return Mat{4, 4, T}(
+        (A[2,3]*A[3,4]*A[4,2] - A[2,4]*A[3,3]*A[4,2] + A[2,4]*A[3,2]*A[4,3] - A[2,2]*A[3,4]*A[4,3] - A[2,3]*A[3,2]*A[4,4] + A[2,2]*A[3,3]*A[4,4]) / determinant,
+        (A[2,4]*A[3,3]*A[4,1] - A[2,3]*A[3,4]*A[4,1] - A[2,4]*A[3,1]*A[4,3] + A[2,1]*A[3,4]*A[4,3] + A[2,3]*A[3,1]*A[4,4] - A[2,1]*A[3,3]*A[4,4]) / determinant,
+        (A[2,2]*A[3,4]*A[4,1] - A[2,4]*A[3,2]*A[4,1] + A[2,4]*A[3,1]*A[4,2] - A[2,1]*A[3,4]*A[4,2] - A[2,2]*A[3,1]*A[4,4] + A[2,1]*A[3,2]*A[4,4]) / determinant,
+        (A[2,3]*A[3,2]*A[4,1] - A[2,2]*A[3,3]*A[4,1] - A[2,3]*A[3,1]*A[4,2] + A[2,1]*A[3,3]*A[4,2] + A[2,2]*A[3,1]*A[4,3] - A[2,1]*A[3,2]*A[4,3]) / determinant,
+
+        (A[1,4]*A[3,3]*A[4,2] - A[1,3]*A[3,4]*A[4,2] - A[1,4]*A[3,2]*A[4,3] + A[1,2]*A[3,4]*A[4,3] + A[1,3]*A[3,2]*A[4,4] - A[1,2]*A[3,3]*A[4,4]) / determinant,
+        (A[1,3]*A[3,4]*A[4,1] - A[1,4]*A[3,3]*A[4,1] + A[1,4]*A[3,1]*A[4,3] - A[1,1]*A[3,4]*A[4,3] - A[1,3]*A[3,1]*A[4,4] + A[1,1]*A[3,3]*A[4,4]) / determinant,
+        (A[1,4]*A[3,2]*A[4,1] - A[1,2]*A[3,4]*A[4,1] - A[1,4]*A[3,1]*A[4,2] + A[1,1]*A[3,4]*A[4,2] + A[1,2]*A[3,1]*A[4,4] - A[1,1]*A[3,2]*A[4,4]) / determinant,
+        (A[1,2]*A[3,3]*A[4,1] - A[1,3]*A[3,2]*A[4,1] + A[1,3]*A[3,1]*A[4,2] - A[1,1]*A[3,3]*A[4,2] - A[1,2]*A[3,1]*A[4,3] + A[1,1]*A[3,2]*A[4,3]) / determinant,
+
+        (A[1,3]*A[2,4]*A[4,2] - A[1,4]*A[2,3]*A[4,2] + A[1,4]*A[2,2]*A[4,3] - A[1,2]*A[2,4]*A[4,3] - A[1,3]*A[2,2]*A[4,4] + A[1,2]*A[2,3]*A[4,4]) / determinant,
+        (A[1,4]*A[2,3]*A[4,1] - A[1,3]*A[2,4]*A[4,1] - A[1,4]*A[2,1]*A[4,3] + A[1,1]*A[2,4]*A[4,3] + A[1,3]*A[2,1]*A[4,4] - A[1,1]*A[2,3]*A[4,4]) / determinant,
+        (A[1,2]*A[2,4]*A[4,1] - A[1,4]*A[2,2]*A[4,1] + A[1,4]*A[2,1]*A[4,2] - A[1,1]*A[2,4]*A[4,2] - A[1,2]*A[2,1]*A[4,4] + A[1,1]*A[2,2]*A[4,4]) / determinant,
+        (A[1,3]*A[2,2]*A[4,1] - A[1,2]*A[2,3]*A[4,1] - A[1,3]*A[2,1]*A[4,2] + A[1,1]*A[2,3]*A[4,2] + A[1,2]*A[2,1]*A[4,3] - A[1,1]*A[2,2]*A[4,3]) / determinant,
+
+        (A[1,4]*A[2,3]*A[3,2] - A[1,3]*A[2,4]*A[3,2] - A[1,4]*A[2,2]*A[3,3] + A[1,2]*A[2,4]*A[3,3] + A[1,3]*A[2,2]*A[3,4] - A[1,2]*A[2,3]*A[3,4]) / determinant,
+        (A[1,3]*A[2,4]*A[3,1] - A[1,4]*A[2,3]*A[3,1] + A[1,4]*A[2,1]*A[3,3] - A[1,1]*A[2,4]*A[3,3] - A[1,3]*A[2,1]*A[3,4] + A[1,1]*A[2,3]*A[3,4]) / determinant,
+        (A[1,4]*A[2,2]*A[3,1] - A[1,2]*A[2,4]*A[3,1] - A[1,4]*A[2,1]*A[3,2] + A[1,1]*A[2,4]*A[3,2] + A[1,2]*A[2,1]*A[3,4] - A[1,1]*A[2,2]*A[3,4]) / determinant,
+        (A[1,2]*A[2,3]*A[3,1] - A[1,3]*A[2,2]*A[3,1] + A[1,3]*A[2,1]*A[3,2] - A[1,1]*A[2,3]*A[3,2] - A[1,2]*A[2,1]*A[3,3] + A[1,1]*A[2,2]*A[3,3]) / determinant
+    )
+end
+
+inv(A::Mat) = typeof(A)(inv(Matrix(A)))
+
+function Random.rand(rng::Random.AbstractRNG, ::Random.SamplerType{<: Mat{R, C, T}}) where {R,C,T}
+    return Mat{R, C, T}(ntuple(i-> rand(rng, T), R*C))
+end
+
+function Random.randn(rng::Random.AbstractRNG, ::Type{<: Mat{R, C, T}}) where {R,C,T}
+    return Mat{R, C, T}(ntuple(i-> randn(rng, T), R*C))
+end
+
+function Base.transpose(a::Mat{R, C, T}) where {R, C, T}
+    idx = CartesianIndices((R, C))
+    data = ntuple(R * C) do i
+        ci, ri = Tuple(idx[i])
+        return a[ri, ci]
+    end
+    Mat{R, C, T}(data)
+end
diff --git a/src/meshes.jl b/src/meshes.jl
index d7d41ffa..9c0bbed5 100644
--- a/src/meshes.jl
+++ b/src/meshes.jl
@@ -1,4 +1,7 @@
 const FaceMesh{Dim,T,Element} = Mesh{Dim,T,Element,<:FaceView{Element}}
+const GLTriangleElement = Triangle{3,Float32}
+const GLTriangleFace = TriangleFace{GLIndex}
+
 
 coordinates(mesh::FaceMesh) = coordinates(getfield(mesh, :simplices))
 faces(mesh::FaceMesh) = faces(getfield(mesh, :simplices))
@@ -14,78 +17,6 @@ function normals(mesh::AbstractMesh)
     return nothing
 end
 
-const GLTriangleElement = Triangle{3,Float32}
-const GLTriangleFace = TriangleFace{GLIndex}
-const PointWithUV{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:uv,),Tuple{Vec{2,T}}}
-const PointWithNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals,),Tuple{Vec{3,T}}}
-const PointWithUVNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals, :uv),
-                                           Tuple{Vec{3,T},Vec{2,T}}}
-const PointWithUVWNormal{Dim,T} = PointMeta{Dim,T,Point{Dim,T},(:normals, :uvw),
-                                            Tuple{Vec{3,T},Vec{3,T}}}
-
-"""
-    TriangleMesh{Dim, T, PointType}
-
-Abstract Mesh with triangle elements of eltype `T`.
-"""
-const TriangleMesh{Dim,T,PointType} = AbstractMesh{TriangleP{Dim,T,PointType}}
-
-"""
-    PlainMesh{Dim, T}
-
-Triangle mesh with no meta information (just points + triangle faces)
-"""
-const PlainMesh{Dim,T} = TriangleMesh{Dim,T,Point{Dim,T}}
-const GLPlainMesh{Dim} = PlainMesh{Dim,Float32}
-const GLPlainMesh2D = GLPlainMesh{2}
-const GLPlainMesh3D = GLPlainMesh{3}
-
-"""
-    UVMesh{Dim, T}
-
-PlainMesh with texture coordinates meta at each point.
-`uvmesh.uv isa AbstractVector{Vec2f}`
-"""
-const UVMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUV{Dim,T}}
-const GLUVMesh{Dim} = UVMesh{Dim,Float32}
-const GLUVMesh2D = UVMesh{2}
-const GLUVMesh3D = UVMesh{3}
-
-"""
-    NormalMesh{Dim, T}
-
-PlainMesh with normals meta at each point.
-`normalmesh.normals isa AbstractVector{Vec3f}`
-"""
-const NormalMesh{Dim,T} = TriangleMesh{Dim,T,PointWithNormal{Dim,T}}
-const GLNormalMesh{Dim} = NormalMesh{Dim,Float32}
-const GLNormalMesh2D = GLNormalMesh{2}
-const GLNormalMesh3D = GLNormalMesh{3}
-
-"""
-    NormalUVMesh{Dim, T}
-
-PlainMesh with normals and uv meta at each point.
-`normalmesh.normals isa AbstractVector{Vec3f}`
-`normalmesh.uv isa AbstractVector{Vec2f}`
-"""
-const NormalUVMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUVNormal{Dim,T}}
-const GLNormalUVMesh{Dim} = NormalUVMesh{Dim,Float32}
-const GLNormalUVMesh2D = GLNormalUVMesh{2}
-const GLNormalUVMesh3D = GLNormalUVMesh{3}
-
-"""
-    NormalUVWMesh{Dim, T}
-
-PlainMesh with normals and uvw (texture coordinates in 3D) meta at each point.
-`normalmesh.normals isa AbstractVector{Vec3f}`
-`normalmesh.uvw isa AbstractVector{Vec3f}`
-"""
-const NormalUVWMesh{Dim,T} = TriangleMesh{Dim,T,PointWithUVWNormal{Dim,T}}
-const GLNormalUVWMesh{Dim} = NormalUVWMesh{Dim,Float32}
-const GLNormalUVWMesh2D = GLNormalUVWMesh{2}
-const GLNormalUVWMesh3D = GLNormalUVWMesh{3}
-
 function decompose_triangulate_fallback(primitive::Meshable; pointtype, facetype)
     positions = decompose(pointtype, primitive)
     faces = decompose(facetype, primitive)
diff --git a/src/offsetintegers.jl b/src/offsetintegers.jl
index 12c567be..c0043c53 100644
--- a/src/offsetintegers.jl
+++ b/src/offsetintegers.jl
@@ -20,7 +20,7 @@ raw(x::Integer) = x
 value(x::OffsetInteger{O,T}) where {O,T} = raw(x) - O
 value(x::Integer) = x
 
-function show(io::IO, oi::OffsetInteger)
+function Base.show(io::IO, oi::OffsetInteger)
     return print(io, "|$(raw(oi)) (indexes as $(value(oi))|")
 end
 
diff --git a/src/precompile.jl b/src/precompile.jl
deleted file mode 100644
index f7b1c52b..00000000
--- a/src/precompile.jl
+++ /dev/null
@@ -1,39 +0,0 @@
-macro warnpcfail(ex::Expr)
-    modl = __module__
-    file = __source__.file === nothing ? "?" : String(__source__.file)
-    line = __source__.line
-    quote
-        $(esc(ex)) || @warn """precompile directive
-     $($(Expr(:quote, ex)))
- failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
-    end
-end
-
-function _precompile_()
-    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
-    @warnpcfail precompile(HyperRectangle{2,Float32}, (Int, Int, Int, Int))
-    @warnpcfail precompile(==, (HyperRectangle{2,Float32}, HyperRectangle{2,Float32}))
-    @warnpcfail precompile(normal_mesh, (Tesselation{3,Float32,Cylinder{3,Float32},1},))
-    @warnpcfail precompile(normal_mesh, (Tesselation{3,Float32,HyperSphere{3,Float32},1},))
-    @warnpcfail precompile(normal_mesh, (HyperSphere{3,Float32},))
-    @warnpcfail precompile(Tuple{typeof(*),SMatrix{4, 4, Float32, 16},HyperRectangle{3, Float32}})   # time: 0.11091917
-    @warnpcfail precompile(Tuple{typeof(coordinates),HyperRectangle{2, Float32},Tuple{Int64, Int64}})   # time: 0.08693867
-    @warnpcfail precompile(union, (HyperRectangle{3, Float32}, HyperRectangle{3, Float32}))
-    @warnpcfail precompile(Tuple{typeof(decompose),Type{Point{2, Float32}},HyperRectangle{2, Float32}})   # time: 0.026609203
-    @warnpcfail precompile(Tuple{Type{HyperRectangle{3, Float32}},HyperRectangle{2, Float32}})   # time: 0.023717888
-    @warnpcfail precompile(Tuple{typeof(+),HyperRectangle{3, Float32},Point{3, Float32}})   # time: 0.006633118
-    @warnpcfail precompile(Tuple{Type{Rect2{T} where T},Float32,Float32,Float32,Float32})   # time: 0.001636267
-    @warnpcfail precompile(Tuple{typeof(*),HyperRectangle{2, Float32},Float32})   # time: 0.001057589
-
-    if Base.VERSION >= v"1.6.0-DEV.1083"
-        @warnpcfail precompile(triangle_mesh, (Polygon{2, Float32, Point2f, LineString{2, Float32, Point2f,
-                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}, false}},
-                                        Vector{LineString{2, Float32, Point2f, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}, false}}}},))
-    else
-        @warnpcfail precompile(triangle_mesh, (Polygon{2, Float32, Point2f, LineString{2, Float32, Point2f,
-                                        Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}}},
-                                        Vector{LineString{2, Float32, Point2f, Base.ReinterpretArray{Line{2, Float32}, 1, Tuple{Point2f, Point2f}, TupleView{Tuple{Point2f, Point2f}, 2, 1, Vector{Point2f}}}}}},))
-    end
-
-    @warnpcfail precompile(split_intersections, (Vector{Point2f},))
-end
diff --git a/src/precompiles.jl b/src/precompiles.jl
new file mode 100644
index 00000000..df41b1e6
--- /dev/null
+++ b/src/precompiles.jl
@@ -0,0 +1,6 @@
+
+function _precompile_()
+    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
+    Point2f(0.5, 0.1) in Triangle(Point2f(0), Point2f(0.5, 1), Point2f(1, 0))
+    decompose(GLTriangleFace, [Point2f(0), Point2f(0.5, 1), Point2f(1, 0)])
+end
diff --git a/src/primitives/cylinders.jl b/src/primitives/cylinders.jl
index 9804cf53..e9467ad2 100644
--- a/src/primitives/cylinders.jl
+++ b/src/primitives/cylinders.jl
@@ -28,24 +28,24 @@ direction(c::Cylinder{N,T}) where {N,T} = (c.extremity .- c.origin) ./ height(c)
 
 function rotation(c::Cylinder{2,T}) where {T}
     d2 = direction(c)
-    u = @SVector [d2[1], d2[2], T(0)]
-    v = @MVector [u[2], -u[1], T(0)]
-    normalize!(v)
-    return hcat(v, u, @SVector T[0, 0, 1])
+    u = Vec{3, T}(d2[1], d2[2], T(0))
+    v = Vec{3, T}(u[2], -u[1], T(0))
+    v = normalize(v)
+    return Mat{3, 3, T}(v..., u..., 0, 0, 1)
 end
 
 function rotation(c::Cylinder{3,T}) where {T}
     d3 = direction(c)
-    u = @SVector [d3[1], d3[2], d3[3]]
+    u = Vec{3, T}(d3[1], d3[2], d3[3])
     if abs(u[1]) > 0 || abs(u[2]) > 0
-        v = @MVector [u[2], -u[1], T(0)]
+        v = Vec{3, T}(u[2], -u[1], T(0))
     else
-        v = @MVector [T(0), -u[3], u[2]]
+        v = Vec{3, T}(T(0), -u[3], u[2])
     end
-    normalize!(v)
-    w = @SVector [u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1],
-                  u[1] * v[2] - u[2] * v[1]]
-    return hcat(v, w, u)
+    v = normalize(v)
+    w = Vec{3, T}(u[2] * v[3] - u[3] * v[2], -u[1] * v[3] + u[3] * v[1],
+                  u[1] * v[2] - u[2] * v[1])
+    return Mat{3, 3, T}(v..., w..., u...)
 end
 
 function coordinates(c::Cylinder{2,T}, nvertices=(2, 2)) where {T}
diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index d8d06348..0ac3a759 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -177,8 +177,8 @@ split(b::Rect, axis, value::Number) = _split(b, axis, value)
 function _split(b::H, axis, value) where {H<:Rect}
     bmin = minimum(b)
     bmax = maximum(b)
-    b1max = setindex(bmax, value, axis)
-    b2min = setindex(bmin, value, axis)
+    b1max = Base.setindex(bmax, value, axis)
+    b2min = Base.setindex(bmin, value, axis)
 
     return H(bmin, b1max - bmin), H(b2min, bmax - b2min)
 end
diff --git a/src/triangulation.jl b/src/triangulation.jl
index 94a10eb6..7bc8f931 100644
--- a/src/triangulation.jl
+++ b/src/triangulation.jl
@@ -68,23 +68,22 @@ Determine if a point is inside of a triangle.
 """
 function Base.in(P::T, triangle::Triangle) where {T<:AbstractPoint}
     A, B, C = coordinates(triangle)
-    a = C .- B
-    b = A .- C
-    c = B .- A
 
-    ap = P .- A
-    bp = P .- B
-    cp = P .- C
+    ap = A .- P
+    bp = B .- P
+    cp = C .- P
 
-    a_bp = a[1] * bp[2] - a[2] * bp[1]
-    c_ap = c[1] * ap[2] - c[2] * ap[1]
-    b_cp = b[1] * cp[2] - b[2] * cp[1]
+    u = cross(bp, cp)
+    v = cross(cp, ap)
+    w = cross(ap, bp)
 
     t0 = zero(eltype(T))
-
-    return ((a_bp >= t0) && (b_cp >= t0) && (c_ap >= t0))
+    dot(u, v) < t0 && return false
+    dot(u, w) < t0 && return false
+    return true
 end
 
+
 function snip(contour::AbstractVector{<:AbstractPoint}, u, v, w, n, V)
     A = contour[V[u]]
     B = contour[V[v]]
diff --git a/src/viewtypes.jl b/src/viewtypes.jl
index b700b6ea..3fbd30bb 100644
--- a/src/viewtypes.jl
+++ b/src/viewtypes.jl
@@ -96,7 +96,7 @@ end
         columns = ntuple(N) do i
             return view(points, ((i - 1) * seglen + 1):(i * seglen))
         end
-        return StructArray{Point{N,T}}(columns)
+        return P.(columns...)
     else
         error("Dim 1 or 2 must be equal to the point dimension!")
     end
diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl
index 7390ec14..fccf527b 100644
--- a/test/geometrytypes.jl
+++ b/test/geometrytypes.jl
@@ -77,7 +77,6 @@ using Test, GeometryBasics
         @test GeometryBasics.faces(m) == faces
         @test GeometryBasics.coordinates(m) ≈ positions
         m = normal_mesh(s)# just test that it works without explicit resolution parameter
-        @test m isa GLNormalMesh
 
         muv = uv_mesh(s)
         @test Rect(Point.(texturecoordinates(muv))) == Rect2f(Vec2f(0), Vec2f(1.0))
diff --git a/test/runtests.jl b/test/runtests.jl
index 956ebf66..19bdf7cd 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,10 +1,9 @@
-using Test, Random, StructArrays, Tables, StaticArrays, OffsetArrays
+using Test, Random, StructArrays, OffsetArrays
 using GeometryBasics
 using LinearAlgebra
 using GeometryBasics: attributes
 
 @testset "GeometryBasics" begin
-
 @testset "algorithms" begin
     cube = Rect(Vec3f(-0.5), Vec3f(1))
     cube_faces = decompose(TriangleFace{Int}, faces(cube))
@@ -39,7 +38,7 @@ end
         @testset "per vertex attributes" begin
             points = rand(Point{3, Float64}, 8)
             tfaces = TetrahedronFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8)]
-            normals = rand(SVector{3, Float64}, 8)
+            normals = rand(Vec{3, Float64}, 8)
             stress = LinRange(0, 1, 8)
             mesh = Mesh(meta(points, normals = normals, stress = stress), tfaces)
 
@@ -111,7 +110,6 @@ end
     end
     @testset "point with metadata" begin
         p = Point(1.1, 2.2)
-        @test p isa AbstractVector{Float64}
         pm = PointMeta(1.1, 2.2; a=1, b=2)
         p1 = Point(2.2, 3.6)
         p2 = [p, p1]
@@ -190,7 +188,6 @@ end
 
     @testset "MetaT{Point}" begin
         p = Point(1.1, 2.2)
-        @test p isa AbstractVector{Float64}
         pm = MetaT(Point(1.1, 2.2); a=1, b=2)
         p1 = Point(2.2, 3.6)
         p2 = [p, p1]
@@ -245,7 +242,7 @@ end
         @testset "per vertex attributes" begin
             points = rand(Point{3, Float64}, 8)
             tfaces = TetrahedronFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8)]
-            normals = rand(SVector{3, Float64}, 8)
+            normals = rand(Vec{3, Float64}, 8)
             stress = LinRange(0, 1, 8)
             mesh_nometa = Mesh(points, tfaces)
             mesh = MetaT(mesh_nometa, normals = normals, stress = stress)
@@ -410,10 +407,6 @@ end
         meshuv = Mesh(meta(points; uv=uv), tfaces)
         meshuvnormal = Mesh(meta(points; normals=normals, uv=uv), tfaces)
 
-        @test mesh isa GLPlainMesh
-        @test meshuv isa GLUVMesh3D
-        @test meshuvnormal isa GLNormalUVMesh3D
-
         t = Tesselation(Rect2f(0, 0, 2, 2), (30, 30))
         m = GeometryBasics.mesh(t, pointtype=Point3f, facetype=QuadFace{Int})
         m2 = GeometryBasics.mesh(m, facetype=QuadFace{GLIndex})
@@ -472,22 +465,22 @@ end
 
     primitive = Sphere(Point3f(0), 1)
     m_normal = normal_mesh(primitive)
-    @test normals(m_normal) isa Vector{Vec3f}
+    @test GeometryBasics.normals(m_normal) isa Vector{Vec3f}
     primitive = Rect2(0, 0, 1, 1)
     m_normal = normal_mesh(primitive)
-    @test normals(m_normal) isa Vector{Vec3f}
+    @test GeometryBasics.normals(m_normal) isa Vector{Vec3f}
     primitive = Rect3(0, 0, 0, 1, 1, 1)
     m_normal = normal_mesh(primitive)
-    @test normals(m_normal) isa Vector{Vec3f}
+    @test GeometryBasics.normals(m_normal) isa Vector{Vec3f}
 
     points = decompose(Point2f, Circle(Point2f(0), 1))
     tmesh = triangle_mesh(points)
-    @test normals(tmesh) == nothing
+    @test GeometryBasics.normals(tmesh) == nothing
 
     m = GeometryBasics.mesh(Sphere(Point3f(0), 1))
-    @test normals(m) == nothing
+    @test GeometryBasics.normals(m) == nothing
     m_normals = pointmeta(m, Normal())
-    @test normals(m_normals) isa Vector{Vec3f}
+    @test GeometryBasics.normals(m_normals) isa Vector{Vec3f}
 
     @test texturecoordinates(m) == nothing
     r2 = Rect2(0.0, 0.0, 1.0, 1.0)
@@ -596,13 +589,11 @@ end
     @test coordinates(m) === decompose(Point{3, Float64}, m)
 
     tmesh = triangle_mesh(m)
-    @test tmesh isa GLPlainMesh
     @test coordinates(tmesh) === decompose(Point3f, tmesh)
 
     nmesh = normal_mesh(m)
-    @test nmesh isa GLNormalMesh
     @test metafree(coordinates(nmesh)) === decompose(Point3f, nmesh)
-    @test normals(nmesh) === decompose_normals(nmesh)
+    @test GeometryBasics.normals(nmesh) === decompose_normals(nmesh)
 
     m = GeometryBasics.mesh(s, pointtype=Point3f)
     @test m isa Mesh{3, Float32}
@@ -715,4 +706,4 @@ end
     include("fixed_arrays.jl")
 end
 
-end  # testset "GeometryBasics"
+end