diff --git a/notebooks/2015-06-25-elasticity-solver-example.ipynb b/notebooks/2015-06-25-elasticity-solver-example.ipynb
index 361a2d5..2c4110f 100644
--- a/notebooks/2015-06-25-elasticity-solver-example.ipynb
+++ b/notebooks/2015-06-25-elasticity-solver-example.ipynb
@@ -36,7 +36,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -47,35 +47,14 @@
        "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)"
       ]
      },
-     "execution_count": 1,
+     "execution_count": 2,
      "metadata": {},
      "output_type": "execute_result"
-    },
-    {
-     "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"
-     ]
     }
    ],
    "source": [
+    "using JuliaFEM\n",
     "using Logging\n",
-    "using ForwardDiff\n",
     "Logging.configure(level=DEBUG)"
    ]
   },
@@ -89,12 +68,12 @@
     "\n",
     "*Design principle 5*: we use 4 space indentation like in Python.\n",
     "\n",
-    "First we write some elementary functions to calculate stiffness matrix."
+    "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}$."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 3,
    "metadata": {
     "collapsed": false
    },
@@ -102,144 +81,56 @@
     {
      "data": {
       "text/plain": [
-       "calc_local_matrices! (generic function with 1 method)"
+       "jacobian (generic function with 1 method)"
       ]
      },
-     "execution_count": 2,
+     "execution_count": 3,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "\"\"\"\n",
-    "Calculate local tangent stiffness matrix and residual force vector\n",
-    "R = T - F for elasticity problem.\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "X : Element coordinates\n",
-    "u : Displacement field\n",
-    "R : Residual force vector\n",
-    "K : Tangent stiffness matrix\n",
-    "basis : Basis functions\n",
-    "dbasis : Derivative of basis functions\n",
-    "lambda : Material parameter\n",
-    "mu : Material parameter\n",
-    "ipoints : integration points\n",
-    "iweights : integration weights\n",
-    "\n",
-    "Returns\n",
-    "-------\n",
-    "None\n",
-    "\n",
-    "Notes\n",
-    "-----\n",
-    "If material parameters are given in list, they are interpolated to gauss\n",
-    "points using shape functions.\n",
-    "\n",
-    "Examples\n",
-    "--------\n",
-    "\n",
-    "\"\"\"\n",
-    "function calc_local_matrices2!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)\n",
-    "    dim, nnodes = size(X)\n",
-    "    I = eye(dim)\n",
-    "    R[:,:] = 0.0\n",
-    "    K[:,:] = 0.0\n",
-    "\n",
-    "    dF = zeros(dim, dim)\n",
-    "\n",
-    "    for m = 1:length(iweights)\n",
-    "        w = iweights[m]\n",
-    "        xi = ipoints[m, :]\n",
-    "        # calculate material parameters\n",
-    "        lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))\n",
-    "        mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))\n",
-    "        Jt = X*dbasis(xi)\n",
-    "        detJ = det(Jt)\n",
-    "        dbasisdX = dbasis(xi)*inv(Jt)\n",
+    "function calc_residual_vector_integrand(el::JuliaFEM.Element, xi)\n",
+    "    # Calculate dN/dX\n",
+    "    dbasisdX = JuliaFEM.get_dbasisdX(el, xi)\n",
     "\n",
+    "    # Calculate residual force vector R(u) = T(u) - F(u)\n",
+    "    function T(u)\n",
+    "        # kinematics\n",
     "        gradu = u*dbasisdX\n",
-    "        F = I + gradu  # Deformation gradient\n",
-    "        E = 1/2*(gradu' + gradu + gradu'*gradu)  # Green-Lagrange strain tensor\n",
-    "        S = lambda*trace(E)*I + 2*mu*E  # PK2 stress tensor\n",
-    "        P = F*S  # PK1 stress tensor\n",
-    "\n",
-    "        R[:,:] += w*P*dbasisdX'*detJ\n",
-    "\n",
-    "        for p = 1:nnodes\n",
-    "            for i = 1:dim\n",
-    "                dF[:,:] = 0.0\n",
-    "                dF[i,:] = dbasisdX[p,:]\n",
-    "                dE = 1/2*(F'*dF + dF'*F)\n",
-    "                dS = lambda*trace(dE)*I + 2*mu*dE\n",
-    "                dP = dF*S + F*dS\n",
-    "                for q = 1:nnodes\n",
-    "                    for j = 1:dim\n",
-    "                        K[dim*(p-1)+i,dim*(q-1)+j] += w*(dP[j,:]*dbasisdX[q,:]')[1]*detJ\n",
-    "                    end\n",
-    "                end\n",
-    "            end\n",
-    "        end\n",
-    "\n",
+    "        F = I + gradu\n",
+    "        E = 1/2*(gradu' + gradu + gradu'*gradu)\n",
+    "        # constitutive equation\n",
+    "        lambda = JuliaFEM.interpolate(el, \"lambda\", xi)\n",
+    "        mu = JuliaFEM.interpolate(el, \"mu\", xi)\n",
+    "        S = lambda*trace(E)*I + 2*mu*E\n",
+    "        P = F*S\n",
+    "        T = P*dbasisdX'\n",
+    "        return T\n",
     "    end\n",
+    "\n",
+    "    u = el.attributes[\"displacement\"]\n",
+    "    R(u) = T(u)\n",
+    "    return R(u) # we skip calculating external force vector for now.\n",
     "end\n",
     "\n",
-    "\"\"\"\n",
-    "Autodiff version.\n",
-    "\"\"\"\n",
-    "function calc_local_matrices!(X, u, R, K, basis, dbasis, lambda_, mu_, ipoints, iweights)\n",
-    "    dim, nnodes = size(X)\n",
-    "    I = eye(dim)\n",
-    "    R[:,:] = 0.0\n",
-    "\n",
-    "    #dF = zeros(dim, dim)\n",
-    "\n",
-    "    function calc_R!(u, R)\n",
-    "        for m = 1:length(iweights)\n",
-    "            w = iweights[m]\n",
-    "            xi = ipoints[m, :]\n",
-    "            # calculate material parameters\n",
-    "            lambda = typeof(lambda_) == Float64 ? lambda_ : dot(lambda_, basis(xi))\n",
-    "            mu = typeof(mu_) == Float64 ? mu_ : dot(mu_, basis(xi))\n",
-    "            Jt = X*dbasis(xi)\n",
-    "            detJ = det(Jt)\n",
-    "            dbasisdX = dbasis(xi)*inv(Jt)\n",
-    "\n",
-    "            gradu = u*dbasisdX\n",
-    "            F = I + gradu  # Deformation gradient\n",
-    "            E = 1/2*(gradu' + gradu + gradu'*gradu)  # Green-Lagrange strain tensor\n",
-    "            S = lambda*trace(E)*I + 2*mu*E  # PK2 stress tensor\n",
-    "            P = F*S  # PK1 stress tensor\n",
-    "\n",
-    "            R[:,:] += w*P*dbasisdX'*detJ\n",
-    "        end\n",
-    "    end\n",
-    "\n",
-    "    # herlper for tangent stiffness matrix\n",
-    "    function R!(u, R)\n",
-    "        R[:] = 0\n",
-    "        calc_R!(reshape(u, dim, nnodes), reshape(R, dim, nnodes))\n",
-    "        #calc_Wext!(reshape(u, 2, 4), reshape(R, 2, 4))\n",
-    "    end\n",
-    "    Jacobian = ForwardDiff.forwarddiff_jacobian(R!, Float64, fadtype=:dual, n=dim*nnodes, m=dim*nnodes)\n",
-    "\n",
-    "    K[:, :] = Jacobian(reshape(u, dim*nnodes))\n",
-    "    R!(reshape(u, dim*nnodes), reshape(R, dim*nnodes))\n",
-    "\n",
-    "end"
+    "calc_residual_vector = JuliaFEM.integrate(calc_residual_vector_integrand)\n",
+    "# calculate partial derivatives of R with respect to field \"displacement\"\n",
+    "calc_tangent_stiffness = JuliaFEM.linearize(calc_residual_vector, \"displacement\")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "That was our geometrically nonlinear elasticity solver. Note how we used automatic differentiation to linearize residual vector.\n",
+    "\n",
     "*Design principle 6*: we test our code. We use FactCheck for testing."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 4,
    "metadata": {
     "collapsed": true
    },
@@ -250,7 +141,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 6,
    "metadata": {
     "collapsed": false,
     "scrolled": false
@@ -267,62 +158,23 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "13-Aug 00:37:26:DEBUG:root:Converged in 6 iterations.\n",
-      "13-Aug 00:37:26:DEBUG:root:solution vector: \n",
+      "17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n",
+      "17-Aug 18:21:25:DEBUG:root:solution vector: \n",
       " [0.0 -0.3991450609547433 -0.07228582695592461 0.0\n",
       " 0.0 -2.1779892317073504 -2.222244754401764 0.0]\n",
-      "13-Aug 00:37:26:DEBUG:root:norm of u: 3.1292483947150047\n",
-      "13-Aug 00:37:27:DEBUG:root:Converged in 6 iterations.\n",
-      "13-Aug 00:37:27:DEBUG:root:solution vector: \n",
+      "17-Aug 18:21:25:DEBUG:root:norm of u: 3.1292483947150047\n",
+      "17-Aug 18:21:25:DEBUG:root:Converged in 6 iterations.\n",
+      "17-Aug 18:21:25:DEBUG:root:solution vector: \n",
       " [0.0 0.7433248532717796 1.048521014723486 0.0\n",
       " 0.0 -2.085766534304891 -1.9606633242166027 0.0]\n",
-      "13-Aug 00:37:27:DEBUG:root:norm of u: 3.1292483947150056\n",
-      "13-Aug 00:37:27:DEBUG:root:Iteration 1\n"
+      "17-Aug 18:21:26:DEBUG:root:norm of u: 3.1292483947150056\n"
      ]
     },
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Array(Float64,(12,12)) 12x12 Array{Float64,2}"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "13-Aug 00:37:27:DEBUG:root:Norm of du: 0.5992228342549063\n",
-      "13-Aug 00:37:27:DEBUG:root:Element displacement: [0.0 -0.02264423092574128 0.022536491965822688 0.0\n",
-      " 0.0 -0.12688379170176511 -0.12679760053383046 0.0]\n",
-      "13-Aug 00:37:27:DEBUG:root:Element displacement: [-0.02264423092574128 -0.029998807330416523 0.030242156525002267 0.022536491965822688\n",
-      " -0.12688379170176511 -0.4041002710283739 -0.40446733271991214 -0.12679760053383046]\n",
-      "13-Aug 00:37:27:DEBUG:root:solution vector: \n",
-      " [-0.02264423092574128 -0.029998807330416523 0.022536491965822688 0.030242156525002267 0.0 0.0\n",
-      " -0.12688379170176511 -0.4041002710283739 -0.12679760053383046 -0.40446733271991214 0.0 0.0]\n",
-      "13-Aug 00:37:27:DEBUG:root:norm of u: 0.5992228342549063\n"
-     ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      ":\n",
-      "  132.8     0.0   23.6    -3.0  -113.6  …   23.6     3.0  -33.2    15.0\n",
-      "    0.0   324.8    3.0    77.6     0.0      -3.0    77.6   15.0   -81.2\n",
-      "   23.6     3.0   66.4   -15.0   -33.2       0.0     0.0    0.0     0.0\n",
-      "   -3.0    77.6  -15.0   162.4    15.0       0.0     0.0    0.0     0.0\n",
-      " -113.6     0.0  -33.2    15.0   132.8     -33.2   -15.0   23.6    -3.0\n",
-      "    0.0  -317.6   15.0   -81.2     0.0  …  -15.0   -81.2    3.0    77.6\n",
-      "  -33.2   -15.0  -56.8     3.0    23.6       0.0     0.0    0.0     0.0\n",
-      "  -15.0   -81.2   -3.0  -158.8     3.0       0.0     0.0    0.0     0.0\n",
-      "   23.6    -3.0    0.0     0.0   -33.2      66.4    15.0  -56.8     3.0\n",
-      "    3.0    77.6    0.0     0.0   -15.0      15.0   162.4   -3.0  -158.8\n",
-      "  -33.2    15.0    0.0     0.0    23.6  …  -56.8    -3.0   66.4   -15.0\n",
-      "   15.0   -81.2    0.0     0.0    -3.0       3.0  -158.8  -15.0   162.4\n",
-      "K norm = 708.0378644377365\n",
-      "du = [-0.02264423092574128 -0.029998807330416523 0.022536491965822688 0.030242156525002267 0.0 0.0\n",
-      " -0.12688379170176511 -0.4041002710283739 -0.12679760053383046 -0.40446733271991214 0.0 0.0]\n",
-      "Out of 3 total facts:"
+      "2 facts verified.\n"
      ]
     },
     {
@@ -331,47 +183,52 @@
        "delayed_handler (generic function with 4 methods)"
       ]
      },
