diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb
index 8a1c2d2..b27ccad 100644
--- a/notebooks/2015-06-25-elasticity-solver-example.ipynb
+++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb
@@ -13,18 +13,36 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:24\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "  likely near /home/jukka/.julia/v0.4/Logging/src/Logging.jl:40\n"
+     ]
+    },
     {
      "data": {
       "text/plain": [
        "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)"
       ]
      },
-     "execution_count": 2,
+     "execution_count": 1,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -40,7 +58,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -51,7 +69,7 @@
        "CPS4"
       ]
      },
-     "execution_count": 3,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -89,7 +107,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 40,
    "metadata": {
     "collapsed": false
    },
@@ -97,11 +115,18 @@
     {
      "data": {
       "text/plain": [
-       "1x8 Array{Float64,2}:\n",
-       " -0.0  -0.0  -0.0  -0.0  -0.0  -0.0  -0.0  -0.0"
+       "8-element Array{Float64,1}:\n",
+       " -189.267      \n",
+       "   26.1073     \n",
+       "    5.7554e-13 \n",
+       "    2.6894e-12 \n",
+       "   -8.73968e-13\n",
+       "   20.0        \n",
+       "  189.267      \n",
+       "  -46.1073     "
       ]
      },
-     "execution_count": 5,
+     "execution_count": 40,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -125,22 +150,43 @@
     "    lambda = young*poisson/((1+poisson)*(1-2*poisson))\n",
     "    lambda = 2*lambda*mu/(lambda + 2*mu)\n",
     "\n",
-    "    invJ = inv(dN*X)\n",
-    "    gradw = dN*invJ\n",
+    "    gradw = dN*inv(dN*X)\n",
+    "    grad(u) = dN*u*inv(dN*(X+u))\n",
+    "    Grad(u) = dN*u*inv(dN*X)\n",
     "\n",
     "    function calc_R(data::Vector)\n",
     "        # create new field similar to displacement and fill it with data\n",
     "        f = similar(u, data)\n",
     "        # kinematics\n",
-    "        gradu = dN*f*invJ\n",
-    "        F = I + gradu  # deformation gradient\n",
-    "        E = 1/2*(F'*F - I)  # strain\n",
+    "        F = I + Grad(f)  # deformation gradient\n",
+    "        C = F'*F\n",
+    "        E = 1/2*(C - I)  # strain\n",
     "        S = 2*mu*E + lambda*trace(E)*I  # stress\n",
-    "        return (F*S*gradw')[:]\n",
+    "        P = F*S\n",
+    "        return (P*gradw')[:]\n",
     "    end\n",
     "\n",
-    "    R = calc_R(u[:])\n",
-    "    Kt = ForwardDiff.jacobian(calc_R, u[:])\n",
+    "    function calc_Wint(data::Vector)\n",
+    "        # create new field similar to displacement and fill it with data\n",
+    "        δu = similar(u, data)\n",
+    "        # kinematics\n",
+    "        F = I + Grad(δu)  # deformation gradient\n",
+    "        C = F'*F\n",
+    "        E = 1/2*(C - I)  # strain\n",
+    "        #E = 1/2*(Grad(δu)' + Grad(δu) + Grad(δu)'*Grad(δu))\n",
+    "        S = 2*mu*E + lambda*trace(E)*I  # stress\n",
+    "        δF = Grad(δu)\n",
+    "        #δɛ = 1/2*(δF' + δF)\n",
+    "        #δE = F'*δɛ*F\n",
+    "        #δE = 1/2*(F'*δF + δF'*F)\n",
+    "        #δE = 1/2*(Grad(δu)' + Grad(δu) + Grad(δu)'*Grad(δu))\n",
+    "        δE = 1/2*E  # ...???? energy 2x\n",
+    "\n",
+    "        return trace(S*δE')\n",
+    "    end\n",
+    "\n",
+    "    R = ForwardDiff.gradient(calc_Wint, u[:])\n",
+    "    Kt = ForwardDiff.hessian(calc_Wint, u[:])\n",
     "    return Kt, -R\n",
     "end\n",
     "\n",
@@ -153,42 +199,13 @@
     "\n",
     "JuliaFEM.has_lhs(equation::CPS4) = true\n",
     "JuliaFEM.has_rhs(equation::CPS4) = true\n",
-    "JuliaFEM.integrate_rhs(equation, 0.0)'"
+    "\n",
+    "JuliaFEM.integrate_rhs(equation, Inf)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "8x8 Array{Float64,2}:\n",
-       "  123.2    15.0    56.8    -3.0   -61.6   -15.0  -118.4     3.0\n",
-       "   15.0   321.2     3.0   158.8   -15.0  -160.6    -3.0  -319.4\n",
-       "   56.8     3.0   123.2   -15.0  -118.4    -3.0   -61.6    15.0\n",
-       "   -3.0   158.8   -15.0   321.2     3.0  -319.4    15.0  -160.6\n",
-       "  -61.6   -15.0  -118.4     3.0   123.2    15.0    56.8    -3.0\n",
-       "  -15.0  -160.6    -3.0  -319.4    15.0   321.2     3.0   158.8\n",
-       " -118.4    -3.0   -61.6    15.0    56.8     3.0   123.2   -15.0\n",
-       "    3.0  -319.4    15.0  -160.6    -3.0   158.8   -15.0   321.2"
-      ]
-     },
-     "execution_count": 6,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "JuliaFEM.integrate_lhs(equation, 0.0)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 41,
    "metadata": {
     "collapsed": false
    },
