diff --git a/CONTRIBUTING.rst b/CONTRIBUTING.rst
index 06b45de..54c212b 100644
--- a/CONTRIBUTING.rst
+++ b/CONTRIBUTING.rst
@@ -2,28 +2,17 @@
 Contributing to JuliaFEM
 ========================
 
-This very important document need to be done.
-
-For now, read 
-
-https://github.com/JuliaLang/julia/blob/master/CONTRIBUTING.md
-
-
-How to contribute
------------------
 Here are the basic steps for contributing to JuliaFEM:
 
 1) Create an account or sign in to `GitHub <https://github.com/>`_.
 
 2) Go to `Git home page <http://git-scm.com/>`_ and download the Git installer. Run the installer to get Git on your computer. It is a version control system used by GitHub. To learn its basics, go through this `Git tutorial <https://try.github.io/levels/1/challenges/1>`_.
 
-3) Install Julia (v0.4+) to your computer. At `Julia readme
-<https://github.com/JuliaLang/julia/blob/master/README.md>`_ you'll find complete instructions for installing it for your platform.
+3) Install Julia (v0.4+) to your computer. At `Julia readme <https://github.com/JuliaLang/julia/blob/master/README.md>`_ you'll find complete instructions for installing it for your platform.
 
 4) Go to the `JuliaFEM GitHub page <https://github.com/JuliaFEM/JuliaFEM.jl>`_. At the top-right corner, press the ``Fork``-button to fork your own copy of JuliaFEM to your repository.
 
-5) Clone JuliaFEM from your repository to your computer. Navigate to the folder you want to clone it to, and type the following command (inserting your GitHub username to its place):
-``git clone https://github.com/your_github_username/JuliaFEM.jl.git``
+5) Clone JuliaFEM from your repository to your computer. Navigate to the folder you want to clone it to, and type the following command (inserting your GitHub username to its place): ``git clone https://github.com/your_github_username/JuliaFEM.jl.git``
 
 6) You can now navigate to JuliaFEM in the folder you chose at step 5. There you'll find the same contents as you see in your GitHub JuliaFEM repository. Now, locate the file you want to modify, open it with your desired text editor, make the changes and save the new version. If you type ``git status``, you'll see that the files you've created or modified are listed under ``untracked files``.
 
@@ -37,48 +26,53 @@ Here are the basic steps for contributing to JuliaFEM:
 
 11) If other contributors ask you to make changes to your pull request, just repeat steps 6-9. Your commits will be updated to your original pull request. Do this until everyone is satisfied and your pull request can be merged to the master branch.
 
+There's also some GUI apps to use git if you don't feel command line comfortable. For OSX and Windows a good application is `SourceTree <https://www.sourcetreeapp.com>`_, for Linux, maybe `SmartGit <http://www.syntevo.com/smartgit/>`_
+
+Developing on local machine
+---------------------------
+
+To set up ready for development, git clone it to your development directory and make symbolic link to julia package directory:
+
+.. code-block:: bash
+
+   cd ~/dev
+   git clone https://github.com/JuliaFEM/JuliaFEM.jl
+   cd ~/.julia/v0.4
+   ln -s ~/dev/JuliaFEM
 
-Developing
-----------
-```bash
-cd ~dev/
-git clone https://github.com/JuliaFEM/JuliaFEM.jl
-cd ~/.julia/v0.4
-ln -s ~/dev/JuliaFEM .
-```
 
 Use of UTF-8 characters in program code
 ---------------------------------------