-     "execution_count": 4,
+     "execution_count": 6,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "facts(\"test solve one element model\") do\n",
-    "    X = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
-    "    F = [0 0; 0 0; 0 -2; 0 0]'\n",
-    "\n",
-    "    # Material properties\n",
-    "    E = 90\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",
-    "\n",
-    "    u = zeros(2, 4)\n",
-    "    du = zeros(2, 4)\n",
-    "    R = zeros(2, 4)\n",
-    "    K = zeros(8, 8)\n",
-    "\n",
+    "function get_test_element()\n",
+    "    # set up one linear quadrangle element\n",
     "    basis(xi) = [\n",
     "        (1-xi[1])*(1-xi[2])/4\n",
     "        (1+xi[1])*(1-xi[2])/4\n",
     "        (1+xi[1])*(1+xi[2])/4\n",
     "        (1-xi[1])*(1+xi[2])/4]\n",
-    "\n",
     "    dbasis(xi) = [-(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",
+    "    ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]'\n",
+    "    iweights = [1.0, 1.0, 1.0, 1.0]\n",
+    "    attributes = Dict()\n",
+    "    e = JuliaFEM.Element(1, [1, 2, 3, 4], basis, dbasis, attributes, ipoints, iweights)\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",
+    "    e.attributes[\"coordinates\"] = [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]'\n",
+    "    e.attributes[\"lambda\"] = la\n",
+    "    e.attributes[\"mu\"] = mu\n",
+    "    e.attributes[\"displacement\"] = [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]'\n",
+    "    return e\n",
+    "end\n",
+    "\n",
+    "facts(\"test solve one element model\") do\n",
+    "\n",
+    "    e = get_test_element()\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",
     "\n",
