From b90ecf1772202d63b1ca4c6193ddcd7b8a68385e Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Sun, 30 Aug 2015 23:20:14 +0300 Subject: [PATCH] some developer guidelines --- .../2015-08-29-developing-juliafem.ipynb | 521 ++++++++++++++++++ 1 file changed, 521 insertions(+) create mode 100644 notebooks/2015-08-29-developing-juliafem.ipynb diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb new file mode 100644 index 0000000..b876d35 --- /dev/null +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -0,0 +1,521 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Developing JuliaFEM\n", + "\n", + "Author(s): Jukka Aho\n", + "\n", + "**Abstract**: General developer notes." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Logger(root,DEBUG,Pipe(open, 0 bytes waiting),root)" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using Logging\n", + "using FactCheck\n", + "Logging.configure(level=DEBUG)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Developing own element\n", + "\n", + "Finite element definition, from [FEniCS-book](https://bitbucket.org/fenics-project/fenics-book/src/7d3a80e7dda0fc279c7964dc6000d57942f11eb3/fenicsbook.cls?at=master) [Ciarlet, 2002]:\n", + "\n", + "- the domain $T$ is a bounded, closed subset of $\\mathbb{R}^d$ (for $d = 1, 2, 3, \\dots$) with nonempty interior and piecewise smooth boundary;\n", + "- the space $\\mathcal{V} = \\mathcal{V}(T)$ is a finite dimensional function space on $T$ of dimension $n$;\n", + "- the set of degrees of freedom (nodes) $\\mathcal{L} = \\{\\ell_1, \\ell_2,\\ldots, \\ell_{n}\\}$ is a basis for the dual space $\\mathcal{V}'$; that is, the space of bounded linear functionals on $\\mathcal{V}$.\n", + "\n", + "We extend this definition so that domain $T$ can also be empty.\n", + "\n", + "Minimum requirements for element:\n", + "- subclass from Element, if not wanting to implement everything by youself\n", + "- define basis and partial derivatives of it, because we need to interpolate over it\n", + "- give connectivity information, how this element is connected to other elements\n", + "- create proper constructor (see example).\n", + "\n", + "Test the element using ``test_element`` function. It it passes, then element implementation should be fine. As an example, we define 4 node quadrilateral element using linear Lagrange basis. We really don't care much how element is implemented as long it's interface is constructed with some rules. The interface is tested using `test_element` and it also gives information how to fix element if something is missing." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM: Element, Equation, set_field, get_field, interpolate, integrate_lhs, integrate_rhs" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here's one basic implementation. The actual element:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "type Quad4 <: Element\n", + " connectivity :: Array{Int, 1}\n", + " fields :: Dict{Any, Any}\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Default constructor, providing connectivity data needed in assembly" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "Quad4" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "Quad4(connectivity) = Quad4(connectivity, Dict{ASCIIString, Any}())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Some basic charasteristics like number of nodes / connectivity points and dimension:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "get_element_dimension (generic function with 2 methods)" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "JuliaFEM.get_number_of_nodes(el::Type{Quad4}) = 4\n", + "JuliaFEM.get_element_dimension(el::Quad4) = 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The most important, it's basis (we probably want to interpolate something with this element):" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "get_dbasisdxi (generic function with 2 methods)" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function JuliaFEM.get_basis(el::Quad4, xi::Array{Float64,1})\n", + " [(1-xi[1])*(1-xi[2])/4\n", + " (1+xi[1])*(1-xi[2])/4\n", + " (1+xi[1])*(1+xi[2])/4\n", + " (1-xi[1])*(1+xi[2])/4]\n", + "end\n", + "\n", + "function JuliaFEM.get_dbasisdxi(el::Quad4, xi::Array{Float64,1})\n", + " [-(1-xi[2])/4.0 -(1-xi[1])/4.0\n", + " (1-xi[2])/4.0 -(1+xi[1])/4.0\n", + " (1+xi[2])/4.0 (1+xi[1])/4.0\n", + " -(1+xi[2])/4.0 (1-xi[1])/4.0]\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we check that everything is well defined:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "30-Aug 23:18:00:INFO:root:number of connectivity points (nodes) in this element: 4\n", + "30-Aug 23:18:00:INFO:root:Constructing element..\n", + "30-Aug 23:18:00:INFO:root:Element dimension: 2\n", + "30-Aug 23:18:01:INFO:root:Setting scalar field [1 2 3 4] to element.\n", + "30-Aug 23:18:01:INFO:root:Interpolating scalar field\n", + "30-Aug 23:18:01:INFO:root:Element Quad4 passed tests.\n" + ] + } + ], + "source": [ + "using JuliaFEM: test_element\n", + "test_element(Quad4)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "If test_element passes, element should be well defined. At least in the sense that it has all necessary things defined ready to be used in JuliaFEM. After building element, one can interpolate things in it. Couple examples:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "el1 = Quad4([1, 2, 3, 4])\n", + "set_field(el1, :temperature, [1 2 3 4])\n", + "set_field(el1, :coordinates, [0.0 0.0 0.0; 10.0 0.0 0.0; 10.0 1.0 0.0; 0.0 1.0 0.0]');\n", + "set_field(el1, :\"heat coefficient\", 1);" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1-element Array{Float64,1}:\n", + " 2.5" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interpolate(el1, :temperature, [0.0, 0.0]) # temperature at the middle poinf of the element, 1/4*(1+2+3+4)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Array{Float64,1}:\n", + " 5.0\n", + " 0.5\n", + " 0.0" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interpolate(el1, :coordinates, [0.0, 0.0]) # midpoint of element" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "1" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "interpolate(el1, :\"heat coefficient\", [0.0, 0.0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Developing own formulation\n", + "\n", + "Let's consider Poisson equation\n", + "\\begin{align}\n", + "\\Delta{u} &= 0 && \\text{on } \\Omega \\\\\n", + "u &= u_0 && \\text{on } \\Gamma_{\\mathrm{D}} \\\\\n", + "\\frac{\\partial u}{\\partial n} &= g && \\text{on } \\Gamma_{\\mathrm{N}}\n", + "\\end{align}\n", + "\n", + "Weak form is, find $u\\in\\mathcal{U}$ such that\n", + "\\begin{equation}\n", + " \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x = \\int_{\\Gamma_{\\mathrm{N}}}g v\\,\\mathrm{d}s \\quad \\forall v\\in\\mathcal{V}.\n", + "\\end{equation}\n", + "\n", + "Minimum requirements for equation: \n", + "- subclass from Equation, if not want to implement from scratch\n", + "- it needs to have lhs and rhs functions\n", + "- default constructor takes the element as input argument\n", + "\n", + "Now we have function `test_equation`, which we can use to test that everything is working as expected. \n", + "\n", + "Again thanks to multiple dispatch, you are free to code your weak form however you want as long as it returns lhs and rhs sides for element dofs. This kind of freedom gives good opportunities to wrap e.g. Fortran code from some other projects. And again we have some suggestions ad following these ideas you get a lot of stuff for free. First we look the left hand side of the equation, that is,\n", + "\\begin{equation}\n", + " \\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "using JuliaFEM: Equation, IntegrationPoint\n", + "\n", + "abstract Heat <: Equation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Our basic data type often looks something like this:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "\"\"\"\n", + "Diffusive heat transfer for 4-node bilinear element.\n", + "\"\"\"\n", + "type DC2D4 <: Heat\n", + " element :: Quad4\n", + " integration_points :: Array{IntegrationPoint, 1}\n", + " global_dofs :: Array{Int64, 1}\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We must provide default constructor which takes element as input argument:" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "DC2D4" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "function DC2D4(el::Quad4)\n", + " integration_points = [\n", + " IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[ 1, 1], 1.0),\n", + " IntegrationPoint(1.0/sqrt(3.0)*[-1, 1], 1.0)]\n", + " set_field(el, \"temperature\", zeros(2, 4))\n", + " DC2D4(el, integration_points, [])\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now the actual implementation for $\\int_{\\Omega}\\nabla u\\cdot\\nabla v\\,\\mathrm{d}x$:" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "get_lhs (generic function with 2 methods)" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "using JuliaFEM: get_element, get_dbasisdX\n", + "\n", + "\"\"\"\n", + "Left hand side defined in integration point\n", + "\"\"\"\n", + "function JuliaFEM.get_lhs(eq::DC2D4, ip)\n", + " el = get_element(eq)\n", + " dNdX = get_dbasisdX(el, ip.xi)\n", + " hc = interpolate(el, :\"temperature heat coefficient\", ip.xi)\n", + " return dNdX*hc*dNdX'\n", + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And that's it. If we want to play with this formulation, we must create element and assign this equation for it:" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "4x4 Array{Float64,2}:\n", + " 4.0 -1.0 -2.0 -1.0\n", + " -1.0 4.0 -1.0 -2.0\n", + " -2.0 -1.0 4.0 -1.0\n", + " -1.0 -2.0 -1.0 4.0" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "el = Quad4([1, 2, 3, 4])\n", + "set_field(el, :coordinates, [0 0; 1 0; 1 1; 0 1]')\n", + "set_field(el, :\"temperature heat coefficient\", 6)\n", + "eq = DC2D4(el)\n", + "integrate_lhs(eq)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Julia 0.4.0-dev", + "language": "julia", + "name": "julia-0.4" + }, + "language_info": { + "name": "julia", + "version": "0.4.0" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +}