diff --git a/Makefile b/Makefile
index 06431a0..dfb9c24 100644
--- a/Makefile
+++ b/Makefile
@@ -59,6 +59,14 @@ doctests:
 test:
 	julia -e 'Pkg.test("JuliaFEM")'
 
+# usage: make test_file FILE=test/test_heat.jl
+# or even: watch make test_file FILE=test/test_heat.jl
+test_file:
+	julia -e 'using JuliaFEM.Test; run_test("$(FILE)"); print_test_statistics()'
+
+test_function:
+	julia -e 'using JuliaFEM.Test; run_test("$(FILE)", :$(FUNCTION)); print_test_statistics()'
+
 #coverage_report:
 #	julia -e 'Pkg.test("JuliaFEM"; coverage=true); cd(Pkg.dir("JuliaFEM")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
 
diff --git a/docs/tutorials/2015-06-14-data-structures.ipynb b/docs/tutorials/2015-06-14-data-structures.ipynb
index e06f96e..89d5bae 100644
--- a/docs/tutorials/2015-06-14-data-structures.ipynb
+++ b/docs/tutorials/2015-06-14-data-structures.ipynb
@@ -8,12 +8,13 @@
     "\n",
     "**Author(s)**: Jukka Aho\n",
     "\n",
-    "**Abstract**: Description of basis data structures. In this notebook the concepts of `Increment`, `TimeStep`, `Field`, `FieldSet`, `SpatialBasis`, `TemporalBasis` are intoduced. With combining these atomic structures one is able to form finite elements and interpolate it's fields in time and spatial domain.\n",
+    "**Abstract**: Description of basis data structures. In this notebook the concepts of `Increment`, `TimeStep`, `DiscreteField`, `ContinuousField`, `FieldSet`, `SpatialBasis` and `TemporalBasis` are intoduced. By combining these atomic structures one is able to form finite elements and interpolate it's fields in time and spatial domain.\n",
     "\n",
     "- `Increment` is the most atomic structure. It's a vector-like object with 1 dimension. Each element in `Increment` can be scalar, vector or tensor (2 or 4 order). It's easy to extend `Increment` to have other data types too.\n",
     "- `TimeStep` is container for increments in certain time $t$.\n",
-    "- `Field` is container for timesteps for a single field.\n",
+    "- `DefaultDiscreteField` is container for timesteps for a single field.\n",
     "- `FieldSet` is container for all fields.\n",
+    "- `DefaultContinuousField` can have continuous field variables.\n",
     "\n",
     "## Revision history\n",
     "\n",
@@ -86,8 +87,8 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 1 entry:\n",
-       "  \"temperature\" => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.Inc…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 1 entry:\n",
+       "  \"temperature\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(1.0,JuliaFE…"
       ]
      },
      "execution_count": 1,
@@ -96,7 +97,7 @@
     }
    ],
    "source": [
-    "using JuliaFEM: Increment, TimeStep, Field, FieldSet\n",
+    "using JuliaFEM: Increment, TimeStep, Field, FieldSet, DefaultDiscreteField\n",
     "\n",
     "fs = FieldSet()\n",
     "i1 = Increment([1, 2, 3])\n",
@@ -105,7 +106,7 @@
     "i3 = Increment([2, 3, 4])\n",
     "i4 = Increment([3, 4, 5])\n",
     "t2 = TimeStep(2.0, Increment[i3, i4])\n",
-    "f1 = Field(TimeStep[t1, t2])\n",
+    "f1 = DefaultDiscreteField(TimeStep[t1, t2])\n",
     "fs[\"temperature\"] = f1\n",
     "fs"
    ]
@@ -127,7 +128,10 @@
     {
      "data": {
       "text/plain": [
-       "([3,4,5],JuliaFEM.Increment{Int64})"
+       "3-element JuliaFEM.Increment{Int64}:\n",
+       " 3\n",
+       " 4\n",
+       " 5"
       ]
      },
      "execution_count": 2,
@@ -136,8 +140,7 @@
     }
    ],
    "source": [
-    "increment = fs[\"temperature\"][end][end]\n",
-    "increment, typeof(increment)"
+    "fs[\"temperature\"][2][2]"
    ]
   },
   {
@@ -188,12 +191,12 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 5 entries:\n",
-       "  \"fourth order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"constant scalar field\"  => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"vector field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"scalar field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"second order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 5 entries:\n",
+       "  \"fourth order tensor fi… => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(…\n",
+       "  \"constant scalar field\"  => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(…\n",
+       "  \"vector field\"           => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(…\n",
+       "  \"scalar field\"           => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(…\n",
+       "  \"second order tensor fi… => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(…"
       ]
      },
      "execution_count": 4,
@@ -228,10 +231,10 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 3 entries:\n",
-       "  \"density\"     => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.Inc…\n",
-       "  \"geometry\"    => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.Inc…\n",
-       "  \"temperature\" => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.Inc…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 3 entries:\n",
+       "  \"density\"     => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaFE…\n",
+       "  \"geometry\"    => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaFE…\n",
+       "  \"temperature\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaFE…"
       ]
      },
      "execution_count": 5,
@@ -287,26 +290,20 @@
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 6 entries:\n",
-       "  \"fourth order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"constant scalar field\"  => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"vector field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"scalar field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"time series\"            => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"second order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…"
-      ]
-     },
-     "execution_count": 7,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: time on last timestep: 1.0\n"
+     ]
     }
    ],
    "source": [
-    "fs2[\"time series\"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5]), (1.0, [1, 1, 1, 1])\n",
-    "fs2[\"time series\"][end].time\n",
-    "fs2"
+    "t0 = TimeStep(0.0, Increment([1, 2, 3, 4]))\n",
+    "t1 = TimeStep(0.5, Increment([2, 3, 4, 5]))\n",
+    "t2 = TimeStep(1.0, Increment([1, 1, 1, 1]))\n",
+    "fs2[\"time series\"] = [t0, t1, t2]\n",
+    "\n",
+    "info(\"time on last timestep: $(fs2[\"time series\"][end].time)\")"
    ]
   },
   {
@@ -324,21 +321,16 @@
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 7 entries:\n",
-       "  \"fourth order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"constant scalar field\"  => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"vector field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"scalar field\"           => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"time series\"            => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"time series 2\"          => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…\n",
-       "  \"second order tensor fi… => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[J…"
-      ]
-     },
-     "execution_count": 8,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: MethodError: `convert` has no method matching convert(::Type{Float64}, ::Array{Int64,1})\nThis may have arisen from a call to the constructor Float64(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n  call{T}(::Type{T}, ::Any)\n  convert(::Type{Float64}, !Matched::Int8)\n  convert(::Type{Float64}, !Matched::Int16)\n  ...\nwhile loading In[8], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: MethodError: `convert` has no method matching convert(::Type{Float64}, ::Array{Int64,1})\nThis may have arisen from a call to the constructor Float64(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n  call{T}(::Type{T}, ::Any)\n  convert(::Type{Float64}, !Matched::Int8)\n  convert(::Type{Float64}, !Matched::Int16)\n  ...\nwhile loading In[8], in expression starting on line 1",
+      "",
+      " in convert at /home/jukka/.julia/v0.4/JuliaFEM/src/fields.jl:205",
+      " in call at essentials.jl:56",
+      " in setindex! at dict.jl:641"
+     ]
     }
    ],
    "source": [
@@ -363,9 +355,9 @@
     {
      "data": {
       "text/plain": [
-       "2-element Array{JuliaFEM.TimeStep{T},1}:\n",
-       " JuliaFEM.Increment[[1,2,3,4]]\n",
-       " JuliaFEM.Increment[[2,3,4,5]]"
+       "2-element Array{JuliaFEM.TimeStep,1}:\n",
+       " JuliaFEM.TimeStep(0.0,JuliaFEM.Increment[[1,2,3,4]])\n",
+       " JuliaFEM.TimeStep(1.0,JuliaFEM.Increment[[2,3,4,5]])"
       ]
      },
      "execution_count": 9,
@@ -375,7 +367,7 @@
    ],
    "source": [
     "T0 = first(fs[\"temperature\"])  # pick first timestep (or to be spesific, last increment of first timestep)\n",
-    "T1 = T0 + 1  # create new increment from old one\n",
+    "T1 = Increment(T0 + 1)  # create new increment from old one\n",
     "timestep = TimeStep(1.0, T1)  # create new timestep at t=1.0\n",
     "push!(fs[\"temperature\"], timestep)  # push to field"
    ]
@@ -533,7 +525,7 @@
    "source": [
     "### Creating continuous and discrete fields\n",
     "\n",
-    "In the last section the concept of fields was demonstrated. The `Field` is actually just a typealias for a `DefaultDiscreteField` and `DefaultDiscreteField` is subtype of `DiscreteField` which is subtype of `AbstractField`:"
+    "In previous section the concept of discrete fields was demonstrated. `DefaultDiscreteField` is subtype of `DiscreteField` which is subtype of `Field`:"
    ]
   },
   {
@@ -555,8 +547,8 @@
     }
    ],
    "source": [
-    "using JuliaFEM: AbstractField, DiscreteField\n",
-    "Field <: DiscreteField <: AbstractField"
+    "using JuliaFEM: DefaultDiscreteField, DiscreteField, Field\n",
+    "DefaultDiscreteField <: DiscreteField <: Field"
    ]
   },
   {
@@ -574,7 +566,7 @@
    },
    "outputs": [],
    "source": [
-    "using JuliaFEM: FieldSet, ContinuousField"
+    "using JuliaFEM: FieldSet, DefaultContinuousField, ContinuousField"
    ]
   },
   {
@@ -587,7 +579,7 @@
     {
      "data": {
       "text/plain": [
-       "MyFunnyField()"
+       "MyContinuousField()"
       ]
      },
      "execution_count": 17,
@@ -596,12 +588,12 @@
     }
    ],
    "source": [
-    "type MyFunnyField <: ContinuousField\n",
+    "type MyContinuousField <: ContinuousField\n",
     "end\n",
-    "function Base.call(f::MyFunnyField, xi::Vector, time::Number)\n",
+    "function Base.call(f::MyContinuousField, xi::Vector, time::Number)\n",
     "    time/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])]\n",
     "end\n",
-    "f = MyFunnyField()"
+    "f = MyContinuousField()"
    ]
   },
   {
@@ -637,8 +629,8 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 1 entry:\n",
-       "  \"basis\" => MyFunnyField()"
+       "Dict{ASCIIString,JuliaFEM.Field} with 1 entry:\n",
+       "  \"basis\" => MyContinuousField()"
       ]
      },
      "execution_count": 19,