-    "    ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
-    "    iweights = [1, 1, 1, 1]\n",
     "    free_dofs = [3, 4, 5, 6]\n",
-    "\n",
     "    for i=1:10\n",
-    "        calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
+    "        R = calc_residual_vector(e)\n",
+    "        K = calc_tangent_stiffness(e)\n",
     "        du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
-    "        u += du\n",
+    "        e.attributes[\"displacement\"] += du\n",
     "        if norm(du) < 1.0e-9\n",
     "            Logging.debug(\"Converged in $i iterations.\")\n",
     "            break\n",
@@ -379,73 +236,36 @@
     "    end\n",
     "\n",
     "    # Tested against Elmer solution\n",
+    "    u = e.attributes[\"displacement\"]\n",
     "    Logging.debug(\"solution vector: \\n $u\")\n",
     "    @fact u[2, 3] --> roughly(-2.222244754401764)\n",
     "    norm1 = norm(u)\n",
     "    Logging.debug(\"norm of u: $(norm(u))\")\n",
     "\n",
-    "    # We rotate model a bit and make sure that L2 norm is same\n",
+    "    # We rotate model a bit and make sure that norm remains same\n",
     "    phi = 30/180*pi\n",
     "    rmat = [\n",
     "        cos(phi) -sin(phi)\n",
     "        sin(phi)  cos(phi)]\n",
-    "    X = rmat*X\n",
+    "    e.attributes[\"coordinates\"] = rmat*e.attributes[\"coordinates\"]\n",
     "    F = rmat*F\n",
-    "    u = zeros(2, 4)\n",
+    "\n",
+    "    e.attributes[\"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",
-    "        calc_local_matrices!(X, u, R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
+    "        R = calc_residual_vector(e)\n",
+    "        K = calc_tangent_stiffness(e)\n",
     "        du[free_dofs] = K[free_dofs, free_dofs] \\ -(R - F)[free_dofs]\n",
-    "        u += du\n",
+    "        e.attributes[\"displacement\"] += du\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",
     "    Logging.debug(\"solution vector: \\n $u\")\n",
     "    Logging.debug(\"norm of u: $(norm(u))\")\n",
     "    @fact norm(u) --> roughly(norm1) \n",
-    "\n",
-    "    # test two element model\n",
-    "    X = [0.0 0.0; 5.0 0.0; 5.0 1.0; 0.0 1.0]'\n",
-    "    u = zeros(2, 6)\n",
-    "    du = zeros(2, 6)\n",
-    "    R = zeros(2, 4)\n",
-    "    K = zeros(8, 8)\n",
-    "    ass1 = [9, 10, 1, 2, 5, 6, 11, 12]\n",
-    "    ass2 = [1, 2, 3, 4, 7, 8, 5, 6]\n",
-    "    free_dofs = collect(1:8)\n",
-    "    F = [0 0; 0 0; 0 0; 0 -0.1; 0 0; 0 0]'\n",
-    "\n",
-    "    A = zeros(12, 12)\n",
-    "    b = zeros(2, 6)\n",
-    "    for i=1:1\n",
-    "        Logging.debug(\"Iteration $i\")\n",
-    "        A[:,:] = 0.0\n",
-    "        b[:] = 0.0\n",
-    "        #Logging.debug(\"Assembling\")\n",
-    "        for ass in (ass1, ass2)\n",
-    "            #Logging.debug(\"ass = $ass, u[ass] = $(u[ass])\")\n",
-    "            calc_local_matrices!(X, u[ass], R, K, basis, dbasis, la, mu, ipoints, iweights)\n",
-    "            A[ass,ass] += K\n",
-    "            b[ass] += R[:]\n",
-    "        end\n",
-    "        dump(round(A, 2))\n",
-    "        println(\"K norm = $(norm(A[free_dofs, free_dofs]))\")\n",
-    "        du[free_dofs] = A[free_dofs, free_dofs] \\ -(b - F)[free_dofs]\n",
-    "        println(\"du = $du\")\n",
-    "        u += du\n",
-    "        Logging.debug(\"Norm of du: $(norm(du))\")\n",
-    "        for ass in (ass1, ass2)\n",
-    "            Logging.debug(\"Element displacement: $(reshape(u[ass], 2, 4))\")\n",
-    "        end\n",
-    "        if norm(du) < 1.0e-9\n",
-    "            Logging.debug(\"Converged in $i iterations.\")\n",
-    "            break\n",
-    "        end\n",
-    "    end\n",
-    "    Logging.debug(\"solution vector: \\n $u\")\n",
-    "    Logging.debug(\"norm of u: $(norm(u))\")\n",
-    "    @pending norm(u) --> :something\n",
     "end"
    ]
   },
