From 7403fb8cb72121b3fade9406bc4b2cbc4e593f51 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Tue, 29 Sep 2015 00:07:56 +0300
Subject: [PATCH] time properly implemented to .. everything.

---
 .../2015-08-29-developing-juliafem.ipynb      | 368 +++++++++++++-----
 src/JuliaFEM.jl                               |   8 +-
 src/elements.jl                               |  95 +++--
 src/equations.jl                              |  24 +-
 src/lagrange.jl                               |  86 ++--
 src/math.jl                                   |  86 +---
 src/types.jl                                  |  14 +-
 test/test_elements.jl                         |  12 +-
 test/test_types.jl                            |   4 +
 9 files changed, 410 insertions(+), 287 deletions(-)

diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb
index 258c2ae..6600a57 100644
--- a/notebooks/2015-08-29-developing-juliafem.ipynb
+++ b/notebooks/2015-08-29-developing-juliafem.ipynb
@@ -13,11 +13,35 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 1,
    "metadata": {
     "collapsed": false
    },
-   "outputs": [],
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n",
+      "WARNING: Base.String is deprecated, use AbstractString instead.\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "Logger(root,DEBUG,PipeEndpoint(open, 0 bytes waiting),root)"
+      ]
+     },
+     "execution_count": 1,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "using Logging\n",
     "using FactCheck\n",
@@ -49,7 +73,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": null,
+   "execution_count": 2,
    "metadata": {
     "collapsed": false
    },
@@ -133,7 +157,7 @@
     {
      "data": {
       "text/plain": [
-       "get_element_dimension (generic function with 8 methods)"
+       "get_element_dimension (generic function with 7 methods)"
       ]
      },
      "execution_count": 5,
@@ -164,21 +188,25 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "28-Sep 12:56:52:INFO:root:Testing element MyQuad4\n",
-      "28-Sep 12:56:52:INFO:root:number of basis functions in this element: 4\n",
-      "28-Sep 12:56:53:INFO:root:Initializing element\n",
-      "28-Sep 12:56:53:INFO:root:Element dimension: 2\n",
-      "28-Sep 12:56:53:INFO:root:Setting scalar field JuliaFEM.Field{Int64}(0.0,1,[1,2,3,4]) to element.\n"
+      "29-Sep 00:06:50:INFO:root:Testing element MyQuad4\n",
+      "29-Sep 00:06:50:INFO:root:number of basis functions in this element: 4\n",
+      "29-Sep 00:06:50:INFO:root:Initializing element\n",
+      "29-Sep 00:06:50:INFO:root:Element dimension: 2\n",
+      "29-Sep 00:06:50:INFO:root:Creating new scalar field JuliaFEM.Field{Array{Int64,1}}(0.0,1,[1,2,3,4])\n",
+      "29-Sep 00:06:50:INFO:root:Pushing field to element.\n",
+      "29-Sep 00:06:50:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
+      "29-Sep 00:06:51:INFO:root:Value: 2.5\n"
      ]
     },
     {
-     "ename": "LoadError",
-     "evalue": "LoadError: MethodError: `convert` has no method matching convert(::Type{Array{JuliaFEM.Field{T},1}}, ::JuliaFEM.Field{Int64})\nThis may have arisen from a call to the constructor Array{JuliaFEM.Field{T},1}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n  call{T}(::Type{T}, ::Any)\n  convert{T}(::Type{Array{T,1}}, !Matched::Range{T})\n  convert{T,S,N}(::Type{Array{T,N}}, !Matched::SubArray{S,N,P<:AbstractArray{T,N},I<:Tuple{Vararg{Union{AbstractArray{T,1},Colon,Int64}}},LD})\n  ...\nwhile loading In[6], in expression starting on line 2",
-     "output_type": "error",
-     "traceback": [
-      "LoadError: MethodError: `convert` has no method matching convert(::Type{Array{JuliaFEM.Field{T},1}}, ::JuliaFEM.Field{Int64})\nThis may have arisen from a call to the constructor Array{JuliaFEM.Field{T},1}(...),\nsince type constructors fall back to convert methods.\nClosest candidates are:\n  call{T}(::Type{T}, ::Any)\n  convert{T}(::Type{Array{T,1}}, !Matched::Range{T})\n  convert{T,S,N}(::Type{Array{T,N}}, !Matched::SubArray{S,N,P<:AbstractArray{T,N},I<:Tuple{Vararg{Union{AbstractArray{T,1},Colon,Int64}}},LD})\n  ...\nwhile loading In[6], in expression starting on line 2",
-      ""
-     ]
+     "data": {
+      "text/plain": [
+       "PipeEndpoint(open, 0 bytes waiting)"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -195,31 +223,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
-   "metadata": {
-    "collapsed": false
-   },
-   "outputs": [
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "24-Sep 21:04:04:INFO:root:Value: [2.5]\n",
-      "24-Sep 21:04:04:INFO:root:Element MyQuad4 passed tests.\n"
-     ]
-    }
-   ],
-   "source": [
-    "using JuliaFEM: set_field, interpolate\n",
-    "el1 = MyQuad4([1, 2, 3, 4])\n",
-    "set_field(el1, :temperature, [1, 2, 3, 4])\n",
-    "set_field(el1, :Geometry, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]);\n",
-    "set_field(el1, :\"heat coefficient\", 1);"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 7,
    "metadata": {
     "collapsed": false
    },
@@ -227,22 +231,59 @@
     {
      "data": {
       "text/plain": [
-       "2.5"
+       "2-element Array{JuliaFEM.Field{T},1}:\n",
+       " JuliaFEM.Field{Int64}(0.0,1,2)\n",
+       " JuliaFEM.Field{Int64}(1.0,1,3)"
       ]
      },
-     "execution_count": 9,
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "29-Sep 00:06:51:INFO:root:Element MyQuad4 passed tests.\n"
+     ]
+    }
+   ],
+   "source": [
+    "using JuliaFEM: new_field!, push_field!, interpolate, dinterpolate\n",
+    "el1 = MyQuad4([1, 2, 3, 4])\n",
+    "new_field!(el1, :temperature, Field(0.0, [0.0, 0.0, 0.0, 0.0]))\n",
+    "push_field!(el1, :temperature, Field(1.0, [1.0, 2.0, 3.0, 4.0]))\n",
+    "new_field!(el1, :Geometry, Field(0.0, Vector[[0.0,0.0,0.0], [10.0,0.0,0.0], [10.0,1.0,0.0], [0.0,1.0,0.0]]))\n",
+    "new_field!(el1, \"heat coefficient\", Field(0.0, 2))\n",
+    "push_field!(el1, \"heat coefficient\", Field(1.0, 3))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "1.25"
+      ]
+     },
+     "execution_count": 8,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
-    "# temperature at the middle poinf of the element, 1/4*(1+2+3+4)\n",
-    "interpolate(el1, :temperature, [0.0, 0.0])"
+    "# temperature at the middle poinf of the element, 1/4*(1+2+3+4) at t=0.5\n",
+    "interpolate(el1, :temperature, [0.0, 0.0], 0.5)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 9,
    "metadata": {
     "collapsed": false
    },