@@ -648,7 +640,7 @@
    ],
    "source": [
     "fs = FieldSet()\n",
-    "fs[\"basis\"] = MyFunnyField()\n",
+    "fs[\"basis\"] = MyContinuousField()\n",
     "fs"
    ]
   },
@@ -679,7 +671,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Naturally we can pass another fields or even fieldsets to continuous field to make fields depend from each other. Here's another example, where `ContinuousField` takes another field and operates it with some function."
+    "Again we have a shortcut to quickly define continuous fields:"
    ]
   },
   {
@@ -692,7 +684,7 @@
     {
      "data": {
       "text/plain": [
-       "call (generic function with 1246 methods)"
+       "6.0"
       ]
      },
      "execution_count": 21,
@@ -701,11 +693,41 @@
     }
    ],
    "source": [
-    "type MyFunnyContinuousField <: ContinuousField\n",
+    "fs[\"continuous field\"] = (xi, t) -> xi[1]*xi[2]*t\n",
+    "fs[\"continuous field\"]([1.0, 2.0], 3.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Naturally we can pass another fields or even fieldsets to continuous field to make fields depend from each other. Here's another example, where `ContinuousField` takes another field and operates it with some function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "call (generic function with 1259 methods)"
+      ]
+     },
+     "execution_count": 22,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "type MyContinuousField2 <: ContinuousField\n",
     "    basis :: Function\n",
     "    discrete_field :: DiscreteField\n",
     "end\n",
-    "function Base.call(field::MyFunnyContinuousField, xi::Vector, time::Number=1.0)\n",
+    "function Base.call(field::MyContinuousField2, xi::Vector, time::Number=1.0)\n",
     "    data = last(field.discrete_field) # get the last timestep last increment\n",
     "    basis = field.basis(xi) # evaluate basis at point ξ.\n",
     "    sum([basis[i]*data[i] for i=1:length(data)]) # sum results\n",
@@ -721,7 +743,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 23,
    "metadata": {
     "collapsed": false
    },
@@ -729,10 +751,10 @@
     {
      "data": {
       "text/plain": [
-       "MyFunnyContinuousField(basis,JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.Increment[[1,2,3,4]]]))"
+       "MyContinuousField2(basis,JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaFEM.Increment[[1,2,3,4]])]))"
       ]
      },
-     "execution_count": 22,
+     "execution_count": 23,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -745,7 +767,7 @@
     "    (1+xi[1])*(1-xi[2]),\n",
     "    (1+xi[1])*(1+xi[2]),\n",
     "    (1-xi[1])*(1+xi[2])]\n",
-    "fs[\"continuous field\"] = MyFunnyContinuousField(basis, fs[\"discrete field\"])"
+    "fs[\"continuous field\"] = MyContinuousField2(basis, fs[\"discrete field\"])"
    ]
   },
   {
@@ -757,7 +779,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 24,
    "metadata": {
     "collapsed": false
    },
@@ -768,7 +790,7 @@
        "(2.5,[1,2,3,4])"
       ]
      },
-     "execution_count": 23,
+     "execution_count": 24,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -786,7 +808,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 25,
    "metadata": {
     "collapsed": false
    },
@@ -794,19 +816,19 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 2 entries:\n",
-       "  \"continuous field\" => MyFunnyContinuousField(basis,JuliaFEM.DefaultDiscreteFi…\n",
-       "  \"discrete field\"   => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFE…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 2 entries:\n",
+       "  \"continuous field\" => MyContinuousField2(basis,JuliaFEM.DefaultDiscreteField(…\n",
+       "  \"discrete field\"   => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…"
       ]
      },
-     "execution_count": 24,
+     "execution_count": 25,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "T0 = last(fs[\"discrete field\"])\n",
-    "push!(fs[\"discrete field\"], TimeStep(1.0, T0 + 1.0))  # push to field\n",
+    "push!(fs[\"discrete field\"], TimeStep(1.0, Increment(T0 + 1.0)))  # push to field\n",
     "fs"
    ]
   },
@@ -819,7 +841,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 26,
    "metadata": {
     "collapsed": false
    },
@@ -827,10 +849,10 @@
     {
      "data": {
       "text/plain": [
-       "(3.5,[2,3,4,5])"
+       "(3.5,[2.0,3.0,4.0,5.0])"
       ]
      },
-     "execution_count": 25,
+     "execution_count": 26,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -848,7 +870,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 27,
    "metadata": {
     "collapsed": false
    },
@@ -859,7 +881,7 @@
        "(2.5,3.5)"
       ]
      },
-     "execution_count": 26,
+     "execution_count": 27,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -877,7 +899,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 28,
    "metadata": {
     "collapsed": false
    },
@@ -888,28 +910,28 @@
        "getindex (generic function with 126 methods)"
       ]
      },
-     "execution_count": 27,
+     "execution_count": 28,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "using JuliaFEM: DiscreteField\n",
-    "type MyFunnyDiscreteField <: DiscreteField\n",
+    "type MyDiscreteField <: DiscreteField\n",
     "    discrete_points :: Vector\n",
     "    continuous_field :: ContinuousField\n",
     "end\n",
-    "Base.length(field::MyFunnyDiscreteField) = length(field.discrete_points)\n",
-    "Base.endof(field::MyFunnyDiscreteField) = endof(field.discrete_points)\n",
-    "Base.last(field::MyFunnyDiscreteField) = Float64[field[i] for i=1:length(field)]\n",
-    "function Base.getindex(field::MyFunnyDiscreteField, idx::Int64)\n",
+    "Base.length(field::MyDiscreteField) = length(field.discrete_points)\n",
+    "Base.endof(field::MyDiscreteField) = endof(field.discrete_points)\n",
+    "Base.last(field::MyDiscreteField) = Float64[field[i] for i=1:length(field)]\n",
+    "function Base.getindex(field::MyDiscreteField, idx::Int64)\n",
     "    field.continuous_field(field.discrete_points[idx])\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 29,
    "metadata": {
     "collapsed": false
    },
@@ -924,14 +946,14 @@
        " 4.24402"
       ]
      },
-     "execution_count": 28,
+     "execution_count": 29,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "discrete_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]\n",
-    "fs[\"discrete field 2\"] = MyFunnyDiscreteField(discrete_points, fs[\"continuous field\"])\n",
+    "fs[\"discrete field 2\"] = MyDiscreteField(discrete_points, fs[\"continuous field\"])\n",
     "last(fs[\"discrete field 2\"])"
    ]
   },
@@ -952,7 +974,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 30,
    "metadata": {
     "collapsed": false
    },
@@ -963,7 +985,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 31,
    "metadata": {
     "collapsed": false
    },
@@ -976,7 +998,7 @@
        " 0.2"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 31,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -996,7 +1018,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 32,
    "metadata": {
     "collapsed": false
    },
@@ -1008,7 +1030,7 @@
        " 0.25  0.25  0.25  0.25"
       ]
      },
-     "execution_count": 31,
+     "execution_count": 32,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1028,7 +1050,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 33,
    "metadata": {
     "collapsed": false
    },
@@ -1039,7 +1061,7 @@
        "([0.0,0.5,1.0,1.5,2.0],[0.0,0.125,0.5,1.125,2.0])"
       ]
      },
-     "execution_count": 32,
+     "execution_count": 33,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1051,28 +1073,6 @@
     "t, x"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 33,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))"
-      ]
-     },
-     "execution_count": 33,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "x2 = tuple(collect(zip(t, x))...)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 34,
@@ -1093,7 +1093,14 @@
     }
    ],
    "source": [
-    "fs[\"particle position\"] = x2\n",
+    "timesteps = TimeStep[]\n",
+    "for (ti, xi) in zip(t, x)\n",
+    "    increment = Increment(xi)\n",
+    "    push!(timesteps, TimeStep(ti, increment))\n",
+    "end\n",
+    "\n",
+    "#fs[\"particle position\"] = t, 1/2*t.^2\n",
+    "fs[\"particle position\"] = timesteps\n",
     "temporalbasis = TemporalBasis((t) -> [1-t, t], (t) -> [-1, 1])\n",
     "call(fs[\"particle position\"], temporalbasis, 1.0)"
    ]
@@ -1105,6 +1112,29 @@
     "It's also possible to take time derivatives. To do so, call `Field` with `TemporalBasis`, time, and additional argument `Val{:derivative}`. Again, same example:"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1-element Array{Float64,1}:\n",
+       " 1.0"
+      ]
+     },
+     "execution_count": 35,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "call(fs[\"particle position\"], temporalbasis, 1.0, Val{:derivative})"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1123,7 +1153,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 36,
    "metadata": {
     "collapsed": false
    },
@@ -1131,12 +1161,12 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.AbstractField} with 2 entries:\n",
-       "  \"geometry\"     => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.In…\n",
-       "  \"displacement\" => JuliaFEM.DefaultDiscreteField(JuliaFEM.TimeStep[JuliaFEM.In…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 2 entries:\n",
+       "  \"geometry\"     => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaF…\n",
+       "  \"displacement\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,JuliaF…"
       ]
      },
-     "execution_count": 35,
+     "execution_count": 36,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1144,13 +1174,15 @@
    "source": [
     "fs = FieldSet()\n",
     "fs[\"geometry\"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]\n",
-    "fs[\"displacement\"] = (0.0, zeros(2, 4)), (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])\n",
+    "u0 = TimeStep(0.0, Increment(zeros(2, 4)))\n",
+    "u1 = TimeStep(1.0, Increment(Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))\n",
+    "fs[\"displacement\"] = [u0, u1]\n",
     "fs"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 37,
    "metadata": {
     "collapsed": false
    },
@@ -1163,7 +1195,7 @@
        " 0.5   "
       ]
      },
-     "execution_count": 36,
+     "execution_count": 37,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1185,7 +1217,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 38,
    "metadata": {
     "collapsed": false
    },
@@ -1198,7 +1230,7 @@
        " 0.0    0.0  "
       ]
      },
-     "execution_count": 37,
+     "execution_count": 38,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1220,7 +1252,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 39,
    "metadata": {
     "collapsed": false
    },
