diff --git a/src/JuliaFEM.jl b/src/JuliaFEM.jl
index 0446e0e..9d83e6c 100644
--- a/src/JuliaFEM.jl
+++ b/src/JuliaFEM.jl
@@ -177,7 +177,9 @@ export @timeit, print_timer
 # using AbaqusReader  # Consolidated into src/readers.jl
 # using AsterReader   # Consolidated into src/readers.jl
 
-# @reexport using HeatTransfer
+# Problem types
+include("problems_heat.jl")
+export Heat, PlaneHeat
 include("problems_elasticity.jl")
 export Elasticity
 include("materials_plasticity.jl")
diff --git a/src/graph/graph_ordering.jl b/src/graph/graph_ordering.jl
index 998de70..57cd931 100644
--- a/src/graph/graph_ordering.jl
+++ b/src/graph/graph_ordering.jl
@@ -5,11 +5,11 @@
 # Consolidated from GraphOrdering.jl
 
 struct GraphOrderingResult
-    perm :: Vector{Int}
-    invperm :: Vector{Int}
-    degrees :: Vector{Int}
-    edge :: Vector{Int}
-    dist :: Vector{Int}
+    perm::Vector{Int}
+    invperm::Vector{Int}
+    degrees::Vector{Int}
+    edge::Vector{Int}
+    dist::Vector{Int}
 end
 
 """
@@ -22,10 +22,10 @@ function bandwidth(G)
     bw = -1
     for (v, adj) in G
         for w in adj
-            bw = max(bw, abs(v-w))
+            bw = max(bw, abs(v - w))
         end
     end
-    return 2*bw + 1
+    return 2 * bw + 1
 end
 
 """
@@ -50,7 +50,7 @@ function symrcm(G, v)
     wrk = zeros(Int, nwrk)
     idx = 2
 
-    for i=1:n
+    for i = 1:n
 
         v = permutation[i]
         adj = G[v]