@@ -256,14 +297,39 @@
        " 0.0"
       ]
      },
-     "execution_count": 10,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "# geometry midpoint of element\n",
-    "interpolate(el1, :Geometry, [0.0, 0.0])"
+    "interpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "3x2 Array{Float64,2}:\n",
+       " 5.0  0.0\n",
+       " 0.0  0.5\n",
+       " 0.0  0.0"
+      ]
+     },
+     "execution_count": 10,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "dinterpolate(el1, :Geometry, [0.0, 0.0], -Inf)"
    ]
   },
   {
@@ -276,7 +342,7 @@
     {
      "data": {
       "text/plain": [
-       "1"
+       "([0.625,0.625,0.625,0.625],[0.75,0.75,0.75,0.75])"
       ]
      },
      "execution_count": 11,
@@ -286,7 +352,7 @@
    ],
    "source": [
     "# interpolating scalar -> scalar.\n",
-    "interpolate(el1, :\"heat coefficient\", [0.0, 0.0])"
+    "interpolate(el1, \"heat coefficient\", [0.0, 0.0], 0.5), interpolate(el1, \"heat coefficient\", [0.0, 0.0], Inf)"
    ]
   },
   {
@@ -304,10 +370,9 @@
    "source": [
     "## Developing own formulation\n",
     "\n",
-    "Let's consider Poisson equation\n",
+    "Let's consider a Laplace 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",
@@ -399,7 +464,7 @@
     "        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),\n",
     "        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]\n",
-    "    set_field(el, :temperature, zeros(2, 4)) # assign new field \"temperature\" to element\n",
+    "    new_field!(el, :temperature) # assign new field \"temperature\" to element\n",
     "    DC2D4(el, integration_points, [])\n",
     "end"
    ]
@@ -435,10 +500,11 @@
     "\"\"\"\n",
     "Left hand side defined in integration point\n",
     "\"\"\"\n",
-    "function JuliaFEM.get_lhs(eq::DC2D4, ip)\n",
+    "function JuliaFEM.get_lhs(eq::DC2D4, ip, t)\n",
     "    el = get_element(eq)\n",
-    "    dNdX = get_dbasisdX(el, ip.xi)\n",
-    "    hc = interpolate(el, :\"temperature thermal conductivity\", ip.xi)\n",
+    "    dNdX = get_dbasisdX(el, ip.xi, t)\n",
+    "    fld = el[\"temperature thermal conductivity\"](t)\n",
+    "    hc = sum(el(ip.xi) * fld)\n",
     "    return dNdX*hc*dNdX'\n",
     "end\n",
     "JuliaFEM.has_lhs(eq::DC2D4) = true"
@@ -476,17 +542,10 @@
    "source": [
     "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
     "el = Quad4([1, 2, 3, 4])\n",
-    "set_field(el, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
-    "set_field(el, :\"temperature thermal conductivity\", 6)\n",
+    "new_field!(el, :Geometry, Field(0.0, Vector[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]))\n",
+    "new_field!(el, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
     "eq = DC2D4(el)\n",
-    "integrate_lhs(eq)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "If rhs or lhs is not defined, integration returns nothing."
+    "integrate_lhs(eq, 1.0)"
    ]
   },
   {
@@ -508,14 +567,14 @@
     }
    ],
    "source": [
-    "isa(integrate_rhs(eq), Void)"
+    "JuliaFEM.has_lhs(eq)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Next heat flux on boundary:"
+    "If rhs or lhs is not defined, integration returns nothing."
    ]
   },
   {
@@ -528,7 +587,7 @@
     {
      "data": {
       "text/plain": [
-       "has_rhs (generic function with 2 methods)"
+       "true"
       ]
      },
      "execution_count": 18,
@@ -536,6 +595,35 @@
      "output_type": "execute_result"
     }
    ],