@@ -1231,7 +1263,7 @@
        "FiniteElement"
       ]
      },
-     "execution_count": 38,
+     "execution_count": 39,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1263,17 +1295,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 40,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
      "ename": "LoadError",
-     "evalue": "LoadError: MethodError: `setindex!` has no method matching setindex!(::FiniteElement, ::Array{Array{Float64,1},1}, ::ASCIIString)\nwhile loading In[39], in expression starting on line 2",
+     "evalue": "LoadError: MethodError: `setindex!` has no method matching setindex!(::FiniteElement, ::Array{Array{Float64,1},1}, ::ASCIIString)\nwhile loading In[40], in expression starting on line 2",
      "output_type": "error",
      "traceback": [
-      "LoadError: MethodError: `setindex!` has no method matching setindex!(::FiniteElement, ::Array{Array{Float64,1},1}, ::ASCIIString)\nwhile loading In[39], in expression starting on line 2",
+      "LoadError: MethodError: `setindex!` has no method matching setindex!(::FiniteElement, ::Array{Array{Float64,1},1}, ::ASCIIString)\nwhile loading In[40], in expression starting on line 2",
       ""
      ]
     }
@@ -1295,6 +1327,15 @@
     "strain_rate = diff(strain)\n",
     "strain_rate([0.0, 0.0], 1.0)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/docs/tutorials/2015-08-29-developing-juliafem.ipynb b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
index d74f407..a485323 100644
--- a/docs/tutorials/2015-08-29-developing-juliafem.ipynb
+++ b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
@@ -47,30 +47,6 @@
     "- variational form, \"principle of minimum potential energy\": there exists some functional  or \"potential function\" $\\Pi$ we are minimizing"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)"
-      ]
-     },
-     "execution_count": 2,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "using Logging\n",
-    "using FactCheck\n",
-    "Logging.configure(level=DEBUG)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -89,13 +65,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 69,
+   "execution_count": 1,
    "metadata": {
     "collapsed": false
    },
    "outputs": [],
    "source": [
-    "using JuliaFEM: Element, Field, FieldSet, Basis"
+    "using JuliaFEM: Element, Basis, FieldSet"
    ]
   },
   {
@@ -107,7 +83,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -116,7 +92,7 @@
     "type MyQuad4 <: Element\n",
     "    connectivity :: Array{Int, 1}\n",
     "    basis :: Basis\n",
-    "    fields :: Dict{ASCIIString, FieldSet}\n",
+    "    fields :: FieldSet\n",
     "end"
    ]
   },
@@ -129,7 +105,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 3,
    "metadata": {
     "collapsed": false
    },
@@ -140,7 +116,7 @@
        "MyQuad4"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -165,7 +141,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 4,
    "metadata": {
     "collapsed": false
    },
@@ -173,10 +149,10 @@
     {
      "data": {
       "text/plain": [
-       "size (generic function with 74 methods)"
+       "size (generic function with 81 methods)"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 4,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -192,6 +168,68 @@
     "Here comes one important thing. We always define our \"things\" so that the first index is dimension, like $(x, y, z)$ or $(\\xi_1, \\xi_2, \\xi_3)$ and second index is basis function number / node id or something similar to that. To motivate this, consider the following example:"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2x4 Array{Int64,2}:\n",
+       " 1  3  5  7\n",
+       " 2  4  6  8"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = [1 2; 3 4; 5 6; 7 8]'"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here, if we consider $x$ as of some field e.g. geometry, our coordinates of first node is $(1, 2)$, second is $(3, 4)$ and so on. Typically on vector field problems the global assembly is something like $(u_1, v_1, u_2, v_2, \\ldots, )$. If fields are defined this way we can easily flatten matrix to vector and back:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "8-element Array{Int64,1}:\n",
+       " 1\n",
+       " 2\n",
+       " 3\n",
+       " 4\n",
+       " 5\n",
+       " 6\n",
+       " 7\n",
+       " 8"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "dim = size(x)\n",
+    "x2 = vec(x)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 7,
@@ -212,68 +250,6 @@
      "output_type": "execute_result"
     }
    ],
-   "source": [
-    "x = [1 2; 3 4; 5 6; 7 8]'"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Here, if we consider $x$ as of some field e.g. geometry, our coordinates of first node is $(1, 2)$, second is $(3, 4)$ and so on. Typically on vector field problems the global assembly is something like $(u_1, v_1, u_2, v_2, \\ldots, )$. If fields are defined this way we can easily flatten matrix to vector and back:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 8,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "8-element Array{Int64,1}:\n",
-       " 1\n",
-       " 2\n",
-       " 3\n",
-       " 4\n",
-       " 5\n",
-       " 6\n",
-       " 7\n",
-       " 8"
-      ]
-     },
-     "execution_count": 8,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "dim = size(x)\n",
-    "x2 = vec(x)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "2x4 Array{Int64,2}:\n",
-       " 1  3  5  7\n",
-       " 2  4  6  8"
-      ]
-     },
-     "execution_count": 9,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
    "source": [
     "reshape(x2, dim)"
    ]
@@ -287,7 +263,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 8,
    "metadata": {
     "collapsed": false
    },
@@ -317,7 +293,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -326,15 +302,22 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "28-Oct 04:26:43:INFO:root:Testing element MyQuad4\n",
-      "28-Oct 04:26:43:INFO:root:element dimension: 2 x 4\n",
-      "28-Oct 04:26:43:INFO:root:Initializing element\n",
-      "28-Oct 04:26:43:INFO:root:basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n",
-      "28-Oct 04:26:43:INFO:root:field val at [0.0,0.0]: 2.5\n",
-      "28-Oct 04:26:43:INFO:root:derivative of basis at [0.0,0.0]: [-0.5 0.5 0.5 -0.5\n",
-      " -0.5 -0.5 0.5 0.5]\n",
-      "28-Oct 04:26:43:INFO:root:field val at [0.0,0.0]: [0.0 2.0]\n",
-      "28-Oct 04:26:43:INFO:root:Element MyQuad4 passed tests.\n"
+      "INFO: Testing element MyQuad4\n",
+      "INFO: element dimension: 2 x 4\n",
+      "INFO: Initializing element\n",
+      "INFO: basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n",
+      "INFO: field val at [0.0,0.0]: 0.0\n"
+     ]
+    },
+    {
+     "ename": "LoadError",
+     "evalue": "LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
+      "",
+      " in call at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:147",
+      " in test_element at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:57"
      ]
     }
    ],
@@ -1291,6 +1274,7 @@
     "    u = basis(\"displacement\", ip, time, variation)\n",
     "    ∇u = dbasis(\"displacement\", ip, time, variation)\n",
     "    F = I + ∇u\n",
+    "    # F = F_e*F_p\n",
     "    b = basis(\"displacement volume load\", ip, time)\n",
     "    E = 1/2*(F'*F - I)\n",
     "    S = λ*trace(E)*I + 2*μ*E\n",
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 6c2524d..a54a662 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -11,6 +11,9 @@ module JuliaFEM
 #using Lexicon
 
 macro debug(msg)
+    if !haskey(ENV, "JuliaFEM_LOG_LEVEL")
+        return :()
+    end
     return :( println("DEBUG: ", $msg) )
 end
 
@@ -32,8 +35,10 @@ function Base.linspace{T<:Array}(X1::T, X2::T, n)
     [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
 end
 
+
+include("fields.jl")  # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
+include("basis.jl")  # interpolation of discrete fields
 include("types.jl")  # type definitions
-#include("interpolate.jl")  # interpolation routines
 
 ### ELEMENTS ###
 include("elements.jl")
@@ -58,6 +63,8 @@ include("solvers.jl")
 include("xdmf.jl")
 include("abaqus_reader.jl")
 
+include("test.jl")
+
 end # module
 
 FEM = JuliaFEM
diff --git a/src/basis.jl b/src/basis.jl
new file mode 100644
index 0000000..df74b1c
--- /dev/null
+++ b/src/basis.jl
@@ -0,0 +1,124 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+abstract AbstractBasis
+
+""" Defined to dimensionless coordinate ξ∈[-1,1]^n. """
+type SpatialBasis <: AbstractBasis
+    basis :: Function
+    dbasisdxi :: Function
+end
+
+typealias Basis SpatialBasis
+
+""" Defined to to interval t∈[0, 1]. """
+type TemporalBasis <: AbstractBasis
+    basis :: Function
+    dbasisdt :: Function
+end
+function TemporalBasis()
+    basis(t) = [1-t, t]
+    dbasis(t) = [-1, 1]
+    return TemporalBasis(basis, dbasis)
+end
+
+function call(b::TemporalBasis, value::Number)
+    b.basis(value)
+end
+
+function call(b::SpatialBasis, value::Vector)
+    b.basis(value)
+end
+
+### INTERPOLATION IN TIME DOMAIN ###
+
+function Base.call(field::Field, basis::TemporalBasis, time)
+    # FieldSet -> Field -> TimeStep -> Increment -> data
+    # special cases, -Inf, +Inf and ~0.0
+    if time > field[end].time
+        return field[end][end]
+    end
+    if (time < field[1].time) || abs(time-field[1].time) < 1.0e-12
+        return field[1][end]
+    end
+    i = length(field)
+    while field[i].time >= time
+        i -= 1
+    end
+    field[i].time == time && return field[i][end]
+    t1 = field[i].time
+    t2 = field[i+1].time
+    inc1 = field[i][end]
+    inc2 = field[i+1][end]
+    # TODO: may there be some reasons for "unphysical" jumps in
+    # fields w.r.t time which should be taken account in some way?
+    # i.e. dt between two fields → 0
+    dt = t2 - t1
+    b = basis.basis((time-t1)/dt)
+    r = Increment[inc1, inc2]
+    return dot(b, r)
+end
+function Base.call(field::DiscreteField, time)
+    return Base.call(field, TemporalBasis(), time)
+end
+
+function Base.call(field::Field, basis::TemporalBasis, time,
+                   derivative::Type{Val{:derivative}})
+    # FieldSet -> Field -> TimeStep -> Increment -> data
+
+    if length(field) == 1
+        # just one timestep, time derivative cannot be evaluated.
+        error("Field length = $(length(field)), cannot evaluate time derivative")
+    end
+
+    function eval_field(i, j)
+        timesteps = TimeStep[field[i], field[j]]
+        increments = Increment[timesteps[1][end], timesteps[2][end]]
+        J = norm(timesteps[2].time - timesteps[1].time)
+        dbasisdt = basis.dbasisdt( (time-timesteps[1].time)/J )
+        return dot(dbasisdt, increments)/J
+    end
+
+    # special cases, +Inf, -Inf, ~0.0
+    if (time > field[end].time) || isapprox(time, field[end].time)
+        return eval_field(endof(field)-1, endof(field))
+    end
+    if (time < field[1].time) || isapprox(time, field[1].time)
+        return eval_field(1, 2)
+    end
+
+    # search for a correct "bin" between time steps
+    i = length(field)
+    while (field[i].time > time) && !isapprox(field[i].time, time)
+        i -= 1
+    end
+
+    if isapprox(field[i].time, time)
+        # This is the hard case, maybe discontinuous time
+        # derivative if linear approximation.
+        # we are on the "mid node" in time axis
+        field1 = eval_field(i-1,i)
+        field2 = eval_field(i,i+1)
+        return 1/2*(field1 + field2)
+    end
+
+    return eval_field(i, i+1)
+
+end
+
+### INTERPOLATION IN SPATIAL DOMAIN ###
+
+function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector)
+    basis = basis.basis(xi)
+    sum([basis[i]*increment[i] for i=1:length(increment)])
+end
+
+function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector,
+                   geometry::Increment, gradient::Type{Val{:gradient}})
+    dbasis = basis.dbasisdxi(xi)
+    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+    grad = inv(J)*dbasis
+    gradf = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
+    return gradf
+end
+
diff --git a/src/elements.jl b/src/elements.jl
index d3c00a3..737600a 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -20,22 +20,22 @@ Raises
 This uses FactCheck and throws exceptions if element is not passing all tests.
 """
 function test_element(element_type)
-    Logging.info("Testing element $element_type")
+    info("Testing element $element_type")
     local element
     dim = nothing
     n = nothing
     try
         dim, n = size(element_type)
     catch
-        Logging.error("Unable to determine element dimensions. Define Base.size(element::Type{$elementtype}) = (dim, nbasis) where dim is spatial dimension of element and nbasis is number of basis functions of element.")
+        error("Unable to determine element dimensions. Define Base.size(element::Type{$elementtype}) = (dim, nbasis) where dim is spatial dimension of element and nbasis is number of basis functions of element.")
     end
-    Logging.info("element dimension: $dim x $n")
+    info("element dimension: $dim x $n")
 
-    Logging.info("Initializing element")
+    info("Initializing element")
     try
         element = element_type(collect(1:n))
     catch
-        Logging.error("""
+        error("""
         Unable to create element with default constructor define function
         $eltype(connectivity) which initializes this element.""")
         return false
