diff --git a/docs/tutorials/2015-08-29-developing-juliafem.ipynb b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
index 27f8ee4..c2912ab 100644
--- a/docs/tutorials/2015-08-29-developing-juliafem.ipynb
+++ b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
@@ -8,24 +8,48 @@
     "\n",
     "Author(s): Jukka Aho\n",
     "\n",
-    "**Abstract**: Developer notes. In this notebook we give general guidelines how to implement things to JuliaFEM.\n",
+    "**Abstract**: Developer notes. In this notebook we give the general guidelines how to implement basic fundamentals to JuliaFEM. The notebook is rapidly updated and should always represent the newest \"style\" how to develop things to JuliaFEM. At this phase of project interfaces are changing rapidly. For this reason we aim to develop couple of simple functions which aim to test that elements, equations, solvers etc. work as expected, in the sense that their interface hasn't any design problems. (Yes, we do believe in testing and follow also other practices proven to be good, like continuous integration). The aim of this document is not to be the guide to \"other contributors\" but is more like an guide line for all contributors, including ourselves, how things should be done.\n",
     "\n",
-    "In short, we have\n",
+    "For demonstrational purposes we use the Poisson equation in our examples to demonstrate the main concepts of the package. Poisson equation has some analogies to solid mechanics. Consider for example $(EAu')' + q = 0$, which is nothing more than a Poisson equation defined in 1d.\n",
+    "\n",
+    "Because this is the second important document in this source repository, any questions or suggestions araising from this notebook are very welcome and desirable and should be adressed to the our issue log: https://github.com/JuliaFEM/JuliaFEM.jl/issues. Easiest way to make world better is to register GitHub (takes 1 min) and drop an issue. And the most important document? It's the \"how to contribute\" guide, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/CONTRIBUTING.rst.\n",
+    "\n",
+    "In this notebook briefly we have defined the following \"objects\":\n",
     "\n",
     "**Element**: object which holds basis functions and fields. One is able to interpolate fields using element. Introduction of new basis functions needs new elements. Typical elements: Lagrange elements, hierarchical elements, etc.\n",
     "\n",
-    "**Equation**: object which defines some field equation which needs to be solve. One element can have several equations but not vice versa. Typical equations: Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc.\n",
+    "**Equation**: object which defines some field equation, like Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc., which needs to be solve. This is the physics we are trying to solve by approximating the equations in element area by some interpolation function provided by element.\n",
     "\n",
     "**Problem**: object which maps equations to elements. For example, HeatProblem which maps Poisson equation to Lagrange elements or ElasticityProblem which maps elasticity equation to Lagrange elements.\n",
     "\n",
-    "**Solver**: object which takes one or more problems, solves them using some strategy (iterative methods, multigrid, direct methods, ...) and updates corresponding fields to elements.\n",
+    "**Solver**: object which takes one or more problems, solves them using some strategy (iterative methods, multigrid, direct methods, ...) and updates corresponding fields to elements. The is the section where any numerical crunching happens and is targeted to high-performance computing. \n",
     "\n",
-    "Moreover, we have (will have) **`test_element`**, **`test_equation`**, **`test_problem`** and **`test_solver`** which can be used to test that implementation has all necessary things defined."
+    "These concepts are quite well separated so that development process can focus just one of thing of interest. In this notebook we will give instructions to all of the sections describing the whole development process from element level to global assembly and solution.\n",
+    "\n",
+    "We have\n",
+    "- [x] **`test_element`**\n",
+    "- [ ] **`test_equation`**\n",
+    "- [ ] **`test_problem`** \n",
+    "- [ ] **`test_solver`**\n",
+    "\n",
+    "which can be used to test that implementation has all necessary things defined.\n",
+    "\n",
+    "\n",
+    "Poisson equation is used in this notebook to demonstrate all phases of development process. It's weak form is: find $u\\in\\mathcal{U}$ such that\n",
+    "\\begin{equation}\n",
+    "    \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x=\\int_{\\Omega}fv\\,\\mathrm{d}x+\\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s\\quad\\forall v\\in\\mathcal{V}.\n",
+    "\\end{equation}\n",
+    "\n",
+    "## Definitions / nomenclature\n",
+    "\n",
+    "- weighted residual form: differential equation is multiplied by a weight function and integrating it.\n",
+    "- weak form, \"principle of virtual work\": typically established by partial integration of a weighted residual form. $\\delta\\mathcal{W}_{\\mathrm{int}}=\\delta\\mathcal{W}_{\\mathrm{ext}}$.\n",
+    "- variational form, \"principle of minimum potential energy\": there exists some functional  or \"potential function\" $\\Pi$ we are minimizing"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -36,7 +60,7 @@
        "Logger(root,DEBUG,Base.PipeEndpoint(open, 0 bytes waiting),root)"
       ]
      },
-     "execution_count": 39,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -65,7 +89,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 4,
    "metadata": {
     "collapsed": false
    },
@@ -83,7 +107,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": 5,
    "metadata": {
     "collapsed": false
    },
@@ -105,7 +129,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 6,
    "metadata": {
     "collapsed": false
    },
@@ -116,7 +140,7 @@
        "MyQuad4"
       ]
      },
-     "execution_count": 42,
+     "execution_count": 6,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -147,7 +171,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 7,
    "metadata": {
     "collapsed": false
    },
@@ -158,7 +182,7 @@
        "get_element_dimension (generic function with 7 methods)"
       ]
      },
-     "execution_count": 43,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -177,7 +201,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 44,
+   "execution_count": 8,
    "metadata": {
     "collapsed": false
    },
@@ -186,7 +210,17 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "20-loka 15:49:17:INFO:root:Testing element MyQuad4\n"
+      "26-Oct 05:30:31:INFO:root:Testing element MyQuad4\n",
+      "26-Oct 05:30:31:INFO:root:number of basis functions in this element: 4\n",
+      "26-Oct 05:30:31:INFO:root:Initializing element\n",
+      "26-Oct 05:30:31:INFO:root:Element dimension: 2\n",
+      "26-Oct 05:30:31:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])\n",
+      "26-Oct 05:30:31:INFO:root:basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n",
+      "26-Oct 05:30:31:INFO:root:field val at [0.0,0.0]: 2.5\n",
+      "26-Oct 05:30:32:INFO:root:derivative of basis at [0.0,0.0]: [-0.5 0.5 0.5 -0.5\n",
+      " -0.5 -0.5 0.5 0.5]\n",
+      "26-Oct 05:30:32:INFO:root:field val at [0.0,0.0]: [0.0 2.0]\n",
+      "26-Oct 05:30:32:INFO:root:Element MyQuad4 passed tests.\n"
      ]
     }
    ],
@@ -204,7 +238,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 45,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -212,23 +246,27 @@
     {
      "data": {
       "text/plain": [
-       "JuliaFEM.FieldSet(symbol(\"heat coefficient\"),JuliaFEM.Field[JuliaFEM.Field{Int64}(0.0,1,2),JuliaFEM.Field{Int64}(1.0,1,3)])"
+       "Dict{Symbol,JuliaFEM.FieldSet} with 4 entries:\n",
+       "  symbol(\"heat coefficien… => JuliaFEM.FieldSet(symbol(\"heat coefficient\"),Juli…\n",
+       "  :geometry                => JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaF…\n",
+       "  :temperature             => JuliaFEM.FieldSet(:temperature,JuliaFEM.Field[Jul…\n",
+       "  :displacement            => JuliaFEM.FieldSet(:displacement,JuliaFEM.Field[Ju…"
       ]
      },
-     "execution_count": 45,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "using JuliaFEM: FieldSet, Field, interpolate, dinterpolate\n",
+    "using JuliaFEM: FieldSet, Field\n",
     "element = MyQuad4([1, 2, 3, 4])\n",
     "\n",
     "geometry_field = Field(0.0, Vector[])  # Create empty field at time t=0.0\n",
-    "push!(geometry_field, [ 0.0, 0.0, 0.0])  # push some values for field\n",
-    "push!(geometry_field, [10.0, 0.0, 0.0])\n",
-    "push!(geometry_field, [10.0, 1.0, 0.0])\n",
-    "push!(geometry_field, [ 0.0, 1.0, 0.0])\n",
+    "push!(geometry_field, [ 0.0, 0.0])  # push some values for field\n",
+    "push!(geometry_field, [ 1.0, 0.0])\n",
+    "push!(geometry_field, [ 1.0, 1.0])\n",
+    "push!(geometry_field, [ 0.0, 1.0])\n",
     "geometry_fieldset = FieldSet(\"geometry\")  # create fieldset \"geometry\"\n",
     "push!(geometry_fieldset, geometry_field)  # add field to fieldset\n",
     "push!(element, geometry_fieldset) # add fieldset to element\n",
@@ -238,15 +276,21 @@
     "push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
     "push!(element, temperature_fieldset)\n",
     "\n",
+    "displacement_fieldset = FieldSet(\"displacement\")\n",
+    "push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))\n",
+    "push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))\n",
+    "push!(element, displacement_fieldset)\n",
+    "\n",
     "heat_coefficient_fieldset = FieldSet(\"heat coefficient\")\n",
     "push!(heat_coefficient_fieldset, Field(0.0, 2))\n",
     "push!(heat_coefficient_fieldset, Field(1.0, 3))\n",
-    "push!(element, heat_coefficient_fieldset)"
+    "push!(element, heat_coefficient_fieldset)\n",
+    "element.fields"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 46,
+   "execution_count": 10,
    "metadata": {
     "collapsed": false
    },
@@ -254,48 +298,26 @@
     {
      "data": {
       "text/plain": [
-       "1.25"
-      ]
-     },
-     "execution_count": 46,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "# temperature at the middle poinf of the element, 1/4*(1+2+3+4) at t=0.5\n",
-    "interpolate(element, \"temperature\", [0.0, 0.0], 0.5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 47,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "3-element Array{Float64,1}:\n",
-       " 5.0\n",
+       "2-element Array{Float64,1}:\n",
        " 0.5\n",
-       " 0.0"
+       " 0.5"
       ]
      },
-     "execution_count": 47,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "# geometry midpoint of element\n",
-    "interpolate(element, \"geometry\", [0.0, 0.0], 0.0)"
+    "using JuliaFEM: get_basis, grad\n",
+    "basis = get_basis(element)\n",
+    "basis(\"geometry\", [0.0, 0.0], 0.0)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 11,
    "metadata": {
     "collapsed": false
    },
@@ -303,25 +325,24 @@
     {
      "data": {
       "text/plain": [
-       "3x2 Array{Float64,2}:\n",
-       " 5.0  0.0\n",
-       " 0.0  0.5\n",
-       " 0.0  0.0"
+       "1x2 Array{Float64,2}:\n",
+       " 0.0  1.0"
       ]
      },
-     "execution_count": 48,
+     "execution_count": 11,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "# interpolate derivatives works too\n",
-    "dinterpolate(element, \"geometry\", [0.0, 0.0], 0.0)"
+    "dbasis = grad(basis)\n",
+    "dbasis(\"temperature\", [0.0, 0.0], 0.5)  # temperature gradient at mid point of element at time t=0.5"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 49,
+   "execution_count": 12,
    "metadata": {
     "collapsed": false
    },
@@ -332,19 +353,19 @@
        "2.5"
       ]
      },
-     "execution_count": 49,
+     "execution_count": 12,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "# interpolating scalar -> scalar.\n",
-    "interpolate(element, \"heat coefficient\", [0.0, 0.0], 0.5)"
+    "basis(\"heat coefficient\", [0.0, 0.0], 0.5)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 50,
+   "execution_count": 13,
    "metadata": {
     "collapsed": false
    },
@@ -355,22 +376,120 @@
        "3"
       ]
      },
-     "execution_count": 50,
+     "execution_count": 13,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "interpolate(element, \"heat coefficient\", [0.0, 0.0], Inf)"
+    "# set time to Inf to get very last field value and -Inf to get first one.\n",
+    "basis(\"heat coefficient\", [0.0, 0.0], Inf)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2x2 Array{Float64,2}:\n",
+       " 0.132813   0.0703125\n",
+       " 0.0703125  0.0078125"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# some continuum mechanics\n",
+    "∇u = dbasis(\"displacement\", [0.0, 0.0], 1.0)\n",
+    "ɛ = 1/2*(∇u + ∇u')\n",
+    "Ω = 1/2*(∇u - ∇u')\n",
+    "X = basis(\"geometry\", [0.0, 0.0], 1.0)\n",
+    "F = I + ∇u\n",
+    "C = F'*F\n",
+    "E = 1/2*(F'*F - I)  # Green-Lagrange strain tensor"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "Basically one could already write local stiffness matrices:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "(\n",
+       "4x4 Array{Float64,2}:\n",
+       "  4.0  -1.0  -2.0  -1.0\n",
+       " -1.0   4.0  -1.0  -2.0\n",
+       " -2.0  -1.0   4.0  -1.0\n",
+       " -1.0  -2.0  -1.0   4.0,\n",
+       "\n",
+       "4x4 Array{Float64,2}:\n",
+       " 4.0  2.0  1.0  2.0\n",
+       " 2.0  4.0  2.0  1.0\n",
+       " 1.0  2.0  4.0  2.0\n",
+       " 2.0  1.0  2.0  4.0,\n",
+       "\n",
+       "4x1 Array{Float64,2}:\n",
+       " 1.0\n",
+       " 1.0\n",
+       " 1.0\n",
+       " 1.0)"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "weights = [1.0, 1.0, 1.0, 1.0]\n",
+    "integration_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]\n",
+    "\n",
+    "M = zeros(4, 4)\n",
+    "K = zeros(4, 4)\n",
+    "f = zeros(4)\n",
+    "k = 6\n",
+    "rho = 36\n",
+    "q = 4\n",
+    "\n",
+    "for (w, xi) in zip(weights, integration_points)\n",
+    "    N = get_basis(element)\n",
+    "    ∇N = grad(N)  # gradient is with respect to \"geometry\" field\n",
+    "    detJ = det(N)(xi)\n",
+    "    M += w*rho*N(xi)'*N(xi)*detJ  # mass matrix\n",
+    "    K += w*k*∇N(xi)'*∇N(xi)*detJ  # stiffness matrix\n",
+    "    f += w*N(xi)'*q*detJ  # force vector\n",
+    "end\n",
+    "K, M, f"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "but we have an higher level interface for that also.\n",
+    "\n",
     "### Summary of developing own elements\n",
     "\n",
-    "Element itself if not calculating anything but only stores fields and basis functions so that the fields can be interpolated. We will provide command `test_element` which will ensure that everything necessary is defined. While lot of things needs to be defined, by subclassing from `Element` most of these are already defined, thanks to multiple dispatch."
+    "Element itself if not calculating anything but only stores fields and basis functions so that the fields can be interpolated. We will provide command `test_element` which will ensure that everything necessary is defined. While a lot of things needs to be defined, by subclassing from `Element` most of these are already defined, thanks to multiple dispatch."
    ]
   },
   {
@@ -379,34 +498,29 @@
    "source": [
     "## Developing own equation\n",
     "\n",
-    "Let's consider a Laplace equation\n",
+    "Let's consider a Poisson equation\n",
     "\\begin{align}\n",
-    "\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\n",
-    "\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}}\n",
+    "-\\nabla \\cdot \\left( k \\nabla u\\right) &= f && \\text{on } \\Omega \\\\\n",
+    "\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}} \\\\\n",
+    "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
     "\\end{align}\n",
     "\n",
     "Weak form is, find $u\\in\\mathcal{U}$ such that\n",
     "\\begin{equation}\n",
-    "    \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Gamma_{\\mathrm{N}}}g v\\,\\mathrm{d}s \\quad \\forall v\\in\\mathcal{V}.\n",
+    "    \\int_{\\Omega}k \\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Omega}fv\\,\\mathrm{d}x + \\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s \\quad\\forall v\\in\\mathcal{V}.\n",
     "\\end{equation}\n",
     "\n",
     "Minimum requirements for equation: \n",
     "- subclass from Equation, if not want to implement from scratch\n",
-    "- it needs to have lhs and rhs functions\n",
     "- provide the name of the unknown field variable trying to solve\n",
     "- default constructor takes the element as input argument\n",
     "\n",
-    "Now we have function `test_equation`, which we can use to test that everything is working as expected. \n",
-    "\n",
-    "Again thanks to multiple dispatch, you are free to code your weak form however you want as long as it returns lhs and rhs sides for element dofs. This kind of freedom gives good opportunities to wrap e.g. Fortran code from some other projects. And again we have some suggestions ad following these ideas you get a lot of stuff for free. First we look the left hand side of the equation, that is,\n",
-    "\\begin{equation}\n",
-    "    \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x\n",
-    "\\end{equation}"
+    "This time we will have function `test_equation`, which we can use to test that everything is working as expected. "
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 51,
+   "execution_count": 16,
    "metadata": {
     "collapsed": false
    },
@@ -417,17 +531,15 @@
        "get_unknown_field_name (generic function with 4 methods)"
       ]
      },