diff --git a/src/mesh/gmsh_wrapper.jl b/src/mesh/gmsh_wrapper.jl
new file mode 100644
index 0000000..46aeee2
--- /dev/null
+++ b/src/mesh/gmsh_wrapper.jl
@@ -0,0 +1,58 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/Gmsh.jl/blob/master/LICENSE
+#
+# Gmsh wrapper - consolidated from Gmsh.jl
+# Provides convenient access to Gmsh API via gmsh_jll
+
+import gmsh_jll
+include(gmsh_jll.gmsh_api)
+import .gmsh
+
+"""
+    gmsh_initialize(argv=String[]; finalize_atexit=true)
+
+Wrapper around `gmsh.initialize` which make sure to only call it if `gmsh` is not already
+initialized. Return `true` if `gmsh.initialize` was called, and `false` if `gmsh` was
+already initialized.
+
+The argument vector `argv` is passed to `gmsh.initialize`. `argv` can be used to pass
+command line options to Gmsh, see [Gmsh documentation for more
+details](https://gmsh.info/doc/texinfo/gmsh.html#index-Command_002dline-options). Note that
+this wrapper prepends the program name to `argv` since Gmsh expects that to be the first
+entry.
+
+If `finalize_atexit` is `true` a Julia exit hook is added, which calls `finalize()`.
+
+**Example**
+```julia
+Gmsh.initialize(["-v", "0"]) # initialize with decreased verbosity
+```
+"""
+function gmsh_initialize(argv=String[]; finalize_atexit=true)
+    if Bool(gmsh.isInitialized())
+        return false
+    end
+    # Prepend a dummy program name in case argv only contains options
+    # see https://gitlab.onelab.info/gmsh/gmsh/-/issues/2112
+    if length(argv) > 0 && startswith(first(argv), "-")
+        argv = pushfirst!(copy(argv), "gmsh")
+    end
+    gmsh.initialize(argv)
+    if finalize_atexit
+        atexit(finalize)
+    end
+    return true
+
+
+"""
+    Gmsh.finalize()
+
+Wrapper around `gmsh.finalize` which make sure to only call it if `gmsh` is initialized.
+"""
+function gmsh_finalize()
+    if Bool(gmsh.isInitialized())
+        gmsh.finalize()
+    end
+
+
+
diff --git a/src/problems_heat.jl b/src/problems_heat.jl
new file mode 100644
index 0000000..6674ccb
--- /dev/null
+++ b/src/problems_heat.jl
@@ -0,0 +1,117 @@
+# This file is a part of JuliaFEM.
+# License is MIT: see https://github.com/JuliaFEM/HeatTransfer.jl/blob/master/LICENSE
+#
+# Heat transfer problem types - consolidated from HeatTransfer.jl
+
+"""
+    Heat
+
+3D heat transfer analysis for JuliaFEM.
+
+# Fields used in formulation
+
+- `thermal conductivity`
+- `heat source`
+- `heat flux`
+- `external temperature`
+- `heat transfer coefficient`
+
+# References
+
+- https://en.wikipedia.org/wiki/Heat_equation
+- https://en.wikipedia.org/wiki/Heat_capacity
+- https://en.wikipedia.org/wiki/Heat_flux
+- 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
+"""
+
+struct PlaneHeat <: FieldProblem end
+struct Heat <: FieldProblem end
+
+get_unknown_field_name(::PlaneHeat) = "temperature"
+get_unknown_field_name(::Heat) = "temperature"
+
+function assemble_elements!(problem::Problem{P}, assembly::Assembly,
+    elements::Vector{Element{M,B}}, time::Float64) where
+{M,B,P<:Union{PlaneHeat,Heat}}
+
+    bi = BasisInfo(B)
+    ndofs = length(bi)
+    Ke = zeros(ndofs, ndofs)
+    fe = zeros(ndofs)
+
+    for element in elements
+        fill!(Ke, 0.0)
+        fill!(fe, 0.0)
+        for ip in get_integration_points(element)
+            J, detJ, N, dN = element_info!(bi, element, ip, time)
+            s = ip.weight * detJ
+            k = element("thermal conductivity", ip, time)
+            Ke += s * k * dN' * dN
+            if haskey(element, "heat source")
+                f = element("heat source", ip, time)
+                fe += s * N' * f
+            end
+        end
+        if haskey(element, "temperature")
+            T = [element("temperature", time)...]
+            fe -= Ke * T
+        end
+        gdofs = get_gdofs(problem, element)
+        add!(assembly.K, gdofs, gdofs, Ke)
+        add!(assembly.f, gdofs, fe)
+    end
+
+end
+
+function assemble_elements!(problem::Problem{PlaneHeat}, assembly::Assembly,
+    elements::Vector{Element{M,B}}, time::Float64) where
+{M,B<:Union{Seg2,Seg3}}
+    return assemble_boundary_elements!(problem, assembly, elements, time)
+end
+
+function assemble_elements!(problem::Problem{Heat}, assembly::Assembly,
+                                    elements::Vector{Element{M,B}}, time::Float64) where
+                                    {M,B<:Union{Tri3,Quad4,Tri6,Quad8,Quad9}}
+    return assemble_boundary_elements!(problem, assembly, elements, time)
+end
+
+function assemble_boundary_elements!(problem::Problem, assembly::Assembly,
+                                     elements::Vector{Element{M,B}}, time::Float64) where {M,B}
+
+    bi = BasisInfo(B)
+    ndofs = length(bi)
+    Ke = zeros(ndofs, ndofs)
+    fe = zeros(ndofs)
+
+    for element in elements
+        fill!(fe, 0.0)
+        fill!(Ke, 0.0)
+        for ip in get_integration_points(element, 2)
+            J, detJ, N, dN = element_info!(bi, element, ip, time)
+            s = ip.weight * detJ
+            if haskey(element, "heat flux")
+                g = element("heat flux", ip, time)
+                fe += s * N' * g
+            end
+            if haskey(element, "heat transfer coefficient") && haskey(element, "external temperature")
+                h = element("heat transfer coefficient", ip, time)
+                Tu = element("external temperature", ip, time)
+                Ke += s * h * N' * N
+                fe += s * N' * h * Tu
+            end
+        end
+        if haskey(element, "temperature")
+            T = [element("temperature", time)...]
+            fe -= Ke * T
+        end
+        gdofs = get_gdofs(problem, element)
+        add!(assembly.K, gdofs, gdofs, Ke)
+        add!(assembly.f, gdofs, fe)
+    end
+
+end
+
+export Heat, PlaneHeat
diff --git a/src/readers/parse_model.jl b/src/readers/parse_model.jl
index 825f97c..69c8ac4 100644
--- a/src/readers/parse_model.jl
+++ b/src/readers/parse_model.jl
@@ -14,78 +14,78 @@ abstract type AbstractBoundaryCondition end
 abstract type AbstractOutputRequest end
 
 mutable struct Mesh