-We have decided not to use them. [See issue #18](https://github.com/JuliaFEM/JuliaFEM.jl/issues/18).
+We have decided not to use them. See issue `#18 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/18>`_.
 
 Supported Julia versions
 ------------------------
-We support Julia versions 0.4+. [See issue #26](https://github.com/JuliaFEM/JuliaFEM.jl/issues/26)
+We support Julia versions 0.4+. See issue `#26 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/26>`_.
 
 Only pull requests to src folder
 --------------------------------
-See [issue #29](https://github.com/JuliaFEM/JuliaFEM.jl/issues/29). This ensures peer review check for contributors and hopefully will decrease the number of merge conflicts. Before making the pull request runn all test: either type `julia> Pkg.test("JuliaFEM")` at REPL or `julia test/runtests.jl` at command line. 
+See issue `#29 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/29>`_. This ensures peer review check for contributors and hopefully will decrease the number of merge conflicts. Before making the pull request runn all test: either type ``julia> Pkg.test("JuliaFEM")`` at REPL or ``julia test/runtests.jl`` at command line. 
 
 New technology should be introduced through notebooks
 -----------------------------------------------------
-[See issue #12](https://github.com/JuliaFEM/JuliaFEM.jl/issues/12). Idea is to introduce new technology as a notebook for the very beginning. Then when it's get mature the notebook will serve functional test for the matter. All notebooks will be included as examples to the documentation. 
+See issue `#12 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/12>`_. Idea is to introduce new technology as a notebook for the very beginning. Then when it's get mature the notebook will serve functional test for the matter. All notebooks will be included as examples to the documentation. 
 
 FactCheck.jl is used to write test for the JuliaFEM.jl package
 --------------------------------------------------------------
-[See issue #27](https://github.com/JuliaFEM/JuliaFEM.jl/issues/27). Use FactCheck.jl package to write the tests. We believe Test Driven Development thus 100 % test coverage is expected. 
+See issue `#27 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/27>`_. Use FactCheck.jl package to write the tests. We believe Test Driven Development thus 100 % test coverage is expected. 
 
 JuliaFEM.jl is using Logging.jl
 -------------------------------
-[See issue #25](https://github.com/JuliaFEM/JuliaFEM.jl/issues/25). We have written a test to check all sources in src folder to find any print statements. Use Logging.jl instead of println().
+See issue `#25 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/25>`_. We have written a test to check all sources in src folder to find any print statements. Use Logging.jl instead of println().
 
 Code indentation
 ----------------
-We use 4 spaces like in Python. [See issue #5](https://github.com/JuliaFEM/JuliaFEM.jl/issues/5)
+We use 4 spaces like in Python. See issue `#5 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/5>`_.
 
 Function docstrings
 -------------------
-We use numpy documentation style in our functions. See https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt. See issue #5.
+We use numpy documentation style in our functions. See `guide <https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt>`_. See issue `#5 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/5>`_.
 
 Documentation
 -------------
-We use restructured text to document this project. Information how to write rst format is described `here <http://sphinx-doc.org/rest.html>`_. `See issue #49 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/49>`_.
+We use restructured text to document this project. Information how to write rst format is described `here <http://sphinx-doc.org/rest.html>`_. See issue `#49 <https://github.com/JuliaFEM/JuliaFEM.jl/issues/49>`_.
diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb
index 59e621c..e9208e6 100644
--- a/notebooks/2015-06-25-elasticity-solver-example.ipynb
+++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb
@@ -36,48 +36,27 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 23,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
-      " in int32 at deprecated.jl:49\n",
-      " in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
-      " in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
-      " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
-      " in anonymous at task.jl:365\n",
-      "while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n",
-      "WARNING: int32(x) is deprecated, use Int32(x) instead.\n",
-      " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
-      " in int32 at deprecated.jl:49\n",
-      " in recv at /Users/jukka/.julia/v0.4/ZMQ/src/ZMQ.jl:617\n",
-      " in recv_ipython at /Users/jukka/.julia/v0.4/IJulia/src/msg.jl:63\n",
-      " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:120\n",
-      " in anonymous at task.jl:365\n",
-      "while loading /Users/jukka/.julia/v0.4/IJulia/src/kernel.jl, in expression starting on line 35\n",
-      "24-Aug 00:56:48:INFO:root:loading types\n",
-      "24-Aug 00:56:48:INFO:root:loading elements\n"
-     ]
-    },
     {
      "data": {
       "text/plain": [
        "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
       ]
      },
-     "execution_count": 1,
+     "execution_count": 23,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
+    "using ForwardDiff\n",
     "using JuliaFEM\n",
-    "using JuliaFEM: Assembly, Element, Quad4, get_integration_points, get_jacobian, get_basis, get_dbasisdxi, get_dbasisdX\n",
+    "using JuliaFEM: Assembly, Element, Equation, Point1, Quad4, IntegrationPoint\n",
+    "using JuliaFEM: get_field, set_field, get_element, linearize, integrate, interpolate, get_dbasisdX\n",
     "using Logging\n",
     "Logging.configure(level=DEBUG)"
    ]
@@ -100,59 +79,6 @@
     "First we construct some type for our element which contains all relevant data. We don't care a much how every element is actually implemented as long as it follows some general rules how the interface is constructed. First we define our element family and it's basis functions, derivatives of them etc. These needs to be defined for each element type only once."
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 2,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_rhs (generic function with 1 method)"
-      ]
-     },
-     "execution_count": 2,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"\n",
-    "Integrate f over element using Gaussian quadrature rules.\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "el::Element\n",
-    "    well defined element\n",
-    "f::Function\n",
-    "    Function to integrate\n",
-    "\"\"\"\n",
-    "function integrate(el::Element, f::Function)\n",
-    "    target = []\n",
-    "    for ip in get_integration_points(el)\n",
-    "        J = get_jacobian(el, ip.xi)\n",
-    "        push!(target, ip.weight*f(el, ip)*det(J))\n",
-    "    end\n",
-    "    return sum(target)\n",
-    "end\n",
-    "\n",
-    "\"\"\"\n",
-    "Return left hand side of the equation Ax = b (i.e. A)\n",
-    "\"\"\"\n",
-    "function get_lhs(el::Element)\n",
-    "    return None\n",
-    "end\n",
-    "\n",
-    "\"\"\"\n",
-    "Return right hand side of the equation Ax = b (i.e. b)\n",
-    "\"\"\"\n",
-    "function get_rhs(el::Element)\n",
-    "    return None\n",
-    "end"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -162,60 +88,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "#abstract Mechanical <: CG"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 4,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_field (generic function with 1 method)"
-      ]
-     },
-     "execution_count": 4,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"\n",
-    "(4-node bilinear) plane stress element\n",
-    "\"\"\"\n",
-    "type CPS4 <: Quad4\n",
-    "    id :: Int\n",
-    "    node_ids :: Array{Int, 1}\n",
-    "    coordinates :: Array{Float64, 2}\n",
-    "    integration_points :: Array{IntegrationPoint, 1}\n",
-    "    attributes :: Dict{ASCIIString, Any}\n",
-    "end\n",
-    "\n",
-    "function get_field(el::Quad4)\n",
-    "    el.attributes[\"displacement\"]\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Each element must provide a standard way how it's initialized. For constructor we need only unique element id and it's connectivity data."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 24,
    "metadata": {
     "collapsed": false
    },
@@ -226,21 +99,29 @@
        "CPS4"
       ]
      },
-     "execution_count": 5,
+     "execution_count": 24,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "function CPS4(element_id, node_ids)\n",
+    "abstract Elasticity <: Equation\n",
+    "\n",
+    "\"\"\"\n",
+    "Plane stress formulation for 4-node bilinear element.\n",
+    "\"\"\"\n",
+    "type CPS4 <: Elasticity\n",
+    "    element :: Quad4\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "end\n",
+    "function CPS4(el::Quad4)\n",
     "    integration_points = [\n",
-    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0, Dict()),\n",
-    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0, Dict()),\n",
-    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0, Dict()),\n",
-    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0, Dict())]\n",
-    "    attributes = Dict(\"displacement\" => zeros(2, 4))\n",
-    "    coordinates = zeros(2, 4)\n",
-    "    CPS4(element_id, node_ids, coordinates, integration_points, attributes)\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
+    "    set_field(el, \"displacement\", zeros(2, 4))\n",
+    "    CPS4(el, integration_points)\n",
     "end"
    ]
   },
@@ -250,12 +131,12 @@
    "source": [
     "Here's the actual basic implementation for mechanical elements\n",
     "\n",
-    "Our task is: for given $\\mathbf{u}$ calculate $\\mathbf{R}(\\mathbf{u}) = \\mathbf{T}(\\mathbf{u}) - \\mathbf{F}(\\mathbf{u})$ and it's partial derivative with respect to $\\mathbf{u}$, i.e. $\\partial \\mathbf{R}(\\mathbf{u}) / \\partial \\mathbf{u}$. Here is our $\\mathbf{T}$:"
+    "Our task is: for given $\\mathbf{u}$ calculate $\\mathbf{R}(\\mathbf{u}) = \\mathbf{T}(\\mathbf{u}) - \\mathbf{F}(\\mathbf{u})$ and it's partial derivative with respect to $\\mathbf{u}$, i.e. $\\partial \\mathbf{R}(\\mathbf{u}) / \\partial \\mathbf{u}$. Because of DRY, we apply this implementation for union of several elements because this works for all 2d and 3d continuum elements."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 25,
    "metadata": {
     "collapsed": false
    },
@@ -266,7 +147,7 @@
        "get_lhs (generic function with 2 methods)"
       ]
      },
-     "execution_count": 6,
+     "execution_count": 25,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -275,8 +156,8 @@
     "\"\"\"\n",
     "Calculate internal nodal forces for continuum element.\n",
     "\"\"\"\n",
-    "function Wint(el::Quad4)\n",
-    "\n",
+    "function Wint(eq::CPS4)\n",
+    "    el = eq.element\n",
     "    dNdX(xi) = get_dbasisdX(el, xi)\n",
     "    # material\n",
     "    lambda(xi) = interpolate(el, \"lambda\", xi)\n",
