diff --git a/docs/tutorials/2015-08-29-developing-juliafem.ipynb b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
index a485323..4de1b5b 100644
--- a/docs/tutorials/2015-08-29-developing-juliafem.ipynb
+++ b/docs/tutorials/2015-08-29-developing-juliafem.ipynb
@@ -8,23 +8,23 @@
     "\n",
     "Author(s): Jukka Aho\n",
     "\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",
+    "**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 aiming 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",
-    "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",
+    "For demonstrational purposes the Poisson equation is used in examples to demonstrate the main concepts of the package. Poisson equation has some analogies to solid mechanics. Consider for example $(EAu')' + q = 0$, what 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",
+    "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 the 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",
+    "In this notebook the following types are briefly described:\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, 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",
+    "**Equation**: object holding some field equation, like Poisson equation $\\Delta u = f$, elasticity equation $\\nabla \\cdot \\sigma = f$, etc., \"solvable\". This is the physics we are solving by approximating the equations in the element area by basis functionsprovided 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",
+    "**Problem**: object which maps equations to elements. For example, HeatProblem maps Poisson equation to Lagrange elements and ElasticityProblem 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. The is the section where any numerical crunching happens and is targeted to high-performance computing. \n",
+    "**Solver**: object 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",
-    "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",
+    "These concepts are quite well separated so that development process can focus to point 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",
@@ -32,7 +32,7 @@
     "- [ ] **`test_problem`** \n",
     "- [ ] **`test_solver`**\n",
     "\n",
-    "which can be used to test that implementation has all necessary things defined.\n",
+    "what 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",
@@ -306,18 +306,12 @@
       "INFO: element dimension: 2 x 4\n",
       "INFO: Initializing element\n",
       "INFO: basis at [0.0,0.0]: [0.25 0.25 0.25 0.25]\n",
-      "INFO: field val at [0.0,0.0]: 0.0\n"
-     ]
-    },
-    {
-     "ename": "LoadError",
-     "evalue": "LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
-     "output_type": "error",
-     "traceback": [
-      "LoadError: MethodError: `inv` has no method matching inv(::Array{Float64,1})\nwhile loading In[9], in expression starting on line 2",
-      "",
-      " in call at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:147",
-      " in test_element at /home/jukka/.julia/v0.4/JuliaFEM/src/elements.jl:57"
+      "INFO: field val at [0.0,0.0]: 2.5\n",
+      "INFO: derivative of basis at [0.0,0.0]:\n",
+      "[-0.5 0.5 0.5 -0.5\n",
+      " -0.5 -0.5 0.5 0.5]\n",
+      "INFO: field val at [0.0,0.0]: [0.0 2.0]\n",
+      "INFO: Element MyQuad4 passed tests.\n"
      ]
     }
    ],
@@ -330,12 +324,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "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:"
+    "If `test_element` passes, elements *interface* should be well defined. After building element, one can interpolate things in it. Here we create three new fields `temperature`, `geometry` and `heat coefficient`, add some values for them in time $t=0.0$ and $t=1.0$ and interpolate. The concept of fields is described in [another notebook](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-06-14-data-structures.ipynb), and we are not particularly focusing on them here. Here's some examples how to store data to element."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 10,
    "metadata": {
     "collapsed": false
    },
@@ -343,51 +337,43 @@
     {
      "data": {
       "text/plain": [
-       "Dict{ASCIIString,JuliaFEM.FieldSet} with 4 entries:\n",
-       "  \"geometry\"         => JuliaFEM.FieldSet(\"geometry\",JuliaFEM.Field[JuliaFEM.Fi…\n",
-       "  \"heat coefficient\" => JuliaFEM.FieldSet(\"heat coefficient\",JuliaFEM.Field[Jul…\n",
-       "  \"displacement\"     => JuliaFEM.FieldSet(\"displacement\",JuliaFEM.Field[JuliaFE…\n",
-       "  \"temperature\"      => JuliaFEM.FieldSet(\"temperature\",JuliaFEM.Field[JuliaFEM…"
+       "Dict{ASCIIString,JuliaFEM.Field} with 4 entries:\n",
+       "  \"geometry\"         => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
+       "  \"heat coefficient\" => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
+       "  \"displacement\"     => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…\n",
+       "  \"temperature\"      => JuliaFEM.DefaultDiscreteField([JuliaFEM.TimeStep(0.0,Ju…"
       ]
      },
-     "execution_count": 12,
+     "execution_count": 10,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "using JuliaFEM: FieldSet, Field\n",
+    "using JuliaFEM: Field, TimeStep, Increment\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])  # 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",
-    "element[\"geometry\"] = geometry_fieldset # add fieldset to element\n",
+    "# create one timestep for t=0\n",
+    "# is equivalent to \n",
+    "# element[\"geometry\"] = (0.0 => Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])\n",
+    "element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
     "\n",
-    "temperature_fieldset = FieldSet(\"temperature\")\n",
-    "push!(temperature_fieldset, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
-    "push!(temperature_fieldset, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
-    "element[\"temperature\"] = temperature_fieldset\n",
+    "# create two timesteps, one for t=0.0 and another for t=1.0\n",
+    "element[\"temperature\"] = (0.0, [0.0, 0.0, 0.0, 0.0]), (1.0, [1.0, 2.0, 3.0, 4.0])\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",
-    "element[\"displacement\"] = displacement_fieldset\n",
+    "# create two timesteps, one is zeros at time t=0.0 and another with t=1.0\n",
+    "element[\"displacement\"] = (\n",
+    "    0.0 => Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]],\n",
+    "    1.0 => Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]])\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",
-    "element[\"heat coefficient\"] = heat_coefficient_fieldset\n",
+    "# create two timesteps, one for t=0.0 and another for t=1.0\n",
+    "element[\"heat coefficient\"] = (0.0, 2), (1.0, 3)\n",
     "element.fields"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 11,
    "metadata": {
     "collapsed": false
    },
@@ -400,7 +386,7 @@
        " 0.5"
       ]
      },
-     "execution_count": 13,
+     "execution_count": 11,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -414,7 +400,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 12,
    "metadata": {
     "collapsed": false
    },
@@ -426,7 +412,7 @@
        " 0.0  1.0"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 12,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -434,7 +420,54 @@
    "source": [
     "# interpolate derivatives works too\n",
     "dbasis = grad(basis)\n",
-    "dbasis(\"temperature\", [0.0, 0.0], 0.5)  # temperature gradient at mid point of element at time t=0.5"
+    "# temperature gradient at mid point of element at time t=0.5\n",
+    "dbasis(\"temperature\", [0.0, 0.0], 0.5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2.5"
+      ]
+     },
+     "execution_count": 13,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# interpolating scalar -> scalar.\n",
+    "basis(\"heat coefficient\", [0.0, 0.0], 0.5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "3"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "# 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)"
    ]
   },
   {
@@ -443,52 +476,6 @@
    "metadata": {
     "collapsed": false
    },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "2.5"
-      ]
-     },
-     "execution_count": 15,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "# interpolating scalar -> scalar.\n",
-    "basis(\"heat coefficient\", [0.0, 0.0], 0.5)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 16,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "3"
-      ]
-     },
-     "execution_count": 16,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "# 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": 17,
-   "metadata": {
-    "collapsed": false
-   },
    "outputs": [
     {
      "data": {
@@ -498,7 +485,7 @@
        " 0.0703125  0.0078125"
       ]
      },
-     "execution_count": 17,
+     "execution_count": 15,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -523,7 +510,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 16,
    "metadata": {
     "collapsed": false
    },
@@ -551,7 +538,7 @@
        " 1.0)"
       ]
      },
-     "execution_count": 18,
+     "execution_count": 16,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -568,12 +555,12 @@
     "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",
+    "    basis = get_basis(element)\n",
+    "    dbasis = grad(basis)  # gradient is with respect to \"geometry\" field\n",
+    "    detJ = det(basis)\n",
+    "    M += w*rho*basis(xi)'*basis(xi)*detJ(xi)  # mass matrix\n",
+    "    K += w*k*dbasis(xi)'*dbasis(xi)*detJ(xi)  # stiffness matrix\n",
+    "    f += w*basis(xi)'*q*detJ(xi)  # force vector\n",
     "end\n",
     "K, M, f"
    ]
@@ -599,7 +586,7 @@
     "\\begin{align}\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",
+    "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}}\n",
     "\\end{align}\n",
     "\n",
     "Weak form is, find $u\\in\\mathcal{U}$ such that\n",
@@ -617,7 +604,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 17,
    "metadata": {
     "collapsed": false
    },
@@ -637,7 +624,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 18,
    "metadata": {
     "collapsed": false
    },
@@ -645,10 +632,10 @@
     {
      "data": {
       "text/plain": [
-       "size (generic function with 76 methods)"
+       "size (generic function with 83 methods)"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -666,7 +653,7 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
     "    if !haskey(element, \"temperature\")\n",
-    "        element[\"temperature\"] = FieldSet()\n",
+    "        element[\"temperature\"] = []\n",
     "    end\n",
     "    DC2D4(element, integration_points)\n",
     "end\n",
@@ -681,7 +668,7 @@
     "    integration_points = [\n",
     "        IntegrationPoint([0.0], 2.0)]\n",
     "    if !haskey(element, \"temperature\")\n",
-    "        element[\"temperature\"] = FieldSet()\n",
+    "        element[\"temperature\"] = []\n",
     "    end\n",
     "    DC2D2(element, integration_points)\n",
     "end\n",
@@ -699,7 +686,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 19,
    "metadata": {
     "collapsed": false
    },
@@ -710,7 +697,7 @@
        "has_force_vector (generic function with 3 methods)"
       ]
      },
-     "execution_count": 21,
+     "execution_count": 19,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -722,23 +709,23 @@
     "function JuliaFEM.get_mass_matrix(equation::DC2D4, ip, time)\n",
     "    element = get_element(equation)\n",
     "    basis = get_basis(element)\n",
-    "    ρ = basis(\"density\")\n",
-    "    return ρ(ip, time) * basis(ip, time)'*basis(ip, time)\n",
+    "    rho = basis(\"density\", ip, time)\n",
+    "    return rho * basis(ip, time)'*basis(ip, time)\n",
     "end\n",
     "\"\"\" 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\")\n",
-    "    return k(ip, time) * dbasis(ip, time)'*dbasis(ip, time)\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\")\n",
-    "    return f(ip, time)*basis(ip, time)'\n",
+    "    f = basis(\"temperature load\", ip, time)\n",
+    "    return f*basis(ip, time)'\n",
     "end\n",
     "JuliaFEM.has_mass_matrix(equation::DC2D4) = true\n",
     "JuliaFEM.has_stiffness_matrix(equation::DC2D4) = true\n",
@@ -748,8 +735,8 @@
     "function JuliaFEM.get_force_vector(equation::DC2D2, ip, time)\n",
     "    element = get_element(equation)\n",
     "    basis = get_basis(element)\n",
-    "    g = basis(\"temperature flux\")\n",
-    "    return g(ip, time)*basis(ip, time)'\n",
+    "    g = basis(\"temperature flux\", ip, time)\n",
+    "    return g*basis(ip, time)'\n",
     "end\n",
     "JuliaFEM.has_force_vector(equation::DC2D2) = true"
    ]
@@ -769,7 +756,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 20,
    "metadata": {
     "collapsed": false
    },
@@ -778,27 +765,23 @@
     "# if no time is given when defining fields it defaults to 0.0.\n",
     "# create volume element\n",
     "element = Quad4([1, 2, 3, 4])\n",
-    "fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
-    "fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n",
-    "fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n",
-    "fieldset4 = FieldSet(\"density\", [Field(36.0)])\n",
-    "element[\"geometry\"] = fieldset1\n",
-    "element[\"temperature thermal conductivity\"] = fieldset2\n",
-    "element[\"temperature load\"] = fieldset3\n",
-    "element[\"density\"] = fieldset4\n",
+    "element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
+    "element[\"temperature thermal conductivity\"] = 6.0\n",
+    "element[\"temperature load\"] = [12.0, 12.0, 12.0, 12.0]\n",
+    "element[\"density\"] = 36.0\n",
     "equation = DC2D4(element)\n",
     "\n",
     "# create boundary element with \n",
     "boundary_element = Seg2([1, 2])\n",
-    "boundary_element[\"geometry\"] = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])])\n",
+    "boundary_element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0]]\n",
     "# linear ramp from 1 to 6 in time 0 to 1\n",
-    "boundary_element[\"temperature flux\"] = FieldSet(\"temperature flux\", [Field(0.0, 0.0), Field(1.0, 6.0)])\n",
+    "boundary_element[\"temperature flux\"] = (0.0 => 0.0, 1.0 => 6.0)\n",
     "boundary_equation = DC2D2(boundary_element);"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 21,
    "metadata": {
     "collapsed": false
    },
@@ -813,7 +796,7 @@
        " 2.0  1.0  2.0  4.0"
       ]
      },
-     "execution_count": 23,
+     "execution_count": 21,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -828,7 +811,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 22,
    "metadata": {
     "collapsed": false
    },
@@ -843,7 +826,7 @@
        " -1.0  -2.0  -1.0   4.0"
       ]
      },
-     "execution_count": 24,
+     "execution_count": 22,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -854,7 +837,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 23,
    "metadata": {
     "collapsed": false
    },
@@ -869,7 +852,7 @@
        " 3.0"
       ]
      },
-     "execution_count": 25,
+     "execution_count": 23,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -887,7 +870,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 26,
+   "execution_count": 24,
    "metadata": {
     "collapsed": false
    },
@@ -900,7 +883,7 @@
        " 1.0"
       ]
      },
-     "execution_count": 26,
+     "execution_count": 24,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -922,7 +905,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 25,
    "metadata": {
     "collapsed": false
    },
@@ -935,14 +918,14 @@
        " 1.0"
       ]
      },
-     "execution_count": 27,
+     "execution_count": 25,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "local_assembly = initialize_local_assembly()\n",
-    "calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\")\n",
+    "calculate_local_assembly!(local_assembly, boundary_equation, \"temperature\", 1.0)\n",
     "b = zeros(4)\n",
     "b[fdofs] = local_assembly.force_vector\n",
     "u = zeros(4)\n",