-     "execution_count": 51,
+     "execution_count": 16,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "using JuliaFEM: Equation, IntegrationPoint, Quad4, get_unknown_field_name\n",
+    "using JuliaFEM: Equation, IntegrationPoint, Quad4, Seg2, get_unknown_field_name\n",
     "\n",
     "abstract Heat <: Equation\n",
-    "\n",
-    "# it's important to define for which fieldset to save unknown values when solving\n",
     "JuliaFEM.get_unknown_field_name(eq::Heat) = symbol(\"temperature\")"
    ]
   },
@@ -435,37 +547,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Our basic data type often looks something like this:"
+    "Basic data type needs to contain element, integration points and global degrees of freedom for global assembly:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 52,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "\"\"\"\n",
-    "Diffusive heat transfer for 4-node bilinear element.\n",
-    "\"\"\"\n",
-    "type DC2D4 <: Heat\n",
-    "    element :: Quad4\n",
-    "    integration_points :: Array{IntegrationPoint, 1}\n",
-    "    global_dofs :: Array{Int64, 1}\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We must provide default constructor which takes element as input argument:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 53,
+   "execution_count": 17,
    "metadata": {
     "collapsed": false
    },
@@ -473,15 +560,21 @@
     {
      "data": {
       "text/plain": [
-       "DC2D4"
+       "size (generic function with 64 methods)"
       ]
      },
-     "execution_count": 53,
+     "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
+    "\"\"\" Diffusive heat transfer for 4-node bilinear element. \"\"\"\n",
+    "type DC2D4 <: Heat\n",
+    "    element :: Quad4\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end\n",
     "function DC2D4(element::Quad4)\n",
     "    integration_points = [\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
@@ -490,19 +583,36 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
     "    push!(element, FieldSet(\"temperature\"))\n",
     "    DC2D4(element, integration_points, [])\n",
-    "end"
+    "end\n",
+    "Base.size(equation::DC2D4) = 4\n",
+    "\n",
+    "\"\"\" Diffusive heat transfer for 2-node linear segment. \"\"\"\n",
+    "type DC2D2 <: Heat\n",
+    "    element :: Seg2\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end\n",
+    "function DC2D2(element::Seg2)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint([0.0], 2.0)]\n",
+    "    push!(element, FieldSet(\"temperature\"))\n",
+    "    DC2D2(element, integration_points, [])\n",
+    "end\n",
+    "Base.size(equation::DC2D2) = 2"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Now the actual implementation for $\\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x$:"
+    "Now the actual implementation for $\\int_{\\Omega}k \\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Omega}fv\\,\\mathrm{d}x + \\int_{\\Gamma_{\\mathrm{N}}}gv\\,\\mathrm{d}s$.\n",
+    "\n",
+    "The low-level way is to write equations for local matrix representation directly:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 54,
+   "execution_count": 18,
    "metadata": {
     "collapsed": false
    },
@@ -510,39 +620,97 @@
     {
      "data": {
       "text/plain": [
-       "has_lhs (generic function with 4 methods)"
+       "has_force_vector (generic function with 3 methods)"
       ]
      },
-     "execution_count": 54,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "using JuliaFEM: get_element, get_dbasisdX, has_lhs, has_rhs\n",
+    "using JuliaFEM: get_basis, get_element\n",
     "\n",
-    "\"\"\"\n",
-    "Left hand side defined in integration point\n",
-    "\"\"\"\n",
-    "function JuliaFEM.get_lhs(eq::DC2D4, ip, t)\n",
-    "    el = get_element(eq)\n",
-    "    dNdX = get_dbasisdX(el, ip.xi, t)\n",
-    "    hc = interpolate(el, \"temperature thermal conductivity\", ip.xi, t)\n",
-    "    return dNdX*hc*dNdX'\n",
+    "\"\"\" Mass matrix for dynamical problems. Not used in this example. \"\"\"\n",
+    "function JuliaFEM.get_mass_matrix(equation::DC2D4, ip, time)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    ρ = basis(\"density\", ip, time)\n",
+    "    return ρ * basis(ip,time)'*basis(ip,time)\n",
     "end\n",
-    "JuliaFEM.has_lhs(eq::DC2D4) = true"
+    "\"\"\" Left hand side defined in integration point. \"\"\"\n",
+    "function JuliaFEM.get_stiffness_matrix(equation::DC2D4, ip, time)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)\n",
+    "    k = basis(\"temperature thermal conductivity\", ip, time)\n",
+    "    return k * dbasis(ip,time)'*dbasis(ip,time)\n",
+    "end\n",
+    "\"\"\" Right hand side defined in integration point. \"\"\"\n",
+    "function JuliaFEM.get_force_vector(equation::DC2D4, ip, time)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    f = basis(\"temperature load\", ip, time)\n",
+    "    return basis(ip,time)'*f\n",
+    "end\n",
+    "JuliaFEM.has_mass_matrix(equation::DC2D4) = true\n",
+    "JuliaFEM.has_stiffness_matrix(equation::DC2D4) = true\n",
+    "JuliaFEM.has_force_vector(equation::DC2D4) = true\n",
+    "\n",
+    "\"\"\" Right hand side defined in integration point. \"\"\"\n",
+    "function JuliaFEM.get_force_vector(equation::DC2D2, ip, time)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    g = basis(\"temperature flux\", ip, time)\n",
+    "    return basis(ip,time)'*g\n",
+    "end\n",
+    "JuliaFEM.has_force_vector(equation::DC2D2) = true"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "And that's it. If we want to play with this formulation, we must create element and assign this equation for it:"
+    "It might at start to feel a bit ackward to pass integration point and time all around as function arguments, so a little reasoning why this is essential might be needed.\n",
+    "\n",
+    "First of all, equations are written in integration or gauss quadrature points for a certain reason. The reason is that for some incremental constitutive models internal variables are needed, which are saved to FieldSets inside integration points. We can always access these variables through `ip.fields` in the same way we can access all other field variables defined inside element via `element.fields`.\n",
+    "\n",
+    "Secondly, starting point for our equation is that everything can be time-dependent. FieldSet is, like name suggests, a set of fields with different times and/or increments. By transferring time parameter inside equation makes it possible to access not only current field values but earlier values too. In this example field variables thermal conductivity $k$, load $f$ and flux $g$ are constant in time but this way they can depend easily from spatial or temporal dimension. Actually let's change the parameter $g$ such a way that it is linear ramp from $g(t) = 6t$.\n",
+    "\n",
+    "If we want to try out this formulation, we must create element and assign this equation for it:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 55,
+   "execution_count": 19,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# create volume element\n",
+    "element = Quad4([1, 2, 3, 4])\n",
+    "fieldset1 = FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
+    "fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)])\n",
+    "fieldset3 = FieldSet(\"temperature load\", [Field(0.0, [12.0, 12.0, 12.0, 12.0])])\n",
+    "fieldset4 = FieldSet(\"density\", [Field(0.0, 1.0)])\n",
+    "push!(element, fieldset1)\n",
+    "push!(element, fieldset2)\n",
+    "push!(element, fieldset3)\n",
+    "push!(element, fieldset4)\n",
+    "equation = DC2D4(element)\n",
+    "\n",
+    "# create boundary element with \n",
+    "boundary_element = Seg2([1, 2])\n",
+    "push!(boundary_element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]])]))\n",
+    "# linear ramp from 1 to 6 in time 0 to 1\n",
+    "push!(boundary_element, FieldSet(\"temperature flux\", [Field(0.0, 0.0), Field(1.0, 6.0)]))\n",
+    "boundary_equation = DC2D2(boundary_element);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
    "metadata": {
     "collapsed": false
    },
@@ -557,28 +725,22 @@
        " -1.0  -2.0  -1.0   4.0"
       ]
      },
-     "execution_count": 55,
+     "execution_count": 20,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
-    "element = Quad4([1, 2, 3, 4])\n",
-    "fieldset1 = FieldSet(\"geometry\")\n",
-    "field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
-    "push!(fieldset1, field1)\n",
-    "fieldset2 = FieldSet(\"temperature thermal conductivity\")\n",
-    "push!(fieldset2, Field(0.0, 6.0))\n",
-    "push!(element, fieldset1)\n",
-    "push!(element, fieldset2)\n",
-    "equation = DC2D4(element)\n",
-    "integrate_lhs(equation, 1.0)"
+    "using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n",
+    "\n",
+    "local_assembly = initialize_local_assembly(equation)\n",
+    "calculate_local_assembly!(local_assembly, equation)\n",
+    "local_assembly.stiffness_matrix"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 56,
+   "execution_count": 21,
    "metadata": {
     "collapsed": false
    },
@@ -586,107 +748,32 @@
     {
      "data": {
       "text/plain": [
-       "true"
+       "4x1 Array{Float64,2}:\n",
+       " 3.0\n",
+       " 3.0\n",
+       " 3.0\n",
+       " 3.0"
       ]
      },
-     "execution_count": 56,
+     "execution_count": 21,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "JuliaFEM.has_lhs(equation)"
+    "local_assembly.force_vector"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If rhs or lhs is not defined, integration returns nothing."
+    "Let's set homogeneous Dirichlet boundary condition ($u=0$) on boundary on $\\Gamma_\\mathrm{D} = \\{(x,y) | 0 \\leq x \\leq 1, y=1\\}$. With source term $f=12$ and $k=6$ we should have constant $T=1$ on free boundary. According to my math we should have $u\\left(x,y\\right)=-\\frac{1}{6}\\left(\\frac{1}{2}fx^{2}-fx\\right)$ which equals to one if $x=1$."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 57,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "true"
-      ]
-     },
-     "execution_count": 57,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "integrate_rhs(equation, 1.0) == nothing"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next heat flux on boundary:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 58,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "has_rhs (generic function with 4 methods)"
-      ]
-     },
-     "execution_count": 58,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "using JuliaFEM: get_basis, Seg2\n",
-    "\n",
-    "\"\"\"\n",
-    "Diffusive heat transfer for 2-node linear segment.\n",
-    "\"\"\"\n",
-    "type DC2D2 <: Heat\n",
-    "    element :: Seg2\n",
-    "    integration_points :: Array{IntegrationPoint, 1}\n",
-    "    global_dofs :: Array{Int64, 1}\n",
-    "end\n",
-    "\n",
-    "function DC2D2(element::Seg2)\n",
-    "    integration_points = [\n",
-    "        IntegrationPoint([0.0], 2.0)]\n",
-    "    push!(element, FieldSet(\"temperature\"))\n",
-    "    DC2D2(element, integration_points, [])\n",
-    "end\n",
-    "\n",
-    "\"\"\"\n",
-    "Right hand side defined in integration point\n",
-    "\"\"\"\n",
-    "function JuliaFEM.get_rhs(eq::DC2D2, ip, t)\n",
-    "    el = get_element(eq)\n",
-    "    h = get_basis(el, ip.xi)\n",
-    "    #f = el[\"temperature flux\"]\n",
-    "    f = interpolate(el, \"temperature flux\", ip.xi, t)\n",
-    "    return h*f\n",
-    "end\n",
-    "JuliaFEM.has_rhs(eq::DC2D2) = true"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 59,
+   "execution_count": 22,
    "metadata": {
     "collapsed": false
    },
@@ -695,21 +782,498 @@
      "data": {
       "text/plain": [
        "2-element Array{Float64,1}:\n",
-       " 50.0\n",
-       " 50.0"
+       " 1.0\n",
+       " 1.0"
       ]
      },
