diff --git a/src/common.jl b/src/common.jl
deleted file mode 100644
index e49c6c3..0000000
--- a/src/common.jl
+++ /dev/null
@@ -1,3 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
diff --git a/src/core.jl b/src/core.jl
deleted file mode 100644
index eccdca1..0000000
--- a/src/core.jl
+++ /dev/null
@@ -1,2 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
diff --git a/src/equations.jl b/src/equations.jl
deleted file mode 100644
index 0851341..0000000
--- a/src/equations.jl
+++ /dev/null
@@ -1,132 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-function get_mass_matrix
-end
-
-function get_stiffness_matrix
-end
-
-function get_force_vector
-end
-
-function get_potential_energy
-end
-
-function get_residual_vector
-end
-
-function has_mass_matrix(problem::Problem, element::Element)
-    default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
-    return method_exists(get_mass_matrix, default_args)
-end
-
-function has_stiffness_matrix(problem::Problem, element::Element)
-    default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
-    return method_exists(get_stiffness_matrix, default_args)
-end
-
-function has_force_vector(problem::Problem, element::Element)
-    default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
-    return method_exists(get_force_vector, default_args)
-end
-
-function has_potential_energy(problem::Problem, element::Element)
-    default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
-    return method_exists(get_potential_energy, default_args)
-end
-
-function has_residual_vector(problem::Problem, element::Element)
-    default_args = Tuple{typeof(problem), typeof(element), IntegrationPoint, Float64}
-    return method_exists(get_residual_vector, default_args)
-end
-
-
-""" Assemble element. """
-function assemble!(assembly::Assembly, problem::Problem, element::Element, time::Number)
-
-    gdofs = get_gdofs(element, problem.dim)
-    unknown_field_name = get_unknown_field_name(problem)
-
-    # 1. if equations are defined we just integrate them, without caring how they are done
-    if has_mass_matrix(problem, element) || has_stiffness_matrix(problem, element) || has_force_vector(problem, element)
-        for ip in get_integration_points(element)
-            w = ip.weight*det(J)
-            if has_mass_matrix(element)
-                add!(assembly.mass_matrix, gdofs, gdofs, w*get_mass_matrix(problem, element, ip, time))
-            end
-            if has_stiffness_matrix(element)
-                add!(assembly.stiffness_matrix, gdofs, gdofs, w*get_stiffness_matrix(problem, element, ip, time))
-            end
-            if has_force_vector(element)
-                add!(assembly.force_vector, gdofs, w*get_force_vector(problem, element, ip, time))
-            end
-        end
-        # external loads -- if any nodal loads is defined add to force vector
-        if haskey(element, "$unknown_field_name nodal load")
-            add!(assembly.force_vector, gdofs, vec(element["$unknown_field_name nodal load"](time)))
-        end
-    end
-
-    # 2. energy form -- user has defined potential energy W -> min!
-    if has_potential_energy(problem, element) && haskey(element, unknown_field_name)
-        field = element[unknown_field_name](time)
-
-        """ Wrapper for potential energy for ForwardDiff. """
-        function calc_W(data::Vector)
-            W = 0.0
-            df = similar(field, data)
-            # integrate potential energy
-            for ip in get_integration_points(element)
-                dw = get_potential_energy(problem, element, ip, time; variation=df)
-                W += ip.weight*dw
-            end
-            # external energy -- if any nodal loads is defined, decrease from potential energy
-            if haskey(element, "$unknown_field_name nodal load")
-                P = element["$unknown_field_name nodal load"](time)
-                W -= dot(vec(P), vec(df))
-            end
-            return W[1]
-        end
-
-        hessian, allresults = ForwardDiff.hessian(calc_W, vec(field), AllResults, cache=autodiffcache)
-        add!(assembly.stiffness_matrix, gdofs, gdofs, hessian)
-        add!(assembly.force_vector, gdofs, -ForwardDiff.gradient(allresults))
-    end
-
-    # 3. virtual work -- user has defined some residual r = p - f = 0
-    if has_residual_vector(problem, element) && haskey(element, unknown_field_name)
-
-        field = DVTI(last(element[unknown_field_name]).data)
-
-        """ Wrapper for virtual work for ForwardDiff. """
-        function calc_R(data::Vector)
-            R = zeros(length(data))
-            df = similar(field, data)
-            gauss_fields = IntegrationPoint[]
-            # integrate residual vector
-            for ip in get_integration_points(element)
-                dr = get_residual_vector(problem, element, ip, time; variation=df)
-                R += ip.weight*dr
-                if ip.changed
-                    push!(gauss_fields, ip)
-                end
-            end
-            # external loads -- if any nodal loads is defined, decrease from residual
-            if haskey(element, "$unknown_field_name nodal load")
-                R -= vec(element["$unknown_field_name nodal load"](time))
-            end
-            #info("return = $R")
-            if length(gauss_fields) != 0
-                update_gauss_fields!(element, gauss_fields, time)
-            end
-            return R
-        end
-
-        jacobian, allresults = ForwardDiff.jacobian(calc_R, vec(field), AllResults, cache=autodiffcache)
-        residual_vector = -ForwardDiff.value(allresults)
-        add!(assembly.stiffness_matrix, gdofs, gdofs, jacobian)
-        add!(assembly.force_vector, gdofs, residual_vector)
-    end
-end
-
diff --git a/src/symbolic.jl b/src/symbolic.jl
deleted file mode 100644
index 4822b73..0000000
--- a/src/symbolic.jl
+++ /dev/null
@@ -1,116 +0,0 @@
-# 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
-
-immutable SymbolicField <: AbstractField
-    name :: ASCIIString
-end
-
-immutable Expression
-    expr :: Expr
-end
-
-function Expression(expression::ASCIIString)
-    return Expression(parse(expression))
-end
-
-function Base.convert(::Type{Field}, name::ASCIIString)
-    return SymbolicField(name)
-end
-
-function Base.convert(::Type{Symbol}, field::SymbolicField)
-    return Symbol(field.name)
-end
-
-function Base.(:*)(n::Number, field::SymbolicField)
-    expr = Expr(:call, :*, n, Symbol(field.name))
-    return Expression(expr)
-end
-
-function grad(field::SymbolicField)
-    expr = Expr(:call, :grad, Symbol(field.name))
-    return Expression(expr)
-end
-
-function diff(field::SymbolicField)
-    expr = Expr(:call, :diff, Symbol(field.name))
-    return Expression(expr)
-end
-
-function Base.call(basis::Basis, expr::Expression, fieldset::FieldSet,
-                   xi::Vector, time::Number)
-    return replace(expr.expr, basis, fieldset, xi, time)
-end
-
-# Unbelievable code. I don't know why or how this works.
-""" Replace symbolic fields with real arrays. """
-function Base.replace(expression::Expr, basis::Basis, fieldset::FieldSet,
-                      xi::Vector, time::Number, data=Dict())
-
-    info("expression = $expression")
-    if expression.head == symbol("'")
-        info("transpose")
-        expr = Expr(:call, :transpose, expression.args...)
-        return replace(expr, basis, fieldset, xi, time, data)
-    end
-    operator = expression.args[1]
-
-    if operator == :diff
-        info("inside diff operator")
-        if !haskey(data, expression)
-            data[expression] = fieldset[string(expression.args[2])](time, Val{:diff})
-        end
-        info("data = $(data[expression])")
-        #return basis(data[expression], xi)
-        return data[expression]
-    end
-    
-    if operator == :grad
-        info("inside gradient operator")
-        info("gradient args: $(expression.args)")
-        field_name = expression.args[2]
-        if isa(field_name, Expr)
-            info("expression inside gradient")
-        end
-        #if startswith(string(field_name), "diff")
-        if !haskey(data, field_name)
-            info("grad: evaluate field $field_name")
-            if isa(field_name, Symbol)
-                data[field_name] = fieldset[string(field_name)](time)
-            elseif isa(field_name, Expr) && (field_name.args[1] == :diff)
-                info("taking time derivative of $(field_name.args[2])")
-                #data[field_name] = fieldset[string(field_name.args[2])](time, Val{:diff})
-                data[field_name] = replace(field_name, basis, fieldset, xi, time, data)
-            end
-        end
-        if !haskey(data, :geometry)
-            info("evaluate geometry")
-            data[:geometry] = fieldset["geometry"](time)
-        end
-        #return Expr(:call, basis, data[:geometry], data[field_name], xi, Val{:gradient})
-        #return :(basis($(data[:geometry]), $(data[field_name]), $xi, Val{:gradient}))
-        info("evaluate gradient")
-        return basis(data[:geometry], data[field_name], xi, Val{:grad})
-    end
-
-    for i in 2:length(expression.args)
-        arg = expression.args[i]
-        if isa(arg, Expr)
-            expression.args[i] = replace(arg, basis, fieldset, xi, time, data)
-        end
-        if isa(arg, Symbol) && haskey(fieldset, string(arg))
-            if !haskey(data, arg)
-                info("evaluate field $arg")
-                data[arg] = fieldset[string(arg)](time)
-            end
-            #expression.args[i] = Expr(:call, basis, data[arg], xi)
-            #return :(basis($(data[arg]), $xi))
-            expression.args[i] = basis(data[arg], xi)
-        end
-    end
-
-    info("new expression: $expression")
-    return expression
-end
-
diff --git a/test/test_symbolic.jl b/test/test_symbolic.jl
deleted file mode 100644
index 04e8e52..0000000
--- a/test/test_symbolic.jl
+++ /dev/null
@@ -1,119 +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.Test
-
-
-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(basis, dbasis)
-end
-
-function get_fieldset()
-    X = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
-    T = Field(
-        (0.0, [0, 0, 0, 0]),
-        (1.0, [1, 2, 3, 4]))
-    u = Field(
-        (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
-        (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
-    # FIXME: this is not working
-    #fieldset = FieldSet("geometry" => X, "temperature" => T, "displacement" => u)
-    fieldset = FieldSet()
-    fieldset["geometry"] = X
-    fieldset["temperature"] = T
-    fieldset["displacement"] = u
-    return fieldset
-end
-
-function test_create_symbolic_field()
-    f = Field("temperature")
-    @test isa(f, Field)
-end
-
-function test_evaluate_symbolic_field()
-    T = Field("temperature")
-    expr = Symbol(T)
-    @test expr == :(temperature)
-end
-
-function test_simple_math()
-    T = Field("temperature")
-    eq = 1/2*T
-    @test eq.expr == :(0.5*temperature)
-end
-
-function test_evaluate_expression()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    T = Field("temperature")
-    expr = 1/2*T
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), 1/2*mean([1, 2, 3, 4]))
-end
-
-function test_evaluate_gradient_of_scalar_field()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    T = Field("temperature")
-    expr = grad(T)
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), [0.0 2.0])
-end
-
-function test_evaluate_gradient_of_vector_field()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    u = Field("displacement")
-    expr = grad(u)
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
-end
-
-function test_evaluate_strain_rate()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    u = Field("displacement")
-    #expr = grad(u)
-    expr = Expression("1/2*(grad(diff(displacement)) + grad(diff(displacement))')")
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), [1.5 0.75; 0.75 2.0])
-end
-
-function test_evaluate_grad_diff()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    u = Field("displacement")
-    #expr = grad(u)
-    expr = Expression("grad(diff(displacement))")
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
-end
-
-function test_grad_diff_simplification()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    u = Field("displacement")
-    expr1 = diff(grad(u))
-    expr2 = Expression("grad(diff(displacement))")
-    info("expr1 = $expr1")
-    info("expr2 = $expr2")
-    @test expr1 == expr2
-end
-
-function test_evaluate_time_derivative()
-    basis = get_basis()
-    fieldset = get_fieldset()
-    T = Field("temperature")
-    expr = diff(T)
-    result = basis(expr, fieldset, [0.0, 0.0], 1.0)
-    @test isapprox(eval(result), mean([1.0, 2.0, 3.0, 4.0]))
-end
-