@@ -953,12 +936,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Some quick notes. Passing time all around is not strictly needed and defaults to $t=+\\inf$, when the most recent field is picked from `FieldSet`. Also when defining `Field` without time, it defaults to $t=0$. It's however good practice to interpolate fields inside functions using time parameter which makes the equations more general. If we want to solve the last equation in time $t=0.5$ we should get $T=0.5$:"
+    "Some quick notes. Passing time all around is not strictly needed and defaults to $t=0$. Also when defining `Field` without time, it defaults to $t=0$. It's however good practice to interpolate fields inside functions using time parameter which makes the equations more general. If we want to solve the last equation in time $t=0.5$ we should get $T=0.5$:"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 28,
+   "execution_count": 26,
    "metadata": {
     "collapsed": false
    },
@@ -971,7 +954,7 @@
        " 0.5"
       ]
      },
-     "execution_count": 28,
+     "execution_count": 26,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1009,7 +992,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 29,
+   "execution_count": 27,
    "metadata": {
     "collapsed": false
    },
@@ -1017,10 +1000,10 @@
     {
      "data": {
       "text/plain": [
-       "size (generic function with 77 methods)"
+       "size (generic function with 84 methods)"
       ]
      },
-     "execution_count": 29,
+     "execution_count": 27,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1038,9 +1021,8 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
     "    if !haskey(element, \"temperature\")\n",
-    "        element[\"temperature\"] = FieldSet()\n",
     "        # Initial configuration needs to be defined if using autodiff\n",
-    "        push!(element[\"temperature\"], Field([0.0, 0.0, 0.0, 0.0]))\n",
+    "        element[\"temperature\"] = zeros(4)\n",
     "    end\n",
     "    DC2D4NL(element, integration_points)\n",
     "end\n",
@@ -1056,7 +1038,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 30,
+   "execution_count": 28,
    "metadata": {
     "collapsed": false
    },
@@ -1067,7 +1049,7 @@
        "has_potential_energy (generic function with 2 methods)"
       ]
      },
-     "execution_count": 30,
+     "execution_count": 28,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1080,8 +1062,9 @@
     "    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*vecdot(∇T, ∇T)\n",
+    "    c = basis(\"temperature nonlinearity coefficient\", ip, time)\n",
+    "    gradT = grad(basis)(\"temperature\", ip, time, variation)\n",
+    "    Wint = 1/2 * k*(1 + c*T) * vecdot(gradT, gradT)\n",
     "    Wext = f*T\n",
     "    return Wint - Wext\n",
     "end\n",
@@ -1097,7 +1080,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 31,
+   "execution_count": 29,
    "metadata": {
     "collapsed": false
    },
@@ -1106,14 +1089,15 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "increment  1, |du| =  1.41421, |r| =    2.121\n",
-      "increment  2, |du| =  0.42426, |r| =    2.121\n",
-      "increment  3, |du| =  0.04657, |r| =    2.121\n",
-      "increment  4, |du| =  0.00057, |r| =    2.121\n",
-      "increment  5, |du| =  0.00000, |r| =    2.121\n",
-      "elapsed time: 0.617104273 seconds\n",
-      "error: 1.5543122344752192e-15\n",
-      "temperature at free end: [0.6666666666666682,0.6666666666666682], should be 0.6666666666666666\n"
+      "increment  1, |du| = 1.4142135624\n",
+      "increment  2, |du| = 0.4242640687\n",
+      "increment  3, |du| = 0.0465655685\n",
+      "increment  4, |du| = 0.0005747960\n",
+      "increment  5, |du| = 0.0000000876\n",
+      "increment  6, |du| = 0.0000000000\n",
+      "increment  7, |du| = 0.0000000000\n",
+      "elapsed time: 1.142529838 seconds\n",
+      "error: 0.0\n"
      ]
     }
    ],
@@ -1121,47 +1105,33 @@
     "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",
-    "    element[\"geometry\"] = fieldset1\n",
-    "    element[\"temperature thermal conductivity\"] = fieldset2\n",
-    "    element[\"temperature load\"] = fieldset3\n",
-    "    element[\"temperature nodal load\"] = fieldset4\n",
+    "    element[\"geometry\"] = Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]\n",
+    "    element[\"temperature thermal conductivity\"] = 6.0\n",
+    "    element[\"temperature load\"] = [0.0, 0.0, 0.0, 0.0]\n",
+    "    element[\"temperature nodal load\"] = [3.0, 3.0, 0.0, 0.0]\n",
+    "    element[\"temperature nonlinearity coefficient\"] = [1.0, 1.0, 1.0, 1.0]\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() # create workspace for local matrices\n",
     "    T = zeros(4) # create workspace for solution vector\n",
-    "    ΔT = zeros(4) # \n",
-    "    fd = [1, 2]\n",
+    "    dT = zeros(4) # \n",
+    "    fd = [1, 2] # free dofs\n",
     "    tic()\n",
     "    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.\n",
-    "    for i=1:5\n",
+    "    for i=1:10\n",
     "        calculate_local_assembly!(la, equation, \"temperature\")  # 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",
+    "        dT[fd] = la.stiffness_matrix[fd,fd] \\ la.force_vector[fd]  # <-- note sign convention, more on this below\n",
+    "        T += dT\n",
+    "        push!(element[\"temperature\"], T) # add new increment to model\n",
     "        # print some convergence information\n",
-    "        @printf(\"increment %2d, |du| = %8.5f, |r| = %8.3f\\n\", i, norm(ΔT), norm(b[fd]))\n",
+    "        @printf(\"increment %2d, |du| = %12.10f\\n\", i, norm(dT))\n",
+    "        err = last(element[\"temperature\"])[1] - 2/3\n",
+    "        isapprox(err, 0.0) && break\n",
     "    end\n",
     "    toc()\n",
-    "    err = element[\"temperature\"][end].values[1] - 2/3\n",
+    "    err = last(element[\"temperature\"])[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()"
    ]
@@ -1177,7 +1147,7 @@
     "\\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",
+    "\\Rightarrow \\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",
@@ -1205,14 +1175,14 @@
     "\\end{equation}\n",
     "where\n",
     "\\begin{equation}\n",
-    "\\delta\\mathcal{W}_{\\mathrm{int}}=\\int\\mathbf{S}:\\delta\\mathbf{E}\\,\\mathrm{d}V_{0}=\\int_{\\Omega_{0}}\\mathbf{P}\\cdot\\mbox{Grad}\\delta\\mathbf{u}\\,\\mathrm{d}V_{0}\n",
+    "\\delta\\mathcal{W}_{\\mathrm{int}}=\\int_{\\Omega_{0}}\\mathbf{S}:\\delta\\mathbf{E}\\,\\mathrm{d}V_{0}=\\int_{\\Omega_{0}}\\mathbf{P}\\cdot\\mbox{Grad}\\delta\\mathbf{u}\\,\\mathrm{d}V_{0}\n",
     "\\end{equation}\n",
     "and $\\boldsymbol{b}_{0}$ is a volume load defined in reference configuration."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 32,
+   "execution_count": 30,
    "metadata": {
     "collapsed": false
    },
@@ -1223,7 +1193,7 @@
        "has_residual_vector (generic function with 3 methods)"
       ]
      },
-     "execution_count": 32,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1248,10 +1218,8 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
     "    if !haskey(element, \"displacement\")\n",
-    "        element[\"displacement\"] = FieldSet()\n",
     "        # initial field must be defined if using autodiff\n",
-    "        u0 = Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]])\n",
-    "        push!(element[\"displacement\"], u0)\n",
+    "        element[\"displacement\"] = zeros(2, 4)\n",
     "    end\n",
     "    CPS4(element, integration_points)\n",
     "end\n",
@@ -1280,7 +1248,6 @@
     "    S = λ*trace(E)*I + 2*μ*E\n",
     "    J = det(F)\n",
     "    σ = J^-1 * F*S*F'\n",
-    "    #push!(ip[\"cauchy stress\"], Field(time, σ))\n",
     "    P = F*S\n",
     "\n",
     "    # residual vector\n",
@@ -1296,12 +1263,12 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "There's a tiny Newton solver which can be used to test single element models to verify formulations. It makes the life a bit easier."
+    "We have a tiny Newton solver which can be used to test single element models to verify formulations. It makes the life a bit easier."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 31,
    "metadata": {
     "collapsed": false
    },
@@ -1312,34 +1279,19 @@
      "text": [
       "displacement at tip: -8.773031198197748\n"
      ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "Success :: (line:-1) :: fact was true\n",
-       "  Expression: disp --> roughly(-8.77303119819776)\n",
-       "    Expected: -8.77303119819776\n",
-       "    Occurred: -8.773031198197748"
-      ]
-     },
-     "execution_count": 33,
-     "metadata": {},
-     "output_type": "execute_result"
     }
    ],
    "source": [
+    "using JuliaFEM.Test\n",
     "using JuliaFEM: get_integration_points, solve!\n",
     "\n",
     "function run_simulation_2()\n",
     "    # create model -- start\n",
-    "    # fields can be created less verbose way. just keep in mind that if time is not defined\n",
-    "    # it defaults to 0.0. fields can be set \"directly\" to the element but it creates a fieldset\n",
-    "    # anyway\n",
     "    element = Quad4([1, 2, 3, 4])\n",
-    "    element[\"geometry\"] = Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]])\n",
-    "    element[\"youngs modulus\"] = Field(500.0)\n",
-    "    element[\"poissons ratio\"] = Field(0.3)\n",
-    "    element[\"displacement volume load\"] = Field(Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]])\n",
+    "    element[\"geometry\"] = Vector[[0.0,0.0], [10.0,0.0], [10.0,1.0], [0.0,1.0]]\n",
+    "    element[\"youngs modulus\"] = 500.0\n",
+    "    element[\"poissons ratio\"] = 0.3\n",
+    "    element[\"displacement volume load\"] = Vector[[0.0,-10.0], [0.0,-10.0], [0.0,-10.0], [0.0,-10.0]]\n",
     "    equation = CPS4(element)\n",
     "    # create model -- end\n",
     "\n",
@@ -1348,7 +1300,7 @@
     "    disp = get_basis(element)(\"displacement\", [1.0, 1.0])[2]\n",
     "    println(\"displacement at tip: $disp\")\n",
     "    # verified using Code Aster.\n",
-    "    @fact disp --> roughly(-8.77303119819776E+00)\n",
+    "    @test isapprox(disp, -8.77303119819776E+00)\n",
     "end\n",
     "run_simulation_2()"
    ]
@@ -1366,31 +1318,25 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 34,
+   "execution_count": 32,
    "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": 34,
-     "metadata": {},
-     "output_type": "execute_result"
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: solution = [1.0000000000000002,1.0000000000000002]\n"
+     ]
     }
    ],
    "source": [
     "using JuliaFEM: LocalAssembly, get_integration_points, initialize_local_assembly!\n",
     "\n",
     "function JuliaFEM.calculate_local_assembly!(assembly::LocalAssembly, equation::DC2D4,\n",
-    "                                            unknown_field_name::ASCIIString, time::Number=Inf,\n",
-    "                                            problem=nothing) # <= notice problem\n",
+    "                                            unknown_field_name::ASCIIString, time::Number=0.0,\n",
+    "                                            problem=nothing)\n",
     "\n",
     "    # zero all workspace. always start with this command, otherwise\n",
     "    # local matrices are added to the old ones.\n",
@@ -1428,7 +1374,9 @@
     "    la = initialize_local_assembly()\n",
     "    calculate_local_assembly!(la, equation, \"temperature\")\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",
+    "    solution = la.stiffness_matrix[free_dofs, free_dofs] \\ la.force_vector[free_dofs]\n",
+    "    info(\"solution = $solution\")\n",
+    "    @test isapprox(solution, [1.0, 1.0])\n",
     "end\n",
     "\n",
     "test_local_assembly()"
@@ -1438,7 +1386,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Note that function takes also problem as input, which allows to access to rest of the equations and elements. From that in the next section."
+    "Note that function takes also problem as input, which allows to access to rest of the equations and elements. More info about that in the next section."
    ]
   },
   {
@@ -1449,12 +1397,12 @@
     "\n",
     "The main purpose of `Problem` is to create a mapping between element types and field equations, so that solver knows how to construct equations. Loosely speaking `Problem` is a set of equations and instructions how they can be created from elements. Every time new element is pushed to problem, element mapping creates corresponding equations ready for integration.\n",
     "\n",
-    "Notice that there is no any limits of how many different equations can be set to one element, so it's totally possible for example solve displacement and temperature equations with same element and even so that problems are loosely coupled, i.e. using temperature field to calculate stresses in mechanical problem."
+    "Notice that there is no any limits of how many \"equations\" can be set to one element, so it's totally possible for example solve displacement and temperature equations using same element and even so that problems are loosely coupled, i.e. using temperature field to calculate stresses in mechanical problem."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 35,
+   "execution_count": 33,
    "metadata": {
     "collapsed": false
    },
@@ -1465,7 +1413,7 @@
        "PlaneHeatProblem"
       ]
      },
-     "execution_count": 35,
+     "execution_count": 33,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1501,7 +1449,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 36,
+   "execution_count": 34,
    "metadata": {
     "collapsed": false
    },