-     "execution_count": 59,
+     "execution_count": 22,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "element = Seg2([1, 2])\n",
-    "push!(element, FieldSet(\"geometry\", [Field(0.0, Vector[[0.0,0.0], [0.0,1.0]])]))\n",
-    "push!(element, FieldSet(\"temperature flux\", [Field(0.0, 100.0)]))\n",
-    "equation = DC2D2(element)\n",
-    "integrate_rhs(equation, 1.0)"
+    "A = local_assembly.stiffness_matrix\n",
+    "b = local_assembly.force_vector\n",
+    "u = zeros(4)\n",
+    "fdofs = [1,2]\n",
+    "u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Notice that we created a separate boundary element which integrates heat flux. Defining boundary elements goes identically with other elements. With flux $g=6$ we should have constant $T=1$ on free boundary. This should equal to $u(x,y) = x$, i.e. 1d bar stretched from other side while other one is fixed."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2-element Array{Float64,1}:\n",
+       " 1.0\n",
+       " 1.0"
+      ]
+     },
+     "execution_count": 23,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "boundary_assembly = JuliaFEM.initialize_local_assembly(boundary_equation)\n",
+    "JuliaFEM.calculate_local_assembly!(boundary_assembly, boundary_equation)\n",
+    "b = zeros(4)\n",
+    "b[fdofs] = boundary_assembly.force_vector\n",
+    "u = zeros(4)\n",
+    "u[fdofs] = A[fdofs, fdofs] \\ b[fdofs]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Some alternative ways to define equations - principle of minimum potential energy\n",
+    "\n",
+    "Fundamentally we are always looking for a solution of systems of equation $Ax = b$, despite the fact that mathematicians doesn't like this \"engineering approach\" at all. For some equations it may be too hard to write left hand side, i.e., stiffness matrix $A$. (At least for me). For this reason JuliaFEM supports automatic differentiation using package called `ForwardDiff`, which makes the analytical linearization unnecessary.\n",
+    "\n",
+    "Let's consider the following functional\n",
+    "\\begin{equation}\n",
+    "\\min\\, J\\left(u\\right)=\\int_{\\Omega}(k+6u)\\left|\\nabla u\\right|^{2}\\,\\mathrm{d}x-Pu,\n",
+    "\\end{equation}\n",
+    "where $P$ contains point loads at the free corners of the domain. This is nothing more but a redefinition of the earlier example problem defined such that the source term $k$, which in previous example was constant, is no more constant, but is replaced a term which depend from field $u$: $q(u) = k + 6u$.\n",
+    "\n",
+    "This functional doesn't actually have any or very little physical meaning, but it's non-linear due to the source term $q(u) = k+6u$ and therefore suits to our needs in demonstration purposes. In mechanical world this could be some sort of elastic spring with a spring \"constant\" depending on displacement (non-linear elasticity problem), which is compressed using point force $P$. The system is still *conservative*, i.e. it has a potential energy and can be solved using standard variational methods, finding a minimum of potential energy in suitable finite set of functions.\n",
+    "\n",
+    "In JuliaFEM we simply define it's potential energy functional $\\Pi(u)$ and let the `ForwardDiff`-package to handle the rest of the hard things. Because the problem is nonlinear, linearization and iterative solution method (=Newton's method) is needed. Accurate solution at the free end is\n",
+    "\\begin{equation}\n",
+    "u = \\frac{2}{3}.\n",
+    "\\end{equation}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "size (generic function with 65 methods)"
+      ]
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "\"\"\" Diffusive heat transfer for 4-node bilinear element, with a nonlinear source term. \"\"\"\n",
+    "type DC2D4NL <: Heat\n",
+    "    element :: Quad4\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end\n",
+    "function DC2D4NL(element::Quad4)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
+    "    push!(element, FieldSet(\"temperature\"))\n",
+    "    # Initial configuration needs to be defined\n",
+    "    push!(element[\"temperature\"], Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
+    "    DC2D4NL(element, integration_points, [])\n",
+    "end\n",
+    "Base.size(equation::DC2D4NL) = 4"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This time we don't write stiffness matrix, etc. explicitly but just write potential energy and let `ForwardDiff` take care of rest of the stuff."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 32,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "has_potential_energy (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 32,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "\"\"\" Calculate a potential Π = Wint - Wext of system. \"\"\"\n",
+    "function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    k = basis(\"temperature thermal conductivity\", ip, time)\n",
+    "    f = basis(\"temperature load\", ip, time)\n",
+    "    T = basis(\"temperature\", ip, time, variation)\n",
+    "    ∇T = grad(basis)(\"temperature\", ip, time, variation)\n",
+    "    Wint = (k + 6*T) * 1/2*∇T*∇T'\n",
+    "    Wext = f*T\n",
+    "    return Wint - Wext\n",
+    "end\n",
+    "JuliaFEM.has_potential_energy(eq::DC2D4NL) = true"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Solving non-linear variational problem is rather complicated thing due to the nature of non-linearity. Here's a heavily commented version of simple Newton iteration."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 34,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "increment  1, |du| =  1.41421, |r| =    4.243\n",
+      "increment  2, |du| =  0.42426, |r| =    4.243\n",
+      "increment  3, |du| =  0.04657, |r| =    4.243\n",
+      "increment  4, |du| =  0.00057, |r| =    4.243\n",
+      "increment  5, |du| =  0.00000, |r| =    4.243\n",
+      "elapsed time: 0.001905274 seconds\n",
+      "error: 1.6653345369377348e-15\n",
+      "temperature at free end: [0.6666666666666683,0.6666666666666683], should be 0.6666666666666666\n"
+     ]
+    }
+   ],
+   "source": [
+    "function run_simulation()\n",
+    "    # create model -- start\n",
+    "    element = Quad4([1, 2, 3, 4])\n",
+    "    fieldset1 = FieldSet(\"geometry\")\n",
+    "    field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
+    "    push!(fieldset1, field1)\n",
+    "    fieldset2 = FieldSet(\"temperature thermal conductivity\")\n",
+    "    push!(fieldset2, Field(0.0, 6.0))\n",
+    "    fieldset3 = FieldSet(\"temperature load\")\n",
+    "    push!(fieldset3, Field(0.0, 0*[12.0, 12.0, 12.0, 12.0]))\n",
+    "    fieldset4 = FieldSet(\"temperature nodal load\")\n",
+    "    push!(fieldset4, Field(0.0, [3.0, 3.0, 0.0, 0.0]))  # <-- P is defined here\n",
+    "    push!(element, fieldset1)\n",
+    "    push!(element, fieldset2)\n",
+    "    push!(element, fieldset3)\n",
+    "    push!(element, fieldset4)\n",
+    "    equation = DC2D4NL(element)\n",
+    "    # create model -- end\n",
+    "\n",
+    "    T0 = element[\"temperature\"][1]  # initial temperature field, we need something to \"variate\"\n",
+    "    la = initialize_local_assembly(equation) # create workspace for local matrices\n",
+    "    T = zeros(4) # create workspace for solution vector\n",
+    "    ΔT = zeros(4) # \n",
+    "    fd = [1, 2]\n",
+    "    tic()\n",
+    "    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n",
+    "    for i=1:5\n",
+    "        la = initialize_local_assembly(equation,la)  # empty workspace -- every iteration should start with this\n",
+    "        calculate_local_assembly!(la, equation)  # calculate local matrices\n",
+    "        ΔT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]  # <-- note sign convention, more on this below\n",
+    "        T = T + ΔT  # add increment to previous value\n",
+    "        # create a new field \"similar\" to field T0 (i.e., same dimension of field variable with new data)\n",
+    "        new_field = similar(T0, T)\n",
+    "        new_field.time = 1.0\n",
+    "        new_field.increment = i\n",
+    "        push!(element[\"temperature\"], new_field)  # add new field to \"temperature\" fieldset of element\n",
+    "        # print some convergence information\n",
+    "        @printf(\"increment %2d, |du| = %8.5f, |r| = %8.3f\\n\", i, norm(ΔT), norm(b[fd]))\n",
+    "    end\n",
+    "    toc()\n",
+    "    err = element[\"temperature\"][end].values[1] - 2/3\n",
+    "    println(\"error: $err\")\n",
+    "    blaa = element[\"temperature\"][end].values[1:2]\n",
+    "    acc = 2/3\n",
+    "    println(\"temperature at free end: $blaa, should be $acc\")\n",
+    "end\n",
+    "run_simulation()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### A short note about sign convention\n",
+    "\n",
+    "When solving a non-linear equations like $\\mathbf{r}\\left(\\mathbf{u}\\right) = 0$ it is necessary to iterate, i.e. solve several linearized solutions using e.g. Newton iterations until residual is small enough. Using truncated Taylor series the solution algorithm is in general\n",
+    "\\begin{equation}\n",
+    "\\mathbf{r}\\left(\\mathbf{u}+\\Delta\\mathbf{u}\\right)=\\mathbf{r}\\left(\\mathbf{u}\\right)+\\frac{\\partial\\mathbf{r}\\left(\\mathbf{u}\\right)}{\\partial\\mathbf{u}}\\Delta\\mathbf{u}+\\cdots\\approx\\mathbf{r}\\left(\\mathbf{u}\\right)+\\mathbf{K}\\Delta\\mathbf{u}=0\n",
+    "\\end{equation}\n",
+    "\\begin{equation}\n",
+    "\\mathbf{K}\\Delta\\mathbf{u}=-\\mathbf{r}\\left(\\mathbf{u}\\right)\n",
+    "\\end{equation}\n",
+    "with a *negative* sign in front of residual vector.\n",
+    "\n",
+    "We do however want to solve something like\n",
+    "\\begin{equation}\n",
+    "\\mathbf{A}\\Delta\\mathbf{x}^{\\left(i+1\\right)}=\\mathbf{b},\\qquad\\mathbf{x}^{\\left(i+1\\right)}=\\mathbf{x}^{\\left(i\\right)}+\\Delta\\mathbf{x}^{\\left(i+1\\right)},\n",
+    "\\end{equation}\n",
+    "where new increment is *added* to the old result. If this doesn't say anything, compare the equation to the 1d Newton iteration (https://en.wikipedia.org/wiki/Newton's_method)\n",
+    "\\begin{equation}\n",
+    "x_{n+1}=x_{n}-\\frac{f\\left(x_{n}\\right)}{f'\\left(x_{n}\\right)}\n",
+    "\\end{equation}\n",
+    "\n",
+    "Our convention is that the negative sign of residual vector is positive. It makes sense because this way the \"normal\" linear system force vector is like it's used to be in school books. Typically the residual vector is defined $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$. In our case the residual vector must be multiplied by $-1$ after the linearization so that the increment can be *added* to old result. In some text books this has already be taken into account by defining right hand side as $\\mathbf{r}=\\mathbf{f}_{\\mathrm{ext}}-\\mathbf{f}_{\\mathrm{int}}$ while linearization is still done for vector $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}$. Quite confusing!\n",
+    "\n",
+    "Just keep in mind that there is always minus somewhere and it should be. In our variational implementation the sign of residual vector is automatically changed so that next increment is  $\\mathbf{A} \\Delta\\mathbf{x}=\\mathbf{b}$. Like in the simple example above.\n",
+    "\n",
+    "### Principle of virtual work\n",
+    "\n",
+    "Defining variational form or \"energy form\", or in other words, using principle of minimum potential energy, requires system to be *conservative*, i.e. there is some energy functional describing the system. Indeed it's a very elegant approach to model e.g. hyperelasticity and other systems without a loss of energy. When system has energy dissipation, principle of minimum potential energy cannot be used for obvious reasons. There is no any \"potential function\" $\\Pi$ which could be defined and variated around it's equilibrium state. Possible situations when this happens include e.g. material plasticity or friction. In these situations the principle of virtual work is useful.\n",
+    "\n",
+    "In this situation we find equation $\\mathbf{r}=\\mathbf{f}_{\\mathrm{int}}-\\mathbf{f}_{\\mathrm{ext}}=0$ which needs to be solved. Again we let `ForwardDiff` to do the linearization of the right hand side. Let's demonstrate this too."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "has_residual_vector (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 37,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using JuliaFEM: get_field\n",
+    "\n",
+    "\"\"\" Diffusive heat transfer for 4-node bilinear element,\n",
+    "with a nonlinear term to be added later.. \"\"\"\n",
+    "type DC2D4NLY <: Heat\n",
+    "    element :: Quad4\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end\n",
+    "function DC2D4NLY(element::Quad4)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
+    "    push!(element, FieldSet(\"temperature\"))\n",
+    "    # Initial configuration needs to be defined\n",
+    "    push!(element[\"temperature\"], Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
+    "    DC2D4NLY(element, integration_points, [])\n",
+    "end\n",
+    "Base.size(equation::DC2D4NLY) = 4\n",
+    "\n",
+    "function JuliaFEM.get_residual_vector(equation::DC2D4NLY, ip, time; variation=nothing)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)\n",
+    "    k = basis(\"temperature thermal conductivity\", ip, time)\n",
+    "    f = basis(\"temperature load\", ip, time)\n",
+    "    T = get_field(basis, \"temperature\", time, variation)\n",
+    "    f_int = k * dbasis(ip,time)'*dbasis(ip,time) * T\n",
+    "    f_ext = f * basis(ip,time)'\n",
+    "    r = f_int[:] - f_ext[:]\n",
+    "    return r\n",
+    "end\n",
+    "JuliaFEM.has_residual_vector(equation::DC2D4NLY) = true"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "increment  "
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "0.020793035"
+      ]
+     },
+     "execution_count": 38,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "1, |du| =  2.82843, |r| =    4.243\n",
+      "increment  2, |du| =  0.00000, |r| =    4.243\n",
+      "increment  3, |du| =  0.00000, |r| =    4.243\n",
+      "increment  4, |du| =  0.00000, |r| =    4.243\n",
+      "increment  5, |du| =  0.00000, |r| =    4.243\n",
+      "temperature field: [2.0000000000000004,2.0000000000000004,0.0,0.0]\n",
+      "elapsed time: 0.020793035 seconds\n"
+     ]
+    }
+   ],
+   "source": [
+    "function run_simulation_2()\n",
+    "    # create model -- start\n",
+    "    element = Quad4([1, 2, 3, 4])\n",
+    "    fieldset1 = FieldSet(\"geometry\")\n",
+    "    field1 = Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
+    "    push!(fieldset1, field1)\n",
+    "    fieldset2 = FieldSet(\"temperature thermal conductivity\")\n",
+    "    push!(fieldset2, Field(0.0, 6.0))\n",
+    "    fieldset3 = FieldSet(\"temperature load\")\n",
+    "    push!(fieldset3, Field(0.0, [12.0, 12.0, 12.0, 12.0]))\n",
+    "    fieldset4 = FieldSet(\"temperature nodal load\")\n",
+    "    push!(fieldset4, Field(0.0, [3.0, 3.0, 0.0, 0.0]))  # <-- P is defined here\n",
+    "    push!(element, fieldset1)\n",
+    "    push!(element, fieldset2)\n",
+    "    push!(element, fieldset3)\n",
+    "    push!(element, fieldset4)\n",
+    "    equation = DC2D4NLY(element)\n",
+    "    # create model -- end\n",
+    "\n",
+    "    T0 = element[\"temperature\"][1]  # initial temperature field, we need something to \"variate\"\n",
+    "    la = initialize_local_assembly(equation) # create workspace for local matrices\n",
+    "    T = zeros(4) # create workspace for solution vector\n",
+    "    ΔT = zeros(4) # \n",
+    "    fd = [1, 2]\n",
+    "    tic()\n",
+    "    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n",
+    "    for i=1:5\n",
+    "        la = initialize_local_assembly(equation,la)  # empty workspace -- every iteration should start with this\n",
+    "        calculate_local_assembly!(la, equation)  # calculate local matrices\n",
+    "        ΔT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]  # <-- note sign convention, more on this below\n",
+    "        T = T + ΔT  # add increment to previous value\n",
+    "        # create a new field \"similar\" to field T0 (i.e., same dimension of field variable with new data)\n",
+    "        new_field = similar(T0, T)\n",
+    "        new_field.time = 1.0\n",
+    "        new_field.increment = i\n",
+    "        push!(element[\"temperature\"], new_field)  # add new field to \"temperature\" fieldset of element\n",
+    "        # print some convergence information\n",
+    "        @printf(\"increment %2d, |du| = %8.5f, |r| = %8.3f\\n\", i, norm(ΔT), norm(b[fd]))\n",
+    "    end\n",
+    "    println(\"temperature field: \",element[\"temperature\"](Inf).values)\n",
+    "    toc()\n",
+    "end\n",
+    "run_simulation_2()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Manually assembling local matrices\n",
+    "\n",
+    "There might still be situations where all the above methods writing field equations are just not enough. This situation can happen for example when calculating tangent stiffness matrix analytically for a nonlinear problem. In this case all matrices are usually calculated at the same time. Or maybe for educational purposes it's important to show how matrices are actually calculated. Or for debugging. \n",
+    "\n",
+    "Anyway, it's possible to override `calculate_local_assembly!` function for your own equation and after that access all the low level stuff. Here's example how to do that:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Success :: (line:-1) :: fact was true\n",
+       "  Expression: la.stiffness_matrix[free_dofs,free_dofs] \\ la.force_vector[free_dofs] --> roughly([1.0,1.0])\n",
+       "    Expected: [1.0,1.0]\n",
+       "    Occurred: [1.0000000000000002,1.0000000000000002]"
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using JuliaFEM: LocalAssembly, get_integration_points\n",
+    "\n",
+    "function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::DC2D4, time::Number=Inf)\n",
+    "    initialize_local_assembly(assembly, equation)  # zero all workspaces\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)\n",
+    "    detJ = det(basis)\n",
+    "    for ip in get_integration_points(equation)\n",
+    "        w = ip.weight * detJ(ip)\n",
+    "        # evaluate fields in integration point\n",
+    "        ρ = basis(\"density\", ip, time)\n",
+    "        k = basis(\"temperature thermal conductivity\", ip, time)\n",
+    "        f = basis(\"temperature load\", ip, time)\n",
+    "        # evaluate basis functions and gradient in integration point\n",
+    "        N = basis(ip, time)\n",
+    "        ∇N = dbasis(ip, time)\n",
+    "        # do assembly\n",
+    "        assembly.mass_matrix += w * ρ*N'*N\n",
+    "        assembly.stiffness_matrix += w * k*∇N'*∇N\n",
+    "        assembly.force_vector += w * (N'*f)[:]\n",
+    "    end\n",
+    "end\n",
+    "\n",
+    "function test_local_assembly()\n",
+    "    # create volume element\n",
+    "    element = Quad4([1, 2, 3, 4])\n",
+    "    fieldset1 = FieldSet(\"geometry\", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
+    "    fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(0.0, 6.0)])\n",
+    "    fieldset3 = FieldSet(\"temperature load\", [Field(0.0, [12.0, 12.0, 12.0, 12.0])])\n",
+    "    fieldset4 = FieldSet(\"density\", [Field(0.0, 1.0)])\n",
+    "    push!(element, fieldset1)\n",
+    "    push!(element, fieldset2)\n",
+    "    push!(element, fieldset3)\n",
+    "    push!(element, fieldset4)\n",
+    "    equation = DC2D4(element)\n",
+    "    la = initialize_local_assembly(equation)\n",
+    "    calculate_local_assembly!(la, equation)\n",
+    "    free_dofs = [1, 2]\n",
+    "    @fact la.stiffness_matrix[free_dofs, free_dofs] \\ la.force_vector[free_dofs] --> roughly([1.0, 1.0])\n",
+    "end\n",
+    "\n",
+    "test_local_assembly()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "To be continued (stuff below this is a little broken at the moment...)"
    ]
   },
   {
@@ -718,12 +1282,12 @@
    "source": [
     "## Defining own problem\n",
     "\n",
-    "- main object: takes a set of elements and maps corresponding field equations to them"
+    "Main objective of \"problem\" is to take a set of elements and map corresponding field equations to them. In this way we can solve several different fields at the same time, like for example temperature + displacement."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 60,
+   "execution_count": 34,
    "metadata": {
     "collapsed": false
    },
@@ -734,7 +1298,7 @@
        "PlaneHeatProblem"
       ]
      },
