diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb
index 3c9077f..d791da8 100644
--- a/notebooks/2015-08-29-developing-juliafem.ipynb
+++ b/notebooks/2015-08-29-developing-juliafem.ipynb
@@ -8,7 +8,19 @@
     "\n",
     "Author(s): Jukka Aho\n",
     "\n",
-    "**Abstract**: General developer notes."
+    "**Abstract**: Developer notes. In this notebook we give general guidelines how to implement things to JuliaFEM.\n",
+    "\n",
+    "In short, we have\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",
+    "\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",
+    "\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."
    ]
   },
   {
@@ -54,21 +66,14 @@
    "source": [
     "## Developing own element\n",
     "\n",
-    "Finite element definition, from [FEniCS-book](https://bitbucket.org/fenics-project/fenics-book/src/7d3a80e7dda0fc279c7964dc6000d57942f11eb3/fenicsbook.cls?at=master) [Ciarlet, 2002]:\n",
-    "\n",
-    "-  the domain $T$ is a bounded, closed subset of $\\mathbb{R}^d$ (for $d = 1, 2, 3, \\dots$) with nonempty interior and piecewise smooth boundary;\n",
-    "- the space $\\mathcal{V} = \\mathcal{V}(T)$ is a finite dimensional function space on $T$ of dimension $n$;\n",
-    "- the set of degrees of freedom (nodes) $\\mathcal{L} = \\{\\ell_1, \\ell_2,\\ldots, \\ell_{n}\\}$ is a basis for the dual space $\\mathcal{V}'$; that is, the space of bounded linear functionals on $\\mathcal{V}$.\n",
-    "\n",
-    "We extend this definition so that domain $T$ can also be empty.\n",
+    "Element contains basis functions and one or several fieldsets.\n",
     "\n",
     "Minimum requirements for element:\n",
-    "- subclass from Element, if not wanting to implement everything by youself\n",
+    "- subclass it from Element, unless want to implement everything by youself\n",
     "- define basis and partial derivatives of it, because we need to interpolate over it\n",
     "- give connectivity information, how this element is connected to other elements\n",
-    "- create proper constructor (see example).\n",
     "\n",
-    "Test the element using ``test_element`` function. It it passes, then element implementation should be fine. As an example, we define 4 node quadrilateral element using linear Lagrange basis. We really don't care much how element is implemented as long it's interface is constructed with some rules. The interface is tested using `test_element` and it also gives information how to fix element if something is missing."
+    "Test the element using ``test_element`` function. It it passes, the element implementation should be fine."
    ]
   },
   {
@@ -79,14 +84,14 @@
    },
    "outputs": [],
    "source": [
-    "using JuliaFEM: Element, Field, Basis"
+    "using JuliaFEM: Element, Field, FieldSet, Basis"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Here's the implementation for element:"
+    "Here's the basic implementation for element:"
    ]
   },
   {
@@ -100,7 +105,7 @@
     "type MyQuad4 <: Element\n",
     "    connectivity :: Array{Int, 1}\n",
     "    basis :: Basis\n",
-    "    fields :: Dict{Symbol, Array{Field, 1}}\n",
+    "    fields :: Dict{Symbol, FieldSet}\n",
     "end"
    ]
   },
@@ -108,7 +113,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Default constructor takes only connectivity information as input argument. Extra arguments may be passed using `args...`"
+    "Default constructor takes only connectivity information as input argument. Most important thing here is to define `Basis` which is used to interpolate fields."
    ]
   },
   {
@@ -130,13 +135,19 @@
     }
    ],
    "source": [
-    "function MyQuad4(connectivity, args...)\n",
-    "    h = Basis((xi) ->\n",
-    "    [(1-xi[1])*(1-xi[2])/4\n",
-    "     (1+xi[1])*(1-xi[2])/4\n",
-    "     (1+xi[1])*(1+xi[2])/4\n",
-    "     (1-xi[1])*(1+xi[2])/4])\n",
-    "    MyQuad4(connectivity, h, Dict())\n",
+    "function MyQuad4(connectivity)\n",
+    "    h(xi) = [\n",
+    "        (1-xi[1])*(1-xi[2])/4\n",
+    "        (1+xi[1])*(1-xi[2])/4\n",
+    "        (1+xi[1])*(1+xi[2])/4\n",
+    "        (1-xi[1])*(1+xi[2])/4]\n",
+    "    dh(xi) = [\n",
+    "        -(1-xi[2])/4.0   -(1-xi[1])/4.0\n",
+    "         (1-xi[2])/4.0   -(1+xi[1])/4.0\n",
+    "         (1+xi[2])/4.0    (1+xi[1])/4.0\n",
+    "        -(1+xi[2])/4.0    (1-xi[1])/4.0]\n",
+    "    basis = Basis(h, dh)\n",
+    "    MyQuad4(connectivity, basis, Dict())\n",
     "end"
    ]
   },
