normal calculation

src/elements.jl

@@ -40,9 +40,11 @@ Which should work if element is defined following some rules. Functions marked w
 
 # These must be implemented for your own element
 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
 
 """
 Create new element with element_name to family element_family
@@ -52,7 +54,7 @@ Examples
 >>> @create_element(Seg2, CG, "2 node linear segment")
 """
 macro create_element(element_name, element_family, element_description)
-    print("Creating element ", element_name, ": ", element_description, "\n")
+#   Logging.debug("Creating element ", element_name, ": ", element_description, "\n")
     eltype = esc(element_name)
     elfam = esc(element_family)
     quote
@@ -108,6 +110,8 @@ End of example.
 
 =#
 
+### LAGRANGE ELEMENTS ###
+
 abstract CG <: Element # Lagrange (continous Galerkin) element family
 
 """
@@ -120,7 +124,7 @@ function calculate_lagrange_basis(P, X)
     for i=1:nbasis
         A[i,:] = P(X[:, i])
     end
-    println("Calculating inverse of A")
+#   Logging.debug("Calculating inverse of A")
     invA = inv(A)'
     basis(xi) = invA*P(xi)
     dbasisdxi = ForwardDiff.jacobian(basis)
@@ -132,7 +136,7 @@ Assign Lagrange basis for element.
 """
 macro create_lagrange_basis(element_name, X, P)
     
-    print("Creating Lagrange basis for element ", element_name, ". ")
+#   Logging.debug("Creating Lagrange basis for element ", element_name, ". ")
     eltype = esc(element_name)
 
     quote
@@ -142,16 +146,17 @@ macro create_lagrange_basis(element_name, X, P)
 
         dim = size($X, 1)
         nbasis = size($X, 2)
-        print("Number of basis functions: ", nbasis, ". ")
-        println("Element dimension: ", dim)
+#       Logging.debug("Number of basis functions: ", nbasis, ". ")
+#       Logging.debug("Element dimension: ", dim)
 
         get_number_of_basis_functions(el::Type{$(esc(element_name))}) = nbasis
+        get_number_of_basis_functions(el::$(esc(element_name))) = nbasis
         get_element_dimension(el::$(esc(element_name))) = dim
 
         basis, dbasisdxi = calculate_lagrange_basis($P, $X)
         get_basis(el::$eltype, xi) = basis(xi)
         get_dbasisdxi(el::$eltype, xi) = dbasisdxi(xi)
-        println("Element ", $element_name, " created.")
+#       Logging.debug("Element ", $element_name, " created.")
     end
 
 end
@@ -160,7 +165,6 @@ end
 
 @create_element(Point1, CG, "1 node point element")
 
-
 # 1d Lagrange elements
 
 @create_element(Seg2, CG, "2 node linear line element")
@@ -177,8 +181,8 @@ end
      -1.0 -1.0 1.0  1.0],
     (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2]])
 
-
 # 3d Lagrange elements