-     "execution_count": 60,
+     "execution_count": 34,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -750,7 +1314,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 61,
+   "execution_count": 35,
    "metadata": {
     "collapsed": false
    },
@@ -761,7 +1325,7 @@
        "get_equation (generic function with 6 methods)"
       ]
      },
-     "execution_count": 61,
+     "execution_count": 35,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -781,7 +1345,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 62,
+   "execution_count": 33,
    "metadata": {
     "collapsed": false
    },
@@ -790,26 +1354,23 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "20-loka 15:49:17:DEBUG:root:total dofs: 4\n"
+      "25-Oct 23:10:39:DEBUG:root:total dofs: 4\n"
      ]
     },
     {
-     "data": {
-      "text/plain": [
-       "4x1 sparse matrix with 2 Float64 entries:\n",
-       "\t[1, 1]  =  300.0\n",
-       "\t[2, 1]  =  300.0"
-      ]
-     },
-     "execution_count": 62,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[33], in expression starting on line 30",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[33], in expression starting on line 30",
+      ""
+     ]
     }
    ],
    "source": [
     "using JuliaFEM: get_connectivity, set_global_dofs!, get_global_dofs\n",
     "using JuliaFEM: add_element!, get_equations, get_matrix_dimension, calculate_global_dofs\n",
-    "using JuliaFEM: assign_global_dofs!, get_lhs, get_rhs\n",
+    "using JuliaFEM: assign_global_dofs!\n",
     "\n",
     "# create elements and add necessary properties like connectivity and geometry\n",
     "el1 = Quad4([2, 3, 4, 5])\n",
@@ -842,22 +1403,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 63,
+   "execution_count": 30,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "2x1 sparse matrix with 2 Float64 entries:\n",
-       "\t[1, 1]  =  100.0\n",
-       "\t[2, 1]  =  100.0"
-      ]
-     },
-     "execution_count": 63,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[30], in expression starting on line 2",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: A not defined\nwhile loading In[30], in expression starting on line 2",
+      ""
+     ]
     }
    ],
    "source": [
@@ -877,7 +1435,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 64,
+   "execution_count": 31,
    "metadata": {
     "collapsed": false
    },
@@ -886,10 +1444,10 @@
      "data": {
       "text/plain": [
        "1-element Array{JuliaFEM.DirichletEquation,1}:\n",
-       " JuliaFEM.DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.FieldSet(symbol(\"reaction force\"),JuliaFEM.Field[]),:geometry=>JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])]))),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)"
+       " JuliaFEM.DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.FieldSet(symbol(\"reaction force\"),JuliaFEM.Field[]),:geometry=>JuliaFEM.FieldSet(:geometry,JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,0,Array{T,1}[[0.0,0.0],[0.0,1.0]])]))),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Symbol,JuliaFEM.FieldSet}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Symbol,JuliaFEM.FieldSet}())],Int64[],fieldval)"
       ]
      },
-     "execution_count": 64,
+     "execution_count": 31,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -907,22 +1465,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 65,
+   "execution_count": 32,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "4x1 sparse matrix with 2 Float64 entries:\n",
-       "\t[3, 1]  =  0.0\n",
-       "\t[4, 1]  =  0.0"
-      ]
-     },
-     "execution_count": 65,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: integrate_lhs not defined\nwhile loading In[32], in expression starting on line 5",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: integrate_lhs not defined\nwhile loading In[32], in expression starting on line 5",
+      "",
+      " in get_lhs at /home/jukka/.julia/v0.4/JuliaFEM/src/problems.jl:119"
+     ]
     }
    ],
    "source": [
@@ -945,24 +1501,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 66,
+   "execution_count": 33,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "4x4 Array{Float64,2}:\n",
-       "  4.0  -1.0  -2.0  -1.0\n",
-       " -1.0   4.0  -1.0  -2.0\n",
-       " -2.0  -1.0   4.0  -1.0\n",
-       " -1.0  -2.0  -1.0   4.0"
-      ]
-     },
-     "execution_count": 66,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[33], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: A not defined\nwhile loading In[33], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -971,24 +1522,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 67,
+   "execution_count": 34,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "4x1 Array{Float64,2}:\n",
-       " 300.0\n",
-       " 300.0\n",
-       "   0.0\n",
-       "   0.0"
-      ]
-     },
-     "execution_count": 67,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: b not defined\nwhile loading In[34], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: b not defined\nwhile loading In[34], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -997,24 +1543,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 68,
+   "execution_count": 35,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "4x4 Array{Float64,2}:\n",
-       " 0.0  0.0  0.0       0.0     \n",
-       " 0.0  0.0  0.0       0.0     \n",
-       " 0.0  0.0  0.333333  0.166667\n",
-       " 0.0  0.0  0.166667  0.333333"
-      ]
-     },
-     "execution_count": 68,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: A2 not defined\nwhile loading In[35], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: A2 not defined\nwhile loading In[35], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -1023,24 +1564,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 69,
+   "execution_count": 36,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "4x1 Array{Float64,2}:\n",
-       " 0.0\n",
-       " 0.0\n",
-       " 0.0\n",
-       " 0.0"
-      ]
-     },
-     "execution_count": 69,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: b2 not defined\nwhile loading In[36], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: b2 not defined\nwhile loading In[36], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -1049,42 +1585,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 70,
+   "execution_count": 37,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "8x8 sparse matrix with 24 Float64 entries:\n",
-       "\t[1, 1]  =  4.0\n",
-       "\t[2, 1]  =  -1.0\n",
-       "\t[3, 1]  =  -2.0\n",
-       "\t[4, 1]  =  -1.0\n",
-       "\t[1, 2]  =  -1.0\n",
-       "\t[2, 2]  =  4.0\n",
-       "\t[3, 2]  =  -1.0\n",
-       "\t[4, 2]  =  -2.0\n",
-       "\t[1, 3]  =  -2.0\n",
-       "\t[2, 3]  =  -1.0\n",
-       "\t⋮\n",
-       "\t[8, 3]  =  0.166667\n",
-       "\t[1, 4]  =  -1.0\n",
-       "\t[2, 4]  =  -2.0\n",
-       "\t[3, 4]  =  -1.0\n",
-       "\t[4, 4]  =  4.0\n",
-       "\t[7, 4]  =  0.166667\n",
-       "\t[8, 4]  =  0.333333\n",
-       "\t[3, 7]  =  0.333333\n",
-       "\t[4, 7]  =  0.166667\n",
-       "\t[3, 8]  =  0.166667\n",
-       "\t[4, 8]  =  0.333333"
-      ]
-     },
-     "execution_count": 70,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: A not defined\nwhile loading In[37], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: A not defined\nwhile loading In[37], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -1093,24 +1606,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 71,
+   "execution_count": 38,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "8x1 sparse matrix with 4 Float64 entries:\n",
-       "\t[1, 1]  =  300.0\n",
-       "\t[2, 1]  =  300.0\n",
-       "\t[7, 1]  =  0.0\n",
-       "\t[8, 1]  =  0.0"
-      ]
-     },
-     "execution_count": 71,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: b not defined\nwhile loading In[38], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: b not defined\nwhile loading In[38], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -1126,28 +1634,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 72,
+   "execution_count": 39,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "data": {
-      "text/plain": [
-       "8x8 Array{Float64,2}:\n",
-       "  4.0  -1.0  -2.0       -1.0       0.0  0.0  0.0       0.0     \n",
-       " -1.0   4.0  -1.0       -2.0       0.0  0.0  0.0       0.0     \n",
-       " -2.0  -1.0   4.0       -1.0       0.0  0.0  0.333333  0.166667\n",
-       " -1.0  -2.0  -1.0        4.0       0.0  0.0  0.166667  0.333333\n",
-       "  0.0   0.0   0.0        0.0       0.0  0.0  0.0       0.0     \n",
-       "  0.0   0.0   0.0        0.0       0.0  0.0  0.0       0.0     \n",
-       "  0.0   0.0   0.333333   0.166667  0.0  0.0  0.0       0.0     \n",
-       "  0.0   0.0   0.166667   0.333333  0.0  0.0  0.0       0.0     "
-      ]
-     },
-     "execution_count": 72,
-     "metadata": {},
-     "output_type": "execute_result"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: Atot not defined\nwhile loading In[39], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: Atot not defined\nwhile loading In[39], in expression starting on line 1",
+      ""
+     ]
     }
    ],
    "source": [
@@ -1156,35 +1655,19 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 73,
+   "execution_count": 40,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Non-zero rows: [1,2,3,4,7,8]"
+     "ename": "LoadError",
+     "evalue": "LoadError: UndefVarError: Atot not defined\nwhile loading In[40], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: UndefVarError: Atot not defined\nwhile loading In[40], in expression starting on line 1",
+      ""
      ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "8x1 Array{Float64,2}:\n",
-       " 100.0\n",
-       " 100.0\n",
-       "   0.0\n",
-       "   0.0\n",
-       "   0.0\n",
-       "   0.0\n",
-       " 600.0\n",
-       " 600.0"
-      ]
-     },
-     "execution_count": 73,
-     "metadata": {},
-     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -1213,7 +1696,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 74,
+   "execution_count": 41,
    "metadata": {
     "collapsed": false
    },
@@ -1221,10 +1704,10 @@
     {
      "data": {
       "text/plain": [
-       "call (generic function with 1270 methods)"
+       "call (generic function with 1288 methods)"
       ]
      },
-     "execution_count": 74,
+     "execution_count": 41,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1316,7 +1799,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 75,
+   "execution_count": 42,
    "metadata": {
     "collapsed": false
    },
@@ -1325,7 +1808,16 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "20-loka 15:49:17:DEBUG:root:total dofs: 4\n"
+      "25-Oct 20:27:59:DEBUG:root:total dofs: 4\n"
+     ]
+    },
+    {
+     "ename": "LoadError",
+     "evalue": "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[42], in expression starting on line 37",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: MethodError: `has_lhs` has no method matching has_lhs(::DC2D4)\nwhile loading In[42], in expression starting on line 37",
+      ""
      ]
     }
    ],
@@ -1379,7 +1871,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 76,
+   "execution_count": 43,
    "metadata": {
     "collapsed": false
    },
@@ -1387,25 +1879,20 @@
     {
      "data": {
       "text/plain": [
-       "Success :: (line:-1) :: fact was true\n",
+       "Error :: (line:-1)\n",
        "  Expression: mean(T) --> roughly(100.0)\n",
-       "    Expected: 100.0\n",
-       "    Occurred: 100.00000000000003"
+       "  UndefVarError: T not defined\n",
+       "   in anonymous at /home/jukka/.julia/v0.4/FactCheck/src/FactCheck.jl:271\n",
+       "   in do_fact at /home/jukka/.julia/v0.4/FactCheck/src/FactCheck.jl:333\n",
+       "   in include_string at loading.jl:266\n",
+       "   in execute_request_0x535c5df2 at /home/jukka/.julia/v0.4/IJulia/src/execute_request.jl:177\n",
+       "   in eventloop at /home/jukka/.julia/v0.4/IJulia/src/IJulia.jl:141\n",
+       "   in anonymous at task.jl:447"
       ]
      },
-     "execution_count": 76,
+     "execution_count": 43,
      "metadata": {},
      "output_type": "execute_result"
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      "Residual norm: 9.845568954283847e-14\n",
-      "Array{T,1}[[0.0,0.0],[0.25,0.0],[0.5,0.0],[0.75,0.0],[1.0,0.0]]\n",
-      "[100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003,100.00000000000003]\n"
-     ]
     }
    ],
    "source": [
diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb
index 02f5773..90a1b4d 100644
--- a/notebooks/2015-06-25-elasticity-solver-example.ipynb
+++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb
@@ -10,7 +10,7 @@
     "\n",
     "**Abstract**: Elasticity equations design notes.\n",
     "\n",
-    "##Weak form\n",
+    "###Weak form\n",
     "\n",
     "Given function spaces\n",
     "\\begin{align}\n",
@@ -20,6 +20,38 @@
     "find $\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{U}}$ such that\n",
     "\\begin{equation}\n",
     "\\delta\\mathcal{W}:=\\int_{\\Omega_{0}}\\rho_{0}\\ddot{\\boldsymbol{u}}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}+\\int_{\\Omega_{0}}\\boldsymbol{S}:\\delta\\boldsymbol{E}\\,\\mathrm{d}V_{0}-\\int_{\\Omega_{0}}\\hat{\\boldsymbol{b}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}V_{0}-\\int_{\\Gamma_{\\sigma}}\\hat{\\boldsymbol{t}}_{0}\\cdot\\delta\\boldsymbol{u}\\,\\mathrm{d}A_{0} =0 \\qquad\\forall\\delta\\boldsymbol{u}\\in\\boldsymbol{\\mathcal{V}}\n",
+    "\\end{equation}\n",
+    "\n",
+    "### Some formulas\n",
+    "\\begin{align}\n",
+    "J & =\\det\\left(F\\right)\\\\\n",
+    "I_{c} & =\\mbox{tr}\\left(C\\right)\\\\\n",
+    "\\mathbf{C} & =\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}\\\\\n",
+    "\\mathbf{F} & =\\mathbf{I}+\\nabla\\mathbf{u}\\\\\n",
+    "\\mathbf{E} & =\\frac{1}{2}\\left(\\mathbf{F}^{\\mathrm{T}}\\mathbf{F}-\\mathbf{I}\\right)\n",
+    "\\end{align}\n",
+    "\n",
+    "### Potential energy\n",
+    "\n",
+    "\\begin{equation}\n",
+    "\\underset{u\\in\\boldsymbol{\\mathcal{U}}}{\\min}\\Pi\\left(\\mathbf{u}\\right)\n",
+    "\\end{equation}\n",
+    "\\begin{equation}\n",
+    "\\Pi\\left(\\mathbf{u}\\right)=\\int_{\\Omega}\\psi\\left(\\mathbf{u}\\right)-\\int_{\\Omega}\\hat{\\mathbf{b}}_{0}\\cdot\\mathbf{u}-\\int_{\\Gamma_{\\sigma}}\\hat{\\mathbf{t}}_{0}\\cdot\\mathbf{u}\\,\\mathrm{d}A_{0}\n",
+    "\\end{equation}\n",
+    "\n",
+    "### Material models\n",
+    "\n",
+    "https://en.wikipedia.org/wiki/Strain_energy_density_function\n",
+    "\n",
+    "Saint-Venant-Kirchhoff model https://en.wikipedia.org/wiki/Hyperelastic_material\n",
+    "\\begin{equation}\n",
+    "\\psi\\left(\\mathbf{E}\\right)=\\frac{\\lambda}{2}\\left[\\mbox{tr}\\left(\\mathbf{E}\\right)\\right]^{2}+\\mu\\mbox{tr}\\left(\\mathbf{E}^2\\right)\n",
+    "\\end{equation}\n",
+    "\n",
+    "neo-Hookean material https://en.wikipedia.org/wiki/Neo-Hookean_solid\n",
+    "\\begin{equation}\n",
+    "\\psi=\\frac{\\mu}{2}\\left(I_{c}-3\\right)-\\mu\\ln\\left(J\\right)+\\frac{\\lambda}{2}\\ln\\left(J\\right)^{2}\n",
     "\\end{equation}"
    ]
   },