@@ -144,7 +155,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Some basic charasteristics like number of basis funcitons and dimension:"
+    "Also some basic information like number of basis functions and element dimension is needed."
    ]
   },
   {
@@ -188,14 +199,14 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "07-Oct 09:28:37:INFO:root:Testing element MyQuad4\n",
-      "07-Oct 09:28:37:INFO:root:number of basis functions in this element: 4\n",
-      "07-Oct 09:28:37:INFO:root:Initializing element\n",
-      "07-Oct 09:28:37:INFO:root:Element dimension: 2\n",
-      "07-Oct 09:28:37:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n",
-      "07-Oct 09:28:37:INFO:root:Pushing field to element.\n",
-      "07-Oct 09:28:38:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
-      "07-Oct 09:28:38:INFO:root:Value: 2.5\n"
+      "09-Oct 01:41:19:INFO:root:Testing element MyQuad4\n",
+      "09-Oct 01:41:19:INFO:root:number of basis functions in this element: 4\n",
+      "09-Oct 01:41:19:INFO:root:Initializing element\n",
+      "09-Oct 01:41:19:INFO:root:Element dimension: 2\n",
+      "09-Oct 01:41:19:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n",
+      "09-Oct 01:41:19:INFO:root:Pushing field to element.\n",
+      "09-Oct 01:41:19:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
+      "09-Oct 01:41:19:INFO:root:Value: 2.5\n"
      ]
     },
     {
@@ -218,7 +229,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "If `test_element` passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:"
+    "If `test_element` passes, elements *interface* should be well defined. After building element, one can interpolate things in it. Here we create three new fieldsets `temperature`, `geometry` and `heat coefficient`, add some values for them in time $t=0.0$ and $t=1.0$ and interpolate:"
    ]
   },
   {
@@ -244,18 +255,21 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "07-Oct 09:28:38:INFO:root:Element MyQuad4 passed tests.\n"
+      "09-Oct 01:41:19:INFO:root:Element MyQuad4 passed tests.\n"
      ]
     }
    ],
    "source": [
-    "using JuliaFEM: new_field!, push_field!, interpolate, dinterpolate\n",
+    "using JuliaFEM: new_fieldset!, add_field!, interpolate, dinterpolate\n",
     "el1 = MyQuad4([1, 2, 3, 4])\n",
-    "new_field!(el1, :temperature, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
-    "push_field!(el1, :temperature, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
-    "new_field!(el1, :Geometry, Field(0.0, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]))\n",
-    "new_field!(el1, \"heat coefficient\", Field(0.0, 2))\n",
-    "push_field!(el1, \"heat coefficient\", Field(1.0, 3))"
+    "new_fieldset!(el1, \"temperature\")\n",
+    "add_field!(el1, \"temperature\", Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
+    "add_field!(el1, \"temperature\", Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
+    "new_fieldset!(el1, \"geometry\")\n",
+    "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]))\n",
+    "new_fieldset!(el1, \"heat coefficient\")\n",
+    "add_field!(el1, \"heat coefficient\", Field(0.0, 2))\n",
+    "add_field!(el1, \"heat coefficient\", Field(1.0, 3))"
    ]
   },
   {
@@ -278,7 +292,7 @@
    ],
    "source": [
     "# temperature at the middle poinf of the element, 1/4*(1+2+3+4) at t=0.5\n",
-    "interpolate(el1, :temperature, [0.0, 0.0], 0.5)"
+    "interpolate(el1, \"temperature\", [0.0, 0.0], 0.5)"
    ]
   },
   {
@@ -304,7 +318,7 @@
    ],
    "source": [
     "# geometry midpoint of element\n",
-    "interpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
+    "interpolate(el1, \"geometry\", [0.0, 0.0], 0.0)"
    ]
   },
   {
@@ -329,7 +343,8 @@
     }
    ],
    "source": [
-    "dinterpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
+    "# interpolate derivatives works too\n",
+    "dinterpolate(el1, \"geometry\", [0.0, 0.0], 0.0)"
    ]
   },
   {
@@ -342,7 +357,7 @@
     {
      "data": {
       "text/plain": [
-       "([0.625,0.625,0.625,0.625],[0.75,0.75,0.75,0.75])"
+       "2.5"
       ]
      },
      "execution_count": 11,
@@ -352,7 +367,29 @@
    ],
    "source": [
     "# interpolating scalar -> scalar.\n",
-    "interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5), interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)"
+    "interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "3"
+      ]
+     },
+     "execution_count": 12,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)"
    ]
   },
   {
@@ -368,7 +405,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Developing own formulation\n",
+    "## Developing own equation\n",
     "\n",
     "Let's consider a Laplace equation\n",
     "\\begin{align}\n",
@@ -384,6 +421,7 @@
     "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",
@@ -396,15 +434,28 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 13,
    "metadata": {
     "collapsed": false
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "get_unknown_field_name (generic function with 1 method)"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "using JuliaFEM: Equation, IntegrationPoint, Quad4\n",
     "\n",
-    "abstract Heat <: Equation"
+    "abstract Heat <: Equation\n",
+    "\n",
+    "get_unknown_field_name(eq::Heat) = symbol(\"temperature\")"
    ]
   },
   {
@@ -416,7 +467,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 14,
    "metadata": {
     "collapsed": false
    },
@@ -441,7 +492,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 15,
    "metadata": {
     "collapsed": false
    },
@@ -452,7 +503,7 @@
        "DC2D4"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 15,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -464,7 +515,7 @@
     "        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",
-    "    new_field!(el, :temperature) # assign new field \"temperature\" to element\n",
+    "    new_fieldset!(el, \"temperature\") # assign new field \"temperature\" to element\n",
     "    DC2D4(el, integration_points, [])\n",
     "end"
    ]
@@ -478,7 +529,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 16,
    "metadata": {
     "collapsed": false
    },
@@ -489,7 +540,7 @@
        "has_lhs (generic function with 2 methods)"
       ]
      },
-     "execution_count": 15,
+     "execution_count": 16,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -503,8 +554,7 @@
     "function JuliaFEM.get_lhs(eq::DC2D4, ip, t)\n",
     "    el = get_element(eq)\n",
     "    dNdX = get_dbasisdX(el, ip.xi, t)\n",
-    "    fld = el[\"temperature thermal conductivity\"](t)\n",
-    "    hc = sum(el(ip.xi) * fld)\n",
+    "    hc = interpolate(el, \"temperature thermal conductivity\", ip.xi, t)\n",
     "    return dNdX*hc*dNdX'\n",
     "end\n",
     "JuliaFEM.has_lhs(eq::DC2D4) = true"
@@ -519,7 +569,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 17,
    "metadata": {
     "collapsed": false
    },
@@ -534,47 +584,20 @@
        " -1.0  -2.0  -1.0   4.0"
       ]
      },
-     "execution_count": 16,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
-    "el = Quad4([1, 2, 3, 4])\n",
-    "new_field!(el, :Geometry, Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]))\n",
-    "new_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
-    "eq = DC2D4(el)\n",
-    "integrate_lhs(eq, 1.0)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 17,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "true"
-      ]
-     },
      "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "JuliaFEM.has_lhs(eq)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If rhs or lhs is not defined, integration returns nothing."
+    "using JuliaFEM: integrate, integrate_lhs, integrate_rhs, new_fieldset!, add_field!\n",
+    "el = Quad4([1, 2, 3, 4])\n",
+    "new_fieldset!(el, \"geometry\")\n",
+    "add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]))\n",
+    "new_fieldset!(el, \"temperature thermal conductivity\")\n",
+    "add_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
+    "eq = DC2D4(el)\n",
+    "integrate_lhs(eq, 1.0)"
    ]
   },
   {
@@ -596,14 +619,14 @@
     }
    ],
    "source": [
-    "isa(integrate_rhs(eq, 1.0), Void)"
+    "JuliaFEM.has_lhs(eq)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Next heat flux on boundary:"
+    "If rhs or lhs is not defined, integration returns nothing."
    ]
   },
   {
@@ -616,7 +639,7 @@
     {
      "data": {
       "text/plain": [
-       "has_rhs (generic function with 2 methods)"
+       "true"
       ]
      },
      "execution_count": 19,
@@ -624,6 +647,35 @@
      "output_type": "execute_result"
     }
    ],
