heat problem

src/JuliaFEM.jl

@@ -9,8 +9,10 @@ module JuliaFEM
 importall Base
 using ForwardDiff
 using JLD
-autodiffcache = ForwardDiffCache()
-# export derivative, jacobian, hessian
+
+#Grad = Val{:Grad}
+#detJ = Val{:detJ}
+#export Grad, detJ
 
 include("common.jl")
 
@@ -169,4 +171,5 @@ module Interfaces
 include("interfaces.jl")
 end
 
+
 end # module

src/assembly.jl

@@ -44,17 +44,30 @@ end
 function assemble_posthook!
 end
 
-function assemble!(problem::Problem, time::Real)
+function assemble!(problem::Problem, time=0.0; auto_initialize=true)
     if !isempty(problem.assembly)
-        warn("problem.assembly is not empty and assembling, are you sure you know what are you doing?")
+        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
     end
-    if method_exists(assemble_prehook!, Tuple{typeof(problem), Real})
+    if isempty(problem.elements)
+        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
+    else
+        first_element = first(problem.elements)
+        unknown_field_name = get_unknown_field_name(problem)
+        if !haskey(first_element, unknown_field_name)
+            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
+            if auto_initialize
+                info("Initializing problem $(problem.name) at time $time automatically.")
+                initialize!(problem, time)
+            end
+        end
+    end
+    if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64})
         assemble_prehook!(problem, time)
     end
     for element in get_elements(problem)
         assemble!(problem.assembly, problem, element, time)
     end
-    if method_exists(assemble_posthook!, Tuple{typeof(problem), Real})
+    if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
         assemble_posthook!(problem, time)
     end
 end

src/elements.jl

@@ -29,10 +29,18 @@ function setindex!(element::Element, data, field_name::ASCIIString)
     element.fields[field_name] = Field(data)
 end
 
+function call(element::Element, field_name::ASCIIString)
+    return element[field_name]
+end
+
 function call(element::Element, field_name::ASCIIString, time)
     return element[field_name](time)
 end
 
+function last(element::Element, field_name::ASCIIString)
+    return last(element[field_name])
+end
+
 function call(element::Element, ip, time)
     return get_basis(element, ip, time)
 end
@@ -80,6 +88,10 @@ function call(element::Element, field::DCTI, ip, time::Float64)
     return field.data
 end
 
+function call(element::Element, field::DCTV, ip, time::Float64)
+    return field(time).data
+end
+
 function call(element::Element, field::CVTV, ip, time::Float64)
     return field(ip, time)
 end
@@ -89,7 +101,9 @@ function call(element::Element, field::Field, ip, time::Float64)
     basis = element(ip, time)
     n = length(element)
     m = length(field_)
-    @assert n == m
+    if n != m
+        error("Error when trying to interpolate field $field at coords $ip and time $time: element length is $n and field length is $m, f = Nᵢfᵢ makes no sense!")
+    end
     return sum([field_[i]*basis[i] for i=1:n])
 end
 
@@ -118,7 +132,7 @@ function update!{K,V}(element::Element, field_name::ASCIIString, data::Pair{Floa
     update!(element, field_name, time => element_data)
 end
 
-function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair{Float64, Union{Real, Vector{Any}}}}...)
+function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair{Float64, Union{Float64, Real, Vector{Any}}}}...)
     for data in datas
         if haskey(element, field_name)
             update!(element[field_name], data)