@@ -456,160 +276,6 @@
     "One element solutions are not particularly interesting so next step is to create function that assembles global matrix from local matrices. Some data types:"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "#type Node\n",
-    "#    id :: Int\n",
-    "#    #elements :: Array{Int64, 1}\n",
-    "#end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 6,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [],
-   "source": [
-    "type Element\n",
-    "    id :: Int\n",
-    "    node_ids :: Array{Int64, 1}\n",
-    "    coordinates :: Array{Float64, 2}\n",
-    "    attributes :: Dict{ASCIIString, Any}\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 7,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "\n",
-      "  Verified: 2\n",
-      "  Pending:  1\n"
-     ]
-    }
-   ],
-   "source": [
-    "type Assembly\n",
-    "    # LHS\n",
-    "    I :: Array{Int64, 1}\n",
-    "    J :: Array{Int64, 1}\n",
-    "    A :: Array{Float64, 1}\n",
-    "    # RHS\n",
-    "    i :: Array{Int64, 1}\n",
-    "    b :: Array{Float64, 1}\n",
-    "    # global dofs for each element\n",
-    "    gdofs :: Dict{Int64, Array{Int64, 1}}\n",
-    "end"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 8,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "get_integration_scheme (generic function with 2 methods)"
-      ]
-     },
-     "execution_count": 8,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
-   "source": [
-    "\"\"\"\n",
-    "Return shape functions and their derivatives for a element.\n",
-    "\n",
-    "Parameters\n",
-    "----------\n",
-    "element::Element\n",
-    "\n",
-    "Returns\n",
-    "-------\n",
-    "tuple (basis, dbasis)\n",
-    "\"\"\"\n",
-    "function get_shape_functions(el::Element)\n",
-    "    ndim, nnodes = size(el.coordinates)\n",
-    "    #Logging.debug(\"ndim = $ndim, nnodes=$nnodes\")\n",
-    "    if (nnodes == 4) & (ndim == 2)\n",
-    "        basis(xi) = [\n",
-    "            (1-xi[1])*(1-xi[2])/4\n",
-    "            (1+xi[1])*(1-xi[2])/4\n",
-    "            (1+xi[1])*(1+xi[2])/4\n",
-    "            (1-xi[1])*(1+xi[2])/4]\n",
-    "        dbasis(xi) = [-(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",
-    "        return basis, dbasis\n",
-    "    elseif (nnodes == 10) & (ndim == 3)\n",
-    "        basis(xi) = [\n",
-    "                (xi[1] + xi[2] + xi[3] - 1)*(2*xi[1] + 2*xi[2] + 2*xi[3] - 1)\n",
-    "                                                          xi[1]*(2*xi[1] - 1)\n",
-    "                                                          xi[2]*(2*xi[2] - 1)\n",
-    "                                                          xi[3]*(2*xi[3] - 1)\n",
-    "                                         -4*xi[1]*(xi[1] + xi[2] + xi[3] - 1)\n",
-    "                                                                4*xi[1]*xi[2]\n",
-    "                                         -4*xi[2]*(xi[1] + xi[2] + xi[3] - 1)\n",
-    "                                         -4*xi[3]*(xi[1] + xi[2] + xi[3] - 1)\n",
-    "                                                                4*xi[1]*xi[3]\n",
-    "                                                                4*xi[2]*xi[3]\n",
-    "        ]\n",
-    "        dbasis(xi) = [\n",
-    "             4*xi[1] + 4*xi[2] + 4*xi[3] - 3  4*xi[1] + 4*xi[2] + 4*xi[3] - 3  4*xi[1] + 4*xi[2] + 4*xi[3] - 3\n",
-    "                                 4*xi[1] - 1                                0                                0\n",
-    "                                           0                      4*xi[2] - 1                                0\n",
-    "                                           0                                0                      4*xi[3] - 1\n",
-    "            -4*(2*xi[1] + xi[2] + xi[3] - 1)                         -4*xi[1]                         -4*xi[1]\n",
-    "                                     4*xi[2]                          4*xi[1]                                0\n",
-    "                                    -4*xi[2] -4*(xi[1] + 2*xi[2] + xi[3] - 1)                         -4*xi[2]\n",
-    "                                    -4*xi[3]                         -4*xi[3] -4*(xi[1] + xi[2] + 2*xi[3] - 1)\n",
-    "                                     4*xi[3]                                0                          4*xi[1]\n",
-    "                                           0                          4*xi[3]                          4*xi[2]\n",
-    "        ]\n",
-    "        return basis, dbasis\n",
-    "    end\n",
-    "    throw(\"Unknown function space, ndim=$ndim, nnodes=$nnodes\")\n",
-    "end\n",
-    "\n",
-    "\"\"\"\n",
-    "\"\"\"\n",
-    "function get_integration_scheme(el::Element, order=2)\n",
-    "    ndim, nnodes = size(el.coordinates)\n",
-    "    if (nnodes == 4) & (order == 2) & (ndim == 2)\n",
-    "        ipoints = 1/sqrt(3)*[-1 -1; 1 -1; 1 1; -1 1]\n",
-    "        iweights = [1, 1, 1, 1]\n",
-    "        return ipoints, iweights\n",
-    "    elseif (nnodes == 10) & (order == 2) & (ndim == 3) # c3d10\n",
-    "        # from code aster documentation\n",
-    "        a = 1/20*(5-sqrt(5))\n",
-    "        b = 1/20*(5+3*sqrt(5))\n",
-    "        ipoints = [a a a; a a b; a b a; b a a]\n",
-    "        iweights = 1/24*[1 1 1 1]\n",
-    "        return ipoints, iweights\n",
-    "    end\n",
-    "end"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 9,
@@ -629,7 +295,7 @@
     }
    ],
    "source": [
-    "function assemble_element!(ass::Assembly, el::Element, io=2)\n",
+    "function assemble_element!(ass::Assembly, el::Element)\n",
     "\n",
     "    # Material properties\n",
     "    E = el.attributes[\"Young\"]\n",
@@ -699,28 +365,28 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "13-Aug 00:37:30:DEBUG:root:Adding nodes to array\n",
-      "13-Aug 00:37:30:DEBUG:root:Creating elements\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 1\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 3.090022136728999\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 2\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 0.32121316021535135\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 3\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 0.040431781939994194\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 4\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 0.0009291101052124042\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 5\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 1.5638899022213175e-7\n",
-      "13-Aug 00:37:30:DEBUG:root:Starting iteration 6\n",
-      "13-Aug 00:37:30:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:30:DEBUG:root:Solution norm = 1.0118539067290854e-14\n",
-      "13-Aug 00:37:30:DEBUG:root:Converged in 6 iterations.\n",
-      "13-Aug 00:37:30:DEBUG:root:Displacement of element = \n",
+      "15-Aug 17:03:55:DEBUG:root:Adding nodes to array\n",
+      "15-Aug 17:03:55:DEBUG:root:Creating elements\n",
+      "15-Aug 17:03:55:DEBUG:root:Starting iteration 1\n",
+      "15-Aug 17:03:55:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:57:DEBUG:root:Solution norm = 3.090022136728999\n",
+      "15-Aug 17:03:57:DEBUG:root:Starting iteration 2\n",
+      "15-Aug 17:03:57:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.32121316021535135\n",
+      "15-Aug 17:03:57:DEBUG:root:Starting iteration 3\n",
+      "15-Aug 17:03:57:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.040431781939994194\n",
+      "15-Aug 17:03:57:DEBUG:root:Starting iteration 4\n",
+      "15-Aug 17:03:57:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:57:DEBUG:root:Solution norm = 0.0009291101052124042\n",
+      "15-Aug 17:03:57:DEBUG:root:Starting iteration 5\n",
+      "15-Aug 17:03:57:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:57:DEBUG:root:Solution norm = 1.5638899022213175e-7\n",
+      "15-Aug 17:03:58:DEBUG:root:Starting iteration 6\n",
+      "15-Aug 17:03:58:DEBUG:root:Assembling\n",
+      "15-Aug 17:03:58:DEBUG:root:Solution norm = 1.0118539067290854e-14\n",
+      "15-Aug 17:03:58:DEBUG:root:Converged in 6 iterations.\n",
+      "15-Aug 17:03:58:DEBUG:root:Displacement of element = \n",
       "[-0.39914506095474334 -0.0722858269559246 0.0 0.0\n",
       " -2.1779892317073504 -2.2222447544017645 0.0 0.0]\n"
      ]