@@ -289,11 +170,11 @@
     "    S(xi, u) = lambda(xi)*trace(E(xi, u))*I + 2*mu(xi)*E(xi, u)\n",
     "    P(xi, u) = F(xi, u)*S(xi, u)\n",
     "    T(xi, u) = P(xi, u)*dNdX(xi)'\n",
-    "    integrate(el, (el, ip) -> T(ip.xi, get_field(el)))\n",
+    "    integrate(eq, (eq, ip) -> T(ip.xi, get_field(get_element(eq), \"displacement\")))\n",
     "end\n",
     "\n",
-    "get_rhs(el::Quad4) = -Wint(el)  # rhs = -R = -(T-F)\n",
-    "get_lhs(el::Quad4) = linearize(Wint, \"displacement\")(el)"
+    "get_rhs(eq::CPS4) = -Wint(eq)  # rhs = -R = -(T-F)\n",
+    "get_lhs(eq::CPS4) = linearize(eq, Wint, \"displacement\")(eq)"
    ]
   },
   {
@@ -301,13 +182,19 @@
    "metadata": {},
    "source": [
     "That was our geometrically nonlinear elasticity solver. Note how we used automatic differentiation to linearize residual vector.\n",
-    "\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
     "*Design principle 6*: we test our code. We use FactCheck for testing."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 26,
    "metadata": {
     "collapsed": true
    },
@@ -318,28 +205,36 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 27,
    "metadata": {
     "collapsed": false
    },
    "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "test rhs\n",
+      "1 fact verified.\n"
+     ]
+    },
     {
      "data": {
       "text/plain": [
-       "get_test_element (generic function with 1 method)"
+       "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 8,
+     "execution_count": 27,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "function get_test_element()\n",
+    "function get_test_equation()\n",
     "    # set up one linear quadrangle element\n",
-    "    element_id = 1\n",
     "    node_ids = [1, 2, 3, 4]\n",
-    "    el = CPS4(element_id, node_ids)\n",
+    "    el = Quad4(node_ids)\n",
+    "    eq = CPS4(el)\n",
     "\n",
     "    E = 90.0\n",
     "    nu = 0.25\n",
@@ -347,15 +242,27 @@
     "    la = E*nu/((1+nu)*(1-2*nu))\n",
     "    la = 2*la*mu/(la + 2*mu)\n",
     "    X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
-    "    set_coordinates(el, X)\n",
-    "    set_material(el, la, mu)\n",
-    "    return el\n",
+    "    #set_coordinates(el, X)\n",
+    "    #set_material(el, la, mu)\n",
+    "    set_field(el, \"coordinates\", X)\n",
+    "    set_field(el, \"lambda\", la)\n",
+    "    set_field(el, \"mu\", mu)\n",
+    "    return eq\n",
+    "end\n",
+    "facts(\"test rhs\") do\n",
+    "    eq = get_test_equation()\n",
+    "    utest = zeros(2, 4)\n",
+    "    utest[1, 3] = 0.5\n",
+    "    set_field(get_element(eq), \"displacement\", utest)\n",
+    "    rhs = get_rhs(eq)\n",
+    "    #println(rhs)\n",
+    "    @fact rhs[2, 4] --> roughly(-8.4)\n",
     "end"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 28,
    "metadata": {
     "collapsed": false,
     "scrolled": false
@@ -372,40 +279,40 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "24-Aug 00:57:07:DEBUG:root:Iteration 1\n",
-      "24-Aug 00:57:10:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Iteration 2\n",
-      "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Iteration 3\n",
-      "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Iteration 4\n",
-      "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Iteration 5\n",
-      "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Iteration 6\n",
-      "24-Aug 00:57:13:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:13:DEBUG:root:Converged in 6 iterations.\n",
-      "24-Aug 00:57:14:DEBUG:root:solution vector: \n",
+      "27-Aug 08:16:52:DEBUG:root:Iteration 1\n",
+      "27-Aug 08:16:53:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Iteration 2\n",
+      "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Iteration 3\n",
+      "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Iteration 4\n",
+      "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Iteration 5\n",
+      "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Iteration 6\n",
+      "27-Aug 08:16:55:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:55:DEBUG:root:Converged in 6 iterations.\n",
+      "27-Aug 08:16:56:DEBUG:root:solution vector: \n",
       " [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n",
       " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
-      "24-Aug 00:57:14:DEBUG:root:norm of u: 3.1292483947150047\n",
-      "24-Aug 00:57:14:DEBUG:root:Iteration 1\n",
-      "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:14:DEBUG:root:Iteration 2\n",
-      "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:14:DEBUG:root:Iteration 3\n",
-      "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:14:DEBUG:root:Iteration 4\n",
-      "24-Aug 00:57:14:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:14:DEBUG:root:Iteration 5\n",
-      "24-Aug 00:57:15:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:15:DEBUG:root:Iteration 6\n",
-      "24-Aug 00:57:15:DEBUG:root:Solving Ax = b\n",
-      "24-Aug 00:57:15:DEBUG:root:Converged in 6 iterations.\n",
-      "24-Aug 00:57:15:DEBUG:root:solution vector: \n",
+      "27-Aug 08:16:56:DEBUG:root:norm of u: 3.1292483947150047\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 1\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 2\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 3\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 4\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 5\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Iteration 6\n",
+      "27-Aug 08:16:56:DEBUG:root:Solving Ax = b\n",
+      "27-Aug 08:16:56:DEBUG:root:Converged in 6 iterations.\n",
+      "27-Aug 08:16:56:DEBUG:root:solution vector: \n",
       " [0.0 1.2578327758133292 1.5202505368695098 0.0\n",
       " 0.0 -1.8223091343697626 -1.6224781337179326 0.0]\n",
-      "24-Aug 00:57:15:DEBUG:root:norm of u: 3.129248394715004\n"
+      "27-Aug 08:16:56:DEBUG:root:norm of u: 3.129248394715004\n"
      ]
     },
     {
@@ -421,7 +328,7 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 9,
+     "execution_count": 28,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -429,7 +336,9 @@
    "source": [
     "facts(\"test solve one element model\") do\n",
     "\n",
-    "    e = get_test_element()\n",
+    "    eq = get_test_equation()\n",
+    "    el = get_element(eq)\n",
+    "    #set_field(el, \"displacement\", zeros(2, 4))\n",
     "    F = [0.0 0.0; 0.0 0.0; 0.0 -2.0; 0.0 0.0]'\n",
     "\n",
     "    du = zeros(2, 4)\n",
@@ -437,13 +346,13 @@
     "    free_dofs = [3, 4, 5, 6]\n",
     "    for i=1:10\n",
     "        Logging.debug(\"Iteration $i\")\n",
-    "        b = get_rhs(e)\n",
-    "        #Logging.debug(\"rhs = $b\")\n",
-    "        A = get_lhs(e)\n",
+    "        b = get_rhs(eq)\n",
+    "        A = get_lhs(eq)\n",
     "        Logging.debug(\"Solving Ax = b\")\n",
     "        du[free_dofs] = A[free_dofs, free_dofs] \\ (b + F)[free_dofs]\n",
     "\n",
-    "        e.attributes[\"displacement\"] += du\n",
+    "        field = get_field(el, \"displacement\") + du\n",
+    "        set_field(el, \"displacement\", field)\n",
     "        if norm(du) < 1.0e-9\n",
     "            Logging.debug(\"Converged in $i iterations.\")\n",
     "            break\n",
@@ -451,7 +360,7 @@
     "    end\n",
     "\n",
     "    # Tested against Elmer solution\n",
-    "    u = e.attributes[\"displacement\"]\n",
+    "    u = get_field(el, \"displacement\")\n",
     "    Logging.debug(\"solution vector: \\n $u\")\n",
     "    @fact u[2, 3] --> roughly(-2.222244754401764)\n",
     "    norm1 = norm(u)\n",
@@ -462,25 +371,26 @@
     "    rmat = [\n",
     "        cos(phi) -sin(phi)\n",
     "        sin(phi)  cos(phi)]\n",
-    "    set_coordinates(e, rmat*get_coordinates(e))\n",
+    "    set_field(el, \"coordinates\", rmat*get_field(el, \"coordinates\"))\n",
     "    F = rmat*F\n",
     "\n",
-    "    e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n",
+    "    set_field(el, \"displacement\", [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]')\n",
     "    du = zeros(2, 4)\n",
     "    for i=1:10\n",
     "        Logging.debug(\"Iteration $i\")\n",
-    "        A = get_lhs(e)\n",
-    "        b = get_rhs(e)\n",
+    "        A = get_lhs(eq)\n",
+    "        b = get_rhs(eq)\n",
     "        Logging.debug(\"Solving Ax = b\")\n",
     "        du[free_dofs] = A[free_dofs, free_dofs] \\ (b + F)[free_dofs]\n",
     "\n",
-    "        e.attributes[\"displacement\"] += du\n",
+    "        field = get_field(el, \"displacement\") + du\n",
+    "        set_field(el, \"displacement\", field)\n",
     "        if norm(du) < 1.0e-9\n",
     "            Logging.debug(\"Converged in $i iterations.\")\n",
     "            break\n",
     "        end\n",
     "    end\n",
-    "    u = e.attributes[\"displacement\"]\n",
+    "    u = get_field(el, \"displacement\")\n",
     "    Logging.debug(\"solution vector: \\n $u\")\n",
     "    Logging.debug(\"norm of u: $(norm(u))\")\n",
     "    @fact norm(u) --> roughly(norm1) \n",
@@ -496,7 +406,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 29,
    "metadata": {
     "collapsed": false
    },
@@ -507,17 +417,19 @@
        "assemble_rhs! (generic function with 1 method)"
       ]
      },
-     "execution_count": 10,
+     "execution_count": 29,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "function assemble_lhs!(ass::Assembly, el::Element)\n",
+    "function assemble_lhs!(ass::Assembly, eq::Equation)\n",
     "\n",
-    "    gdofs = ass.gdofs[el.id]\n",
+    "    el::Element = get_element(eq)\n",
+    "    elid = get_field(el, \"id\")\n",
+    "    gdofs = ass.gdofs[elid]\n",
     "\n",
-    "    A = get_lhs(el)\n",
+    "    A = get_lhs(eq)\n",
     "    if !(A == None)\n",
     "        ii, jj = size(A)\n",
     "        for i=1:ii\n",
@@ -530,11 +442,14 @@
     "    end\n",
     "end\n",
     "\n",
-    "function assemble_rhs!(ass::Assembly, el::Element)\n",
+    "function assemble_rhs!(ass::Assembly, eq::Equation)\n",
+    "\n",
+    "    el::Element = get_element(eq)\n",
+    "    elid = get_field(el, \"id\")\n",
+    "    gdofs = ass.gdofs[elid]\n",
     "\n",
-    "    gdofs = ass.gdofs[el.id]\n",
     "    \n",
-    "    b = get_rhs(el)\n",
+    "    b = get_rhs(eq)\n",
     "    if !(b == None)\n",
     "        for i=1:length(b)\n",
     "            push!(ass.i, gdofs[i])\n",
@@ -546,7 +461,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 30,
    "metadata": {
     "collapsed": false
    },
@@ -554,129 +469,30 @@
     {
      "data": {
       "text/plain": [
-       "get_field (generic function with 2 methods)"
+       "get_lhs (generic function with 2 methods)"
       ]
      },
-     "execution_count": 11,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "function update_field(el::Element, du)\n",
-    "    el.attributes[\"displacement\"][:] += du\n",
-    "end\n",
-    "function get_field(el::Element)\n",
-    "    return el.attributes[\"displacement\"]\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "We also need to construct our element in somehow \"standard\" way. My proposal is: element id and node ids (connectivity) information. Of course other fields must also be provided. Here's example for CPS4 element:"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Almost done. Also we must define interface how to set coordinates, material properties etc. for element."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "set_attribute (generic function with 1 method)"
-      ]
-     },
-     "execution_count": 12,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "function set_attribute(el::Element, field, value)\n",
-    "    el.attributes[field] = value\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "This might look a bit cumbersome at this point, but notice that everything is very general so far. We really don't have to define all this stuff when creating new elements if we follow some general construct. To demonstrate that, we define 0-dimensional \"point force\" element. Keep in mind, we need to have some functions defined how to operate using element, but they are already defined to superclass of element. So we need to only take care of construcor, rhs and (maybe) lhs. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "using JuliaFEM: CG\n",
-    "abstract Point0 <: CG\n",
     "\"\"\"\n",
     "1-node point force element for plane stress problems.\n",
     "\"\"\"\n",
-    "type CPS1 <: Point0\n",
-    "    id :: Int\n",
-    "    node_ids :: Array{Int, 1}\n",
-    "    attributes :: Dict{ASCIIString, Any}\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 14,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_lhs (generic function with 3 methods)"
-      ]
-     },
-     "execution_count": 14,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"\n",
-    "Constructor of CPS1\n",
-    "\"\"\"\n",
-    "function CPS1(element_id, node_ids)\n",
-    "    attributes = Dict(\n",
-    "        \"displacement\" => zeros(2, 1),\n",
-    "        \"displacement nodal load\" => zeros(2, 1))\n",
-    "    CPS1(element_id, node_ids, attributes)\n",
+    "type CPS1 <: Elasticity\n",
+    "    element :: Point1\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
     "end\n",