+   "source": [
+    "isa(integrate_rhs(eq, 1.0), Void)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next heat flux on boundary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "has_rhs (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 19,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
    "source": [
     "using JuliaFEM: get_basis, Seg2\n",
     "\n",
@@ -551,25 +639,25 @@
     "function DC2D2(el::Seg2)\n",
     "    integration_points = [\n",
     "        IntegrationPoint([0], 2.0)]\n",
-    "    set_field(el, :temperature, zeros(2, 1))\n",
+    "    new_field!(el, :temperature)\n",
     "    DC2D2(el, integration_points, [])\n",
     "end\n",
     "\n",
     "\"\"\"\n",
     "Right hand side defined in integration point\n",
     "\"\"\"\n",
-    "function JuliaFEM.get_rhs(eq::DC2D2, ip)\n",
+    "function JuliaFEM.get_rhs(eq::DC2D2, ip, t)\n",
     "    el = get_element(eq)\n",
-    "    N = get_basis(el, ip.xi)\n",
-    "    f = interpolate(el, :\"temperature flux\", ip.xi)\n",
-    "    return f*N\n",
+    "    ϕ = get_basis(el)\n",
+    "    f = el[\"temperature flux\"]\n",
+    "    return ϕ(ip.xi)*f(t)\n",
     "end\n",
     "JuliaFEM.has_rhs(eq::DC2D2) = true"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 20,
    "metadata": {
     "collapsed": false
    },
@@ -582,17 +670,17 @@
        " 50.0"
       ]
      },
-     "execution_count": 19,
+     "execution_count": 20,
      "metadata": {},
      "output_type": "execute_result"
     }
    ],
    "source": [
     "el = Seg2([1, 2])\n",
-    "set_field(el, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
-    "set_field(el, :\"temperature flux\", 100.0)\n",
+    "new_field!(el, :Geometry, Field(0.0, Vector[[0.0,0.0], [0.0,1.0]]))\n",
+    "new_field!(el, \"temperature flux\", Field(0.0, 100.0))\n",
     "eq = DC2D2(el)\n",
-    "integrate_rhs(eq)"
+    "integrate_rhs(eq, 1.0)"
    ]
   },
   {
@@ -607,7 +695,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 21,
    "metadata": {
     "collapsed": false
    },
@@ -618,7 +706,7 @@
        "PlaneHeatProblem"
       ]
      },
-     "execution_count": 20,
+     "execution_count": 21,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -627,14 +715,14 @@
     "using JuliaFEM: Problem, get_equation, get_dimension\n",
     "\n",
     "type PlaneHeatProblem <: Problem\n",
-    "    equations :: Array{Any, 1}\n",
+    "    equations :: Array{Equation, 1}\n",
     "end\n",
     "PlaneHeatProblem() = PlaneHeatProblem([])"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 22,
    "metadata": {
     "collapsed": false
    },
@@ -645,7 +733,7 @@
        "get_equation (generic function with 3 methods)"
       ]
      },
-     "execution_count": 21,
+     "execution_count": 22,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -674,7 +762,7 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "24-Sep 21:04:17:DEBUG:root:Problem (matrix) dimension: 4\n"
+      "29-Sep 00:06:53:DEBUG:root:Problem (matrix) dimension: 4\n"
      ]
     },
     {
@@ -687,7 +775,7 @@
        " -2.0  -1.0   4.0  -1.0\n",
        " -1.0  -2.0  -1.0   4.0,\n",
        "\n",
-       "[50.0,50.0,0.0,0.0])"
+       "[300.0,300.0,0.0,0.0])"
       ]
      },
      "execution_count": 23,
