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

notebooks/2015-08-29-developing-juliafem.ipynb

@@ -66,7 +66,7 @@
    },
    "outputs": [],
    "source": [
-    "using JuliaFEM: Element, Equation, set_field, get_field, interpolate, integrate_lhs, integrate_rhs"
+    "using JuliaFEM: Element"
    ]
   },
   {
@@ -208,12 +208,13 @@
      "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"
+      "31-Aug 20:05:30:INFO:root:number of connectivity points (nodes) in this element: 4\n",
+      "31-Aug 20:05:31:INFO:root:Constructing element..\n",
+      "31-Aug 20:05:31:INFO:root:Element dimension: 2\n",
+      "31-Aug 20:05:31:INFO:root:Setting scalar field [1 2 3 4] to element.\n",
+      "31-Aug 20:05:31:INFO:root:Interpolating scalar field at [0.0,0.0]\n",
+      "31-Aug 20:05:31:INFO:root:Value: [2.5]\n",
+      "31-Aug 20:05:31:INFO:root:Element Quad4 passed tests.\n"
      ]
     }
    ],
@@ -237,6 +238,7 @@
    },
    "outputs": [],
    "source": [
+    "using JuliaFEM: set_field, interpolate\n",
     "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",
@@ -253,8 +255,7 @@
     {
      "data": {
       "text/plain": [
-       "1-element Array{Float64,1}:\n",
-       " 2.5"
+       "2.5"
       ]
      },
      "execution_count": 9,
@@ -414,7 +415,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))\n",
+    "    set_field(el, :temperature, zeros(2, 4))\n",
     "    DC2D4(el, integration_points, [])\n",
     "end"
    ]
@@ -453,7 +454,7 @@
     "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",
+    "    hc = interpolate(el, :\"temperature thermal conductivity\", ip.xi)\n",
     "    return dNdX*hc*dNdX'\n",
     "end"
    ]
@@ -488,13 +489,131 @@
     }
    ],
    "source": [
+    "using JuliaFEM: integrate, integrate_lhs, integrate_rhs\n",
     "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",
+    "set_field(el, :\"temperature thermal conductivity\", 6)\n",
     "eq = DC2D4(el)\n",
     "integrate_lhs(eq)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next heat flux on boundary:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "31-Aug 20:05:35:INFO:root:number of connectivity points (nodes) in this element: 2\n",
+      "31-Aug 20:05:35:INFO:root:Constructing element..\n",
+      "31-Aug 20:05:35:INFO:root:Element dimension: 1\n",
+      "31-Aug 20:05:35:INFO:root:Setting scalar field [1 2] to element.\n",
+      "31-Aug 20:05:35:INFO:root:Interpolating scalar field at [0.0]\n",
+      "31-Aug 20:05:35:INFO:root:Value: [1.5]\n",
+      "31-Aug 20:05:35:INFO:root:Element Seg2 passed tests.\n"
+     ]
+    }
+   ],
+   "source": [
+    "type Seg2 <: Element\n",
+    "    connectivity :: Array{Int, 1}\n",
+    "    fields :: Dict{Any, Any}\n",
+    "end\n",
+    "Seg2(connectivity) = Seg2(connectivity, Dict{ASCIIString, Any}())\n",
+    "JuliaFEM.get_number_of_nodes(el::Type{Seg2}) = 2\n",
+    "JuliaFEM.get_element_dimension(el::Seg2) = 1\n",
+    "JuliaFEM.get_basis(el::Seg2, xi::Array{Float64,1}) = [0.5*(1-xi[1]), 0.5*(1+xi[1])]\n",
+    "JuliaFEM.get_dbasisdxi(el::Seg2, xi::Array{Float64,1}) = [-0.5 0.5]'\n",
+    "test_element(Seg2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "get_rhs (generic function with 2 methods)"
+      ]
+     },
+     "execution_count": 29,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "using JuliaFEM: get_basis\n",
+    "\n",
+    "\"\"\"\n",
+    "Diffusive heat transfer for 4-node bilinear element.\n",
+    "\"\"\"\n",
+    "type DC2D2 <: Heat\n",
+    "    element :: Seg2\n",
+    "    integration_points :: Array{IntegrationPoint, 1}\n",
+    "    global_dofs :: Array{Int64, 1}\n",
+    "end\n",
+    "\n",
+    "function DC2D2(el::Seg2)\n",
+    "    integration_points = [\n",
+    "        IntegrationPoint([0], 2.0)]\n",
+    "    set_field(el, :temperature, zeros(2, 1))\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",
+    "    el = get_element(eq)\n",
+    "    N = get_basis(el, ip.xi)\n",
+    "    f = interpolate(el, :\"temperature flux\", ip.xi)\n",
+    "    return f*N\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "2-element Array{Float64,1}:\n",
+       " 50.0\n",
+       " 50.0"
+      ]
+     },
+     "execution_count": 30,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "el = Seg2([1, 2])\n",
+    "set_field(el, :coordinates, [0.0 0.0; 0.0 1.0]')\n",
+    "set_field(el, :\"temperature flux\", 100.0)\n",
+    "eq = DC2D2(el)\n",
+    "integrate_rhs(eq)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,