@@ -132,6 +146,12 @@ function update!(element::Element, field_name::ASCIIString, datas::Union{Real, V
     end
 end
 
+function update!(element::Element, field_name::ASCIIString, datas::Pair...)
+    for data in datas
+        update!(element, field_name, data)
+    end
+end
+
 function update!(element::Element, field_name::ASCIIString, data::Pair{Float64, Vector{Any}})
     if haskey(element, field_name)
         update!(element[field_name], data)
@@ -200,11 +220,12 @@ function update!(elements::Vector, field_name::ASCIIString, data)
     end
 end
 
+dbasis_cache = ForwardDiff.jacobian
 """ Evaluate partial derivatives of basis functions using ForwardDiff. """
 function get_dbasis(element::Element, ip, time)
     xi = isa(ip, IP) ? ip.coords : ip
     basis(xi) = vec(get_basis(element, xi, time))
-    return ForwardDiff.jacobian(basis, xi, cache=autodiffcache)'
+    return ForwardDiff.jacobian(basis, xi)'
 end
 
 """ Check existence of field. """

src/heat.jl

@@ -22,23 +22,59 @@ where
     f = temperature load                    defined on volume elements
     g = temperature flux                    defined on boundary elements
 
+Parameters
+----------
+temperature thermal conductivity
+temperature load
+temperature flux
+
+thermal conductivity
+heat source
+heat flux
+heat transfer coefficient
+external temperature
+
+
+Formulations
+------------
+1D, 2D, 3D
+
 References
 ----------
 https://en.wikipedia.org/wiki/Heat_equation
-
+https://en.wikipedia.org/wiki/Heat_capacity
+https://en.wikipedia.org/wiki/Thermal_conduction
+https://en.wikipedia.org/wiki/Thermal_conductivity
+https://en.wikipedia.org/wiki/Thermal_diffusivity
+https://en.wikipedia.org/wiki/Volumetric_heat_capacity
 """
 type Heat <: FieldProblem
+    formulation :: ASCIIString
+end
+
+function Heat()
+    return Heat("3D")
 end
 
 function get_unknown_field_name(problem::Problem{Heat})
     return "temperature"
 end
 
-function get_unknown_field_type(problem::Problem{Heat})
-    return Float64
+function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
+    formulation = Val{Symbol(problem.properties.formulation)}
+    assemble!(assembly, problem, element, time, formulation)
 end
 
-function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
+# 3d heat problems
+
+function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
+    info("Unknown element type $E for 3d heat problem!")
+end
+
+typealias Heat3DVolumeElements Union{Tet4, Tet10, Hex8}
+typealias Heat3DSurfaceElements Union{Tri3, Tri6, Quad4}
+
+function assemble!{E<:Heat3DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
     gdofs = get_gdofs(problem, element)
     field_name = get_unknown_field_name(problem)
     nnodes = length(element)
@@ -57,6 +93,85 @@ function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element,
             f = element("$field_name load", ip, time)
             fq += w*N'*f
         end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+function assemble!{E<:Heat3DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("3D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        if haskey(element, "$field_name flux")
+            q = element("$field_name flux", ip, time)
+            fq += w*N'*q
+        end
+        if haskey(element, "$field_name heat transfer coefficient")
+            h = element("$field_name heat transfer coefficient", ip, time)
+            Tu = element("$field_name external temperature", ip, time)
+            K += w*h*N'*N
+            fq += w*N'*h*Tu
+        end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+# 2d heat problems
+
+function assemble!{E}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    info("Unknown element type $E for 2d heat problem!")
+end
+
+typealias Heat2DVolumeElements Union{Tri3, Tri6, Quad4}
+typealias Heat2DSurfaceElements Union{Seg2, Seg3}
+
+function assemble!{E<:Heat2DVolumeElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
+        if haskey(element, "$field_name thermal conductivity")
+            dN = element(ip, time, Val{:Grad})
+            k = element("$field_name thermal conductivity", ip, time)
+            K += w*k*dN'*dN
+        end
+        if haskey(element, "$field_name load")
+            f = element("$field_name load", ip, time)
+            fq += w*N'*f
+        end
+    end
+    T = vec(element[field_name](time))
+    fq -= K*T
+    add!(assembly.K, gdofs, gdofs, K)
+    add!(assembly.f, gdofs, fq)
+end
+
+function assemble!{E<:Heat2DSurfaceElements}(assembly::Assembly, problem::Problem{Heat}, element::Element{E}, time, ::Type{Val{Symbol("2D")}})
+    gdofs = get_gdofs(problem, element)
+    field_name = get_unknown_field_name(problem)
+    nnodes = length(element)
+    K = zeros(nnodes, nnodes)
+    fq = zeros(nnodes)
+    for ip in get_integration_points(element)
+        detJ = element(ip, time, Val{:detJ})
+        w = ip.weight*detJ
+        N = element(ip, time)
         if haskey(element, "$field_name flux")
             g = element("$field_name flux", ip, time)
             fq += w*N'*g

src/lagrange_macro.jl

@@ -95,17 +95,36 @@ end
     (xi) -> [1.0, xi[1], xi[2], xi[1]^2, xi[2]^2, xi[1]*xi[2]])
 
 @create_lagrange_element(Quad4, "4 node bilinear quadrangle element",
-    [-1.0  1.0 1.0 -1.0
-     -1.0 -1.0 1.0  1.0],
+    [-1.0  1.0  1.0 -1.0
+     -1.0 -1.0  1.0  1.0],
     (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2]])
 
+@create_lagrange_element(Quad9, "9 node bilinear quadrangle element",
+    [-1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0
+     -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0],
+    (xi) -> [1.0,     xi[1],   xi[2],         xi[1]*xi[2],
+             xi[1]^2, xi[2]^2, xi[1]^2*xi[2], xi[1]*xi[2]^2])
+
 # 3d Lagrange elements
 
 @create_lagrange_element(Hex8, "8 node hexahedra",
     [-1.0  1.0  1.0 -1.0 -1.0  1.0 1.0 -1.0
      -1.0 -1.0  1.0  1.0 -1.0 -1.0 1.0  1.0
      -1.0 -1.0 -1.0 -1.0  1.0  1.0 1.0  1.0],
-    (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2], xi[3], xi[1]*xi[3], xi[2]*xi[3], xi[1]*xi[2]*xi[3]])
+    (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2], xi[3],
+             xi[1]*xi[3], xi[2]*xi[3], xi[1]*xi[2]*xi[3]])
+
+#=
+@create_lagrange_element(Hex20, "20 node hexahedra",
+    [
+        -1.0  1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0
+        -1.0 -1.0  1.0  1.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0 -1.0 -1.0  1.0  1.0 -1.0  0.0  1.0  0.0
+        -1.0 -1.0 -1.0 -1.0  1.0  1.0  1.0  1.0 -1.0 -1.0 -1.0 -1.0  0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0
+    ],
+    (xi) -> [1.0, xi[1], xi[2], xi[1]*xi[2], xi[3], xi[1]*xi[3], xi[2]*xi[3], xi[1]*xi[2]*xi[3]
+             x[1]^2, 
+    ])
+=#
 
 @create_lagrange_element(Tet4, "4 node tetrahedron",
     [0.0 1.0 0.0 0.0

src/preprocess_aster_reader.jl

@@ -309,21 +309,17 @@ function get_element_sets(med::MEDFile, mesh_name)
     return es
 end
 
-# hex8 nodes rotating cw first in yz plane then x+1
-
 global const med_elmap = Dict{Symbol, Vector{Int}}(
+    :PO1 => [1],
     :SE2 => [1, 2],
     :SE3 => [1, 2, 3],
     :TR3 => [1, 2, 3],
-    :TR6 => [1, 2, 3, 4, 5, 6],
     :QU4 => [1, 2, 3, 4],
-    :HE8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
     :TE4 => [3, 2, 1, 4],
-    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8],
-    :PO1 => [1]
-#   :T10 => [3, 4, 1, 2, 10, 8, 7, 6, 9, 5]
-#   :T10 => [5, 9, 6, 7, 8, 10, 2, 1, 4, 3]
-)
+    :TR6 => [1, 2, 3, 4, 5, 6],
+    :QU8 => [1, 2, 3, 4, 5, 6, 7, 8],
+    :HE8 => [4, 8, 7, 3, 1, 5, 6, 2],  # ..?
+    :T10 => [3, 2, 1, 4, 6, 5, 7, 10, 9, 8])
 
 function get_connectivity(med::MEDFile, elsets, mesh_name)
     elsets[0] = :OTHER

src/problems.jl

@@ -145,7 +145,7 @@ end
     take last known value and set it as a initial quess for next
     time increment.
 """
-function initialize!(problem::Problem, time::Real)
+function initialize!(problem::Problem, time=0.0)
     field_name = get_unknown_field_name(problem)
     field_dim = get_unknown_field_dimension(problem)
     for element in get_elements(problem)
@@ -289,8 +289,12 @@ function get_parent_field_name{P<:BoundaryProblem}(problem::Problem{P})
     return problem.parent_field_name
 end
 
-function push!(problem::Problem, element)
-    push!(problem.elements, element)
+function push!(problem::Problem, elements...)
+    push!(problem.elements, elements...)
+end
+
+function push!(problem::Problem, elements::Vector)
+    push!(problem.elements, elements...)
 end
 
 function get_gdofs(element::Element, dim::Int)

src/solvers.jl

@@ -511,13 +511,18 @@ function call(solver::Solver{Linear}; F=nothing, show_info=true, return_factoriz
 end
 
 """ Convenience function to call linear solver. """
-function LinearSolver(problems...)
+function LinearSolver(problems::Problem...)
     solver = Solver(Linear, "default linear solver")
     if length(problems) != 0
         push!(solver, problems...)
     end
     return solver
 end
+function LinearSolver(name::ASCIIString, problems::Problem...)
+    solver = LinearSolver(problems...)
+    solver.name = name
+    return solver
+end
 
 ### End of linear quasistatic solver