@@ -700,11 +788,11 @@
     "\n",
     "# create elements and add necessary properties like connectivity and geometry\n",
     "el1 = Quad4([1, 2, 3, 4])\n",
-    "set_field(el1, :Geometry, Vector[[0,0], [1,0], [1,1], [0,1]])\n",
-    "set_field(el1, :\"temperature thermal conductivity\", 6)\n",
+    "new_field!(el1, :Geometry, Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))\n",
+    "new_field!(el1, \"temperature thermal conductivity\", Field(0.0, 6.0))\n",
     "el2 = Seg2([1, 2])\n",
-    "set_field(el2, :Geometry, Vector[[0.0,0.0], [0.0,1.0]])\n",
-    "set_field(el2, :\"temperature flux\", 100.0)\n",
+    "new_field!(el2, :Geometry, Field(0.0, Vector[[0.0, 0.0], [0.0, 1.0]]))\n",
+    "new_field!(el2, \"temperature flux\", Field(1.0, 600.0))\n",
     "\n",
     "problem = PlaneHeatProblem()\n",
     "add_element!(problem, el1)\n",
@@ -716,20 +804,108 @@
     "n = get_matrix_dimension(problem)\n",
     "\n",
     "# integrate and assembly\n",
+    "t = 1.0\n",
     "A = zeros(n, n)\n",
     "b = zeros(n)\n",
     "for eq in get_equations(problem)\n",
     "    dofs = get_global_dofs(eq)\n",
     "    if has_lhs(eq)\n",
-    "        A[dofs, dofs] += integrate_lhs(eq)\n",
+    "        A[dofs, dofs] += integrate_lhs(eq, t)\n",
     "    end\n",
     "    if has_rhs(eq)\n",
-    "        b[dofs] += integrate_rhs(eq)\n",
+    "        b[dofs] += integrate_rhs(eq, t)\n",
     "    end\n",
     "end\n",
     "A, b"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2-element Array{Float64,1}:\n",
+       " 100.0\n",
+       " 100.0"
+      ]
+     },
+     "execution_count": 24,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fdofs = [1, 2]\n",
+    "A[fdofs, fdofs] \\ b[fdofs]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We still need to consider Dirichlet boundary conditions:\n",
+    "\\begin{align}\n",
+    "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n",
+    "\\end{align}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2x4 Array{Float64,2}:\n",
+       " 0.0  0.0  0.333333  0.166667\n",
+       " 0.0  0.0  0.166667  0.333333"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "M = 1/6*[0 0 2 1; 0 0 1 2]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "6-element Array{Float64,1}:\n",
+       " 100.0        \n",
+       " 100.0        \n",
+       "   1.36187e-14\n",
+       "   1.77636e-15\n",
+       " 600.0        \n",
+       " 600.0        "
+      ]
+     },
+     "execution_count": 26,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "[A M'; M zeros(2, 2)] \\ [b; 0; 0]"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -742,13 +918,15 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 0.5.0-dev",
+   "display_name": "Julia 0.4.0-rc2",
    "language": "julia",
-   "name": "julia-0.5"
+   "name": "julia-0.4"
   },
   "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
    "name": "julia",
-   "version": "0.5.0"
+   "version": "0.4.0"
   }
  },
  "nbformat": 4,
diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index c1e2db4..6f57e0e 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -8,7 +8,13 @@ using Logging
 @Logging.configure(level=DEBUG)
 
 include("types.jl") # type definitions
-include("elements.jl") # elements
+
+### ELEMENTS ###
+include("elements.jl")
+include("lagrange.jl") # Lagrange elements
+#include("hierarchical.jl") # P-elements
+
+
 include("equations.jl") # formulations
 include("problems.jl") # problems
 
diff --git a/src/elements.jl b/src/elements.jl
index 23ed479..cae93f4 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -90,16 +90,13 @@ End of example.
 
 # These must be implemented for your own element
 get_number_of_basis_functions(el::Type{Element}) = nothing
-get_number_of_basis_functions(el::Element) = nothing
-get_element_dimension(el::Element) = nothing
-get_dbasisdxi(el::Element, xi) = nothing
-get_connectivity(el::Element) = el.connectivity
+get_element_dimension(el::Type{Element}) = nothing
 
 ### LAGRANGE ELEMENTS ###
-include("lagrange.jl")
+#include("lagrange.jl")
 
 ### HIERARCHICAL P-ELEMENTS ###
-include("hierarchical.jl")
+#include("hierarchical.jl")
 
 ### COMMON ELEMENT ROUTINES ###
 
@@ -144,14 +141,15 @@ function test_element(eltype)
 
     # try to interpolate some scalar field
     fld = Field(0.0, collect(1:n))
-    Logging.info("Pushing scalar field $fld to element.")
+    Logging.info("Creating new scalar field $fld")
+    Logging.info("Pushing field to element.")
     new_field!(el, :field1)
     push_field!(el, :field1, fld)
     @fact el[:field1][1] --> fld
 
     mid = zeros(dim)
     try