+   "source": [
+    "isa(integrate_rhs(eq, 1.0), Void)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next heat flux on boundary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "has_rhs (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 20,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "using JuliaFEM: get_basis, Seg2\n",
     "\n",
@@ -639,7 +691,7 @@
     "function DC2D2(el::Seg2)\n",
     "    integration_points = [\n",
     "        IntegrationPoint([0], 2.0)]\n",
-    "    new_field!(el, :temperature)\n",
+    "    new_fieldset!(el, \"temperature\")\n",
     "    DC2D2(el, integration_points, [])\n",
     "end\n",
     "\n",
@@ -648,51 +700,14 @@
     "\"\"\"\n",
     "function JuliaFEM.get_rhs(eq::DC2D2, ip, t)\n",
     "    el = get_element(eq)\n",
-    "    ϕ = get_basis(el)\n",
-    "    f = el[\"temperature flux\"]\n",
-    "    return ϕ(ip.xi)*f(t)\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": 20,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "2-element Array{Float64,1}:\n",
-       " 50.0\n",
-       " 50.0"
-      ]
-     },
-     "execution_count": 20,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "el = Seg2([1, 2])\n",
-    "new_field!(el, :Geometry, Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n",
-    "new_field!(el, \"temperature flux\", Field(0.0, 100.0))\n",
-    "eq = DC2D2(el)\n",
-    "integrate_rhs(eq, 1.0)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## Defining own problem\n",
-    "\n",
-    "- takes a set of elements and maps corresponding equations for them\n",
-    "- problem should have all required information in order to be solvable"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 21,
@@ -703,7 +718,9 @@
     {
      "data": {
       "text/plain": [
-       "PlaneHeatProblem"
+       "2-element Array{Float64,1}:\n",
+       " 50.0\n",
+       " 50.0"
       ]
      },
      "execution_count": 21,
@@ -711,6 +728,43 @@
      "output_type": "execute_result"
     }
    ],
+   "source": [
+    "el = Seg2([1, 2])\n",
+    "new_fieldset!(el, \"geometry\")\n",
+    "add_field!(el, \"geometry\", Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n",
+    "new_fieldset!(el, \"temperature flux\")\n",
+    "add_field!(el, \"temperature flux\", Field(0.0, 100.0))\n",
+    "eq = DC2D2(el)\n",
+    "integrate_rhs(eq, 1.0)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Defining own problem\n",
+    "\n",
+    "- main object: takes a set of elements and maps corresponding field equations to them"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "PlaneHeatProblem"
+      ]
+     },
+     "execution_count": 22,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "using JuliaFEM: Problem, get_equation, get_dimension\n",
     "\n",
@@ -722,7 +776,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 23,
    "metadata": {
     "collapsed": false
    },
@@ -733,7 +787,7 @@
        "get_equation (generic function with 3 methods)"
       ]
      },
-     "execution_count": 22,
+     "execution_count": 23,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -753,7 +807,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 24,
    "metadata": {
     "collapsed": false
    },
@@ -762,7 +816,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "07-Oct 09:28:41:DEBUG:root:total dofs: 4\n"
+      "09-Oct 01:41:22:DEBUG:root:total dofs: 4\n"
      ]
     },
     {
@@ -780,7 +834,7 @@
        "\t[2, 1]  =  300.0"
       ]
      },
-     "execution_count": 23,
+     "execution_count": 24,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -791,11 +845,15 @@
     "\n",
     "# create elements and add necessary properties like connectivity and geometry\n",
     "el1 = Quad4([2, 3, 4, 5])\n",
-    "new_field!(el1, :Geometry, Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
-    "new_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
+    "new_fieldset!(el1, \"geometry\")\n",
+    "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
+    "new_fieldset!(el1, \"temperature thermal conductivity\")\n",
+    "add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
     "el2 = Seg2([2, 3])\n",
-    "new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
-    "new_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
+    "new_fieldset!(el2, \"geometry\")\n",
+    "add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
+    "new_fieldset!(el2, \"temperature flux\")\n",
+    "add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
     "\n",
     "problem = PlaneHeatProblem()\n",
     "add_element!(problem, el1)\n",
@@ -891,7 +949,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 25,
    "metadata": {
     "collapsed": false
    },
@@ -904,7 +962,7 @@
        "\t[2, 1]  =  100.0"
       ]
      },
-     "execution_count": 24,
+     "execution_count": 25,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -926,7 +984,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 26,
    "metadata": {
     "collapsed": false
    },
@@ -937,7 +995,7 @@
        "DirichletProblem"
       ]
      },
-     "execution_count": 25,
+     "execution_count": 26,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -953,7 +1011,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 27,
    "metadata": {
     "collapsed": false
    },
@@ -964,7 +1022,7 @@
        "get_equation (generic function with 4 methods)"
       ]
      },
-     "execution_count": 26,
+     "execution_count": 27,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -972,6 +1030,8 @@
    "source": [
     "abstract DirichletBC <: Equation\n",
     "\n",
+    "get_unknown_field_name(eq::DirichletBC) = symbol(\"reaction force\")\n",
+    "\n",
     "\"\"\"\n",
     "Dirichlet boundary condition element for 2 node line segment\n",
     "(plane problems).\n",
@@ -986,7 +1046,7 @@
     "    integration_points = [\n",
     "        IntegrationPoint([-sqrt(1/3)], 1.0),\n",
     "        IntegrationPoint([+sqrt(1/3)], 1.0)]\n",
-    "    new_field!(el, \"reaction force\")\n",
+    "    new_fieldset!(el, \"reaction force\")\n",
     "    fieldval(X, t) = 0.0\n",
     "    DBC2D2(el, integration_points, [], fieldval)\n",
     "end\n",
@@ -999,10 +1059,10 @@
     "\n",
     "function JuliaFEM.get_rhs(eq::DBC2D2, ip, t)\n",
     "    el = get_element(eq)\n",
-    "    h = get_basis(el)\n",
+    "    h = get_basis(el, ip.xi)\n",
     "    f = eq.fieldval\n",
-    "    X = interpolate(el, :Geometry, ip.xi, t)\n",
-    "    return h(ip.xi)*f(X, t)\n",
+    "    X = interpolate(el, \"geometry\", ip.xi, t)\n",
+    "    return h*f(X, t)\n",
     "end\n",
     "JuliaFEM.has_lhs(eq::DBC2D2) = true\n",
     "JuliaFEM.has_rhs(eq::DBC2D2) = true\n",
@@ -1010,42 +1070,43 @@
     "JuliaFEM.get_equation(pr::Type{DirichletProblem}, el::Type{Seg2}) = DBC2D2"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 27,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "1-element Array{JuliaFEM.Equation,1}:\n",
-       " DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(:Geometry=>JuliaFEM.Field[JuliaFEM.Field{Array{Array{T,1},1}}(0.0,1,Array{T,1}[[0.0,0.0],[0.0,1.0]])],symbol(\"reaction force\")=>JuliaFEM.Field[])),[JuliaFEM.IntegrationPoint([-0.5773502691896257],1.0,Dict{Any,Any}()),JuliaFEM.IntegrationPoint([0.5773502691896257],1.0,Dict{Any,Any}())],Int64[],fieldval)"
-      ]
-     },
-     "execution_count": 27,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "# create elements and add necessary properties like connectivity and geometry\n",
-    "el3 = Seg2([4, 5])\n",
-    "new_field!(el3, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
-    "\n",
-    "bc1 = DirichletProblem()\n",
-    "#get_equation(typeof(bc1), typeof(el3))\n",
-    "#DBC2D2(el3)\n",
-    "add_element!(bc1, el3)"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 28,
    "metadata": {
     "collapsed": false
    },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1-element Array{JuliaFEM.Equation,1}:\n",
