diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 7f4573e..07cfebf 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -7,7 +7,8 @@ VERSION < v"0.4-" && using Docile
 using Lexicon
 
 include("types.jl") # type definitions
-include("interpolation.jl")
+include("math.jl") # basic mathematical operations
+
 include("elasticity_solver.jl")
 include("xdmf.jl")
 include("abaqus_reader.jl")
diff --git a/src/interpolation.jl b/src/interpolation.jl
deleted file mode 100644
index e0441ac..0000000
--- a/src/interpolation.jl
+++ /dev/null
@@ -1,53 +0,0 @@
-# This file is a part of JuliaFEM.
-# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
-
-"""
-Interpolate field variable using basis functions f for point ip.
-This function tries to be as general as possible and allows interpolating
-lot of different fields.
-
-Parameters
-----------
-field :: Array{Number, dim}
-  Field variable
-basis :: Function
-  Basis functions
-ip :: Array{Number, 1}
-  Point to interpolate
-"""
-function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip)
-    result = dot(field, basis(ip))
-    return result
-end
-function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
-    m, n = size(field)
-    bip = basis(ip)
-    tmp = size(bip)
-    if length(tmp) == 1
-        ndim = 1
-        nnodes = tmp[1]
-    else
-        ndim, nnodes = size(bip)
-    end
-    if ndim == 1
-        if n == nnodes
-            result = field * bip
-        elseif m == nnodes
-            result = field' * bip
-        end
-    else
-      if n == nnodes
-        result = bip' * field
-    elseif m == nnodes
-      result = bip' * field'
-    end
-    end
-    if length(result) == 1
-        result = result[1]
-    end
-    return result
-end
-function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
-    return interpolate(e.attributes[field], derivative ? e.dbasis : e.basis, x)
-end
-
diff --git a/src/math.jl b/src/math.jl
new file mode 100644
index 0000000..c4c45e6
--- /dev/null
+++ b/src/math.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
+
+using ForwardDiff
+
+"""
+Interpolate field variable using basis functions f for point ip.
+This function tries to be as general as possible and allows interpolating
+lot of different fields.
+
+Parameters
+----------
+field :: Array{Number, dim}
+  Field variable
+basis :: Function
+  Basis functions
+ip :: Array{Number, 1}
+  Point to interpolate
+"""
+function interpolate(field::Float64, basis::Function, ip::Array{Float64,1})
+    # dummy function, unable to interpolate scalar value!
+    return field
+end
+function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip)
+    result = dot(field, basis(ip))
+    return result
+end
+function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
+    m, n = size(field)
+    bip = basis(ip)
+    tmp = size(bip)
+    if length(tmp) == 1
+        ndim = 1
+        nnodes = tmp[1]
+    else
+        ndim, nnodes = size(bip)
+    end
+    if ndim == 1
+        if n == nnodes
+            result = field * bip
+        elseif m == nnodes
+            result = field' * bip
+        end
+    else
+      if n == nnodes
+        result = bip' * field
+    elseif m == nnodes
+      result = bip' * field'
+    end
+    end
+    if length(result) == 1
+        result = result[1]
+    end
+    return result
+end
+function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
+    return interpolate(e.attributes[field], derivative ? e.dbasis : e.basis, x)
+end
+
+
+function get_basis(el::Element, xi)
+    return el.basis(xi)
+end
+
+function get_dbasisdX(el::Element, xi)
+    J = interpolate(el, "coordinates", xi; derivative=true)
+    dbasisdX = el.dbasis(xi)*inv(J)
+    return dbasisdX
+end
+
+"""
+Linearize function f w.r.t some given field, i.e. calculate dR/du
+
+Parameters
+----------
+f::Function
+    (possibly) nonlinear function to linearize
+field::ASCIIString
+    field variable
+"""
+function linearize(f::Function, field::ASCIIString)
+    function jacobian(el::Element, xi)
+        dim, nnodes = size(el.attributes[field])
+        function helper!(x, y)
+            orig = copy(el.attributes[field])
+            el.attributes[field] = reshape(x, dim, nnodes)
+            y[:] = f(el, xi)
+            el.attributes[field] = copy(orig)
+        end
+        jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
+        return jac(el.attributes[field][:])
+    end
+    return jacobian
+end
+
+
+"""
+Integrate f over element using Gaussian quadrature rules.
+
+Parameters
+----------
+el::Element
+    well defined element
+f::Function
+    Function to integrate
+target::ASCIIString
+    Where to save result (el.attributes[target])
+"""
+function integrate!(el::Element, f::Function, target::ASCIIString)
+    # set target to zero
+    el.attributes[target][:] = 0.0
+    for m = 1:length(el.iweights)
+        w = el.iweights[m]
+        xi = el.ipoints[:, m]
+        J = interpolate(el, "coordinates", xi; derivative=true)
+        el.attributes[target] += w*f(el, xi)*det(J)
+    end
+end
+
+
diff --git a/src/types.jl b/src/types.jl
index 6ddedfd..05d38e9 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -4,7 +4,7 @@
 
 type Element
     id :: Int
-    element_type :: Int
+#   element_type :: Int
     node_ids :: Array{Int, 1}
     basis :: Function
     dbasis :: Function
diff --git a/test/test_interpolation.jl b/test/test_math.jl
similarity index 100%
rename from test/test_interpolation.jl
rename to test/test_math.jl