-    nodes :: Dict{Int, Vector{Float64}}
-    node_sets :: Dict{String, Vector{Int}}
-    elements :: Dict{Int, Vector{Int}}
-    element_types :: Dict{Int, Symbol}
-    element_sets :: Dict{String, Vector{Int}}
-    surface_sets :: Dict{String, Vector{Tuple{Int, Symbol}}}
-    surface_types :: Dict{String, Symbol}
+    nodes::Dict{Int,Vector{Float64}}
+    node_sets::Dict{String,Vector{Int}}
+    elements::Dict{Int,Vector{Int}}
+    element_types::Dict{Int,Symbol}
+    element_sets::Dict{String,Vector{Int}}
+    surface_sets::Dict{String,Vector{Tuple{Int,Symbol}}}
+    surface_types::Dict{String,Symbol}
 end
 
-function Mesh(d::Dict{String, Dict})
+function Mesh(d::Dict{String,Dict})
     return Mesh(d["nodes"], d["node_sets"], d["elements"],
-                d["element_types"], d["element_sets"],
-                d["surface_sets"], d["surface_types"])
+        d["element_types"], d["element_sets"],
+        d["surface_sets"], d["surface_types"])
 end
 
 mutable struct Model
-    path :: String
-    name :: String
-    mesh :: Mesh
-    materials :: Dict{Symbol, AbstractMaterial}
-    properties :: Vector{AbstractProperty}
-    boundary_conditions :: Vector{AbstractBoundaryCondition}
-    steps :: Vector{AbstractStep}
+    path::String
+    name::String
+    mesh::Mesh
+    materials::Dict{Symbol,AbstractMaterial}
+    properties::Vector{AbstractProperty}
+    boundary_conditions::Vector{AbstractBoundaryCondition}
+    steps::Vector{AbstractStep}
 end
 
 mutable struct SolidSection <: AbstractProperty
-    element_set :: Symbol
-    material_name :: Symbol
+    element_set::Symbol
+    material_name::Symbol
 end
 
 mutable struct Material <: AbstractMaterial
-    name :: Symbol
-    properties :: Vector{AbstractMaterialProperty}
+    name::Symbol
+    properties::Vector{AbstractMaterialProperty}
 end
 
 mutable struct Elastic <: AbstractMaterialProperty
-    E :: Float64
-    nu :: Float64
+    E::Float64
+    nu::Float64
 end
 
 mutable struct Step <: AbstractStep
-    kind :: Union{Symbol, Nothing} # STATIC, ... (was Nullable{Symbol} in Julia 0.x)
-    boundary_conditions :: Vector{AbstractBoundaryCondition}
-    output_requests :: Vector{AbstractOutputRequest}
+    kind::Union{Symbol,Nothing} # STATIC, ... (was Nullable{Symbol} in Julia 0.x)
+    boundary_conditions::Vector{AbstractBoundaryCondition}
+    output_requests::Vector{AbstractOutputRequest}
 end
 
 mutable struct BoundaryCondition <: AbstractBoundaryCondition