@@ -1509,10 +1457,10 @@
     {
      "data": {
       "text/plain": [
-       "call (generic function with 1319 methods)"
+       "call (generic function with 1271 methods)"
       ]
      },
-     "execution_count": 36,
+     "execution_count": 34,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1541,7 +1489,7 @@
     "    Cu       = g\n",
     "\n",
     "\"\"\"\n",
-    "function call(solver::SimpleSolver, time::Number=Inf)\n",
+    "function call(solver::SimpleSolver, time::Number=0.0)\n",
     "    p1, p2 = get_problems(solver)\n",
     "\n",
     "    ga1 = initialize_global_assembly(p1)\n",
@@ -1559,22 +1507,21 @@
     "    b = [b1; b2]\n",
     "\n",
     "    # solve problem\n",
-    "    nz = unique(rowvals(A))  # here we remove any zero rows\n",
+    "    nz = unique(rowvals(A))  # remove any zero rows\n",
     "    x = zeros(b)\n",
     "    x[nz] = lufact(A[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",
+    "    x1 = full(x[1:length(b1)])\n",
+    "    x2 = full(x[length(b1)+1:end])\n",
     "\n",
     "    # update field for elements in problem 1\n",
     "    for equation in get_equations(p1)\n",
     "        element = get_element(equation)\n",
     "        field_name = get_unknown_field_name(p1)\n",
     "        gdofs = get_gdofs(p1, equation)\n",
-    "        element_solution = full(x1[gdofs])\n",
-    "        field = Field(time, element_solution)\n",
-    "        push!(element[field_name], field)\n",
+    "        local_solution = reshape(x1[gdofs], size(equation))\n",
+    "        push!(element[field_name], local_solution)\n",
     "    end\n",
     "\n",
     "    # update field for elements in problem 2 (Dirichlet boundary)\n",
@@ -1582,9 +1529,8 @@
     "        element = get_element(equation)\n",
     "        field_name = get_unknown_field_name(p2)\n",
     "        gdofs = get_gdofs(p2, equation)\n",
-    "        element_solution = full(x2[gdofs])\n",
-    "        field = Field(time, element_solution)\n",
-    "        push!(element[field_name], field)\n",
+    "        local_solution = reshape(x2[gdofs], size(equation))\n",
+    "        push!(element[field_name], local_solution)\n",
     "    end\n",
     "end"
    ]
@@ -1600,7 +1546,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 35,
    "metadata": {
     "collapsed": false
    },
@@ -1609,25 +1555,12 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "28-Oct 04:27:34:INFO:root:assembling problem for temperature\n",
-      "28-Oct 04:27:34:INFO:root:dimension of unknown field: 1, problem dofs: 4\n",
-      "28-Oct 04:27:35:INFO:root:assembling problem for reaction force\n",
-      "28-Oct 04:27:35:INFO:root:dimension of unknown field: 1, problem dofs: 4\n",
-      "28-Oct 04:27:35:INFO:root:Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n"
+      "INFO: assembling problem for temperature\n",
+      "INFO: dimension of unknown field: 1, problem dofs: 4\n",
+      "INFO: assembling problem for reaction force\n",
+      "INFO: dimension of unknown field: 1, problem dofs: 4\n",
+      "INFO: Temperature at point X = [0.5,0.0] is T = 100.00000000000003\n"
      ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "Success :: (line:-1) :: fact was true\n",
-       "  Expression: T --> roughly(100.0)\n",
-       "    Expected: 100.0\n",
-       "    Occurred: 100.00000000000003"
-      ]
-     },
-     "execution_count": 37,
-     "metadata": {},
-     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -1638,8 +1571,6 @@
     "\"\"\"\n",
     "function get_heatproblem()\n",
     "    el1 = Quad4([1, 2, 3, 4])\n",
-    "    # these might look like normal values but believe me, they\n",
-    "    # are fields with temporal and spatial dimension\n",
     "    el1[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
     "    el1[\"temperature thermal conductivity\"] = 6.0\n",
     "    el1[\"temperature load\"] = zeros(4)\n",
@@ -1647,17 +1578,7 @@
     "\n",
     "    el2 = Seg2([1, 2])\n",
     "    el2[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0]]\n",
-    "    # Boundary load, linear ramp 0 -> 600 at time 0 -> 1\n",
-    "    # yet another simplification, if field is given as a tuple,\n",
-    "    # multiple fields are created. there is 1 second time step between\n",
-    "    # each field. So the following is basically same as\n",
-    "    # fieldset = FieldSet(\"temperature flux\")\n",
-    "    # field1 = Field(0.0, 0.0)\n",
-    "    # field2 = Field(1.0, 600.0)\n",
-    "    # push!(fieldset, field1)\n",
-    "    # push!(fieldset, field2)\n",
-    "    # element[\"temperature flux\"] = fieldset\n",
-    "    el2[\"temperature flux\"] = (0.0, 600.0)\n",
+    "    el2[\"temperature flux\"] = (0.0 => 0.0, 1.0 => 600.0)\n",
     "\n",
     "    problem1 = PlaneHeatProblem()\n",
     "    push!(problem1, el1)\n",
@@ -1687,28 +1608,20 @@
     "    push!(solver, problem2)\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 = [0.0, -1.0]\n",
     "    el2 = get_element(problem1.equations[2])\n",
     "    basis = get_basis(el2)\n",
     "    X = basis(\"geometry\", xi, 1.0)\n",
-    "    T = basis(\"temperature\", xi, 1.0)\n",
-    "    Logging.info(\"Temperature at point X = $X is T = $T\")\n",
-    "    @fact T --> roughly(100.0)\n",
+    "    T = basis(\"temperature\", xi, 1.0)[1]\n",
+    "    info(\"Temperature at point X = $X is T = $T\")\n",
+    "    @test isapprox(T, 100.0)\n",
     "end\n",
     "\n",
     "main()"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "In this notebook the basic instructions how to develop JuliaFEM has been given. The most imporant concepts has been considered; how to develop own element with own basis, several ways how to define own equation, and how to finally assemble and calculate the problem using solver. Any comments and/or discussion about technical details, theory, programming, or from life in general is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1717,99 +1630,328 @@
     "\n",
     "In last section of this tutorial we consider some of the more advanced things which may araise when developing own models.\n",
     "\n",