@@ -43,8 +75,9 @@
    ],
    "source": [
     "using ForwardDiff\n",
-    "using JuliaFEM: Quad4, Field, FieldSet, IntegrationPoint\n",
-    "using JuliaFEM: interpolate, get_element, get_dbasisdX, dinterpolate\n",
+    "using JuliaFEM: Quad4, Field, FieldSet, IntegrationPoint, Equation, LocalAssembly\n",
+    "using JuliaFEM: get_element, get_basis, grad, get_integration_points\n",
+    "using JuliaFEM: initialize_local_assembly, calculate_local_assembly!\n",
     "using Logging\n",
     "using FactCheck\n",
     "Logging.configure(level=DEBUG)"
@@ -60,7 +93,7 @@
     {
      "data": {
       "text/plain": [
-       "CPS4"
+       "size (generic function with 63 methods)"
       ]
      },
      "execution_count": 3,
@@ -69,10 +102,12 @@
     }
    ],
    "source": [
-    "abstract Elasticity <: JuliaFEM.Equation\n",
+    "abstract Elasticity <: Equation\n",
     "abstract PlaneElasticity <: Elasticity\n",
     "abstract PlaneStressElasticity <: PlaneElasticity\n",
     "\n",
+    "JuliaFEM.get_unknown_field_name(equation::Elasticity) = \"displacement\"\n",
+    "\n",
     "\"\"\" Plane stress formulation for 4-node bilinear element. \"\"\"\n",
     "type CPS4 <: PlaneStressElasticity\n",
     "    element :: Quad4\n",
@@ -87,7 +122,9 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
     "    push!(element, FieldSet(\"displacement\"))\n",
     "    CPS4(element, integration_points, [])\n",
-    "end"
+    "end\n",
+    "\n",
+    "JuliaFEM.size(eq::CPS4) = 8\n"
    ]
   },
   {
@@ -101,7 +138,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 25,
    "metadata": {
     "collapsed": false
    },
@@ -109,74 +146,83 @@
     {
      "data": {
       "text/plain": [
-       "has_rhs (generic function with 4 methods)"
+       "has_potential_energy (generic function with 2 methods)"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 25,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "function get_lhs_and_rhs(equation::CPS4, ip, time)\n",
-    "    # boilerplate code start\n",
+    "\"\"\"\n",
+    "Calculate internal energy of system. This can be\n",
+    "used to define own material models.\n",
+    "\n",
+    "Parameters\n",
+    "----------\n",
+    "equation\n",
+    "    field equation we are solving\n",
+    "ip\n",
+    "    integration point, which can be used to access fields\n",
+    "time\n",
+    "    current time\n",
+    "F\n",
+    "    deformation gradient\n",
+    "\n",
+    "Returns\n",
+    "-------\n",
+    "Internal energy of system.\n",
+    "\"\"\"\n",
+    "function calculate_internal_energy(equation::Equation, ip::IntegrationPoint, time::Number, F::Matrix)\n",
     "    element = get_element(equation)\n",
-    "    geometry = element[\"geometry\"](time)\n",
-    "    basis = FEM.get_basis(element)\n",
-    "    dbasis = FEM.diff(basis)(ip.xi)\n",
-    "    grad(u) = dbasis*u*inv(dbasis*geometry)\n",
-    "    # boilerplate code end -- replace with a macro?\n",
+    "    basis = get_basis(element)\n",
     "\n",
-    "    # interpolate fields in temporal dimension\n",
-    "    young = element[\"youngs modulus\"](time)\n",
-    "    poisson = element[\"poissons ratio\"](time)\n",
-    "    displacement = element[\"displacement\"](time)\n",
-    "\n",
-    "    # interpolate material in spatial dimension\n",
-    "    young = interpolate(basis, young, ip)\n",
-    "    poisson = interpolate(basis, poisson, ip)\n",
+    "    # material parameters\n",
+    "    young = basis(\"youngs modulus\", ip, time)\n",
+    "    poisson = basis(\"poissons ratio\", ip, time)\n",
     "    mu = young/(2*(1+poisson))\n",
     "    lambda = young*poisson/((1+poisson)*(1-2*poisson))\n",
-    "    lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d\n",
-    "\n",
-    "    function W(data::Vector)\n",
-    "        Wint = 0.0\n",
-    "        # create new field u, similar to field displacement, and fill it with data\n",
-    "        u = similar(displacement, data)\n",
-    "        F = I + grad(u)  # deformation gradient\n",
-    "        E = 1/2*(F'*F - I)  # strain\n",
-    "        S = 2*mu*E + lambda*trace(E)*I  # stress\n",
-    "        Wint += 1/2*trace(S*E')\n",
-    "\n",
-    "        Wext = 0.0\n",
-    "        # any volume load?\n",
-    "        if haskey(element, \"displacement volume load\")\n",
-    "            b = interpolate(element, \"displacement volume load\", ip, time)\n",
-    "            δu = interpolate(basis, u, ip)\n",
-    "            Wext += dot(b, δu)\n",
-    "        end\n",
-    "        return Wint - Wext\n",
+    "    if isa(equation, PlaneStressElasticity)\n",
+    "        lambda = 2*lambda*mu/(lambda + 2*mu)  # <- correction for 2d\n",
     "    end\n",
     "\n",
-    "    R = ForwardDiff.gradient(W, displacement[:])\n",
-    "    K = ForwardDiff.hessian(W, displacement[:])\n",
-    "    return K, -R\n",
+    "    E = 1/2*(F'*F - I)  # strain\n",
+    "    Wint = 1/2*lambda*trace(E)^2 + mu*trace(E*E')\n",
+    "    # alternative way to calculate this:\n",
+    "    # S = lambda*trace(E)*I + 2*mu*E\n",
+    "    # Wint = 1/2*trace(S*E')\n",
+    "    return Wint\n",
     "end\n",
     "\n",
-    "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n",
-    "    get_lhs_and_rhs(equation, ip, time)[1]\n",
-    "end\n",
-    "function JuliaFEM.get_rhs(equation::CPS4, ip, time)\n",
-    "    get_lhs_and_rhs(equation, ip, time)[2]\n",
+    "function JuliaFEM.get_potential_energy(equation::Elasticity, ip, time; variation=nothing)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)\n",
+    "\n",
+    "    u = basis(\"displacement\", ip, time, variation)\n",
+    "    ∇u = dbasis(\"displacement\", ip, time, variation)\n",
+    "    F = I + ∇u  # deformation gradient\n",
+    "\n",
+    "    # internal energy\n",
+    "    Wint = calculate_internal_energy(equation, ip, time, F)\n",
+    "\n",
+    "    # external energy -- any volume load?\n",
+    "    Wext = 0.0\n",
+    "    if haskey(element, \"displacement volume load\")\n",
+    "        b = basis(\"displacement volume load\", ip, time)\n",
+    "        Wext += dot(b, u)\n",
+    "    end\n",
+    "\n",
+    "    return Wint - Wext\n",
     "end\n",
     "\n",
-    "JuliaFEM.has_lhs(equation::CPS4) = true\n",
-    "JuliaFEM.has_rhs(equation::CPS4) = true"
+    "JuliaFEM.has_potential_energy(equation::CPS4) = true"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 26,
    "metadata": {
     "collapsed": false
    },
@@ -185,19 +231,8 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "testing primary field with point load versus code aster solution\n",
-      "increment  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",
-      "increment  7, norm =  0.00000, Wint =   83.109, Wext =   83.109, |Wint-Wext| =  0.00000\n",
-      "increment  8, norm =  0.00000, Wint =   83.109, Wext =   83.109, |Wint-Wext| =  0.00000\n",
-      "increment  9, norm =  0.00000, Wint =   83.109, Wext =   83.109, |Wint-Wext| =  0.00000\n",
-      "increment 10, norm =  0.00000, Wint =   83.109, Wext =   83.109, |Wint-Wext| =  0.00000\n",
-      "elapsed time: 4.257759666 seconds\n",
-      "1 fact verified.\n"
+      "testing primary field with nodal load versus code aster solution\n",
+      "increment "
      ]
     },
     {
@@ -206,50 +241,65 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 26,
      "metadata": {},
      "output_type": "execute_result"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      " 1, |du| =  5.77653\n",
+      "increment  2, |du| =  0.99988\n",
+      "increment  3, |du| =  0.28354\n",
+      "increment  4, |du| =  0.07071\n",
+      "increment  5, |du| =  0.00082\n",
+      "increment  6, |du| =  0.00000\n",
+      "increment  7, |du| =  0.00000\n",
+      "increment  8, |du| =  0.00000\n",
+      "increment  9, |du| =  0.00000\n",
+      "increment 10, |du| =  0.00000\n",
+      "elapsed time: 0.03522289 seconds\n",
+      "1 fact verified.\n"
+     ]
     }
    ],
    "source": [
-    "facts(\"testing primary field with point load versus code aster solution\") do\n",
+    "facts(\"testing primary field with nodal load versus code aster solution\") do\n",
     "    element = Quad4([1, 2, 3, 4])\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(\"youngs modulus\", [Field(0.0, 500.0)]))\n",
     "    push!(element, FieldSet(\"poissons ratio\", [Field(0.0, 0.3)]))\n",
+    "    push!(element, FieldSet(\"displacement nodal load\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, -20.0], [0.0, 0.0]])]))\n",
     "    equation = CPS4(element)\n",
-    "\n",
+    "    \n",
     "    u0 = Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n",
     "    push!(element[\"displacement\"], u0)\n",
     "    u = zeros(8)\n",
     "    du = zeros(8)\n",
     "    fd = [3, 4, 5, 6]\n",
-    "    f = zeros(8)\n",
-    "    f[6] = -20.0\n",
+    "    la = initialize_local_assembly(equation)\n",
     "    tic()\n",
     "    for i=1:10\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",
+    "        la = initialize_local_assembly(equation, la)\n",
+    "        calculate_local_assembly!(la, equation)\n",
+    "        du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]\n",
     "        u += du\n",
     "        new_field = similar(u0, u)\n",
     "        new_field.time = 1.0\n",
     "        new_field.increment = i\n",
     "        push!(element[\"displacement\"], new_field)\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",
+    "        @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n",
     "    end\n",
     "    toc()\n",
     "    # verified using Code Aster.\n",
-    "    @fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-4.15546385452579E+00)\n",
+    "    @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-4.15546385452579E+00)\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 27,
    "metadata": {
     "collapsed": false
    },
@@ -259,7 +309,7 @@
      "output_type": "stream",
      "text": [
       "testing primary field with volume load versus code aster solution\n",
-      "increment  "
+      "increment "
      ]
     },
     {
@@ -268,7 +318,7 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 27,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -276,17 +326,17 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "1, norm = 14.44128, Wint = -509.167, Wext =    0.000, |Wint-Wext| = 509.16667\n",
-      "increment  2, norm =  4.01742, Wint = 4335.728, Wext =    0.000, |Wint-Wext| = 4335.72801\n",
-      "increment  3, norm =  1.54645, Wint =  705.913, Wext =    0.000, |Wint-Wext| = 705.91292\n",
-      "increment  4, norm =  1.12361, Wint =   60.492, Wext =    0.000, |Wint-Wext| = 60.49208\n",
-      "increment  5, norm =  0.79119, Wint =   -1.486, Wext =    0.000, |Wint-Wext| =  1.48555\n",
-      "increment  6, norm =  0.12733, Wint =    5.743, Wext =    0.000, |Wint-Wext| =  5.74331\n",
-      "increment  7, norm =  0.00725, Wint =    0.080, Wext =    0.000, |Wint-Wext| =  0.07997\n",
-      "increment  8, norm =  0.00001, Wint =    0.000, Wext =    0.000, |Wint-Wext| =  0.00045\n",
-      "increment  9, norm =  0.00000, Wint =    0.000, Wext =    0.000, |Wint-Wext| =  0.00000\n",
-      "increment 10, norm =  0.00000, Wint =    0.000, Wext =    0.000, |Wint-Wext| =  0.00000\n",
-      "elapsed time: 0.068575305 seconds\n",
+      " 1, |du| = 14.44128\n",
+      "increment  2, |du| =  4.01742\n",
+      "increment  3, |du| =  1.54645\n",
+      "increment  4, |du| =  1.12361\n",
+      "increment  5, |du| =  0.79119\n",
+      "increment  6, |du| =  0.12733\n",
+      "increment  7, |du| =  0.00725\n",
+      "increment  8, |du| =  0.00001\n",
+      "increment  9, |du| =  0.00000\n",
+      "increment 10, |du| =  0.00000\n",
+      "elapsed time: 0.004752839 seconds\n",
       "1 fact verified.\n"
      ]
     }
@@ -306,26 +356,22 @@
     "    u = zeros(8)\n",
     "    du = zeros(8)\n",
     "    fd = [3, 4, 5, 6]\n",
-    "    f = zeros(8)\n",
-    "    f[6] = -20.0*0\n",
+    "    la = initialize_local_assembly(equation)\n",
     "    tic()\n",
     "    for i=1:10\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",
+    "        la = initialize_local_assembly(equation, la)\n",
+    "        calculate_local_assembly!(la, equation)\n",
+    "        du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]\n",
     "        u += du\n",
     "        new_field = similar(u0, u)\n",
     "        new_field.time = 1.0\n",
     "        new_field.increment = i\n",
     "        push!(element[\"displacement\"], new_field)\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",
+    "        @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n",
     "    end\n",
     "    toc()\n",
     "    # verified using Code Aster.\n",
-    "    @fact interpolate(element, \"displacement\", [1.0, 1.0], Inf)[2] --> roughly(-8.77303119819776E+00)\n",
+    "    @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-8.77303119819776E+00)\n",
     "end"
    ]
   },
@@ -333,104 +379,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Method 2, Voigt notation, analytical linearization"
+    "### Method 2, Voigt notation, analytical linearization\n",
+    "\n",
+    "I think I don't need to mention which way is more elegant. Here's the linearization of system is done manually anyway, maybe it has better performance."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 170,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_rhs (generic function with 6 methods)"
-      ]
-     },
-     "execution_count": 170,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "function get_lhs_and_rhs(equation::CPS4, ip, time)\n",
-    "    element = get_element(equation)\n",
-    "\n",
-    "    # fields\n",
-    "    X = element[\"geometry\"](time)\n",
-    "    u = element[\"displacement\"](time)\n",
-    "    young = element[\"young\"](time)\n",
-    "    poisson = element[\"poisson\"](time)\n",
-    "\n",
-    "    # material\n",
-    "    N = FEM.get_basis(element)\n",
-    "    young = interpolate(N, young, ip)\n",
-    "    poisson = interpolate(N, poisson, ip)\n",
-    "\n",
-    "    dN = FEM.diff(N)(ip.xi)\n",
-    "    invJ = inv(dN*X)\n",
-    "    dNdX = dN*invJ\n",
-    "\n",
-    "    # kinematics\n",
-    "    gradu = dN*u*invJ\n",
-    "    F = I + gradu  # deformation gradient\n",
-    "    E = 1/2*(F'*F - I)  # GL strain tensor\n",
-    "    E = [E[1,1], E[2,2], E[1,2]*2]  # go to Voigt\n",
-    "\n",
-    "    # constitutive equations\n",
-    "    D = young/(1-poisson^2) * [1 poisson 0; poisson 1 0; 0 0 1/2*(1-poisson)]\n",
-    "    S = D*E\n",
-    "    T = zeros(4, 4)\n",
-    "    T[1,1] = S[1]\n",
-    "    T[2,2] = S[2]\n",
-    "    T[1,2] = T[2,1] = S[3]\n",
-    "    T[3:4,3:4] = T[1:2,1:2]\n",
-    "\n",
-    "    # linear part\n",
-    "    B_L = zeros(3, 8)\n",
-    "    for i=1:4\n",
-    "        B_L[1, 2*(i-1)+1] = F[1,1]*dNdX[i,1]\n",
-    "        B_L[1, 2*(i-1)+2] = F[2,1]*dNdX[i,1]\n",
-    "        B_L[2, 2*(i-1)+1] = F[1,2]*dNdX[i,2]\n",
-    "        B_L[2, 2*(i-1)+2] = F[2,2]*dNdX[i,2]\n",
-    "        B_L[3, 2*(i-1)+1] = F[1,1]*dNdX[i,2] + F[1,2]*dNdX[i,1]\n",
-    "        B_L[3, 2*(i-1)+2] = F[2,1]*dNdX[i,2] + F[2,2]*dNdX[i,1]\n",
-    "    end\n",
-    "    K_L = B_L'*D*B_L\n",
-    "\n",
-    "    # nonlinear part\n",
-    "    B_NL = zeros(4, 8)\n",
-    "    for i=1:4\n",
-    "        B_NL[1, 2*(i-1)+1] = dNdX[i,1]\n",
-    "        B_NL[2, 2*(i-1)+1] = dNdX[i,2]\n",
-    "        B_NL[3, 2*(i-1)+2] = dNdX[i,1]\n",
-    "        B_NL[4, 2*(i-1)+2] = dNdX[i,2]\n",
-    "    end\n",
-    "    K_NL = B_NL'*T*B_NL\n",
-    "\n",
-    "    fint = B_L'*S\n",
-    "\n",
-    "    R = fint\n",
-    "    Kt = K_L + K_NL\n",
-    "\n",
-    "    print('.')\n",
-    "    \n",
-    "    return Kt, -R\n",
-    "end\n",
-    "\n",
-    "function JuliaFEM.get_lhs(equation::CPS4, ip, time)\n",
-    "    get_lhs_and_rhs(equation, ip, time)[1]\n",
-    "end\n",
-    "function JuliaFEM.get_rhs(equation::CPS4, ip, time)\n",
-    "    get_lhs_and_rhs(equation, ip, time)[2]\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 172,
+   "execution_count": 24,
    "metadata": {
     "collapsed": false
    },