+
 @create_element(Tet10, CG, "10 node quadratic tetrahedron")
 @create_lagrange_basis(Tet10,
     [0.0 1.0 0.0 0.0 0.5 0.5 0.0 0.0 0.5 0.0
@@ -187,8 +191,12 @@ end
     (xi) -> [    1.0,   xi[1],       xi[2],       xi[3],     xi[1]^2,
              xi[2]^2, xi[3]^2, xi[1]*xi[2], xi[2]*xi[3], xi[3]*xi[1]])
 
-# Common element routines
 
+### HIERARCHICAL ELEMENTS ###
+
+include("hierarchical.jl")
+
+# Common element routines
 
 """
 Test routine for element.
@@ -257,13 +265,13 @@ function test_element(eltype)
 end
 
 """
-Get jacobian of element evaluated at point ξ on element.
+Get jacobian of element evaluated at point ξ on element in reference configuration.
 
 Notes
 -----
 This function assumes that element has field :geometry defined.
 """
-function get_jacobian(el::Element, xi)
+function get_Jacobian(el::Element, xi)
     dbasisdxi = get_dbasisdxi(el, xi)
     X = get_field(el, :geometry)
     J = X*dbasisdxi
@@ -271,15 +279,38 @@ function get_jacobian(el::Element, xi)
 end
 
 """
-Evaluate partial derivatives of basis function w.r.t
-material description X, i.e. dbasis/dX
+Get jacobian of element evaluated at point ξ on element in current configuration.
+
+Notes
+-----
+This function assumes that element has fields :geometry and :displacement defined.
+"""
+function get_jacobian(el::Element, xi)
+    dbasisdxi = get_dbasisdxi(el, xi)
+    X = get_field(el, :geometry)
+    u = get_field(el, :displacement)
+    j = (X+u)*dbasisdxi
+    return j
+end
+
+"""
+Evaluate partial derivatives of basis, dbasis/dX
 """
 function get_dbasisdX(el::Element, xi)
     dbasisdxi = get_dbasisdxi(el, xi)
-    J = get_jacobian(el, xi)
+    J = get_Jacobian(el, xi)
     dbasisdxi*inv(J)
 end
 
+"""
+Evaluate partial derivatives of basis, dbasis/dx
+"""
+function get_dbasisdx(el::Element, xi)
+    dbasisdxi = get_dbasisdxi(el, xi)
+    j = get_jacobian(el, xi)
+    dbasisdxi*inv(j)
+end
+
 """ Set field variable. """
 function set_field(el::Element, field_name, field_value)
     el.fields[field_name] = field_value
@@ -290,18 +321,44 @@ function get_field(el::Element, field_name)
     el.fields[field_name]
 end
 
-""" 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)
-    result = zeros(dim)
-    for i=1:nnodes
-        result += basis[i]*f[:,i]
-    end
-    if dim == 1
-        return result[1]
-    else
-        return result
-    end
+"""Evaluate some field in point ξ on element using basis functions.
+
+Parameters
+----------
+el :: Element
+field :: Union{ASCIIString, Symbol}
+xi :: Vector
+
+Returns
+-------
+Scalar, Vector, Tensor, depending on what is type of field to interpolate.
+
+Notes
+-----
+This has another version which returns multiple values for set of coordinates {ξᵢ}.
+dinterpolate returns derivatives.
+
+Examples
+--------
+>>> field = [1.0, 2.0, 3.0, 4.0]
+>>> set_field(el, :temperature, field)
+>>> interpolate(el, :temperature, [0.0, 0.0])
+15.0
+"""
+function interpolate(el::Element, field, xi::Vector)
+    field = get_field(el, field)
+    sum(get_basis(el, xi) .* field)
+end
+function interpolate(el::Element, field, xis::Array{Vector, 1})
+    field = get_field(el, field)
+    interpolate_(xi) = sum(get_basis(el, xi) .* field)
+    map(interpolate_, xis)
+end
+function dinterpolate(el::Element, field, xi::Vector)
+    field = get_field(el, field)
+    dbasis = get_dbasisdxi(el, xi)
+    if isa(dbasis, Vector)
+        return sum(dbasis .* field)
+    end
+    return sum([fld[i]*g[i,:] for i in 1:length(fld)])
 end

src/equations.jl

@@ -42,7 +42,7 @@ function get_detJ(eq::Equation, ip::IntegrationPoint)
     get_detJ(el, ip)
 end
 function get_detJ(el::Element, ip::IntegrationPoint)
-    J = get_jacobian(el, ip.xi)
+    J = get_Jacobian(el, ip.xi)
     n, m = size(J)
     if n != m # for manifolds
         return norm(J)

src/hierarchical.jl

@@ -3,17 +3,36 @@
 
 # some preliminary code for constructing hierarchical elements
 