-        f = get_basis(el)(mid)
+        get_basis(el)(mid)
     catch
         Logging.error("""
         Unable to evaluate basis, define function 'get_basis' for
@@ -167,11 +165,12 @@ function test_element(eltype)
 
     Logging.info("Interpolating scalar field at $mid")
     f(field, xi, t) = el(xi)*el[field](t)
-    i = f(:field, mid, 0.0)
+    i = f(:field1, mid, 0.0)
     Logging.info("Value: $i")
     Logging.info("Element $eltype passed tests.")
 end
 
+get_connectivity(el::Element) = el.connectivity
 
 """
 Get basis functions of element.
@@ -180,6 +179,32 @@ get_basis(el::Element) = el.basis
 get_basis(el::Element, xi::Vector) = el.basis(xi)
 Base.call(el::Element, xi::Vector) = el.basis(xi)
 
+"""
+Get partial derivatives of basis functions of element.
+"""
+get_dbasisdxi(el::Element) = el.basis.dbasisdxi
+get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi)
+
+"""
+Interpolate field on element.
+"""
+function interpolate(el::Element, field::Symbol, xi::Vector, t::Number)
+    get_basis(el, xi)*el[field](t)
+end
+function interpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
+    interpolate(el, Symbol(field), xi, t)
+end
+
+"""
+Interpolate derivative of field on element.
+"""
+function dinterpolate(el::Element, field::Symbol, xi::Vector, t::Number)
+    get_dbasisdxi(el, xi)*el[field](t)
+end
+function dinterpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
+    dinterpolate(el, Symbol(field), xi, t)
+end
+
 """
 Get jacobian of element evaluated at point ξ on element in reference configuration.
 
@@ -188,6 +213,7 @@ Parameters
 el::Element
 xi::Vector
 geometry_field::Any, optional
+time::Number, optional, default=0.0
 
 Returns
 -------
@@ -198,8 +224,8 @@ Notes
 -----
 Big "J" comes from reference (undeformed) configuration.
 """
-function get_Jacobian(el::Element, xi, geometry_field=:Geometry)
-    dinterpolate(el, geometry_field, xi)
+function get_Jacobian(el::Element, xi, t, geometry_field=:Geometry)
+    dinterpolate(el, geometry_field, xi, t)
 end
 
 
@@ -210,11 +236,11 @@ Notes
 -----
 Small "j" comes from current (deformed) configuration.
 """
-function get_jacobian(el::Element, xi, geometry_field=:Geometry, displacement_field=:displacement)
+function get_jacobian(el::Element, xi, t, geometry_field=:Geometry, displacement_field=:displacement)
     dbasisdxi = get_dbasisdxi(el, xi)
-    X = get_field(el, geometry_field)
-    u = get_field(el, displacement_field)
-    j = (X+u)*dbasisdxi
+    X = get_field(el, geometry_field)(t)
+    u = get_field(el, displacement_field)(t)
+    j = dbasisdxi*(X+u)
     return j
 end
 
@@ -222,9 +248,9 @@ end
 """
 Evaluate partial derivatives of basis, dbasis/dX
 """
-function get_dbasisdX(el::Element, xi)
+function get_dbasisdX(el::Element, xi, t)
     dbasisdxi = get_dbasisdxi(el, xi)
-    J = get_Jacobian(el, xi)
+    J = get_Jacobian(el, xi, t)
     dbasisdxi*inv(J)
 end
 
@@ -232,32 +258,48 @@ end
 """
 Evaluate partial derivatives of basis, dbasis/dx
 """
-function get_dbasisdx(el::Element, xi)
+function get_dbasisdx(el::Element, xi, t)
     dbasisdxi = get_dbasisdxi(el, xi)
-    j = get_jacobian(el, xi)
+    j = get_jacobian(el, xi, t)
     dbasisdxi*inv(j)
 end
 
+
 """ Create new empty field of some type. """
-function new_field!(el::Element, field_name)
+function new_field!(el::Element, field_name::Symbol)
     el.fields[field_name] = Field[]
 end
+function new_field!(el::Element, field_name::Symbol, field::Field)
+    new_field!(el, field_name)
+    push_field!(el, field_name, field)
+end
+function new_field!(el::Element, field_name::ASCIIString, field::Field)
+    new_field!(el, Symbol(field_name), field)
+end
 
 
 """ Push to existing set field of fields. """
-function push_field!(el::Element, field_name, field::Field)
+function push_field!(el::Element, field_name::Symbol, field::Field)
     push!(el.fields[field_name], field)
 end
+function push_field!(el::Element, field_name::ASCIIString, field::Field)
+    push_field!(el, Symbol(field_name), field)
+end
 
 
 """ Get field variable. """
-function get_field(el::Element, field_name)
+function get_field(el::Element, field_name::Symbol)
     el.fields[field_name]
 end
-function Base.getindex(el::Element, field_name)
-    el.fields[field_name]
+function get_field(el::Element, field_name::ASCIIString)
+    el.fields[Symbol(field_name)]
+end
+function Base.getindex(el::Element, field_name::Union{ASCIIString, Symbol})
+    get_field(el, field_name)
 end
 
+
+#=
 """
 Evaluate some field in point ξ on element using basis functions.
 