@@ -51,15 +51,15 @@ function test_element(element_type)
     dbasis = grad(basis)
     mid = zeros(dim)
     val1 = basis(mid, 0.0)
-    Logging.info("basis at $mid: $val1")
+    info("basis at $mid: $val1")
     val2 = basis("field1", mid, 0.0)
-    Logging.info("field val at $mid: $val2")
+    info("field val at $mid: $val2")
     val3 = dbasis(mid, 0.0)
-    Logging.info("derivative of basis at $mid: $val3")
+    info("derivative of basis at $mid: $val3")
     val4 = dbasis("field1", mid, 0.0)
-    Logging.info("field val at $mid: $val4")
+    info("field val at $mid: $val4")
 
-    Logging.info("Element $element_type passed tests.")
+    info("Element $element_type passed tests.")
 end
 
 """ Get FieldSet from element. """
@@ -75,7 +75,8 @@ Examples
 JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
 """
 function Base.setindex!(element::Element, field_data, field_name)
-    element.fields[field_name] = field_data
+    #element.fields[field_name] = field_data
+    setindex!(element.fields, field_data, field_name)
 end
 
 function get_connectivity(el::Element)
diff --git a/src/interpolate.jl b/src/interpolate.jl
deleted file mode 100644
index 7910521..0000000
--- a/src/interpolate.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: Basis, Field, FieldSet, diff
-
-
-"""
-Interpolate field u using basis N in point xi.
-"""
-function interpolate{T}(N::Basis, u::Field{Vector{T}}, xi::Array{Float64,1})
-    N(xi)*u
-end
-"""
-Interpolate field u using basis N in set of points xi. Convenient function.
-"""
-function interpolate{T}(N::Basis, u::Field{Vector{T}}, xis::Array{Array{Float64,1},1})
-    T[N(xi)*u for xi in xis]
-end
-function interpolate{T}(N::Basis, u::Field{T}, xi::Array{Float64,1})
-    u.values
-end
-
-"""
-Interpolate a field from fieldset for some time t.
-"""
-function interpolate(fields::FieldSet, t::Number)
-    if length(fields) == 0
-        throw("Empty set of fields: $fields")
-    end
-    if t <= fields[1].time
-        return Field(t, fields[1].values)
-    end
-    if t >= fields[end].time
-        return Field(t, fields[end].values)
-    end
-    i = length(fields)
-    while fields[i].time >= t
-        i -= 1
-    end
-    if fields[i].time == t
-        return fields[i]
-    end
-    #Logging.debug("doing linear interpolation between fields $i and $(i+1)")
-    f1 = fields[i]
-    t1 = f1.time
-    f2 = fields[i+1]
-    t2 = f2.time
-    dt = t2 - t1
-    nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
-    f = Field(t, nw)
-    return f
-end
-
-function call(fieldset::FieldSet, time::Number)
-    interpolate(fieldset, time)
-end
-
-function interpolate(basis::Basis, field::Field, ip::IntegrationPoint)
-    interpolate(basis, field, ip.xi)
-end
-
-#function dinterpolate(basis::Basis, u::Field, xi::Array{Float64, 1})
-#    basis.dbasisdxi(xi)*u
-#end
-
diff --git a/src/types.jl b/src/types.jl
index a4552c9..bf38ee1 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,273 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
-
-#abstract AbstractField{T,N} <: AbstractArray{T,N}
-
-abstract AbstractField
-
-abstract DiscreteField <: AbstractField
-
-abstract ContinuousField <: AbstractField
-abstract TimeContinuousField <: ContinuousField
-abstract SpatialContinuousField <: ContinuousField
-abstract TimeAndSpatialContinuousField <: ContinuousField
-# should we introduce time and spatial discontinuous fields
-# for discontinuous galerkin?
-
-### DEFAULT DISCRETE FIELD ###
-
-# 1. Increment
-
-# FIXME: This should be Vector.
-#typealias Increment Vector
-type Increment{T} <: AbstractVector{T}
-    data :: Vector{T}
-end
-Base.size(increment::Increment) = Base.size(increment.data)
-Base.linearindexing(::Type{Increment}) = Base.LinearFast()
-Base.getindex(increment::Increment, i::Int) = increment.data[i]
-Base.setindex!(increment::Increment, v, i::Int) = (increment.data[i] = v)
-Base.similar{T}(increment::Increment, ::Type{T}) = Increment(similar(increment.data))
-Base.dot(v::Number, i::Increment) = v*i
-
-function Base.convert(::Type{Increment}, data::Number)
-    Increment([data])
-end
-function Base.convert{T}(::Type{Increment}, data::Array{T, 2})
-    Increment([data[:,i] for i=1:size(data, 2)])
-end
-function Base.convert{T}(::Type{Increment}, data::Array{T, 3})
-    Increment([data[:,:,i] for i=1:size(data, 3)])
-end
-function Base.convert{T}(::Type{Increment}, data::Array{T, 4})
-    Increment([data[:,:,:,i] for i=1:size(data, 4)])
-end
-function Base.convert{T}(::Type{Increment}, data::Array{T, 5})
-    Increment([data[:,:,:,:,i] for i=1:size(data, 5)])
-end
-function Base.zeros(::Type{Increment}, dims...)
-    Increment(zeros(dims...))
-end
-function Base.vec(increment::Increment)
-    [increment.data...;]
-end
-function Base.similar{T}(increment::Increment{Vector{T}}, data::Vector{T})
-    Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
-end
-
-# 2. TimeStep
-
-type TimeStep{T} <: AbstractVector{T}
-    time :: Float64
-    increments :: Vector{T}
-end
-Base.size(timestep::TimeStep) = Base.size(timestep.increments)
-Base.linearindexing(::Type{TimeStep}) = Base.LinearFast()
-Base.getindex(timestep::TimeStep, i::Int) = timestep.increments[i]
-
-function Base.convert(::Type{TimeStep}, time::Number, increment::Increment)
-    TimeStep(time, Increment[increment])
-end
-
-function Base.push!(timestep::TimeStep, increment::Increment)
-    push!(timestep.increments, increment)
-end
-
-# 3. DefaultDiscreteField
-
-type DefaultDiscreteField <: DiscreteField
-    timesteps :: Vector{TimeStep}
-end
-Base.size(field::DefaultDiscreteField) = Base.size(field.timesteps)
-Base.linearindexing(::Type{DefaultDiscreteField}) = Base.LinearFast()
-Base.getindex(field::DefaultDiscreteField, i::Int) = field.timesteps[i]
-Base.length(field::DefaultDiscreteField) = length(field.timesteps)
-Base.endof(field::DefaultDiscreteField) = endof(field.timesteps)
-Base.first(field::DefaultDiscreteField) = field[1][end]
-Base.last(field::DefaultDiscreteField) = field[end][end]
-function Base.push!(field::DefaultDiscreteField, timestep::TimeStep)
-    push!(field.timesteps, timestep)
-end
-
-
-typealias Field DefaultDiscreteField
-
-### CONTINUOUS FIELDS ###
-
-
-# fix print_matrix
-#function Base.print_matrix(::Base.AbstractIOBuffer, field::ContinuousField, args...)
-    # TODO: anything nice to print?
-#end
-
-### FIELDSET ###
-
-typealias FieldSet Dict{ASCIIString, AbstractField}
-
-"""Quicky add discrete field to fieldset.
-
-Examples
---------
->>> fs = FieldSet()
->>> fs["myfield"] = [1, 2, 3, 4]
-"""
-function Base.convert(::Type{AbstractField}, data::Union{Array, Number})
-    increment = Increment(data)
-    timestep = TimeStep(0.0, Increment[increment])
-    field = DefaultDiscreteField(TimeStep[timestep])
-    return field
-end
-
-""" Quicky add several time steps at once in tuple.
-
-Examples
---------
->>> fs = FieldSet()
->>> fs["myfield"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5])
-
-or
-
->>> fs["myfield"] = [1, 2, 3, 4], [2, 3, 4, 5]
-
-"""
-function Base.convert(::Type{AbstractField}, data::Tuple)
-    timesteps = TimeStep[]
-    for (i, timestep) in enumerate(data)
-        if isa(timestep, Tuple)
-            push!(timesteps, TimeStep(Float64(timestep[1]), Increment(timestep[2])))
-        else
-            push!(timesteps, TimeStep(Float64(i-1), Increment(timestep)))
-        end
-    end
-    return DefaultDiscreteField(timesteps)
-end
-
-### BASIS ###
-
-abstract AbstractBasis
-
-""" Defined to dimensionless coordinate ξ∈[-1,1]^n. """
-type SpatialBasis <: AbstractBasis
-    basis :: Function
-    dbasisdxi :: Function
-end
-
-typealias Basis SpatialBasis
-
-""" Defined to to interval t∈[0, 1]. """
-type TemporalBasis <: AbstractBasis
-    basis :: Function
-    dbasisdt :: Function
-end
-function TemporalBasis()
-    basis(t) = [1-t, t]
-    dbasis(t) = [-1, 1]
-    return TemporalBasis(basis, dbasis)
-end
-
-function call(b::TemporalBasis, value::Number)
-    b.basis(value)
-end
-
-function call(b::SpatialBasis, value::Vector)
-    b.basis(value)
-end
-
-### INTERPOLATION IN TIME DOMAIN ###
-
-function Base.call(field::Field, basis::TemporalBasis, time)
-    # FieldSet -> Field -> TimeStep -> Increment -> data
-    # special cases, -Inf, +Inf and ~0.0
-    if time > field[end].time
-        return field[end][end]
-    end
-    if (time < field[1].time) || abs(time-field[1].time) < 1.0e-12
-        return field[1][end]
-    end
-    i = length(field)
-    while field[i].time >= time
-        i -= 1
-    end
-    field[i].time == time && return field[i][end]
-    t1 = field[i].time
-    t2 = field[i+1].time
-    inc1 = field[i][end]
-    inc2 = field[i+1][end]
-    # TODO: may there be some reasons for "unphysical" jumps in
-    # fields w.r.t time which should be taken account in some way?
-    # i.e. dt between two fields → 0
-    dt = t2 - t1
-    b = basis.basis((time-t1)/dt)
-    r = Increment[inc1, inc2]
-    return dot(b, r)
-end
-function Base.call(field::DiscreteField, time)
-    return Base.call(field, TemporalBasis(), time)
-end
-
-function Base.call(field::Field, basis::TemporalBasis, time,
-                   derivative::Type{Val{:derivative}})
-    # FieldSet -> Field -> TimeStep -> Increment -> data
-
-    if length(field) == 1
-        # just one timestep, time derivative cannot be evaluated.
-        error("Field length = $(length(field)), cannot evaluate time derivative")
-    end
-
-    function eval_field(i, j)
-        timesteps = TimeStep[field[i], field[j]]
-        increments = Increment[timesteps[1][end], timesteps[2][end]]
-        J = norm(timesteps[2].time - timesteps[1].time)
-        dbasisdt = basis.dbasisdt( (time-timesteps[1].time)/J )
-        return dot(dbasisdt, increments)/J
-    end
-
-    # special cases, +Inf, -Inf, ~0.0
-    if (time > field[end].time) || isapprox(time, field[end].time)
-        return eval_field(endof(field)-1, endof(field))
-    end
-    if (time < field[1].time) || isapprox(time, field[1].time)
-        return eval_field(1, 2)
-    end
-
-    # search for a correct "bin" between time steps
-    i = length(field)
-    #while field[i].time >= time + 1.0e-12
-    while (field[i].time > time) && !isapprox(field[i].time, time)
-        i -= 1
-    end
-
-    if isapprox(field[i].time, time)
-        # This is the hard case, maybe discontinuous time
-        # derivative if linear approximation.
-        # we are on the "mid node" in time axis
-        field1 = eval_field(i-1,i)
-        field2 = eval_field(i,i+1)
-        return 1/2*(field1 + field2)
-    end
-
-    return eval_field(i, i+1)
-
-end
-
-### INTERPOLATION IN SPATIAL DOMAIN ###
-
-function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector)
-    basis = basis.basis(xi)
-    sum([basis[i]*increment[i] for i=1:length(increment)])
-end
-
-function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector,
-                   geometry::Increment, gradient::Type{Val{:gradient}})
-    dbasis = basis.dbasisdxi(xi)
-    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
-    grad = inv(J)*dbasis
-    gradf = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
-    return gradf
-end
-
 ### INTEGRATIONPOINT ###
 
 """