-    "function get_rhs(el::CPS1)\n",
-    "    return el.attributes[\"displacement nodal load\"]\n",
+    "function CPS1(el::Point1)\n",
+    "    integration_points = []\n",
+    "    set_field(el, \"displacement\", zeros(2, 1))\n",
+    "    set_field(el, \"displacement nodal load\", zeros(2, 1))\n",
+    "    CPS1(el, integration_points)\n",
     "end\n",
-    "function get_lhs(el::CPS1)\n",
-    "    return None\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "And that's basically everything needs to be defined. During assembling get_lhs returns nothing and does not assemble anything to stiffness matrix (of course in case of follower point force the direction depents on the normal ..) and get_rhs adds simple point force pointing to some direction."
+    "get_rhs(eq::CPS1) = get_field(get_element(eq), \"displacement nodal load\")\n",
+    "get_lhs(eq::CPS1) = None"
    ]
   },
   {
@@ -691,7 +507,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 31,
    "metadata": {
     "collapsed": false
    },
@@ -700,45 +516,45 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "one element assembly\n"
+      "two element assembly\n"
      ]
     },
     {
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "24-Aug 01:06:39:DEBUG:root:Starting iteration 1\n",
-      "24-Aug 01:06:39:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:39:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:39:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 3.0900221367289986\n",
-      "24-Aug 01:06:40:DEBUG:root:Starting iteration 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.3212131602153472\n",
-      "24-Aug 01:06:40:DEBUG:root:Starting iteration 3\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.04043178193999703\n",
-      "24-Aug 01:06:40:DEBUG:root:Starting iteration 4\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 0.0009291101052105739\n",
-      "24-Aug 01:06:40:DEBUG:root:Starting iteration 5\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 1.5638899027804743e-7\n",
-      "24-Aug 01:06:40:DEBUG:root:Starting iteration 6\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 1\n",
-      "24-Aug 01:06:40:DEBUG:root:Assembling element 2\n",
-      "24-Aug 01:06:40:DEBUG:root:Solution norm = 1.0464940956129567e-14\n",
-      "24-Aug 01:06:40:DEBUG:root:Converged in 6 iterations.\n",
-      "24-Aug 01:06:40:DEBUG:root:Displacement of element = \n",
+      "27-Aug 08:17:00:DEBUG:root:Starting iteration 1\n",
+      "27-Aug 08:17:00:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:00:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:00:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 3.0900221367289986\n",
+      "27-Aug 08:17:01:DEBUG:root:Starting iteration 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.3212131602153472\n",
+      "27-Aug 08:17:01:DEBUG:root:Starting iteration 3\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.04043178193999703\n",
+      "27-Aug 08:17:01:DEBUG:root:Starting iteration 4\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 0.0009291101052105739\n",
+      "27-Aug 08:17:01:DEBUG:root:Starting iteration 5\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 1.5638899027804743e-7\n",
+      "27-Aug 08:17:01:DEBUG:root:Starting iteration 6\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 1\n",
+      "27-Aug 08:17:01:DEBUG:root:Assembling lhs and rhs from equation 2\n",
+      "27-Aug 08:17:01:DEBUG:root:Solution norm = 1.0464940956129567e-14\n",
+      "27-Aug 08:17:01:DEBUG:root:Converged in 6 iterations.\n",
+      "27-Aug 08:17:01:DEBUG:root:Displacement of element = \n",
       "[0.0 -0.39914506095474334 -0.0722858269559246 0.0\n",
       " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n"
      ]