-    "### Boundary element access to parent element + overriding equations in problems\n",
-    "\n",
-    "This kind of situation might happen when one is almost happy for some problem setting, but would like to change just one or two equations from it. For example boundary equation is not satisfying all the requirements and one would like to test something new. \n",
-    "\n",
-    "### Accessing integration points\n",
-    "\n",
-    "### Fields as a function of something.\n",
-    "- statistical variables\n",
-    "- field dependent from another field\n",
-    "- field dependent from time or spatial domain etc.\n"
+    "### Boundary element access to parent element + overriding equations in problems\n"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 36,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "has_potential_energy (generic function with 3 methods)"
+      ]
+     },
+     "execution_count": 36,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "\"\"\" 2-node radiation boundary element. \"\"\"\n",
+    "type DC2D2NL <: Heat\n",
+    "    element :: Seg2\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "end\n",
+    "\n",
+    "function DC2D2NL(element::Seg2)\n",
+    "    # https://pomax.github.io/bezierinfo/legendre-gauss.html\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint([0.0], 128/225),\n",
+    "        IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),\n",
+    "        IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),\n",
+    "        IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),\n",
+    "        IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900)\n",
+    "    ]\n",
+    "    if !haskey(element, \"temperature\")\n",
+    "        element[\"temperature\"] = zeros(2)\n",
+    "    end\n",
+    "    DC2D2NL(element, integration_points)\n",
+    "end\n",
+    "\n",
+    "function Base.size(equation::DC2D2NL)\n",
+    "    return (1, 2)\n",
+    "end\n",
+    "\n",
+    "function JuliaFEM.get_potential_energy(equation::DC2D2NL, ip, time; variation=nothing)\n",
+    "    element = get_element(equation)\n",
+    "    basis = get_basis(element)\n",
+    "    T = basis(\"temperature\", ip, time, variation)[1]\n",
+    "    T_ext = basis(\"temperature external\", ip, time)[1]\n",
+    "    coeff = basis(\"temperature coefficient\", ip, time)[1]\n",
+    "    q0 = coeff*(T_ext^4 - T^4)\n",
+    "    Wint = 0.0\n",
+    "    Wext = q0*T\n",
+    "    return Wint - Wext\n",
+    "end\n",
+    "\n",
+    "function JuliaFEM.has_potential_energy(equation::DC2D2NL)\n",
+    "    return true\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Defining new equation mapping to old problem is one line command. Here's the equation `DC2D2` is replaced to $\\rightarrow$ `DC2DCNL`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 37,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "INFO: residual norm: 0.7071067811865476\n",
+      "INFO: residual norm: 2.2097084800638154e-9\n",
+      "INFO: residual norm: 2.09345658541016e-16\n",
+      "INFO: Temperature at the midpoint of boundary element: 0.4999999984375\n"
+     ]
+    }
+   ],
+   "source": [
+    "function run_radiation_model()\n",
+    "\n",
+    "    # create elements\n",
+    "    N = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
+    "    element1 = Quad4([1, 2, 3, 4])\n",
+    "    element1[\"geometry\"] = Vector[N[1], N[2], N[3], N[4]]\n",
+    "    element1[\"temperature thermal conductivity\"] = 6.0\n",
+    "    element1[\"temperature load\"] = [0.0, 0.0, 0.0, 0.0]\n",
+    "    element1[\"density\"] = 36.0\n",
+    "    element1[\"temperature nonlinearity coefficient\"] = [0.0, 0.0, 0.0, 0.0]\n",
+    "    element2 = Seg2([1, 2])\n",
+    "    element2[\"geometry\"] = Vector[N[1], N[2]]\n",
+    "    #element2[\"emissivity\"] = 0.5\n",
+    "    element2[\"temperature coefficient\"] = 3.0e-8 # ~ 5.7e-8 * 0.5\n",
+    "    element2[\"temperature external\"] = 100.0\n",
+    "\n",
+    "    # set initial conditions\n",
+    "    element1[\"temperature\"] = zeros(4)\n",
+    "    element2[\"temperature\"] = zeros(2)\n",
+    "    # create equations\n",
+    "\n",
+    "    equation1 = DC2D4NL(element1)\n",
+    "    equation2 = DC2D2NL(element2)\n",
+    "    #equations = [equation1, equation2]\n",
+    "\n",
+    "    # create workspace for local matrices\n",
+    "    la1 = initialize_local_assembly(equation1)\n",
+    "    la2 = initialize_local_assembly(equation2)\n",
+    "\n",
+    "    #solve!(problem, free_dofs; max_iterations=5, dump_matrices=true)\n",
+    "\n",
+    "    T = zeros(4)\n",
+    "    dT = zeros(4)\n",
+    "    free_dofs = [1, 2]\n",
+    "    K = zeros(4,4)\n",
+    "    f = zeros(4)\n",
+    "    for i=1:20\n",
+    "        fill!(K, 0.0)\n",
+    "        fill!(f, 0.0)\n",
+    "        calculate_local_assembly!(la1, equation1, \"temperature\")\n",
+    "        calculate_local_assembly!(la2, equation2, \"temperature\")\n",
+    "\n",
+    "        K += la1.stiffness_matrix\n",
+    "        f += la1.force_vector\n",
+    "        #info(\"f(linear) = \\n$(la1.force_vector)\")\n",
+    "        #info(\"K(linear) = \\n$(la1.stiffness_matrix[free_dofs, free_dofs])\")\n",
+    "\n",
+    "        K[free_dofs, free_dofs] += la2.stiffness_matrix\n",
+    "        f[free_dofs] += la2.force_vector\n",
+    "        #info(\"K(nonlinear) = \\n$(la2.stiffness_matrix)\")\n",
+    "        #info(\"f(nonlinear) = \\n$(la2.force_vector)\")\n",
+    "\n",
+    "        dT[free_dofs] = K[free_dofs, free_dofs] \\ f[free_dofs]\n",
+    "        T += dT\n",
+    "\n",
+    "        push!(element1[\"temperature\"], T)\n",
+    "        push!(element2[\"temperature\"], T[free_dofs])\n",
+    "        T_mid = get_basis(element2)(\"temperature\", [0.0])[1]\n",
+    "        info(\"residual norm: $(norm(dT))\")\n",
+    "        norm(dT) < 1.0e-12 && break\n",
+    "    end\n",
+    "    #solve!(problem, free_dofs; max_iterations=50)\n",
+    "\n",
+    "    T_mid = get_basis(element2)(\"temperature\", [0.0])[1]\n",
+    "    info(\"Temperature at the midpoint of boundary element: $T_mid\")\n",
+    "end\n",
+    "\n",
+    "run_radiation_model()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Accessing integration points\n",
+    "\n",
+    "It's possible to save internal variables to integration points. This is described [constitutive modelling tutorial](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-11-03-constitutive-modelling-using-juliafem.ipynb).\n",
+    "\n",
+    "### Statistical fields\n",
+    "\n",
+    "It has been already discussed in [data structures](https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/docs/tutorials/2015-06-14-data-structures.ipynb) how to make fields dependent from other fields. Another idea is to make field dependent from external resources (data from logging system / internet / random.org). For example, sampling data from normal distribution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "first (generic function with 13 methods)"
+      ]
+     },
+     "execution_count": 38,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "type RandomField <: JuliaFEM.DiscreteField\n",
+    "    mu :: Float64\n",
+    "    std :: Float64\n",
+    "end\n",
+    "Base.first(field::RandomField) = Increment(randn(4).*field.std^2 + field.mu)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "4-element JuliaFEM.Increment{Float64}:\n",
+       "  9.18727\n",
+       " 10.2647 \n",
+       "  7.93236\n",
+       " 12.0819 "
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "RandomField(10.0, 2.0)(0.0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 40,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "function calc_stochastic(N=500)\n",
+    "    solutions = Vector[]\n",
+    "    element = Quad4([1, 2, 3, 4])\n",
+    "    element[\"geometry\"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]\n",
+    "    element[\"temperature load\"] = [12.0, 12.0, 12.0, 12.0]\n",
+    "    element[\"density\"] = 36.0\n",
+    "    element[\"temperature thermal conductivity\"] = RandomField(6.0, 1.0)\n",
+    "    equation = DC2D4(element)\n",
+    "    la = initialize_local_assembly()\n",
+    "    fd = [1,2]\n",
+    "    for i=1:N\n",
+    "        calculate_local_assembly!(la, equation, \"temperature\", 0.0)\n",
+    "        push!(solutions, la.stiffness_matrix[fd, fd] \\ la.force_vector[fd])\n",
+    "    end\n",
+    "    return solutions\n",
+    "end\n",
+    "solutions = calc_stochastic();"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2-element Array{Float64,1}:\n",
+       " 0.997586\n",
+       " 0.998368"
+      ]
+     },
+     "execution_count": 41,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "mean(solutions)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 42,
    "metadata": {
     "collapsed": true
    },
    "outputs": [],
    "source": [
-    "using JuliaFEM: get_default_integration_points\n",
-    "\"\"\" 2-node radiation boundary element. \"\"\"\n",
-    "type DC2D2RAD <: Heat\n",
-    "    element :: Seg2\n",
-    "    integration_points :: Array{IntegrationPoint, 1}\n",
-    "end\n",
-    "function DC2D2RAD(element::Seg2)\n",
-    "    integration_points = [\n",
-    "        IntegrationPoint([0.0], 2.0)]\n",
-    "    if !haskey(element, \"temperature\")\n",
-    "        element[\"temperature\"] = FieldSet()\n",
-    "        push!(element[\"temperature\"], Field([0.0, 0.0]))\n",
-    "    end\n",
-    "    DC2D2RAD(element, integration_points)\n",
-    "end\n",
-    "Base.size(equation::DC2D2RAD) = (1, 2)\n",
-    "\n",
-    "\"\"\" Calculate potential energy caused by radiation.\n",
-    "https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_constant\n",
-    "\"\"\"\n",
-    "function JuliaFEM.get_potential_energy(equation::DC2D2RAD, ip, time; variation=nothing)\n",
-    "    element = get_element(equation)\n",
-    "    basis = get_basis(element)\n",
-    "    eps = basis(\"emissivity\", ip, time)\n",
-    "    #sig = basis(\"stefan-boltzmann constant\", ip, time)\n",
-    "    sig = 5.670367e-8  # i guess stefan-boltzmann constant is constant ;)\n",
-    "    T = basis(\"temperature\", ip, time, variation)\n",
-    "    T_ext = basis(\"temperature external\", ip, time)\n",
-    "    q = eps*sig*((T_ext+273.15)^4 - (T+273.15)^4)\n",
-    "    println(ForwardDiff.value(q*T))\n",
-    "    return q\n",
-    "end\n",
-    "JuliaFEM.has_potential_energy(equation::DC2D2RAD) = true\n",
-    "\n",
-    "Defining new equation mapping to old problem is one line command. Here we replace `DC2D2` $\\rightarrow$ `DC2DCRAD`\n",
-    "\n",
-    "function run_radiation_model()\n",
-    "    # this is the same as before\n",
-    "    element = Quad4([1, 2, 3, 4])\n",
-    "    fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
-    "    fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n",
-    "    fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n",
-    "    fieldset4 = FieldSet(\"density\", [Field(36.0)])\n",
-    "    push!(element, fieldset1)\n",
-    "    push!(element, fieldset2)\n",
-    "    push!(element, fieldset3)\n",
-    "    push!(element, fieldset4)\n",
-    "\n",
-    "    # create boundary element\n",
-    "    boundary_element = Seg2([1, 2])\n",
-    "    push!(boundary_element, FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])]))\n",
-    "    push!(boundary_element, FieldSet(\"emissivity\", [Field(0.5)]))\n",
-    "    push!(boundary_element, FieldSet(\"temperature external\", [Field(20.0)]))\n",
-    "\n",
-    "    # set initial conditions\n",
-    "    push!(element, FieldSet(\"temperature\", [Field([0.0, 0.0, 0.0, 0.0])]))\n",
-    "    push!(boundary_element, FieldSet(\"temperature\", [Field([0.0, 1.0])]))\n",
-    "    \n",
-    "    # create problem, change element mapping\n",
-    "    problem = PlaneHeatProblem()\n",
-    "    problem[Seg2] = DC2D2RAD  # Seg2 was previous DC2D2\n",
-    "    push!(problem, element)\n",
-    "    push!(problem, boundary_element)\n",
-    "\n",
-    "    # run our \"unit test solver\"\n",
-    "    free_dofs = [3, 4]\n",
-    "    solve!(problem, free_dofs; max_iterations=4, dump_matrices=true)\n",
-    "    basis = get_basis(boundary_element)\n",
-    "    T = basis(\"temperature\", [0.0])\n",
-    "    println(\"Temperature at the midpoint of element: $T\")\n",
-    "end\n",
-    "\n",
-    "run_radiation_model()"
+    "using PyPlot"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAADUCAYAAAAcN2ODAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X1UVHX+B/D3DMIwIpIMIIMco9UfiS7Es2lmGavhimXrilH+2FTQNRIXc43M8mxsWCoLu9LxKX/ppsuKT7s+lG2amlYrKGPH3wbJnvUJHBVEEXF4/v7+8DeT08A4I3fmDvl+ncPJvvd7733PJ/XTvXMfFEIIASIiIuqUUu4AREREroyNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKNkoiIiIr2CiJiIisYKMkIiKygo2SiIjICpsbZWNjI5YsWYLExET4+vpCqVRi48aNVtdpbW3F0KFDoVQqkZeX1+mc9evXIywsDGq1GqGhoSgsLLTvExARETmQzY2ypqYGOTk5+O677xAZGQkAUCgUVtdZuXIlLly40OXcNWvWID09HeHh4SgsLMSIESOQmZmJZcuW2fMZiIiIHMbmRhkUFIRLly7hzJkzWL58+V3nX7lyBTk5OcjOzu50ucFgwBtvvIGkpCQUFxdj5syZ2LhxI1588UXk5OTg+vXrtn8KIiIiB7G5UXp4eCAgIAAAYMsLR7KzszFkyBC8+OKLnS4/ePAg6urq8PLLL5uNZ2RkoLGxEXv37rU1GhERkcM45GKekpIS/PnPf0ZBQUGXc3Q6HQAgNjbWbDw6OhpKpRInT550RDQiIiK7SN4ohRCYO3cunn/+eQwfPrzLeXq9Hm5ubvDz8zMb9/DwgEajwcWLF6WORkREZLdeUm9ww4YN+N///V/s2LHD6jyDwQAPD49Ol6lUKhgMBovx2tpafPrppwgJCYFarZYkLxER9TwGgwFnz57F008/bXHAJTVJG+WNGzfw+uuvY+HChRgwYIDVuWq1Gi0tLZ0ua2pq6rQRfvrpp5g2bZokWYmIqOfbtGlTl9fCSEXSRrlixQq0trYiOTkZZ8+eBQBUVVUBAOrq6nD27FkMGDAA7u7u0Gq1aG9vR21trdn/DbS0tKCurg5BQUEW2w8JCQFwuzBhYWFSRneIrKws5Ofnyx3DJj0pK9Cz8jKrYzCrY/SUrOXl5Zg2bZqpLziSpI3ywoULuHbtGoYNG2axLDc3F7m5uTh58iQiIiJM92KWlpZi/PjxpnnHjx9HR0eHafmdjEeZYWFhiI6OljK6Q/j4+PSInEDPygr0rLzM6hjM6hg9KSsAp3wNJ2mjzMzMxHPPPWc2dvnyZcyePRvTp0/Hs88+a+r+CQkJ8PX1xapVq8wa5apVq+Dl5YUJEyZIGY2IiOie2NUoCwsLcf36ddMVqbt27cL58+cB3G6SUVFRiIqKMlvHeAp22LBheOaZZ0zjnp6eyMnJQUZGBpKTkzFu3DgcOXIEmzdvRm5uLh544IHufC4iIiJJ2NUo8/LycO7cOQC3H0m3c+dO7NixAwqFAqmpqejbt69dO58zZw7c3d2Rl5eHXbt2YeDAgSgoKEBmZqZd2yGS0q1bt1BRUSHZ9urr61FWVtatbQwZMgS9e/eWKBER2cOuRnnmzBm7dxASEoKOjo4ul6elpSEtLc3u7fYEKSkpckewWU/KCjg2b0VFBWJiYiTdZne3d+LECad8b9STfh8wq2P0pKzOohC2PI/ORZSVlSEmJsZpf2nQ/cn4+wzYBEDuq6vLAUzj73miH3BmP5D8gQNEPx5hANiciO53fHEzERGRFWyUREREVrBREhERWcFGSUREZAUbJRERkRVslERERFawURIREVlhc6NsbGzEkiVLkJiYCF9fXyiVSmzcuNFsjhACGzZswDPPPIOBAweiT58+CA8PxzvvvIPm5uZOt7t+/XqEhYVBrVYjNDQUhYWF3ftEREREErK5UdbU1CAnJwffffed6RVYCoXCbE5jYyNmzJiBq1evYs6cOfjjH/+I+Ph4LFmyxOwNIUZr1qxBeno6wsPDUVhYiBEjRiAzMxPLli3r5sciIiKShs1P5gkKCsKlS5cQEBCAEydOIC4uzmKOSqXCV199hUcffdQ0NnPmTISEhGDJkiU4cOAAEhISAAAGgwFvvPEGkpKSUFxcbJrb0dGBnJwczJo1i28QISIi2dl8ROnh4YGAgAAAt0+xdsbd3d2sSRpNmjQJAMzeyHDw4EHU1dXh5ZdfNpubkZGBxsZG7N2719ZoREREDuOUi3kuXboEAPDz8zON6XQ6AEBsbKzZ3OjoaCiVSpw8edIZ0YiIiKxySqNctmwZfHx8zL6n1Ov1cHNzM2uewO0jV41GY3o5NBERkZwc/vaQ3NxcHDhwAKtWrTJ7sbPBYICHh0en66hUKhgMBkdHIyIiuiuHNsotW7bgzTffRFpaGmbPnm22TK1Wo6WlpdP1mpqaoFaru9xuVlYWfHx8zMZSUlL4wlEioh+hoqIiFBUVmY3V19c7bf8Oa5SfffYZUlNTkZSUhNWrV1ss12q1aG9vR21trdnp15aWFtTV1SEoKKjLbefn5/MltkRE94nODoS+f8G64znkO8pjx47hueeeQ3x8PIqLi6FUWu4mKioKAFBaWmo2fvz4cXR0dJju1SQiIpKT5I2yvLwcEyZMwE9+8hPs2bMHKpWq03lPPfUUfH19sWrVKrPxVatWwcvLCxMmTJA6GhERkd3sOvVaWFiI69evm65I3bVrF86fPw8AyMzMhEKhwNNPP43r169j4cKF2L17t9n6gwcPNt1n6enpiZycHGRkZCA5ORnjxo3DkSNHsHnzZuTm5vJhA0RE5BLsapR5eXk4d+4cgNuPr9u5cyd27NgBhUKB1NRUdHR0oKqqCgqFAtnZ2Rbrv/TSS2YPJJgzZw7c3d2Rl5eHXbt2YeDAgSgoKEBmZmY3PxYREZE07GqUZ86cueucjo4OuwKkpaUhLS3NrnXox+fWrVtmT26SU3l5udwRiMiFOPw+SiJbVFRUOO0KNiIie7BRkovZBCBM5gwfA3hT5gxE5CrYKMnFhAGQ+x5Znnolou855VmvREREPRUbJRERkRVslERERFawURIREVnBRklERGQFGyUREZEVNjfKxsZGLFmyBImJifD19YVSqcTGjRs7nVteXo7ExER4e3tDo9EgNTUVtbW1nc5dv349wsLCoFarERoaisLCwnv7JERERA5gc6OsqalBTk4OvvvuO9MrsBQKhcW8qqoqjB49Gv/5z3+wdOlSLFiwAHv37sXYsWPR2tpqNnfNmjVIT09HeHg4CgsLMWLECGRmZmLZsmXd/FhERETSsPmBA0FBQbh06RICAgJw4sQJxMXFdTovNzcXBoMBOp0OwcHBAID4+HiMHTsWGzZsQHp6OgDAYDDgjTfeQFJSEoqLiwEAM2fOREdHB3JycjBr1iy+QYSIiGRn8xGlh4cHAgICAABCiC7nbd++HUlJSaYmCQAJCQkIDQ01NUQAOHjwIOrq6vDyyy+brZ+RkYHGxkbs3bvX5g9BRETkKJJezFNdXY2amhrExsZaLIuLi4NOpzP9u/HXP5wbHR0NpVKJkydPShmNiIjonkjaKPV6PQBAq9VaLNNqtairqzN9T6nX6+Hm5gY/Pz+zeR4eHtBoNKaXQxMREclJ0kZpMBgAACqVymKZp6en2RyDwQAPD49Ot6NSqUzziIiI5CTp20PUajUAoLm52WJZU1OT2Ry1Wo2WlpZOt9PU1GSa15msrCz4+PiYjaWkpCAlJeWechMRkesqKipCUVGR2Vh9fb3T9i9pozSecjWegr2TXq+HRqOBu7u7aW57eztqa2vNTr+2tLSgrq4OQUFBXe4nPz8f0dFyv4qJiIicobMDobKyMqe97F3SU68DBgyAv78/SktLLZaVlJSY7r8EgKioKACwmHv8+HF0dHSYzSUiIpKL5I+wmzx5Mvbs2YOqqirT2IEDB1BZWYkpU6aYxp566in4+vpi1apVZuuvWrUKXl5emDBhgtTRiIiI7GbXqdfCwkJcv37ddEXqrl27cP78eQBAZmYm+vbti0WLFmHr1q0YM2YM5s2bh4aGBixfvhwRERGYPn26aVuenp7IyclBRkYGkpOTMW7cOBw5cgSbN29Gbm4uHzZAREQuwa5GmZeXh3PnzgG4/fi6nTt3YseOHVAoFEhNTUXfvn0RHByMw4cPY/78+cjOzoZKpUJSUhLy8vJM308azZkzB+7u7sjLy8OuXbswcOBAFBQUIDMzU7pPSERE1A12NcozZ87YNG/o0KHYt2+fTXPT0tKQlpZmTwwiIiKn4Wu2iIiIrGCjJCIisoKNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKNkoiIiIr2CiJiIiskPQ1W0bHjx9HTk4OSktLUV9fj4EDB+KFF17AggULzN4zWV5ejqysLHz55Zfw8PDAhAkT8Ic//MHstVtE97fbLzAvLy+XOYe5IUOGoHfv3nLHIHIKyRvlqVOnMGrUKAQFBeE3v/kNfH198dVXX2HJkiU4ceIE/va3vwEAqqqqMHr0aPTr1w9Lly5FQ0MDVqxYgVOnTqGkpMTiubBE96ezAIBp06bJG+MHTpw4wXfC0n1D8kb517/+FS0tLdi7dy/CwsIA3H6ea0dHB/785z+jvr4ePj4+yM3NhcFggE6nQ3BwMAAgPj4eY8eOxYYNG5Ceni51NKIebBOAMLlDACgH4FpNm8jRJG+UxlOrAQEBZuOBgYFwc3ODh4cHAGD79u1ISkoyNUkASEhIQGhoKIqLi9koicyEAeARHJEcJL+YZ8aMGejfvz9mzpyJb775BhcuXMCWLVuwevVqZGZmQq1Wo7q6GjU1NYiNjbVYPy4uDjqdTupYRERE90TyI8qgoCB8+eWX+PnPf46oqCjT+OLFi/H2228DAPR6PQBAq9VarK/ValFXV4fW1lZ+T0lERLKTvFFevnwZ48ePBwCsW7cOGo0Ge/bswTvvvIP+/fsjIyMDBsPtK/lUKpXF+p6engAAg8HARklERLKTvFHm5OSguroap0+fRlBQEABg0qRJ6OjowGuvvYaUlBTT95jNzc0W6zc1NQGA2W0kP5SVlQUfHx+zsZSUFKSkpEj1MYiIyEUUFRWhqKjIbKy+vt5p+5e8UR49ehRRUVGmJmk0ceJEbNiwASdPnsTDDz8M4PtTsHfS6/XQaDRWjybz8/N5aToR0X2iswOhsrIyxMTEOGX/kl/M09raivb29k7HAaCtrQ0DBgyAv78/SktLLeaVlJQgMjJS6lhERET3RPJGGR0djbKyMlRWVpqNFxUVwc3NDREREQCAyZMnY8+ePaiqqjLNOXDgACorKzFlyhSpYxEREd0TyU+9/va3v8X27dvx+OOP45VXXoGvry/27NmDffv2IT09HYGBgQCARYsWYevWrRgzZgzmzZuHhoYGLF++HBEREZg+fbrUsYiIiO6J5I0yIiIChw4dwpIlS7B8+XI0NTXhJz/5CXJzc7Fw4ULTvODgYBw+fBjz589HdnY2VCoVkpKSkJeXx6tdneDWrVuoqKiQO4aJqz3LlIjIyCEPRY+Pj8cnn3xy13lDhw7Fvn37HBGB7qKiosJpX4QTEfVkDmmU1JO4yjNEPwbwptwhiIgssFHe91zlGaI89UpErokvbiYiIrKCjZKIiMgKNkoiIiIr2CiJiIisYKMkIiKygo2SiIjICoc1yrKyMjzzzDPQaDTw8vJCeHg4Vq5caTanvLwciYmJ8Pb2hkajQWpqKmprax0ViYiIyG4OuY/yH//4ByZOnIiYmBi89dZb6NOnD/7973+jurraNKeqqgqjR49Gv379sHTpUjQ0NGDFihU4deoUSkpK+Bg7IiJyCZI3yhs3biA1NRUTJ07Etm3bupyXm5sLg8EAnU6H4OBgALcffTd27Fhs2LAB6enpUkcjIiKym+SnXv/yl7/gypUreOeddwAAjY2N6OjosJi3fft2JCUlmZokACQkJCA0NBTFxcVSxyIiIronkjfK/fv3o2/fvrhw4QIefvhheHt7w8fHBy+//DKam5sBANXV1aipqUFsbKzF+nFxcdDpdFLHIiIiuieSN8rKykq0tbVh0qRJGD9+PHbs2IEZM2Zg9erVpvdM6vV6AIBWq7VYX6vVoq6uDq2trVJHIyIispvk31HevHkTt27dwpw5c1BQUAAAmDRpElpaWrBmzRq8/fbbMBgMAACVSmWxvqenJwDAYDDwgh4iIpKd5EeUarUaAJCSkmI2bvz3f/7zn6Y5xlOxd2pqajLbDhERkZwkP6IMCgrCt99+i/79+5uNBwQEAACuXbuGoKAgAN+fgr2TXq+HRqOxejSZlZUFHx8fs7GUlBSL5kxERD1fUVERioqKzMbq6+udtn/JG2VsbCz279+Pqqoq/Nd//Zdp/OLFiwAAf39/BAUFwd/fH6WlpRbrl5SUIDIy0uo+8vPzER3tCu9QJCIiR+vsQKisrAwxMTFO2b/kp16Tk5MBAOvXrzcb/+CDD+Du7o4nn3wSADB58mTs2bMHVVVVpjkHDhxAZWUlpkyZInUsIiKieyL5EWVkZCRmzJiB//mf/0FbWxtGjx6NQ4cOYdu2bVi0aBECAwMBAIsWLcLWrVsxZswYzJs3Dw0NDVi+fDkiIiJMV8cSERHJzSGPsFu9ejUGDhyIDz/8EDt37kRISAgKCgqQmZlpmhMcHIzDhw9j/vz5yM7OhkqlQlJSEvLy8ni1KxERuQyHNMpevXrhrbfewltvvWV13tChQ7Fv3z5HRCAiIpIEX7NFRERkhUOOKInox+r2w0LKy8tlzvG9IUOGoHfv3nLHoB8xNkoissNZAMC0adPkjXGHEydO8HYxcig2SiK6B5sAhMmcoRyA6zRs+vFioySiexAGgEdxdH/gxTxERERWsFESERFZwUZJRERkBRslERGRFU5plO+88w6USiXCw8MtlpWXlyMxMRHe3t7QaDRITU1FbW2tM2IRERHdlcOveq2qqkJubi68vLygUCgslo0ePRr9+vXD0qVL0dDQgBUrVuDUqVMoKSnhM1+JiEh2Dm+UCxYswMiRI9HW1mZxpJibmwuDwQCdTofg4GAAQHx8PMaOHYsNGzYgPT3d0fGIiIiscuip1y+++ALbt29HQUEBhBAWR5Tbt29HUlKSqUkCQEJCAkJDQ1FcXOzIaERERDZxWKNsb2/H3LlzkZ6ejmHDhlksr66uRk1NDWJjYy2WxcXFQafTOSoaERGRzRx26nX16tU4f/48Pv/8806X6/V6AIBWq7VYptVqUVdXh9bWVn5PSUREsnLIEeXVq1dN76PUaDSdzjEYbr+FQKVSWSzz9PQ0m0NERCQXhzTKxYsXw8/PD3Pnzu1yjlqtBgA0NzdbLGtqajKbQ0REJBfJT71WVlZi3bp1KCgoQFVVlWm8qakJLS0tOHfuHPr27Ws65Wo8BXsnvV4PjUbT5WnXrKws+Pj4mI2lpKQgJSVFwk9CRESuoKioCEVFRWZj9fX1Ttu/5I2yuroaHR0dyMzMRGZmpsXyhx56CL/5zW/whz/8Af7+/igtLbWYU1JSgsjIyC73kZ+fz/fPERHdJzo7ECorK0NMTIxT9i95owwPD8fOnTvNbgURQmDx4sW4efMm/vjHP2LQoEEAgMmTJ2Pjxo2oqqoy3SJy4MABVFZW4tVXX5U6GhERkd0kb5QajQbPPvusxXh+fj4A4JlnnjGNLVq0CFu3bsWYMWMwb948NDQ0YPny5YiIiMD06dOljkZERGQ3pz0UXaFQWDxwIDg4GIcPH8agQYOQnZ2NFStWICkpCZ999hlvCyEiIpfg8EfYGR08eLDT8aFDh2Lfvn3OikFERGQXvmaLiIjICjZKIiIiK9goiYiIrGCjJCIisoKNkoiIyAo2SiIiIivYKImIiKxgoyQiIrKCjZKIiMgKyRtlaWkpXnnlFQwbNgx9+vTBgw8+iKlTp6KystJibnl5ORITE+Ht7Q2NRoPU1FTU1tZKHYmIiOieSf4Iu/feew9ff/01pkyZgoiICOj1ehQWFiI6Ohr//Oc/MWzYMABAVVUVRo8ejX79+mHp0qVoaGjAihUrcOrUKZSUlPBZr0RE5BIkb5Svvvoq4uLi0KvX95ueOnUqwsPD8e677+Kjjz4CAOTm5sJgMECn05lesRUfH4+xY8diw4YNSE9PlzoaERGR3SQ/9TpixAizJgkAgwcPxtChQ1FRUWEa2759O5KSkkxNEgASEhIQGhqK4uJiqWMRERHdE6dczCOEwOXLl+Hn5wcAqK6uRk1NDWJjYy3mxsXFQafTOSMWERHRXTmlUW7evBkXL17E1KlTAQB6vR4AoNVqLeZqtVrU1dWhtbXVGdGIiIiscnijrKioQEZGBkaOHIlf/epXAACDwQAAUKlUFvM9PT3N5hAREcnJoY3y0qVLmDBhAvr164dt27ZBoVAAANRqNQCgubnZYp2mpiazOURERHKS/KpXo/r6eowfPx43btzAkSNHEBgYaFpmPOVqPAV7J71eD41GY/X2kKysLPj4+JiNpaSkICUlRaL0RETkKoqKilBUVGQ2Vl9f77T9O6RRNjU1YeLEifj3v/+N/fv3Y8iQIWbLBwwYAH9/f5SWllqsW1JSgsjISKvbz8/PR3R0tKSZiYjINXV2IFRWVoaYmBin7F/yU6/t7e2YOnUqjh07hq1bt2L48OGdzps8eTL27NmDqqoq09iBAwdQWVmJKVOmSB2LiIjonjjkgQO7d+/GxIkTUVtbi02bNpktnzZtGgBg0aJF2Lp1K8aMGYN58+ahoaEBy5cvR0REBKZPny51LCIionsieaP85ptvoFAosHv3buzevdtsmUKhMDXK4OBgHD58GPPnz0d2djZUKhWSkpKQl5fHx9cRkQ1uXxlfXl4ucw5zQ4YMQe/eveWOQRKSvFEePHjQ5rlDhw7Fvn37pI5ARPeFswC+P0vlKk6cOMFrKH5kHHbVKxGRc2wCECZ3CADlAFyraZM02CiJqIcLA8AjOHIcvriZiIjICjZKIiIiK3jqlYhIEq53FS6vwJUGG6UT3bp1y+ydnHJypT/MRD8OZwG41lW4vAJXGmyUTlRRUeG0Ry4RkVxc4SpcXoErJTZKWbjCH6SPAbwpcwaiHyNehftjw0YpC1f4g8RTr0REtpD1qtfm5ma89tprCAoKQu/evfHoo49i//79ckYiIiIyI2ujfOmll5Cfn4///u//xp/+9Ce4ubnh5z//Ob788ks5Y0nmh+9Pc209KSvQs/Iyq2N8JXcAO/Scuvasv7ecQ7ZGWVJSgi1btuDdd9/Fe++9h7S0NHz++ed48MEHsXDhQrliSapn/YbrSVmBnpWXWR3ja7kD2KHn1LVn/b3lHLI1ym3btqFXr16YNWuWaUylUmHmzJn4+uuvUV1dLVc0IiIiE9ku5tHpdAgNDUWfPn3MxuPi4gAAJ0+exIABA7q1j7///e/46iv5Ts98++23eO2110z/funSJdmyENH95N4fflBfX4+ysjKpAwHouQ9AkK1R6vV6aLVai3Hj2MWLF7u9j/379+P999eiV68Hu72te9HWdgn5+TuN/4bW1jOy5CCi+81ZAPf+8ANH3e/dUx+AIFujNBgMUKlUFuOenp6m5Z2tA9j+f0lXrlyBUhmA1tYF3UjaHe+jtTXj/399HYDx6PJjyH97hvGCKWOWKgCb5YtjkeduHJnX3ix3052sUme5m7tldXYea+r+/5+ukAWwXhtn//kyZpkJwPKAxLotAKZKGwd6AOslfSKYcVud9QrJCZkMGzZM/OxnP7MY/9e//iUUCoVYu3atxbJNmzYJAPzhD3/4wx/+CABi06ZNDu9Xsh1RarXaTk+v6vV6AEBQUJDFsqeffhqbNm1CSEgI1Gq1wzMSEZFrMhgMOHv2LJ5++mmH70u2RhkVFYVDhw6hoaEB3t7epvFjx44BACIjIy3W8fPzw4svvui0jERE5Loee+wxp+xHtttDfvnLX6K9vR1r1641jTU3N+PDDz/Eo48+2u0rXomIiKQg2xFlfHw8pkyZgtdffx1XrlzBoEGDsHHjRpw/fx4ffvihXLGIiIjMKIQQQq6dNzc3480338SmTZtw7do1PPLII8jJycHYsWPlikRERGRG1kZJRETk6mR9KDoREZGrc3qj7M6rtfbv34+EhAQEBATA29sbjzzyCFauXImOjg6LuV999RVGjRoFLy8vaLVazJs3D42NjS6X9cknn4RSqbT4GT9+vNOyfvbZZ6Za+fr6YsqUKTh37lync+Wuq61ZpaprY2MjlixZgsTERPj6+kKpVGLjxo02r3/9+nXMmjUL/v7+6NOnD5566inodLpO53a3ts7KKkVtu5P10qVLyM7OxpgxY+Dt7Q2lUonDhw93OV/OutqTVe66HjhwADNmzEBoaCi8vLwwaNAgpKend/noTTnrak9WSf4ucPidmj/w/PPPC3d3d7Fw4UKxbt06MXLkSOHu7i6OHj1qdb1PPvlEKBQKER4eLgoKCsTatWvFpEmThEKhEPPmzTObq9PphKenp4iJiRFr1qwRixcvFp6enmL8+PEul/WJJ54QAwcOFJs3bzb7OXjwoFOy7t69WyiVShEfHy9Wrlwpfv/73wt/f38RHBwsampqzObKXVd7skpV1zNnzgiFQiFCQkLEmDFjhEKhEBs3brRp3fb2djFy5EjRp08f8fbbb4v3339fDBs2TPTt21dUVlaazZWits7KKkVtu5P14MGDQqFQiIcffliMHDlSKBQKcfjw4U7nyl1Xe7LKXdeYmBgxaNAgkZ2dLdavXy8WLVok+vbtKwIDA8WlS5fM5spdV3uySlFXpzbKY8eOCYVCIfLy8kxjTU1NYvDgwWLkyJFW133hhReEp6enuHbtmtn4E088IXx8fMzGxo8fLwYMGCAaGhpMYx988IFQKBTiH//4h0tlfeKJJ0R4eLhNmRyRdejQoSI0NFS0traaxr755hvh5uYmXn31VbO5ctfVnqxS1FUIIZqbm8Xly5eFEEIcP37crj/MW7ZsEQqFQmzfvt00VlNTI/r16yfnUErKAAAIPklEQVReeOEFs7lS1NZZWaWobXeyNjQ0mP5sbd261Wrzkbuu9mSVu65HjhyxGPviiy+EQqEQixcvNhuXu672ZJWirk499dqdV2up1WqoVCr4+PiYjQcGBpo9jf7GjRvYv38/pk2bZvZmktTUVPTp0wfFxcUuk9VICIH29nbcvHnTpmxSZa2rq0N5eTmee+459Or1/Z1CERERGDJkCP7617+axuSuqz1ZjbpbVwDw8PBAQECAaXv22LZtGwIDA/GLX/zCNObn54fk5GT8/e9/R2trKwDpauuMrEbdrW13svbp0wcPPPDAXee5Ql1tzWokZ11HjRplMfb444/D19cXFRUVpjFXqKutWY26W1enNkpbXq3Vlblz56KjowOzZ89GRUUFzp07h9WrV2Pnzp14/fXXTfNOnTqFtrY2xMbGmq3v7u6OyMjILr8fkiOr0enTp+Hl5YW+fftCq9XirbfeQltbm005u5O1ubkZADp9HGDv3r2h1+tx5coVAPLX1Zasly9fNhvvbl27S6fTdfqmhLi4ONy6dQunT58GIF1tnZHVSO7a2sIV6movV6vrzZs30dDQAD8/P9OYq9a1s6xG3a2rUx840J1Xaz3yyCP4/PPPMXHiRHzwwQcAADc3N7z//vtmRyfGZ8V2tp/AwEAcPXrUZbICwODBg5GQkIDw8HA0NjZi69at+P3vf4/Tp093epQkZdb+/fvjgQcesKjJ1atX8e233wIAqqurERAQIHtdbcl68eJF9O/fH4A0de0uvV6PJ5980mL8zs86bNgwyWrbHbZmBVyjtrZwhbrawxXrWlBQgNbWVkyd+v3bRFy1rp1lBaSpq1Mb5b28WsuooqICEyZMwIMPPojly5fD09MTf/nLX/DKK6+gf//+ePbZZ8220dV+bH0lizOyAjA1UqMXX3wRs2fPxrp165CVlYXhw4c7LKtSqcTs2bPx3nvvYdGiRZg+fTpu3LiBhQsXorW1FUII07py19WerIA0de2upqYmmz6rVLXtDluzAq5RW1u4Ql3t4Wp1/eKLL/C73/0OU6dONfufKFesa1dZAWnq6tRTr2q12nQK7U5NTU2m5V1ZsGABevXqhUOHDmHatGn45S9/iR07dmDUqFHIyMgw3XZh3EZX+7H17dqOzNre3m5136+++iqA25dAOzrr22+/jZkzZ2LZsmV4+OGHERcXBw8PD8ycORMATKdIXaGutmbtir117S5bP6tUte2O7vx3AZxfW1u4Ql27S666VlRU4LnnnkNERIRFo3G1ulrL2hV76+rURnkvr9YyOnr0KJ566imL/wgTJ07ExYsXcfbsWdM+7tzmD/djbR/OytrVPYpGwcHBAG5fwOLorO7u7li3bh0uXryII0eO4PTp0/jkk09w/fp1uLm5YfDgwaZ93LnNH+7HGXW1NWtX7K1rd9n6WaWqbXd0578L4Pza2sIV6tpdctT1woULGDduHPr164ePP/4YXl5eZstdqa53y9oVe+vq1EYZFRWF06dPo6GhwWzc2qu1jNra2jo9EjNejWf8YvanP/0pevXqhdLSUrN5LS0tOHnypNV9ODtrV/7zn/8AAPz9/R2e1SggIACPPfYYBg8ejPb2dhw6dAjDhw83NXtXqKutWbtib127KzIyEmVlZRZX9B07dgxeXl4IDQ0FIF1tnZG1K86urS1coa7d5ey6Xr16FePGjUNrays+/fRT03f+d3KVutqStSt217VbN5fYyXgP3YoVK0xjxnvoRowYYRrT6/WivLzc7H65UaNGCY1GI65evWoaa2trEzExMcLHx0e0tbWZxsePHy+CgoI6vcfn008/dZmsN27cEE1NTWb77ejoEFOnThVKpVLodDqHZ+3Mu+++KxQKhdixY4fZuNx1tTWrVHX9odLS0i7v9eosr/HexG3btpnGampqxAMPPCBSUlLM1peits7I6oja2pv1TrbcRylnXW3N6gp1vXnzpoiPjxc+Pj6irKzM6rblrqutWaWqq9OfzJOcnGx6KsuaNWvEyJEjhYeHh9kNpL/61a+EQqEQ586dM4198sknQqlUisGDB4tly5aJP/3pT2LEiBFCoVCI3Nxcs32UlZUJT09PER0dLVatWiXeeOMNoVarRWJioktlPXjwoAgMDBTz588X77//vlixYoV47LHHhEKhEL/+9a+dkvWjjz4SkyZNEvn5+WLt2rUiOTlZKBQKMWvWLIt9yF1XW7NKWVchhFi5cqXIyckRc+bMEQqFQkyePFnk5OSInJwcUV9f32Xe9vZ2MWLECOHt7W32tBsfHx9x+vRps31IVVtHZ5WytveaVQhhmvf8888LhUIhZs6caRpzpbramtUV6vrss8+a8n300UdmP3/7299cqq62ZpWqrk5vlE1NTeK3v/2t0Gq1wtPTUwwfPtziSQ4vvfSSUCqVFr/h9u3bJx5//HHh5eUlVCqVeOSRR8TatWs73c/Ro0fFY489JtRqtejfv7+YO3euuHnzpktlPXPmjEhOThYPPfSQUKvVwsvLS8TFxXX5mRyRtaSkRDzxxBPC19dXqNVqERUVZXX/ctbV1qxS1lUIIUJCQoRCoRAKhUIolUqhVCpNvzbm6+r3wbVr10RaWprw8/MTXl5eYsyYMeLEiROd7keK2jo6q5S17U5W47w7/2n89Q/JXVdbsrpCXUNCQszy3fnz0EMPWexHzrramlWquvI1W0RERFbwNVtERERWsFESERFZwUZJRERkBRslERGRFWyUREREVrBREhERWcFGSUREZAUbJRERkRVslERERFawURIREVnBRklERGQFGyUREZEV/wep8Wkcr6uWHQAAAABJRU5ErkJggg==",
+      "text/plain": [
+       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x7f28c0d89590>)"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "figure(figsize=(5, 2))\n",
+    "plt[:hist](Float64[d[1] for d in solutions]);"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "What is left is to code own MCMC sampler and start doing Bayesian FEM."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this notebook the basic instructions how to develop JuliaFEM has been given. The most imporant concepts has been considered; how to develop own element with own basis, several ways how to define own equation, and how to finally assemble and calculate the problem using solver. Any comments and/or discussion about technical details, theory, programming etc. is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues"
    ]
   }
  ],