diff --git a/test/runtests.jl b/test/runtests.jl
index c1c801e..7d7d09e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -5,11 +5,13 @@
 
 
 using JuliaFEM
-using FactCheck
-using Logging
-@Logging.configure(level=DEBUG)
+using JuliaFEM.Test
 
+#using FactCheck
+#using Logging
+#@Logging.configure(level=DEBUG)
 
+#=
 facts("Testing if somebody used print, println(), @sprint in src directory") do
   # TODO: make better reqular expression. Currently it will match all print words
   lines_with_print = Dict()
@@ -77,6 +79,7 @@ facts("Looking the [src,test] folders *.jl files header information") do
   @fact files_no_license => isempty out_str
 end
 
+#=
 test_files = readdir(Pkg.dir("JuliaFEM")*"/test")
 for test_file in test_files
   Logging.info("checking is $test_file is real test file")
@@ -85,6 +88,7 @@ for test_file in test_files
     include(test_file)
   end
 end
+=#
 
 # Keep this at the end of this file (include statements above this)
 @Logging.configure(level=DEBUG)
@@ -92,3 +96,33 @@ end
 for dic in FactCheck.getstats()
   @debug(dic[1], ": ",dic[2])
 end
+=#
+
+
+### NEW STYLE OF TESTING
+
+using JuliaFEM.Test
+
+function run_tests()
+
+    for test_file in readdir(Pkg.dir("JuliaFEM")*"/test")
+        info("checking is $test_file is real test file")
+        if (startswith(test_file, "test_")) & (endswith(test_file, ".jl"))
+            run_test(test_file)
+        end
+    end
+
+    passed, failed, errors, critical = print_test_statistics()
+
+    # at the very end throw error if something is failed
+    if failed + errors critical != 0
+        error("Some tests has failed. Fix them. Now.")
+        exit(1)
+    else
+        info("""All tests has passed \o/ .""")
+        exit(0)
+    end
+
+end
+
+run_tests()
diff --git a/test/test_basis.jl b/test/test_basis.jl
new file mode 100644
index 0000000..99c6367
--- /dev/null
+++ b/test/test_basis.jl
@@ -0,0 +1,120 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module BasisTests
+
+using JuliaFEM.Test
+
+using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
+
+
+"""basic continuum interpolations"""
+function test_basic_interpolations()
+
+    element = Quad4([1, 2, 3, 4])
+
+    element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0])
+    element["displacement"] = (
+        Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]],
+        Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
+
+    # from my old home works
+    basis = get_basis(element)
+    dbasis = grad(basis)
+    @test isapprox(basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0),  [9/16, 1/2])
+    gradu = dbasis("displacement", [0.0, 0.0], 1.0)
+    epsilon = 1/2*(gradu + gradu')
+    rotation = 1/2*(gradu - gradu')
+    X = basis("geometry", [0.0, 0.0], 1.0)
+    k = 0.25
+    epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
+    rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
+    @test isapprox(epsilon, epsilon_wanted)
+    @test isapprox(rotation, rotation_wanted)
+    F = I + gradu
+    @test isapprox(F, [X[2]*k+1 X[1]*k; 0 1])
+    C = F'*F
+    @test isapprox(C, [(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)
+    @test isapprox(E, [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)
+    @test isapprox(U, [1.24235 0.13804; 0.13804 1.02149])
+end
+
+function test_interpolation_in_temporal_basis()
+    info("testing interpolation on temporal basis")
+    temporalbasis = TemporalBasis((t) -> [1-t, t], (t) -> [-1, 1])
+    @test temporalbasis(0.2) == [0.8, 0.2]
+    i1 = Increment([0.0])
+    i2 = Increment([1.0])
+    i3 = Increment([2.0])
+    t1 = TimeStep(0.0, Increment[i1])
+    t2 = TimeStep(2.0, Increment[i2])
+    t3 = TimeStep(4.0, Increment[i3])
+    field = Field(TimeStep[t1, t2, t3])
+    @test call(field, temporalbasis, -Inf) == [0.0]
+    @test call(field, temporalbasis,  0.0) == [0.0]
+    @test call(field, temporalbasis,  1.0) == [0.5]
+    @test call(field, temporalbasis,  2.0) == [1.0]
+    @test call(field, temporalbasis,  3.0) == [1.5]
+    @test call(field, temporalbasis,  4.0) == [2.0]
+    @test call(field, temporalbasis, +Inf) == [2.0]
+    @test call(field, temporalbasis, +Inf, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis, -Inf, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis,  0.0, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis,  0.5, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis,  1.0, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis,  1.5, Val{:derivative}) == [0.5]
+    @test call(field, temporalbasis,  2.0, Val{:derivative}) == [0.5]
+    fs = FieldSet()
+
+    t = collect(linspace(0, 2, 5))
+    x = 1/2*t.^2
+    x2 = tuple(collect(zip(t, x))...)
+    # => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
+    fs["particle"] = x2
+    position = call(fs["particle"], temporalbasis, 1.0)[1]
+    @test isapprox(position, 0.50)
+    velocity = call(fs["particle"], temporalbasis, 2.0, Val{:derivative})[1]
+    @test isapprox(velocity, (2.0-1.125)/0.5) # = 1.75
+    velocity = call(fs["particle"], temporalbasis, 1.0, Val{:derivative})[1]
+    v1 = (0.500 - 0.125)/0.5
+    v2 = (1.125 - 0.500)/0.5
+    info("v1 = $v1, v2 = $v2")
+    info(mean([v1, v2]))
+    @test isapprox(velocity, mean([v1, v2])) # = 1.00
+
+    # FIXME, returns wrong type.
+    @test isa(position, Increment) == true
+    @test isa(velocity, Increment) == true
+end
+
+function test_interpolation_in_spatial_basis()
+    info("testing interpolation on spatial basis")
+    basis(xi) = 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])]'
+    dbasis(xi) = 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])]
+    spatialbasis = SpatialBasis(basis, dbasis)
+    @test spatialbasis.basis([0.0, 0.0]) == 1/4*[1 1 1 1]
+
+    fs = FieldSet()
+    fs["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
+    fs["displacement"] = (0.0, zeros(2, 4)), (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
+
+    X = call(last(fs["geometry"]), spatialbasis, [0.0, 0.0])
+    u = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0])
+    x = X+u
+    @test X ≈ 1/2*[1, 1]
+    @test x ≈ [9/16, 1/2]
+
+    gradu = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0], last(fs["geometry"]), Val{:gradient})
+    @test gradu ≈ [0.125 0.125; 0.0 0.0]
+end
+
+end
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index 7a6ad38..c33910f 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -1,13 +1,13 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+module ElasticityTests
 