@@ -756,44 +572,47 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 17,
+     "execution_count": 31,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "facts(\"one element assembly\") do\n",
+    "facts(\"two element assembly\") do\n",
     "    # set up element 1\n",
-    "    element_id = 1\n",
-    "    node_ids = [1, 2, 3, 4]\n",
-    "    el1 = CPS4(element_id, node_ids)\n",
+    "    el1 = Quad4([1, 2, 3, 4])\n",
+    "    eq1 = CPS4(el1)\n",
     "    # assign properties to element, e.g. coordinates, material properties, ...\n",
     "    E = 90.0\n",
     "    nu = 0.25\n",
     "    mu = E/(2*(1+nu))\n",
     "    la = E*nu/((1+nu)*(1-2*nu))\n",
     "    la = 2*la*mu/(la + 2*mu)\n",
-    "    set_coordinates(el1, [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n",
-    "    set_attribute(el1, \"lambda\", la)\n",
-    "    set_attribute(el1, \"mu\", mu)\n",
+    "    set_field(el1, \"coordinates\", [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')\n",
+    "    set_field(el1, \"lambda\", la)\n",
+    "    set_field(el1, \"mu\", mu)\n",
+    "    set_field(el1, \"id\", 1)\n",
     "\n",
     "    # set up element 2\n",
-    "    el2 = CPS1(2, [3])  # Create nodal point force element with id 2 for node 3\n",
-    "    set_attribute(el2, \"displacement nodal load\", [0.0, -2.0])\n",
+    "    el2 = Point1([3])\n",
+    "    eq2 = CPS1(el2)  # Create nodal point force element with id 2 for node 3\n",
+    "    set_field(el2, \"displacement nodal load\", [0.0, -2.0])\n",
+    "    set_field(el2, \"id\", 2)\n",
     "\n",
     "    elements = [el1, el2]\n",
+    "    equations = [eq1, eq2]\n",
     "\n",
     "    for i=1:10\n",
     "        Logging.debug(\"Starting iteration $i\")\n",
     "        Logging.debug(\"Assembling\")\n",
     "        ass = Assembly()\n",
-    "        ass.gdofs[el1.id] = [1, 2, 3, 4, 5, 6, 7, 8]\n",
-    "        ass.gdofs[el2.id] = [5, 6]\n",
+    "        ass.gdofs[get_field(el1, \"id\")] = [1, 2, 3, 4, 5, 6, 7, 8]\n",
+    "        ass.gdofs[get_field(el2, \"id\")] = [5, 6]\n",
     "\n",
-    "        for (j, el) in enumerate(elements)\n",
-    "            Logging.debug(\"Assembling element $j\")\n",
-    "            assemble_lhs!(ass, el)\n",
-    "            assemble_rhs!(ass, el)\n",
+    "        for (j, eq) in enumerate(equations)\n",
+    "            Logging.debug(\"Assembling lhs and rhs from equation $j\")\n",
+    "            assemble_lhs!(ass, eq)\n",
+    "            assemble_rhs!(ass, eq)\n",
     "        end\n",
     "\n",
     "        # (Dirichlet) boundary conditions \"handled\"\n",
@@ -809,15 +628,18 @@
     "\n",
     "        # update solution back to elements\n",
     "        for el in elements\n",
-    "            eldu = du[ass.gdofs[el.id]]\n",
-    "            update_field(el, eldu)\n",
+    "            eldu = du[ass.gdofs[get_field(el, \"id\")]]\n",
+    "            #tmp = get_field(el, \"displacement\")\n",
+    "            #set_field(el, \"displacement\", tmp+eldu)\n",
+    "            #update_field(el, eldu)\n",
+    "            el.fields[\"displacement\"][:] += eldu\n",
     "        end\n",
     "        if norm(du) < 1.0e-9\n",
     "            Logging.debug(\"Converged in $i iterations.\")\n",
     "            break\n",
     "        end\n",
     "    end    \n",