@@ -197,18 +214,19 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "increment 1, norm = 3.0900221367289094\n"
+      "increment  "
      ]
     },
     {
      "data": {
       "text/plain": [
-       "2-element Array{Float64,1}:\n",
-       " -0.399145\n",
-       " -2.17799 "
+       "Success :: (line:-1) :: fact was true\n",
+       "  Expression: (interpolate(element,\"displacement\",[1.0,1.0],Inf))[2] --> roughly(-4.15546385452579)\n",
+       "    Expected: -4.15546385452579\n",
+       "    Occurred: -4.155463854525783"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 41,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -216,23 +234,21 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "increment 2, norm = 0.32121316021534246\n",
-      "increment 3, norm = 0.040431781940025155\n",
-      "increment 4, norm = 0.0009291101052125905\n",
-      "increment 5, norm = 1.5638899099729665e-7\n",
-      "increment 6, norm = 1.0201033584846351e-14\n",
-      "increment 7, norm = 1.3679680786658739e-15\n",
-      "increment 8, norm = 2.463346396682122e-15\n",
-      "increment 9, norm = 9.485831749489368e-16\n",
-      "increment 10, norm = 1.3372608112297287e-15\n"
+      "1, norm =  5.77653, Wint =    0.000, Wext =   81.525, |Wint-Wext| = 81.52516\n",
+      "increment  2, norm =  0.99988, Wint =  173.895, Wext =   79.037, |Wint-Wext| = 94.85823\n",
+      "increment  3, norm =  0.28354, Wint =   87.084, Wext =   82.175, |Wint-Wext| =  4.90953\n",
+      "increment  4, norm =  0.07071, Wint =   82.548, Wext =   83.100, |Wint-Wext| =  0.55141\n",
+      "increment  5, norm =  0.00082, Wint =   83.116, Wext =   83.109, |Wint-Wext| =  0.00644\n",
+      "increment  6, norm =  0.00000, Wint =   83.109, Wext =   83.109, |Wint-Wext| =  0.00000\n",
+      "83.10927709051566\n"
      ]
     }
    ],
    "source": [
     "element = Quad4([1, 2, 3, 4])\n",
-    "push!(element, FieldSet(\"young\", [Field(0.0, 90.0)]))\n",
+    "push!(element, FieldSet(\"young\", [Field(0.0, 500.0)]))\n",
     "push!(element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])]))\n",
-    "push!(element, FieldSet(\"poisson\", [Field(0.0, 0.25)]))\n",
+    "push!(element, FieldSet(\"poisson\", [Field(0.0, 0.3)]))\n",
     "equation = CPS4(element)\n",
     "\n",
     "u0 = Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n",
@@ -241,8 +257,8 @@
     "du = zeros(8)\n",
     "fd = [3, 4, 5, 6]\n",
     "f = zeros(8)\n",
-    "f[6] = -2.0\n",
-    "for i=1:10\n",
+    "f[6] = -20.0\n",
+    "for i=1:6\n",
     "    A = JuliaFEM.integrate_lhs(equation, 1.0)\n",
     "    b = JuliaFEM.integrate_rhs(equation, 1.0)\n",
     "    du[fd] = A[fd,fd] \\ (b[fd]+f[fd])\n",
@@ -251,9 +267,14 @@
     "    new_field.time = 1.0\n",
     "    new_field.increment = i\n",
     "    push!(element[\"displacement\"], new_field)\n",
-    "    println(\"increment $i, norm = $(norm(du))\")\n",
+    "    Wint = (-b'*u)[1]\n",
+    "    Wext = (f'*u)[1]\n",
+    "    W = abs(Wint-Wext)\n",
+    "    @printf(\"increment %2d, norm = %8.5f, Wint = %8.3f, Wext = %8.3f, |Wint-Wext| = %8.5f\\n\", i, norm(du), Wint, Wext, W)\n",
     "end\n",
-    "interpolate(element, \"displacement\", [1.0, -1.0], 1.0)"
+    "println((f'*u)[1])\n",
+    "@fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-4.15546385452579E+00)\n",
+    "#interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2], (f'*u)[1]"
    ]
   },
   {
@@ -265,7 +286,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 170,
    "metadata": {
     "collapsed": false
    },
@@ -273,11 +294,10 @@
     {
      "data": {
       "text/plain": [
-       "1x8 Array{Float64,2}:\n",
-       " -0.0  -0.0  -0.0  -0.0  -0.0  -0.0  -0.0  -0.0"
+       "get_rhs (generic function with 6 methods)"
       ]
      },
-     "execution_count": 7,
+     "execution_count": 170,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -342,6 +362,9 @@
     "\n",
     "    R = fint\n",
     "    Kt = K_L + K_NL\n",
+    "\n",
+    "    print('.')\n",
+    "    \n",
     "    return Kt, -R\n",
     "end\n",
     "\n",
@@ -350,14 +373,56 @@
     "end\n",
     "function JuliaFEM.get_rhs(equation::CPS4, ip, time)\n",
     "    get_lhs_and_rhs(equation, ip, time)[2]\n",
-    "end\n",
-    "\n",
-    "JuliaFEM.integrate_rhs(equation, 0.0)'"
+    "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 172,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      ".."
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "8x8 Array{Float64,2}:\n",
+       "  119.461      13.486      61.8587   …   -26.0246  -113.693       2.80625\n",
+       "   13.486     309.576      15.6443      -148.518     -3.10973  -316.14   \n",
+       "   61.8587     15.6443    130.663        -53.2776   -60.0524      1.57406\n",
+       "    9.73233   155.082      36.0593      -303.877      1.51267  -163.427  \n",
+       "  -67.6267    -26.0206   -132.47          63.0312    56.539      10.2937 \n",
+       "  -26.0246   -148.518     -53.2776   …   299.439     16.271     152.956  \n",
+       " -113.693      -3.10973   -60.0524        16.271    117.207     -14.674  \n",
+       "    2.80625  -316.14        1.57406      152.956    -14.674     326.611  "
+      ]
+     },
+     "execution_count": 172,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      ".."
+     ]
+    }
+   ],
+   "source": [
+    "JuliaFEM.integrate_lhs(equation, Inf)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 169,
    "metadata": {
     "collapsed": false
    },
@@ -365,24 +430,16 @@
     {
      "data": {
       "text/plain": [
-       "8x8 Array{Float64,2}:\n",
-       "  123.2    15.0    56.8    -3.0   -61.6   -15.0  -118.4     3.0\n",
-       "   15.0   321.2     3.0   158.8   -15.0  -160.6    -3.0  -319.4\n",
-       "   56.8     3.0   123.2   -15.0  -118.4    -3.0   -61.6    15.0\n",
-       "   -3.0   158.8   -15.0   321.2     3.0  -319.4    15.0  -160.6\n",
-       "  -61.6   -15.0  -118.4     3.0   123.2    15.0    56.8    -3.0\n",
-       "  -15.0  -160.6    -3.0  -319.4    15.0   321.2     3.0   158.8\n",
-       " -118.4    -3.0   -61.6    15.0    56.8     3.0   123.2   -15.0\n",
-       "    3.0  -319.4    15.0  -160.6    -3.0   158.8   -15.0   321.2"
+       "2.927103942720631"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 169,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "JuliaFEM.integrate_lhs(equation, 0.0)"
+    "(1/2*u'*JuliaFEM.integrate_lhs(equation, 0.0)*u)[1]"
    ]
   },
   {
diff --git a/src/types.jl b/src/types.jl
index 1a2a93b..8e080a7 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -13,7 +13,7 @@ type Field{T}
 end
 """ Initialize field. """
 function Field(time, values)
-    Field(time, 1, values)
+    Field(time, 0, values)
 end
 """ Get length of a field (number of basis functions in practice). """
 function Base.length(f::Field)