diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index b876d35..a447e6f 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/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, diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl index 45ace6a..a3a5f3b 100644 --- a/src/JuliaFEM.jl +++ b/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") diff --git a/src/elements.jl b/src/elements.jl index 72ba279..829fac0 100644 --- a/src/elements.jl +++ b/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 ### diff --git a/src/equations.jl b/src/equations.jl new file mode 100644 index 0000000..04037a6 --- /dev/null +++ b/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 + diff --git a/src/math.jl b/src/math.jl index 6a4f0e5..623b53d 100644 --- a/src/math.jl +++ b/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. """ diff --git a/src/types.jl b/src/types.jl index c75a291..66a79fc 100644 --- a/src/types.jl +++ b/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}