diff --git a/REQUIRE b/REQUIRE
index 6771323..841e852 100644
--- a/REQUIRE
+++ b/REQUIRE
@@ -7,5 +7,6 @@ Formatting
 Logging
 TimerOutputs
 AbaqusReader
-FEMQuad
 AsterReader
+FEMBasis
+FEMQuad
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index ccd78ed..354e986 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -19,6 +19,12 @@ import Base: getindex, setindex!, convert, length, size, isapprox, similar,
              start, first, next, done, last, endof, vec, ==, +, -, *, /, haskey, copy,
              push!, isempty, empty!, append!, sparse, full, read
 
+using FEMBasis
+using FEMBasis: AbstractBasis
+using FEMQuad
+using AbaqusReader
+using AsterReader
+
 using Logging
 
 Logging.configure(level=INFO)
@@ -45,7 +51,7 @@ export AbstractPoint, Point, IntegrationPoint, IP, Node
 
 ### ELEMENTS ###
 include("elements.jl") # common element routines
-export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
+export Node, Element, update!, get_connectivity, get_basis,
        get_dbasis, inside, get_local_coordinates, get_element_type,
        filter_by_element_type, get_element_id
 
@@ -62,11 +68,6 @@ export Poi1,
        Wedge6,
        Hex8, Hex20, Hex27
 
-include("elements_nurbs.jl")
-export NSeg, NSurf, NSolid, is_nurbs
-
-#include("hierarchical.jl") # P-elements
-
 include("integrate.jl")  # default integration points for elements
 export get_integration_points
 
diff --git a/src/elements.jl b/src/elements.jl
index cb8cabf..80e01c1 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -1,9 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-abstract type AbstractElement end
-
-type Element{E<:AbstractElement}
+type Element{E<:AbstractBasis}
     id :: Int
     connectivity :: Vector{Int}
     integration_points :: Vector{IP}
