From fb402a733865e58b4858e7fdbcc4ceee082feb7d Mon Sep 17 00:00:00 2001 From: Jukka Aho Date: Wed, 8 Jun 2016 23:42:36 +0300 Subject: [PATCH] fixed heat assembly + test. add several values to element at once by using update!(element, "field", 0 => 1.0, 1 => 2.0, ..., N => NN) --- src/elements.jl | 12 +++++---- src/heat.jl | 30 ++++++++++----------- test/test_heat.jl | 67 ++++++++++++++++++++++------------------------- 3 files changed, 53 insertions(+), 56 deletions(-) diff --git a/src/elements.jl b/src/elements.jl index e433def..9f40439 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -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 diff --git a/src/heat.jl b/src/heat.jl index 4e34743..0bf04fb 100644 --- a/src/heat.jl +++ b/src/heat.jl @@ -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 diff --git a/test/test_heat.jl b/test/test_heat.jl index 2f43968..6d3921b 100644 --- a/test/test_heat.jl +++ b/test/test_heat.jl @@ -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