@@ -844,48 +510,48 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "13-Aug 00:37:31:DEBUG:root:Creating elements\n",
-      "13-Aug 00:37:31:INFO:root:create_ldof2gdofmap: dofs per node: 2\n",
-      "13-Aug 00:37:31:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],5=>[9,10],6=>[11,12],1=>[1,2])\n",
-      "13-Aug 00:37:32:INFO:root:solve!: dofs per node: 2\n",
-      "13-Aug 00:37:32:DEBUG:root:Problem size = 12\n",
-      "13-Aug 00:37:32:DEBUG:root:Starting iteration 1\n",
-      "13-Aug 00:37:32:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:32:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:32:DEBUG:root:Added 5 Lagrange multipliers to model\n",
-      "13-Aug 00:37:32:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:32:DEBUG:root:Solving system of equations. Total size = 16\n",
-      "13-Aug 00:37:32:DEBUG:root:Solution norm du = 0.6015838690633517\n",
-      "13-Aug 00:37:32:DEBUG:root:Starting iteration 2\n",
-      "13-Aug 00:37:32:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:32:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:32:DEBUG:root:Added 5 Lagrange multipliers to model\n",
-      "13-Aug 00:37:32:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:32:DEBUG:root:Solving system of equations. Total size = 16\n",
-      "13-Aug 00:37:32:DEBUG:root:Solution norm du = 0.013420417380980414\n",
-      "13-Aug 00:37:33:DEBUG:root:Starting iteration 3\n",
-      "13-Aug 00:37:33:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
-      "13-Aug 00:37:33:DEBUG:root:Solution norm du = 0.00032202957936854873\n",
-      "13-Aug 00:37:33:DEBUG:root:Starting iteration 4\n",
-      "13-Aug 00:37:33:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
-      "13-Aug 00:37:33:DEBUG:root:Solution norm du = 9.906677094801476e-8\n",
-      "13-Aug 00:37:33:DEBUG:root:Starting iteration 5\n",
-      "13-Aug 00:37:33:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:33:DEBUG:root:Added 5 Lagrange multipliers to model\n",
-      "13-Aug 00:37:33:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:33:DEBUG:root:Solving system of equations. Total size = 16\n",
-      "13-Aug 00:37:33:DEBUG:root:Solution norm du = 5.027624635820698e-15\n",
-      "13-Aug 00:37:33:DEBUG:root:Converged in 5 iterations.\n",
-      "13-Aug 00:37:33:DEBUG:root:Displacement of element = \n",
+      "15-Aug 17:03:59:DEBUG:root:Creating elements\n",
+      "15-Aug 17:03:59:INFO:root:create_ldof2gdofmap: dofs per node: 2\n",
+      "15-Aug 17:04:00:DEBUG:root:Dict(4=>[7,8],2=>[3,4],3=>[5,6],5=>[9,10],6=>[11,12],1=>[1,2])\n",
+      "15-Aug 17:04:04:INFO:root:solve!: dofs per node: 2\n",
+      "15-Aug 17:04:04:DEBUG:root:Problem size = 12\n",
+      "15-Aug 17:04:04:DEBUG:root:Starting iteration 1\n",
+      "15-Aug 17:04:04:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n",
+      "15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.6015838690633517\n",
+      "15-Aug 17:04:05:DEBUG:root:Starting iteration 2\n",
+      "15-Aug 17:04:05:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:05:DEBUG:root:Solving system of equations. Total size = 16\n",
+      "15-Aug 17:04:05:DEBUG:root:Solution norm du = 0.013420417380980414\n",
+      "15-Aug 17:04:05:DEBUG:root:Starting iteration 3\n",
+      "15-Aug 17:04:05:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:05:DEBUG:root:Added 5 Lagrange multipliers to model\n",
+      "15-Aug 17:04:05:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
+      "15-Aug 17:04:06:DEBUG:root:Solution norm du = 0.00032202957936854873\n",
+      "15-Aug 17:04:06:DEBUG:root:Starting iteration 4\n",
+      "15-Aug 17:04:06:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n",
+      "15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
+      "15-Aug 17:04:06:DEBUG:root:Solution norm du = 9.906677094801476e-8\n",
+      "15-Aug 17:04:06:DEBUG:root:Starting iteration 5\n",
+      "15-Aug 17:04:06:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:06:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:06:DEBUG:root:Added 5 Lagrange multipliers to model\n",
+      "15-Aug 17:04:06:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:06:DEBUG:root:Solving system of equations. Total size = 16\n",
+      "15-Aug 17:04:06:DEBUG:root:Solution norm du = 5.027624635820698e-15\n",
+      "15-Aug 17:04:06:DEBUG:root:Converged in 5 iterations.\n",
+      "15-Aug 17:04:06:DEBUG:root:Displacement of element = \n",
       "[-0.02442313597467864 -0.039356000063335075 0.021097993207232584 0.020877031423993653\n",
       " -0.12673626841485705 -0.40433021969759375 -0.40656320177872923 -0.1275776940913048]\n"
      ]
