Merge branch 'master' of https://github.com/JuliaFEM/JuliaFEM.jl

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>`_.

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
+}

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
 

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
+
+

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