diff --git a/src/elements.jl b/src/elements.jl index 11a2530..a6535f9 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -2,16 +2,115 @@ # License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md using FactCheck +using ForwardDiff abstract Element -export get_number_of_nodes +#= ELEMENT DEFINITIONS -get_number_of_nodes(el::Type{Element}) = -1 -get_element_dimension(el::Element) = -1 -get_dbasisdx(el::Element) = nothing -get_basis(el::Element) = nothing +Each element must have +1. Connectivity information. How element is connected to other elements. + This is typically node ids in Lagrange elements. +2. Ability to store fields, in array of shape dim × nnodes, where dim is + dimension of field and nnodes is number of nodes of element. Note that + this is always 2d array. +3. Default constructor which takes connectivity as argument. +4. Basis functions and derivative of basis functions. + +These rules probably will change, but there's a test_element function which +tests element and that it obeys current rules. If test_element passes, +everything should be ok. I use Quad4 as an example element here. +Several functions are inherited from Element abstract type: + +- get_connectivity +- get_number_of_basis_functions +- get_element_dimension * +- get_basis * +- get_dbasisdxi +- get_dbasisdX +- get_field +- set_field +- interpolate +- ... + +Which should work if element is defined following some rules. Functions marked with asterisk * are the ones which must necessarily to implement by your own. + +=# + +# These must be implemented for your own element +get_element_dimension(el::Element) = nothing +get_basis(el::Element, xi) = nothing +""" Return partial derivatives of basis using ForwardDiff """ +function get_dbasisdxi(el::Element, xi) + f(xi) = get_basis(el, xi) + ForwardDiff.jacobian(f, xi) +end +function get_number_of_basis_functions(el::Type{Element}) + Logging.info("You really should define get_number_of_basis_functions") + # this is hack, evaluate basis in some point and return length of vector. + bf = get_basis([0.0, 0.0, 0.0]) + length(bf) +end + +abstract CG <: Element # Lagrange (continous Galerkin) element family + +""" +4 node bilinear quadrangle element + +X = [-1.0 1.0 1.0 -1.0 + -1.0 -1.0 1.0 1.0] + +P(xi) = [1.0, xi[1], xi[2], xi[1]*xi[2]] + +dP(xi) = [0.0 1.0 0.0 xi[2] + 0.0 0.0 1.0 xi[1]] +""" +type Quad4 <: CG + connectivity :: Array{Int, 1} + fields :: Dict{Any, Any} +end + +""" Default contructor. """ +Quad4(connectivity) = Quad4(connectivity, Dict{Any, Any}()) + +""" Return number of basis functions of this element. """ +get_number_of_basis_functions(el::Type{Quad4}) + +""" Return element dimension (length of xi vector). """ +get_element_dimension(el::Type{Quad4}) = 2 + +""" Return basis functions for this element (xi dim = 2, functions = 4). """ +function get_basis(el::Quad4, xi) + [(1-xi[1])*(1-xi[2])/4 + (1+xi[1])*(1-xi[2])/4 + (1+xi[1])*(1+xi[2])/4 + (1-xi[1])*(1+xi[2])/4] +end + +""" Return partial derivatives of basis functions. """ +function get_dbasisdxi(el::Quad4, xi) + [-(1-xi[2])/4.0 -(1-xi[1])/4.0 + (1-xi[2])/4.0 -(1+xi[1])/4.0 + (1+xi[2])/4.0 (1+xi[1])/4.0 + -(1+xi[2])/4.0 (1-xi[1])/4.0] +end +# TODO: create_lagrange_element(:Quad4, X, P, dP) + +# Common element routines + +""" +Test routine for element. + +Parameters +---------- +eltype::Type{Element} + Element to test + +Raises +------ +This uses FactCheck and throws exception if element is not passing. +""" function test_element(eltype) local el n = get_number_of_nodes(eltype) @@ -26,6 +125,7 @@ function test_element(eltype) Logging.error("""Unable to create element with default constructor define function $eltype(connectivity) which initializes this element. """) + return false end dim = get_element_dimension(el) Logging.info("Element dimension: $dim") @@ -58,14 +158,16 @@ function test_element(eltype) Logging.info("Element $eltype passed tests.") end - """ -Get jacobian of element evaluated at point xi +Get jacobian of element evaluated at point ξ on element. + +Notes +----- +This function assumes that element has field :geometry defined. """ function get_jacobian(el::Element, xi) dbasisdxi = get_dbasisdxi(el, xi) X = get_field(el, :geometry) - #J = interpolate(X, dbasisdxi, xi)' J = X*dbasisdxi return J end @@ -80,31 +182,18 @@ function get_dbasisdX(el::Element, xi) dbasisdxi*inv(J) end -""" -Set field variable. -""" +""" Set field variable. """ function set_field(el::Element, field_name, field_value) el.fields[field_name] = field_value end -""" -Get field variable. -""" +""" Get field variable. """ 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::Union(ASCIIString, Symbol), xi::Array{Float64,1}) +""" Evaluate some field in point ξ on element using basis functions. """ +function interpolate(el::Element, field, xi) f = get_field(el, field) basis = get_basis(el, xi) dim, nnodes = size(f) @@ -119,10 +208,7 @@ function interpolate(el::Element, field::Union(ASCIIString, Symbol), xi::Array{F end end -### Lagrange family ### - -abstract CG <: Element # Lagrange element family - +#= """ Create new Lagrange element @@ -169,8 +255,11 @@ function create_lagrange_element(element_name, X, P, dP) end +=# + # 0d Lagrange elements +#= """ 1 node point element """ @@ -211,49 +300,10 @@ end #create_lagrange_element(:Seg3, X, P, dP) # 2d Lagrange elements - -""" -4 node bilinear quadrangle element -""" -type Quad4 <: CG - node_ids :: Array{Int, 1} - fields :: Dict{ASCIIString, Any} -end -function Quad4(node_ids) - fields = Dict{ASCIIString, Any}() - Quad4(node_ids, fields) -end -function get_basis(el::Quad4, xi) - [(1-xi[1])*(1-xi[2])/4 - (1+xi[1])*(1-xi[2])/4 - (1+xi[1])*(1+xi[2])/4 - (1-xi[1])*(1+xi[2])/4] -end -function get_dbasisdxi(el::Quad4, xi) - [-(1-xi[2])/4.0 -(1-xi[1])/4.0 - (1-xi[2])/4.0 -(1+xi[1])/4.0 - (1+xi[2])/4.0 (1+xi[1])/4.0 - -(1+xi[2])/4.0 (1-xi[1])/4.0] -end -#X = [ -# -1.0 -1.0 -# 1.0 -1.0 -# 1.0 1.0 -# -1.0 1.0] -#P = (xi) -> [ -# 1.0 -# xi[1] -# xi[2] -# xi[1]*xi[2]] -#dP = (xi) -> [ -# 0.0 0.0 -# 1.0 0.0 -# 0.0 1.0 -# xi[2] xi[1]] -#create_lagrange_element(:Quad4, X, P, dP) +=# # 3d Lagrange elements - +#= """ 10 node quadratic tethahedron """ @@ -261,6 +311,7 @@ type Tet10 <: CG node_ids :: Array{Int, 1} fields :: Dict{ASCIIString, Any} end +=# # X = [ # 0.0 0.0 0.0 # 1.0 0.0 0.0 diff --git a/src/equations.jl b/src/equations.jl index 04037a6..2f5ae63 100644 --- a/src/equations.jl +++ b/src/equations.jl @@ -21,7 +21,9 @@ type IntegrationPoint end IntegrationPoint(xi, weight) = IntegrationPoint(xi, weight, Dict{ASCIIString, Any}()) +has_lhs(eq::Equation) = false get_lhs(eq::Equation, xi) = nothing +has_rhs(eq::Equation) = false get_rhs(eq::Equation, xi) = nothing get_element(eq::Equation) = eq.element get_integration_points(eq::Equation) = eq.integration_points @@ -29,8 +31,8 @@ get_integration_points(eq::Equation) = eq.integration_points 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) +integrate_lhs(eq::Equation) = has_lhs(eq) ? integrate(eq, get_lhs) : nothing +integrate_rhs(eq::Equation) = has_rhs(eq) ? integrate(eq, get_rhs) : nothing """ Return determinant of Jacobian for numerical integration.