Glassvision/JuliaFEM.jl
1
0
Fork 0
mirror of https://github.com/JuliaFEM/JuliaFEM.jl.git synced 2026-08-28 15:21:02 +00:00
Code Issues Packages Projects Releases Wiki Activity

fixed heat assembly + test. add several values to element at once by using update!(element, "field", 0 => 1.0, 1 => 2.0, ..., N => NN)

Browse Source
This commit is contained in:
Jukka Aho
2016-06-08 23:42:36 +03:00
parent c94a9e4407
commit fb402a7338
3 changed files with 53 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/elements.jl
+7 -5
View File
@@ -115,11 +115,13 @@ function update!(element::Element, field_name::ASCIIString, data::Dict)
element[field_name] = [data[i] for i in get_connectivity(element)]
end
function update!(element::Element, field_name::ASCIIString, data::Union{Real, Vector, Pair})
if haskey(element, field_name)
update!(element[field_name], data)
else
element[field_name] = data
function update!(element::Element, field_name::ASCIIString, datas::Union{Real, Vector, Pair}...)
for data in datas
if haskey(element, field_name)
update!(element[field_name], data)
else
element[field_name] = data
end
end
end
src/heat.jl
+14 -16
View File
@@ -3,18 +3,14 @@
# Heat problems
abstract HeatProblem <: AbstractProblem
function HeatProblem(dim::Int=1, elements=[])
return Problem{HeatProblem}(dim, elements)
type Heat <: FieldProblem
end
function get_unknown_field_name{P<:HeatProblem}(::Type{P})
function get_unknown_field_name(problem::Problem{Heat})
return "temperature"
end
function get_unknown_field_type{P<:HeatProblem}(::Type{P})
# scalar field
function get_unknown_field_type(problem::Problem{Heat})
return Float64
end
@@ -42,29 +38,31 @@ References
https://en.wikipedia.org/wiki/Heat_equation
"""
function assemble!(assembly::Assembly, problem::Problem{HeatProblem}, element::Element, time::Number)
function assemble!(assembly::Assembly, problem::Problem{Heat}, element::Element, time=0.0)
gdofs = get_gdofs(problem, element)
gdofs = get_gdofs(element, problem.dim)
for ip in get_integration_points(element)
w = ip.weight
J = get_jacobian(element, ip, time)
detJ = element(ip, time, Val{:detJ})
w = ip.weight*detJ
N = element(ip, time)
if haskey(element, "density")
rho = element("density", ip, time)
add!(assembly.mass_matrix, gdofs, gdofs, w*rho*N'*N*det(J))
add!(assembly.M, gdofs, gdofs, w*rho*N'*N)
end
if haskey(element, "temperature thermal conductivity")
dN = element(ip, time, Val{:grad})
dN = element(ip, time, Val{:Grad})
k = element("temperature thermal conductivity", ip, time)
add!(assembly.stiffness_matrix, gdofs, gdofs, w*k*dN'*dN*det(J))
add!(assembly.K, gdofs, gdofs, w*k*dN'*dN)
end
if haskey(element, "temperature load")
f = element("temperature load", ip, time)
add!(assembly.force_vector, gdofs, w*N'*f*det(J))
add!(assembly.f, gdofs, w*N'*f)
end
if haskey(element, "temperature flux")
g = element("temperature flux", ip, time)
add!(assembly.force_vector, gdofs, w*N'*g*norm(J))
add!(assembly.f, gdofs, w*N'*g)
end
end
end
test/test_heat.jl
+32 -35
View File
@@ -1,63 +1,60 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
# unit tests for heat equations
using JuliaFEM
using JuliaFEM.Test
module HeatTests # always wrap tests to module ending with "Tests"
using JuliaFEM.Test # always use JuliaFEM.Test, not Base.Test
@testset "test one element heat problem" begin
using JuliaFEM.Core: Seg2, Quad4, HeatProblem, assemble
function test_one_element() # always start test function with name test_
X = Dict{Int, Vector{Float64}}(
1 => [0.0,0.0],
2 => [1.0,0.0],
3 => [1.0,1.0],
4 => [0.0,1.0])
# volume element
element = Quad4([1, 2, 3, 4])
element = Element(Quad4, [1, 2, 3, 4])
element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
element["temperature thermal conductivity"] = 6.0
element["temperature load"] = [12.0, 12.0, 12.0, 12.0]
element["density"] = 36.0
update!(element, "geometry", X)
update!(element, "temperature thermal conductivity", 6.0)
update!(element, "temperature load", [12.0, 12.0, 12.0, 12.0])
update!(element, "density", 36.0)
# boundary element
boundary_element = Seg2([1, 2])
boundary_element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
boundary_element = Element(Seg2, [1, 2])
update!(boundary_element, "geometry", X)
# linear ramp from 0 to 6 in time 0 to 1
boundary_element["temperature flux"] = (0.0 => 0.0, 1.0 => 6.0)
update!(boundary_element, "temperature flux", 0.0 => 0.0, 1.0 => 6.0)
problem = HeatProblem()
push!(problem, element)
push!(problem, boundary_element)
problem = Problem(Heat, "one element heat problem", 1)
push!(problem, element, boundary_element)
# Set constant source f=12 with k=6. Accurate solution is
# T=1 on free boundary, u(x,y) = -1/6*(1/2*f*x^2 - f*x)
assembly = assemble(problem, 0.0)
fdofs = [1, 2]
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
assemble!(problem, 0.0)
A = full(problem.assembly.K)
b = full(problem.assembly.f)
info("stiffness matrix = \n$(round(A, 3))")
@test isapprox(A, [
A_expected = [
4.0 -1.0 -2.0 -1.0
-1.0 4.0 -1.0 -2.0
-2.0 -1.0 4.0 -1.0
-1.0 -2.0 -1.0 4.0
])
-1.0 -2.0 -1.0 4.0]
@test isapprox(A[fdofs, fdofs] \ b[fdofs], [1.0, 1.0])
@test isapprox(A, A_expected)
free_dofs = [1, 2]
@test isapprox(A[free_dofs, free_dofs] \ b[free_dofs], [1.0, 1.0])
# Set constant flux g=6 on boundary. Accurate solution is
# u(x,y) = x which equals T=1 on boundary.
# at time t=1.0 all loads should be on.
assembly = assemble(problem, 1.0)
A = full(assembly.stiffness_matrix)
b = full(assembly.force_vector)
T = A[fdofs, fdofs] \ b[fdofs]
info("T = $T")
empty!(problem.assembly)
assemble!(problem, 1.0)
A = full(problem.assembly.K)
b = full(problem.assembly.f)
T = A[free_dofs, free_dofs] \ b[free_dofs]
@test isapprox(T, [2.0, 2.0])
end
end