-using Base.Test
-
+using JuliaFEM.Test
 using JuliaFEM: Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem
 
 
-function run()
+function test_elasticity_one_element()
     element = Quad4([1, 2, 3, 4])
     element["geometry"] = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
     element["youngs modulus"] = 500.0
@@ -24,4 +24,5 @@ function run()
     @test disp ≈ -8.77303119819776
 end
 
-run()
+
+end
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 2f7cdb6..06ad936 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -1,18 +1,22 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using FactCheck
-using JuliaFEM: Element, Basis, Field, FieldSet, FunctionSpace
+module ElementTests
+
+using JuliaFEM.Test
+
+using JuliaFEM: Element, Basis, Field, FieldSet, FunctionSpace, test_element
 
 """ Prototype element
 
 This should always pass test_element if everything is ok.
 """
 type MockElement <: Element
-    connectivity :: Array{Int, 1}
+    connectivity :: Vector{Int}
     basis :: Basis
-    fields :: Dict{ASCIIString, FieldSet}
+    fields :: FieldSet
 end
+
 function MockElement(connectivity)
 
     h(xi) = 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])]
@@ -24,33 +28,26 @@ function MockElement(connectivity)
     basis = Basis(h, dh)
     MockElement(connectivity, basis, Dict())
 end
+
 Base.size(element::Type{MockElement}) = (2, 4)
 
-using JuliaFEM: test_element
-facts("test test_element against mock element") do
+"""test test_element against mock element"""
+function test_mockelement()
     test_element(MockElement)
 end
 
-
-facts("test adding fieldsets and fields to element") do
+""" test adding fieldsets and fields to element"""
+function test_add_fields_to_element()
     el = MockElement([1, 2, 3, 4])
- 
-    fieldset = JuliaFEM.FieldSet("geometry")
-    field1 = JuliaFEM.Field(0.0, [0.0, 0.0, 0.0, 0.0])
-    push!(fieldset, field1)
-    field2 = JuliaFEM.Field(1.0, [1.0, 1.0, 1.0, 1.0])
-    push!(fieldset, field2)
-
-    el["geometry"] = fieldset
-    fields = el["geometry"]
-    @fact length(fields) --> 2
-    @fact fields[1] --> field1
-    @fact fields[2] --> field2
+    el["geometry"] = [0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0]
+    field = el["geometry"]
+    @test length(field) == 2  # two time steps
 end
 
+#=
 facts("interpolation of fields in some function space") do
 
-    element = MockElement([1, 2, 3, 4])
+    el = MockElement([1, 2, 3, 4])
     fieldset1 = FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])
     fieldset2 = FieldSet("constant scalar field", [Field(0.0, 1.0)])
     fieldset3 = FieldSet("scalar field", [Field(0.0, [1.0, 2.0, 3.0, 4.0])])
@@ -79,4 +76,6 @@ facts("interpolation of fields in some function space") do
     @fact v("vector field 3", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8, 9+10+11+12]
     @fact v("tensor field 1", xi, t) --> 1/4*[1+2+3+4 5+6+7+8; 9+10+11+12 13+14+15+16]
 end
+=#
 
+end
diff --git a/test/test_fields.jl b/test/test_fields.jl
index 5632dc5..1557989 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -2,115 +2,222 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 
-module TypesTests
+module FieldTests
 
-using JuliaFEM: Increment, TimeStep, AbstractField, DefaultDiscreteField, FieldSet
+using JuliaFEM: Increment, TimeStep, Field, DefaultDiscreteField, FieldSet
 using JuliaFEM: TemporalBasis, SpatialBasis, ContinuousField, DiscreteField
-using JuliaFEM: Field
+using JuliaFEM: DefaultContinuousField
 
+using JuliaFEM.Test
 
-using Base.Test
-
+""" testing Increment """
 function test_increment()
-    info("testing Increment")
-    # testing Increment
+
+    # constant increment
+    I3 = Increment(1)
+    @test isa(I3, Increment)
+    @test length(I3) == 1
+    # FIXME: i don't like 1-length arrays
+    @test I3 == [1]
+
+    # two increments with vector data
     I1 = Increment([1, 2, 3])
     I2 = Increment([2, 3, 4])
-    @test dot(I1, I2) == 20
-    @test dot([1,2,3], I2) == 20
-    @test dot(I1, [2,3,4]) == 20
+    @test length(I1) == 3
+    @test length(I2) == 3
+
+    # basic math
     @test 1/2*(I1+I2) == [1.5, 2.5, 3.5]
     @test I1 + 1 == [2, 3, 4]
     @test I1 - 1 == [0, 1, 2]
     @test I1*3 == [3, 6, 9]
     @test I1+I2 == [3, 5, 7]
+
+    # dot product
+    @test dot(I1, I2) == 20
+    @test dot([1,2,3], I2) == 20
+    @test dot(I1, [2,3,4]) == 20
+    @test dot([1, 2], Increment[I1, I2])
+
+    # similarity
     f = zeros(Increment, 2, 4)
     @test length(f) == 4
-
     g = similar(f, ones(8))
     @test typeof(f) == typeof(g)
     @test length(f) == length(g)
 
-    # promotion of increment
-    @test typeof(I1+1) == typeof(I1)
-    @test typeof(I1-1) == typeof(I1)
-    @test typeof(I1*3) == typeof(I1)
+    # vec
+    @test vec(g) == ones(8)
+
+    # promotion
+    # FIXME: how to do promotion so that modified increment is still increment?
+    @test isa(I1+1, Increment)
+    @test isa(I1-1, Increment)
+    @test isa(3*I1, Increment)
+    @test isa(1/2*I1, Increment)
+    @test isa(I1+I2, Increment)
+    @test isa(I1-I2, Increment)
 
-    # FIXME
-    #@test typeof(I1) == typeof(I1+I2)
-    #@test typeof(I1/2) == typeof(I1)
-    #@test typeof(1/2*S1) == typeof(I1)
 end
-test_increment()
 
+""" testing TimeStep """
 function test_timestep()
-    info("testing TimeStep")
+
+    info("test_timestep(): create empty timestep")
+    ts = TimeStep()
+    @test length(ts) == 0
+    @test ts.time == 0.0
+
+    info("create timestep with two increments")
     i1 = Increment([1, 2, 3])
     i2 = Increment([2, 3, 4])
-    i3 = Increment([2, 3, 4])
-    i4 = Increment([3, 4, 5])
-    t1 = TimeStep(1.0, Increment[i1, i2])
-    t2 = TimeStep(2.0, Increment[i3, i4])
-    @test length(t1) == length(t2) == 2
-    t3 = TimeStep(3.0, i1+1)
+    increments = Increment[i1, i2]
+    ts = TimeStep(1.0, increments)
+    @test length(ts) == 2
+
+    info("create timestep with scalar value")
+    ts = TimeStep(1)
+    @test length(ts) == 1
+    @test ts.time == 0.0
+    @test isa(ts[1], Increment)
+    @test ts[1] == [1]
+
+    info("create timestep compactly for time t=0.0")
+    ts = TimeStep([1, 2, 3])
+    @test length(ts) == 1
+    @test ts.time == 0.0
+    @test isa(ts[1], Increment)
+    @test ts[1] == [1, 2, 3]
+
+    info("create timestep compactly, add three increments compactly for time t=0.0")
+    ts = TimeStep(1, 2, 3)
+    @test length(ts) == 3
+    @test ts.time == 0.0
+    @test isa(ts[1], Increment)
+
+    info("create timestep compactly, add two increments compactly for time t=0.0")
+    ts = TimeStep([1, 2, 3], [2, 3, 4])
+    @test length(ts) == 2
+    @test ts.time == 0.0
+    @test isa(ts[1], Increment)
+    @test ts[1] == [1, 2, 3]
+    @test ts[2] == [2, 3, 4]
+
+    info("create standard timesteps")
+    info(TimeStep(0.5, Increment([1, 2])))
+
+    @test TimeStep(0.5, [1, 2]).time == 0.5
+    @test TimeStep(0.5, [1, 2]) == [1, 2]
+    @test TimeStep(0.5, 1).time == 0.5
+    @test TimeStep(0.5, 1) == [1]
+
 end
-test_timestep()
-
-function test_watta_fak()
-    # TODO: this test will fail if Increment is typealiased to Vector 
-    fs = FieldSet()
-    fs["discrete field"] = [1, 2, 3, 4]
-    T0 = last(fs["discrete field"])
-    info("last discrete field: $T0, ", typeof(T0))
-    T1 = T0 + 1
-    info("adding 1 to discrete field: $T1, ", typeof(T1))
-    ts = TimeStep(1.0, T1)
-    info("creating time step: $ts, ", typeof(ts))
-    push!(fs["discrete field"], ts)
-    info("last discrete field = ", last(fs["discrete field"]))
-
-    info("fieldset: $fs")
-
-    @test last(fs["discrete field"]) == [2, 3, 4, 5]
-end
-test_watta_fak()
 
+""" testing DefaultDiscreteField """
 function test_default_discrete_field()
-    info("testing DefaultDiscreteField")
+
+    info("test_default_discrete_field(): the traditional way")
     i1 = Increment([1, 2, 3])
     i2 = Increment([2, 3, 4])
