diff --git a/src/FEMBase.jl b/src/FEMBase.jl index 8708666..b705c8c 100644 --- a/src/FEMBase.jl +++ b/src/FEMBase.jl @@ -22,7 +22,7 @@ end export info, debug using FEMBasis -using FEMBasis: AbstractBasis, jacobian +using FEMBasis: AbstractBasis, jacobian, get_reference_element_coordinates import FEMBasis: interpolate using FEMQuad: get_quadrature_points include("fields.jl") @@ -48,6 +48,7 @@ export SparseMatrixCOO, SparseVectorCOO, Node, BasisInfo, export is_field_problem, is_boundary_problem, get_elements, get_connectivity, get_parent_field_name, get_reference_coordinates, + get_reference_element_coordinates, get_assembly, get_nonzero_rows, get_nonzero_columns, eval_basis!, get_basis, get_dbasis, grad!, assemble_mass_matrix!, get_local_coordinates, inside, diff --git a/src/assembly.jl b/src/assembly.jl index 5e6f6a7..376dad3 100644 --- a/src/assembly.jl +++ b/src/assembly.jl @@ -11,73 +11,98 @@ function isapprox(a1::Assembly, a2::Assembly) return T end -function assemble_prehook!(::Problem, time) end +function assemble_prehook!{T<:Number}(::Problem, ::T) end -function assemble_posthook!(::Problem, time) end +function assemble_posthook!{T<:Number}(::Problem, ::T) end -function assemble!(problem::Problem, time=0.0; auto_initialize=true) - if !isempty(problem.assembly) - warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?") - end - 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 - assembly = get_assembly(problem) - assemble_prehook!(problem, time) - for (element_type, elements) in group_by_element_type(get_elements(problem)) - assemble!(problem, assembly, elements, time) - end - assemble_posthook!(problem, time) - return true -end - -function assemble!{P}(::Assembly, ::Problem{P}, element::Element, time::Float64) +# will be deprecated +function assemble!{P}(::Assembly, ::Problem{P}, ::Element, ::Any) warn("One must define assemble! function for problem of type $P. ", "Not doing anything.") return nothing end +# will be deprecated +function assemble!{P}(assembly::Assembly, problem::Problem{P}, + elements::Vector{Element}, time) + warn("This is default assemble! function. Decreased performance can be ", + "expected without preallocation of memory. One should implement ", + "`assemble_elements!(problem, assembly, elements, time)` function.") + for element in elements + assemble!(assembly, problem, element, time) + end + return nothing +end + """ - assemble!(problem, assembly, elements, time) + assemble_elements!(problem, assembly, elements, time) + +Assemble elements for problem. -This should be overridden with own custom assemble operation. +This should be overridden with own `assemble_elements!`-implementation. """ -function assemble!{E}(problem::Problem, assembly::Assembly, - elements::Vector{Element{E}}, time::Float64) +function assemble_elements!{E}(problem::Problem, assembly::Assembly, + elements::Vector{Element{E}}, time) elements2 = convert(Vector{Element}, elements) assemble!(assembly, problem, elements2, time) end -function assemble!{P}(assembly::Assembly, problem::Problem{P}, - elements::Vector{Element}, time) - warn("This is default assemble! function. Decreased performance can be ", - "expected without preallocation of memory.") - for element in elements - assemble!(assembly, problem, element, time) +function assemble!(problem::Problem, time) + + assemble_prehook!(problem, time) + elements = get_elements(problem) + assembly = get_assembly(problem) + + if !isempty(assembly) + warn("Assembling elements for problem $(problem.name): problem.assembly ", + "is not empty before assembling. This is probably causing unexpected ", + "results. To remove old assembly, use `empty!(problem.assembly)`") + assemble_posthook!(problem, time) + return nothing end + + if isempty(elements) + warn("Assembling elements for problem $(problem.name): problem.elements ", + "is empty, there is no elements in problem. Before assembling ", + "problem, elements must be added using ", + "`add_elements!(problem, elements)`.") + assemble_posthook!(problem, time) + return nothing + end + + first_element = first(elements) + unknown_field_name = get_unknown_field_name(problem) + if !haskey(first_element, unknown_field_name) + warn("Assembling elements for problem $(problem.name): seems that ", + "problem is uninitialized. To initialize problem, use ", + "`initialize!(problem, time)`.") + info("Initializing problem $(problem.name) at time $time automatically.") + initialize!(problem, time) + end + + for (element_type, elements) in group_by_element_type(elements) + assemble_elements!(problem, assembly, elements, time) + end + assemble_posthook!(problem, time) return nothing end -function assemble_mass_matrix!(problem::Problem, time) +function assemble!(problem::Problem) + warn("assemble!(problem) will be deprecated. Use assemble!(problem, time)") + assemble!(problem, 0.0) +end + +function assemble_mass_matrix!(problem::Problem, time::Float64) if !isempty(problem.assembly.M) - info("Mass matrix for $(problem.name) is already assembled, skipping assemble routine") - return + info("Mass matrix for $(problem.name) is already assembled, ", + "not assembling.") + return nothing end elements = get_elements(problem) for (element_type, elements) in group_by_element_type(get_elements(problem)) assemble_mass_matrix!(problem::Problem, elements, time) end - return + return nothing end function assemble_mass_matrix!{Basis}(problem::Problem, elements::Vector{Element{Basis}}, time) @@ -137,7 +162,7 @@ function assemble_mass_matrix!(problem::Problem, elements::Vector{Element{Tet10} -6 -4 -6 -4 16 16 8 16 32 16 -6 -6 -4 -4 8 16 16 16 16 32] - function is_CM(element::Element{Tet10}, X; rtol=1.0e-6) + function is_CM(::Element{Tet10}, X; rtol=1.0e-6) isapprox(X[5], 1/2*(X[1]+X[2]); rtol=rtol) || return false isapprox(X[6], 1/2*(X[2]+X[3]); rtol=rtol) || return false isapprox(X[7], 1/2*(X[3]+X[1]); rtol=rtol) || return false diff --git a/src/elements.jl b/src/elements.jl index 84b3ed3..757858d 100644 --- a/src/elements.jl +++ b/src/elements.jl @@ -149,14 +149,14 @@ function interpolate(element::Element, field_name::String, time::Float64) end end -function get_basis{B}(element::Element{B}, ip, time) +function get_basis{B}(element::Element{B}, ip, ::Any) T = typeof(first(ip)) N = zeros(T, 1, length(element)) eval_basis!(B, N, tuple(ip...)) return N end -function get_dbasis{B}(element::Element{B}, ip, time) +function get_dbasis{B}(element::Element{B}, ip, ::Any) T = typeof(first(ip)) dN = zeros(T, size(element)...) eval_dbasis!(B, dN, tuple(ip...)) @@ -221,39 +221,30 @@ end function interpolate(element::Element, field_name, ip, time) field = element[field_name] - return interpolate_(element, field, ip, time) + return interpolate_field(element, field, ip, time) end -@lintpragma("Ignore unused ip") -@lintpragma("Ignore unused element") -function interpolate_(element::Element, field::DCTI, ip, time::Float64) - return field.data +function interpolate_field{T<:Union{DCTI,DCTV}}(::Element, field::T, ::Any, time) + return interpolate_field(field, time) end -function interpolate_(element::Element, field::DCTV, ip, time::Float64) - return interpolate_(field, time) +function interpolate_field{T<:Union{DVTV,DVTI}}(element::Element, field::T, ip, time) + data = interpolate_field(field, time) + basis = get_basis(element, ip, time) + N = length(basis) + return sum(data[i]*basis[i] for i=1:N) end -function interpolate_(element::Element, field::CVTV, ip, time::Float64) - return field(ip, time) +function interpolate_field{T<:Union{DVTVd,DVTId}}(element::Element, field::T, ip, time) + data = interpolate_field(field, time) + basis = element(ip, time) + N = length(element) + c = get_connectivity(element) + return sum(data[c[i]]*basis[i] for i=1:N) end -function interpolate_{F<:AbstractField}(element::Element, field::F, ip, time::Float64) - field_ = interpolate(field, time) - basis = element(ip, time) - n = length(element) - if isa(field_, Tuple) - m = length(field_) - 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) - else - c = get_connectivity(element) - return sum(field_[c[i]]*basis[i] for i=1:n) - end +function interpolate_field(::Element, field::CVTV, ip, time) + return field(ip, time) end function size(element::Element, dim) diff --git a/src/elements_lagrange.jl b/src/elements_lagrange.jl index 2b2228a..073fbe2 100644 --- a/src/elements_lagrange.jl +++ b/src/elements_lagrange.jl @@ -8,24 +8,23 @@ using FEMBasis type Poi1 <: AbstractBasis end -function get_basis(element::Element{Poi1}, ip, time) +function get_basis(::Element{Poi1}, ::Any, ::Any) return [1] end -function get_dbasis(element::Element{Poi1}, ip, time) +function get_dbasis(::Element{Poi1}, ::Any, ::Any) return [0] end -function (element::Element{Poi1})(ip, time::Float64, ::Type{Val{:detJ}}) +function (::Element{Poi1})(::Any, ::Float64, ::Type{Val{:detJ}}) return 1.0 end -function get_integration_order(element::Poi1) +function get_integration_order(::Poi1) return 1 end -@lintpragma("Ignore unused order") -function get_integration_points(element::Poi1, order::Int64) +function get_integration_points(::Poi1, ::Int64) return [ (1.0, (0.0, )) ] end @@ -51,7 +50,7 @@ function inside(::Union{Type{Tri3}, Type{Tri6}, Type{Tri7}, Type{Tet4}, Type{Tet return all(xi .>= 0.0) && (sum(xi) <= 1.0) end -function get_reference_coordinates{B}(element::Element{B}) +function get_reference_coordinates{B}(::Element{B}) return get_reference_element_coordinates(B) end diff --git a/src/fields.jl b/src/fields.jl index 35bad06..4ec5a4b 100644 --- a/src/fields.jl +++ b/src/fields.jl @@ -1,6 +1,11 @@ # This file is a part of JuliaFEM. # License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE +""" + AbstractField + +Abstract supertype for all fields in JuliaFEM. +""" abstract type AbstractField end function length{F<:AbstractField}(f::F) @@ -11,109 +16,103 @@ function size{F<:AbstractField}(f::F) return size(f.data) end -function =={F<:AbstractField}(f::F, y) - return ==(f.data, y) +function =={F<:AbstractField}(x::F, y) + return ==(x.data, y) +end + +function =={F<:AbstractField}(x, y::F) + return ==(x, y.data) +end + +function =={F<:AbstractField}(x::F, y::F) + return ==(x.data, y.data) end function getindex{F<:AbstractField}(f::F, i::Int64) return getindex(f.data, i) end -# DISCRETE, CONSTANT, TIME INVARIANT +function interpolate_field(field, ::Any) + return field.data +end + +function update_field!(field, data) + field.data = data +end """ - DCTI(T) + DCTI{T} <: AbstractField Discrete, constant, time-invariant field. + This field is constant in both spatial direction and time direction, i.e. df/dX = 0 and df/dt = 0. + +# Example + +```jldoctest +julia> DCTI(1) +FEMBase.DCTI{Int64}(1) +``` """ type DCTI{T} <: AbstractField data :: T end -""" - update!(f::DCTI, data) - -Update new value to field. -""" -function update!(f::DCTI, data) - f.data = data +function getindex(field::DCTI, ::Int64) + return field.data end -@lintpragma("Ignore unused time") """ - interpolate_(f::DCTI, time) + DVTI{N,T} <: AbstractField -Interpolate constant, time-invariant DCTI field in time direction. That is -trivially only the data itself. -""" -function interpolate_(f::DCTI, time) - return f.data -end - -@lintpragma("Ignore unused i") -""" - getindex(f::DCTI, i) +Discrete, variable, time-invariant field. -Return `i`th item from constant field. Trivially the data itself. -""" -function getindex(f::DCTI, i::Int64) - return f.data -end +This is constant in time direction, but not in spatial direction, i.e. df/dt = 0 +but df/dX != 0. The basic structure of data is `Tuple`, and it is implicitly +assumed that length of field matches to the number of shape functions, so that +interpolation in spatial direction works. -## DISCRETE, VARIABLE, TIME INVARIANT +# Example -""" Discrete, variable, time-invariant field. This is constant in time direction, -but not in spatial direction, i.e. df/dt = 0 but df/dX != 0. The basic structure -of data is Vector, and it is implicitly assumed that length of field matches to -the number of shape functions, so that interpolation in spatial direction works. +```jldoctest +julia> DVTI(1, 2, 3) +FEMBase.DVTI{3,Int64}((1, 2, 3)) +``` """ type DVTI{N,T} <: AbstractField data :: NTuple{N,T} end -""" - update!(f::DVTI, data) - -Update new value to field. -""" -function update!(f::DVTI, data) - f.data = data +function DVTI(data...) + return DVTI(data) end """ - interpolate_(f::DVTI, time) + DCTV{T} <: AbstractField -Interpolate variable, time-invariant DVTI field in time direction. -""" -function interpolate_{N,T}(f::DVTI{N,T}, time) - return f.data -end +Discrete, constant, time variant field. This type of field can change in time +direction but not in spatial direction. -""" - interpolate(a, b) +# Example -A helper function for interpolate routines. Given iterables `a` and `b`, -calculate c = aᵢbᵢ. Length of `a` can be less than `b`, but not vice versa. -""" -function interpolate(a, b) - @assert length(a) <= length(b) - return sum(a[i]*b[i] for i=1:length(a)) -end +Field having value 5 at time 0.0 and value 10 at time 1.0: -## DISCRETE, CONSTANT, TIME VARIANT +```jldoctest +julia> DCTV(0.0 => 5, 1.0 => 10) +FEMBase.DCTV{Int64}(Pair{Float64,Int64}[0.0=>5, 1.0=>10]) +``` +""" type DCTV{T} <: AbstractField data :: Vector{Pair{Float64,T}} end -""" - update!(f::DCTV, time => data) +function DCTV{T}(data::Pair{Float64,T}...) + return DCTV(collect(data)) +end -Update new value to field. -""" -function update!{T}(f::DCTV, data::Pair{Float64, T}) +function update_field!{T}(f::DCTV, data::Pair{Float64, T}) if isapprox(last(f.data).first, data.first) f.data[end] = data else @@ -121,23 +120,7 @@ function update!{T}(f::DCTV, data::Pair{Float64, T}) end end -function DCTV{T}(a::Pair{Float64,T}, b::Pair{Float64,T}) - return DCTV([a, b]) -end - -""" - interpolate_(f::DCTV, time) - -Interpolate constant time variant DCTV field in time direction. - -First algorithm checks that is time out of range, i.e. time is smaller than -time of first frame or larger than last frame. If that is the case, return -first or last frame. Secondly algorithm finds is given time exact match to -time of some frame and return that frame. At last, we find correct bin so -that t0 < time < t1 and use linear interpolation. - -""" -function interpolate_(field::DCTV, time) +function interpolate_field(field::DCTV, time) time < first(field.data).first && return first(field.data).second time > last(field.data).first && return last(field.data).second for i=reverse(1:length(field)) @@ -155,18 +138,28 @@ function interpolate_(field::DCTV, time) end end -## DISCRETE, VARIABLE, TIME VARIANT +""" + DVTV{N,T} <: AbstractField + +Discrete, variable, time variant field. The most general discrete field can +change in both temporal and spatial direction. +# Example + +```jldoctest +julia> DVTV(0.0 => (1, 2), 1.0 => (2, 3)) +FEMBase.DVTV{2,Int64}(Pair{Float64,Tuple{Int64,Int64}}[0.0=>(1, 2), 1.0=>(2, 3)]) +``` +""" type DVTV{N,T} <: AbstractField data :: Vector{Pair{Float64,NTuple{N,T}}} end -""" - update!(f::DVTV, time => data) +function DVTV{N,T}(data::Pair{Float64,NTuple{N,T}}...) + return DVTV(collect(data)) +end -Update new value to field. -""" -function update!{N,T}(f::DVTV, data::Pair{Float64, NTuple{N,T}}) +function update_field!{N,T}(f::DVTV, data::Pair{Float64, NTuple{N,T}}) if isapprox(last(f.data).first, data.first) f.data[end] = data else @@ -174,22 +167,7 @@ function update!{N,T}(f::DVTV, data::Pair{Float64, NTuple{N,T}}) end end -function DVTV{N,T}(a::Pair{Float64,NTuple{N,T}}, b::Pair{Float64,NTuple{N,T}}) - return DVTV([a, b]) -end - -""" - interpolate_(f::DVTV, time) - -Interpolate variable, time variant DVTV field in time direction. - -First algorithm checks that is time out of range, i.e. time is smaller than -time of first frame or larger than last frame. If that is the case, return -first or last frame. Secondly algorithm finds is given time exact match to -time of some frame and return that frame. At last, we find correct bin so -that t0 < time < t1 and use linear interpolation. -""" -function interpolate_{N,T}(field::DVTV{N,T}, time) +function interpolate_field{N,T}(field::DVTV{N,T}, time) time < first(field.data).first && return first(field.data).second time > last(field.data).first && return last(field.data).second for i=reverse(1:length(field)) @@ -208,9 +186,16 @@ function interpolate_{N,T}(field::DVTV{N,T}, time) end """ - CVTV(f) + CVTV <: AbstractField Continuous, variable, time variant field. + +# Example + +```jldoctest +julia> f = CVTV((xi,t) -> xi*t) +FEMBase.CVTV(#1) +``` """ type CVTV <: AbstractField data :: Function @@ -229,21 +214,7 @@ type DVTId{T} <: AbstractField data :: Dict{Int64, T} end -""" - interpolate_(f::DVTId, time) - -Interpolate DVTId, returns trivially the content as this is time invariant field. -""" -function interpolate_(f::DVTId, time) - return f.data -end - -""" - update!(field::DVTId, data::Dict) - -Update data to field. -""" -function update!{T}(field::DVTId{T}, data::Dict{Int64, T}) +function update_field!{T}(field::DVTId{T}, data::Dict{Int64, T}) merge!(field.data, data) end @@ -260,18 +231,7 @@ function DVTVd{T}(data::Pair{Float64,Dict{Int64,T}}...) return DVTVd(collect(data)) end -""" - interpolate_(f::DVTVd, time) - -Interpolate variable, time variant DVTVd dictionary field in time direction. - -# Notes -First check that is outside of range -> extrapolate -Secondly check is "exact match" in time -At last, find the correct bin and use linear interpolation - -""" -function interpolate_{T}(field::DVTVd{T}, time) +function interpolate_field{T}(field::DVTVd{T}, time) time < first(field.data).first && return first(field.data).second time > last(field.data).first && return last(field.data).second for i=reverse(1:length(field)) @@ -293,12 +253,7 @@ function interpolate_{T}(field::DVTVd{T}, time) end end -""" - update!(f::DCTVd, time => data) - -Update new value to dictionary field. -""" -function update!{T}(f::DVTVd, data::Pair{Float64,Dict{Int64,T}}) +function update_field!{T}(f::DVTVd, data::Pair{Float64,Dict{Int64,T}}) if isapprox(last(f.data).first, data.first) f.data[end] = data else @@ -306,148 +261,53 @@ function update!{T}(f::DVTVd, data::Pair{Float64,Dict{Int64,T}}) end end -""" - field(x) - -Create new discrete, constant, time invariant field from value `x`. - -# Example - -```@example abc -f = field(1.0) -``` - -```@example abc -using FEMBase -f = field(1.0) -``` - -```jldoctest -julia> field(1.0) -FEMBase.DCTI{Float64}(1.0) -``` - -""" -function field(x) - return DCTI(x) +function new_field(data) + return DCTI(data) end -""" - field(x::NTuple{N,T}) - -Create new discrete, variable, time invariant field from tuple `x`. - -# Example - -```julia -f = field( (1.0, 2.0) ) -``` -""" -function field{N,T}(data::NTuple{N,T}) +function new_field(data...) return DVTI(data) end -""" - field(x::Pair{Float64,T}) - -Create new discrete, constant, time variant field from pair `x`. - -# Example +function new_field{N,T}(data::NTuple{N,T}) + return DVTI(data) +end -```julia -f = field(1.0=>1.0, 2.0=>2.0) -``` -""" -function field{T}(data::Pair{Float64,T}...) +function new_field{T}(data::Pair{Float64,T}...) return DCTV(collect(data)) end -""" - field(x::Pair{Float64,NTuple{N,T}) - -Create new discrete, variable, time variant field from pair `x`. - -# Example - -```julia -f = field(1.0=>(1.0,2.0), 2.0=>(2.0,3.0)) -``` -""" -function field{N,T}(data::Pair{Float64,NTuple{N,T}}...) +function new_field{N,T}(data::Pair{Float64,NTuple{N,T}}...) return DVTV(collect(data)) end -""" - field(x::Function) - -Create new, continuous, variable, time variant field from function `x`. - -# Example - -```julia -f = field( (xi,t) -> xi[1]*t ) -``` -""" -function field(data::Function) +function new_field(data::Function) return CVTV(data) end -""" - field(x::Pair{Int64,T}) - -Create new discrete, variable, time invariant dictionary field from `x`. - -# Example - -```julia -f = field(1 => 1.0, 2 => 2.0) -``` -""" -function field{T}(data::Pair{Int64, T}...) +function new_field{T}(data::Pair{Int64, T}...) return DVTId(Dict(data)) end -""" - field(x::Pair{Float64, NTuple{N, Pair{Int64, T}}}) - -Create new discrete, variable, time variant dictionary field from `x`. - -# Example - -```julia -X1 = (1 => 1.0, 2 => 2.0) -X2 = (1 => 2.0, 2 => 3.0) -f = field(0.0 => X1, 1.0 => X2) -``` -""" -function field{N,T}(data::Pair{Float64, NTuple{N, Pair{Int64, T}}}...) +function new_field{N,T}(data::Pair{Float64, NTuple{N, Pair{Int64, T}}}...) return DVTVd(collect(t => Dict(d) for (t, d) in data)) end -""" - field(x::Dict) - -Create new discrete, variable, time invariant dictionary field from dictionary `x`. - -# Example - -```julia -X = Dict(1 => [0.0,0.0], 2 => [1.0,0.0]) -f = field(X) -``` -""" -function field{T}(data::Dict{Int64,T}) +function new_field{T}(data::Dict{Int64,T}) return DVTId(data) end +function new_field{T}(data::Pair{Float64, Dict{Int64, T}}...) + return DVTVd(collect(data)) +end + """ - field(t::Float64 => x::Dict, ...) + field(x) -Create new discrete, variable, time variant dictionary field from pair of time -and dictionary. +Create new field. Field type is deduced from data type. """ -function field{T}(data::Pair{Float64, Dict{Int64, T}}...) - return DVTVd(collect(data)) +function field(data...) + return new_field(data...) end """ @@ -457,8 +317,9 @@ Interpolate field in time direction. # Examples -For time invariant fields `DCTI`, `DVTI`, `DVTId` solution is trivially the -data inside field as fields does not depend from time: +For time invariant fields [`DCTI`](@ref), [`DVTI`](@ref), [`DVTId`](@ref) +solution is trivially the data inside field as fields does not depend from +the time: ```jldoctest julia> a = field(1.0) @@ -486,11 +347,38 @@ Dict{Int64,Float64} with 2 entries: 1 => 1.0 ``` -For time invariant fields DCTI, DVTI, DVTId trivial solution is returned. For time variant fields DCTV, DVTV, DVTVd linear interpolation is performed. +# Other notes + +First algorithm checks that is time out of range, i.e. time is smaller than +time of first frame or larger than last frame. If that is the case, return +first or last frame. Secondly algorithm finds is given time exact match to +time of some frame and return that frame. At last, we find correct bin so +that t0 < time < t1 and use linear interpolation. + """ function interpolate{F<:AbstractField}(field::F, time) - return interpolate_(field, time) + return interpolate_field(field, time) +end + +""" + interpolate(a, b) + +A helper function for interpolate routines. Given iterables `a` and `b`, +calculate c = aᵢbᵢ. Length of `a` can be less than `b`, but not vice versa. +""" +function interpolate(a, b) + @assert length(a) <= length(b) + return sum(a[i]*b[i] for i=1:length(a)) +end + +""" + update!(field, data) + +Update new value to field. +""" +function update!{F<:AbstractField}(field::F, data) + update_field!(field, data) end diff --git a/src/integrate.jl b/src/integrate.jl index e727116..44e9990 100644 --- a/src/integrate.jl +++ b/src/integrate.jl @@ -51,6 +51,6 @@ for (E, R) in integration_rule_mapping end # All good codes needs a special case. Here we have it: Poi1 -function get_integration_points(element::Poi1) +function get_integration_points(::Poi1) [ (1.0, (0.0,) ) ] end diff --git a/src/problems.jl b/src/problems.jl index b4fc164..37e669e 100644 --- a/src/problems.jl +++ b/src/problems.jl @@ -223,7 +223,9 @@ function initialize!(problem::Problem, element::Element, time::Float64) if field_dim == 1 # scalar field empty_field = tuple(zeros(nnodes)...) else # vector field - empty_field = tuple([zeros(field_dim) for i=1:nnodes]...) + # FIXME: the most effective way to do + # ([0.0,0.0], [0.0,0.0], ..., [0.0,0.0]) ? + empty_field = tuple(map((x)->zeros(field_dim)*x, 1:nnodes)...) end # initialize primary field diff --git a/src/sparse.jl b/src/sparse.jl index d3374af..e6002ff 100644 --- a/src/sparse.jl +++ b/src/sparse.jl @@ -57,12 +57,25 @@ function sparse(A::SparseMatrixCOO, n::Int, m::Int, f::Function; tol=1.0e-12) return B end +# For backward compatibility, will be removed function push!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64) push!(A.I, I) push!(A.J, J) push!(A.V, V) end +function add!(A::SparseMatrixCOO, I::Int, J::Int, V::Float64) + push!(A.I, I) + push!(A.J, J) + push!(A.V, V) +end + +function add!(A::SparseMatrixCOO, I::Int, V::Float64) + push!(A.I, I) + push!(A.J, 1) + push!(A.V, V) +end + function empty!(A::SparseMatrixCOO) empty!(A.I) empty!(A.J) @@ -107,11 +120,10 @@ function add!(A::SparseMatrixCOO, dofs1::Vector{Int}, dofs2::Vector{Int}, data:: n, m = size(data) for j=1:m for i=1:n - push!(A.I, dofs1[i]) - push!(A.J, dofs2[j]) + add!(A, dofs1[i], dofs2[j], data[i,j]) end end - append!(A.V, vec(data)) + return nothing end """ Add sparse matrix of CSC to COO. """ diff --git a/test/test_elements.jl b/test/test_elements.jl index 0a27a89..c668ac9 100644 --- a/test/test_elements.jl +++ b/test/test_elements.jl @@ -25,7 +25,7 @@ end el = Element(Quad4, [1, 2, 3, 4]) pr = Problem(Dummy, "problem", 2) add_elements!(pr, [el]) - assemble!(pr) + assemble!(pr, 0.0) @test true end @@ -343,12 +343,6 @@ end @test isapprox(element("test", (0.0, 0.0), 0.5), 0.5) end -@testset "interpolation failure" begin - el = Element(Seg2, [1, 2]) - update!(el, "g", (1.0, 2.0, 3.0)) - @test_throws Exception el("g", [0.0], 0.0) -end - @testset "get integration points for Quad4" begin el = Element(Quad4, [1, 2, 3, 4]) ips = get_integration_points(el) diff --git a/test/test_fields.jl b/test/test_fields.jl index fac294d..94ca3a8 100644 --- a/test/test_fields.jl +++ b/test/test_fields.jl @@ -44,6 +44,7 @@ end @test interpolate(a, 0.0) == (1, 2) update!(a, (2,3)) @test a == (2,3) + @test (2,3) == a # vector field b = DVTI(([1,2], [2,3])) @@ -60,6 +61,9 @@ end @test interpolate(c, 0.0) == ([1 2; 3 4], [2 3; 4 5]) update!(c, ([2 3; 4 5], [5 6; 7 8])) @test c == ([2 3; 4 5], [5 6; 7 8]) + + d = DVTI(2, 3) + @test a == d end @testset "DCTV field" begin @@ -124,6 +128,7 @@ end @test isa(field(1.0), DCTI) @test isa(field(1.0 => 1.0), DCTV) @test isa(field((1.0,2.0)), DVTI) + @test isa(field(1, 2), DVTI) @test isa(field(1.0 => (1.0,2.0)), DVTV) @test isa(field((xi,t) -> xi[1]*t), CVTV) @test isa(field(1 => [1.0, 2.0], 10 => [2.0, 3.0]), DVTId) diff --git a/test/test_problems.jl b/test/test_problems.jl index fabcb0c..a8be814 100644 --- a/test/test_problems.jl +++ b/test/test_problems.jl @@ -8,7 +8,7 @@ using FEMBase: get_assembly, get_elements using Base.Test import FEMBase: get_unknown_field_dimension, get_unknown_field_name -import FEMBase: get_formulation_type, assemble! +import FEMBase: get_formulation_type, assemble_elements! type P1 <: FieldProblem A :: Bool @@ -26,9 +26,17 @@ function P2() return P2(:incremental) end +type P3 <: FieldProblem +end + get_unknown_field_name(p::P1) = "P1" get_formulation_type(p::Problem{P2}) = p.properties.formulation +@testset "test get_unknown_field_name" begin + p = Problem(P3, "P3", 1) + @test get_unknown_field_name(p) == "N/A" +end + @testset "test creating new problems" begin p1 = Problem(P1, "P1", 1) p2 = Problem(P2, "P2", 2, "p") @@ -130,10 +138,10 @@ end @test isapprox(e2("lambda", (0.0,), 0.0), [0.0, 0.0]) end -function assemble!{E}(problem::Problem{P1}, - assembly::Assembly, - elements::Vector{Element{E}}, - time::Float64) +function assemble_elements!{E}(problem::Problem{P1}, + assembly::Assembly, + elements::Vector{Element{E}}, + time::Float64) info("Assembling elements of kind $E") bi = BasisInfo(E) @@ -156,7 +164,7 @@ function assemble!{E}(problem::Problem{P1}, end -@testset "test assemble test problem" begin +@testset "test assemble field problem" begin el1 = Element(Quad4, [1, 2, 3, 4]) el2 = Element(Tri3, [3, 2, 5]) X = Dict(1 => [0.0, 0.0], @@ -181,3 +189,72 @@ end 0.0 0.0 -3.0 0.0 3.0] @test isapprox(K, K_expected) end + +type DirBC <: BoundaryProblem +end + +function assemble_elements!{E}(problem::Problem{DirBC}, + assembly::Assembly, + elements::Vector{Element{E}}, + time::Float64) + + name = get_parent_field_name(problem) + dim = get_unknown_field_dimension(problem) + + data = Dict{Int64,Float64}() + for element in elements + for i=1:dim + haskey(element, "$name $dim") || continue + gdofs = get_gdofs(problem, element) + ldofs = gdofs[i:dim:end] + xis = get_reference_element_coordinates(E) + for (ldof, xi) in zip(ldofs, xis) + data[ldof] = interpolate(element, "$name $dim", xi, time) + end + end + end + + for (dof, val) in data + add!(assembly.C1, dof, dof, 1.0) + add!(assembly.C2, dof, dof, 1.0) + add!(assembly.g, dof, val) + end + + return nothing + +end + +@testset "test assemble boundary problem" begin + el1 = Element(Seg2, [1, 2]) + el2 = Element(Seg2, [2, 3]) + X = Dict(1 => [0.0, 0.0], + 2 => [1.0, 0.0], + 3 => [1.0, 1.0], + 4 => [0.0, 1.0], + 5 => [2.0, 1.0]) + elements = [el1, el2] + update!(elements, "geometry", X) + T0 = Dict(1 => 0.5, 2 => 0.0, 3 => 1.0) + update!(elements, "temperature 1", T0) + problem = Problem(DirBC, "boundary", 1, "temperature") + add_elements!(problem, elements) + time = 0.0 + assemble!(problem, time) + C = full(problem.assembly.C1) + g = full(problem.assembly.g) + println(C) + println(g) + C_expected = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] + g_expected = [0.5, 0.0, 1.0] + @test isapprox(C, C_expected) + @test isapprox(g, g_expected) +end + +@testset "test deprecated assembly procedure" begin + problem = Problem(P3, "P3", 1) + elements = [Element(Seg2, [1, 2])] + add_elements!(problem, elements) + assemble!(problem.assembly, problem, first(elements), 0.0) + assemble!(problem.assembly, problem, problem.elements, 0.0) + assemble_elements!(problem, problem.assembly, elements, 0.0) +end