@@ -1136,7 +802,7 @@
       "\n",
       "WARNING: deprecated syntax \"{a=>b, ...}\" at /Users/jukka/.julia/v0.4/JuliaFEM/src/abaqus_reader.jl:32.\n",
       "Use \"Dict{Any,Any}(a=>b, ...)\" instead.\n",
-      "13-Aug 00:37:37:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n",
+      "15-Aug 17:04:20:INFO:root:Registered handlers: Any[\"ELEMENT\",\"NODE\",\"NSET\"]\n",
       "WARNING: beginswith is deprecated, use startswith instead.\n",
       " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
       " in beginswith at deprecated.jl:30\n",
@@ -1155,8 +821,8 @@
       " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
       " in anonymous at task.jl:365\n",
       "while loading In[13], in expression starting on line 3\n",
-      "13-Aug 00:37:38:DEBUG:root:Found NODE section\n",
-      "13-Aug 00:37:38:DEBUG:root:Found ELEMENT section\n",
+      "15-Aug 17:04:21:DEBUG:root:Found NODE section\n",
+      "15-Aug 17:04:22:DEBUG:root:Found ELEMENT section\n",
       "WARNING: integer(s::AbstractString) is deprecated, use parse(Int,s) instead.\n",
       " in depwarn at /Applications/Julia-0.4.0-dev-539c818c4e.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
       " in integer at deprecated.jl:49\n",
@@ -1169,14 +835,14 @@
       " in eventloop at /Users/jukka/.julia/v0.4/IJulia/src/IJulia.jl:123\n",
       " in anonymous at task.jl:365\n",
       "while loading In[13], in expression starting on line 3\n",