+       " DBC2D2(JuliaFEM.Seg2([4,5],JuliaFEM.Basis(basis,j),Dict(symbol(\"reaction force\")=>JuliaFEM.Field[],: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)"
+      ]
+     },
+     "execution_count": 28,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# create elements and add necessary properties like connectivity and geometry\n",
+    "el3 = Seg2([4, 5])\n",
+    "new_fieldset!(el3, \"geometry\")\n",
+    "add_field!(el3, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
+    "\n",
+    "bc1 = DirichletProblem()\n",
+    "#get_equation(typeof(bc1), typeof(el3))\n",
+    "#DBC2D2(el3)\n",
+    "add_element!(bc1, el3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "data": {
@@ -1055,7 +1116,7 @@
        "\t[4, 1]  =  0.0"
       ]
      },
-     "execution_count": 28,
+     "execution_count": 29,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1080,7 +1141,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 30,
    "metadata": {
     "collapsed": false
    },
@@ -1095,7 +1156,7 @@
        " -1.0  -2.0  -1.0   4.0"
       ]
      },
-     "execution_count": 29,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1106,7 +1167,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 31,
    "metadata": {
     "collapsed": false
    },
@@ -1121,7 +1182,7 @@
        "   0.0"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 31,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1132,7 +1193,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 32,
    "metadata": {
     "collapsed": false
    },
@@ -1147,7 +1208,7 @@
        " 0.0  0.0  0.166667  0.333333"
       ]
      },
-     "execution_count": 31,
+     "execution_count": 32,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1158,7 +1219,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 33,
    "metadata": {
     "collapsed": false
    },
@@ -1173,7 +1234,7 @@
        " 0.0"
       ]
      },
-     "execution_count": 32,
+     "execution_count": 33,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1184,7 +1245,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 34,
    "metadata": {
     "collapsed": false
    },
@@ -1217,7 +1278,7 @@
        "\t[4, 8]  =  0.333333"
       ]
      },
-     "execution_count": 33,
+     "execution_count": 34,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1228,7 +1289,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 35,
    "metadata": {
     "collapsed": false
    },
@@ -1243,7 +1304,7 @@
        "\t[8, 1]  =  0.0"
       ]
      },
-     "execution_count": 34,
+     "execution_count": 35,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1300,7 +1361,7 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Unique rows: [1,2,3,4,7,8]\n"
+      "Non-zero rows: [1,2,3,4,7,8]\n"
      ]
     },
     {
@@ -1324,7 +1385,7 @@
    ],
    "source": [
     "r = unique(rowvals(Atot))\n",
-    "println(\"Unique rows: $r\")\n",
+    "println(\"Non-zero rows: $r\")\n",
     "xtot = zeros(btot)\n",
     "F = lufact(Atot[r,r])\n",
     "s = full(btot[r])\n",
@@ -1332,6 +1393,282 @@
     "full(xtot)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "## Developing own solver\n",
+    "\n",
+    "Last part. Defining own solver.\n",
+    "\n",
+    "- takes a set of problems (typically main problem + boundary problems)\n",
+    "- solves them, updates fields"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "add_problem! (generic function with 1 method)"
+      ]
+     },
+     "execution_count": 38,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "abstract Solver\n",
+    "\n",
+    "\"\"\"\n",
+    "Simple solver for educational purposes.\n",
+    "\"\"\"\n",
+    "type SimpleSolver <: Solver\n",
+    "    problems\n",
+    "end\n",
+    "\n",
+    "\"\"\"\n",
+    "Default initializer\n",
+    "\"\"\"\n",
+    "function SimpleSolver()\n",
+    "    SimpleSolver(Problem[])\n",
+    "end\n",
+    "\n",
+    "function add_problem!(solver::Solver, problem::Problem)\n",
+    "    push!(solver.problems, problem)\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "call (generic function with 1253 methods)"
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "get_problems(s::Solver) = s.problems\n",
+    "\"\"\"\n",
+    "Call solver to solve a set of problems.\n",
+    "\n",
+    "This is simple serial solver for demonstration purposes. It handles the most\n",
+    "common situation, i.e., some main field problem and it's Dirichlet boundary.\n",
+    "\"\"\"\n",
+    "function call(solver::SimpleSolver, t)\n",
+    "    problems = get_problems(solver)\n",
+    "    problem1 = problems[1]\n",
+    "    problem2 = problems[2]\n",
+    "\n",
+    "    # calculate order of degrees of freedom in global matrix\n",
+    "    # and set the ordering to problems\n",
+    "    dofmap = calculate_global_dofs(problem1)\n",
+    "    assign_global_dofs!(problem1, dofmap)\n",
+    "    assign_global_dofs!(problem2, dofmap)\n",
+    "\n",
+    "    # assemble problem 1\n",
+    "    A1 = sparse(get_lhs(problem1, t)...)\n",
+    "    b1 = sparsevec(get_rhs(problem1, t)..., size(A1, 1))\n",
+    "\n",
+    "    # assemble problem 2\n",
+    "    A2 = sparse(get_lhs(problem2, t)...)\n",
+    "    b2 = sparsevec(get_rhs(problem2, t)..., size(A2, 1))\n",
+    "    \n",
+    "    # make one monolithic assembly\n",
+    "    A = [A1 A2; A2' zeros(A2)]\n",
+    "    b = [b1; b2]\n",
+    "\n",
+    "    # solve problem\n",
+    "    nz = unique(rowvals(A))\n",
+    "    x = zeros(b)\n",
+    "    x[nz] = lufact(Atot[nz,nz]) \\ full(b[nz])\n",
+    "\n",
+    "    # get \"problem-wise\" solution vectors\n",
+    "    x1 = x[1:length(b1)]\n",
+    "    x2 = x[length(b1)+1:end]\n",
+    "\n",
+    "    # check residual\n",
+    "    R1 = A1*x1 - b1\n",
+    "    R2 = A2*x2 - b2\n",
+    "    println(\"Residual norm: $(norm(R1+R2))\")\n",
+    "\n",
+    "    # update field for elements in problem 1\n",
+    "    for equation in get_equations(problem1)\n",
+    "        gdofs = get_global_dofs(equation)\n",
+    "        element = get_element(equation)\n",
+    "        field_name = get_unknown_field_name(equation)  # field we are solving\n",
+    "        field = Field(t, full(x1[gdofs])[:])\n",
+    "        add_field!(element, field_name, field)\n",
+    "    end\n",
+    "\n",
+    "    # update field for elements in problem 2\n",
+    "    for equation in get_equations(problem2)\n",
+    "        gdofs = get_global_dofs(equation)\n",
+    "        element = get_element(equation)\n",
+    "        field_name = get_unknown_field_name(equation)\n",
+    "        field = Field(t, full(x2[gdofs]))\n",
+    "        add_field!(element, field_name, field)\n",
+    "    end\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Summary\n",
+    "\n",
+    "Next we collect all things together and solve Poisson equation in domain $\\Omega = [0,1]\\times[0,1]$ with some Neumann boundary on $\\Gamma_1$ and Dirichlet boundary on $\\Gamma_2$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5-element Array{Array{Float64,1},1}:\n",
+       " [-1.0,0.0]\n",
+       " [-0.5,0.5]\n",
+       " [0.0,1.0] \n",
+       " [0.5,1.5] \n",
+       " [1.0,2.0] "
+      ]
+     },
+     "execution_count": 40,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "\"\"\"\n",
+    "Simple linspace extension to multidimensional values.\n",
+    "\"\"\"\n",
+    "function Base.linspace(X1, X2, n)\n",
+    "    [1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]\n",
+    "end\n",
+    "\n",
+    "linspace([-1.0, 0.0], [1.0, 2.0], 5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "09-Oct 01:41:25:DEBUG:root:total dofs: 4\n"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "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": [
+    "using JuliaFEM: get_fieldset\n",
+    "\n",
+    "# Define Problem 1:\n",
+    "# - Field function: Laplace equation Δu=0 in Ω={u∈R²|(x,y)∈[0,1]×[0,1]}\n",
+    "# - Neumann boundary on Γ₁={0<=x<=1, y=0}, ∂u/∂n=600 on Γ₁\n",
+    "\n",
+    "el1 = Quad4([1, 2, 3, 4])\n",
+    "new_fieldset!(el1, \"geometry\")\n",
+    "add_field!(el1, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
+    "new_fieldset!(el1, \"temperature thermal conductivity\")\n",
+    "add_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
+    "\n",
+    "el2 = Seg2([1, 2])\n",
+    "new_fieldset!(el2, \"geometry\")\n",
+    "add_field!(el2, \"geometry\", Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0]]))\n",
+    "\n",
+    "# Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n",
+    "new_fieldset!(el2, \"temperature flux\")\n",
+    "add_field!(el2, \"temperature flux\", Field(0.0, 0.0))\n",
+    "add_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
+    "\n",
+    "problem1 = PlaneHeatProblem()\n",
+    "add_element!(problem1, el1)\n",
+    "add_element!(problem1, el2)\n",
+    "\n",
+    "# Define Problem 2:\n",
+    "# - Dirichlet boundary Γ₂={0<=x<=1, y=1}, u=0 on Γ₂\n",
+    "\n",
+    "el3 = Seg2([3, 4])\n",
+    "new_fieldset!(el3, \"geometry\")\n",
+    "add_field!(el3, \"geometry\", Field(0.0, Vector[[1.0, 1.0], [0.0, 1.0]]))\n",
+    "\n",
+    "problem2 = DirichletProblem()\n",
+    "add_element!(problem2, el3)\n",
+    "\n",
+    "# Create a solver for a set of problems\n",
+    "solver = SimpleSolver()\n",
+    "add_problem!(solver, problem1)\n",
+    "add_problem!(solver, problem2)\n",
+    "\n",
+    "# Solve problem at time t=1.0 and update fields\n",
+    "call(solver, 1.0)\n",
+    "\n",
+    "# Postprocess.\n",
+    "# Interpolate temperature field along boundary of Γ₁ at time t=1.0\n",
+    "xi = linspace([-1.0], [1.0], 5)\n",
+    "X = interpolate(el2, \"geometry\", xi, 1.0)\n",
+    "T = interpolate(el2, \"temperature\", xi, 1.0)\n",
+    "println(X)\n",
+    "println(T)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## TODO\n",
+    "\n",
+    "- dynamics\n",
+    "- different basis for trial and test (Petrov-Galerkin)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Frequently asked questions\n",
+    "\n",
+    "Any question about data structures, better ideas and improvements are very welcome."
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 6f57e0e..de5b55b 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -7,7 +7,9 @@ using Lexicon
 using Logging
 @Logging.configure(level=DEBUG)
 