+    t1 = TimeStep(1.0, Increment[i1, i2])
     i3 = Increment([2, 3, 4])
     i4 = Increment([3, 4, 5])
-    t1 = TimeStep(1.0, Increment[i1, i2])
     t2 = TimeStep(2.0, Increment[i3, i4])
     timesteps = TimeStep[t1, t2]
     f1 = DefaultDiscreteField(timesteps)
     @test length(f1) == 2
-    @test isa(f1, AbstractField) == true
-end
-test_default_discrete_field()
+    @test isa(f1, Field)
+    @test f1[1][1] == [1, 2, 3]
+    @test f1[1][2] == [2, 3, 4]
+    @test f1[2][1] == [2, 3, 4]
+    @test f1[2][2] == [3, 4, 5]
+    @test f1[1].time == 1.0
+    @test f1[2].time == 2.0
 
+    info("test_default_discrete_field(): quick way, this creates one timestep with vector value")
+    f1 = DefaultDiscreteField([1, 2, 3])
+    info("f1 = $f1")
+    @test isa(f1[1], TimeStep)
+    @test isa(f1[1][1], Increment)
+    @test f1[1][1] == [1, 2, 3]
+    @test f1[1].time == 0.0
+
+    info("test_default_discrete_field(): quick way, two timesteps with constant value")
+    f1 = DefaultDiscreteField(1, 2)
+    @test length(f1) == 2
+    @test isa(f1[1], TimeStep)
+    @test isa(f1[2], TimeStep)
+    @test isa(f1[1][1], Increment)
+    @test isa(f1[2][1], Increment)
+    @test f1[1][1] == [1]
+    @test f1[2][1] == [2]
+    @test f1[1].time == 0.0
+    @test f1[2].time == 1.0
+
+    info("test_default_discrete_field(): quick way, one timestep with scalar value")
+    f1 = DefaultDiscreteField(1)
+    @test length(f1) == 1
+    @test isa(f1[1], TimeStep)
+    @test isa(f1[1][1], Increment)
+    @test f1[1][1] == [1]
+    @test f1[1].time == 0.0
+
+    info("test_default_discrete_field(): quick way, two timesteps with vector value")
+    f1 = DefaultDiscreteField([1, 2, 3], [3, 4, 5])
+    @test length(f1) == 2
+    @test isa(f1[1], TimeStep)
+    @test isa(f1[2], TimeStep)
+    @test isa(f1[1][1], Increment)
+    @test isa(f1[2][1], Increment)
+    @test f1[1][1] == [1, 2, 3]
+    @test f1[2][1] == [3, 4, 5]
+    @test f1[1].time == 0.0
+    @test f1[2].time == 1.0
+
+    info("test_default_discrete_field(): quick way, set time vector also")
+    f1 = DefaultDiscreteField( (0.5, [1, 2, 3]), (1.0, [3, 4, 5]) )
+    @test isa(f1[1], TimeStep)
+    @test isa(f1[2], TimeStep)
+    @test isa(f1[1][1], Increment)
+    @test isa(f1[2][1], Increment)
+    @test f1[1][1] == [1, 2, 3]
+    @test f1[2][1] == [3, 4, 5]
+    @test f1[1].time == 0.5
+    @test f1[2].time == 1.0
+
+end
+
+""" testing DefaultContinuousField """
+function test_default_continuous_field()
+
+    function myfield(xi::Vector, time::Float64)
+        time/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
+
+    f = DefaultContinuousField(myfield)
+    @test f([0.0, 0.0], 1.0) == [0.25 0.25 0.25 0.25]
+    
+end
+
+""" testing FieldSet """
 function test_fieldset()
-    i1 = Increment([1, 2, 3])
-    i2 = Increment([2, 3, 4])
-    i3 = Increment([2, 3, 4])
-    i4 = Increment([3, 4, 5])
-    t1 = TimeStep(1.0, Increment[i1, i2])
-    t2 = TimeStep(2.0, Increment[i3, i4])
-    timesteps = TimeStep[t1, t2]
-    f1 = DefaultDiscreteField(timesteps)
-    info("testing adding discrete field to FieldSet")
+
+    info("test_fieldset(): testing adding discrete field to FieldSet")
     fs = FieldSet()
-    fs["temperature"] = f1
+    fs["temperature"] = DefaultDiscreteField([1, 2, 3])
     @test length(fs) == 1
 
-    info("testing adding discrete fields quickly")
-    # the easy way
+    info("test_fieldset(): testing adding discrete fields quickly")
     fs2 = FieldSet()
     fs2["temperature"] = [1, 2, 3, 4]
     @test fs2["temperature"][end][end] == [1, 2, 3, 4]
     @test last(fs2["temperature"]) == [1, 2, 3, 4]
 
+    info("test_fieldset(): testing adding all kind of discrete fields")
     fs2 = FieldSet()
     fs2["constant scalar field"] = 1
     fs2["scalar field"] = [1, 2, 3, 4]
@@ -120,45 +227,61 @@ function test_fieldset()
     timestep = fs2["vector field"][end]
     @test timestep.time == 0.0
 
-    info("testing adding timesteps")
-    # add another timestep
+    info("test_fieldset(): testing adding timesteps")
     fs = FieldSet()
     fs["temperature"] = [1, 2, 3, 4]
     T0 = last(fs["temperature"])  # last increment of last field
-    info("last temperature = $T0")
-    T1 = T0 + 1
-    @test typeof(T0) == typeof(T1)
+    info("last temperature T0 = $T0")
+    T1 = Increment(T0 + 1)
+    info("typeof T1 = $(typeof(T1))")
     timestep = TimeStep(1.0, Increment[T1])  # new list of increments for timestep
     push!(fs["temperature"], timestep)
     T2 = last(fs["temperature"])
-    info("last temperature = $T2")
+    info("last temperature T2 = $T2")
     @test last(fs["temperature"]) == [2, 3, 4, 5]
-    # or more easily
+
+    info("test_fieldset(): testing adding timesteps compactly")
     timestep = TimeStep(2.0, T1)
     push!(fs["temperature"], timestep)
     @test length(fs["temperature"].timesteps) == 3
 
-    info("test adding several time steps at once")
+    info("test_fieldset(): test adding several time steps at once without time vector")
     fs3 = FieldSet()
-    fs3["time series 1"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5]), (1.0, [1, 1, 1, 1])
-    @test fs3["time series 1"][end].time == 1.0
-    fs3["time series 2"] = [1, 2, 3, 4], [2, 3, 4, 5], [1, 1, 1, 1]
-    @test fs3["time series 2"][end].time == 2.0
-end
-test_fieldset()
+    fs3["time series 2"] = [1, 2, 3, 4], [2, 3, 4, 5]
+    info(fs3)
+    @test fs3["time series 2"][1].time == 0.0
+    @test fs3["time series 2"][2].time == 1.0
+    @test fs3["time series 2"][1][end] == [1, 2, 3, 4]
+    @test fs3["time series 2"][2][end] == [2, 3, 4, 5]
+
+    info("test_fieldset(): test adding several time steps at once with time vector")
+    fs3 = FieldSet()
+    fs3["time series 1"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5])
+    @test fs3["time series 1"][1].time == 0.0
+    @test fs3["time series 1"][2].time == 0.5
+    @test fs3["time series 1"][1][end] == [1, 2, 3, 4]
+    @test fs3["time series 1"][2][end] == [2, 3, 4, 5]
+
+    info("test_fieldset(): adding continuous field")
+    fs = FieldSet()
+    fs["continuous field"] = (xi, t) -> xi[1]*xi[2]*t
+    @test fs["continuous field"]([1.0, 2.0], 3.0) == 6.0
 
-type MyFunnyContinuousField <: ContinuousField
-    basis :: Function
-    discretefield :: DiscreteField
 end
-function Base.call(field::MyFunnyContinuousField, xi::Vector, time::Number=1.0)
-    data = last(field.discretefield) # get the last timestep last increment
+
+type MyContinuousField <: ContinuousField
+    basis :: Function
+    discrete_field :: DiscreteField
+end
+function Base.call(field::MyContinuousField, xi::Vector, time::Number=1.0)
+    data = last(field.discrete_field) # get the last timestep last increment
     info("data = $data, typeof data = $(typeof(data))")
     basis = time*field.basis(xi) # evaluate basis at point ξ.
     sum([basis[i]*data[i] for i=1:length(data)]) # sum results
 end
+
+""" testing ContinuousField """
 function test_continuous_field()
-    info("testing continuous field")
     fs = FieldSet()
     fs["discrete field"] = [1, 2, 3, 4]
     basis(xi) = 1/4*[
@@ -166,7 +289,7 @@ function test_continuous_field()
         (1+xi[1])*(1-xi[2]),
         (1+xi[1])*(1+xi[2]),
         (1-xi[1])*(1+xi[2])]
-    fs["continuous field"] = MyFunnyContinuousField(basis, fs["discrete field"])
+    fs["continuous field"] = MyContinuousField(basis, fs["discrete field"])
     @test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(1+2+3+4)
     T0 = last(fs["discrete field"])
     T1 = T0 + 1.0
@@ -174,21 +297,20 @@ function test_continuous_field()
     push!(fs["discrete field"], TimeStep(1.0, T0+1.0))
     @test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(2+3+4+5)
 end
-test_continuous_field()
 
-type MyFunnyDiscreteField <: DiscreteField
+type MyDiscreteField <: DiscreteField
     discrete_points :: Vector
-    continuousfield :: ContinuousField
+    continuous_field :: ContinuousField
 end
-Base.length(field::MyFunnyDiscreteField) = length(field.discrete_points)
-Base.endof(field::MyFunnyDiscreteField) = endof(field.discrete_points)
-Base.last(field::MyFunnyDiscreteField) = Float64[field[i] for i=1:length(field)]
-function Base.getindex(field::MyFunnyDiscreteField, idx::Int64)
-    field.continuousfield(field.discrete_points[idx])
+Base.length(field::MyDiscreteField) = length(field.discrete_points)
+Base.endof(field::MyDiscreteField) = endof(field.discrete_points)
+Base.last(field::MyDiscreteField) = Float64[field[i] for i=1:length(field)]
+function Base.getindex(field::MyDiscreteField, idx::Int64)
+    field.continuous_field(field.discrete_points[idx])
 end
 
+""" testing DiscreteField """
 function test_discrete_field()