@@ -17,10 +15,23 @@ Examples
 --------
 julia> element = Element(Tri3, [1, 2, 3])
 """
-function Element{E<:AbstractElement}(::Type{E}, connectivity::Vector{Int})
+function Element{E<:AbstractBasis}(::Type{E}, connectivity::Vector{Int})
     return Element{E}(-1, connectivity, [], Dict(), E())
 end
 
+"""
+    length(element::Element)
+
+Return the number of nodes in element.
+"""
+function length{B}(element::Element{B})
+    return length(B)
+end
+
+function size{B}(element::Element{B})
+    return size(B)
+end
+
 function getindex(element::Element, field_name::AbstractString)
     return element.fields[field_name]
 end
diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl
index d3f3d42..df79789 100644
--- a/src/elements_lagrange.jl
+++ b/src/elements_lagrange.jl
@@ -1,66 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-global const ELEMENT_DESCRIPTIONS = Dict(
-    "Poi1" => "1 node discrete point element",
-    "Seg2" => "2 node linear segment/line element",
-    "Seg3" => "3 node quadratic segment/line element",
-    "Tri3" => "3 node linear triangle element",
-    "Tri6" => "6 node quadratic triangle element",
-    "Tri7" => "7 node quadratic triangle element (has middle node)",
-    "Quad4" => "4 node linear quadrangle element",
-    "Quad8" => "8 node quadratic quadrangle element (Serendip)",
-    "Quad9" => "9 node quadratic quadrangle element",
-    "Tet4" => "4 node linear tetrahedral element",
-    "Tet10" => "10 node quadratic tetrahedral element",
-    "Pyr5" => "5 node linear pyramid element",
-    "Wedge6" => "6 node linear prismatic element (wedge)",
-    "Wedge15" => "15 node quadratic prismatic element (wedge)",
-    "Hex8" => "8 node linear hexahedral element",
-    "Hex20" => "20 node biquadratic hexahedral element",
-    "Hex27" => "27 node quadratic hexahedral element")
+using FEMBasis
 
-global const ELEMENT_SIZES = Dict(
-    "Poi1" => (0, 1),
-    "Seg2" => (1, 2),
-    "Seg3" => (1, 3),
-    "Tri3" => (2, 3),
-    "Tri6" => (2, 6),
-    "Tri7" => (2, 7),
-    "Quad4" => (2, 4),
-    "Quad8" => (2, 8),
-    "Quad9" => (2, 9),
-    "Tet4" => (3, 4),
-    "Tet10" => (3, 10),
-    "Pyr5" => (3,5),
-    "Wedge6" => (3, 6),
-    "Wedge15" => (3, 15),
-    "Hex8" => (3, 8),
-    "Hex20" => (3, 20),
-    "Hex27" => (3, 27))
-
-""" Return description line of element. """
-function description{T}(element::Element{T})
-    element_type = last(split("$T", '.'))
-    return get(ELEMENT_DESCRIPTIONS, element_type, "Unknown element description")
-end
-
-""" Return size of element, i.e. tuple (n, m) where n is dimension of element
-(0, 1, 2, 3) and m is number of nodes. """
-function size{T}(element::Element{T})
-    element_type = last(split("$T", '.'))
-    return ELEMENT_SIZES[element_type]
-end
-
-""" Return length of element, i.e. number of nodes. """
-function length{T}(element::Element{T})
-    return size(element)[end]
-end
-
-### 0d element
-
-
-type Poi1 <: AbstractElement
+#    "Poi1" => (0, 1),
+"1 node discrete point element",
+type Poi1 <: AbstractBasis
 end
 
 function get_basis(element::Element{Poi1}, ip, time)
@@ -83,541 +28,32 @@ function get_integration_points(element::Poi1, order::Int64)
     return [ (1.0, [] ) ]
 end
 
-function get_reference_coordinates(::Type{Poi1})
+function size(::Type{Poi1})
+    return (0, 1)
+end
+
+function length(::Type{Poi1})
+    return 1
+end
+
+function FEMBasis.get_reference_element_coordinates(::Type{Poi1})
     Vector{Float64}[[0.0]]
 end
 
-### 1d elements
-
-type Seg2 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Seg2})
-    Vector{Float64}[
-        [-1.0], # N1
-        [ 1.0]] # N2
-end
-
-function get_interpolation_polynomial(::Type{Seg2}, xi)
-    [1.0 xi[1]]
-end
-
-function get_interpolation_polynomial(::Type{Seg2}, xi, ::Type{Val{:partial_derivatives}})
-    [0.0 1.0]
-end
-
-#
-
-type Seg3 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Seg3})
-    Vector{Float64}[
-        [-1.0], # N1
-        [ 1.0], # N2
-        [ 0.0]] # N3
-end
-
-function get_interpolation_polynomial(::Type{Seg3}, xi)
-    [1.0 xi[1] xi[1]^2]
-end
-
-function get_interpolation_polynomial(::Type{Seg3}, xi, ::Type{Val{:partial_derivatives}})
-    [0.0 1.0 2.0*xi[1]]
-end
-
-### 2d elements
-
-type Tri3 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Tri3})
-    Vector{Float64}[
-        [0.0, 0.0], # N1
-        [1.0, 0.0], # N2
-        [0.0, 1.0]] # N3
-end
-
-function get_interpolation_polynomial(::Type{Tri3}, xi)
-    [
-        1 xi[1] xi[2]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Tri3}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0.0 1.0 0.0
-        0.0 0.0 1.0
-    ]
-end
-
-#
-
-type Tri6 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Tri6})
-    Vector{Float64}[
-        [0.0, 0.0], # N1
-        [1.0, 0.0], # N2
-        [0.0, 1.0], # N3
-        [0.5, 0.0], # N4
-        [0.5, 0.5], # N5
-        [0.0, 0.5]] # N6
-end
-
-function get_interpolation_polynomial(::Type{Tri6}, xi)
-    [
-        1 xi[1] xi[2] xi[1]^2 xi[1]*xi[2] xi[2]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Tri6}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  1  0  2*xi[1]  xi[2]  0
-        0  0  1  0        xi[1]  2*xi[2]
-    ]
-end
-
-#
-
-type Tri7 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Tri7})
-    Vector{Float64}[
-        [0.0, 0.0], # N1
-        [1.0, 0.0], # N2
-        [0.0, 1.0], # N3
-        [0.5, 0.0], # N4
-        [0.5, 0.5], # N5
-        [0.0, 0.5], # N6
-        [1/3, 1/3]] # N7
-end
-
-function get_interpolation_polynomial(::Type{Tri7}, xi)
-    [
-        1 xi[1] xi[2] xi[1]^2 xi[1]*xi[2] xi[2]^2 xi[1]^2*xi[2]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Tri7}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  1  0  2*xi[1]  xi[2]  0        2*xi[1]*xi[2]^2
-        0  0  1  0        xi[1]  2*xi[2]  2*xi[1]^2*xi[2]
-    ]
-end
-
-#
-
-type Quad4 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Quad4})
-    Vector{Float64}[
-        [-1.0, -1.0], # N1
-        [ 1.0, -1.0], # N2
-        [ 1.0,  1.0], # N3
-        [-1.0,  1.0]] # N4
-end
-
-function get_interpolation_polynomial(::Type{Quad4}, xi)
-    [
-        1.0 xi[1] xi[2] xi[1]*xi[2]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Quad4}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0 1 0 xi[2]
-        0 0 1 xi[1]      
-    ]
-end
-
-#
-
-type Quad8 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Quad8})
-    Vector{Float64}[
-        [-1.0, -1.0], # N1
-        [ 1.0, -1.0], # N2
-        [ 1.0,  1.0], # N3
-        [-1.0,  1.0], # N4
-        [ 0.0, -1.0], # N5
-        [ 1.0,  0.0], # N6
-        [ 0.0,  1.0], # N7
-        [-1.0,  0.0]] # N8
-end
-
-function get_interpolation_polynomial(::Type{Quad8}, xi)
-    [
-        1 xi[2] xi[1] xi[2]^2 xi[1]*xi[2] xi[1]^2 xi[1]*xi[2]^2 xi[1]^2*xi[2]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Quad8}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  1  0        xi[2]  2*xi[1]  xi[2]^2        2*xi[1]*xi[2]
-        0  1  0  2*xi[2]  xi[1]  0        2*xi[1]*xi[2]  xi[1]^2
-    ]
-end
-
-#
-
-type Quad9 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Quad9})
-    Vector{Float64}[
-        [-1.0, -1.0], # N1
-        [ 1.0, -1.0], # N2
-        [ 1.0,  1.0], # N3
-        [-1.0,  1.0], # N4
-        [ 0.0, -1.0], # N5
-        [ 1.0,  0.0], # N6
-        [ 0.0,  1.0], # N7
-        [-1.0,  0.0], # N8
-        [ 0.0,  0.0]] # N9
-end
-
-function get_interpolation_polynomial(::Type{Quad9}, xi)
-    [
-        1 xi[2] xi[1] xi[2]^2 xi[1]*xi[2] xi[1]^2 xi[1]*xi[2]^2 xi[1]^2*xi[2] xi[1]^2*xi[2]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Quad9}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  1  0        xi[2]  2*xi[1]  xi[2]^2        2*xi[1]*xi[2]  2*xi[1]*xi[2]^2
-        0  1  0  2*xi[2]  xi[1]  0        2*xi[1]*xi[2]  xi[1]^2        2*xi[1]^2*xi[2]
-    ]
-end
-
-### 3d elements
-
-type Tet4 <: AbstractElement
+function get_basis{B}(element::Element{B}, ip, time)
+    T = typeof(first(ip))
+    N = zeros(T, 1, length(B))
+    eval_basis!(B, N, tuple(ip...))
+    return N
 end
 
-function get_reference_coordinates(::Type{Tet4})
-    Vector{Float64}[
-        [0.0, 0.0, 0.0], # N1
-        [1.0, 0.0, 0.0], # N2
-        [0.0, 1.0, 0.0], # N3
-        [0.0, 0.0, 1.0]] # N4
+function get_dbasis{B}(element::Element{B}, ip, time)
+    T = typeof(first(ip))
+    dN = zeros(T, size(B)...)
+    eval_dbasis!(B, dN, tuple(ip...))
+    return dN
 end
 
-function get_interpolation_polynomial(::Type{Tet4}, xi)
-    [
-        1.0 xi[1] xi[2] xi[3]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Tet4}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0.0 1.0 0.0 0.0
-        0.0 0.0 1.0 0.0
-        0.0 0.0 0.0 1.0
-    ]
-end
-
-#
-
-type Tet10 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Tet10})
-    Vector{Float64}[
-        [0.0, 0.0, 0.0], # N1
-        [1.0, 0.0, 0.0], # N2
-        [0.0, 1.0, 0.0], # N3
-        [0.0, 0.0, 1.0], # N4
-        [0.5, 0.0, 0.0], # N5
-        [0.5, 0.5, 0.0], # N6
-        [0.0, 0.5, 0.0], # N7
-        [0.0, 0.0, 0.5], # N8
-        [0.5, 0.0, 0.5], # N9
-        [0.0, 0.5, 0.5]] # N10
-end
-
-function get_interpolation_polynomial(::Type{Tet10}, xi)
-    [
-        1.0 xi[3] xi[2] xi[1] xi[3]^2 xi[2]*xi[3] xi[2]^2 xi[1]*xi[3] xi[1]*xi[2] xi[1]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Tet10}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  0  1  0        0      0        xi[3]  xi[2]  2*xi[1]
-        0  0  1  0  0        xi[3]  2*xi[2]  0      xi[1]  0      
-        0  1  0  0  2*xi[3]  xi[2]  0        xi[1]  0      0          
-    ]
-end
-
-#
-
-type Pyr5 <: AbstractElement
-end
-
-
-function get_reference_coordinates(::Type{Pyr5})
-    Vector{Float64}[
-        [-1.0,-1.0,-1.0], # N1
-        [ 1.0,-1.0,-1.0], # N2
-        [ 1.0, 1.0,-1.0], # N3
-        [-1.0, 1.0,-1.0], # N4
-        [ 0.0, 0.0, 1.0]] # N5
-end
-
-function get_interpolation_polynomial(::Type{Pyr5}, xi)
-    [
-        1.0/8.0*(1.0-1.0*xi[1])*(1.0-1.0*xi[2])*(1.0-1.0*xi[3])
-        1.0/8.0*(1.0+1.0*xi[1])*(1.0-1.0*xi[2])*(1.0-1.0*xi[3])
-        1.0/8.0*(1.0+1.0*xi[1])*(1.0+1.0*xi[2])*(1.0-1.0*xi[3])
-        1.0/8.0*(1.0-1.0*xi[1])*(1.0+1.0*xi[2])*(1.0-1.0*xi[3])
-        1.0/2.0*(1.0+xi[3])
-    ]'
-end
-
-function get_interpolation_polynomial(::Type{Pyr5}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        -0.125*(1.0-xi[2])*(1.0-xi[3])   0.125*(1.0-xi[2])*(1.0-xi[3])   0.125*(1.0+xi[2])*(1.0-xi[3])  -0.125*(1.0+xi[2])*(1.0-xi[3])  0.0
-        -0.125*(1.0-xi[1])*(1.0-xi[3])  -0.125*(1.0+xi[1])*(1.0-xi[3])   0.125*(1.0+xi[1])*(1.0-xi[3])   0.125*(1.0-xi[1])*(1.0-xi[3])  0.0
-        -0.125*(1.0-xi[1])*(1.0-xi[2])  -0.125*(1.0+xi[1])*(1.0-xi[2])  -0.125*(1.0+xi[1])*(1.0+xi[2])  -0.125*(1.0-xi[1])*(1.0+xi[2])  0.5
-    ]
-end
-
-#
-
-type Wedge6 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Wedge6})
-    Vector{Float64}[
-        [0.0, 0.0, -1.0], # N1
-        [1.0, 0.0, -1.0], # N2
-        [0.0, 1.0, -1.0], # N3
-        [0.0, 0.0,  1.0], # N4
-        [1.0, 0.0,  1.0], # N5
-        [0.0, 1.0,  1.0]] # N6
-end
-
-function get_interpolation_polynomial(::Type{Wedge6}, x)
-    [
-        1 x[1] x[2] x[3] x[1]*x[3] x[2]*x[3]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Wedge6}, x, ::Type{Val{:partial_derivatives}})
-    [
-        0  1  0  0  x[3]  0   
-        0  0  1  0  0     x[3]
-        0  0  0  1  x[1]  x[2]
-    ]
-end
-
-#
-
-type Wedge15 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Wedge15})
-    Vector{Float64}[
-        [0.0, 0.0, -1.0], # N1
-        [1.0, 0.0, -1.0], # N2
-        [0.0, 1.0, -1.0], # N3
-        [0.0, 0.0,  1.0], # N4
-        [1.0, 0.0,  1.0], # N5
-        [0.0, 1.0,  1.0], # N6
-        [0.5, 0.0, -1.0], # N7
-        [0.5, 0.5, -1.0], # N8
-        [0.0, 0.5, -1.0], # N9
-        [0.5, 0.0,  1.0], # N10
-        [0.5, 0.5,  1.0], # N11
-        [0.0, 0.5,  1.0], # N12
-        [0.0, 0.0,  0.0], # N13
-        [1.0, 0.0,  0.0], # N14
-        [0.0, 1.0,  0.0]] # N15
-end
-
-function get_interpolation_polynomial(::Type{Wedge15}, x)
-    [
-        1 x[1] x[1]^2 x[2] x[1]*x[2] x[2]^2 x[3] x[1]*x[3] x[1]^2*x[3] x[2]*x[3] x[1]*x[2]*x[3] x[2]^2*x[3] x[3]^2 x[1]*x[3]^2 x[2]*x[3]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Wedge15}, x, ::Type{Val{:partial_derivatives}})
-    [
-        0 1 2*x[1] 0 x[2] 0 0 x[3] 2*x[1]*x[3] 0 x[2]*x[3] 0 0 x[3]^2 0
-        0 0 0 1 x[1] 2*x[2] 0 0 0 x[3] x[1]*x[3] 2*x[2]*x[3] 0 0 x[3]^2
-        0 0 0 0 0 0 1 x[1] x[1]^2 x[2] x[1]*x[2] x[2]^2 2*x[3] 2*x[1]*x[3] 2*x[2]*x[3]
-    ]
-end
-
-#
-
-type Hex8 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Hex8})
-    Vector{Float64}[
-        [-1.0, -1.0, -1.0], # N1
-        [ 1.0, -1.0, -1.0], # N2
-        [ 1.0,  1.0, -1.0], # N3
-        [-1.0,  1.0, -1.0], # N4
-        [-1.0, -1.0,  1.0], # N5
-        [ 1.0, -1.0,  1.0], # N6
-        [ 1.0,  1.0,  1.0], # N7
-        [-1.0,  1.0,  1.0]] # N8
-end
-
-function get_interpolation_polynomial(::Type{Hex8}, xi)
-    [
-        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Hex8}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]
-        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]
-        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]
-    ]
-end
-
-#
-
-type Hex20 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Hex20})
-    Vector{Float64}[
-        [-1.0, -1.0, -1.0], # N1
-        [ 1.0, -1.0, -1.0], # N2
-        [ 1.0,  1.0, -1.0], # N3
-        [-1.0,  1.0, -1.0], # N4
-        [-1.0, -1.0,  1.0], # N5
-        [ 1.0, -1.0,  1.0], # N6
-        [ 1.0,  1.0,  1.0], # N7
-        [-1.0,  1.0,  1.0], # N8
-        [ 0.0, -1.0, -1.0], # N9
-        [ 1.0,  0.0, -1.0], # N10
-        [ 0.0,  1.0, -1.0], # N11
-        [-1.0,  0.0, -1.0], # N12
-        [-1.0, -1.0,  0.0], # N13
-        [ 1.0, -1.0,  0.0], # N14
-        [ 1.0,  1.0,  0.0], # N15
-        [-1.0,  1.0,  0.0], # N16
-        [ 0.0, -1.0,  1.0], # N17
-        [ 1.0,  0.0,  1.0], # N18
-        [ 0.0,  1.0,  1.0], # N19
-        [-1.0,  0.0,  1.0]] # N20
-end
-
-function get_interpolation_polynomial(::Type{Hex20}, xi)
-    [
-        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3] xi[3]^2 xi[2]^2 xi[1]^2 xi[2]*xi[3]^2 xi[2]^2*xi[3] xi[1]*xi[3]^2 xi[1]*xi[2]^2 xi[1]^2*xi[3] xi[1]^2*xi[2] xi[1]*xi[2]*xi[3]^2 xi[1]*xi[2]^2*xi[3] xi[1]^2*xi[2]*xi[3]
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Hex20}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]  0        0        2*xi[1]  0              0              xi[3]^2        xi[2]^2        2*xi[1]*xi[3]  2*xi[1]*xi[2]  xi[2]*xi[3]^2        xi[2]^2*xi[3]        2*xi[1]*xi[2]*xi[3]
-        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]  0        2*xi[2]  0        xi[3]^2        2*xi[2]*xi[3]  0              2*xi[1]*xi[2]  0              xi[1]^2        xi[1]*xi[3]^2        2*xi[1]*xi[2]*xi[3]  xi[1]^2*xi[3]
-        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]  2*xi[3]  0        0        2*xi[2]*xi[3]  xi[2]^2        2*xi[1]*xi[3]  0              xi[1]^2        0              2*xi[1]*xi[2]*xi[3]  xi[1]*xi[2]^2        xi[1]^2*xi[2]
-    ]
-end
-
-###
-
-type Hex27 <: AbstractElement
-end
-
-function get_reference_coordinates(::Type{Hex27})
-    Vector{Float64}[
-        [-1.0, -1.0, -1.0], # N1
-        [ 1.0, -1.0, -1.0], # N2
-        [ 1.0,  1.0, -1.0], # N3
-        [-1.0,  1.0, -1.0], # N4
-        [-1.0, -1.0,  1.0], # N5
-        [ 1.0, -1.0,  1.0], # N6
-        [ 1.0,  1.0,  1.0], # N7
-        [-1.0,  1.0,  1.0], # N8
-        [ 0.0, -1.0, -1.0], # N9
-        [ 1.0,  0.0, -1.0], # N10
-        [ 0.0,  1.0, -1.0], # N11
-        [-1.0,  0.0, -1.0], # N12
-        [-1.0, -1.0,  0.0], # N13
-        [ 1.0, -1.0,  0.0], # N14
-        [ 1.0,  1.0,  0.0], # N15
-        [-1.0,  1.0,  0.0], # N16
-        [ 0.0, -1.0,  1.0], # N17
-        [ 1.0,  0.0,  1.0], # N18
-        [ 0.0,  1.0,  1.0], # N19
-        [-1.0,  0.0,  1.0], # N20
-        [ 0.0,  0.0, -1.0], # N21
-        [ 0.0, -1.0,  0.0], # N22
-        [ 1.0,  0.0,  0.0], # N23
-        [ 0.0,  1.0,  0.0], # N24
-        [-1.0,  0.0,  0.0], # N25
-        [ 0.0,  0.0,  1.0], # N26
-        [ 0.0,  0.0,  0.0]] # N27
-end
-
-function get_interpolation_polynomial(::Type{Hex27}, xi)
-    [
-        1 xi[3] xi[2] xi[1] xi[2]*xi[3] xi[1]*xi[3] xi[1]*xi[2] xi[1]*xi[2]*xi[3] xi[3]^2 xi[2]^2 xi[1]^2 xi[2]*xi[3]^2 xi[2]^2*xi[3] xi[1]*xi[3]^2 xi[1]*xi[2]^2 xi[1]^2*xi[3] xi[1]^2*xi[2] xi[2]^2*xi[3]^2 xi[1]*xi[2]*xi[3]^2 xi[1]*xi[2]^2*xi[3] xi[1]^2*xi[3]^2 xi[1]^2*xi[2]*xi[3] xi[1]^2*xi[2]^2 xi[1]*xi[2]^2*xi[3]^2 xi[1]^2*xi[2]*xi[3]^2 xi[1]^2*xi[2]^2*xi[3] xi[1]^2*xi[2]^2*xi[3]^2
-    ]
-end
-
-function get_interpolation_polynomial(::Type{Hex27}, xi, ::Type{Val{:partial_derivatives}})
-    [
-        0  0  0  1  0      xi[3]  xi[2]  xi[2]*xi[3]  0        0        2*xi[1]  0              0              xi[3]^2        xi[2]^2        2*xi[1]*xi[3]  2*xi[1]*xi[2]  0                xi[2]*xi[3]^2        xi[2]^2*xi[3]        2*xi[1]*xi[3]^2  2*xi[1]*xi[2]*xi[3]  2*xi[1]*xi[2]^2  xi[2]^2*xi[3]^2        2*xi[1]*xi[2]*xi[3]^2  2*xi[1]*xi[2]^2*xi[3]  2*xi[1]*xi[2]^2*xi[3]^2
-        0  0  1  0  xi[3]  0      xi[1]  xi[1]*xi[3]  0        2*xi[2]  0        xi[3]^2        2*xi[2]*xi[3]  0              2*xi[1]*xi[2]  0              xi[1]^2        2*xi[2]*xi[3]^2  xi[1]*xi[3]^2        2*xi[1]*xi[2]*xi[3]  0                xi[1]^2*xi[3]        2*xi[1]^2*xi[2]  2*xi[1]*xi[2]*xi[3]^2  xi[1]^2*xi[3]^2        2*xi[1]^2*xi[2]*xi[3]  2*xi[1]^2*xi[2]*xi[3]^2
-        0  1  0  0  xi[2]  xi[1]  0      xi[1]*xi[2]  2*xi[3]  0        0        2*xi[2]*xi[3]  xi[2]^2        2*xi[1]*xi[3]  0              xi[1]^2        0              2*xi[2]^2*xi[3]  2*xi[1]*xi[2]*xi[3]  xi[1]*xi[2]^2        2*xi[1]^2*xi[3]  xi[1]^2*xi[2]        0                2*xi[1]*xi[2]^2*xi[3]  2*xi[1]^2*xi[2]*xi[3]  xi[1]^2*xi[2]^2        2*xi[1]^2*xi[2]^2*xi[3]
-    ]
-end
-
-###
-
-macro create_basis(T)
-    quote
-        T = $T
-        global get_basis, get_dbasis, length, size
-        X = get_reference_coordinates(T)
-        nbasis = length(X)
-        A = zeros(nbasis, nbasis)
-        for i=1:nbasis
-            A[i,:] = get_interpolation_polynomial(T, X[i])
-        end
-        invA = inv(A)
-        function get_basis(element::Element{$T}, ip, time)
-            return get_interpolation_polynomial($T, ip)*invA
-        end
-        function get_dbasis(element::Element{$T}, ip, time)
-            return get_interpolation_polynomial($T, ip, Val{:partial_derivatives})*invA
-        end
-    end
-end
-
-@create_basis Seg2
-@create_basis Seg3
-@create_basis Tri3
-@create_basis Tri6
-@create_basis Tri7
-@create_basis Quad4
-@create_basis Quad8
-@create_basis Quad9
-@create_basis Tet4
-@create_basis Tet10
-@create_basis Pyr5
-@create_basis Wedge6
-@create_basis Wedge15
-@create_basis Hex8
-@create_basis Hex20
-@create_basis Hex27
-
 function inside(::Union{Type{Seg2}, Type{Seg3}, Type{Quad4}, Type{Quad8},
                         Type{Quad9}, Type{Pyr5}, Type{Hex8}, Type{Hex20},
                         Type{Hex27}}, xi)
@@ -628,7 +64,7 @@ function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet
     return all(xi .>= 0.0) && (sum(xi) <= 1.0)
 end
 
-function get_reference_coordinates{E}(element::Element{E})
-    get_reference_coordinates(E)
+function get_reference_coordinates{B}(element::Element{B})
+    return get_reference_element_coordinates(B)
 end
 
diff --git a/src/elements_nurbs.jl b/src/elements_nurbs.jl
deleted file mode 100644
index 7124557..0000000
--- a/src/elements_nurbs.jl
+++ /dev/null
@@ -1,149 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using ForwardDiff
-# TODO: evaluate partial derivatives of basis functions without forwarddiff
-
-""" NURBS segment. """
-type NSeg <: AbstractElement
-    order :: Int
-    knots :: Vector{Float64}
-    weights :: Vector{Float64}
-end
-
-function NSeg()
-    NSeg(1,
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(4))
-end
-
-type NSurf <: AbstractElement
-    order_u :: Int
-    order_v :: Int
-    knots_u :: Vector{Float64}
-    knots_v :: Vector{Float64}
-    weights :: Matrix{Float64}
-end
-
-function NSurf()
-    NSurf(1, 1,
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(2, 2))
-end
-
-type NSolid <: AbstractElement
-    order_u :: Int
-    order_v :: Int
-    order_w :: Int
-    knots_u :: Vector{Float64}
-    knots_v :: Vector{Float64}
-    knots_w :: Vector{Float64}
-    weights :: Array{Float64, 3}
-end
-
-function NSolid()
-    NSolid(1, 1, 1,
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        [-1.0, -1.0, 1.0, 1.0],
-        ones(2, 2, 2))
-end
-
-function NURBS(i, p, u, t)
-    p == 0 && return t[i] <= u <= t[i+1] ? 1.0 : 0.0
-    anom = u-t[i]
-    adenom = t[i+p]-t[i]
-    a = isapprox(adenom, 0.0) ? 0.0 : anom/adenom
-    bnom = t[i+p+1]-u
-    bdenom = t[i+p+1]-t[i+1]
-    b = isapprox(bdenom, 0.0) ? 0.0 : bnom/bdenom
-    result = a*NURBS(i,p-1,u,t) + b*NURBS(i+1,p-1,u,t)
-    return result
-end
-
-function get_basis(element::Element{NSeg}, xi::Vector, time)
-    pu = element.properties.order
-    tu = element.properties.knots
-    w = element.properties.weights
-    nu = length(tu)-pu-1
-    u = xi[1]
-    N = vec([w[j]*NURBS(j,pu,u,tu) for j=1:nu])'
-    return N/sum(N)
-end
-
-function get_basis(element::Element{NSurf}, xi::Vector, time)
-    pu = element.properties.order_u
-    pv = element.properties.order_v
-    tu = element.properties.knots_u
-    tv = element.properties.knots_v
-    w = element.properties.weights
-    nu = length(tu)-pu-1
-    nv = length(tv)-pv-1
-    u, v = xi
-    N = vec([w[i,j]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv) for i=1:nu, j=1:nv])'
-    return N / sum(N)
-end
-
-function get_basis(element::Element{NSolid}, xi::Vector, time)
-    pu = element.properties.order_u
-    pv = element.properties.order_v
-    pw = element.properties.order_w
-    tu = element.properties.knots_u
-    tv = element.properties.knots_v
-    tw = element.properties.knots_w
-    weights = element.properties.weights
-    nu = length(tu)-pu-1
-    nv = length(tv)-pv-1
-    nw = length(tw)-pw-1
-    u, v, w = xi
-    N = vec([weights[i,j,k]*NURBS(i,pu,u,tu)*NURBS(j,pv,v,tv)*NURBS(k,pw,w,tw) for i=1:nu, j=1:nv, k=1:nw])'
-    return N / sum(N)
-end
-
-# TODO: evaluate partial derivatives of basis functions without forwarddiff
-""" Evaluate partial derivatives of basis functions using ForwardDiff. """
-function get_dbasis{E<:Union{NSeg, NSurf, NSolid}}(element::Element{E}, ip, time)
-    xi = isa(ip, IP) ? ip.coords : ip
-    basis(xi) = vec(get_basis(element, xi, time))
-    return ForwardDiff.jacobian(basis, xi)'
-end
-
-function length(element::Element{NSeg})
-    nu = length(element.properties.knots) - element.properties.order - 1
-    return nu
-end
-
-function size(element::Element{NSeg})
-    return (1, length(element))
-end
-
-function length(element::Element{NSurf})
-    nu = length(element.properties.knots_u) - element.properties.order_u - 1
-    nv = length(element.properties.knots_v) - element.properties.order_v - 1
-    return nu*nv
-end
-
-function size(element::Element{NSurf})
-    return (2, length(element))
-end
-
-function length(element::Element{NSolid})
-    nu = length(element.properties.knots_u) - element.properties.order_u - 1
-    nv = length(element.properties.knots_v) - element.properties.order_v - 1
-    nw = length(element.properties.knots_w) - element.properties.order_w - 1
-    return nu*nv*nw
-end
-
-function size(element::Element{NSolid})
-    return (3, length(element))
-end
-
-function is_nurbs(element::Element)
-    return false
-end
-
-function is_nurbs{E<:Union{NSeg, NSurf, NSolid}}(element::Element{E})
-    return true
-end
-
diff --git a/src/fields.jl b/src/fields.jl
index 504e356..aac9857 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -416,7 +416,7 @@ function (field::CCTV)(xi::Vector, time::Number)
     return field.data(time)
 end
 
-function (field::CVTV)(xi::Vector, time::Number)
+function (field::CVTV)(xi, time)
     return field.data(xi, time)
 end
 
diff --git a/src/problems_contact_3d.jl b/src/problems_contact_3d.jl
index 85455a9..a23838c 100644
--- a/src/problems_contact_3d.jl
+++ b/src/problems_contact_3d.jl
@@ -122,7 +122,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri3}, time
     Q3 = create_rotation_matrix(slave_element, time)
 
     # project slave nodes to auxiliary plane (x0, Q)
-    xi = mean(get_reference_coordinates(slave_element))
+    xi = get_mean_xi(slave_element)
     N = vec(get_basis(slave_element, xi, time))
     x0 = N*X1
     n0 = N*n1
@@ -280,7 +280,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
             #la = sub_slave_element("lambda", time)
 
             # create auxiliary plane
-            xi = mean(get_reference_coordinates(sub_slave_element))
+            xi = get_mean_xi(sub_slave_element)
             N = vec(get_basis(sub_slave_element, xi, time))
             x0 = N*X1
             n0 = N*n1
@@ -356,7 +356,7 @@ function assemble!(problem::Problem{Contact}, slave_element::Element{Tri6}, time
         n1 = sub_slave_element("normal", time)
             
         # create auxiliary plane
-        xi = mean(get_reference_coordinates(sub_slave_element))
+        xi = get_mean_xi(sub_slave_element)
         N = vec(get_basis(sub_slave_element, xi, time))
         x0 = N*X1
         n0 = N*n1
diff --git a/src/problems_dirichlet.jl b/src/problems_dirichlet.jl
index b9ff261..1519556 100644
--- a/src/problems_dirichlet.jl
+++ b/src/problems_dirichlet.jl
@@ -56,7 +56,7 @@ function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
             for i=1:field_dim
                 haskey(element, field_name*" $i") || continue
                 ldofs = gdofs[i:field_dim:end]
-                xis = get_reference_coordinates(typeof(element.properties))
+                xis = get_reference_coordinates(element)
                 vals = Float64[]
                 for xi in xis
                     g = element(field_name*" $i", xi, time)
diff --git a/src/problems_mortar.jl b/src/problems_mortar.jl
index 29e9cea..fc89010 100644
--- a/src/problems_mortar.jl
+++ b/src/problems_mortar.jl
@@ -146,7 +146,7 @@ function diagnose_interface(problem::Problem{Mortar}, time::Float64)
         n1 = Field([normals[j] for j in slave_element_nodes])
 
         # project slave nodes to auxiliary plane (x0, Q)
-        xi = mean(get_reference_coordinates(slave_element))
+        xi = get_mean_xi(slave_element)
         N = vec(get_basis(slave_element, xi, time))
         x0 = N*X1
         n0 = N*n1
diff --git a/src/problems_mortar_3d.jl b/src/problems_mortar_3d.jl
index cfb5dca..88dde0d 100644
--- a/src/problems_mortar_3d.jl
+++ b/src/problems_mortar_3d.jl
@@ -103,11 +103,11 @@ function approx_in{T}(q::T, P::Vector{T}; rtol=1.0e-4, atol=0.0)
     return false
 end
 
-function get_polygon_clip(xs, xm, n)
+function get_polygon_clip{T}(xs::Vector{T}, xm::Vector{T}, n::T)
     # objective: search does line xm1 - xm2 clip xs
     nm = length(xm)
     ns = length(xs)
-    P = Vector[]
+    P = T[]
 
     # 1. test is master point inside slave, if yes, add to clip
     for i=1:nm
@@ -330,6 +330,17 @@ function split_quadratic_elements(elements::Vector, time::Float64)
     return new_elements
 end
 
+function get_mean_xi(element::Element)
+    xi = zeros(2)
+    coords = get_reference_coordinates(element)
+    for (xi1,xi2) in coords
+        xi[1] += xi1
+        xi[2] += xi2
+    end
+    xi /= length(coords)
+    return xi
+end
+
 """ Assemble linear surface element to problem.
 
 Dual basis is constructed such that partially integrated slave segments are taken into account in a proper way.