diff --git a/src/basis.jl b/src/basis.jl
index 5343198..62cacec 100644
--- a/src/basis.jl
+++ b/src/basis.jl
@@ -45,12 +45,14 @@ function Base.call(field::DiscreteField, time::Number,
 
     # special cases, only 1 timestep defined or time = -Inf -> return first ts
     if (length(field) == 1) || (time == -Inf)
-        return field[1][end]
+        #return field[1][end]
+        return first(field)
     end
 
     # special case, time = +Inf -> return last ts
     if time == +Inf
-        return field[end][end]
+        #return field[end][end]
+        return last(field)
     end
 
     # very likely we are always near some defined timestep, usually field
diff --git a/src/dirichlet.jl b/src/dirichlet.jl
index c63360d..ba9ae5a 100644
--- a/src/dirichlet.jl
+++ b/src/dirichlet.jl
@@ -35,7 +35,7 @@ function DBC2D2(element::Seg2)
         IntegrationPoint([-sqrt(1/3)], 1.0),
         IntegrationPoint([+sqrt(1/3)], 1.0)]
     if !haskey(element, "reaction force")
-        element["reaction force"] = FieldSet()
+        element["reaction force"] = zeros(1, 2)
     end
     DBC2D2(element, integration_points)
 end
diff --git a/src/elements.jl b/src/elements.jl
index 737600a..2acd923 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -42,9 +42,9 @@ function test_element(element_type)
     end
 
     # try to interpolate some scalar field