src/JuliaFEM.jl

@@ -7,10 +7,9 @@ using Lexicon
 using Logging
 @Logging.configure(level=DEBUG)
 
-Logging.info("loading types")
 include("types.jl") # type definitions
-Logging.info("loading elements")
 include("elements.jl") # elements
+include("equations.jl") # formulations
 include("math.jl") # basic mathematical operations
 
 include("elasticity_solver.jl")

src/elements.jl

@@ -1,16 +1,70 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
+using FactCheck
+
 abstract Element
 
-get_element(eq::Equation) = eq.element
+export get_number_of_nodes
+
+get_number_of_nodes(el::Type{Element}) = -1
+get_element_dimension(el::Element) = -1
+get_dbasisdx(el::Element) = nothing
+get_basis(el::Element) = nothing
+
+function test_element(eltype)
+    local el
+    n = get_number_of_nodes(eltype)
+    Logging.info("number of connectivity points (nodes) in this element: $n")
+    @fact n --> not(-1) """Unable to determine number of nodes for $eltype
+    define a function 'get_number_of_nodes' which returns the number of nodes for this element."""
+
+    Logging.info("Constructing element..")
+    try
+        el = eltype(collect(1:n))
+    catch
+        Logging.error("""Unable to create element with default constructor
+        define function $eltype(connectivity) which initializes this element.
+        """)
+    end
+    dim = get_element_dimension(el)
+    Logging.info("Element dimension: $dim")
+    @fact dim --> not(-1) """Unable to get element dimension
+    define function 'get_element_dimension' which return the dimension of this element (1, 2, 3)"""
+
+    # try to interpolate some scalar field
+    fld = collect(1:n)'
+    Logging.info("Setting scalar field $fld to element.")
+    set_field(el, "field1", fld)
+    @fact get_field(el, "field1") --> fld
+
+    try
+        get_basis(el, zeros(dim))
+    catch
+        Logging.error("""Unable to evaluate basis, define function 'get_basis' for this element.
+        """)
+    end
+    try
+        get_dbasisdxi(el, zeros(dim))
+    catch
+        Logging.error("""Unable to evaluate partial derivatives of basis, define function 'get_dbasisdxi' for this element.
+        """)
+    end
+
+    xi = zeros(dim)
+    Logging.info("Interpolating scalar field at $xi")
+    i = interpolate(el, "field1", zeros(dim))
+    Logging.info("Value: $i")
+    Logging.info("Element $eltype passed tests.")
+end
+
 
 """
 Get jacobian of element evaluated at point xi
 """
 function get_jacobian(el::Element, xi)
     dbasisdxi = get_dbasisdxi(el, xi)
-    X = get_field(el, "coordinates")
+    X = get_field(el, :coordinates)
     #J = interpolate(X, dbasisdxi, xi)'
     J = X*dbasisdxi
     return J
@@ -40,11 +94,29 @@ function get_field(el::Element, field_name)
     el.fields[field_name]
 end
 
+#"""
+#Evaluate field in point xi using basis functions.
+#"""
+#function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})
+#    (get_basis(el, xi)'*get_field(el, field))'
+#end
+
 """
 Evaluate field in point xi using basis functions.
 """