@@ -439,24 +395,17 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      ".."
+      "testing primary field with nodal load versus code aster solution\n",
+      "increment "
      ]
     },
     {
      "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  "
+       "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 172,
+     "execution_count": 24,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -464,70 +413,143 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      ".."
+      " 1, |du| =  5.77653\n",
+      "increment  2, |du| =  1.05083\n",
+      "increment  3, |du| =  0.39176\n",
+      "increment  4, |du| =  0.21527\n",
+      "increment  5, |du| =  0.16285\n",
+      "increment  6, |du| =  0.13115\n",
+      "increment  7, |du| =  0.10774\n",
+      "increment  8, |du| =  0.08994\n",
+      "increment  9, |du| =  0.07617\n",
+      "increment 10, |du| =  0.06533\n",
+      "elapsed time: 0.004335419 seconds\n",
+      "  Failure :: (line:-1) :: fact was false\n",
+      "    Expression: ((get_basis(element))(\"displacement\",[1.0,1.0]))[2] --> roughly(-4.15546385452579)\n",
+      "      Expected: -5.091745430627231 ≅ -4.15546385452579\n",
+      "Out of 1 total fact:\n",
+      "  Failed:   1\n"
      ]
     }
    ],
    "source": [
-    "JuliaFEM.integrate_lhs(equation, Inf)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 169,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "2.927103942720631"
-      ]
-     },
-     "execution_count": 169,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "(1/2*u'*JuliaFEM.integrate_lhs(equation, 0.0)*u)[1]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 30,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_lhs (generic function with 2 methods)"
-      ]
-     },
-     "execution_count": 30,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"\n",
-    "1-node point force element for plane stress problems.\n",
-    "\"\"\"\n",
-    "type CPS1 <: Elasticity\n",
-    "    element :: Point1\n",
+    "\"\"\" Plane stress formulation for 4-node bilinear element, manual formulation. \"\"\"\n",
+    "type CPS4M <: PlaneStressElasticity\n",
+    "    element :: Quad4\n",
     "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
     "end\n",
-    "function CPS1(el::Point1)\n",
-    "    integration_points = []\n",
-    "    set_field(el, \"displacement\", zeros(2, 1))\n",
-    "    set_field(el, \"displacement nodal load\", zeros(2, 1))\n",
-    "    CPS1(el, integration_points)\n",
+    "function CPS4M(element::Quad4)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
+    "    push!(element, FieldSet(\"displacement\"))\n",
+    "    CPS4M(element, integration_points, [])\n",
     "end\n",
-    "get_rhs(eq::CPS1) = get_field(get_element(eq), \"displacement nodal load\")\n",
-    "get_lhs(eq::CPS1) = None"
+    "\n",
+    "JuliaFEM.size(eq::CPS4M) = 8\n",
+    "\n",
+    "function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::CPS4M, time::Number=Inf)\n",
+    "    initialize_local_assembly(assembly, equation)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)\n",
+    "    detJ = det(basis)\n",
+    "    ndofs = size(equation)\n",
+    "    nnodes = round(Int, ndofs/2)\n",
+    "\n",
+    "    B_L = zeros(3, ndofs)\n",
+    "    B_NL = zeros(4, ndofs)\n",
+    "\n",
+    "    for ip in get_integration_points(equation)\n",
+    "\n",
+    "        fill!(B_L, 0.0)\n",
+    "        fill!(B_NL, 0.0)\n",
+    "\n",
+    "        u = basis(\"displacement\", ip, time)\n",
+    "        ∇u = dbasis(\"displacement\", ip, time)\n",
+    "        young = basis(\"youngs modulus\", ip, time)\n",
+    "        poisson = basis(\"poissons ratio\", ip, time)\n",
+    "\n",
+    "        # kinematics\n",
+    "        F = I + ∇u  # deformation gradient\n",
+    "        E = 1/2*(F'*F - I)  # GL strain tensor\n",
+    "        E = [E[1,1], E[2,2], E[1,2]*2]  # go to Voigt\n",
+    "\n",
+    "        # constitutive equations -- calculate stress\n",
+    "        D = young/(1-poisson^2) * [1 poisson 0; poisson 1 0; 0 0 1/2*(1-poisson)]\n",
+    "        S = D*E\n",
+    "        T = zeros(4, 4)\n",
+    "        T[1,1] = S[1]\n",
+    "        T[2,2] = S[2]\n",
+    "        T[1,2] = T[2,1] = S[3]\n",
+    "        T[3:4,3:4] = T[1:2,1:2]\n",
+    "\n",
+    "        dNdX = dbasis(ip, time)\n",
+    "\n",
+    "        # linear part\n",
+    "        for i=1:nnodes\n",
+    "            B_L[1, 2*(i-1)+1] = F[1,1]*dNdX[1,i]\n",
+    "            B_L[1, 2*(i-1)+2] = F[2,1]*dNdX[1,i]\n",
+    "            B_L[2, 2*(i-1)+1] = F[1,2]*dNdX[2,i]\n",
+    "            B_L[2, 2*(i-1)+2] = F[2,2]*dNdX[2,i]\n",
+    "            B_L[3, 2*(i-1)+1] = F[1,1]*dNdX[2,i] + F[1,2]*dNdX[1,i]\n",
+    "            B_L[3, 2*(i-1)+2] = F[2,1]*dNdX[2,i] + F[2,2]*dNdX[1,i]\n",
+    "        end\n",
+    "\n",
+    "        # nonlinear part\n",
+    "        for i=1:nnodes\n",
+    "            B_NL[1, 2*(i-1)+1] = dNdX[1,i]\n",
+    "            B_NL[2, 2*(i-1)+1] = dNdX[2,i]\n",
+    "            B_NL[3, 2*(i-1)+2] = dNdX[1,i]\n",
+    "            B_NL[4, 2*(i-1)+2] = dNdX[2,i]\n",
+    "        end\n",
+    "\n",
+    "        s = ip.weight*detJ(ip)\n",
+    "        assembly.stiffness_matrix += s*B_L'*D*B_L + s*B_NL'*T*B_NL\n",
+    "        assembly.force_vector += -s*B_L'*S\n",
+    "\n",
+    "    end\n",
+    "\n",
+    "    if haskey(element, \"displacement nodal load\")\n",
+    "        assembly.force_vector += element[\"displacement nodal load\"](time)[:]\n",
+    "    end\n",
+    "\n",
+    "end\n",
+    "\n",
+    "facts(\"testing primary field with nodal load versus code aster solution\") do\n",
+    "    element = Quad4([1, 2, 3, 4])\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(\"youngs modulus\", [Field(0.0, 500.0)]))\n",
+    "    push!(element, FieldSet(\"poissons ratio\", [Field(0.0, 0.3)]))\n",
+    "    push!(element, FieldSet(\"displacement nodal load\", [Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, -20.0], [0.0, 0.0]])]))\n",
+    "    equation = CPS4M(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",
+    "    push!(element[\"displacement\"], u0)\n",
+    "    u = zeros(8)\n",
+    "    du = zeros(8)\n",
+    "    \n",
+    "    fd = [3, 4, 5, 6]\n",
+    "    la = initialize_local_assembly(equation)\n",
+    "    tic()\n",
+    "    for i=1:10\n",
+    "        la = initialize_local_assembly(equation, la)\n",
+    "        calculate_local_assembly!(la, equation)\n",
+    "        du[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]\n",
+    "        u += du\n",
+    "        new_field = similar(u0, u)\n",
+    "        new_field.time = 1.0\n",
+    "        new_field.increment = i\n",
+    "        push!(element[\"displacement\"], new_field)\n",
+    "        @printf(\"increment %2d, |du| = %8.5f\\n\", i, norm(du))\n",
+    "    end\n",
+    "    toc()\n",
+    "    # verified using Code Aster. Noh, toimi se eilen.\n",
+    "    @fact get_basis(element)(\"displacement\", [1.0, 1.0])[2] --> roughly(-4.15546385452579E+00)\n",
+    "end"
    ]
   },
   {
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index b5db2a5..2314de1 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -10,6 +10,9 @@ using Lexicon
 using Logging
 @Logging.configure(level=DEBUG)
 
+using ForwardDiff
+autodiffcache = ForwardDiffCache()
+
 """ Simple linspace extension to arrays.
 
 Examples
diff --git a/src/elements.jl b/src/elements.jl
index f882ba3..d4fe7b9 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -8,20 +8,18 @@ Related notebooks
 2015-08-29-developing-juliafem.ipynb
 =#
 
-using JuliaFEM: interpolate
 using FactCheck
 using ForwardDiff
 
-
 abstract Element
 
 """ Get FieldSet from element. """
-function Base.getindex(element::Element, field_name::Union{Symbol, ASCIIString})
+function Base.getindex(element::Element, field_name)
     element.fields[symbol(field_name)]
 end
 
 """ Add new FieldSet to element. """
-function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name::Union{Symbol, ASCIIString})
+function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name)
     fieldset.name = symbol(fieldset_name)
     element.fields[fieldset.name] = fieldset
 end
@@ -126,103 +124,150 @@ function test_element(element_type)
     fieldset = FieldSet("field1")
     push!(fieldset, field)
     push!(element, fieldset)
-    @fact element["field1"][1] --> fld
+    push!(element, FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])]))
 
     # evaluate basis functions at middle point of element
     mid = zeros(dim)
     try
-        get_basis(element)(mid)
+        basis = get_basis(element)
+        val1 = basis(mid, 0.0)
+        Logging.info("basis at $mid: $val1")
+        val2 = basis("field1", mid, 0.0)
+        Logging.info("field val at $mid: $val2")
     catch
         Logging.error("""
         Unable to evaluate basis, define function 'get_basis' for
         this element.""")
     end
     try
-        get_dbasisdxi(element)(mid)
+        basis = get_basis(element)
+        dbasis = grad(basis)
+        val3 = dbasis(mid, 0.0)
+        Logging.info("derivative of basis at $mid: $val3")
+        val4 = dbasis("field1", mid, 0.0)
+        Logging.info("field val at $mid: $val4")
     catch
         Logging.error("""
         Unable to evaluate partial derivatives of basis,
         define function 'get_dbasisdxi' for this element.""")
     end
 
-    Logging.info("Interpolating scalar field at $mid")
-    i = interpolate(element, "field1", mid, 0.0)
-    Logging.info("Value: $i")
     Logging.info("Element $element_type passed tests.")
 end
 
 
-get_connectivity(el::Element) = el.connectivity
-
-""" Get basis functions of element. """
-get_basis(el::Element) = el.basis
-get_basis(el::Element, xi::Vector) = el.basis(xi)
-Base.call(el::Element, xi::Vector) = el.basis(xi)
-
-""" Get partial derivatives of basis functions of element. """
-get_dbasisdxi(el::Element) = el.basis.dbasisdxi
-get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi)
-get_dbasisdxi(el::Element, ip::IntegrationPoint) = el.basis.dbasisdxi(ip.xi)
-
-""" Interpolate field on element. """
-function interpolate(element::Element, field_name, xi::Vector, time::Number)
-    fieldset = element[field_name]
-    field = interpolate(fieldset, time)
-    basis = get_basis(element)
-    interpolate(basis, field, xi)
-end
-function interpolate(element::Element, field_name, ip::IntegrationPoint, time::Number)
-    interpolate(element, field_name, ip.xi, time)
+function get_connectivity(el::Element)
+    el.connectivity
 end
 
-""" Interpolate derivative of field on element. """
-function dinterpolate(element::Element, field_name, xi::Vector, time::Number)
-    fieldset = element[field_name]
-    field = interpolate(fieldset, time)
-    basis = get_basis(element)
-    dinterpolate(basis, field, xi)
+type MixedFunctionSpace
+    element1 :: Element
+    element2 :: Element
 end
-function dinterpolate(element::Element, field_name, ip::IntegrationPoint, time::Number)
-    dinterpolate(element, field_name, ip.xi, time)
+
+type FunctionSpace
+    element :: Element
 end
 
+type GradientFunctionSpace
+    element :: Element
+end
+
+function grad(u::FunctionSpace)
+    GradientFunctionSpace(u.element)
+end
+
+""" Evaluate field on element function space. """
+function call(u::FunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
+    f = !isa(variation, Void) ? variation : u.element[field_name](t)
+    if length(f) == 1
+        return f.values
+    end
+    h = u.element.basis.basis(xi)
+    return h*f
+end
+
+""" If basis is called without a field, return basis functions evaluated at that point. """
+function call(u::FunctionSpace, xi::Vector, t::Number=Inf)
+    return u.element.basis.basis(xi)'
+end
+
+""" Evaluate gradient of field on element function space. """
+function call(gradu::GradientFunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
+    f = !isa(variation, Void) ? variation : gradu.element[field_name](t)
+    X = gradu.element["geometry"](t)
+    b = gradu.element.basis.dbasisdxi(xi)
+    return b*f*inv(b*X)
+end
+
+""" If gradient of basis is called without a field, return "empty" gradient evaluated at that point. """
+function call(gradu::GradientFunctionSpace, xi::Vector, t::Number=Inf)
+    X = gradu.element["geometry"](t)
+    b = gradu.element.basis.dbasisdxi(xi)
+    return (b*inv(b*X))'
+end
+
+# on-line functions to get api more easy to use, ip -> xi.ip
+call(u::FunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t)
+call(u::FunctionSpace, ip::IntegrationPoint) = call(u, ip.xi)
+call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number) = call(u, ip.xi, t)
+call(u::GradientFunctionSpace, ip::IntegrationPoint) = call(u, ip.xi)
+
+""" Return field from function space. """
+function get_field(u::FunctionSpace, field_name, time=Inf)
+    return u.element[field_name](time)
+end
+
+""" Return field from function space. """
+function get_field(u::FunctionSpace, field_name, time=Inf, variation=nothing)
+    return !isa(variation, Void) ? variation : u.element[field_name](time)
+end
+
+""" Return fieldset from function space. """
+function get_fieldset(u::FunctionSpace, field_name)
+    return u.element[field_name]
+end
+
+# i think these will be the most called functions.
+call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+
+function jacobian(u::FunctionSpace, xi, t)
+    u.element.basis.dbasisdxi(xi)*u.element["geometry"](t)
+end
+
+function jacobian(u::FunctionSpace, ip::IntegrationPoint, t::Number)
+    jacobian(u, ip.xi, t)
+end
+
+function jacobian(u::FunctionSpace, xi)
+    jacobian(u, xi, Inf)
+end
+
+function LinAlg.det(u::FunctionSpace)
+    function detJ(args...)
+        J = jacobian(u, args...)
+        m, n = size(J)
+        return m == n ? det(J) : norm(J)
+    end
+    return detJ
+end
+
+function get_basis(element::Element)
+    return FunctionSpace(element)
+end
+
+Base.(:+)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) + v(args...)
+Base.(:-)(u::FunctionSpace, v::FunctionSpace) = (args...) -> u(args...) - v(args...)
+Base.(:+)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) + v(args...)
+Base.(:-)(u::GradientFunctionSpace, v::GradientFunctionSpace) = (args...) -> u(args...) - v(args...)
+
+
 """ Check does fieldset exist. """
 function Base.haskey(element::Element, what)
     haskey(element.fields, symbol(what))
 end
 
-"""
-Get jacobian of element evaluated at point ξ on element in reference configuration.
-
-Parameters
-----------
-element :: Element
-xi :: Vector
-    spatial coordinate
-time :: Float64
-    temporal coordinate
-geometry_field :: optional
-
-Returns
--------
-Vector or Matrix
-    depending on element dimension
-
-"""
-function get_jacobian(element::Element, xi, time, geometry_field="geometry")
-    dinterpolate(element, geometry_field, xi, time)
-end
-
-""" Evaluate partial derivatives of basis, dbasis/dX, at some time t"""
-function get_dbasisdX(el::Element, xi, t)
-    dbasisdxi = get_dbasisdxi(el, xi)
-    J = get_jacobian(el, xi, t)
-    dbasisdxi*inv(J)
-end
-
-
-
-
 
 # FIXME: These two needs integration -- maybe not in elements.jl ..?
 """
diff --git a/src/equations.jl b/src/equations.jl
index 8df27ec..5ee0de6 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -3,68 +3,174 @@
 
 abstract Equation
 
+abstract Assembly
+
+""" Local element assembly. """
+type LocalAssembly <: Assembly
+    ndofs :: Int
+    mass_matrix :: Matrix
+    stiffness_matrix :: Matrix
+    force_vector :: Matrix
+    potential_energy# :: Union{Array, Float64}
+    residual_vector :: Vector
+end
+
+function LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector::Matrix)
+    LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector[:])
+end
+
+""" Initialize workspace for local assembly. """
+function LocalAssembly(equation::Equation)
+    ndofs = size(equation)
+    mass_matrix = zeros(ndofs, ndofs)
+    stiffness_matrix = zeros(ndofs, ndofs)
+    force_vector = zeros(ndofs, 1)
+    potential_energy = 0.0
+    residual_vector = zeros(ndofs)
+    return LocalAssembly(ndofs, mass_matrix, stiffness_matrix, force_vector,
+                         potential_energy, residual_vector)
+end
+
+function initialize_local_assembly(equation::Equation)
+    LocalAssembly(equation)
+end
+
+function initialize_local_assembly(equation::Equation, assembly::LocalAssembly)
+    if size(equation) != assembly.ndofs
+        # if problem size changes, automatically initialize new work space
+        return initialize_local_assembly(equation)
+    end
+    # otherwise, empty workspace ready for next iteration
+    fill!(assembly.mass_matrix, 0.0)
+    fill!(assembly.stiffness_matrix, 0.0)
+    fill!(assembly.force_vector, 0.0)
+    assembly.potential_energy = 0.0
+    fill!(assembly.residual_vector, 0.0)
+    return assembly
+end
+function initialize_local_assembly(assembly::LocalAssembly, equation::Equation)
+    initialize_local_assembly(equation, assembly)
+end
+
 function get_unknown_field_name(equation::Equation)
     eqtype = typeof(equation)
     error("define get_unknown_field_name for this equation type $eqtype")
 end
 