-    element["field1"] = Field(0.0, collect(1:n))
+    element["field1"] = Field(collect(1:n))
     # TODO: how to parametrize this?
-    element["geometry"] = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    element["geometry"] = Field(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
     basis = get_basis(element)
@@ -55,7 +55,7 @@ function test_element(element_type)
     val2 = basis("field1", mid, 0.0)
     info("field val at $mid: $val2")
     val3 = dbasis(mid, 0.0)
-    info("derivative of basis at $mid: $val3")
+    info("derivative of basis at $mid:\n$val3")
     val4 = dbasis("field1", mid, 0.0)
     info("field val at $mid: $val4")
 
@@ -64,135 +64,130 @@ end
 
 """ Get FieldSet from element. """
 function Base.getindex(element::Element, field_name)
-    element.fields[field_name]
+    return element.fields[field_name]
 end
 
-"""Add new FieldSet to element.
+"""Add new Field to element.
 
 Examples
 --------
->>> element["geometry"] = [1, 2, 3, 4]
-JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
+>>> element["temperature"] = [1, 2, 3, 4]
+>>> element["temperature"] = (0.0, [0, 0, 0, 0]), (1.0, [1, 2, 3, 4])
+>>> element["temperature"] = (0.0 => [0, 0, 0, 0], 1.0 => [1, 2, 3, 4])
 """
 function Base.setindex!(element::Element, field_data, field_name)
-    #element.fields[field_name] = field_data
     setindex!(element.fields, field_data, field_name)
 end
 
+function Base.setindex!(element::Element, field_data::Tuple, field_name)
+    field = Field()
+    for (time, data) in field_data
+        ts = TimeStep(time, Increment[Increment(data)])
+        push!(field, ts)
+    end
+    element[field_name] = field
+end
+
 function get_connectivity(el::Element)
-    el.connectivity
+    return el.connectivity
 end
 
 abstract AbstractFunctionSpace
 
 type FunctionSpace <: AbstractFunctionSpace
-    element :: Element
+    basis :: Basis
+    fields :: FieldSet
 end
 
 type GradientFunctionSpace <: AbstractFunctionSpace
-    element :: Element
-end
-
-type MixedFunctionSpace <: AbstractFunctionSpace
-    element1 :: Element
-    element2 :: Element
+    basis :: Basis
+    fields :: FieldSet
 end
 
 function get_basis(element::Element)
-    return FunctionSpace(element)
+    return FunctionSpace(element.basis, element.fields)
 end
 
 function get_dbasis(element::Element)
-    return GradientFunctionSpace(element)
+    return GradientFunctionSpace(element.basis, element.fields)
 end
 
 function grad(u::FunctionSpace)
-    return 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.data[1]
-    end
-    h = u.element.basis.basis(xi)
-    #@debug("vec(h) = $(vec(h)), size(h) = $(size(vec(h)))")
-    #@debug("f = $f, size(f) = $(size(f))")
-    #return dot(vec(h), f)
-    return sum(vec(h).*f)
+    return GradientFunctionSpace(u.basis, u.fields)
 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)
-    dN = gradu.element.basis.dbasisdxi(xi)
-    J = sum([dN[:,i]*X[i]' for i=1:length(X)])
-    grad = inv(J)*dN
-    gradf = sum([grad[:,i]*f[i]' for i=1:length(f)])'
-    return gradf
+function call(u::FunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
+    return u.basis(xi)
 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)
-    dN = gradu.element.basis.dbasisdxi(xi)
-    J = sum([dN[:,i]*X[i]' for i=1:length(X)])
-    grad = inv(J)*dN
-    return grad
+function call(gradu::GradientFunctionSpace, xi::Union{Vector, IntegrationPoint}, t::Number=0.0)
+    geometry = gradu.fields["geometry"](t)
+    gradu.basis(geometry, xi, Val{:grad})
 end
 
+""" Evaluate field on element function space. """
+function call(u::FunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
+    field = !isa(variation, Void) ? variation : u.fields[field_name](t)
+    if length(field) == 1
+        return field.data[1]
+    end
+    u.basis(field, xi)
+end
+
+""" Evaluate gradient of field on element function space. """
+function call(gradu::GradientFunctionSpace, field_name, xi::Union{Vector, IntegrationPoint}, t::Number=0.0, variation=nothing)
+    field = !isa(variation, Void) ? variation : gradu.fields[field_name](t)
+    geometry = gradu.fields["geometry"](t)
+    gradu.basis(geometry, field, xi, Val{:grad})
+end
+
+
 # on-line functions to get api more easy to use, ip -> xi.ip
-call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
-call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
+#call(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
+#call(u::GradientFunctionSpace, ip::IntegrationPoint, t::Number=Inf) = call(u, ip.xi, t)
 # i think these will be the most called functions.
-call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
-call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=Inf, variation=nothing) = call(u, field_name, ip.xi, t, variation)
-call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
-call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
+#call(u::FunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+#call(u::GradientFunctionSpace, field_name, ip::IntegrationPoint, t::Number=0.0, variation=nothing) = call(u, field_name, ip.xi, t, variation)
+#call(u::FunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
+#call(u::GradientFunctionSpace, field_name) = (args...) -> call(u, field_name, args...)
 
 """ Return a field from function space. """
-function get_field(u::FunctionSpace, field_name, time=Inf)
-    return u.element[field_name](time)
+function get_field(u::FunctionSpace, field_name, time::Number=0.0)
+    return u.fields[field_name](time)
 end
 
 """ Return a 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)