-include("types.jl") # type definitions
+include("types.jl")  # type definitions
+
+include("interpolate.jl")  # interpolation routines
 
 ### ELEMENTS ###
 include("elements.jl")
diff --git a/src/elements.jl b/src/elements.jl
index 9c54896..d8a75b3 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -8,9 +8,11 @@ Related notebooks
 2015-08-29-developing-juliafem.ipynb
 =#
 
+using JuliaFEM: interpolate
 using FactCheck
 using ForwardDiff
 
+
 abstract Element
 
 #= ELEMENT DEFINITIONS
@@ -143,9 +145,10 @@ function test_element(eltype)
     fld = Field(0.0, collect(1:n))
     Logging.info("Creating new scalar field $fld")
     Logging.info("Pushing field to element.")
-    new_field!(el, :field1)
-    push_field!(el, :field1, fld)
-    @fact el[:field1][1] --> fld
+    new_fieldset!(el, "field1")
+    add_field!(el, "field1", fld)
+    fieldset = get_fieldset(el, "field1")
+    @fact fieldset[1] --> fld
 
     mid = zeros(dim)
     try
@@ -164,8 +167,9 @@ function test_element(eltype)
     end
 
     Logging.info("Interpolating scalar field at $mid")
-    f(field, xi, t) = el(xi)*el[field](t)
-    i = f(:field1, mid, 0.0)
+    #f(field, xi, t) = el(xi)*el[field](t)
+    #i = f(:field1, mid, 0.0)
+    i = interpolate(el, "field1", mid, 0.0)
     Logging.info("Value: $i")
     Logging.info("Element $eltype passed tests.")
 end
@@ -188,21 +192,22 @@ get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi)
 """
 Interpolate field on element.
 """
-function interpolate(el::Element, field::Symbol, xi::Vector, t::Number)
-    get_basis(el, xi)*el[field](t)
-end
-function interpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
-    interpolate(el, Symbol(field), xi, t)
+function interpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number)
+    fieldset = get_fieldset(el, symbol(field_name))
+    field = interpolate(fieldset, t)
+    basis = get_basis(el)
+    interpolate(basis, field, xi)
 end
 
 """
 Interpolate derivative of field on element.
 """
-function dinterpolate(el::Element, field::Symbol, xi::Vector, t::Number)
-    get_dbasisdxi(el, xi)*el[field](t)
-end
-function dinterpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
-    dinterpolate(el, Symbol(field), xi, t)
+function dinterpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number)
+    #get_dbasisdxi(el, xi)*el[field](t)
+    fieldset = get_fieldset(el, symbol(field_name))
+    field = interpolate(fieldset, t)
+    basis = get_basis(el)
+    dinterpolate(basis, field, xi)
 end
 
 """
@@ -210,155 +215,64 @@ Get jacobian of element evaluated at point ξ on element in reference configurat
 
 Parameters
 ----------
-el::Element
-xi::Vector
-geometry_field::Any, optional
-time::Number, optional, default=0.0
+el :: Element
+xi :: Vector
+geometry_field :: Any, optional
+time :: Number
 
 Returns
 -------
 Vector or Matrix
     depending on element type
 
-Notes
------
-Big "J" comes from reference (undeformed) configuration.
 """
-function get_Jacobian(el::Element, xi, t, geometry_field=:Geometry)
+function get_jacobian(el::Element, xi, t, geometry_field=symbol("geometry"))
     dinterpolate(el, geometry_field, xi, t)
 end
 