+abstract Hierarchical <: Element
+
+bin(n, k) = prod([(n + 1 - i)/i for i=1:k])
+
 """
 Return Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]
 
-Maybe slow version?
+Parameters
+----------
+n :: Int
+    order of polynomial
+
+Returns
+-------
+function
+    Legendgre polynomial of order n in interval ξ ∈ [-1, 1]
 """
-function get_legendre_polynomial_2(n)
-    bin(n, k) = prod([(n + 1 - i)/i for i=1:k])
-    P(xi) = sum([2^n*xi.^k*bin(n, k)*bin(1/2*(n+k-1), n) for k=0:n])
+function get_legendre_polynomial(n::Int)
+    P(xi) = 2^n*sum([xi.^k*bin(n, k)*bin(1/2*(n+k-1), n) for k=0:n])
     P
 end
 
+"""
+Return derivative of Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]
+"""
+function get_legendre_polynomial_derivative(n::Int)
+    dP(xi) = 2^n*sum([k*xi.^(k-1)*bin(n, k)*bin(1/2*(n+k-1), n) for k=0:n])
+    dP
+end
+
 """
 Return Legendgre polynomial of order n to inverval ξ ∈ [1, 1].
 
@@ -32,7 +51,7 @@ Notes
 Uses Bonnet's recursion formula. See
 https://en.wikipedia.org/wiki/Legendre_polynomials
 """
-function get_legendre_polynomial(n)
+function get_legendre_polynomial_recursive(n)
     if n == 0
         P(xi) = 1
     elseif n == 1
@@ -48,11 +67,11 @@ end
 """
 Return derivative of Legendgre polynomial of order n to inverval ξ ∈  [-1, 1]
 """
-function get_legendre_polynomial_derivative(n)
+function get_legendre_polynomial_derivative_recursive(n)
     if n == 0
-        P(xi) = 0*xi
+        P(xi) = 0
     elseif n == 1
-        P(xi) = 0*xi + 1
+        P(xi) = 1
     else
         Pm1 = get_legendre_polynomial_derivative(n-1)
         Pm2 = get_legendre_polynomial_derivative(n-2)
@@ -83,9 +102,9 @@ Return derivative of hierarchical shape function of order N
 """
 function get_hierarchial_basis_derivative(n)
     if n == 1
-        dN(xi) = 1/2*(0*xi - 1)
+        dN(xi) = -1/2
     elseif n == 2
-        dN(xi) = 1/2*(0*xi + 1)
+        dN(xi) = 1/2
     else
         j = n-1
         Pj = get_legendre_polynomial_derivative(j)
@@ -95,3 +114,47 @@ function get_hierarchial_basis_derivative(n)
     return dN
 end
 
+"""
+Set degree of hierarchical element
+"""
+function set_degree(el::Hierarchical, degree)
+    el.degree = degree
+end
+
+"""
+Get degree of hierarchical element
+"""
+function get_degree(el::Hierarchical)
+    el.degree
+end
+
+"""
+Hierarchical 1d segment element.
+"""
+type PSeg <: Hierarchical
+    connectivity :: Array{Int, 1}
+    fields :: Dict{Any, Any}
+    degree :: Int
+end
+PSeg(connectivity) = PSeg(connectivity, Dict{Any,Any}(), 1)
+get_number_of_basis_functions(el::Type{PSeg}) = 2
+get_number_of_basis_functions(el::PSeg) = 2 + el.degree - 1
+get_element_dimension(el::Type{PSeg}) = 1
+function get_basis(el::PSeg, xi)
+    m = get_number_of_basis_functions(el)
+    out = zeros(m)
+    for n=1:m
+        N = get_hierarchial_basis(n)
+        out[n] = N(xi[1])
+    end
+    return out
+end
+function get_dbasisdxi(el::PSeg, xi)
+    m = get_number_of_basis_functions(el)
+    out = zeros(m)
+    for n=1:m
+        dN = get_hierarchial_basis_derivative(n)
+        out[n] = dN(xi[1])
+    end
+    return out
+end