diff --git a/src/FEMBase.jl b/src/FEMBase.jl
index d45f8a9..2a9fcbc 100644
--- a/src/FEMBase.jl
+++ b/src/FEMBase.jl
@@ -36,5 +36,6 @@ export Poi1, Seg2, Seg3, Tri3, Tri6, Tri7, Quad4, Quad8, Quad9,
 export update!, add_elements!, add!, get_gdofs, group_by_element_type,
        get_unknown_field_name, get_unknown_field_dimension,
        get_integration_points, initialize!, assemble!
+export DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV
 
 end
diff --git a/src/assembly.jl b/src/assembly.jl
index 2220905..09072fd 100644
--- a/src/assembly.jl
+++ b/src/assembly.jl
@@ -44,6 +44,11 @@ function assemble!(problem::Problem, time=0.0; auto_initialize=true)
     return true
 end
 
+function assemble!{P}(assembly::Assembly, problem::Problem{P}, element::Element, time)
+    warn("One must define assemble! function for problem of type $P. Not doing anything.")
+    return nothing
+end
+
 function assemble!(assembly::Assembly, problem::Problem, elements::Vector{Element}, time)
     warn("assemble!() this is default assemble operation, decreased performance can be expected without preallocation of memory!")
     for element in elements
diff --git a/src/problems.jl b/src/problems.jl
index f53e31c..f2eae3a 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -395,6 +395,11 @@ function get_unknown_field_dimension(problem::Problem)
     return problem.dimension
 end
 
+function get_unknown_field_name{P<:AbstractProblem}(::Type{P})
+    warn("The name of unknown field (e.g. displacement, temperature, ...) of the problem type must be given by defining function `get_unknown_field_name`")
+    return "N/A"
+end
+
 """ Return the name of the unknown field of this problem. """
 function get_unknown_field_name{P}(problem::Problem{P})
     return get_unknown_field_name(P)
diff --git a/test/runtests.jl b/test/runtests.jl
index 0ba8309..e1e4244 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,6 +6,13 @@ using TimerOutputs
 const to = TimerOutput()
 
 test_files = String[]
+push!(test_files, "test_add_elements.jl")
+push!(test_files, "test_common_failures.jl")
+push!(test_files, "test_elements.jl")
+push!(test_files, "test_elements_2.jl")
+push!(test_files, "test_fields.jl")
+push!(test_files, "test_integration_points.jl")
+push!(test_files, "test_sparse.jl")
 
 @testset "FEMBase.jl" begin
     for fn in test_files