-    "    disp = get_field(el1)\n",
+    "    disp = get_field(el1, \"displacement\")\n",
     "    Logging.debug(\"Displacement of element = \\n$disp\")\n",
     "    @fact norm(disp) --> roughly(3.1292483947150043)\n",
     "end"
@@ -832,7 +654,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 12,
    "metadata": {
     "collapsed": false
    },
@@ -852,7 +674,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 13,
    "metadata": {
     "collapsed": false
    },
@@ -863,7 +685,7 @@
        "MPC"
       ]
      },
-     "execution_count": 19,
+     "execution_count": 13,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -879,7 +701,26 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "abstract Problem\n",
+    "\n",
+    "type PlaneStressProblem <: Problem\n",
+    "    equations :: Array{Equation, 1}\n",
+    "    boundary_conditions :: Array{BoundaryCondition, 1}\n",
+    "\n",
+    "#    dofmap :: Dict{Int64, Array{Int64,1}}()  # a dict node_id : (dof1, dof2, ...)\n",
+    "#    solver_parameters :: SolverParameters\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
    "metadata": {
     "collapsed": false
    },
@@ -890,7 +731,7 @@
        "create_ldof2gdofmap (generic function with 1 method)"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 17,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -923,40 +764,9 @@
     "end"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "Next, we want to solve a _problem_. It's yet another container and if it's defined somewhat standard way default solver can solve it. Nothing stops user to write his/hers own solver."