@@ -296,8 +338,6 @@ function interpolate(el::Element, field, xis::Array{Vector, 1})
     map(interpolate_, xis)
 end
 
-"""
-"""
 function dinterpolate(el::Element, field, xi::Number)
     dinterpolate(el, field, [xi])
 end
@@ -309,12 +349,13 @@ function dinterpolate(el::Element, field, xi::Vector)
     end
     return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
 end
+=#
 
 """
 calculate "local" normals in elements, in a way that
 n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1]
 """
-function calculate_normals!(el::Element, field_name=:Normals)
+function calculate_normals!(el::Element, t, field_name=:Normals)
     new_field!(el, field_name, Vector)
     for xi in Vector[[-1.0], [1.0]]
         t = dinterpolate(el, :Geometry, xi)
diff --git a/src/equations.jl b/src/equations.jl
index 1676955..89a40b0 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -32,23 +32,23 @@ get_integration_points(eq::Equation) = eq.integration_points
 get_connectivity(eq::Equation) = get_connectivity(get_element(eq))
 get_basis(eq::Equation, ip::IntegrationPoint) = get_basis(get_element(eq), ip.xi)
 get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq), ip.xi)
-interpolate(eq::Equation, field::Union(ASCIIString, Symbol), ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
-integrate_lhs(eq::Equation) = has_lhs(eq) ? integrate(eq, get_lhs) : nothing
-integrate_rhs(eq::Equation) = has_rhs(eq) ? integrate(eq, get_rhs) : nothing
+interpolate(eq::Equation, field::Union{ASCIIString, Symbol}, ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
+integrate_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
+integrate_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
 
 
 """
 Return determinant of Jacobian for numerical integration.
 """
-function get_detJ(eq::Equation, ip::IntegrationPoint)
+function get_detJ(eq::Equation, ip::IntegrationPoint, t::Float64)
     el = get_element(eq)
-    get_detJ(el, ip)
+    get_detJ(el, ip, t)
 end
-function get_detJ(el::Element, ip::IntegrationPoint)
-    J = get_detJ(el, ip.xi)
+function get_detJ(el::Element, ip::IntegrationPoint, t::Float64)
+    get_detJ(el, ip.xi, t)
 end
-function get_detJ(el::Element, xi::Vector)
-    J = get_Jacobian(el, xi)
+function get_detJ(el::Element, xi::Vector, t::Float64)
+    J = get_Jacobian(el, xi, t)
     s = size(J)
     return s[1] == s[2] ? det(J) : norm(J)
 end
@@ -63,10 +63,10 @@ eq::Equation
 f::Function
     Function to integrate
 """
-function integrate(eq::Equation, f::Function)
+function integrate(eq::Equation, f::Function, t::Float64)
     target = []
     for ip in get_integration_points(eq)
-        push!(target, ip.weight*f(eq, ip)*get_detJ(eq, ip))
+        push!(target, ip.weight*f(eq, ip, t)*get_detJ(eq, ip, t))
     end
     return sum(target)
 end
@@ -80,4 +80,4 @@ function set_global_dofs!(eq::Equation, dofs)
 end
 
 # Equations for heat problems
-include("heat_equations.jl")
+#include("heat_equations.jl")
diff --git a/src/lagrange.jl b/src/lagrange.jl
index 5ee8115..92901c4 100644
--- a/src/lagrange.jl
+++ b/src/lagrange.jl
@@ -5,31 +5,9 @@
 
 abstract CG <: Element
 
-"""
-Create new element with element_name to family element_family
-
-Examples
---------
->>> @create_element(Seg2, CG, "2 node linear segment")
-"""
-macro create_element(element_name, element_family, element_description)
-#   Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
-    eltype = esc(element_name)
-    elfam = esc(element_family)
-    quote
-        global get_element_description
-        type $eltype <: $elfam
-            connectivity :: Array{Int, 1}
-            fields :: Dict{Any, Any}
-        end
-        $eltype(connectivity) = $eltype(connectivity, Dict{Any, Any}())
-        get_element_description(el::Type{$eltype}) = $element_description
-    end
-end
-
 """
 Given polynomial P and coordinates of reference element, calculate
-Lagrange basis function and partial derivatives.
+Lagrange basis functions
 """
 function calculate_lagrange_basis(P, X)
     dim, nbasis = size(X)
@@ -40,70 +18,66 @@ function calculate_lagrange_basis(P, X)
 #   Logging.debug("Calculating inverse of A")
     invA = inv(A)'
     basis(xi) = invA*P(xi)
-    dbasisdxi = ForwardDiff.jacobian(basis)
-    basis, dbasisdxi
+    basis
 end
 
 """
-Assign Lagrange basis for element.
+Create new Lagrange element
+
+Examples
+--------
+>>> @create_lagrange_element(Seg2, "2 node linear segment", X, P)
 """
-macro create_lagrange_basis(element_name, X, P)
-
-#   Logging.debug("Creating Lagrange basis for element ", element_name, ". ")
+macro create_lagrange_element(element_name, element_description, X, P)
+#   Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
     eltype = esc(element_name)
-
     quote
-
+        global get_element_description
         global get_number_of_basis_functions, get_element_dimension
-        global get_basis, get_dbasisdxi
-
         dim = size($X, 1)
         nbasis = size($X, 2)
-#       Logging.debug("Number of basis functions: ", nbasis, ". ")
-#       Logging.debug("Element dimension: ", dim)
-
-        get_number_of_basis_functions(el::Type{$(esc(element_name))}) = nbasis
-        get_number_of_basis_functions(el::$(esc(element_name))) = nbasis
-        get_element_dimension(el::$(esc(element_name))) = dim
-
-        basis, dbasisdxi = calculate_lagrange_basis($P, $X)
-        get_basis(el::$eltype, xi) = basis(xi)
-        get_dbasisdxi(el::$eltype, xi) = dbasisdxi(xi)
-#       Logging.debug("Element ", $element_name, " created.")
+        h = calculate_lagrange_basis($P, $X)
+        type $eltype <: CG
+            connectivity :: Array{Int, 1}
+            basis :: Basis
+            fields :: Dict{Symbol, Array{Field, 1}}
+        end
+        function $eltype(connectivity, args...)
+            $eltype(connectivity, Basis(h), Dict())
+        end
+        get_element_description(el::Type{$eltype}) = $element_description
+        get_number_of_basis_functions(el::Type{$eltype}) = nbasis
+        get_element_dimension(el::Type{$eltype}) = dim
     end
-
 end
 
 # 0d Lagrange element
 
-@create_element(Point1, CG, "1 node point element")
+#@create_element(Point1, CG, "1 node point element")
 
 # 1d Lagrange elements
 
-@create_element(Seg2, CG, "2 node linear line element")
-@create_lagrange_basis(Seg2, [-1.0 1.0], (xi) -> [1.0, xi[1]])
+@create_lagrange_element(Seg2, "2 node linear line element",
+    [-1.0 1.0], (xi) -> [1.0, xi[1]])
 
-@create_element(Seg3, CG, "3 node quadratic line element")
-@create_lagrange_basis(Seg3, [-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
+@create_lagrange_element(Seg3, "3 node quadratic line element",
+    [-1.0 1.0 0.0], (xi) -> [1.0, xi[1], xi[1]^2])
 
 # 2d Lagrange elements
 
-@create_element(Tri3, CG, "3 node bilinear triangle element")
-@create_lagrange_basis(Tri3,
+@create_lagrange_element(Tri3, "3 node bilinear triangle element",
     [0.0 1.0 0.0
      0.0 0.0 1.0],
     (xi) -> [1.0, xi[1], xi[2]])
 
-@create_element(Quad4, CG, "4 node bilinear quadrangle element")
-@create_lagrange_basis(Quad4,
+@create_lagrange_element(Quad4, "4 node bilinear quadrangle element",
     [-1.0  1.0 1.0 -1.0
      -1.0 -1.0 1.0  1.0],
     (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2]])
 
 # 3d Lagrange elements
 
-@create_element(Tet10, CG, "10 node quadratic tetrahedron")
-@create_lagrange_basis(Tet10,
+@create_lagrange_element(Tet10, "10 node quadratic tetrahedron",
     [0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
      0.0 0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5
      0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.5 0.5 0.5],
diff --git a/src/math.jl b/src/math.jl
index 8326d78..c11c414 100644
--- a/src/math.jl
+++ b/src/math.jl
@@ -1,69 +1,8 @@
 # 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, ...
-
 using ForwardDiff
 
-#export interpolate, integrate, linearize
-
-"""
-Interpolate field variable using basis functions f for point ip.
-This function tries to be as general as possible and allows interpolating
-lot of different fields.
-
-Parameters
-----------
-field :: Array{Number, dim}
-  Field variable
-basis :: Function
-  Basis functions
-ip :: Array{Number, 1}
-  Point to interpolate
-"""
-function interpolate(field::Float64, basis::Function, ip::Array{Float64,1})
-    # dummy function, unable to interpolate scalar value!
-    return field
-end
-function interpolate{T<:Real}(field::Array{T,1}, basis::Function, ip)
-    result = dot(field, basis(ip))
-    return result
-end
-function interpolate{T<:Real}(field::Array{T,2}, basis::Function, ip)
-    m, n = size(field)
-    bip = basis(ip)
-    tmp = size(bip)
-    if length(tmp) == 1
-        ndim = 1
-        nnodes = tmp[1]
-    else
-        ndim, nnodes = size(bip)
-    end
-    if ndim == 1
-        if n == nnodes
-            result = field * bip
-        elseif m == nnodes
-            result = field' * bip
-        end
-    else
-      if n == nnodes
-        result = bip' * field
-    elseif m == nnodes
-      result = bip' * field'
-    end
-    end
-    if length(result) == 1
-        result = result[1]
-    end
-    return result
-end
-#function interpolate(e::Element, field::ASCIIString, x::Array{Float64,1}; derivative=false)
-#    basis = derivative ? get_dbasisdxi(e) : get_basis(e)
-#    return interpolate(e.attributes[field], basis, x)
-#end
-
-
-
 """
 Linearize function f w.r.t some given field, i.e. calculate dR/du
 
@@ -78,7 +17,6 @@ Returns
 -------
 Array{Float64, 2}
     jacobian / "tangent stiffness matrix"
-
 """
 function linearize(f::Function, el::Element, field::ASCIIString)
     dim, nnodes = size(el.attributes[field])
@@ -92,6 +30,7 @@ function linearize(f::Function, el::Element, field::ASCIIString)
     return jac(el.attributes[field][:])
 end
 
+
 """
 This version returns another function which can be then evaluated against field
 """
@@ -111,6 +50,7 @@ function linearize(f::Function, field::ASCIIString)
     return jacobian
 end
 
+
 """
 In-place version, no additional garbage collection.
 """
@@ -128,8 +68,6 @@ function linearize!(f::Function, el::Element, field::ASCIIString, target::ASCIIS
 end
 
 
-
-
 """
 This version returns a function which must be operated with element e
 """
@@ -157,25 +95,6 @@ function integrate!(f::Function, el::Element, target)
     end
 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)
@@ -194,4 +113,3 @@ function linearize(eq::Equation, f::Function, field::ASCIIString)
     return jacobian
 end
 
-
diff --git a/src/types.jl b/src/types.jl
index 455baa3..6d03c90 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -10,7 +10,7 @@ using ForwardDiff
 type Field{T}
     time :: Float64
     increment :: Int64
-    values :: Array{T, 1}
+    values :: T
 end
 
 """ Initialize field. """
@@ -25,8 +25,16 @@ Base.length(f::Field) = length(f.values)
 Base.getindex(f::Field, i::Int64) = f.values[i]
 
 """ Interpolate field h(ξ)*f = x*f """
-Base.(:*)(x::Array{Float64, 1}, f::Field) = sum(x .* f.values)
-Base.(:*)(x::Array{Float64, 2}, f::Field) = sum([f[i]*x[i,:] for i in 1:length(f)])
+function interpolate{T}(x::Vector, f::Field{Vector{T}})
+    sum([f[i]*x[i] for i in 1:length(f)])
+end
+function interpolate{T}(x::Matrix, f::Field{Vector{T}})
+    sum([f[i]*x[i,:] for i in 1:length(f)])
+end
+function interpolate(x::Vector, f::Field)
+    f.values*x
+end
+Base.(:*)(x::Union{Vector, Matrix}, f::Field) = interpolate(x, f)
 
 """ Interpolate field (h*f)(ξ) """
 Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
diff --git a/test/test_elements.jl b/test/test_elements.jl
index ea962aa..7353438 100644
--- a/test/test_elements.jl
+++ b/test/test_elements.jl
@@ -2,13 +2,7 @@
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
 using FactCheck
+using JuliaFEM: test_element
 
-facts("test set and interpolate field variable") do
-    el = Quad4(1, [1, 2, 3, 4])
-    set_coordinates(el, [0.0 0.0; 10.0 0.0; 10.0 1.0; 0.0 1.0]')
-    set_field(el, "displacement", [0.0 0.0; 0.0 0.0; 0.5 0.0; 0.0 0.0]'')
-    fval = interpolate(el, "displacement", [0.0, 1.0])
-    Logging.debug(fval)
-    @fact fval --> roughly([0.25 0.0]')
-end
-
+using JuliaFEM: Quad4
+test_element(Quad4)
diff --git a/test/test_types.jl b/test/test_types.jl
index 86c3d2e..a94b88b 100644
--- a/test/test_types.jl
+++ b/test/test_types.jl
@@ -57,4 +57,8 @@ facts("test fields and interpolation") do
    @fact diff(h)([0.0, 0.0])*X --> [0.5 0.0; 0.0 0.5]
    @fact (diff(h)*X)([0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
 
+   b = Basis((xi) -> [1/2*(1-xi[1]), 1/2*(1+xi[1])])
+   f = Field(0.0, 100.0)
+   @fact b(0.0) * f --> [50.0, 50.0]
+
 end