diff --git a/test/test_add_elements.jl b/test/test_add_elements.jl
new file mode 100644
index 0000000..0641c0d
--- /dev/null
+++ b/test/test_add_elements.jl
@@ -0,0 +1,16 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+type Dummy <: FEMBase.FieldProblem
+end
+
+@testset "add elements to problem" begin
+    problem = Problem(Dummy, "test", 2)
+    element = Element(Quad4, [1, 2, 3, 4])
+    elements = [element]
+    add_elements!(problem, elements)
+    @test problem.elements[1] == element
+end
diff --git a/test/test_common_failures.jl b/test/test_common_failures.jl
new file mode 100644
index 0000000..0af229e
--- /dev/null
+++ b/test/test_common_failures.jl
@@ -0,0 +1,15 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+type Dummy <: FieldProblem end
+
+@testset "no unknown field name or assemble!-function defined" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    pr = Problem(Dummy, "problem", 2)
+    add_elements!(pr, [el])
+    assemble!(pr)
+    @test true
+end
diff --git a/test/test_elements.jl b/test/test_elements.jl
new file mode 100644
index 0000000..6e241e1
--- /dev/null
+++ b/test/test_elements.jl
@@ -0,0 +1,75 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using Base.Test
+
+using FEMBase: group_by_element_type
+
+@testset "add time dependent field to element" begin
+    el = Element(Seg2, [1, 2])
+    u1 = Vector{Float64}[[0.0, 0.0], [0.0, 0.0]]
+    u2 = Vector{Float64}[[1.0, 1.0], [1.0, 1.0]]
+    update!(el, "displacement", 0.0 => u1)
+    update!(el, "displacement", 1.0 => u2)
+    @test length(el["displacement"]) == 2
+    @test isapprox(el("displacement", [0.0], 0.0), [0.0, 0.0])
+    @test isapprox(el("displacement", [0.0], 0.5), [0.5, 0.5])
+    @test isapprox(el("displacement", [0.0], 1.0), [1.0, 1.0])
+    el2 = Element(Poi1, [1])
+    update!(el2, "force 1", 0.0 => 1.0)
+end
+
+@testset "add CVTV field to element" begin
+    el = Element(Seg2, [1, 2])
+    f(xi, time) = xi[1]*time
+    update!(el, "my field", f)
+    v = el("my field", [1.0], 2.0)
+    @test isapprox(v, 2.0)
+end
+
+@testset "add DCTI to element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    update!(el, "displacement load", DCTI([4.0, 8.0]))
+    @test isa(el["displacement load"], DCTI)
+    @test !isa(el["displacement load"].data, DCTI)
+    update!(el, "displacement load 2", [4.0, 8.0])
+    @test isa(el["displacement load 2"], DCTI)
+    update!(el, "temperature", [1.0, 2.0, 3.0, 4.0])
+    @test isa(el["temperature"], DVTI)
+    @test isapprox(el("displacement load", [0.0, 0.0], 0.0), [4.0, 8.0])
+end
+
+@testset "interpolate DCTI from element" begin
+    el = Element(Seg2, [1, 2])
+    update!(el, "foobar", 1.0)
+    fb = el("foobar", [0.0], 0.0)
+    @test isa(fb, Float64)
+    @test isapprox(fb, 1.0)
+end
+
+@testset "add elements to elements" begin
+    el1 = Element(Seg2, [1, 2])
+    el2 = Element(Seg2, [3, 4])
+    update!(el1, "master elements", [el2])
+    lst = el1("master elements", 0.0)
+    @test isa(lst, Vector)
+end
+
+@testset "extend basis" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    expected = [
+        0.25 0.00 0.25 0.00 0.25 0.00 0.25 0.00
+        0.00 0.25 0.00 0.25 0.00 0.25 0.00 0.25]
+    @test isapprox(el([0.0, 0.0], 0.0, 2), expected)
+end
+
+@testset "group elements" begin
+    e1 = Element(Seg2, [1, 2])
+    e2 = Element(Quad4, [1, 2, 3, 4])
+    elements = [e1, e2]
+    r = group_by_element_type(elements)
+    @test length(r) == 2
+    @test first(r[Element{Seg2}]) == e1
+    @test first(r[Element{Quad4}]) == e2
+end
diff --git a/test/test_elements_2.jl b/test/test_elements_2.jl
new file mode 100644
index 0000000..f929125
--- /dev/null
+++ b/test/test_elements_2.jl
@@ -0,0 +1,54 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using FEMBase: get_local_coordinates, inside
+using Base.Test
+
+@testset "inverse isoparametric mapping" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    X1 = el("geometry", [0.1, 0.2], time)
+    xi = get_local_coordinates(el, X1, time)
+    X2 = el("geometry", xi, time)
+    info("X1 = $X1, X2 = $X2")
+    @test isapprox(X1, X2)
+end
+
+@testset "inside of linear element" begin
+    el = Element(Quad4, [1, 2, 3, 4])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [1.0, 1.0],
+        4 => [0.0, 1.0])
+    update!(el, "geometry", X)
+    time = 0.0
+    @test inside(el, [0.5, 0.5], time) == true
+    @test inside(el, [1.0, 0.5], time) == true
+    @test inside(el, [1.0, 1.0], time) == true
+    @test inside(el, [1.01, 1.0], time) == false
+    @test inside(el, [1.0, 1.01], time) == false
+end
+
+@testset "inside of quadratic element" begin
+    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
+    X = Dict{Int64, Vector{Float64}}(
+        1 => [0.0, 0.0],
+        2 => [1.0, 0.0],
+        3 => [0.0, 1.0],
+        4 => [0.5, 0.2],
+        5 => [0.8, 0.6],
+        6 => [-0.2, 0.5])
+    update!(el, "geometry", X)
+    p = [0.94, 0.3] # visually checked to be inside
+    @test inside(el, p, 0.0) == true
+    p = [-0.2, 0.8] # visually checked to be outside
+    @test inside(el, p, 0.0) == false
+end
diff --git a/test/test_fields.jl b/test/test_fields.jl
new file mode 100644
index 0000000..0f8718f
--- /dev/null
+++ b/test/test_fields.jl
@@ -0,0 +1,202 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using FEMBase: Field
+using Base.Test
+
+@testset "discrete, constant, time invariant field" begin
+    @test DCTI(0.0).data == 0.0
+    @test isa(Field(0.0), DCTI)
+    f = DCTI(0.0)
+    update!(f, 1.0)
+    @test isapprox(f, DCTI(1.0))
+    @test isapprox(f, 1.0)
+    @test 2*f == 2.0 # multiply by constant
+    @test f(1.0) == 1.0 # time interpolation
+    @test isapprox(reshape([2.0],1,1)*f, 2.0) # wanted behavior?
+end
+
+@testset "discrete, variable, time invariant field" begin
+    @test DVTI([1.0, 2.0]).data == [1.0, 2.0]
+    @test isa(Field([1.0, 2.0]), DVTI)
+    
+    f = DVTI(zeros(2))
+    update!(f, [2.0, 3.0])
+    @test isapprox(f.data, [2.0, 3.0])
+    @test length(f) == 2
+
+    # slicing
+    @test isapprox(f[1], 2.0)
+    @test isapprox(f[[1, 2]], [2.0, 3.0])
+
+    # boolean comparison and multiplying by a constant
+    @test f == DVTI([2.0, 3.0])
+    @test isapprox(2*f, [4.0, 6.0])
+
+    f3 = 2*f
+    @test isa(f3, DVTI)
+    @test f3+f == 3*f
+    @test f3-f == f
+
+    # spatial interpolation
+    N = [1.0, 2.0]
+    @test isapprox(N*f, 8.0)
+
+    # time interpolation
+    @test isapprox(f(1.0), [2.0, 3.0])
+
+    # spatial interpolation of vector valued variable field
+    f2 = DVTI(Vector[[1.0, 2.0], [3.0, 4.0]])
+    @test isapprox(f2[1], [1.0, 2.0])
+    @test isapprox(f2[2], [3.0, 4.0])
+    @test length(f2) == 2
+    @test isapprox(N*f2, [1.0, 2.0] + [6.0, 8.0])
+
+    # iteration of DVTI field
+    s = zeros(2)
+    for j in f2
+        s += j
+    end
+    @test isapprox(s, [4.0, 6.0])
+
+    @test vec(f2) == [1.0, 2.0, 3.0, 4.0]
+    @test isapprox([1.0 2.0]*f, [8.0]'')
+
+    new_data = [2.0, 3.0, 4.0, 5.0]
+    f4 = similar(f2, new_data)
+    @test isa(f4, DVTI)
+    @test isapprox(f4.data[1], [2.0, 3.0])
+    @test isapprox(f4.data[2], [4.0, 5.0])
+end
+
+@testset "discrete, constant, time-variant field" begin
+    f = Field(0.0 => 1.0)
+    @test isa(f, DCTV)
+    @test last(f).time == 0.0
+    @test last(f).data == 1.0
+    update!(f, 0.0 => 2.0)
+    @test last(f).time == 0.0
+    @test last(f).data == 2.0
+    @test length(f) == 1
+    update!(f, 1.0 => 3.0)
+    @test last(f).time == 1.0
+    @test last(f).data == 3.0
+    @test length(f) == 2
+
+    @testset "interpolation in time direction" begin
+        @test isa(f(0.0), DCTI) # converts to time-invariant after time interpolation
+        @test isapprox(f(-1.0), 2.0)
+        @test isapprox(f(0.0), 2.0)
+        @test isapprox(f(0.5), 2.5)
+        @test isapprox(f(1.0), 3.0)
+        @test isapprox(f(2.0), 3.0)
+    end
+
+    # create several time steps at once
+    f = DCTV(0.0 => 1.0, 1.0 => 2.0)
+    @test isapprox(f(0.5), 1.5)
+
+end
+
+@testset "discrete, variable, time-variant field" begin
+    f = Field(0.0 => [1.0, 2.0])
+    @test isa(f, DVTV)
+    @test last(f).time == 0.0
+    @test last(f).data == [1.0, 2.0]
+    update!(f, 0.0 => [2.0, 3.0])
+    @test last(f).time == 0.0
+    @test last(f).data == [2.0, 3.0]
+    @test length(f) == 1
+    update!(f, 1.0 => [3.0, 4.0])
+    @test last(f).time == 1.0
+    @test last(f).data == [3.0, 4.0]
+    @test length(f) == 2
+    
+    @testset "interpolation in time direction" begin
+        @test isa(f(0.0), DVTI) # converts to time-invariant after time interpolation
+        @test isapprox(f(-1.0), [2.0, 3.0])
+        @test isapprox(f(0.0), [2.0, 3.0])
+        @test isapprox(f(0.5), [2.5, 3.5])
+        @test isapprox(f(1.0), [3.0, 4.0])
+        @test isapprox(f(2.0), [3.0, 4.0])
+    end
+
+    # create several time steps at once
+    f = DVTV(0.0 => [1.0, 2.0], 1.0 => [2.0, 3.0])
+    @test isapprox(f(0.5), [1.5, 2.5])
+end
+
+@testset "continuous, constant, time-invariant field" begin
+    f = Field(() -> 2.0)
+    @test isapprox(f([1.0], 2.0), 2.0)
+
+end
+
+@testset "continuous, constant, time variant field" begin
+    f = Field((time::Float64) -> 2.0*time)
+    @test isapprox(f([1.0], 2.0), 4.0)
+
+end
+
+@testset "continuous, variable, time invariant field" begin
+    f = Field((xi::Vector) -> sum(xi))
+    @test isapprox(f([1.0, 2.0], 2.0), 3.0)
+end
+
+@testset "continuous, variable, time variant field" begin
+    f = Field((xi::Vector, t::Float64) -> xi[1]*t)
+    @test isapprox(f([1.0], 2.0), 2.0)
+end
+
+@testset "unknown function argument for continuous field" begin
+    @test_throws ErrorException Field((a, b, c) -> a*b*c)
+end
+
+@testset "dictionary fields" begin
+    f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0])
+    f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0])
+    f = Field(0.0 => f1, 1.0 => f2)
+    @test isa(f, DVTV)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f(1.0)[2], [1.0, 1.0])
+
+    f = Field(0.0 => f1)
+    update!(f, 1.0 => f2)
+    @test isa(f, DVTV)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f(1.0)[2], [1.0, 1.0])
+
+    f = Field(f1)
+    @test isapprox(f(0.0)[1], [0.0, 0.0])
+    @test isapprox(f[1], [0.0, 0.0])
+
+    f = Field(f1)
+    @test isa(f, DVTI)
+end
+
+@testset "test interpolation of discrete constant time-variant field" begin
+    f = DCTV(0.0 => 0.0, 1.0 => 1.0)
+    @test isapprox(f(-1.0), DCTI(0.0))
+    @test isapprox(f( 0.0), DCTI(0.0))
+    @test isapprox(f( 0.3), DCTI(0.3))
+    @test isapprox(f( 0.5), DCTI(0.5))
+    @test isapprox(f( 0.9), DCTI(0.9))
+    @test isapprox(f( 1.0), DCTI(1.0))
+    @test isapprox(f( 1.5), DCTI(1.0))
+
+    f2 = DCTV(0.0 => 0.0, 0.25 => -0.1, 0.50 => -0.1)
+    @test isapprox(f2(0.0), DCTI(0.0))
+    @test isapprox(f2(0.25), DCTI(-0.1))
+    @test isapprox(f2(0.50), DCTI(-0.1))
+    @test isapprox(f2(0.35), DCTI(-0.1))
+end
+
+@testset "default fields" begin
+    @test DCTI().data == nothing
+    @test Field().data == nothing
+    @test DCTI(1) == DCTI(1)
+    @test DCTI(1) == 1
+    @test length(DCTI(1)) == 1
+    @test length(DVTV()) == 0
+end
diff --git a/test/test_integration_points.jl b/test/test_integration_points.jl
new file mode 100644
index 0000000..679ed08
--- /dev/null
+++ b/test/test_integration_points.jl
@@ -0,0 +1,16 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using FEMBase: IP
+using Base.Test
+
+@testset "test integration point" begin
+    a = sqrt(1.0/3.0)
+    ip = IP(1, 1.0, (a, -a))
+    strain = [1.0 2.0; 3.0 4.0]
+    update!(ip, "strain", 0.0 => strain)
+    @test isapprox(ip("strain", 0.0), strain)
+    @test isapprox(ip("strain"), strain)
+end
+
diff --git a/test/test_sparse.jl b/test/test_sparse.jl
new file mode 100644
index 0000000..e9063c6
--- /dev/null
+++ b/test/test_sparse.jl
@@ -0,0 +1,49 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE
+
+using FEMBase
+using FEMBase: SparseMatrixCOO, SparseVectorCOO, optimize!,
+               resize_sparse, resize_sparsevec
+using Base.Test
+
+@testset "Add to SparseMatrixCOO" begin
+    A = SparseMatrixCOO()
+    A2 = reshape(collect(1:9), 3, 3)
+    add!(A, sparse(A2))
+    @test isapprox(full(A), full(A2))
+end
+
+@testset "Add to SparseVectorCOO" begin
+    b = SparseVectorCOO()
+    b2 = collect(1:3)
+    add!(b, sparse(b2))
+    @test isapprox(full(b), full(b2))
+end
+
+@testset "Failure to add data to sparse vector due dimensino mismatch" begin
+    b = SparseVectorCOO()
+    @test_throws ErrorException add!(b, [1, 2], [1.0, 2.0, 3.0])
+end
+
+@testset "Test combining of SparseMatrixCOO" begin
+    k = convert(Matrix{Float64}, reshape(collect(1:9), 3, 3))
+    dofs1 = [1, 2, 3]
+    dofs2 = [2, 3, 4]
+    A = SparseMatrixCOO()
+    add!(A, dofs1, dofs1, k)
+    add!(A, dofs2, dofs2, k)
+    A1 = full(A)
+    optimize!(A)
+    A2 = full(A)
+    @test isapprox(A1, A2)
+end
+
+@testset "resize of sparse matrix and sparse vector" begin
+    A = sparse(rand(3, 3))
+    B = resize_sparse(A, 4, 4)
+    @test size(B) == (4, 4)
+    a = sparse(rand(3))
+    b = resize_sparsevec(a, 4)
+    @test size(b) == (4, )
+end
+