From e4d9cbb6c9f7d8cd2afedd0d34c62d59ac95bb23 Mon Sep 17 00:00:00 2001
From: Jukka Aho <jukka.aho@kapsi.fi>
Date: Mon, 14 Sep 2015 16:15:36 +0300
Subject: [PATCH] added heat equations.

---
 src/elements.jl       | 29 ++++++++++++++++------
 src/equations.jl      | 12 ++++++++-
 src/heat_equations.jl | 58 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 90 insertions(+), 9 deletions(-)
 create mode 100644 src/heat_equations.jl

diff --git a/src/elements.jl b/src/elements.jl
index c635fd9..c6c4213 100644
--- a/src/elements.jl
+++ b/src/elements.jl
@@ -135,7 +135,7 @@ end
 Assign Lagrange basis for element.
 """
 macro create_lagrange_basis(element_name, X, P)
-    
+
 #   Logging.debug("Creating Lagrange basis for element ", element_name, ". ")
     eltype = esc(element_name)
 
@@ -271,11 +271,14 @@ 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 = X*dbasisdxi
-    return J
+#function get_Jacobian(el::Element, xi)
+#    dbasisdxi = get_dbasisdxi(el, xi)
+#    X = get_field(el, :geometry)
+#    J = X*dbasisdxi
+#    return J
+#end
+function get_Jacobian(el::Element, xi, geometry_field=:geometry)
+    dinterpolate(el, geometry_field, xi)
 end
 
 """
@@ -316,6 +319,16 @@ 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[]
+end
+
+""" Push to existing field. """
+function push_field!(el::Element, field_name, field_value)
+    push!(el.fields[field_name], field_value)
+end
+
 """ Get field variable. """
 function get_field(el::Element, field_name)
     el.fields[field_name]
@@ -355,10 +368,10 @@ function interpolate(el::Element, field, xis::Array{Vector, 1})
     map(interpolate_, xis)
 end
 function dinterpolate(el::Element, field, xi::Vector)
-    field = get_field(el, field)
+    fld = 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)])
+    return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
 end
diff --git a/src/equations.jl b/src/equations.jl
index 1c46979..ba07347 100644
--- a/src/equations.jl
+++ b/src/equations.jl
@@ -57,7 +57,7 @@ Integrate f over element
 Parameters
 ----------
 eq::Equation
-    
+
 f::Function
     Function to integrate
 """
@@ -69,3 +69,13 @@ function integrate(eq::Equation, f::Function)
     return sum(target)
 end
 
+function get_global_dofs(eq::Equation)
+    eq.global_dofs
+end
+
+function set_global_dofs(eq::Equation, dofs)
+    eq.global_dofs = dofs
+end
+
+# Equations for heat problems
+include("heat_equations.jl")
diff --git a/src/heat_equations.jl b/src/heat_equations.jl
new file mode 100644
index 0000000..c09d93a
--- /dev/null
+++ b/src/heat_equations.jl
@@ -0,0 +1,58 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
+
+abstract Heat <: Equation
+
+"""
+Diffusive heat transfer for 4-node bilinear element.
+"""
+type DC2D4 <: Heat
+    element :: Quad4
+    integration_points :: Array{IntegrationPoint, 1}
+    global_dofs :: Array{Int64, 1}
+end
+function DC2D4(el::Quad4)
+    integration_points = [
+        IntegrationPoint(1.0/sqrt(3.0)*[-1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1, -1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[ 1,  1], 1.0),
+        IntegrationPoint(1.0/sqrt(3.0)*[-1,  1], 1.0)]
+    set_field(el, :temperature, zeros(2, 4))
+    DC2D4(el, integration_points, [])
+end
+function get_lhs(eq::DC2D4)
+    function get_lhs_(eq, ip)
+        el = get_element(eq)
+        xi = ip.xi
+        dNdX = get_dbasisdX(el, xi)'
+        hc = interpolate(el, :"temperature heat coefficient", xi)
+        return dNdX'*hc*dNdX
+    end
+    integrate(eq, get_lhs_)
+end
+
+
+"""
+Diffusive heat transfer for 2-node linear segment.
+"""
+type DC2D2 <: Heat
+    element :: Seg2
+    integration_points :: Array{IntegrationPoint, 1}
+    global_dofs :: Array{Int64, 1}
+end
+function DC2D2(el::Seg2)
+    integration_points = [IntegrationPoint([0.0], 1.0)]
+    set_field(el, :temperature, zeros(2, 1))
+    set_field(el, :"temperature flux", zeros(2, 1))
+    DC2D2(el, integration_points, [])
+end
+function get_rhs(eq::DC2D2)
+    function get_rhs_(eq, ip)
+        el = get_element(eq)
+        xi = ip.xi
+        N = get_basis(el, xi)
+        f = interpolate(el, :"temperature flux", xi)
+        return f*N
+    end
+    integrate(eq, get_rhs_)
+end
\ No newline at end of file