-   ]
-  },
   {
    "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": [
-    "type SolverParameters\n",
-    "    max_iterations :: Int64\n",
-    "    eps :: Float64\n",
-    "end\n",
-    "SolverParameters() = SolverParameters(10, 0.05)\n",
-    "\n",
-    "abstract Problem\n",
-    "abstract ElasticityProblem <: Problem\n",
-    "type PlaneStressProblem <: ElasticityProblem\n",
-    "    elements :: Array{Element, 1}\n",
-    "    dofmap :: Dict{Int64, Array{Int64,1}}()  # a dict node_id : (dof1, dof2, ...)\n",
-    "    dirichlet_bcs :: Array{BoundaryCondition, 1}\n",
-    "    solver_parameters :: SolverParameters\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 18,
    "metadata": {
     "collapsed": false,
     "scrolled": false
@@ -973,53 +783,53 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "24-Aug 01:07:29:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n",
-      "24-Aug 01:07:29:INFO:root:solve!: dofs per node: 2\n",
-      "24-Aug 01:07:29:DEBUG:root:Problem size = 8\n",
-      "24-Aug 01:07:29:DEBUG:root:Starting iteration 1\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:29:DEBUG:root:Solution norm du = 3.0900221367289444\n",
-      "24-Aug 01:07:29:DEBUG:root:Starting iteration 2\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.32121316021534796\n",
-      "24-Aug 01:07:29:DEBUG:root:Starting iteration 3\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.040431781940014504\n",
-      "24-Aug 01:07:29:DEBUG:root:Starting iteration 4\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:29:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:29:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:29:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:29:DEBUG:root:Solution norm du = 0.0009291101052065917\n",
-      "24-Aug 01:07:30:DEBUG:root:Starting iteration 5\n",
-      "24-Aug 01:07:30:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:30:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:30:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:30:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:30:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:30:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n",
-      "24-Aug 01:07:30:DEBUG:root:Starting iteration 6\n",
-      "24-Aug 01:07:30:DEBUG:root:Assembling lhs\n",
-      "24-Aug 01:07:30:DEBUG:root:Assembling rhs\n",
-      "24-Aug 01:07:30:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "24-Aug 01:07:30:DEBUG:root:Added 4 Lagrange multipliers to model\n",
-      "24-Aug 01:07:30:DEBUG:root:Solving system of equations. Total size = 12\n",
-      "24-Aug 01:07:30:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n",
-      "24-Aug 01:07:30:DEBUG:root:Converged in 6 iterations.\n",
-      "24-Aug 01:07:30:DEBUG:root:Displacement on upper right = \n",
+      "24-Aug 18:34:38:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],1=>[1,2])\n",
+      "24-Aug 18:34:38:INFO:root:solve!: dofs per node: 2\n",
+      "24-Aug 18:34:39:DEBUG:root:Problem size = 8\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 1\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 3.0900221367289444\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 2\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.32121316021534796\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 3\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.040431781940014504\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 4\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 0.0009291101052065917\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 5\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 1.5638899136781228e-7\n",
+      "24-Aug 18:34:39:DEBUG:root:Starting iteration 6\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling lhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Assembling rhs\n",
+      "24-Aug 18:34:39:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "24-Aug 18:34:39:DEBUG:root:Added 4 Lagrange multipliers to model\n",
+      "24-Aug 18:34:39:DEBUG:root:Solving system of equations. Total size = 12\n",
+      "24-Aug 18:34:39:DEBUG:root:Solution norm du = 1.0913504694802626e-14\n",
+      "24-Aug 18:34:39:DEBUG:root:Converged in 6 iterations.\n",
+      "24-Aug 18:34:39:DEBUG:root:Displacement on upper right = \n",
       "[-0.07228582695592467\n",
       " -2.222244754401765]\n"
      ]
@@ -1037,7 +847,7 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 22,
+     "execution_count": 18,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -2020,6 +1830,157 @@
    "source": [
     "d.Image(\"/tmp/piston.png\")"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 43,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "using ForwardDiff"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 51,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "f(x::Vector) = sum(sin, x) + prod(tan, x) * sum(sqrt, x);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 52,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5-element Array{Float64,1}:\n",
+       " 0.903617\n",
+       " 0.97582 \n",
+       " 0.36872 \n",
+       " 0.233261\n",
+       " 0.213665"
+      ]
+     },
+     "execution_count": 52,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "x = rand(5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 53,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2.5477827220775753"
+      ]
+     },
+     "execution_count": 53,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "f(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 54,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "5-element Array{Float64,1}:\n",
+       " 0.906762\n",
+       " 0.860405\n",
+       " 1.35153 \n",
+       " 1.59159 \n",
+       " 1.64713 "
+      ]
+     },
+     "execution_count": 54,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "g = ForwardDiff.gradient(f);\n",
+    "g(x)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 55,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Function"
+      ]
+     },
+     "execution_count": 55,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "typeof(g)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 56,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "ename": "LoadError",
+     "evalue": "LoadError: MethodError: `grad` has no method matching grad(::Function)\nwhile loading In[56], in expression starting on line 1",
+     "output_type": "error",
+     "traceback": [
+      "LoadError: MethodError: `grad` has no method matching grad(::Function)\nwhile loading In[56], in expression starting on line 1",
+      ""
+     ]
+    }
+   ],
+   "source": [
+    "j = ForwardDiff.jacobian(g)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb
new file mode 100644
index 0000000..b876d35
--- /dev/null
+++ b/notebooks/2015-08-29-developing-juliafem.ipynb
@@ -0,0 +1,521 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Developing JuliaFEM\n",
+    "\n",
+    "Author(s): Jukka Aho\n",
+    "\n",
+    "**Abstract**: General developer notes."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using Logging\n",
+    "using FactCheck\n",
+    "Logging.configure(level=DEBUG)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Developing own element\n",
+    "\n",
+    "Finite element definition, from [FEniCS-book](https://bitbucket.org/fenics-project/fenics-book/src/7d3a80e7dda0fc279c7964dc6000d57942f11eb3/fenicsbook.cls?at=master) [Ciarlet, 2002]:\n",
+    "\n",
+    "-  the domain $T$ is a bounded, closed subset of $\\mathbb{R}^d$ (for $d = 1, 2, 3, \\dots$) with nonempty interior and piecewise smooth boundary;\n",
+    "- the space $\\mathcal{V} = \\mathcal{V}(T)$ is a finite dimensional function space on $T$ of dimension $n$;\n",
+    "- the set of degrees of freedom (nodes) $\\mathcal{L} = \\{\\ell_1, \\ell_2,\\ldots, \\ell_{n}\\}$ is a basis for the dual space $\\mathcal{V}'$; that is, the space of bounded linear functionals on $\\mathcal{V}$.\n",
+    "\n",
+    "We extend this definition so that domain $T$ can also be empty.\n",
+    "\n",
+    "Minimum requirements for element:\n",
+    "- subclass from Element, if not wanting to implement everything by youself\n",
+    "- define basis and partial derivatives of it, because we need to interpolate over it\n",
+    "- give connectivity information, how this element is connected to other elements\n",
+    "- create proper constructor (see example).\n",
+    "\n",
+    "Test the element using ``test_element`` function. It it passes, then element implementation should be fine. As an example, we define 4 node quadrilateral element using linear Lagrange basis. We really don't care much how element is implemented as long it's interface is constructed with some rules. The interface is tested using `test_element` and it also gives information how to fix element if something is missing."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "using JuliaFEM: Element, Equation, set_field, get_field, interpolate, integrate_lhs, integrate_rhs"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here's one basic implementation. The actual element:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "type Quad4 <: Element\n",
+    "    connectivity :: Array{Int, 1}\n",
+    "    fields :: Dict{Any, Any}\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Default constructor, providing connectivity data needed in assembly"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Quad4"
+      ]
+     },
+     "execution_count": 4,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "Quad4(connectivity) = Quad4(connectivity, Dict{ASCIIString, Any}())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Some basic charasteristics like number of nodes / connectivity points and dimension:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "get_element_dimension (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "JuliaFEM.get_number_of_nodes(el::Type{Quad4}) = 4\n",
+    "JuliaFEM.get_element_dimension(el::Quad4) = 2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The most important, it's basis (we probably want to interpolate something with this element):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "get_dbasisdxi (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "function JuliaFEM.get_basis(el::Quad4, xi::Array{Float64,1})\n",
+    "    [(1-xi[1])*(1-xi[2])/4\n",
+    "     (1+xi[1])*(1-xi[2])/4\n",
+    "     (1+xi[1])*(1+xi[2])/4\n",
+    "     (1-xi[1])*(1+xi[2])/4]\n",
+    "end\n",
+    "\n",
+    "function JuliaFEM.get_dbasisdxi(el::Quad4, xi::Array{Float64,1})\n",
+    "    [-(1-xi[2])/4.0    -(1-xi[1])/4.0\n",
+    "      (1-xi[2])/4.0    -(1+xi[1])/4.0\n",
+    "      (1+xi[2])/4.0     (1+xi[1])/4.0\n",
+    "     -(1+xi[2])/4.0     (1-xi[1])/4.0]\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we check that everything is well defined:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "30-Aug 23:18:00:INFO:root:number of connectivity points (nodes) in this element: 4\n",
+      "30-Aug 23:18:00:INFO:root:Constructing element..\n",
+      "30-Aug 23:18:00:INFO:root:Element dimension: 2\n",
+      "30-Aug 23:18:01:INFO:root:Setting scalar field [1 2 3 4] to element.\n",
+      "30-Aug 23:18:01:INFO:root:Interpolating scalar field\n",
+      "30-Aug 23:18:01:INFO:root:Element Quad4 passed tests.\n"
+     ]
+    }
+   ],
+   "source": [
+    "using JuliaFEM: test_element\n",
+    "test_element(Quad4)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If test_element passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "el1 = Quad4([1, 2, 3, 4])\n",
+    "set_field(el1, :temperature, [1 2 3 4])\n",
+    "set_field(el1, :coordinates, [0.0 0.0 0.0; 10.0 0.0 0.0; 10.0 1.0 0.0; 0.0 1.0 0.0]');\n",
+    "set_field(el1, :\"heat coefficient\", 1);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1-element Array{Float64,1}:\n",
+       " 2.5"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interpolate(el1, :temperature, [0.0, 0.0]) # temperature at the middle poinf of the element, 1/4*(1+2+3+4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "3-element Array{Float64,1}:\n",
+       " 5.0\n",
+       " 0.5\n",
+       " 0.0"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interpolate(el1, :coordinates, [0.0, 0.0]) # midpoint of element"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "interpolate(el1, :\"heat coefficient\", [0.0, 0.0])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Developing own formulation\n",
+    "\n",
+    "Let's consider Poisson equation\n",
+    "\\begin{align}\n",
+    "\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\n",
+    "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
+    "\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}}\n",
+    "\\end{align}\n",
+    "\n",
+    "Weak form is, find $u\\in\\mathcal{U}$ such that\n",
+    "\\begin{equation}\n",
+    "    \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Gamma_{\\mathrm{N}}}g v\\,\\mathrm{d}s \\quad \\forall v\\in\\mathcal{V}.\n",
+    "\\end{equation}\n",
+    "\n",
+    "Minimum requirements for equation: \n",
+    "- subclass from Equation, if not want to implement from scratch\n",
+    "- it needs to have lhs and rhs functions\n",
+    "- default constructor takes the element as input argument\n",
+    "\n",
+    "Now we have function `test_equation`, which we can use to test that everything is working as expected. \n",
+    "\n",
+    "Again thanks to multiple dispatch, you are free to code your weak form however you want as long as it returns lhs and rhs sides for element dofs. This kind of freedom gives good opportunities to wrap e.g. Fortran code from some other projects. And again we have some suggestions ad following these ideas you get a lot of stuff for free. First we look the left hand side of the equation, that is,\n",
+    "\\begin{equation}\n",
+    "    \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x\n",
+    "\\end{equation}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "using JuliaFEM: Equation, IntegrationPoint\n",
+    "\n",
+    "abstract Heat <: Equation"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Our basic data type often looks something like this:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "\"\"\"\n",
+    "Diffusive heat transfer for 4-node bilinear element.\n",
+    "\"\"\"\n",
+    "type DC2D4 <: Heat\n",
+    "    element :: Quad4\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We must provide default constructor which takes element as input argument:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "DC2D4"
+      ]
+     },
+     "execution_count": 14,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "function DC2D4(el::Quad4)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
+    "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
+    "    set_field(el, \"temperature\", zeros(2, 4))\n",
+    "    DC2D4(el, integration_points, [])\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now the actual implementation for $\\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x$:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "get_lhs (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using JuliaFEM: get_element, get_dbasisdX\n",
+    "\n",
+    "\"\"\"\n",
+    "Left hand side defined in integration point\n",
+    "\"\"\"\n",
+    "function JuliaFEM.get_lhs(eq::DC2D4, ip)\n",
+    "    el = get_element(eq)\n",
+    "    dNdX = get_dbasisdX(el, ip.xi)\n",
+    "    hc = interpolate(el, :\"temperature heat coefficient\", ip.xi)\n",
+    "    return dNdX*hc*dNdX'\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "And that's it. If we want to play with this formulation, we must create element and assign this equation for it:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "4x4 Array{Float64,2}:\n",
+       "  4.0  -1.0  -2.0  -1.0\n",
+       " -1.0   4.0  -1.0  -2.0\n",
+       " -2.0  -1.0   4.0  -1.0\n",
+       " -1.0  -2.0  -1.0   4.0"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "el = Quad4([1, 2, 3, 4])\n",
+    "set_field(el, :coordinates, [0 0; 1 0; 1 1; 0 1]')\n",
+    "set_field(el, :\"temperature heat coefficient\", 6)\n",
+    "eq = DC2D4(el)\n",
+    "integrate_lhs(eq)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 0.4.0-dev",
+   "language": "julia",
+   "name": "julia-0.4"
+  },
+  "language_info": {
+   "name": "julia",
+   "version": "0.4.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}
diff --git a/src/elements.jl b/src/elements.jl
index 3371115..72ba279 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -3,13 +3,16 @@
 
 abstract Element
 
+get_element(eq::Equation) = eq.element
+
 """
 Get jacobian of element evaluated at point xi
 """
 function get_jacobian(el::Element, xi)
-    dbasisdxi(xi) = get_dbasisdxi(el, xi)
+    dbasisdxi = get_dbasisdxi(el, xi)
     X = get_field(el, "coordinates")
-    J = interpolate(X, dbasisdxi, xi)'
+    #J = interpolate(X, dbasisdxi, xi)'
+    J = X*dbasisdxi
     return J
 end
 
@@ -103,6 +106,10 @@ type Point1 <: CG
     node_ids :: Array{Int, 1}
     fields :: Dict{ASCIIString, Any}
 end
+function Point1(node_ids)
+    fields = Dict{ASCIIString, Any}()
+    Point1(node_ids, fields)
+end
 
 # 1d Lagrange elements
 
diff --git a/src/math.jl b/src/math.jl
index 5983178..6a4f0e5 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -1,9 +1,7 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-"""
-This module contains math stuff, including interpolation, integration, linearization, ...
-"""
+## This module contains math stuff, including interpolation, integration, linearization, ...
 
 using ForwardDiff
 
@@ -99,15 +97,16 @@ This version returns another function which can be then evaluated against field
 """
 function linearize(f::Function, field::ASCIIString)
     function jacobian(el::Element, args...)
-        dim, nnodes = size(el.attributes[field])
+        fld = get_field(el, field)
+        dim, nnodes = size(fld)
         function helper!(x, y)
-            orig = copy(el.attributes[field])
-            el.attributes[field] = reshape(x, dim, nnodes)
+            orig = copy(fld)
+            set_field(el, field,  reshape(x, dim, nnodes))
             y[:] = f(el, args...)
-            el.attributes[field] = copy(orig)
+            set_field(el, field, copy(orig))
         end
         jac = ForwardDiff.forwarddiff_jacobian(helper!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)
-        return jac(el.attributes[field][:])
+        return jac(fld[:])
     end
     return jacobian
 end
@@ -184,3 +183,63 @@ function integrate!(f::Function, el::Element, target)
     end
 end
 
+get_integration_points(eq::Equation) = eq.integration_points
+
+"""
+Integrate f over element using Gaussian quadrature rules.
+
+Parameters
+----------
+el::Element
+    well defined element
+f::Function
+    Function to integrate
+"""
+function integrate(eq::Equation, f::Function)
+    target = []
+    for ip in get_integration_points(eq)
+        J = get_jacobian(eq.element, ip.xi)
+        push!(target, ip.weight*f(eq, ip)*det(J))
+    end
+    return sum(target)
+end
+
+
+"""
+Evaluate field in point xi using basis functions.
+"""
+function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})
+    f = get_field(el, field)
+    if !isa(f, Array)
+        # This is scalar, nothing to interpolate
+        return f
+    end
+    basis = get_basis(el, xi)
+    dim, nnodes = size(f)
+    result = zeros(dim)
+    for i=1:nnodes
+        result += basis[i]*f[:,i]
+    end
+    return result
+end
+
+
+function linearize(eq::Equation, f::Function, field::ASCIIString)
+    function jacobian(eq::Equation, args...)
+        el = get_element(eq)
+        fld = get_field(el, field)
+        dim, nnodes = size(fld)
+        function helper(x::Vector)
+            orig = copy(fld)
+            set_field(el, field,  reshape(x, dim, nnodes))
+            y = f(eq, args...)
+            set_field(el, field, orig)
+            return y[:]
+        end
+        jac = ForwardDiff.jacobian(helper)
+        return jac(fld[:])
+    end
+    return jacobian
+end
+
+
diff --git a/src/types.jl b/src/types.jl
index 2ef95f1..c75a291 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -1,7 +1,8 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-export IntegrationPoint, Assembly
+abstract Equation
+
 
 """
 Integration point
@@ -19,7 +20,7 @@ type IntegrationPoint
     weight :: Float64
     attributes :: Dict{ASCIIString, Any}
 end
-IntegrationPoint(xi, weight) = IntegrationPoint(xi, weight, Dict{AsciiString, Any}())
+IntegrationPoint(xi, weight) = IntegrationPoint(xi, weight, Dict{ASCIIString, Any}())
 
 type Assembly
     # LHS
diff --git a/test_elements.jl b/test/test_elements.jl
similarity index 100%
rename from test_elements.jl
rename to test/test_elements.jl