-
-has_lhs(eq::Equation) = false
-get_lhs(eq::Equation, xi) = nothing
-has_rhs(eq::Equation) = false
-get_rhs(eq::Equation, xi) = nothing
-get_element(eq::Equation) = eq.element
-get_integration_points(eq::Equation) = eq.integration_points
-
-# couple convenient functions -- could make weak form definition easier
-get_connectivity(eq::Equation) = get_connectivity(get_element(eq))
-get_basis(eq::Equation, ip::IntegrationPoint) = get_basis(get_element(eq), ip.xi)
-get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq), ip.xi)
-interpolate(eq::Equation, field::Union{ASCIIString, Symbol}, ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
-integrate_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
-integrate_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
-get_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
-get_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
+has_mass_matrix(equation::Equation) = false
+get_mass_matrix(equation::Equation, ip, time) = nothing
+has_stiffness_matrix(equation::Equation) = false
+get_stiffness_matrix(equation::Equation, ip, time) = nothing
+has_force_vector(equation::Equation) = false
+get_force_vector(equation::Equation, ip, time) = nothing
+has_residual_vector(equation::Equation) = false
+get_residual_vector(equation::Equation, ip, time) = nothing
+has_potential_energy(equation::Equation) = false
+get_potential_energy(equation::Equation, ip, time) = nothing
+get_element(equation::Equation) = equation.element
+get_number_of_dofs(equation::Equation) = nothing
+get_integration_points(equation::Equation) = equation.integration_points
 
 
-"""
-Return determinant of Jacobian for numerical integration.
-"""
-function get_detJ(eq::Equation, ip::IntegrationPoint, t::Float64)
-    el = get_element(eq)
-    get_detJ(el, ip, t)
-end
-function get_detJ(el::Element, ip::IntegrationPoint, t::Float64)
-    get_detJ(el, ip.xi, t)
-end
-function get_detJ(el::Element, xi::Vector, t::Float64)
-    J = get_jacobian(el, xi, t)
-    s = size(J)
-    return s[1] == s[2] ? det(J) : norm(J)
-end
+""" Return a local assembly for element. """
+function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation, time::Number=Inf)
 
-"""
-Integrate f over element
+    initialize_local_assembly(assembly, equation) # zero all
 
-Parameters
-----------
-eq::Equation
+    element = get_element(equation)
+    basis = get_basis(element)
+    detJ = det(basis)
+    field_name = get_unknown_field_name(equation)
 
-f::Function
-    Function to integrate
-"""
-function integrate(eq::Equation, f::Function, t::Float64)
-    target = []
-    for ip in get_integration_points(eq)
-        push!(target, ip.weight*f(eq, ip, t)*get_detJ(eq, ip, t))
+    # 1. if equations are defined we just integrate them
+    if has_mass_matrix(equation) || has_stiffness_matrix(equation) || has_force_vector(equation)
+        for ip in get_integration_points(equation)
+            s = ip.weight*detJ(ip)
+            if has_mass_matrix(equation)
+                assembly.mass_matrix += s*get_mass_matrix(equation, ip, time)
+            end
+            if has_stiffness_matrix(equation)
+                assembly.stiffness_matrix += s*get_stiffness_matrix(equation, ip, time)
+            end
+            if has_force_vector(equation)
+                assembly.force_vector += s*get_force_vector(equation, ip, time)[:]
+            end
+            # external loads -- if any nodal loads is defined add to force vector
+            if haskey(element, "$field_name nodal load")
+                assembly.force_vector += element["$field_name nodal load"](time)[:]
+            end
+        end
     end
-    return sum(target)
+
+    # 2. variational / energy form - user has defined some potential energy / variational form
+    if has_potential_energy(equation)
+        field_name = get_unknown_field_name(equation)
+        element = get_element(equation)
+        field = element[field_name](time)
+        function potential_energy(data::Vector)
+            # calculate potential energy for some setting. this is needed by forwarddiff
+            assembly.potential_energy = 0.0
+            df = similar(field, data)
+            # integrate potential energy
+            for ip in get_integration_points(equation)
+                dw = get_potential_energy(equation, ip, time; variation=df)
+                assembly.potential_energy += ip.weight * dw * detJ(ip)
+            end
+            # external energy -- if any nodal loads is defined, decrease from potential energy
+            if haskey(element, "$field_name nodal load")
+                P = element["$field_name nodal load"](time)
+                assembly.potential_energy -= dot(P[:], df[:])
+            end
+            if isa(assembly.potential_energy, Array)
+                return assembly.potential_energy[1]
+            end
+            return assembly.potential_energy
+        end
+        hessian, allresults = ForwardDiff.hessian(potential_energy, field[:],
+                                                  AllResults, cache=autodiffcache)
+        assembly.stiffness_matrix += hessian
+        assembly.force_vector -= ForwardDiff.gradient(allresults)  # <--- minus explained in tutorial
+        assembly.potential_energy = ForwardDiff.value(allresults)
+    end
+
+    # 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v
+    if has_residual_vector(equation)
+        field_name = get_unknown_field_name(equation)
+        element = get_element(equation)
+        field = element[field_name](time)
+        function residual_vector(data::Vector)
+            fill!(assembly.residual_vector, 0.0)
+            df = similar(field, data)
+            # integrate W
+            for ip in get_integration_points(equation)
+                dr = get_residual_vector(equation, ip, time; variation=df)
+                assembly.residual_vector += ip.weight*dr*detJ(ip)
+            end
+            # external loads -- if any nodal loads is defined, remove from residual
+            if haskey(element, "$field_name nodal load")
+                assembly.residual_vector -= element["$field_name nodal load"](time)[:]
+            end
+            return assembly.residual_vector
+        end
+        jacobian, allresults = ForwardDiff.jacobian(residual_vector, field[:],
+                                                    AllResults, cache=autodiffcache)
+        assembly.stiffness_matrix += jacobian
+        assembly.force_vector -= ForwardDiff.value(allresults)  # <-- minus explained in tutorial
+    end
+
 end
 
+function calculate_local_assembly!(equation::Equation, assembly::LocalAssembly, time::Number=Inf)
+    calculate_local_assembly!(assembly, equation)
+end
+
+
+""" Get global degrees of freedom for this element. """
 function get_global_dofs(eq::Equation)
     eq.global_dofs
 end
 
+""" Set global degrees of freedom for this element. """
 function set_global_dofs!(eq::Equation, dofs)
     eq.global_dofs = dofs
 end
diff --git a/src/interpolate.jl b/src/interpolate.jl
index b156823..c1dbfda 100644
--- a/src/interpolate.jl
+++ b/src/interpolate.jl
@@ -59,8 +59,7 @@ function interpolate(basis::Basis, field::Field, ip::IntegrationPoint)
     interpolate(basis, field, ip.xi)
 end
 
-function dinterpolate(N::Basis, u::Field, xi::Array{Float64, 1})
-    dN = diff(N)
-    dN(xi)*u
+function dinterpolate(basis::Basis, u::Field, xi::Array{Float64, 1})
+    basis.dbasisdxi(xi)*u
 end
 
diff --git a/src/types.jl b/src/types.jl
index 8e080a7..f9330f5 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -70,12 +70,17 @@ JuliaFEM.Field{Array{Array{T,1},1}}(0.5,1,Array{T,1}[[1.0,1.0],[1.0,1.0]])
 
 """
 function Base.similar(field::Field, data::Vector)
+    fdim = round(Int, length(data)/length(field))  # dimension of field variable
+    if fdim == 1
+        new_field = Field(field.time, data)
+        return new_field
+    end
     new_field = Field(field.time, similar(field.values))
-    data = reshape(data, round(Int, length(data)/length(field)), length(field))
+    data = reshape(data, fdim, length(field))
     for i=1:length(new_field)
         new_field.values[i] = data[:,i]
     end
-    new_field
+    return new_field
 end
 
 
@@ -117,7 +122,6 @@ function Base.endof(fieldset::FieldSet)
 end
 
 
-
 """ Basis function. """
 type Basis
     basis :: Function
@@ -128,9 +132,9 @@ function Basis(basis)
     Basis(basis, ForwardDiff.jacobian(basis))
 end
 """ Get partial derivative of basis function. """
-diff(h::Basis) = h.dbasisdxi
-derivative(h::Basis) = h.dbasisdxi
-
+function grad(basis::Basis)
+    (ip) -> basis.dbasisdxi(ip.xi)
+end
 
 
 """
@@ -158,7 +162,10 @@ end
 
 # convenient functions -- maybe this is not correct place for them
 """ Evaluate basis function in point ξ. """