-    info("testing discrete field")
     fs = FieldSet()
     fs["discrete field"] = [1, 2, 3, 4]
     basis(xi) = 1/4*[
@@ -196,98 +318,35 @@ function test_discrete_field()
         (1+xi[1])*(1-xi[2]),
         (1+xi[1])*(1+xi[2]),
         (1-xi[1])*(1+xi[2])]
-    fs["continuous field"] = MyFunnyContinuousField(basis, fs["discrete field"])
+    fs["continuous field"] = MyContinuousField(basis, fs["discrete field"])
     discrete_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]
-    fs["discrete field 2"] = MyFunnyDiscreteField(discrete_points, fs["continuous field"])
+    fs["discrete field 2"] = MyDiscreteField(discrete_points, fs["continuous field"])
     @test last(fs["discrete field 2"]) ≈ [
         1.7559830641437073,
         2.0893163974770410,
         2.9106836025229590,
         3.2440169358562922]
 end
-test_discrete_field()
 
-
-function test_interpolation_in_temporal_basis()
-    info("testing interpolation on temporal basis")
-    temporalbasis = TemporalBasis((t) -> [1-t, t], (t) -> [-1, 1])
-    @test temporalbasis(0.2) == [0.8, 0.2]
-    i1 = Increment([0.0])
-    i2 = Increment([1.0])
-    i3 = Increment([2.0])
-    t1 = TimeStep(0.0, Increment[i1])
-    t2 = TimeStep(2.0, Increment[i2])
-    t3 = TimeStep(4.0, Increment[i3])
-    field = Field(TimeStep[t1, t2, t3])
-    @test call(field, temporalbasis, -Inf) == [0.0]
-    @test call(field, temporalbasis,  0.0) == [0.0]
-    @test call(field, temporalbasis,  1.0) == [0.5]
-    @test call(field, temporalbasis,  2.0) == [1.0]
-    @test call(field, temporalbasis,  3.0) == [1.5]
-    @test call(field, temporalbasis,  4.0) == [2.0]
-    @test call(field, temporalbasis, +Inf) == [2.0]
-    @test call(field, temporalbasis, +Inf, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis, -Inf, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis,  0.0, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis,  0.5, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis,  1.0, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis,  1.5, Val{:derivative}) == [0.5]
-    @test call(field, temporalbasis,  2.0, Val{:derivative}) == [0.5]
-    fs = FieldSet()
-
-    t = collect(linspace(0, 2, 5))
-    x = 1/2*t.^2
-    x2 = tuple(collect(zip(t, x))...)
-    # => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
-    fs["particle"] = x2
-    position = call(fs["particle"], temporalbasis, 1.0)[1]
-    @test position ≈ 0.50
-    velocity = call(fs["particle"], temporalbasis, 2.0, Val{:derivative})[1]
-    @test velocity ≈ (2.0-1.125)/0.5 # = 1.75
-    velocity = call(fs["particle"], temporalbasis, 1.0, Val{:derivative})[1]
-    v1 = (0.500 - 0.125)/0.5
-    v2 = (1.125 - 0.500)/0.5
-    info("v1 = $v1, v2 = $v2")
-    info(mean([v1, v2]))
-    @test velocity ≈ mean([v1, v2]) # = 1.00
-
-    # FIXME, returns wrong type.
-    #=
-    @test isa(position, Increment) == true
-    @test isa(velocity, Increment) == true
-    =#
+function test_field_conversion()
+    i1 = Increment([1, 2, 3])
+    i2 = Increment([2, 3, 4])
+    t1 = TimeStep(1.0, Increment[i1, i2])
+    i3 = Increment([2, 3, 4])
+    i4 = Increment([3, 4, 5])
+    t2 = TimeStep(2.0, Increment[i3, i4])
+    timesteps = TimeStep[t1, t2]
+    info("timesteps = $timesteps")
+    f1 = Field(timesteps)
+    info("field = $f1")
+    @test length(f1) == 2
+    @test isa(f1, Field)
+    @test f1[1][1] == [1, 2, 3]
+    @test f1[1][2] == [2, 3, 4]
+    @test f1[2][1] == [2, 3, 4]
+    @test f1[2][2] == [3, 4, 5]
+    @test f1[1].time == 1.0
+    @test f1[2].time == 2.0
 end
-test_interpolation_in_temporal_basis()
-
-function test_interpolation_in_spatial_basis()
-    info("testing interpolation on spatial basis")
-    basis(xi) = 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])]'
-    dbasis(xi) = 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])]
-    spatialbasis = SpatialBasis(basis, dbasis)
-    @test spatialbasis.basis([0.0, 0.0]) == 1/4*[1 1 1 1]
-
-    fs = FieldSet()
-    fs["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
-    fs["displacement"] = (0.0, zeros(2, 4)), (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
-
-    X = call(last(fs["geometry"]), spatialbasis, [0.0, 0.0])
-    u = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0])
-    x = X+u
-    @test X ≈ 1/2*[1, 1]
-    @test x ≈ [9/16, 1/2]
-
-    gradu = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0], last(fs["geometry"]), Val{:gradient})
-    @test gradu ≈ [0.125 0.125; 0.0 0.0]
-end
-test_interpolation_in_spatial_basis()
-
-
-println("test_fields.jl: all test passing.")
 
 end
diff --git a/test/test_global_assembly.jl b/test/test_global_assembly.jl
index 6862590..c48aced 100644
--- a/test/test_global_assembly.jl
+++ b/test/test_global_assembly.jl
@@ -1,11 +1,14 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+module GlobalAssemblyTests
+
+using JuliaFEM.Test
 using JuliaFEM: Quad4, Seg2, FieldSet, Field, PlaneHeatProblem
 using JuliaFEM: initialize_global_assembly, calculate_global_assembly!
-using FactCheck
 
-facts("assemble a simple two element problem and solve") do 
+"""assemble a simple two element problem and solve"""
+function test_asssembly()
     el1 = Quad4([1, 2, 3, 4])
     el1["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
     el1["temperature thermal conductivity"] = 6.0
@@ -28,5 +31,7 @@ facts("assemble a simple two element problem and solve") do
     A = lufact(global_assembly.stiffness_matrix[free_dofs, free_dofs])
     b = full(global_assembly.force_vector)[free_dofs]
     u = A \ b
-    @fact u --> roughly([101.0, 101.0])
+    @test isapprox(u, roughly([101.0, 101.0]))
+end
+
 end
diff --git a/test/test_heat.jl b/test/test_heat.jl
index 0fb5b4d..56f5e1a 100644
--- a/test/test_heat.jl
+++ b/test/test_heat.jl
@@ -3,10 +3,17 @@
 
 # unit tests for heat equations
 
-using FactCheck
-using JuliaFEM: Seg2, Quad4, Field, FieldSet, DC2D4, initialize_local_assembly, calculate_local_assembly!, DC2D2
+module HeatTests  # always wrap tests to module ending with "Tests"
 
-facts("tests on [0x1]x[0x1] domain") do
+using JuliaFEM.Test  # always use JuliaFEM.Test, not Base.Test
+
+using JuliaFEM: Seg2, Quad4, Field, FieldSet, DC2D4,
+                initialize_local_assembly, calculate_local_assembly!,
+                DC2D2
+
+
+"tests on [0x1]x[0x1] domain"
+function test_one_element()  # always start test function with name test_
 
     # volume element
     element = Quad4([1, 2, 3, 4])
@@ -29,7 +36,7 @@ facts("tests on [0x1]x[0x1] domain") do
     fdofs = [1, 2]
     A = la.stiffness_matrix
     b = la.force_vector
-    @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0])
+    @test isapprox(A[fdofs, fdofs] \ b[fdofs], [1.0, 1.0])
 
     # Set constant flux g=6 on boundary. Accurate solution is
     # u(x,y) = x which equals T=1 on boundary.
@@ -37,7 +44,8 @@ facts("tests on [0x1]x[0x1] domain") do
 
     calculate_local_assembly!(la, boundary_equation, "temperature")
     b = la.force_vector
-    @fact A[fdofs, fdofs] \ b[fdofs] --> roughly([1.0, 1.0])
+    @test isapprox(A[fdofs, fdofs] \ b[fdofs], [1.0, 1.0])  # always use @test to test things.
 
 end
 
+end
diff --git a/test/test_interpolation.jl b/test/test_interpolation.jl
deleted file mode 100644
index ef39013..0000000
--- a/test/test_interpolation.jl
+++ /dev/null
@@ -1,39 +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: get_basis, grad, FieldSet, Field, Quad4
-using FactCheck
-
-
-facts("basic continuum interpolations") do
-
-    element = Quad4([1, 2, 3, 4])
-
-    element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
-    element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0])
-    element["displacement"] = (
-        Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]],
-        Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
-
-    # from my old home works
-    basis = get_basis(element)
-    dbasis = grad(basis)
-    @fact basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0) --> [9/16, 1/2]
-    gradu = dbasis("displacement", [0.0, 0.0], 1.0)
-    epsilon = 1/2*(gradu + gradu')
-    rotation = 1/2*(gradu - gradu')
-    X = basis("geometry", [0.0, 0.0], 1.0)
-    k = 0.25
-    epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
-    rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
-    @fact epsilon --> roughly(epsilon_wanted)
-    @fact rotation --> roughly(rotation_wanted)
-    F = I + gradu
-    @fact F --> [X[2]*k+1 X[1]*k; 0 1]
-    C = F'*F
-    @fact C --> [(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)
-    @fact E --> [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)
-    #@fact U --> roughly([1.24235 0.13804; 0.13804 1.02149])
-end
diff --git a/test/test_math.jl b/test/test_math.jl
index 37caf1b..1d97002 100644
--- a/test/test_math.jl
+++ b/test/test_math.jl
@@ -1,10 +1,15 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using JuliaFEM: interpolate
+module MathTests
+using JuliaFEM.Test
 
-using FactCheck
+function test_math()
+    @test 1+1 == 2
+end
 
+#using JuliaFEM: interpolate
+#=
 facts("test interpolation of different field variables") do
     N(xi) = [
       (1-xi[1])*(1-xi[2])/4
@@ -31,3 +36,6 @@ facts("test interpolation of different field variables") do
     @fact interpolate(F5, dNdξ, [0.0, 0.0]) --> [5.0 0.0; 0.0 0.5]
     @fact interpolate(F6, N, [0.0, 0.0]) --> 36
 end
+=#
+
+end