-
-"""
-Get jacobian of element evaluated at point ξ on element in current configuration.
-
-Notes
------
-Small "j" comes from current (deformed) configuration.
-"""
-function get_jacobian(el::Element, xi, t, geometry_field=:Geometry, displacement_field=:displacement)
-    dbasisdxi = get_dbasisdxi(el, xi)
-    X = get_field(el, geometry_field)(t)
-    u = get_field(el, displacement_field)(t)
-    j = dbasisdxi*(X+u)
-    return j
-end
-
-
 """
 Evaluate partial derivatives of basis, dbasis/dX
 """
 function get_dbasisdX(el::Element, xi, t)
     dbasisdxi = get_dbasisdxi(el, xi)
-    J = get_Jacobian(el, xi, t)
+    J = get_jacobian(el, xi, t)
     dbasisdxi*inv(J)
 end
 
 
-"""
-Evaluate partial derivatives of basis, dbasis/dx
-"""
-function get_dbasisdx(el::Element, xi, t)
-    dbasisdxi = get_dbasisdxi(el, xi)
-    j = get_jacobian(el, xi, t)
-    dbasisdxi*inv(j)
+""" Create new empty set of fields for element. """
+function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString})
+    el.fields[symbol(field_name)] = FieldSet()
+end
+function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field)
+    new_fieldset!(el, symbol(field_name))
+    add_field!(el, symbol(field_name), field)
 end
 
 
-""" Create new empty field of some type. """
-function new_field!(el::Element, field_name::Symbol)
-    el.fields[field_name] = Field[]
-end
-function new_field!(el::Element, field_name::Symbol, field::Field)
-    new_field!(el, field_name)
-    push_field!(el, field_name, field)
-end
-function new_field!(el::Element, field_name::ASCIIString, field::Field)
-    new_field!(el, Symbol(field_name), field)
-end
-function new_field!(el::Element, field_name::ASCIIString)
-    new_field!(el, Symbol(field_name))
+""" Add new field to fieldset of element. """
+function add_field!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field)
+    push!(el.fields[symbol(field_name)], field)
 end
 
 
-""" Push to existing set field of fields. """
-function push_field!(el::Element, field_name::Symbol, field::Field)
-    push!(el.fields[field_name], field)
+""" Get fieldset. """
+function get_fieldset(el::Element, field_name::Union{Symbol, ASCIIString})
+    el.fields[symbol(field_name)]
 end
-function push_field!(el::Element, field_name::ASCIIString, field::Field)
-    push_field!(el, Symbol(field_name), field)
+""" Get fieldset, convenient function. """
+function Base.getindex(el::Element, field_name::Union{Symbol, ASCIIString})
+    get_fieldset(el, field_name)
 end
 
 
-""" Get field variable. """
-function get_field(el::Element, field_name::Symbol)
-    el.fields[field_name]
-end
-function get_field(el::Element, field_name::ASCIIString)
-    el.fields[Symbol(field_name)]
-end
-function Base.getindex(el::Element, field_name::Union{ASCIIString, Symbol})
-    get_field(el, field_name)
-end
 
 
-#=
-"""
-Evaluate some field in point ξ on element using basis functions.
-
-Parameters
-----------
-el :: Element
-field :: Any
-xi :: Vector
-
-Returns
--------
-Scalar, Vector, Tensor, depending on what is type of field to interpolate.
-
-Notes
------
-This has another version which returns multiple values for set of coordinates {ξᵢ}.
-dinterpolate returns derivatives.
-
-Examples
---------
->>> field = [1.0, 2.0, 3.0, 4.0]
->>> set_field(el, :temperature, field)
->>> interpolate(el, :temperature, [0.0, 0.0])
-15.0
-"""
-function interpolate(el::Element, field, xi::Number)
-    interpolate(el, field, [xi])
-end
-function interpolate(el::Element, field, xi::Vector)
-    field = get_field(el, field)
-    sum(get_basis(el, xi) .* field)
-end
-function interpolate(el::Element, field, xis::Array{Vector, 1})
-    field = get_field(el, field)
-    interpolate_(xi) = sum(get_basis(el, xi) .* field)
-    map(interpolate_, xis)
-end
-
-function dinterpolate(el::Element, field, xi::Number)
-    dinterpolate(el, field, [xi])
-end
-function dinterpolate(el::Element, field, xi::Vector)
-    fld = get_field(el, field)
-    dbasis = get_dbasisdxi(el, xi)
-    if isa(dbasis, Vector)
-        return sum(dbasis .* fld)
-    end
-    return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
-end
-=#
-
 """
 calculate "local" normals in elements, in a way that
 n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1]
 """
-function calculate_normals!(el::Element, t, field_name=:Normals)
+function calculate_normals!(el::Element, t, field_name=symbol("normals"))
     new_field!(el, field_name, Vector)
     for xi in Vector[[-1.0], [1.0]]
         t = dinterpolate(el, :Geometry, xi)
@@ -371,7 +285,7 @@ end
 """
 Alter normal field such that normals of adjacent elements are averaged.
 """
-function average_normals!(elements, normal_field=:Normals)
+function average_normals!(elements, normal_field=symbol("normals"))
     d = Dict()
     for el in elements
         c = get_connectivity(el)
diff --git a/src/equations.jl b/src/equations.jl
index 89a40b0..7acd838 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -35,6 +35,8 @@ get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq),
 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
 
 
 """
@@ -48,7 +50,7 @@ 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)
+    J = get_jacobian(el, xi, t)
     s = size(J)
     return s[1] == s[2] ? det(J) : norm(J)
 end
@@ -79,5 +81,3 @@ function set_global_dofs!(eq::Equation, dofs)
     eq.global_dofs = dofs
 end
 
-# Equations for heat problems
-#include("heat_equations.jl")
diff --git a/src/interpolate.jl b/src/interpolate.jl
new file mode 100644
index 0000000..f4c66ee
--- /dev/null
+++ b/src/interpolate.jl
@@ -0,0 +1,57 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+using JuliaFEM: Basis, Field, FieldSet, diff
+
+
+"""
+Interpolate field u using basis N in point xi.
+"""
+function interpolate{T}(N::Basis, u::Field{Vector{T}}, xi::Array{Float64,1})
+    N(xi)*u
+end
+"""
+Interpolate field u using basis N in set of points xi. Convenient function.
+"""
+function interpolate{T}(N::Basis, u::Field{Vector{T}}, xis::Array{Array{Float64,1},1})
+    T[N(xi)*u for xi in xis]
+end
+function interpolate{T}(N::Basis, u::Field{T}, xi::Array{Float64,1})
+    u.values
+end
+
+"""
+Interpolate a field from fieldset for some time t.
+"""
+function interpolate(fields::FieldSet, t::Number)
+    if length(fields) == 0
+        throw("Empty set of fields.")
+    end
+    if t <= fields[1].time
+        return Field(t, fields[1].values)
+    end
+    if t >= fields[end].time
+        return fields[end]
+    end
+    i = length(fields)
+    while fields[i].time >= t
+        i -= 1
+    end
+    if fields[i].time == t
+        return fields[i]
+    end
+    #Logging.debug("doing linear interpolation between fields $i and $(i+1)")
+    f1 = fields[i]
+    t1 = f1.time
+    f2 = fields[i+1]
+    t2 = f2.time
+    dt = t2 - t1
+    nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
+    f = Field(t, nw)
+    return f
+end
+
+function dinterpolate(N::Basis, u::Field, xi::Array{Float64, 1})
+    dN = diff(N)
+    dN(xi)*u
+end
diff --git a/src/problems.jl b/src/problems.jl
index 914df2a..bae79d6 100644
--- a/src/problems.jl
+++ b/src/problems.jl
@@ -10,9 +10,9 @@ get_equation(pr::Type{Problem}, el::Type{Element}) = nothing
 """
 Add new element to problem
 """
-function add_element!(pr::Problem, el::Element)
-    eq = get_equation(typeof(pr), typeof(el))
-    push!(pr.equations, eq(el))
+function add_element!(problem::Problem, element::Element)
+    equation = get_equation(typeof(problem), typeof(element))
+    push!(problem.equations, equation(element))
 end
 
 """