-      "13-Aug 00:37:39:DEBUG:root:120 elements found\n",
-      "13-Aug 00:37:40:INFO:root:Creating ELSET Body1\n",
-      "13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
-      "13-Aug 00:37:40:DEBUG:root:Creating node set SUPPORT\n",
-      "13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
-      "13-Aug 00:37:40:DEBUG:root:Creating node set LOAD\n",
-      "13-Aug 00:37:40:DEBUG:root:Found NSET section\n",
-      "13-Aug 00:37:40:DEBUG:root:Creating node set TOP\n"
+      "15-Aug 17:04:27:DEBUG:root:120 elements found\n",
+      "15-Aug 17:04:27:INFO:root:Creating ELSET Body1\n",
+      "15-Aug 17:04:27:DEBUG:root:Found NSET section\n",
+      "15-Aug 17:04:28:DEBUG:root:Creating node set SUPPORT\n",
+      "15-Aug 17:04:28:DEBUG:root:Found NSET section\n",
+      "15-Aug 17:04:28:DEBUG:root:Creating node set LOAD\n",
+      "15-Aug 17:04:28:DEBUG:root:Found NSET section\n",
+      "15-Aug 17:04:28:DEBUG:root:Creating node set TOP\n"
      ]
     },
     {
@@ -1214,47 +880,47 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "13-Aug 00:37:41:DEBUG:root:Creating elements\n",
-      "13-Aug 00:37:41:INFO:root:create_ldof2gdofmap: dofs per node: 3\n",
-      "13-Aug 00:37:41:INFO:root:solve!: dofs per node: 3\n",
-      "13-Aug 00:37:41:DEBUG:root:Problem size = 894\n",
-      "13-Aug 00:37:41:DEBUG:root:Starting iteration 1\n",
-      "13-Aug 00:37:41:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:43:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:43:DEBUG:root:Added 28 Lagrange multipliers to model\n",
-      "13-Aug 00:37:43:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:43:DEBUG:root:Solving system of equations. Total size = 921\n",
-      "13-Aug 00:37:43:DEBUG:root:Solution norm du = 54.19642700242575\n",
-      "13-Aug 00:37:43:DEBUG:root:Starting iteration 2\n",
-      "13-Aug 00:37:43:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:45:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:45:DEBUG:root:Added 28 Lagrange multipliers to model\n",
-      "13-Aug 00:37:45:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:45:DEBUG:root:Solving system of equations. Total size = 921\n",
-      "13-Aug 00:37:45:DEBUG:root:Solution norm du = 1.68729144400063\n",
-      "13-Aug 00:37:45:DEBUG:root:Starting iteration 3\n",
-      "13-Aug 00:37:45:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:47:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:47:DEBUG:root:Added 28 Lagrange multipliers to model\n",
-      "13-Aug 00:37:47:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:47:DEBUG:root:Solving system of equations. Total size = 921\n",
-      "13-Aug 00:37:47:DEBUG:root:Solution norm du = 0.04175080278208098\n",
-      "13-Aug 00:37:47:DEBUG:root:Starting iteration 4\n",
-      "13-Aug 00:37:47:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:48:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:48:DEBUG:root:Added 28 Lagrange multipliers to model\n",
-      "13-Aug 00:37:48:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:48:DEBUG:root:Solving system of equations. Total size = 921\n",
-      "13-Aug 00:37:49:DEBUG:root:Solution norm du = 2.4176157836844963e-5\n",
-      "13-Aug 00:37:49:DEBUG:root:Starting iteration 5\n",
-      "13-Aug 00:37:49:DEBUG:root:Assembling\n",
-      "13-Aug 00:37:50:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
-      "13-Aug 00:37:50:DEBUG:root:Added 28 Lagrange multipliers to model\n",
-      "13-Aug 00:37:50:DEBUG:root:Adding Neumann boundary conditions\n",
-      "13-Aug 00:37:50:DEBUG:root:Solving system of equations. Total size = 921\n",
-      "13-Aug 00:37:50:DEBUG:root:Solution norm du = 1.4448231751500831e-11\n",
-      "13-Aug 00:37:50:DEBUG:root:Converged in 5 iterations.\n",
-      "13-Aug 00:37:50:INFO:root:Maximum absolute displacement in y direction: 5.245400568184194\n"
+      "15-Aug 17:04:30:DEBUG:root:Creating elements\n",
+      "15-Aug 17:04:31:INFO:root:create_ldof2gdofmap: dofs per node: 3\n",
+      "15-Aug 17:04:31:INFO:root:solve!: dofs per node: 3\n",
+      "15-Aug 17:04:31:DEBUG:root:Problem size = 894\n",
+      "15-Aug 17:04:31:DEBUG:root:Starting iteration 1\n",
+      "15-Aug 17:04:31:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:36:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:36:DEBUG:root:Added 28 Lagrange multipliers to model\n",
+      "15-Aug 17:04:36:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:36:DEBUG:root:Solving system of equations. Total size = 921\n",
+      "15-Aug 17:04:36:DEBUG:root:Solution norm du = 72.87727091053921\n",
+      "15-Aug 17:04:36:DEBUG:root:Starting iteration 2\n",
+      "15-Aug 17:04:36:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:40:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:40:DEBUG:root:Added 28 Lagrange multipliers to model\n",
+      "15-Aug 17:04:40:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:40:DEBUG:root:Solving system of equations. Total size = 921\n",
+      "15-Aug 17:04:40:DEBUG:root:Solution norm du = 2.737285770100854\n",
+      "15-Aug 17:04:40:DEBUG:root:Starting iteration 3\n",
+      "15-Aug 17:04:40:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:44:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:44:DEBUG:root:Added 28 Lagrange multipliers to model\n",
+      "15-Aug 17:04:44:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:44:DEBUG:root:Solving system of equations. Total size = 921\n",
+      "15-Aug 17:04:44:DEBUG:root:Solution norm du = 0.07997112801214978\n",
+      "15-Aug 17:04:44:DEBUG:root:Starting iteration 4\n",
+      "15-Aug 17:04:44:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:48:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:48:DEBUG:root:Added 28 Lagrange multipliers to model\n",
+      "15-Aug 17:04:48:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:48:DEBUG:root:Solving system of equations. Total size = 921\n",
+      "15-Aug 17:04:48:DEBUG:root:Solution norm du = 6.406557430235748e-5\n",
+      "15-Aug 17:04:48:DEBUG:root:Starting iteration 5\n",
+      "15-Aug 17:04:48:DEBUG:root:Assembling\n",
+      "15-Aug 17:04:52:DEBUG:root:Adding Dirichlet boundary conditions using Lagrange multipliers\n",
+      "15-Aug 17:04:52:DEBUG:root:Added 28 Lagrange multipliers to model\n",
+      "15-Aug 17:04:52:DEBUG:root:Adding Neumann boundary conditions\n",
+      "15-Aug 17:04:52:DEBUG:root:Solving system of equations. Total size = 921\n",
+      "15-Aug 17:04:52:DEBUG:root:Solution norm du = 6.870072007793185e-11\n",
+      "15-Aug 17:04:52:DEBUG:root:Converged in 5 iterations.\n",
+      "15-Aug 17:04:52:INFO:root:Maximum absolute displacement in y direction: 6.9925206227884695\n"
      ]
     }
    ],
@@ -1367,8 +1033,8 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "13-Aug 00:37:51:INFO:root:Number of nodes in model: 298\n",
-      "13-Aug 00:37:51:INFO:root:Number of elements in model: 120\n"
+      "15-Aug 17:04:54:INFO:root:Number of nodes in model: 298\n",
+      "15-Aug 17:04:54:INFO:root:Number of elements in model: 120\n"
      ]
     }
    ],