+function get_field(u::FunctionSpace, field_name, time::Number=0.0, variation=nothing)
+    return !isa(variation, Void) ? variation : u.fields[field_name](time)
 end
 
-""" Return a fieldset from function space. """
+""" Return a field from function space. """
 function get_fieldset(u::FunctionSpace, field_name)
-    return u.element[field_name]
+    return u.fields[field_name]
 end
 
 """ Get a determinant of element in point ξ. """
-function LinAlg.det(u::FunctionSpace, xi::Vector, t::Number=Inf)
-    X = u.element["geometry"](t)
-    dN = u.element.basis.dbasisdxi(xi)
+function LinAlg.det(u::FunctionSpace, xi::Vector, time::Number=0.0)
+    X = u.fields["geometry"](time)
+    dN = u.basis.dbasisdxi(xi)
     J = sum([dN[:,i]*X[i]' for i=1:length(X)])
     m, n = size(J)
     return m == n ? det(J) : norm(J)
 end
-function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, t::Number=Inf)
-    LinAlg.det(u, ip.xi, t)
+function LinAlg.det(u::FunctionSpace, ip::IntegrationPoint, time::Number=0.0)
+    LinAlg.det(u, ip.xi, time)
 end
 function LinAlg.det(u::FunctionSpace)
     return (args...) -> det(u, args...)
 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...)
 
-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. """
+""" Check does field exist. """
 function Base.haskey(element::Element, what)
     haskey(element.fields, what)
 end
diff --git a/src/equations.jl b/src/equations.jl
index d53d8fa..c62bf02 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -57,10 +57,10 @@ function initialize_local_assembly!(assembly::LocalAssembly, equation::Equation)
 end
 
 has_mass_matrix(equation::Equation) = false
-function get_mass_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+function get_mass_matrix(equation::Equation, ip, time=0.0, problem=nothing)
     get_mass_matrix(equation, ip, time)
 end
-function get_mass_matrix(equation::Equation, ip, time=Inf)
+function get_mass_matrix(equation::Equation, ip, time=0.0)
     get_mass_matrix(equation, ip)
 end
 function get_mass_matrix(equation::Equation, ip)
@@ -68,10 +68,10 @@ function get_mass_matrix(equation::Equation, ip)
 end
 
 has_stiffness_matrix(equation::Equation) = false
-function get_stiffness_matrix(equation::Equation, ip, time=Inf, problem=nothing)
+function get_stiffness_matrix(equation::Equation, ip, time=0.0, problem=nothing)
     get_stiffness_matrix(equation, ip, time)
 end
-function get_stiffness_matrix(equation::Equation, ip, time=Inf)
+function get_stiffness_matrix(equation::Equation, ip, time=0.0)
     get_stiffness_matrix(equation, ip)
 end
 function get_stiffness_matrix(equation::Equation, ip)
@@ -79,10 +79,10 @@ function get_stiffness_matrix(equation::Equation, ip)
 end
 
 has_force_vector(equation::Equation) = false
-function get_force_vector(equation::Equation, ip, time=Inf, problem=nothing)
+function get_force_vector(equation::Equation, ip, time=0.0, problem=nothing)
     get_force_vector(equation, ip, time)
 end
-function get_force_vector(equation::Equation, ip, time=Inf)
+function get_force_vector(equation::Equation, ip, time=0.0)
     get_force_vector(equation, ip)
 end
 function get_force_vector(equation::Equation, ip)
@@ -90,10 +90,10 @@ function get_force_vector(equation::Equation, ip)
 end
 
 has_residual_vector(equation::Equation) = false
-function get_residual_vector(equation::Equation, ip, time=Inf, problem=nothing)
+function get_residual_vector(equation::Equation, ip, time=0.0, problem=nothing)
     get_residual_vector(equation, ip, time)
 end
-function get_residual_vector(equation::Equation, ip, time=Inf)
+function get_residual_vector(equation::Equation, ip, time=0.0)
     get_residual_vector(equation, ip)
 end
 function get_residual_vector(equation::Equation, ip)
@@ -101,10 +101,10 @@ function get_residual_vector(equation::Equation, ip)
 end
 
 has_potential_energy(equation::Equation) = false
-function get_potential_energy(equation::Equation, ip, time=Inf, problem=nothing)
+function get_potential_energy(equation::Equation, ip, time=0.0, problem=nothing)
     get_potential_energy(equation, ip, time)
 end
-function get_potential_energy(equation::Equation, ip, time=Inf)
+function get_potential_energy(equation::Equation, ip, time=0.0)
     get_potential_energy(equation, ip)
 end
 function get_potential_energy(equation::Equation, ip)
@@ -117,7 +117,7 @@ get_integration_points(equation::Equation) = equation.integration_points
 
 """ Return a local assembly for element. """
 function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
-                                   unknown_field_name::ASCIIString, time::Number=Inf,
+                                   unknown_field_name::ASCIIString, time::Number=0.0,
                                    problem=nothing)
 
     initialize_local_assembly!(assembly, equation) # zero all
@@ -174,6 +174,7 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
         assembly.stiffness_matrix += hessian
         assembly.force_vector -= ForwardDiff.gradient(allresults)  # <--- minus explained in tutorial
         assembly.potential_energy = ForwardDiff.value(allresults)
+        #info("potential energy of system: $(assembly.potential_energy)")
     end
 
     # 3. virtual work form - user has defined residual vector δW_int(u,δu) + δW_ext(u,δu) = 0 ∀ v
diff --git a/src/fields.jl b/src/fields.jl
index b64f81e..bffdd4f 100644
--- a/src/fields.jl
+++ b/src/fields.jl
@@ -214,6 +214,15 @@ function Base.push!(field::DefaultDiscreteField, timestep::TimeStep)
     push!(field.timesteps, timestep)
 end
 
+function Base.push!(field::DefaultDiscreteField, data::Union{Vector, Matrix})
+    push!(field[end], Increment(data))
+end
+
+function Base.push!(field::DefaultDiscreteField, data::Pair)
+    ts = TimeStep(data[1], Increment(data[2]))
+    push!(field, ts)
+end
+
 """Quickly create fields.
 
 Examples
@@ -287,3 +296,7 @@ function Base.convert(::Type{ContinuousField}, data::Function)
     return convert(DefaultContinuousField, data)
 end
 
+function Base.length(::Field)
+    return 1
+end
+
diff --git a/src/integrate.jl b/src/integrate.jl
index 5b8497f..f5de31e 100644
--- a/src/integrate.jl
+++ b/src/integrate.jl
@@ -16,3 +16,48 @@ function get_default_integration_points(element::Seg2)
         IntegrationPoint([0.0], 2.0)
     ]
 end