diff --git a/src/types.jl b/src/types.jl
index 6d03c90..dcd1f0d 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -5,89 +5,55 @@
 
 using ForwardDiff
 
-
-""" Field. """
+""" Field is a fundamental type which holds some values in some time t """
 type Field{T}
     time :: Float64
     increment :: Int64
     values :: T
 end
-
 """ Initialize field. """
 function Field(time, values)
     Field(time, 1, values)
 end
-
 """ Get length of a field (number of basis functions in practice). """
-Base.length(f::Field) = length(f.values)
-
+function Base.length(f::Field)
+    length(f.values)
+end
 """ Get field discrete value at point i. """
-Base.getindex(f::Field, i::Int64) = f.values[i]
-
-""" Interpolate field h(ξ)*f = x*f """
-function interpolate{T}(x::Vector, f::Field{Vector{T}})
+function Base.getindex(f::Field, i::Int64)
+    f.values[i]
+end
+""" Multiply field with some constant k. """
+function Base.(:*)(k::Number, f::Field)
+    Field(f.time, k*f.values)
+end
+""" Multiply field with some vector x. """
+function Base.(:*)(x::Vector, f::Field)
+    @assert length(x) == length(f)
     sum([f[i]*x[i] for i in 1:length(f)])
 end
-function interpolate{T}(x::Matrix, f::Field{Vector{T}})
+""" Multiply field with some matrix x. """
+# function Base.(:*){T}(x::Matrix, f::Field{Vector{T}})
+function Base.(:*)(x::Matrix, f::Field)
     sum([f[i]*x[i,:] for i in 1:length(f)])
 end
-function interpolate(x::Vector, f::Field)
-    f.values*x
-end
-Base.(:*)(x::Union{Vector, Matrix}, f::Field) = interpolate(x, f)
-
-""" Interpolate field (h*f)(ξ) """
-Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
-
-""" Multiply field with some constant k. """
-Base.(:*)(k::Float64, f::Field) = Field(f.time, k*f.values)
-
 """ Sum two fields. """
 function Base.(:+)(f1::Field, f2::Field)
     @assert(f1.time == f2.time, "Cannot add fields: time mismatch, $(f1.time) != $(f2.time)")
     Field(f1.time, f1.values + f2.values)
 end
 
-""" Interpolate from set of fields x*[f1, f2] where x is evaluated basis. """
+"""
+FieldSet is array of fields, each field maybe having different time and/or increment.
+"""
+typealias FieldSet Array{Field, 1}
+""" Multiply fieldset with some vector x. """
 Base.(:*)(x::Array{Float64, 1}, f::Array{Field}) = sum(x .* f)
 
-""" Interpolate from set of fields with basis b, i.e. f(t) = b(t)*[f1, f2] """
-Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
-
-"""
-Interpolate a field from finite set of fields some time t ∈ R.
-"""
-function call(fields :: Array{Field, 1}, t::Number)
-    if length(fields) == 0
-        throw("Empty set of fields.")
-    end
-    if t <= fields[1].time
-        return Field(t, fields[1].values)
-    end
-    if t >= fields[end].time
-        return fields[end]
-    end
-    i = length(fields)
-    while fields[i].time >= t
-        i -= 1
-    end
-    if fields[i].time == t
-        return fields[i]
-    end
-    #Logging.debug("doing linear interpolation between fields $i and $(i+1)")
-    f1 = fields[i]
-    t1 = f1.time
-    f2 = fields[i+1]
-    t2 = f2.time
-    dt = t2 - t1
-    nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
-    f = Field(t, nw)
-    return f
+""" Add new field to fieldset. """
+function add_field!(fs::FieldSet, field::Field)
+    push!(fs, field)
 end
-function call(field::Field, t::Float64)
-    Field(t, field.increment, field.values)
-end
-
 
 
 """ Basis function. """
@@ -95,21 +61,23 @@ type Basis
     basis :: Function
     dbasisdxi :: Function
 end
-
 """ Constructor of basis function. """
 function Basis(basis)
     Basis(basis, ForwardDiff.jacobian(basis))
 end
-
-""" Interpolate field f using basis b. """
-Base.(:*)(b::Basis, f::Field) = (x) -> b(x)*f
-Base.(:*)(b::Basis, f::Array{Field}) = (t) -> b(t)*f
-
-""" Evaluate basis function in point ξ. """
-call(b::Basis, xi) = b.basis(xi)
-
 """ Get partial derivative of basis function. """
-∂(h::Basis) = h.dbasisdxi
 diff(h::Basis) = h.dbasisdxi
 derivative(h::Basis) = h.dbasisdxi
 
+
+# convenient functions
+""" Evaluate basis function in point ξ. """
+call(b::Basis, xi) = b.basis(xi)
+#""" 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] """
+#Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
+#""" Interpolate field f using basis b. """
+#Base.(:*)(b::Basis, f::Field) = (x) -> b(x)*f
+#Base.(:*)(b::Basis, f::Array{Field}) = (t) -> b(t)*f
+
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 7353438..ec162d1 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -2,7 +2,49 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 using FactCheck
-using JuliaFEM: test_element
+using JuliaFEM: Element, Basis, FieldSet
+
+# prototype element
+type MockElement <: Element
+    connectivity :: Array{Int, 1}
+    basis :: Basis
+    fields :: Dict{Symbol, FieldSet}
+end
+function MockElement(connectivity)
+    h(xi) = [
+        (1-xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1+xi[2])/4
+        (1-xi[1])*(1+xi[2])/4]
+    dh(xi) = [
+        -(1-xi[2])/4.0   -(1-xi[1])/4.0
+         (1-xi[2])/4.0   -(1+xi[1])/4.0
+         (1+xi[2])/4.0    (1+xi[1])/4.0
+        -(1+xi[2])/4.0    (1-xi[1])/4.0]
+    basis = Basis(h, dh)
+    MockElement(connectivity, basis, Dict())
+end
+JuliaFEM.get_number_of_basis_functions(el::Type{MockElement}) = 4
+JuliaFEM.get_element_dimension(el::Type{MockElement}) = 2
+
+
+using JuliaFEM: test_element
+facts("test test_element against mock element") do
+    test_element(MockElement)
+end
+
+
+using JuliaFEM: new_fieldset!, add_field!, Field, get_fieldset
+facts("test adding fieldsets and fields to element") do
+    el = MockElement([1, 2, 3, 4])
+    fieldset = new_fieldset!(el, "geometry")
+    field1 = Field(0.0, [0.0, 0.0, 0.0, 0.0])
+    add_field!(el, "geometry", field1)
+    field2 = Field(1.0, [1.0, 1.0, 1.0, 1.0])
+    add_field!(fieldset, field2)
+    fields = get_fieldset(el, "geometry")
+    @fact length(fields) --> 2
+    @fact fields[1] --> field1
+    @fact fields[2] --> field2
+end
 