-call(b::Basis, xi) = b.basis(xi)
+call(b::Basis, xi::Vector) = b.basis(xi)
+call(b::Basis, ip::IntegrationPoint) = b.basis(ip.xi)
+Base.(:*)(basis::Basis, fs::FieldSet) = (xi, t) -> basis(xi)*fs(t)
+
 #""" Interpolate field (h*f)(ξ) """
 #Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
 #""" Interpolate from set of fields with basis b, i.e. f(t) = b(t)*[f1, f2] """
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 170357a..c3b1f48 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -2,7 +2,7 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 using FactCheck
-using JuliaFEM: Element, Basis, FieldSet
+using JuliaFEM: Element, Basis, Field, FieldSet, FunctionSpace
 
 """ Prototype element
 
@@ -51,3 +51,35 @@ facts("test adding fieldsets and fields to element") do
     @fact fields[2] --> field2
 end
 
+facts("interpolation of fields in some function space") do
+
+    element = MockElement([1, 2, 3, 4])
+    fieldset1 = FieldSet("geometry", [Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])
+    fieldset2 = FieldSet("constant scalar field", [Field(0.0, 1.0)])
+    fieldset3 = FieldSet("scalar field", [Field(0.0, [1.0, 2.0, 3.0, 4.0])])
+    fieldset4 = FieldSet("vector field 1", [Field(0.0, Vector[[1.0], [2.0], [3.0], [4.0]])])
+    fieldset5 = FieldSet("vector field 2", [Field(0.0, Vector[[1.0, 5.0], [2.0, 6.0], [3.0, 7.0], [4.0, 8.0]])])
+    fieldset6 = FieldSet("vector field 3", [Field(0.0, Vector[[1.0, 5.0, 9.0], [2.0, 6.0, 10.0], [3.0, 7.0, 11.0], [4.0, 8.0, 12.0]])])
+    fieldset7 = FieldSet("tensor field 1", [Field(0.0, Matrix[[1.0 5.0; 9.0 13.0], [2.0 6.0; 10.0 14.0], [3.0 7.0; 11.0 15.0], [4.0 8.0; 12.0 16.0]])])
+
+    push!(element, fieldset1)
+    push!(element, fieldset2)
+    push!(element, fieldset3)
+    push!(element, fieldset4)
+    push!(element, fieldset5)
+    push!(element, fieldset6)
+    push!(element, fieldset7)
+
+    xi = [0.0, 0.0]
+    t = 0.0
+    u = FunctionSpace(element)
+    v = FunctionSpace(element)
+
+    @fact v("constant scalar field", xi, t) --> 1.0
+    @fact v("scalar field", xi, t) --> 1/4*(1+2+3+4)
+    @fact v("vector field 1", xi, t) --> [1/4*(1+2+3+4)]
+    @fact v("vector field 2", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8]
+    @fact v("vector field 3", xi, t) --> 1/4*[1+2+3+4, 5+6+7+8, 9+10+11+12]
+    @fact v("tensor field 1", xi, t) --> 1/4*[1+2+3+4 5+6+7+8; 9+10+11+12 13+14+15+16]
+end
+
diff --git a/test/test_interpolation.jl b/test/test_interpolation.jl
new file mode 100644
index 0000000..1b4dade
--- /dev/null
+++ b/test/test_interpolation.jl
@@ -0,0 +1,50 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
+using FactCheck
+
+element = Quad4([1, 2, 3, 4])
+
+geometry_field = Field(0.0, Vector[])  # Create empty field at time t=0.0
+push!(geometry_field, [ 0.0, 0.0])  # push some values for field
+push!(geometry_field, [ 1.0, 0.0])
+push!(geometry_field, [ 1.0, 1.0])
+push!(geometry_field, [ 0.0, 1.0])
+geometry_fieldset = FieldSet("geometry")  # create fieldset "geometry"
+push!(geometry_fieldset, geometry_field)  # add field to fieldset
+push!(element, geometry_fieldset) # add fieldset to element
+
+temperature_fieldset = FieldSet("temperature")
+push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))
+push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))
+push!(element, temperature_fieldset)
+
+displacement_fieldset = FieldSet("displacement")
+push!(displacement_fieldset, Field(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))
+push!(displacement_fieldset, Field(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]]))
+push!(element, displacement_fieldset)
+
+facts("basic continuum interpolations") do
+    # from my old home works
+    basis = get_basis(element)
+    dbasis = grad(basis)
+    @fact basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0) --> [9/16, 1/2]
+    gradu = dbasis("displacement", [0.0, 0.0], 1.0)
+    epsilon = 1/2*(gradu + gradu')
+    rotation = 1/2*(gradu - gradu')
+    X = basis("geometry", [0.0, 0.0], 1.0)
+    k = 0.25
+    epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
+    rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
+    @fact epsilon --> roughly(epsilon_wanted)
+    @fact rotation --> roughly(rotation_wanted)
+    F = I + gradu
+    @fact F --> [X[2]*k+1 X[1]*k; 0 1]
+    C = F'*F
+    @fact C --> [(X[2]*k+1)^2  (X[2]*k+1)*X[1]*k; (X[2]*k+1)*X[1]*k  X[1]^2*k^2+1]
+    E = 1/2*(F'*F - I)
+    @fact E --> [1/2*(X[2]*k + 1)^2-1/2  1/2*(X[2]*k+1)*X[1]*k; 1/2*(X[2]*k + 1)*X[1]*k  1/2*X[1]^2*k^2]
+    U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I)
+    #@fact U --> roughly([1.24235 0.13804; 0.13804 1.02149])
+end
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson1.comm b/verification/2015-10-09-verification-of-poisson-equation/poisson1.comm
new file mode 100644
index 0000000..e22d842
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson1.comm
@@ -0,0 +1,39 @@
+# heat source f=12
+
+DEBUT()
+
+MAIL = LIRE_MAILLAGE()
+
+MO = AFFE_MODELE(
+    MAILLAGE=MAIL,
+    AFFE = _F(MAILLE='E1', PHENOMENE='THERMIQUE', MODELISATION='PLAN'))
+
+MAT = DEFI_MATERIAU(
+    THER = _F(LAMBDA = 6.0))
+
+CHMAT = AFFE_MATERIAU(
+    MAILLAGE = MAIL,
+    AFFE = _F(MAILLE = 'E1', MATER = MAT))
+
+BC = AFFE_CHAR_THER( # Dirichlet boundary condition on 0 <= X <= 1, Y = 1
+    MODELE = MO,
+    TEMP_IMPO = (_F(NOEUD = ('N3','N4'), TEMP=0)))
+
+# Heat source f = 12
+LO = AFFE_CHAR_THER(
+    MODELE = MO,
+    SOURCE = _F(MAILLE='E1', SOUR=12.0))
+
+RESU = THER_LINEAIRE(
+    MODELE=MO,
+    CHAM_MATER=CHMAT,
+    EXCIT=(
+        _F(CHARGE=BC),
+        _F(CHARGE=LO)))
+
+IMPR_RESU(
+    MODELE = MO,
+    FORMAT = 'RESULTAT',
+    RESU = _F(RESULTAT = RESU))
+
+FIN()
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson1.mail b/verification/2015-10-09-verification-of-poisson-equation/poisson1.mail
new file mode 100644
index 0000000..2df78f1
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson1.mail
@@ -0,0 +1,13 @@
+
+ COOR_2D
+ N1          0.0 0.0
+ N2          1.0 0.0
+ N3          1.0 1.0
+ N4          0.0 1.0
+ FINSF
+
+ QUAD4
+ E1         N1 N2 N3 N4
+ FINSF
+
+ FIN
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson1.resu b/verification/2015-10-09-verification-of-poisson-equation/poisson1.resu
new file mode 100644
index 0000000..d41d6dd
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson1.resu
@@ -0,0 +1,92 @@
+                                                                                
+                                                                                
+              -- CODE_ASTER -- VERSION : EXPLOITATION (stable) --               
+                                                                                
+                          Version 11.4.0 du 05/06/2013                          
+                         Copyright EDF R&D 1991 - 2015                          
+                                                                                
+                    Exécution du : Fri Oct 23 23:35:34 2015                     
+                       Nom de la machine : jukka-desktop                        
+                              Architecture : 64bit                              
+                          Type de processeur : x86_64                           
+                Système d'exploitation : Linux 3.13.0-66-generic                
+                        Langue des messages : en (UTF-8)                        
+                                                                                
+   
+   !------------------------------------------------------------------------------------!
+   ! <A> <SUPERVIS2_2>                                                                  !
+   !                                                                                    !
+   !     Vous utilisez une vieille version de Code_Aster.                               !
+   !                                                                                    !
+   !     En mettant à jour votre version, vous bénéficierez des dernières améliorations !
+   !     apportées au code depuis 15 mois.                                              !
+   !     Si vous avez des développements privés, vous risquez d'avoir un travail        !
+   !     important de portage si vous ne suivez pas les mises à jour.                   !
+   !                                                                                    !
+   !                                                                                    !
+   ! Ceci est une alarme. Si vous ne comprenez pas le sens de cette                     !
+   ! alarme, vous pouvez obtenir des résultats inattendus !                             !
+   !------------------------------------------------------------------------------------!
+   
+                           Parallélisme MPI : inactif                           
+                          Parallélisme OpenMP : actif                           
+                        Nombre de processus utilisés : 1                        
+                      Version de la librairie HDF5 : 1.8.8                      
+                      Version de la librairie MED : 3.0.6                       
+                          Librairie MUMPS : installée                           
+                    Version de la librairie SCOTCH : 5.1.10                     
+                  Mémoire limite pour l'exécution : 4096.00 Mo                  
+                   consommée par l'initialisation : 197.46 Mo                   
+                  par les objets du jeu de commandes : 1.15 Mo                  
+                 reste pour l'allocation dynamique : 3897.26 Mo                 
+                Taille limite des fichiers d'échange : 48.00 Go                 
+                                                                                
+
+ --------------------------------------------------------------------------------
+ ASTER 11.04.00 CONCEPT RESU CALCULE LE 23/10/2015 A 23:35:34 DE TYPE EVOL_THER  
+
+
+ ======>
+
+ ------>
+ CHAMP AUX NOEUDS DE NOM SYMBOLIQUE  TEMP                                        
+ NUMERO D'ORDRE: 0 INST:   0.00000000000000E+00                                  
+ NOEUD                 TEMP    
+ N1        1.00000000000000E+00
+ N2        1.00000000000000E+00
+ N3        1.11022302462516E-16
+ N4        1.11022302462516E-16
+ 
+ <I> <FIN> FERMETURE DE LA BASE "GLOBALE" EFFECTUEE.
+
+ <FIN> Arrêt normal dans "FIN".
+ <I> <FIN> ARRET NORMAL DANS "FIN" PAR APPEL A "JEFINI".
+  
+ <I> <FIN> MEMOIRE JEVEUX MINIMALE REQUISE POUR L'EXECUTION :                      20.87 Mo
+ <I> <FIN> MEMOIRE JEVEUX OPTIMALE REQUISE POUR L'EXECUTION :                      27.32 Mo
+ <I> <FIN> MAXIMUM DE MEMOIRE UTILISEE PAR LE PROCESSUS LORS DE L'EXECUTION :     226.60 Mo
+
+ ********************************************************************************
+ * COMMAND                  :       USER :     SYSTEM :   USER+SYS :    ELAPSED *
+ ********************************************************************************
+ * init (jdc)               :       0.15 :       0.02 :       0.17 :       0.17 *
+ *  . compile               :       0.00 :       0.00 :       0.00 :       0.00 *
+ *  . exec_compile          :       0.05 :       0.01 :       0.06 :       0.05 *
+ *  . report                :       0.00 :       0.00 :       0.00 :       0.01 *
+ *  . build                 :       0.00 :       0.00 :       0.00 :       0.00 *
+ * DEBUT                    :       0.01 :       0.02 :       0.03 :       0.03 *
+ * LIRE_MAILLAGE            :       0.00 :       0.00 :       0.00 :       0.01 *
+ * AFFE_MODELE              :       0.00 :       0.00 :       0.00 :       0.00 *
+ * DEFI_MATERIAU            :       0.00 :       0.00 :       0.00 :       0.00 *
+ * AFFE_MATERIAU            :       0.01 :       0.00 :       0.01 :       0.00 *
+ * AFFE_CHAR_THER           :       0.00 :       0.00 :       0.00 :       0.01 *
+ * AFFE_CHAR_THER           :       0.00 :       0.00 :       0.00 :       0.00 *
+ * THER_LINEAIRE            :       0.01 :       0.00 :       0.01 :       0.01 *
+ * IMPR_RESU                :       0.01 :       0.00 :       0.01 :       0.00 *
+ * FIN                      :       0.01 :       0.01 :       0.02 :       0.02 *
+ *  . part Superviseur      :       0.17 :       0.04 :       0.21 :       0.21 *
+ *  . part Fortran          :       0.04 :       0.01 :       0.05 :       0.04 *
+ ********************************************************************************
+ * TOTAL_JOB                :       0.21 :       0.05 :       0.26 :       0.25 *
+ ********************************************************************************
+
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson2.comm b/verification/2015-10-09-verification-of-poisson-equation/poisson2.comm
new file mode 100644
index 0000000..223ef5f
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson2.comm
@@ -0,0 +1,40 @@
+# Flux q=6 on free boundary
+# B1 = boundary element
+
+DEBUT()
+
+MAIL = LIRE_MAILLAGE()
+
+MO = AFFE_MODELE(
+    MAILLAGE=MAIL,
+    AFFE = _F(MAILLE=('B1', 'E1'), PHENOMENE='THERMIQUE', MODELISATION='PLAN'))
+
+MAT = DEFI_MATERIAU(
+    THER = _F(LAMBDA = 6.0))
+
+CHMAT = AFFE_MATERIAU(
+    MAILLAGE = MAIL,
+    AFFE = _F(MAILLE = 'E1', MATER = MAT))
+
+BC = AFFE_CHAR_THER( # Dirichlet boundary condition on 0 <= X <= 1, Y = 1
+    MODELE = MO,
+    TEMP_IMPO = (_F(NOEUD = ('N3','N4'), TEMP=0)))
+
+# Heat flux on free boundary
+LO = AFFE_CHAR_THER(
+    MODELE = MO,
+    FLUX_REP = _F(MAILLE='B1', FLUN=6))
+
+RESU = THER_LINEAIRE(
+    MODELE=MO,
+    CHAM_MATER=CHMAT,
+    EXCIT=(
+        _F(CHARGE=BC),
+        _F(CHARGE=LO)))
+
+IMPR_RESU(
+    MODELE = MO,
+    FORMAT = 'RESULTAT',
+    RESU = _F(RESULTAT = RESU))
+
+FIN()
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson2.mail b/verification/2015-10-09-verification-of-poisson-equation/poisson2.mail
new file mode 100644
index 0000000..acc66e2
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson2.mail
@@ -0,0 +1,17 @@
+
+ COOR_2D
+ N1          0.0 0.0
+ N2          1.0 0.0
+ N3          1.0 1.0
+ N4          0.0 1.0
+ FINSF
+
+ QUAD4
+ E1         N1 N2 N3 N4
+ FINSF
+
+ SEG2
+ B1         N1 N2
+ FINSF
+
+ FIN
diff --git a/verification/2015-10-09-verification-of-poisson-equation/poisson2.resu b/verification/2015-10-09-verification-of-poisson-equation/poisson2.resu
new file mode 100644
index 0000000..a814c1d
--- /dev/null
+++ b/verification/2015-10-09-verification-of-poisson-equation/poisson2.resu
@@ -0,0 +1,92 @@
+                                                                                
+                                                                                
+              -- CODE_ASTER -- VERSION : EXPLOITATION (stable) --               
+                                                                                
+                          Version 11.4.0 du 05/06/2013                          
+                         Copyright EDF R&D 1991 - 2015                          
+                                                                                
+                    Exécution du : Fri Oct 23 23:52:48 2015                     
+                       Nom de la machine : jukka-desktop                        
+                              Architecture : 64bit                              
+                          Type de processeur : x86_64                           
+                Système d'exploitation : Linux 3.13.0-66-generic                
+                        Langue des messages : en (UTF-8)                        
+                                                                                
+   
+   !------------------------------------------------------------------------------------!
+   ! <A> <SUPERVIS2_2>                                                                  !
+   !                                                                                    !
+   !     Vous utilisez une vieille version de Code_Aster.                               !
+   !                                                                                    !
+   !     En mettant à jour votre version, vous bénéficierez des dernières améliorations !
+   !     apportées au code depuis 15 mois.                                              !
+   !     Si vous avez des développements privés, vous risquez d'avoir un travail        !
+   !     important de portage si vous ne suivez pas les mises à jour.                   !
+   !                                                                                    !
+   !                                                                                    !
+   ! Ceci est une alarme. Si vous ne comprenez pas le sens de cette                     !
+   ! alarme, vous pouvez obtenir des résultats inattendus !                             !
+   !------------------------------------------------------------------------------------!
+   
+                           Parallélisme MPI : inactif                           
+                          Parallélisme OpenMP : actif                           
+                        Nombre de processus utilisés : 1                        
+                      Version de la librairie HDF5 : 1.8.8                      
+                      Version de la librairie MED : 3.0.6                       
+                          Librairie MUMPS : installée                           
+                    Version de la librairie SCOTCH : 5.1.10                     
+                  Mémoire limite pour l'exécution : 4096.00 Mo                  
+                   consommée par l'initialisation : 197.46 Mo                   
+                  par les objets du jeu de commandes : 1.14 Mo                  
+                 reste pour l'allocation dynamique : 3897.26 Mo                 
+                Taille limite des fichiers d'échange : 48.00 Go                 
+                                                                                
+
+ --------------------------------------------------------------------------------
+ ASTER 11.04.00 CONCEPT RESU CALCULE LE 23/10/2015 A 23:52:48 DE TYPE EVOL_THER  
+
+
+ ======>
+
+ ------>
+ CHAMP AUX NOEUDS DE NOM SYMBOLIQUE  TEMP                                        
+ NUMERO D'ORDRE: 0 INST:   0.00000000000000E+00                                  
+ NOEUD                 TEMP    
+ N1        1.00000000000000E+00
+ N2        1.00000000000000E+00
+ N3        5.55111512312578E-17
+ N4        5.55111512312578E-17
+ 
+ <I> <FIN> FERMETURE DE LA BASE "GLOBALE" EFFECTUEE.
+
+ <FIN> Arrêt normal dans "FIN".
+ <I> <FIN> ARRET NORMAL DANS "FIN" PAR APPEL A "JEFINI".
+  
+ <I> <FIN> MEMOIRE JEVEUX MINIMALE REQUISE POUR L'EXECUTION :                      21.00 Mo
+ <I> <FIN> MEMOIRE JEVEUX OPTIMALE REQUISE POUR L'EXECUTION :                      27.32 Mo
+ <I> <FIN> MAXIMUM DE MEMOIRE UTILISEE PAR LE PROCESSUS LORS DE L'EXECUTION :     226.59 Mo
+
+ ********************************************************************************
+ * COMMAND                  :       USER :     SYSTEM :   USER+SYS :    ELAPSED *
+ ********************************************************************************
+ * init (jdc)               :       0.16 :       0.01 :       0.17 :       0.17 *
+ *  . compile               :       0.00 :       0.00 :       0.00 :       0.00 *
+ *  . exec_compile          :       0.04 :       0.01 :       0.05 :       0.06 *
+ *  . report                :       0.01 :       0.00 :       0.01 :       0.00 *
+ *  . build                 :       0.00 :       0.00 :       0.00 :       0.00 *
+ * DEBUT                    :       0.02 :       0.01 :       0.03 :       0.04 *
+ * LIRE_MAILLAGE            :       0.00 :       0.00 :       0.00 :       0.00 *
+ * AFFE_MODELE              :       0.00 :       0.00 :       0.00 :       0.00 *
+ * DEFI_MATERIAU            :       0.00 :       0.00 :       0.00 :       0.00 *
+ * AFFE_MATERIAU            :       0.00 :       0.00 :       0.00 :       0.01 *
+ * AFFE_CHAR_THER           :       0.01 :       0.00 :       0.01 :       0.00 *
+ * AFFE_CHAR_THER           :       0.00 :       0.00 :       0.00 :       0.00 *
+ * THER_LINEAIRE            :       0.01 :       0.01 :       0.02 :       0.02 *
+ * IMPR_RESU                :       0.00 :       0.00 :       0.00 :       0.00 *
+ * FIN                      :       0.01 :       0.01 :       0.02 :       0.03 *
+ *  . part Superviseur      :       0.18 :       0.03 :       0.21 :       0.23 *
+ *  . part Fortran          :       0.03 :       0.02 :       0.05 :       0.04 *
+ ********************************************************************************
+ * TOTAL_JOB                :       0.21 :       0.05 :       0.26 :       0.27 *
+ ********************************************************************************
+