diff --git a/.gitignore b/.gitignore
index 1795503..392e664 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,9 @@
 docs/build/html
 *.swp
 *.lnk
+*.mess
+*.log
+*.aux
+*.out
+*.tex
+*.export
diff --git a/src/basis.jl b/src/basis.jl
index df74b1c..f1571df 100644
--- a/src/basis.jl
+++ b/src/basis.jl
@@ -1,88 +1,219 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-abstract AbstractBasis
+abstract Basis <: ContinuousField
 
-""" Defined to dimensionless coordinate ξ∈[-1,1]^n. """
-type SpatialBasis <: AbstractBasis
+### ELEMENT BASIS
+
+""" This is the normal "user defined" basis functions familiar from school books. """
+type ElementBasis <: Basis
     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)
+function Basis(basis::Function, dbasisdxi::Function)
+    return ElementBasis(basis, dbasisdxi)
 end
 
-function call(b::TemporalBasis, value::Number)
-    b.basis(value)
+function Base.call(basis::ElementBasis, xi::Vector, time::Number=0.0)
+    basis.basis(xi)  # passing time does not make much sense actually for this...
 end
 
-function call(b::SpatialBasis, value::Vector)
-    b.basis(value)
+""" Interpolate increment in spatial domain using ElementBasis. """
+function Base.call(basis::ElementBasis, increment::Increment, xi::Vector)
+    basis = basis.basis(xi)
+    sum([basis[i]*increment[i] for i=1:length(increment)])
 end
 
-### INTERPOLATION IN TIME DOMAIN ###
+### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD
 
-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
+""" Here we add field we are wanting to interpolate with ElementBasis. """
+type ElementFieldBasis <: Basis
+    element_basis :: ElementBasis
+    field :: DiscreteField
+    time_extrapolation :: Symbol
+    time_interpolation :: Symbol
+end
+
+function Basis(basis::Function, dbasisdxi::Function, field::DiscreteField,
+               time_extrapolation=:linear, time_interpolation=:linear)
+    element_basis = ElementBasis(basis, dbasisdxi)
+    return ElementFieldBasis(element_basis, field, time_extrapolation,
+                             time_interpolation)
+end
+
+function Base.call(basis::ElementFieldBasis, xi::Vector, time::Number)
+    increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation)
+    return basis.element_basis(increment, xi)
+end
+
+""" Interpolate discrete field in time domain. """
+function Base.call(field::DiscreteField, time::Number,
+                   time_extrapolation::Symbol=:linear,
+                   time_interpolation::Symbol=:linear)
+
+    # special cases, only 1 timestep defined or time = -Inf -> return first ts
+    if (length(field) == 1) || (time == -Inf)
         return field[1][end]
     end
+
+    # special case, time = +Inf -> return last ts
+    if time == +Inf
+        return field[end][end]
+    end
+
+    # very likely we are always near some defined timestep, usually field
+    # defined only on t = 0.0, test neighbourhood for timesteps
+    for i=1:length(field)
+        if isapprox(field[i].time, time)
+            return field[i][end]
+        end
+    end
+
+    # special case: out of time domain in positive direction, very likely
+    # to happen in incremental constitutive models
+    if time > field[end].time
+        if time_extrapolation == :constant
+            # constant time extrapolation, return last field
+            return field[end][end]
+        else
+            # multiple fields, pick last and second last and do linear interpolation
+            f1 = field[end-1]
+            f2 = field[end]
+            dt = abs(f2.time - f1.time)
+            i1 = f1[end]
+            i2 = f2[end]
+            di = i2 - i1
+            increment = Increment(i2 + di./dt * (time-f2.time))
+            return increment
+        end
+    end
+
+    # special case: out of time domain in negative direction
+    if time < field[1].time
+        if time_extrapolation == :constant
+            # constant time extrapolation, return first field
+            return field[1][end]
+        else
+            # multiple fields, pick first and second and do linear interpolation
+            f1 = field[1]
+            f2 = field[2]
+            dt = abs(f2.time - f1.time)
+            i1 = f1[end]
+            i2 = f2[end]
+            di = i2 - i1
+            increment = Increment(i1 - di./dt * (f1.time - time))
+            return increment
+        end
+    end
+
+    # find correct bin and perform interpolation
     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)
+
+    if time_interpolation == :linear
+        t1 = field[i].time
+        t2 = field[i+1].time
+        inc1 = field[i][end]
+        inc2 = field[i+1][end]
+        dt = abs(t2 - t1)
+        t = (time-t1)/dt
+        increment = Increment((1-t)*inc1 + t*inc2)
+        return increment
+    end
+
+    if time_interpolation == :constant
+        # nearest neightbour interpolation, i.e. pick nearest defined field
+        t1 = field[i].time
+        t2 = field[i+1].time
+        dt1 = abs(t1-time)
+        dt2 = abs(t2-time)
+        if dt1 < dt2
+            return field[i][end]
+        else
+            return field[i+1][end]
+        end
+    end
+
 end
 
