diff --git a/notebooks/2015-08-29-developing-juliafem.ipynb b/notebooks/2015-08-29-developing-juliafem.ipynb index b2e99d7..258c2ae 100644 --- a/notebooks/2015-08-29-developing-juliafem.ipynb +++ b/notebooks/2015-08-29-developing-juliafem.ipynb @@ -13,22 +13,11 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": { "collapsed": false }, - "outputs": [ - { - "data": { - "text/plain": [ - "Logger(root,DEBUG,PipeEndpoint(open, 0 bytes waiting),root)" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "using Logging\n", "using FactCheck\n", @@ -60,20 +49,20 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "using JuliaFEM: Element" + "using JuliaFEM: Element, Field, Basis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Here's one basic implementation. The actual element:" + "Here's the implementation for element:" ] }, { @@ -86,7 +75,8 @@ "source": [ "type MyQuad4 <: Element\n", " connectivity :: Array{Int, 1}\n", - " fields :: Dict{Any, Any}\n", + " basis :: Basis\n", + " fields :: Dict{Symbol, Array{Field, 1}}\n", "end" ] }, @@ -94,7 +84,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Default constructor, providing connectivity data needed in assembly" + "Default constructor takes only connectivity information as input argument. Extra arguments may be passed using `args...`" ] }, { @@ -116,7 +106,14 @@ } ], "source": [ - "MyQuad4(connectivity) = MyQuad4(connectivity, Dict{Any, Any}())" + "function MyQuad4(connectivity, args...)\n", + " h = Basis((xi) ->\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", + " MyQuad4(connectivity, h, Dict())\n", + "end" ] }, { @@ -136,7 +133,7 @@ { "data": { "text/plain": [ - "get_element_dimension (generic function with 7 methods)" + "get_element_dimension (generic function with 8 methods)" ] }, "execution_count": 5, @@ -146,48 +143,7 @@ ], "source": [ "JuliaFEM.get_number_of_basis_functions(el::Type{MyQuad4}) = 4\n", - "JuliaFEM.get_element_dimension(el::MyQuad4) = 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 7 methods)" - ] - }, - "execution_count": 6, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "function JuliaFEM.get_basis(el::MyQuad4, xi)\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::MyQuad4, xi)\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" + "JuliaFEM.get_element_dimension(el::Type{MyQuad4}) = 2" ] }, { @@ -199,7 +155,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": 6, "metadata": { "collapsed": false }, @@ -208,23 +164,21 @@ "name": "stderr", "output_type": "stream", "text": [ - "24-Sep 21:04:04:INFO:root:Testing element MyQuad4\n", - "24-Sep 21:04:04:INFO:root:number of basis functions in this element: 4\n", - "24-Sep 21:04:04:INFO:root:Initializing element\n", - "24-Sep 21:04:04:INFO:root:Element dimension: 2\n", - "24-Sep 21:04:04:INFO:root:Setting scalar field [1 2 3 4] to element.\n", - "24-Sep 21:04:04:INFO:root:Interpolating scalar field at [0.0,0.0]\n" + "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" ] }, { - "data": { - "text/plain": [ - "PipeEndpoint(open, 0 bytes waiting)" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" + "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", + "" + ] } ], "source": [ @@ -775,6 +729,15 @@ "end\n", "A, b" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/src/elements.jl b/src/elements.jl index 28cb335..23ed479 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -1,6 +1,13 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md +#= +Related notebooks +----------------- + +2015-08-29-developing-juliafem.ipynb +=# + using FactCheck using ForwardDiff @@ -8,6 +15,8 @@ abstract Element #= ELEMENT DEFINITIONS +TODO: rewrite instructions after notebook. + Each element must have 1. Connectivity information. How element is connected to other elements. @@ -83,7 +92,6 @@ End of example. get_number_of_basis_functions(el::Type{Element}) = nothing get_number_of_basis_functions(el::Element) = nothing get_element_dimension(el::Element) = nothing -get_basis(el::Element, xi) = nothing get_dbasisdxi(el::Element, xi) = nothing get_connectivity(el::Element) = el.connectivity @@ -128,41 +136,50 @@ function test_element(eltype) return false end - dim = get_element_dimension(el) + dim = get_element_dimension(eltype) Logging.info("Element dimension: $dim") @fact dim --> not(nothing) """ 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 + fld = Field(0.0, collect(1:n)) + Logging.info("Pushing scalar field $fld to element.") + new_field!(el, :field1) + push_field!(el, :field1, fld) + @fact el[:field1][1] --> fld + mid = zeros(dim) try - get_basis(el, zeros(dim)) + f = get_basis(el)(mid) catch Logging.error(""" Unable to evaluate basis, define function 'get_basis' for this element.""") end try - get_dbasisdxi(el, zeros(dim)) + get_dbasisdxi(el)(mid) 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("Interpolating scalar field at $mid") + f(field, xi, t) = el(xi)*el[field](t) + i = f(:field, mid, 0.0) Logging.info("Value: $i") Logging.info("Element $eltype passed tests.") end +""" +Get basis functions of element. +""" +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 jacobian of element evaluated at point ξ on element in reference configuration. @@ -221,22 +238,15 @@ function get_dbasisdx(el::Element, xi) dbasisdxi*inv(j) end - -""" Set field variable. """ -function set_field(el::Element, field_name, field_value) - el.fields[field_name] = field_value -end - - """ Create new empty field of some type. """ -function new_field!(el::Element, field_name, field_type) - el.fields[field_name] = field_type[] +function new_field!(el::Element, field_name) + el.fields[field_name] = Field[] end -""" Push to existing field. """ -function push_field!(el::Element, field_name, field_value) - push!(el.fields[field_name], field_value) +""" Push to existing set field of fields. """ +function push_field!(el::Element, field_name, field::Field) + push!(el.fields[field_name], field) end @@ -244,7 +254,9 @@ end function get_field(el::Element, field_name) el.fields[field_name] end - +function Base.getindex(el::Element, field_name) + el.fields[field_name] +end """ Evaluate some field in point ξ on element using basis functions. diff --git a/src/types.jl b/src/types.jl index eed1943..455baa3 100644 --- a/src/types.jl +++ b/src/types.jl @@ -49,9 +49,12 @@ Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld """ Interpolate a field from finite set of fields some time t ∈ R. """ -function call(fields :: Array{Field, 1}, t::Float64) +function call(fields :: Array{Field, 1}, t::Number) + if length(fields) == 0 + throw("Empty set of fields.") + end if t <= fields[1].time - return fields[1] + return Field(t, fields[1].values) end if t >= fields[end].time return fields[end]