-function interpolate(el::Element, field::ASCIIString, xi::Array{Float64,1})
-    (get_basis(el, xi)'*get_field(el, field))'
+function interpolate(el::Element, field::Union(ASCIIString, Symbol), xi::Array{Float64,1})
+    f = get_field(el, field)
+    basis = get_basis(el, xi)
+    dim, nnodes = size(f)
+    result = zeros(dim)
+    for i=1:nnodes
+        result += basis[i]*f[:,i]
+    end
+    if dim == 1
+        return result[1]
+    else
+        return result
+    end
 end
 
 ### Lagrange family ###

src/equations.jl

@@ -0,0 +1,69 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+abstract Equation
+
+"""
+Integration point
+
+xi :: Array{Float64, 1}
+    (dimensionless) coordinates of integration point
+weight :: Float64
+    Integration weight
+attributes :: Dict{Any, Any}
+    This is used to save internal variables of IP needed e.g. for incremental
+    material models.
+"""
+type IntegrationPoint
+    xi :: Array{Float64, 1}
+    weight :: Float64
+    attributes :: Dict{Any, Any}
+end
+IntegrationPoint(xi, weight) = IntegrationPoint(xi, weight, Dict{ASCIIString, Any}())
+
+get_lhs(eq::Equation, xi) = nothing
+get_rhs(eq::Equation, xi) = nothing
+get_element(eq::Equation) = eq.element
+get_integration_points(eq::Equation) = eq.integration_points
+# couple convenient functions -- could make weak form definition easier
+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) = integrate(eq, get_lhs)
+integrate_rhs(eq::Equation) = integrate(eq, get_rhs)
+
+"""
+Return determinant of Jacobian for numerical integration.
+"""
+function get_detJ(eq::Equation, ip::IntegrationPoint)
+    el = get_element(eq)
+    get_detJ(el, ip)
+end
+function get_detJ(el::Element, ip::IntegrationPoint)
+    J = get_jacobian(el, ip.xi)
+    n, m = size(J)
+    if n != m # for manifolds
+        return norm(J)
+    else
+        return det(J)
+    end
+end
+
+"""
+Integrate f over element
+
+Parameters
+----------
+eq::Equation
+    
+f::Function
+    Function to integrate
+"""
+function integrate(eq::Equation, f::Function)
+    target = []
+    for ip in get_integration_points(eq)
+        push!(target, ip.weight*f(eq, ip)*get_detJ(eq, ip))
+    end
+    return sum(target)
+end
+

src/math.jl

@@ -129,32 +129,6 @@ end
 
 
 
-"""
-Integrate f over element using Gaussian quadrature rules.
-
-Parameters
-----------
-el::Element
-    well defined element
-f::Function
-    Function to integrate
-"""
-function integrate(f::Function, el::Element)
-    target = []
-    for ip in el.integration_points
-        J = interpolate(el, "coordinates", ip.xi; derivative=true)
-        push!(target, ip.weight*f(el, ip)*det(J))
-    end
-    return sum(target)
-end
-#function integrate(f::Function, integration_points::Array{IntegrationPoint, 1}, Xargs...)
-#    target = []
-#    for ip in integration_points
-#        J = interpolate(el, "coordinates", ip.xi; derivative=true)
-#        push!(target, ip.weight*f(ip, args...)*det(J))
-#    end
-#    return sum(target)
-#end
 
 """
 This version returns a function which must be operated with element e
@@ -183,28 +157,6 @@ 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.
 """

src/types.jl

@@ -1,27 +1,6 @@
 # This file is a part of JuliaFEM.
 # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
 
-abstract Equation
-
-
-"""
-Integration point
-
-xi :: Array{Float64, 1}
-    (dimensionless) coordinates of integration point
-weight :: Float64
-    Integration weight
-attributes :: Dict{ASCIIString, Any}
-    This is used to save internal variables of IP needed e.g. for incremental
-    material models.
-"""
-type IntegrationPoint
-    xi :: Array{Float64, 1}
-    weight :: Float64
-    attributes :: Dict{ASCIIString, Any}
-end
-IntegrationPoint(xi, weight) = IntegrationPoint(xi, weight, Dict{ASCIIString, Any}())
-
 type Assembly
     # LHS
     I :: Array{Int64, 1}