-function Base.call(field::Field, basis::TemporalBasis, time,
-                   derivative::Type{Val{:derivative}})
+### ELEMENT GRADIENT BASIS = ELEMENT BASIS + GEOMETRY
+
+""" Gradient of ElementBasis, needs geometry information. """
+type ElementGradientBasis <: Basis
+    element_basis :: ElementBasis
+    geometry :: DiscreteField
+    time_extrapolation :: Symbol
+    time_interpolation :: Symbol
+end
+
+function ElementGradientBasis(element_basis::ElementBasis, geometry::DiscreteField)
+    return ElementGradientBasis(element_basis, geometry, :linear, :linear)
+end
+
+function Base.call(basis::ElementGradientBasis, xi::Vector, time::Number=0.0)
+    dbasis = basis.element_basis.dbasisdxi(xi)
+    geometry = basis.geometry(time, basis.time_extrapolation, basis.time_interpolation)
+    J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
+    grad = inv(J)*dbasis
+    return grad
+end
+
+### ELEMENT FIELD GRADIENT BASIS = ELEMENT GRADIENT BASIS + FIELD
+
+""" Gradient of ElementFieldBasis, needs field to interpolate. """
+type ElementFieldGradientBasis <: Basis
+    element_gradient_basis :: ElementGradientBasis
+    field :: DiscreteField
+    time_extrapolation :: Symbol
+    time_interpolation :: Symbol
+end
+
+function ElementFieldGradientBasis(element_gradient_basis::ElementGradientBasis,
+                                   field::DiscreteField)
+    return ElementFieldGradientBasis(element_gradient_basis, field, :linear, :linear)
+end
+
+function Base.call(basis::ElementFieldGradientBasis, xi::Vector, time::Number=0.0)
+    grad = basis.element_gradient_basis(xi, time)
+    increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation)
+    gradf = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
+    return gradf
+end
+
+### INTERPOLATION IN TIME DOMAIN ###
+
+function Base.call(field::DiscreteField, time::Number,
+                   derivative::Type{Val{:derivative}},
+                   time_extrapolation::Symbol=:linear,
+                   time_interpolation::Symbol=:linear)
+
     # FieldSet -> Field -> TimeStep -> Increment -> data
 
+    time_extrapolation == :linear || error("$time_extrapolation not implemented")
+    time_interpolation == :linear || error("$time_interpolation not implemented")
+
     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
+        t1 = field[i]
+        t2 = field[j]
+        J = abs(t2.time - t1.time)
+        t = (time-t1.time)/J
+        result = 1/J*((1-t)*t1[end] + t*t2[end])
+        return Increment(result)
     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
@@ -106,19 +237,3 @@ function Base.call(field::Field, basis::TemporalBasis, time,
 
 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/fields.jl b/src/fields.jl
index 7570634..b64f81e 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -131,7 +131,7 @@ end
 # FIXME: having some serious problems here to get tuple form working.
 
 # 3. DefaultDiscreteField
-type DefaultDiscreteField <: DiscreteField
+immutable DefaultDiscreteField <: DiscreteField
     timesteps :: Vector{TimeStep}
 #=
     function DefaultDiscreteField(data::Array)
@@ -170,6 +170,26 @@ function Base.size(field::DefaultDiscreteField)
     return size(field.timesteps)
 end
 
+function Base.length(field::DefaultDiscreteField)
+    return length(field.timesteps)
+end
+
+function Base.start(::DefaultDiscreteField)
+    return 1
+end
+
+function Base.next(field::DefaultDiscreteField, state)
+    return (field[state+1], state+1)
+end
+
+function Base.done(field::DefaultDiscreteField, state)
+    return state > length(field)
+end
+
+function eltype(::Type{DefaultDiscreteField})
+    return TimeStep
+end
+
 function Base.linearindexing(::Type{DefaultDiscreteField})
     return LinearFast()
 end
@@ -178,10 +198,6 @@ function Base.getindex(field::DefaultDiscreteField, i::Int)
     return field.timesteps[i]
 end
 
-function Base.length(field::DefaultDiscreteField)
-    return length(field.timesteps)
-end
-
 function Base.endof(field::DefaultDiscreteField)
     return endof(field.timesteps)
 end
diff --git a/src/types.jl b/src/types.jl
index bf38ee1..5bb67c7 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -23,6 +23,6 @@ function IntegrationPoint(xi, weight)
     IntegrationPoint(xi, weight, Dict())
 end
 
-call(b::SpatialBasis, ip::IntegrationPoint) = b.basis(ip.xi)
+call(N::ElementBasis, ip::IntegrationPoint) = N(ip.xi)
 
 
diff --git a/test/test_basis.jl b/test/test_basis.jl
index 99c6367..5029080 100644
--- a/test/test_basis.jl
+++ b/test/test_basis.jl
@@ -5,9 +5,163 @@ module BasisTests
 
 using JuliaFEM.Test
 
-using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
+using JuliaFEM
+using JuliaFEM: Basis, ElementGradientBasis, ElementFieldGradientBasis, Field
+
+function get_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])]
+
+    return basis, dbasis
+end
 
 
+function test_basic_interpolation()
+    basis, dbasis = get_basis()
+    b = Basis(basis, dbasis)
+    @test b([0.0, 0.0]) == 1/4*[1 1 1 1]
+    @test b([0.0, 0.0], 1.0) == 1/4*[1 1 1 1]
+end
+
+function test_basic_interpolation_of_field()
+    # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+    temperature = Field(
+        (0.0, [0.0, 0.0, 0.0, 0.0]),
+        (1.0, [1.0, 2.0, 3.0, 4.0]))
+    basis, dbasis = get_basis()
+    b = Basis(basis, dbasis, temperature)
+    T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+    @test b([0.0, 0.0], 0.0) == T([0.5, 0.5], 0.0)
+    @test b([0.0, 0.0], 0.6) == T([0.5, 0.5], 0.6)
+    @test b([0.0, 0.0], 1.0) == T([0.5, 0.5], 1.0) 
+end
+
+function test_linear_time_extrapolation_of_field()
+    temperature = Field(
+        (0.0, [0.0, 0.0, 0.0, 0.0]),
+        (1.0, [1.0, 2.0, 3.0, 4.0]))
+    basis, dbasis = get_basis()
+    b = Basis(basis, dbasis, temperature, :linear)
+    T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+    @test b([0.0, 0.0], -1.0) == T([0.5, 0.5], -1.0)
+    @test b([0.0, 0.0],  3.0) == T([0.5, 0.5],  3.0)
+end
+
+function test_constant_time_extrapolation_of_field()
+    temperature = Field(
+        (0.0, [0.0, 0.0, 0.0, 0.0]),
+        (1.0, [1.0, 2.0, 3.0, 4.0]))
+    basis, dbasis = get_basis()
+    b = Basis(basis, dbasis, temperature, :constant)
+    T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
+    @test b([0.0, 0.0], -1.0) == T([0.5, 0.5], 0.0)
+    @test b([0.0, 0.0],  3.0) == T([0.5, 0.5], 1.0)
+end
+
+function test_time_extrapolation_of_field_with_single_timestep()
+    temperature = Field([1.0, 2.0, 3.0, 4.0])
+    basis, dbasis = get_basis()
+    b = Basis(basis, dbasis, temperature)
+    @test b([0.0, 0.0], 1.0) == mean([1.0, 2.0, 3.0, 4.0])
+end
+
+function test_gradient_interpolation_empty_gradient()
+    X = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
+    geometry = Field(X)
+#   P(X) = [1.0, X[1], X[2], X[1]*X[2]]
+#   basis2, dbasis2 = JuliaFEM.calculate_lagrange_basis(P, X)
+    basis, dbasis = get_basis()
+    N = Basis(basis, dbasis)
+    dN = ElementGradientBasis(N, geometry)
+    @test dN([0.0, 0.0]) == 1/2*[-1 1 1 -1; -1 -1 1 1]
+#   @test dN([0.0, 0.0]) == dbasis2([0.5, 0.5])
+end
+
+function test_gradient_interpolation_of_scalar_field()
+    # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
+    geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
+
+    temperature = Field([1, 2, 3, 4])
+    basis, dbasis = get_basis()
+    N = Basis(basis, dbasis)
+
+    dN = ElementGradientBasis(N, geometry)
+    dT = ElementFieldGradientBasis(dN, temperature)
+    dT_expected(X) = [1-2*X[2] 3-2*X[1]]
+ 
+    @test dT([0.0, 0.0]) == dT_expected([0.5, 0.5])
+end
+
+function test_interpolation_of_vector_field()
+    # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
+    geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
+    displacement = Field(
+        (0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]),
+        (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]))
+
+    basis, dbasis = get_basis()
+    X = Basis(basis, dbasis, geometry)
+    u = Basis(basis, dbasis, displacement)
+    u_expected(X,t) = [1/4*t*X[1]*X[2], 0]
+    # x = X + u
+    x = X([0.0, 0.0], 1.0) + u([0.0, 0.0], 1.0)
+    @test isapprox(x, [9/16, 1/2])
+    @test isapprox(u([0.0, 0.0], 1.0), u_expected([0.5, 0.5], 1.0))
+end
+
+function test_interpolation_of_gradient_of_vector_field()
+    # in unit square, u(X) = t*[X[1]*X[2]/4, X[1]*(X[1]+X[2])/2]
+    # => u_i,j = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
+    geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
+    displacement = Field(
+        (0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]]),
+        (1.0, Vector[[0.0, 0.0], [0.0, 0.5], [0.25, 1.0], [0.0, 0.0]]))
+
+    basis, dbasis = get_basis()
+    N = Basis(basis, dbasis)
+    dN = ElementGradientBasis(N, geometry)
+    dU = ElementFieldGradientBasis(dN, displacement)
+    dU_expected(X, t) = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
+
+    @test isapprox(dU([0.0, 0.0], 1.0), dU_expected([0.5, 0.5], 1.0))
+end
+
+# TODO: how on earth make this work without some serious spaghetti code
+function test_time_derivative_gradient_interpolation_of_field()
+    # in unit square, u(X) = t*[X[1]*X[2]/4, X[1]*(X[1]+X[2])/2]
+    # => u_i,j = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
+    # => d(u_i,j)/dt = [X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
+    geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
+    displacement = Field(
+        (0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]]),
+        (1.0, Vector[[0.0, 0.0], [0.0, 0.5], [0.25, 1.0], [0.0, 0.0]]))
+
+    basis, dbasis = get_basis()
+    N = Basis(basis, dbasis)
+    # wanted
+    #u = Basis(basis, dbasis, displacement)
+    #L = grad(diff(u))
+    #D = 1/2*(L + L')
+    #@text isapprox(D([0.0, 0.0], 1.0), ...)
+
+    xi = [0.0, 0.0]
+    time = 1.0
+    grad = ElementGradientBasis(N, geometry)(xi, time)
+    increment = displacement(time, Val{:derivative}, :linear, :linear)
+    diffgradu = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
+    diffgradu_expected(X, t) = [X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
+    @test diffgradu == diffgradu_expected([0.5, 0.5], 1.0)
+end
+
+#=
 """basic continuum interpolations"""
 function test_basic_interpolations()
 
@@ -42,10 +196,9 @@ function test_basic_interpolations()
     @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])