+
+function line3()
+    [
+        IntegrationPoint([0.0], 8/9),
+        IntegrationPoint([-sqrt(3/5)], 5/9),
+        IntegrationPoint([+sqrt(3/5)], 5/9)
+    ]
+end
+
+function line5()
+    [
+        IntegrationPoint([-1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900),
+        IntegrationPoint([-1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
+        IntegrationPoint([0.0], 128/225),
+        IntegrationPoint([ 1/3*sqrt(5 - 2*sqrt(10/7))], (322+13*sqrt(70))/900),
+        IntegrationPoint([ 1/3*sqrt(5 + 2*sqrt(10/7))], (322-13*sqrt(70))/900)
+    ]
+end
+
+    #integration_points = [
+    #    IntegrationPoint([ 0.0000000000000000], 0.5688888888888889),
+    #    IntegrationPoint([-0.5384693101056831], 0.4786286704993665),
+    #    IntegrationPoint([ 0.5384693101056831], 0.4786286704993665),
+    #    IntegrationPoint([-0.9061798459386640], 0.2369268850561891),
+    #    IntegrationPoint([ 0.9061798459386640], 0.2369268850561891)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([+sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+    #    IntegrationPoint([-sqrt(3/7 - 2/7*sqrt(6/5))], (18+sqrt(30))/36)
+    #    IntegrationPoint([+sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+    #    IntegrationPoint([-sqrt(3/7 + 2/7*sqrt(6/5))], (18-sqrt(30))/36)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([0.0], 8/9),
+    #    IntegrationPoint([-sqrt(3/5)], 5/9),
+    #    IntegrationPoint([+sqrt(3/5)], 5/9)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([-sqrt(1/3)], 1)
+    #    IntegrationPoint([+sqrt(1/3)], 1)
+    #]
+    #integration_points = [
+    #    IntegrationPoint([0.0], 2)
+    #]
+
diff --git a/src/solvers.jl b/src/solvers.jl
index e7d98d6..580e90b 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -10,10 +10,10 @@ Solve field equations for single element with some dofs fixed. This can be used
 to test nonlinear element formulations.
 """
 function solve!(equation::Equation, unknown_field_name::ASCIIString,
-                 free_dofs::Array{Int, 1}, time::Number=Inf; 
+                 free_dofs::Array{Int, 1}, time::Number=0.0; 
                  max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
     element = get_element(equation)
-    x0 = element[unknown_field_name](-Inf)
+    x0 = element[unknown_field_name](0.0)
     x = zeros(prod(size(equation)))
     dx = fill!(similar(x), 0.0)
     la = initialize_local_assembly()
@@ -27,15 +27,10 @@ function solve!(equation::Equation, unknown_field_name::ASCIIString,
         end
         dx[free_dofs] = A \ b
         x += dx
-        new_field = similar(x0, x)
-        new_field.time = time
-        new_field.increment = i
-        push!(element[unknown_field_name], new_field)
-        if norm(dx) < tolerance
-            return
-        end
+        push!(element[unknown_field_name], reshape(x, size(equation)))
+        norm(dx) < tolerance && return
     end
-    Logging.err("Did not converge in $max_iterations iterations")
+    error("Did not converge in $max_iterations iterations")
 end
 
 """
@@ -44,15 +39,18 @@ to test nonlinear element formulations. Dirichlet boundary is assumed to be homo
 and degrees of freedom are eliminated. So if boundary condition is known in nodal
 points and everything is zero this should be quite good.
 """
-function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf;
+function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=1.0;
                 max_iterations::Int=10, tolerance::Float64=1.0e-12, dump_matrices::Bool=false)
+    info("start solver")
     ga = initialize_global_assembly(problem)
     x = zeros(ga.ndofs)
     dx = fill!(similar(x), 0.0)
     field_name = get_unknown_field_name(problem)
     dim = get_unknown_field_dimension(problem)
     for i=1:max_iterations
+        info("calculate global assembly")
         calculate_global_assembly!(ga, problem)
+        info("done")
         A = ga.stiffness_matrix[free_dofs, free_dofs]
         b = ga.force_vector[free_dofs]
         if dump_matrices
@@ -60,20 +58,17 @@ function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf;
             dump(full(b)')
         end
         dx[free_dofs] = lufact(A) \ full(b)
+        info("Difference in solution norm: $(norm(dx))")
         x += dx
         for equation in get_equations(problem)
             element = get_element(equation)
-            conn = get_connectivity(element)
-            gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
-            old_field = element[field_name](Inf)
-            new_field = similar(old_field, full(x[gdofs]))
-            push!(element[field_name][end], new_field)
-        end
-        if norm(dx) < tolerance
-            return
+            gdofs = get_gdofs(problem, equation)
+            data = reshape(full(x[gdofs]), size(equation))
+            push!(element[field_name], data)
         end
+        norm(dx) < tolerance && return
     end
-    Logging.err("Did not converge in $max_iterations iterations")
+    error("Did not converge in $max_iterations iterations")
 end
 
 """ Add new problem to solver. """
diff --git a/src/types.jl b/src/types.jl
index 7e3dea4..9fa6c78 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -31,3 +31,19 @@ end
 function Base.convert(::Type{Number}, ip::IntegrationPoint)
     return ip.xi
 end
+
+function Base.call(basis::Basis, ip::IntegrationPoint)
+    return basis(ip.xi)
+end
+
+function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint)
+    return call(basis, increment, ip.xi)
+end
+
+function Base.call(basis::Basis, increment::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+    return call(basis, increment, ip.xi, Val{:grad})
+end
+
+function Base.call(basis::Basis, geometry::Increment, field::Increment, ip::IntegrationPoint, ::Type{Val{:grad}})
+    return call(basis, geometry, field, ip.xi, Val{:grad})
+end
diff --git a/test/test_dirichlet.jl b/test/test_dirichlet.jl
new file mode 100644
index 0000000..146f20d
--- /dev/null
+++ b/test/test_dirichlet.jl
@@ -0,0 +1,18 @@
+
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module TestAutoDiffWeakForm
+
+using JuliaFEM.Test
+using JuliaFEM
+using JuliaFEM: Seg2, DirichletProblem
+
+function test_dirichlet_problem()
+    element = Seg2([3, 4])
+    element["geometry"] = Vector[[1.0, 1.0], [0.0, 1.0]]
+    problem = DirichletProblem(1)
+    push!(problem, element)
+end
+
+end
diff --git a/test/test_elasticity.jl b/test/test_elasticity.jl
index c33910f..0704c98 100644
--- a/test/test_elasticity.jl
+++ b/test/test_elasticity.jl
@@ -4,7 +4,9 @@
 module ElasticityTests
 
 using JuliaFEM.Test
-using JuliaFEM: Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem
+using JuliaFEM: Quad4, Field, FieldSet, CPS4,
+                get_basis, solve!,
+                PlaneStressElasticityProblem
 
 
 function test_elasticity_one_element()
@@ -21,7 +23,7 @@ function test_elasticity_one_element()
     disp = get_basis(element)("displacement", [1.0, 1.0])[2]
     info("displacement at tip: $disp")
     # verified using Code Aster.
-    @test disp ≈ -8.77303119819776
+    @test isapprox(disp, -8.77303119819776)
 end
 
 
diff --git a/test/test_elements.jl b/test/test_elements.jl
index 06ad936..8e465d1 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -19,7 +19,11 @@ end
 
 function MockElement(connectivity)
 
-    h(xi) = 1/4*[(1-xi[1])*(1-xi[2])   (1+xi[1])*(1-xi[2])   (1+xi[1])*(1+xi[2])   (1-xi[1])*(1+xi[2])]
+    h(xi) = 1/4*[
+        (1-xi[1])*(1-xi[2])
+        (1+xi[1])*(1-xi[2])
+        (1+xi[1])*(1+xi[2])
+        (1-xi[1])*(1+xi[2])]'
 
     dh(xi) = 1/4*[
         -(1-xi[2])    (1-xi[2])   (1+xi[2])  -(1+xi[2])
@@ -39,9 +43,47 @@ end
 """ test adding fieldsets and fields to element"""
 function test_add_fields_to_element()
     el = MockElement([1, 2, 3, 4])
-    el["geometry"] = [0.0, 0.0, 0.0, 0.0], [1.0, 1.0, 1.0, 1.0]
+    #geometry = Field([0.0, 0.0, 0.0, 0.0])
+    el["geometry"] = Field([0.0, 0.0, 0.0, 0.0])
+    @test el["geometry"][1].time == 0.0
+    @test last(el["geometry"]) == [0.0, 0.0, 0.0, 0.0]
+    el["geometry"] = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
+    @test last(el["geometry"])[3] == [1.0, 1.0]
+    el["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    @test last(el["geometry"])[3] == [1.0, 1.0]
+    el["geometry"] = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
+    @test last(el["geometry"])[3] == [1.0, 1.0]
+    el["geometry"] = (0.0, [0.0, 0.0, 0.0, 0.0]), (1.0, [1.0, 1.0, 1.0, 1.0])
     field = el["geometry"]
     @test length(field) == 2  # two time steps
+    el["boundary flux"] = (0.0, 0.0), (1.0, 6.0)
+end
+
+function test_add_fields_to_element_2()
+    el = MockElement([1, 2, 3, 4])
+    el["data"] = (0.0 => [1, 2], 1.0 => [2, 3])
+    @test length(el["data"]) == 2
+    @test el["data"][1].time == 0.0
+    @test el["data"][2].time == 1.0
+    @test last(el["data"][1]) == [1, 2]
+    @test last(el["data"][2]) == [2, 3]
+end
+
+function test_add_data_to_element_using_push()
+    el = MockElement([1, 2, 3, 4])
+    el["data"] = [0, 0, 0, 0]
+   
+    push!(el["data"], [1, 2, 3, 4])
+    @test length(el["data"]) == 1
+    @test length(el["data"][1]) == 2
+    @test el["data"][1].time == 0.0
+
+    push!(el["data"], 1.0 => [2, 3, 4, 5])  # creates new timestep at t=1.0
+    push!(el["data"], [3, 4, 5, 6])  # adds new increment data to last timestep
+    @test length(el["data"]) == 2
+    @test length(el["data"][2]) == 2
+    @test el["data"][2].time == 1.0
+    
 end
 
 #=
diff --git a/test/test_potential_energy.jl b/test/test_potential_energy.jl
new file mode 100644
index 0000000..68f45ae
--- /dev/null
+++ b/test/test_potential_energy.jl
@@ -0,0 +1,197 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module ElementTests
+
+using JuliaFEM.Test
+
+using JuliaFEM
+using JuliaFEM: Equation, Quad4, IntegrationPoint, initialize_local_assembly,
+                get_element, get_basis, grad, calculate_local_assembly!,
+                PlaneHeatProblem, Seg2, HeatEquation, Problem, solve!
+
+
+""" Diffusive heat transfer for 4-node bilinear element, with a nonlinear source term. """
+type DC2D4NL <: Equation
+    element :: Quad4
+    integration_points :: Array{IntegrationPoint, 1}
+end
+
+function DC2D4NL(element::Quad4, initial_temperature=zeros(4))
+    integration_points = [
+        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]
+    if !haskey(element, "temperature")
+        element["temperature"] = initial_temperature
+    end
+    DC2D4NL(element, integration_points)
+end
+
+function Base.size(equation::DC2D4NL)
+    return (1, 4)
+end
+
+""" Nonlinear flux term. """
+type DC2D2NL <: Equation
+    element :: Seg2
+    integration_points :: Array{IntegrationPoint, 1}
+end
+
+function DC2D2NL(element::Seg2, initial_temperature=zeros(2))
+    #integration_points = [
+    #    IntegrationPoint([0.0], 2.0)]
+    integration_points = JuliaFEM.line5()
+    if !haskey(element, "temperature")
+        element["temperature"] = initial_temperature
+    end
+    DC2D2NL(element, integration_points)
+end
+
+function Base.size(equation::DC2D2NL)
+    return (1, 2)
+end
+
+""" Calculate a potential Π = Wint - Wext of system. """
+function JuliaFEM.get_potential_energy(equation::DC2D4NL, ip, time; variation=nothing)
+    element = get_element(equation)
+    basis = get_basis(element)
+    k = basis("temperature thermal conductivity", ip, time)
+    f = basis("temperature load", ip, time)
+    T = basis("temperature", ip, time, variation)
+    c = basis("temperature nonlinearity coefficient", ip, time)
+    gradT = grad(basis)("temperature", ip, time, variation)
+    Wint = (k + c*T) * 1/2*vecdot(gradT, gradT)
+    #Wint = k*1/2*vecdot(gradT, gradT)
+    Wext = f*T
+    #Wext = 0.0
+    return Wint - Wext
+end
+
+function JuliaFEM.has_potential_energy(eq::DC2D4NL)
+    return true
+end
+
+function JuliaFEM.get_potential_energy(equation::DC2D2NL, ip, time; variation=nothing)
+    element = get_element(equation)
+    basis = get_basis(element)
+    T = basis("temperature", ip, time, variation)[1]
+    Wint = 0.0
+    sig = 5.7e-8
+    eps = basis("emissivity", ip, time)[1]
+    T_ext = basis("temperature external", ip, time)[1]
+    q0 = eps*sig*((T_ext+273.15)^4 - (T+273.15)^4)
+    Wext = q0*T
+    W = Wint - Wext
+    return W
+end
+
+function JuliaFEM.has_potential_energy(eq::DC2D2NL)
+    return true
+end
+
+function test_potential_energy_method()
+
+    # create model -- start
+    element = Quad4([1, 2, 3, 4])
+    element["geometry"] = Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
+    element["temperature thermal conductivity"] = 6.0
+    element["temperature load"] = [0.0, 0.0, 0.0, 0.0]
+    element["temperature nodal load"] = [3.0, 3.0, 0.0, 0.0]
+    element["temperature nonlinearity coefficient"] = [6.0, 6.0, 6.0, 6.0]
+    equation = DC2D4NL(element)
+    # create model -- end
+
+    la = initialize_local_assembly() # create workspace for local matrices
+    T = zeros(4) # create workspace for solution vector
+    dT = zeros(4) # 
+    fd = [1, 2] # free dofs
+    tic()
+    # start loops, in principle solve ∂r(u)/∂uΔu = -r(u) and update.
+    for i=1:10
+        calculate_local_assembly!(la, equation, "temperature")  # calculate local matrices
+        dT[fd] = la.stiffness_matrix[fd,fd] \ la.force_vector[fd]
+        T += dT
+        push!(element["temperature"], T) # add new increment to model
+        info("T = $T")
+        @printf("increment %2d, |du| = %8.5f\n", i, norm(dT))
+        err = last(element["temperature"])[1] - 2/3
+        isapprox(err, 0.0) && break
+    end
+    toc()
+    err = last(element["temperature"])[1] - 2/3
+    info("error: $err")
+    @test isapprox(err, 0.0)
+end
+
+type TestProblem <: Problem
+    unknown_field_name :: ASCIIString
+    unknown_field_dimension :: Int
+    equations :: Array{Equation, 1}
+    element_mapping :: Dict{DataType, DataType}
+end
+
+function TestProblem(equations=[])
+    element_mapping = Dict(
+        Quad4 => DC2D4NL,
+        Seg2 => DC2D2NL)
+    TestProblem("temperature", 1, equations, element_mapping)
+end
+
+function test_potential_energy_method_2()
+
+    # create model -- start
+    N = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
+    element = Quad4([1, 2, 3, 4])
+    element["geometry"] = Vector[N[1], N[2], N[3], N[4]]
+    element["temperature thermal conductivity"] = 6.0
+    element["temperature load"] = [0.0, 0.0, 0.0, 0.0]
+    element["temperature nonlinearity coefficient"] = [0.0, 0.0, 0.0, 0.0]
+    #equation1 = DC2D4NL(element, initial_temperature=ones(4))
+    equation1 = DC2D4NL(element)
+
+    boundary_element = Seg2([1, 2])
+    boundary_element["geometry"] = Vector[N[1], N[2]]
+    boundary_element["emissivity"] = 0.5
+    boundary_element["temperature external"] = 10.0
+    #equation2 = DC2D2NL(boundary_element, initial_temperature=ones(4))
+    equation2 = DC2D2NL(boundary_element)
+    # create model -- end
+
+    element["temperature"] = ones(4)
+    boundary_element["temperature"] = ones(2)
+
+    equations = [equation1, equation2]
+    la = initialize_local_assembly() # create workspace for local matrices
+    T = zeros(4) # create workspace for solution vector
+    dT = zeros(4) # 
+    fd = [1, 2] # free dofs
+    info("equation 1")
+    calculate_local_assembly!(la, equation1, "temperature")
+    info("stiffness matrix: $(la.stiffness_matrix)")
+#   info("force vector: $(la.force_vector)")
+    info("equation 2")
+    calculate_local_assembly!(la, equation2, "temperature")
+#   info("stiffness matrix: $(la.stiffness_matrix)")
+    info("force vector: $(la.force_vector)")
+
+    info("Creating problem")
+    #problem = PlaneHeatProblem("temperature", 1, equations, Dict())
+    problem = TestProblem(equations)
+
+    free_dofs = [1, 2]
+    tic()
+    solve!(problem, free_dofs; max_iterations=10)
+    toc()
+    temp = get_basis(boundary_element)("temperature", [0.0])[1]
+    info("temperature = $temp")
+    #err = last(element["temperature"])[1] - 2/3
+    #info("error: $err")
+    # 0.3888756709834147 tulee jostakin syysta...
+    # tai -0.39411350336960116
+    info(boundary_element["temperature"])
+    @test isapprox(temp, 2.93509690572300E+00)  # tested using Code Aster
+end
+
+end
diff --git a/test/test_randomfields.jl b/test/test_randomfields.jl
new file mode 100644
index 0000000..a93ab11
--- /dev/null
+++ b/test/test_randomfields.jl
@@ -0,0 +1,34 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module RandomFieldTests
+
+using JuliaFEM
+using JuliaFEM: DiscreteField, Field, Increment, Quad4
+using JuliaFEM.Test
+
+type RandomField <: DiscreteField
+    mu :: Float64
+    std :: Float64
+end
+
+Base.first(field::RandomField) = Increment(randn(2, 4).*field.std^2 + field.mu)
+
+
+function test_interpolate_in_time()
+    r = RandomField(10.0, 0.0)
+    f = Increment(ones(2, 4)*10.0)
+    @test r(0.0) == f
+    @test r(-Inf) == f
+    @test r(+Inf) == f
+    @test r(1.0) == f
+end
+
+function test_interpolate_in_spatial_domain()
+    basis = Quad4([1, 2, 3, 4]).basis
+    r = RandomField(10.0, 0.0)
+    feval = basis(r(0.0), [0.0, 0.0])
+    @test feval == [10.0, 10.0]
+end
+
+end
diff --git a/test/test_virtual_work.jl b/test/test_virtual_work.jl
new file mode 100644
index 0000000..5d1acef
--- /dev/null
+++ b/test/test_virtual_work.jl
@@ -0,0 +1,81 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+module TestAutoDiffWeakForm
+
+using JuliaFEM.Test
+using JuliaFEM
+using JuliaFEM: Quad4, Equation, IntegrationPoint,
+                solve!, get_field, get_element, get_basis,
+                grad
+
+""" Plane stress formulation for 4-node bilinear element. """
+type CPS4 <: Equation
+    element :: Quad4
+    integration_points :: Array{IntegrationPoint, 1}
+end
+
+function CPS4(element::Quad4)
+    integration_points = [
+        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]
+    if !haskey(element, "displacement")
+        # initial field must be defined if using autodiff
+        element["displacement"] = zeros(2, 4)
+    end
+    CPS4(element, integration_points)
+end
+
+JuliaFEM.size(eq::CPS4) = (2, 4)
+
+function JuliaFEM.get_residual_vector(equation::CPS4, ip, time; variation=nothing)
+    element = get_element(equation)
+    basis = get_basis(element)
+    dbasis = grad(basis)
+
+    # material parameters
+    E = basis("youngs modulus", ip, time)
+    nu = basis("poissons ratio", ip, time)
+    mu = E/(2*(1+nu))
+    la = E*nu/((1+nu)*(1-2*nu))
+    la = 2*la*mu/(la + 2*mu)  # <- correction for 2d
+
+    # elasticity formulation
+    u = basis("displacement", ip, time, variation)
+    gradu = dbasis("displacement", ip, time, variation)
+    F = I + gradu
+    b = basis("displacement volume load", ip, time)
+    E = 1/2*(F'*F - I)
+    S = la*trace(E)*I + 2*mu*E
+    P = F*S
+
+    # residual vector
+    r_int = P*dbasis(ip,time)
+    r_ext = b*basis(ip,time)
+    r = r_int - r_ext
+    return vec(r)
+end
+
+JuliaFEM.has_residual_vector(equation::CPS4) = true
+
+function test_residual_form()
+    # create model -- start
+    element = Quad4([1, 2, 3, 4])
+    element["geometry"] = Vector[[0.0,0.0], [10.0,0.0], [10.0,1.0], [0.0,1.0]]
+    element["youngs modulus"] = 500.0
+    element["poissons ratio"] = 0.3
+    element["displacement volume load"] = Vector[[0.0,-10.0], [0.0,-10.0], [0.0,-10.0], [0.0,-10.0]]
+    equation = CPS4(element)
+    # create model -- end
+
+    free_dofs = [3, 4, 5, 6]
+    solve!(equation, "displacement", free_dofs)  # launch a newton solver for single element
+    disp = get_basis(element)("displacement", [1.0, 1.0])[2]
+    println("displacement at tip: $disp")
+    # verified using Code Aster.
+    @test isapprox(disp, -8.77303119819776E+00)
+end
+
+end