-    kind :: Symbol # BOUNDARY, CLOAD, DLOAD, DSLOAD, ...
-    data :: Vector
-    options :: Dict
+    kind::Symbol # BOUNDARY, CLOAD, DLOAD, DSLOAD, ...
+    data::Vector
+    options::Dict
 end
 
 mutable struct OutputRequest <: AbstractOutputRequest
-    kind :: Symbol # NODE, EL, SECTION, ...
-    data :: Vector
-    options :: Dict
-    target :: Symbol # PRINT, FILE
+    kind::Symbol # NODE, EL, SECTION, ...
+    data::Vector
+    options::Dict
+    target::Symbol # PRINT, FILE
 end
 
 ### Utility functions to parse ABAQUS .inp file to data model
 
 mutable struct Keyword
-    name :: String
-    options :: Vector{Union{String, Pair}}
+    name::String
+    options::Vector{Union{String,Pair}}
 end
 
 mutable struct AbaqusReaderState
-    section :: Union{Keyword, Nothing}  # was Nullable{Keyword}
-    material :: Union{AbstractMaterial, Nothing}  # was Nullable
-    property :: Union{AbstractProperty, Nothing}  # was Nullable
-    step :: Union{AbstractStep, Nothing}  # was Nullable
-    data :: Vector{String}
+    section::Union{Keyword,Nothing}  # was Nullable{Keyword}
+    material::Union{AbstractMaterial,Nothing}  # was Nullable
+    property::Union{AbstractProperty,Nothing}  # was Nullable
+    step::Union{AbstractStep,Nothing}  # was Nullable
+    data::Vector{String}
 end
 
 function get_data(state::AbaqusReaderState)
@@ -161,7 +161,7 @@ function maybe_close_section!(model, state)
     isnull(state.section) && return
     section_name = state.section.name
     @debug("Close section: $section_name")
-    args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}}
+    args = Tuple{Model,AbaqusReaderState,Type{Val{Symbol(section_name)}}}
     if hasmethod(close_section!, args)
         close_section!(model, state, Val{Symbol(section_name)})
     else
@@ -175,7 +175,7 @@ function maybe_open_section!(model, state)
     section_name = state.section.name
     section_options = state.section.options
     @debug("New section: $section_name with options $section_options")
-    args = Tuple{Model, AbaqusReaderState, Type{Val{Symbol(section_name)}}}
+    args = Tuple{Model,AbaqusReaderState,Type{Val{Symbol(section_name)}}}
     if hasmethod(open_section!, args)
         open_section!(model, state, Val{Symbol(section_name)})
     else
@@ -259,12 +259,12 @@ BOUNDARY = register_abaqus_keyword("BOUNDARY")
 CLOAD = register_abaqus_keyword("CLOAD")
 DLOAD = register_abaqus_keyword("DLOAD")
 DSLOAD = register_abaqus_keyword("DSLOAD")
-const BOUNDARY_CONDITIONS = Union{BOUNDARY, CLOAD, DLOAD, DSLOAD}
+const BOUNDARY_CONDITIONS = Union{BOUNDARY,CLOAD,DLOAD,DSLOAD}
 
 NODE_PRINT = register_abaqus_keyword("NODE PRINT")
 EL_PRINT = register_abaqus_keyword("EL PRINT")
 SECTION_PRINT = register_abaqus_keyword("SECTION PRINT")
-const OUTPUT_REQUESTS = Union{NODE_PRINT, EL_PRINT, SECTION_PRINT}
+const OUTPUT_REQUESTS = Union{NODE_PRINT,EL_PRINT,SECTION_PRINT}
 
 ## Properties