@@ -1498,31 +1164,31 @@
      "data": {
       "text/plain": [
        "Dict{Any,Any} with 298 entries:\n",
-       "  288 => [-0.11672879201601734,-4.945796660856531,0.34641446977739976]\n",
-       "  11  => [-0.10351195894967631,-4.756564623966246,0.3419475317395125]\n",
-       "  158 => [-0.004781967238459207,-0.15800625216143013,-0.03641750463465448]\n",
-       "  215 => [-0.009679264390442352,-3.1283115653077243,0.1060066992444004]\n",
-       "  134 => [-0.012153531002752209,-2.568061461984137,0.1358780416251414]\n",
-       "  160 => [0.045841800099574996,-0.3952996715282528,-0.04337647899598485]\n",
-       "  29  => [-0.019610173840903154,-0.13106872577808154,-0.04187494739380516]\n",
-       "  131 => [0.04911611452200464,-2.5861023637514955,0.15290704689462717]\n",
-       "  249 => [-0.18589774363886435,-4.932931637644884,0.27789363131184]\n",
-       "  207 => [0.05588675730587902,-2.563471323701039,0.0669746466494972]\n",
-       "  173 => [-0.017826777186855432,-3.4323824666471117,0.08923268461041232]\n",
-       "  289 => [-0.08168560138254058,-4.782381469728047,0.4566974667444831]\n",
-       "  74  => [-0.013759506382279562,-2.9054222587693794,0.1389413161160111]\n",
-       "  201 => [0.006183182231407451,-2.8511450493652664,0.14659955727801022]\n",
-       "  176 => [0.04229935999959756,-2.349387641785469,0.09405299762367403]\n",
-       "  57  => [-0.0022029054234204426,-3.062804411311773,0.14846103116563106]\n",
-       "  31  => [0.005313697788644969,-0.25237461591643895,-0.030307540603490793]\n",
-       "  285 => [-0.10859185872935265,-4.832314881155342,0.35832378735042353]\n",
-       "  70  => [0.042811232565748765,-2.7950866978216338,0.1421452814276251]\n",
-       "  33  => [-0.016450709743293015,-0.379272230318473,0.027549328283933208]\n",
-       "  252 => [-0.14470395638682446,-0.9284887677850233,0.043594035381458736]\n",
-       "  114 => [-0.14842327883584733,-0.3532980011491899,-0.018641651848308006]\n",
-       "  165 => [-0.04346921888336824,-3.2745129973662452,0.07316785266950793]\n",
-       "  96  => [0.1501458469524447,-3.8471140528652787,0.1489845291162686]\n",
-       "  133 => [0.09045793517058681,-2.3950191788705957,0.12799724500423504]\n",
+       "  288 => [-0.21784897897079505,-6.64256844921929,0.18339325827190953]\n",
+       "  11  => [-0.17342410820059667,-6.343125806991628,0.17544367431766225]\n",
+       "  158 => [-0.0007723334905839533,-0.2321904799732306,-0.03883691791045858]\n",
+       "  215 => [-0.03299227384130259,-4.2236276356937505,-0.0017757488281830807]\n",
+       "  134 => [-0.06336705102703495,-3.3332586850719337,0.06555932581531249]\n",
+       "  160 => [0.03205303832663507,-0.536361456950779,-0.07291283371223893]\n",
+       "  29  => [-0.01620524134284992,-0.19455583841608232,-0.03838421619811468]\n",
+       "  131 => [0.04015396763870081,-3.5101776818531523,0.10045744211179368]\n",
+       "  249 => [-0.2491004521409732,-6.570813533283954,0.06837890174506163]\n",
+       "  207 => [0.03974407267311282,-3.48588780058722,0.007971846403163977]\n",
+       "  173 => [-0.04472757590353555,-4.582185697159958,-0.04434348606400347]\n",
+       "  289 => [-0.16494613158592017,-6.360857515886982,0.33435404817411213]\n",
+       "  74  => [-0.03449870789750817,-3.947284664599636,0.0692691027146262]\n",
+       "  201 => [-0.00128347996939122,-3.8527236749372054,0.054216444201382656]\n",
+       "  176 => [-0.004729440266412542,-3.132546258410676,0.0036604388995965386]\n",
+       "  57  => [-0.021453552968802154,-4.134921994863808,0.04888988010136468]\n",
+       "  31  => [-0.005045439296021616,-0.3424143474221696,-0.0365148544462663]\n",
+       "  285 => [-0.19371941654927557,-6.446502168659889,0.18022053743770028]\n",
+       "  70  => [0.006277500845425116,-3.802694347992476,0.06251589929767469]\n",
+       "  33  => [-0.02778798411401373,-0.5091869603980016,0.03538255963543149]\n",
+       "  252 => [-0.11915625386204923,-1.2182594899567352,-0.025975085489623087]\n",
+       "  114 => [-0.14643901527350212,-0.4628626545501056,-0.0057156433415159625]\n",
+       "  165 => [-0.07462160313837857,-4.426674712264929,-0.04391022455570823]\n",
+       "  96  => [0.0541064545893336,-5.057801826440659,-0.01669155772495867]\n",
+       "  133 => [0.01615115174975327,-3.140610327091724,0.028951801068022566]\n",
        "  ⋮   => ⋮"
       ]
      },
@@ -1552,9 +1218,9 @@
      "data": {
       "text/plain": [
        "3x298 Array{Float64,2}:\n",
-       " -0.174259   -0.175081   0.0359566   0.0279021   0.00925611   0.0379208  …  -0.0248959  -0.116472     0.191156   0.193018  -0.0153311  -0.123282  0.0\n",
-       " -0.978534   -1.32327   -2.42646    -2.50175    -2.72095     -2.59357       -3.20696    -0.707364    -1.76025   -1.33439   -2.22191    -1.93398   0.0\n",
-       "  0.0535874   0.111758   0.125464    0.115461    0.145374     0.1251         0.0902774   0.00418633   0.368958   0.387435   0.185441    0.230546  0.0"
+       " -0.202475   -0.17437     0.0131267   -0.0340868  -0.0252162   0.0521758  …  -0.054781    -0.169903   0.186936   0.102666  -0.0970747  -0.164397  0.0\n",
+       " -1.34204    -1.66404    -3.20587     -3.27673    -3.68732    -3.36789       -4.33822     -1.03943   -2.27466   -1.67467   -2.91213    -2.52353   0.0\n",
+       "  0.0211856   0.0324162   0.00213878   0.0377156   0.0683502   0.030081      -0.00684333  -0.035509   0.304867   0.26033    0.113802    0.152433  0.0"
       ]
      },
      "execution_count": 20,
@@ -1602,7 +1268,7 @@
     {
      "data": {
       "text/plain": [
-       "28215"
+       "28330"
       ]
      },
      "execution_count": 22,
diff --git a/src/math.jl b/src/math.jl
index 5f0864b..b839cb6 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -75,7 +75,7 @@ Return partial derivatives of shape functions w.r.t X using chain rule.
 """
 function get_dbasisdX(el::Element, xi)
     J = interpolate(el, "coordinates", xi; derivative=true)
-    dbasisdX = el.dbasis(xi)*inv(J)
+    dbasisdX = el.dbasis(xi)*inv(J')
     return dbasisdX
 end