-using JuliaFEM: Quad4
-test_element(Quad4)
diff --git a/test/test_types.jl b/test/test_types.jl
index 90c90ee..936ac84 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -1,18 +1,11 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-using JuliaFEM: Basis, Field, get_field, diff
+using JuliaFEM: Basis, Field, FieldSet, interpolate, dinterpolate
 using FactCheck
 
-facts("test fields and interpolation") do
-
-    # simple interpolation in domain [-1, 1]
-    N = Basis((ξ) -> [0.5*(1.0-ξ[1]), 0.5*(1.0+ξ[1])])
-    u = Field(0.0, [0.0, 1.0])
-    @fact N([0.0])*u --> 0.5
-    @fact (N*u)([0.0]) --> 0.5
-
-    # multiply of field with constant
+facts("test fields") do
+    # multiple field with some constant
     u1 = Field(0.0, [0.0, 1.0])
     u2 = 3.0*u1
     @fact u1.time --> 0.0
@@ -24,48 +17,68 @@ facts("test fields and interpolation") do
     u2 = Field(0.0, [1.0, 2.0])
     u3 = u1 + u2
     @fact u3.values --> [1.0, 3.0]
+end
 
-    # interpolation between two fields in time domain
+facts("test interpolation of fields") do
+
+    # interpolation of field in spatial domain
+    N = Basis((xi) -> [0.5*(1.0-xi[1]), 0.5*(1.0+xi[1])])
+    u = Field(0.0, [0.0, 1.0])
+    @fact interpolate(N, u, [0.0]) --> 0.5
+
+    # interpolation of fieldset in time domain
     u1 = Field(0.0, [0.0, 1.0])
-    u2 = Field(0.0, [1.0, 2.0])
-    t = Basis((t) -> [1-t, t])
-    u = Field[u1, u2]
-    u2 = (t*u)(0.5)
-    @fact u2.values --> [0.5, 1.5]
+    u2 = Field(1.0, [1.0, 2.0])
+    u = FieldSet([u1, u2])
+    @fact interpolate(u, 0.5).values --> [0.5, 1.5]
+    @fact interpolate(u, 0.5).time --> 0.5
 
-    # interpolation in set of fields is defined for every time value
+    # interpolation of fieldset is defined for every time value:
     u1 = Field(0.0, [0.0, 0.0])
     u2 = Field(1.0, [1.0, 2.0])
     u3 = Field(2.0, [0.5, 1.5])
-    u = Field[u1, u2, u3]
-    @fact u(-1.0).values --> [0.0, 0.0]  # "out of range -" -> first known value
-    @fact u(0.0).values --> [0.0, 0.0]
-    @fact u(1.0).values --> [1.0, 2.0]
-    @fact u(2.0).values --> [0.5, 1.5]
-    @fact u(3.0).values --> [0.5, 1.5]  # "out of range +" -> last known value
-    @fact u(0.5).values --> [0.5, 1.0]
-    @fact u(1.5).values --> [0.75, 1.75]
+    u = FieldSet([u1, u2, u3])
+    @fact interpolate(u, -1.0).values --> [0.0, 0.0]  # "out of range -" -> first known value
+    @fact interpolate(u, 0.0).values --> [0.0, 0.0]
+    @fact interpolate(u, 1.0).values --> [1.0, 2.0]
+    @fact interpolate(u, 2.0).values --> [0.5, 1.5]
+    @fact interpolate(u, 3.0).values --> [0.5, 1.5]  # "out of range +" -> last known value
+    @fact interpolate(u, 0.5).values --> [0.5, 1.0]
+    @fact interpolate(u, 1.5).values --> [0.75, 1.75]
     # use Inf to get very first or last value of field
-    @fact u(-Inf).values --> [0.0, 0.0]
-    @fact u(+Inf).values --> [0.75, 1.75]
+    @fact interpolate(u, -Inf).values --> [0.0, 0.0]
+    @fact interpolate(u, +Inf).values --> [0.5, 1.5]
 
     # multidimensional interpolation with and without derivatives
-    X = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
-    h = Basis((xi) ->
-    [(1-xi[1])*(1-xi[2])/4
-     (1+xi[1])*(1-xi[2])/4
-     (1+xi[1])*(1+xi[2])/4
-     (1-xi[1])*(1+xi[2])/4])
-   # midpoint of field
-   @fact (h*X)([0.0, 0.0]) --> [0.5, 0.5]
-   @fact h([0.0, 0.0])*X --> [0.5, 0.5]
-   # derivatives of field at midpoint
-   @fact diff(h)([0.0, 0.0])*X --> [0.5 0.0; 0.0 0.5]
-   @fact (diff(h)*X)([0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
+    h(xi) = [
+        (1-xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1-xi[2])/4
+        (1+xi[1])*(1+xi[2])/4
+        (1-xi[1])*(1+xi[2])/4]
+    dh(xi) = [
+        -(1-xi[2])/4.0   -(1-xi[1])/4.0
+         (1-xi[2])/4.0   -(1+xi[1])/4.0
+         (1+xi[2])/4.0    (1+xi[1])/4.0
+        -(1+xi[2])/4.0    (1-xi[1])/4.0]
+    N = Basis(h, dh)
 
-   # multiplying scalar field with a vector -> vector
-   b = Basis((xi) -> [1/2*(1-xi[1]), 1/2*(1+xi[1])])
-   f = Field(0.0, 100.0)
-   @fact b(0.0) * f --> [50.0, 50.0]
+    X = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+
+    # get midpoint of field in spatial domain
+    @fact interpolate(N, X, [0.0, 0.0]) --> [0.5, 0.5]
+    # derivatives of field at midpoint
+    @fact dinterpolate(N, X, [0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
+
+    # interpolate of scalar field -> scalar
+    H = Field(0.0, 6.0)
+    @fact interpolate(N, H, [0.0, 0.0]) --> 6.0
+
+    # multiplying scalar field with a vector -> vector
+    # this is actually not so good idea...
+    #h(xi) = [1/2*(1-xi[1]), 1/2*(1+xi[1])]
+    #dh(xi) = [-1/2 1/2]'
+    #N = Basis(h, dh)
+    #f = Field(0.0, 100.0)
+    #@fact interpolate(N, f, [0.0]) --> [50.0, 50.0]
 
 end