@@ -90,31 +243,4 @@ function test_interpolation_in_temporal_basis()
     @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_fields.jl b/test/test_fields.jl
index 8ec933c..df4c871 100644
--- a/test/test_fields.jl
+++ b/test/test_fields.jl
@@ -6,8 +6,7 @@ module FieldTests
 
 using JuliaFEM
 using JuliaFEM: Increment, TimeStep, Field, DefaultDiscreteField, FieldSet
-using JuliaFEM: TemporalBasis, SpatialBasis, ContinuousField, DiscreteField
-using JuliaFEM: DefaultContinuousField
+using JuliaFEM: ContinuousField, DiscreteField, DefaultContinuousField
 
 using JuliaFEM.Test
 
@@ -189,7 +188,9 @@ function test_default_discrete_field_quick_way_two_timesteps_with_vector_value()
 end
 
 function test_default_discrete_field_quick_way_set_time_vector_also()
-    f1 = DefaultDiscreteField( (0.5, [1, 2, 3]), (1.0, [3, 4, 5]) )
+    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)
@@ -200,6 +201,20 @@ function test_default_discrete_field_quick_way_set_time_vector_also()
     @test f1[2].time == 1.0
 end
 
+function test_default_discrete_field_for_loop()
+    field = DefaultDiscreteField(
+        (0.5, [1, 2, 3]),
+        (1.0, [3, 4, 5]),
+        (1.5, [4, 5, 6]))
+    timesteps = [ts for ts in field]
+    @test timesteps[1].time == 0.5
+    @test timesteps[2].time == 1.0
+    @test timesteps[3].time == 1.5
+    @test timesteps[1][end] == [1, 2, 3]
+    @test timesteps[2][end] == [3, 4, 5]
+    @test timesteps[3][end] == [4, 5, 6]
+end
+
 function test_default_continuous_field()
 
     function myfield(xi::Vector, time::Float64)
@@ -244,13 +259,10 @@ function test_adding_timesteps()
     fs = FieldSet()
     fs["temperature"] = [1, 2, 3, 4]
     T0 = last(fs["temperature"])  # last increment of last field
-    @debug("last temperature T0 = $T0")
     T1 = Increment(T0 + 1)
-    @debug("typeof T1 = $(typeof(T1))")
     timestep = TimeStep(1.0, Increment[T1])  # new list of increments for timestep
     push!(fs["temperature"], timestep)
     T2 = last(fs["temperature"])
-    @debug("last temperature T2 = $T2")
     @test length(fs["temperature"]) == 2
     @test last(fs["temperature"]) == [2, 3, 4, 5]
     @test fs["temperature"][end].time == 1.0