@@ -357,7 +368,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
     n1 = slave_element("normal", time)
 
     # project slave nodes to auxiliary plane (x0, Q)
-    xi = mean(get_reference_coordinates(slave_element))
+    xi = get_mean_xi(slave_element)
     first_slave_element && debug("midpoint xi = $xi")
     N = vec(get_basis(slave_element, xi, time))
     x0 = N*X1
@@ -465,7 +476,7 @@ function assemble!{E<:Union{Tri3, Quad4}}(problem::Problem{Mortar}, slave_elemen
         all_cells = get_cells(P, C0)
         for cell in all_cells
             virtual_element = Element(Tri3, Int[])
-            update!(virtual_element, "geometry", cell)
+            virtual_element.fields["geometry"] = DVTI(cell)
 
             # 5. loop integration point of integration cell
             for ip in get_integration_points(virtual_element, 3)
@@ -579,7 +590,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
             n1 = sub_slave_element("normal", time)
 
             # create auxiliary plane
-            xi = mean(get_reference_coordinates(sub_slave_element))
+            xi = get_mean_xi(sub_slave_element)
             first_slave_element && debug("midpoint xi = $xi")
             N = vec(get_basis(sub_slave_element, xi, time))
             x0 = N*X1
@@ -662,7 +673,7 @@ function assemble!{E<:Union{Tri6}}(problem::Problem{Mortar}, slave_element::Elem
         n1 = sub_slave_element("normal", time)
 
         # create auxiliary plane
-        xi = mean(get_reference_coordinates(sub_slave_element))
+        xi = get_mean_xi(sub_slave_element)
         first_slave_element && debug("midpoint xi = $xi")
         N = vec(get_basis(sub_slave_element, xi, time))
         x0 = N*X1
diff --git a/src/types.jl b/src/types.jl
index d65e65e..a784c99 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -33,6 +33,18 @@ function (point::Point)(field_name, time=0.0)
     point.fields[field_name](time).data
 end
 
+function start(point::Point)
+    return start(point.coords)
+end
+
+function done(point::Point, i)
+    return done(point.coords, i)
+end
+
+function next(point::Point, i)
+    return next(point.coords, i)
+end
+
 function update!{T}(point::Point, field_name, val::Pair{Float64, T})
     if haskey(point, field_name)
         update!(point[field_name], val)
diff --git a/test/test_basis.jl b/test/test_basis.jl
deleted file mode 100644
index 6fec06c..0000000
--- a/test/test_basis.jl
+++ /dev/null
@@ -1,287 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using JuliaFEM.Testing
-importall Base
-import JuliaFEM: get_basis, get_dbasis
-
-type TestElement <: AbstractElement
-end
-
-function get_basis(element::Element{TestElement}, xi, time)
-    1/4*[
-        (1-xi[1])*(1-xi[2])
-        (1+xi[1])*(1-xi[2])
-        (1+xi[1])*(1+xi[2])
-        (1-xi[1])*(1+xi[2])]'
-end
-
-function get_dbasis(element::Element{TestElement}, xi, time)
-    1/4*[
-        -(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
-        -(1-xi[1])   -(1+xi[1])   (1+xi[1])   (1-xi[1])]
-end
-
-function length(element::Element{TestElement})
-    return 4
-end
-
-function size(element::Element{TestElement})
-    return (2, 4)
-end
-
-function get_element()
-    element = Element(TestElement, [1, 2, 3, 4])
-    X = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [1.0, 0.0],
-        3 => [1.0, 1.0],
-        4 => [0.0, 1.0])
-    T = Dict{Int64, Float64}(
-        1 => 1.0,
-        2 => 2.0,
-        3 => 3.0,
-        4 => 4.0)
-    u1 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [0.0, 0.0],
-        3 => [1/4, 0.0],
-        4 => [0.0, 0.0])
-    u2 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0,  0.0],
-        2 => [1.0, -1.0],
-        3 => [2.0,  3.0],
-        4 => [0.0,  0.0])
-    update!(element, "geometry", X)
-    update!(element, "temperature", T)
-    update!(element, "displacement1", u1)
-    update!(element, "displacement2", u2)
-    return element
-end
-
-@testset "spatial interpolation in basis" begin
-    element = get_element()
-    @test isapprox(element([0.0, 0.0], 0.0), 1/4*[1 1 1 1])
-    @test isapprox(element([0.0, 0.0], 1.0), 1/4*[1 1 1 1])
-end
-
-@testset "gradient of shape functions" begin
-    element = get_element()
-    grad = element([0.0, 0.0], 0.0, Val{:Grad})
-    @test isapprox(grad, 1/2*[-1 1 1 -1; -1 -1 1 1])
-end
-
-@testset "interpolation of scalar field in spatial domain" begin
-    # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
-    element = get_element()
-    T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
-    T_interpolated = element("temperature", [0.0, 0.0], 0.0)
-    @test isapprox(T_interpolated, T_known([0.5, 0.5]))
-end
-
-@testset "interpolation of gradient of scalar field in spatial domain" begin
-    # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
-    element = get_element()
-    gradT = element("temperature", [0.0, 0.0], 0.0, Val{:Grad})
-    gradT_expected(X) = [1-2*X[2] 3-2*X[1]]
-    @test isapprox(gradT, gradT_expected([0.5, 0.5]))
-end
-
-@testset "test interpolation of vector field" begin
-    # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
-    element = get_element()
-    u = element("displacement1", [0.0, 0.0], 0.0)
-#   x = X+u
-    u_expected(X) = [1/4*X[1]*X[2], 0]
-#   @test isapprox(x, [9/16, 1/2])
-    @test isapprox(u, u_expected([0.5, 0.5]))
-end
-
-@testset "interpolation of gradient of vector_field" begin
-    # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
-    # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    element = get_element()
-#    displacement = Field(
-#        (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
-#        (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
-    gradu = element("displacement2", [0.0, 0.0], 0.0, Val{:Grad})
-    gradu_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    @test isapprox(gradu, gradu_expected([0.5, 0.5]))
-end
-
-#= TODO: Fix test
-@testset "linear time extrapolation of field" begin
-    #T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
-    T = DVTV()
-    update!(T, 0.0 => [0.0, 0.0, 0.0, 0.0])
-    update!(T, 1.0 => [1.0, 2.0, 3.0, 4.0])
-    @test T(-1.0) == -1.0*[1.0, 2.0, 3.0, 4.0]
-    @test T( 3.0) ==  3.0*[1.0, 2.0, 3.0, 4.0]
-    # when going to \pm infinity, return the last one.
-    @test T(-Inf) ==  0.0*[1.0, 2.0, 3.0, 4.0]
-    @test T(+Inf) ==  1.0*[1.0, 2.0, 3.0, 4.0]
-end
-=#
-
-#= TODO: Fix test
-@testset "constant time extrapolation of field" begin
-    #T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
-    T = DVTV()
-    update!(T, 0.0 => [0.0, 0.0, 0.0, 0.0])
-    update!(T, 1.0 => [1.0, 2.0, 3.0, 4.0])
-    @test isapprox(T(-1.0, Val{:constant}), [0.0, 0.0, 0.0, 0.0])
-    @test isapprox(T( 3.0, Val{:constant}), [1.0, 2.0, 3.0, 4.0])
-end
-=#
-
-#= TODO: Fix test
-@testset "time extrapolation of field with only one timestep" begin
-    T = DVTV()
-    update!(T, 0.0 => [1.0, 2.0, 3.0, 4.0])
-    @test isapprox(T(1.0), [1.0, 2.0, 3.0, 4.0])
-end
-=#
-
-#= TODO: Fix test
-@testset "interpolation in temporal direction" begin
-    field = DCTV()
-    update!(field, 0.0 => 0.0)
-    update!(field, 2.0 => 1.0)
-    update!(field, 4.0 => 2.0)
-    @test isapprox(field(-Inf), 0.0)
-    @test isapprox(field( 0.0), 0.0)
-    @test isapprox(field( 1.0), 0.5)
-    @test isapprox(field( 2.0), 1.0)
-    @test isapprox(field( 3.0), 1.5)
-    @test isapprox(field( 4.0), 2.0)
-    @test isapprox(field(+Inf), 2.0)
-end
-=#
-
-#= TODO: Fix test
-@testset "time derivative interpolation in temporal basis in constant velocity" begin
-    field = DCTV()
-    update!(field, 0.0 => 0.0)
-    update!(field, 2.0 => 1.0)
-    update!(field, 4.0 => 2.0)
-    @test isapprox(field(+Inf, Val{:diff}), 0.5)
-    @test isapprox(field(-Inf, Val{:diff}), 0.5)
-    @test isapprox(field( 0.0, Val{:diff}), 0.5)
-    @test isapprox(field( 0.5, Val{:diff}), 0.5)
-    @test isapprox(field( 1.0, Val{:diff}), 0.5)
-    @test isapprox(field( 1.5, Val{:diff}), 0.5)
-    @test isapprox(field( 2.0, Val{:diff}), 0.5)
-end
-=#
-
-#= TODO: Fix test
-@testset "time derivative interpolation in temporal basis in variable velocity" begin
-    pos = DCTV()
-    for ti in linspace(0, 2, 5)
-        update!(pos, ti => 1/2*ti^2)
-    end
-    # => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
-    velocity = pos(1.0, Val{:diff})
-    v1 = (0.500 - 0.125)/0.5
-    v2 = (1.125 - 0.500)/0.5
-    @test isapprox(velocity, mean([v1, v2])) # = 1.00
-    velocity = pos(2.0, Val{:diff})
-    @test isapprox(velocity, (2.0-1.125)/0.5) # = 1.75
-end
-=#
-
-function test_time_derivative_gradient_interpolation_of_field()
-    # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
-    # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    # => d(u_i,j)/dt = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    X = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [1.0, 0.0],
-        3 => [1.0, 1.0],
-        4 => [0.0, 1.0])
-    u1 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0,  0.0],
-        2 => [0.5, -0.5],
-        3 => [1.0,  1.5],
-        4 => [0.0,  0.0])
-    u2 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0,  0.0],
-        2 => [1.5, -1.5],
-        3 => [3.0,  4.5],
-        4 => [0.0,  0.0])
-    element = Element(TestElement, [1, 2, 3, 4])
-    update!(element, "geometry", X)
-    update!(element, "displacement", 0.5 => u1)
-    update!(element, "displacement", 1.5 => u2)
-
-    xi = [0.0, 0.0]
-    time = 1.2
-    diffgradu = element("displacement", xi, time, Val{:diff}, Val{:Grad})
-    diffgradu_expected(X, t) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
-    @test diffgradu == diffgradu_expected([0.5, 0.5], 1.2)
-end
-
-@testset "some continuum mechanics interpolations" begin
-    X = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [1.0, 0.0],
-        3 => [1.0, 1.0],
-        4 => [0.0, 1.0])
-    u1 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [0.0, 0.0],
-        3 => [0.0, 0.0],
-        4 => [0.0, 0.0])
-    u2 = Dict{Int64, Vector{Float64}}(
-        1 => [0.0, 0.0],
-        2 => [0.0, 0.0],
-        3 => [1/4, 0.0],
-        4 => [0.0, 0.0])
-    element = Element(Quad4, [1, 2, 3, 4])
-    update!(element, "geometry", X)
-    update!(element, "displacement", 0.0 => u1)
-    update!(element, "displacement", 1.0 => u2)
-
-    # from my old home works
-    X = element("geometry", [0.0, 0.0], 1.0)
-    u = element("displacement", [0.0, 0.0], 1.0)
-    x = X + u
-    x_expected = [9/16, 1/2]
-    gradu = element("displacement", [0.0, 0.0], 1.0, Val{:Grad})
-    epsilon = 1/2*(gradu + gradu')
-    rotation = 1/2*(gradu - gradu')
-    k = 0.25
-    epsilon_expected = [
-        X[2]*k      1/2*X[1]*k
-        1/2*X[1]*k  0]
-    rotation_expected = [
-        0          k/2*X[1]
-        -k/2*X[1]         0]
-    F = I + gradu
-    F_expected = [
-        X[2]*k+1 X[1]*k
-        0        1]
-    C = F'*F
-    C_expected = [
-        (X[2]*k+1)^2       (X[2]*k+1)*X[1]*k
-        (X[2]*k+1)*X[1]*k  X[1]^2*k^2+1]
-    E = 1/2*(F'*F - I)
-    E_expected = [
-        1/2*(X[2]*k + 1)^2-1/2   1/2*(X[2]*k+1)*X[1]*k
-        1/2*(X[2]*k + 1)*X[1]*k  1/2*X[1]^2*k^2]
-    U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I)
-#   U_expected = [1.24235 0.13804; 0.13804 1.02149]
-
-    @test isapprox(x, x_expected)
-    @test isapprox(epsilon, epsilon_expected)
-    @test isapprox(rotation, rotation_expected)
-    @test isapprox(F, F_expected)
-    @test isapprox(C, C_expected)
-    @test isapprox(E, E_expected)
-    # TODO: Fix test
-#   @test isapprox(U, U_expected)
-end
-
-
-
diff --git a/test/test_define_new_element.jl b/test/test_define_new_element.jl
deleted file mode 100644
index 8c12b7f..0000000
--- a/test/test_define_new_element.jl
+++ /dev/null
@@ -1,65 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using JuliaFEM.Testing
-importall Base
-import JuliaFEM: get_basis, get_dbasis, get_integration_points
-
-type MyQuad4 <: AbstractElement
-end
-
-function get_basis(element::Element{MyQuad4}, ip, time)
-    1/4*[(1-ip[1])*(1-ip[2]) (1+ip[1])*(1-ip[2]) (1+ip[1])*(1+ip[2]) (1-ip[1])*(1+ip[2])]
-end
-
-function get_dbasis(element::Element{MyQuad4}, ip, time)
-    1/4*[-(1-ip[2])    (1-ip[2])   (1+ip[2])  -(1+ip[2])
-         -(1-ip[1])   -(1+ip[1])   (1+ip[1])   (1-ip[1])]
-end
-
-function get_integration_points(element::MyQuad4)
-    [
-        (1.0, 1.0/sqrt(3.0)*[-1, -1]),
-        (1.0, 1.0/sqrt(3.0)*[ 1, -1]),
-        (1.0, 1.0/sqrt(3.0)*[ 1,  1]),
-        (1.0, 1.0/sqrt(3.0)*[-1,  1])
-    ]
-end
-
-function length(element::Element{MyQuad4})
-    return 4
-end
-
-function size(element::Element{MyQuad4})
-    return (2, 4)
-end
-
-@testset "test new element" begin
-    el = Element(MyQuad4, Int[])
-    el["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
-    el["displacement"] = Vector{Float64}[[0.0,0.0], [0.0,0.0], [1.0,0.0], [0.0,0.0]]
-    @test isapprox(el("geometry", [0.0, 0.0], 0.0), [0.5, 0.5])
-    @test isapprox(el("displacement", [0.0, 0.0], 0.0), [0.25, 0.0])
-    el["temperature thermal conductivity"] = 6.0
-    dim = length(el)
-    K = zeros(dim, dim)
-    A = 0.0
-    time = 0.0
-    for ip in get_integration_points(el)
-        dN = el(ip, time, Val{:Grad})
-        detJ = el(ip, time, Val{:detJ})
-        w = ip.weight*detJ
-        c = el("temperature thermal conductivity", ip, time)
-        K += w*c*dN'*dN
-        A += w
-    end
-    @test isapprox(A, 1.0)
-    K_expected = [
-        4.0  -1.0  -2.0  -1.0
-       -1.0   4.0  -1.0  -2.0
-       -2.0  -1.0   4.0  -1.0
-       -1.0  -2.0  -1.0   4.0]
-    @test isapprox(K, K_expected)
-end
-
diff --git a/test/test_elements_nurbs.jl b/test/test_elements_nurbs.jl
deleted file mode 100644
index 328ee2a..0000000
--- a/test/test_elements_nurbs.jl
+++ /dev/null
@@ -1,28 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using JuliaFEM.Testing
-
-@testset "NSeg interpolate" begin
-    element = Element(NSeg, [1, 2])
-    @test element([0.0], 0.0) == [0.5 0.5]
-    @test size(element) == (1, 2)
-    @test is_nurbs(element)
-    element2 = Element(Seg2, [1, 2])
-    @test !is_nurbs(element2)
-end
-
-@testset "NSurf interpolate" begin
-    element = Element(NSurf, [1, 2, 3, 4])
-    @test element([0.0, 0.0], 0.0) == [0.25 0.25 0.25 0.25]
-    @test size(element) == (2, 4)
-    @test is_nurbs(element)
-end
-
-@testset "NSolid interpolate" begin
-    element = Element(NSolid, [1, 2, 3, 4, 5, 6, 7, 8])
-    @test element([0.0, 0.0, 0.0], 0.0) == [0.125 0.125 0.125 0.125 0.125 0.125 0.125 0.125]
-    @test size(element) == (3, 8)
-    @test is_nurbs(element)
-end
diff --git a/test/test_lagrange.jl b/test/test_lagrange.jl
deleted file mode 100644
index 8fcb931..0000000
--- a/test/test_lagrange.jl
+++ /dev/null
@@ -1,89 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-using JuliaFEM
-using JuliaFEM.Testing
-
-ALL_ELEMENTS = [
-    Seg2, Seg3,
-    Tri3, Tri6, Tri7,
-    Quad4, Quad8, Quad9,
-    Tet4, Tet10,
-    Wedge6,
-    Hex8, Hex20, Hex27
-]
-
-info("basic data for elements implemented so far:")
-for element_type in [Poi1; ALL_ELEMENTS]
-    element = Element(element_type, Int[])
-    element_length = length(element)
-    element_size = size(element)
-    element_description = description(element)
-    info("Element $element_type, description = $element_description, length = $element_length, size = $element_size")
-end
-
-
-ALL_ELEMENTS_NODES = [
-    [1,2], [1,2,3],
-    [1,2,3], [1,2,3,4,5,6], [1,2,3,4,5,6,7],
-    [1,2,3,4], [1,2,3,4,5,6,7,8], [1,2,3,4,5,6,7,8,9],
-    [1,2,3,4], [1,2,3,4,5,6,7,8,9,10],
-    [1,2,3,4,5,6],
-    [1,2,3,4,5,6,7,8],
-    [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,15,17,18,19,20],
-    [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,15,17,18,19,20,
-    21,22,23,24,25,26,27]
-]
-
-@testset "Evaluating basis" begin
-    for (T, nod) in zip(ALL_ELEMENTS,ALL_ELEMENTS_NODES)
-        el = Element(T,nod)
-        nnodes = length(el)
-        for (i, X) in enumerate(get_reference_coordinates(T))
-            Ni = vec(el(X))
-            expected = zeros(nnodes)
-            expected[i] = 1.0
-            @test isapprox(Ni, expected)
-        end
-    end
-end
-
-function get_volume{T<:AbstractElement}(::Type{T},nodes)
-    X = get_reference_coordinates(T)
-    element = Element(T,nodes)
-    update!(element, "geometry", X)
-    V = 0.0
-    for ip in get_integration_points(element)
-        V += ip.weight*element(ip, 0.0, Val{:detJ})
-    end
-    return V
-end
-
-RESULTS = [2.0, 2.0, 0.5, 0.5, 0.5, 2.0^2, 2.0^2,
-           2.0^2, 1/6, 1/6, 1.0, 2.0^3, 2.0^3, 2.0^3,]
-
-@testset "Calculate reference element length/area/volume" begin
-    for (T, nod, res) in zip(ALL_ELEMENTS,ALL_ELEMENTS_NODES,
-                             RESULTS)
-        @test isapprox(get_volume(T,nod),  res)
-    end
-end
-
-SIZES = [(1,2), (1,3), (2,3), (2,6), (2,7), (2,4),
-         (2,8), (2,9), (3,4), (3,10), (3,6), (3,8),
-         (3,20), (3,27)]
-
-@testset "element size" begin
-    for (T, nod, res) in zip(ALL_ELEMENTS,ALL_ELEMENTS_NODES, SIZES)
-        @test size(Element(T,nod)) == res
-    end
-end
-
-@testset "element length" begin
-    for i in 1:length(ALL_ELEMENTS)
-        typ = ALL_ELEMENTS[i]
-        vec = ALL_ELEMENTS_NODES[i]
-        el = Element(typ,vec)
-        @test length(el) == length(vec)
-    end
-end
diff --git a/test/test_solvers.jl b/test/test_solvers.jl
index f1a188d..ab66fde 100644
--- a/test/test_solvers.jl
+++ b/test/test_solvers.jl
@@ -41,7 +41,7 @@ using JuliaFEM.Testing
 
     # Postprocess.
     # Interpolate temperature field along boundary of Γ₁ at time t=1.0
-    xi = [0.0, -1.0]
+    xi = (0.0, )
     X = el2("geometry", xi, 1.0)
     T = el2("temperature", xi, 1.0)
     info("Temperature at point X = $X is T = $T")