mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
data structures iteration #3
This commit is contained in:
Jukka Aho
File diff suppressed because it is too large
Load Diff
@@ -89,13 +89,13 @@
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"execution_count": 69,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM: Element, Field, FieldSet, Basis, Quad4"
|
||||
"using JuliaFEM: Element, Field, FieldSet, Basis"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1724,6 +1724,109 @@
|
||||
"source": [
|
||||
"In this notebook the basic instructions how to develop JuliaFEM has been given. The most imporant concepts has been considered; how to develop own element with own basis, several ways how to define own equation, and how to finally assemble and calculate the problem using solver. Any comments and/or discussion about technical details, theory, programming, or from life in general is very desirable; our issue log is in address https://github.com/JuliaFEM/JuliaFEM.jl/issues"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Advanced stuff\n",
|
||||
"\n",
|
||||
"In last section of this tutorial we consider some of the more advanced things which may araise when developing own models.\n",
|
||||
"\n",
|
||||
"### Boundary element access to parent element + overriding equations in problems\n",
|
||||
"\n",
|
||||
"This kind of situation might happen when one is almost happy for some problem setting, but would like to change just one or two equations from it. For example boundary equation is not satisfying all the requirements and one would like to test something new. \n",
|
||||
"\n",
|
||||
"### Accessing integration points\n",
|
||||
"\n",
|
||||
"### Fields as a function of something.\n",
|
||||
"- statistical variables\n",
|
||||
"- field dependent from another field\n",
|
||||
"- field dependent from time or spatial domain etc.\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"using JuliaFEM: get_default_integration_points\n",
|
||||
"\"\"\" 2-node radiation boundary element. \"\"\"\n",
|
||||
"type DC2D2RAD <: Heat\n",
|
||||
" element :: Seg2\n",
|
||||
" integration_points :: Array{IntegrationPoint, 1}\n",
|
||||
"end\n",
|
||||
"function DC2D2RAD(element::Seg2)\n",
|
||||
" integration_points = [\n",
|
||||
" IntegrationPoint([0.0], 2.0)]\n",
|
||||
" if !haskey(element, \"temperature\")\n",
|
||||
" element[\"temperature\"] = FieldSet()\n",
|
||||
" push!(element[\"temperature\"], Field([0.0, 0.0]))\n",
|
||||
" end\n",
|
||||
" DC2D2RAD(element, integration_points)\n",
|
||||
"end\n",
|
||||
"Base.size(equation::DC2D2RAD) = (1, 2)\n",
|
||||
"\n",
|
||||
"\"\"\" Calculate potential energy caused by radiation.\n",
|
||||
"https://en.wikipedia.org/wiki/Stefan%E2%80%93Boltzmann_constant\n",
|
||||
"\"\"\"\n",
|
||||
"function JuliaFEM.get_potential_energy(equation::DC2D2RAD, ip, time; variation=nothing)\n",
|
||||
" element = get_element(equation)\n",
|
||||
" basis = get_basis(element)\n",
|
||||
" eps = basis(\"emissivity\", ip, time)\n",
|
||||
" #sig = basis(\"stefan-boltzmann constant\", ip, time)\n",
|
||||
" sig = 5.670367e-8 # i guess stefan-boltzmann constant is constant ;)\n",
|
||||
" T = basis(\"temperature\", ip, time, variation)\n",
|
||||
" T_ext = basis(\"temperature external\", ip, time)\n",
|
||||
" q = eps*sig*((T_ext+273.15)^4 - (T+273.15)^4)\n",
|
||||
" println(ForwardDiff.value(q*T))\n",
|
||||
" return q\n",
|
||||
"end\n",
|
||||
"JuliaFEM.has_potential_energy(equation::DC2D2RAD) = true\n",
|
||||
"\n",
|
||||
"Defining new equation mapping to old problem is one line command. Here we replace `DC2D2` $\\rightarrow$ `DC2DCRAD`\n",
|
||||
"\n",
|
||||
"function run_radiation_model()\n",
|
||||
" # this is the same as before\n",
|
||||
" element = Quad4([1, 2, 3, 4])\n",
|
||||
" fieldset1 = FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])])\n",
|
||||
" fieldset2 = FieldSet(\"temperature thermal conductivity\", [Field(6.0)])\n",
|
||||
" fieldset3 = FieldSet(\"temperature load\", [Field([12.0, 12.0, 12.0, 12.0])])\n",
|
||||
" fieldset4 = FieldSet(\"density\", [Field(36.0)])\n",
|
||||
" push!(element, fieldset1)\n",
|
||||
" push!(element, fieldset2)\n",
|
||||
" push!(element, fieldset3)\n",
|
||||
" push!(element, fieldset4)\n",
|
||||
"\n",
|
||||
" # create boundary element\n",
|
||||
" boundary_element = Seg2([1, 2])\n",
|
||||
" push!(boundary_element, FieldSet(\"geometry\", [Field(Vector[[0.0, 0.0], [1.0, 0.0]])]))\n",
|
||||
" push!(boundary_element, FieldSet(\"emissivity\", [Field(0.5)]))\n",
|
||||
" push!(boundary_element, FieldSet(\"temperature external\", [Field(20.0)]))\n",
|
||||
"\n",
|
||||
" # set initial conditions\n",
|
||||
" push!(element, FieldSet(\"temperature\", [Field([0.0, 0.0, 0.0, 0.0])]))\n",
|
||||
" push!(boundary_element, FieldSet(\"temperature\", [Field([0.0, 1.0])]))\n",
|
||||
" \n",
|
||||
" # create problem, change element mapping\n",
|
||||
" problem = PlaneHeatProblem()\n",
|
||||
" problem[Seg2] = DC2D2RAD # Seg2 was previous DC2D2\n",
|
||||
" push!(problem, element)\n",
|
||||
" push!(problem, boundary_element)\n",
|
||||
"\n",
|
||||
" # run our \"unit test solver\"\n",
|
||||
" free_dofs = [3, 4]\n",
|
||||
" solve!(problem, free_dofs; max_iterations=4, dump_matrices=true)\n",
|
||||
" basis = get_basis(boundary_element)\n",
|
||||
" T = basis(\"temperature\", [0.0])\n",
|
||||
" println(\"Temperature at the midpoint of element: $T\")\n",
|
||||
"end\n",
|
||||
"\n",
|
||||
"run_radiation_model()"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
File diff suppressed because one or more lines are too long
+9
-5
@@ -6,10 +6,14 @@ This is JuliaFEM -- Finite Element Package
|
||||
"""
|
||||
module JuliaFEM
|
||||
|
||||
using Logging
|
||||
@Logging.configure(level=DEBUG)
|
||||
#using Logging
|
||||
#@Logging.configure(level=DEBUG)
|
||||
#using Lexicon
|
||||
|
||||
macro debug(msg)
|
||||
return :( println("DEBUG: ", $msg) )
|
||||
end
|
||||
|
||||
using Lexicon
|
||||
using ForwardDiff
|
||||
autodiffcache = ForwardDiffCache()
|
||||
|
||||
@@ -24,12 +28,12 @@ Examples
|
||||
[1.0]
|
||||
|
||||
"""
|
||||
function Base.linspace(X1, X2, n)
|
||||
function Base.linspace{T<:Array}(X1::T, X2::T, n)
|
||||
[1/2*(1-ti)*X1 + 1/2*(1+ti)*X2 for ti in linspace(-1, 1, n)]
|
||||
end
|
||||
|
||||
include("types.jl") # type definitions
|
||||
include("interpolate.jl") # interpolation routines
|
||||
#include("interpolate.jl") # interpolation routines
|
||||
|
||||
### ELEMENTS ###
|
||||
include("elements.jl")
|
||||
|
||||
+2
-2
@@ -49,8 +49,8 @@ function calculate_global_assembly!(assembly::GlobalAssembly, problem::Problem,
|
||||
unknown_field_name = get_unknown_field_name(problem)
|
||||
initialize_global_assembly!(assembly, problem) # zero all
|
||||
dim, ndofs = size(problem)
|
||||
Logging.info("assembling problem for $unknown_field_name")
|
||||
Logging.info("dimension of unknown field: $dim, problem dofs: $ndofs")
|
||||
info("assembling problem for $unknown_field_name")
|
||||
info("dimension of unknown field: $dim, problem dofs: $ndofs")
|
||||
local_assembly = initialize_local_assembly()
|
||||
for (i, equation) in enumerate(get_equations(problem))
|
||||
calculate_local_assembly!(local_assembly, equation, unknown_field_name, time)
|
||||
|
||||
+1
-2
@@ -121,8 +121,7 @@ end
|
||||
function CPS4(element::Quad4)
|
||||
integration_points = get_default_integration_points(element)
|
||||
if !haskey(element, "displacement")
|
||||
element["displacement"] = FieldSet()
|
||||
push!(element["displacement"], Field(Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]))
|
||||
element["displacement"] = zeros(2, 4)
|
||||
end
|
||||
CPS4(element, integration_points)
|
||||
end
|
||||
|
||||
+7
-47
@@ -69,56 +69,13 @@ end
|
||||
|
||||
"""Add new FieldSet to element.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> field = Field(0.0, [1, 2, 3, 4])
|
||||
>>> fieldset = FieldSet("geometry", Field[field])
|
||||
>>> element["geometry"] = fieldset
|
||||
JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
|
||||
"""
|
||||
function Base.setindex!(element::Element, fieldset::FieldSet, fieldset_name)
|
||||
fieldset.name = fieldset_name
|
||||
element.fields[fieldset.name] = fieldset
|
||||
end
|
||||
|
||||
"""Add new FieldSet to element.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> field = Field(0.0, [1, 2, 3, 4])
|
||||
>>> element["geometry"] = field
|
||||
JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
|
||||
"""
|
||||
function Base.setindex!(element::Element, field::Field, fieldset_name)
|
||||
element[fieldset_name] = FieldSet(field)
|
||||
end
|
||||
|
||||
"""Add new FieldSet to element.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> element["geometry"] = [1, 2, 3, 4]
|
||||
JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("geometry"=>JuliaFEM.FieldSet("geometry",JuliaFEM.Field[JuliaFEM.Field{Array{Int64,1}}(0.0,0,[1,2,3,4])])))
|
||||
"""
|
||||
function Base.setindex!(element::Element, field_data::Union{Number, Array}, fieldset_name)
|
||||
element[fieldset_name] = Field(field_data)
|
||||
end
|
||||
|
||||
"""Add new FieldSet to element.
|
||||
|
||||
Notes
|
||||
-----
|
||||
This last version takes tuple and each cell in tuple is converted to new field.
|
||||
Time in field is 0.0, 1.0, ..., n
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> element["load"] = (1, 2)
|
||||
JuliaFEM.Quad4([1,2,3,4],JuliaFEM.Basis(basis,dbasisdxi),Dict("load"=>JuliaFEM.FieldSet("load",JuliaFEM.Field[JuliaFEM.Field{Int64}(0.0,0,1),JuliaFEM.Field{Int64}(1.0,0,2)])))
|
||||
"""
|
||||
function Base.setindex!(element::Element, field_data::Tuple, fieldset_name)
|
||||
fields = Field[Field(Float64(i-1), field) for (i,field) in enumerate(field_data)]
|
||||
element.fields[fieldset_name] = FieldSet(fieldset_name, fields)
|
||||
function Base.setindex!(element::Element, field_data, field_name)
|
||||
element.fields[field_name] = field_data
|
||||
end
|
||||
|
||||
function get_connectivity(el::Element)
|
||||
@@ -156,10 +113,13 @@ end
|
||||
function call(u::FunctionSpace, field_name, xi::Vector, t::Number=Inf, variation=nothing)
|
||||
f = !isa(variation, Void) ? variation : u.element[field_name](t)
|
||||
if length(f) == 1
|
||||
return f.values
|
||||
return f.data[1]
|
||||
end
|
||||
h = u.element.basis.basis(xi)
|
||||
return dot(vec(h), f)
|
||||
#@debug("vec(h) = $(vec(h)), size(h) = $(size(vec(h)))")
|
||||
#@debug("f = $f, size(f) = $(size(f))")
|
||||
#return dot(vec(h), f)
|
||||
return sum(vec(h).*f)
|
||||
end
|
||||
|
||||
""" If basis is called without a field, return basis functions evaluated at that point. """
|
||||
|
||||
+4
-1
@@ -182,7 +182,10 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::Equation,
|
||||
field = element[unknown_field_name](time)
|
||||
function residual_vector(data::Vector)
|
||||
fill!(assembly.residual_vector, 0.0)
|
||||
df = similar(field, data)
|
||||
#@debug("field: $field, length = $(size(field))")
|
||||
#@debug("data: $data, size = $(size(data))")
|
||||
#df = similar(field, data)
|
||||
df = Increment(reshape(data, size(equation)...))
|
||||
# integrate W
|
||||
for ip in get_integration_points(equation)
|
||||
dr = get_residual_vector(equation, ip, time; variation=df)
|
||||
|
||||
+4
-4
@@ -46,8 +46,8 @@ function calculate_local_assembly!(assembly::LocalAssembly, equation::HeatEquati
|
||||
w = ip.weight * detJ(ip)
|
||||
N = basis(ip, time)
|
||||
if haskey(element, "density")
|
||||
ρ = basis("density", ip, time)
|
||||
assembly.mass_matrix += w * ρ*N'*N
|
||||
rho = basis("density", ip, time)
|
||||
assembly.mass_matrix += w * rho*N'*N
|
||||
end
|
||||
if haskey(element, "temperature thermal conductivity")
|
||||
dN = dbasis(ip, time)
|
||||
@@ -92,7 +92,7 @@ end
|
||||
function DC2D4(element::Quad4)
|
||||
integration_points = get_default_integration_points(element)
|
||||
if !haskey(element, "temperature")
|
||||
element["temperature"] = FieldSet()
|
||||
element["temperature"] = zeros(4)
|
||||
end
|
||||
DC2D4(element, integration_points)
|
||||
end
|
||||
@@ -106,7 +106,7 @@ end
|
||||
function DC2D2(element::Seg2)
|
||||
integration_points = get_default_integration_points(element)
|
||||
if !haskey(element, "temperature")
|
||||
element["temperature"] = FieldSet()
|
||||
element["temperature"] = zeros(2)
|
||||
end
|
||||
DC2D2(element, integration_points)
|
||||
end
|
||||
|
||||
+4
-3
@@ -37,13 +37,14 @@ macro create_lagrange_element(element_name, element_description, X, P)
|
||||
type $eltype <: CG
|
||||
connectivity :: Array{Int, 1}
|
||||
basis :: Basis
|
||||
fields :: Dict{ASCIIString, FieldSet}
|
||||
fields :: FieldSet
|
||||
end
|
||||
function $eltype(connectivity, args...)
|
||||
$eltype(connectivity, Basis(basis, dbasisdxi), Dict())
|
||||
$eltype(connectivity, Basis(basis, dbasisdxi), FieldSet())
|
||||
end
|
||||
get_element_description(el::Type{$eltype}) = $element_description
|
||||
Base.size(el::Type{$eltype}) = Base.size($X)
|
||||
Base.size(element::Type{$eltype}) = Base.size($X)
|
||||
Base.size(element::$eltype) = Base.size($X)
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
+1
-3
@@ -67,9 +67,7 @@ function solve!(problem::Problem, free_dofs::Array{Int, 1}, time::Number=Inf;
|
||||
gdofs = vec(vcat([dim*conn'-i for i=dim-1:-1:0]...))
|
||||
old_field = element[field_name](Inf)
|
||||
new_field = similar(old_field, full(x[gdofs]))
|
||||
new_field.time = time
|
||||
new_field.increment = i
|
||||
push!(element[field_name], new_field)
|
||||
push!(element[field_name][end], new_field)
|
||||
end
|
||||
if norm(dx) < tolerance
|
||||
return
|
||||
|
||||
+233
-125
@@ -3,150 +3,272 @@
|
||||
|
||||
# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
|
||||
using ForwardDiff
|
||||
#abstract AbstractField{T,N} <: AbstractArray{T,N}
|
||||
|
||||
""" Field is a fundamental data type which holds some values in some time t """
|
||||
type Field{T}
|
||||
time :: Number
|
||||
increment :: Int64
|
||||
values :: T
|
||||
abstract AbstractField
|
||||
|
||||
abstract DiscreteField <: AbstractField
|
||||
|
||||
abstract ContinuousField <: AbstractField
|
||||
abstract TimeContinuousField <: ContinuousField
|
||||
abstract SpatialContinuousField <: ContinuousField
|
||||
abstract TimeAndSpatialContinuousField <: ContinuousField
|
||||
# should we introduce time and spatial discontinuous fields
|
||||
# for discontinuous galerkin?
|
||||
|
||||
### DEFAULT DISCRETE FIELD ###
|
||||
|
||||
# 1. Increment
|
||||
|
||||
# FIXME: This should be Vector.
|
||||
#typealias Increment Vector
|
||||
type Increment{T} <: AbstractVector{T}
|
||||
data :: Vector{T}
|
||||
end
|
||||
""" Initialize field. """
|
||||
function Field(time, values)
|
||||
Field(time, 0, values)
|
||||
Base.size(increment::Increment) = Base.size(increment.data)
|
||||
Base.linearindexing(::Type{Increment}) = Base.LinearFast()
|
||||
Base.getindex(increment::Increment, i::Int) = increment.data[i]
|
||||
Base.setindex!(increment::Increment, v, i::Int) = (increment.data[i] = v)
|
||||
Base.similar{T}(increment::Increment, ::Type{T}) = Increment(similar(increment.data))
|
||||
Base.dot(v::Number, i::Increment) = v*i
|
||||
|
||||
function Base.convert(::Type{Increment}, data::Number)
|
||||
Increment([data])
|
||||
end
|
||||
function Field(values)
|
||||
Field(0.0, 0, values)
|
||||
function Base.convert{T}(::Type{Increment}, data::Array{T, 2})
|
||||
Increment([data[:,i] for i=1:size(data, 2)])
|
||||
end
|
||||
""" Get length of a field (number of basis functions in practice). """
|
||||
function Base.length(f::Field)
|
||||
length(f.values)
|
||||
function Base.convert{T}(::Type{Increment}, data::Array{T, 3})
|
||||
Increment([data[:,:,i] for i=1:size(data, 3)])
|
||||
end
|
||||
""" Push value to field. """
|
||||
function Base.push!(f::Field, value)
|
||||
push!(f.values, value)
|
||||
function Base.convert{T}(::Type{Increment}, data::Array{T, 4})
|
||||
Increment([data[:,:,:,i] for i=1:size(data, 4)])
|
||||
end
|
||||
""" Get field discrete value at point i. """
|
||||
function Base.getindex(f::Field, i::Int64)
|
||||
f.values[i]
|
||||
function Base.convert{T}(::Type{Increment}, data::Array{T, 5})
|
||||
Increment([data[:,:,:,:,i] for i=1:size(data, 5)])
|
||||
end
|
||||
""" Multiply field with some constant k. """
|
||||
function Base.(:*)(k::Number, f::Field)
|
||||
Field(f.time, k*f.values)
|
||||
function Base.zeros(::Type{Increment}, dims...)
|
||||
Increment(zeros(dims...))
|
||||
end
|
||||
function Base.vec(increment::Increment)
|
||||
[increment.data...;]
|
||||
end
|
||||
function Base.similar{T}(increment::Increment{Vector{T}}, data::Vector{T})
|
||||
Increment(reshape(data, round(Int, length(data)/length(increment)), length(increment)))
|
||||
end
|
||||
|
||||
""" Inner product of field and vector x. """
|
||||
function Base.dot(x::Vector, f::Field)
|
||||
@assert length(x) == length(f)
|
||||
sum([f[i]*x[i] for i in 1:length(f)])
|
||||
# 2. TimeStep
|
||||
|
||||
type TimeStep{T} <: AbstractVector{T}
|
||||
time :: Float64
|
||||
increments :: Vector{T}
|
||||
end
|
||||
Base.size(timestep::TimeStep) = Base.size(timestep.increments)
|
||||
Base.linearindexing(::Type{TimeStep}) = Base.LinearFast()
|
||||
Base.getindex(timestep::TimeStep, i::Int) = timestep.increments[i]
|
||||
|
||||
function Base.convert(::Type{TimeStep}, time::Number, increment::Increment)
|
||||
TimeStep(time, Increment[increment])
|
||||
end
|
||||
|
||||
function Base.size(field::Field)
|
||||
(length(field.values[1]), length(field.values))
|
||||
function Base.push!(timestep::TimeStep, increment::Increment)
|
||||
push!(timestep.increments, increment)
|
||||
end
|
||||
|
||||
#""" Multiply field with some matrix x. """
|
||||
# function Base.(:*){T}(x::Matrix, f::Field{Vector{T}})
|
||||
#function Base.(:*)(x::Matrix, f::Field)
|
||||
# sum([f[i]*x[:,i]' for i in 1:length(f)])
|
||||
# 3. DefaultDiscreteField
|
||||
|
||||
type DefaultDiscreteField <: DiscreteField
|
||||
timesteps :: Vector{TimeStep}
|
||||
end
|
||||
Base.size(field::DefaultDiscreteField) = Base.size(field.timesteps)
|
||||
Base.linearindexing(::Type{DefaultDiscreteField}) = Base.LinearFast()
|
||||
Base.getindex(field::DefaultDiscreteField, i::Int) = field.timesteps[i]
|
||||
Base.length(field::DefaultDiscreteField) = length(field.timesteps)
|
||||
Base.endof(field::DefaultDiscreteField) = endof(field.timesteps)
|
||||
Base.first(field::DefaultDiscreteField) = field[1][end]
|
||||
Base.last(field::DefaultDiscreteField) = field[end][end]
|
||||
function Base.push!(field::DefaultDiscreteField, timestep::TimeStep)
|
||||
push!(field.timesteps, timestep)
|
||||
end
|
||||
|
||||
|
||||
typealias Field DefaultDiscreteField
|
||||
|
||||
### CONTINUOUS FIELDS ###
|
||||
|
||||
|
||||
# fix print_matrix
|
||||
#function Base.print_matrix(::Base.AbstractIOBuffer, field::ContinuousField, args...)
|
||||
# TODO: anything nice to print?
|
||||
#end
|
||||
|
||||
""" Sum two fields. """
|
||||
function Base.(:+)(f1::Field, f2::Field)
|
||||
@assert(f1.time == f2.time, "Cannot add fields: time mismatch, $(f1.time) != $(f2.time)")
|
||||
Field(f1.time, f1.values + f2.values)
|
||||
end
|
||||
### FIELDSET ###
|
||||
|
||||
""" Return data from field as a long array.
|
||||
typealias FieldSet Dict{ASCIIString, AbstractField}
|
||||
|
||||
"""Quicky add discrete field to fieldset.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> f = Field(0.0, Vector[[1.0, 2.0], [3.0, 4.0]])
|
||||
>>> f[:]
|
||||
[1.0, 2.0, 3.0, 4.0]
|
||||
|
||||
>>> fs = FieldSet()
|
||||
>>> fs["myfield"] = [1, 2, 3, 4]
|
||||
"""
|
||||
function Base.getindex(field::Field, c::Colon)
|
||||
[field.values...;]
|
||||
end
|
||||
function Base.vec(field::Field)
|
||||
[field.values...;]
|
||||
function Base.convert(::Type{AbstractField}, data::Union{Array, Number})
|
||||
increment = Increment(data)
|
||||
timestep = TimeStep(0.0, Increment[increment])
|
||||
field = DefaultDiscreteField(TimeStep[timestep])
|
||||
return field
|
||||
end
|
||||
|
||||
""" Return field similar to input but with new data in it.
|
||||
""" Quicky add several time steps at once in tuple.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> f = Field(0.5, Vector[[1.0, 2.0], [3.0, 4.0]])
|
||||
>>> similar(f, ones(4))
|
||||
JuliaFEM.Field{Array{Array{T,1},1}}(0.5,1,Array{T,1}[[1.0,1.0],[1.0,1.0]])
|
||||
>>> fs = FieldSet()
|
||||
>>> fs["myfield"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5])
|
||||
|
||||
or
|
||||
|
||||
>>> fs["myfield"] = [1, 2, 3, 4], [2, 3, 4, 5]
|
||||
|
||||
"""
|
||||
function Base.similar(field::Field, data::Vector)
|
||||
fdim = round(Int, length(data)/length(field)) # dimension of field variable
|
||||
if fdim == 1
|
||||
new_field = Field(field.time, data)
|
||||
return new_field
|
||||
function Base.convert(::Type{AbstractField}, data::Tuple)
|
||||
timesteps = TimeStep[]
|
||||
for (i, timestep) in enumerate(data)
|
||||
if isa(timestep, Tuple)
|
||||
push!(timesteps, TimeStep(Float64(timestep[1]), Increment(timestep[2])))
|
||||
else
|
||||
push!(timesteps, TimeStep(Float64(i-1), Increment(timestep)))
|
||||
end
|
||||
end
|
||||
new_field = Field(field.time, similar(field.values))
|
||||
data = reshape(data, fdim, length(field))
|
||||
for i=1:length(new_field)
|
||||
new_field.values[i] = data[:,i]
|
||||
end
|
||||
return new_field
|
||||
return DefaultDiscreteField(timesteps)
|
||||
end
|
||||
|
||||
### BASIS ###
|
||||
|
||||
abstract AbstractBasis
|
||||
|
||||
|
||||
|
||||
""" FieldSet is set of fields, each field can have different time and/or increment. """
|
||||
type FieldSet
|
||||
name :: ASCIIString
|
||||
fields :: Array{Field, 1}
|
||||
end
|
||||
""" Initializer for FieldSet. """
|
||||
function FieldSet(field_name::ASCIIString)
|
||||
FieldSet(field_name, [])
|
||||
end
|
||||
function FieldSet()
|
||||
FieldSet("unknown field", [])
|
||||
end
|
||||
function FieldSet(fields::Array{Field, 1})
|
||||
FieldSet("unknown field", fields)
|
||||
end
|
||||
""" Add new field to fieldset. """
|
||||
function Base.push!(fs::FieldSet, field::Field)
|
||||
push!(fs.fields, field)
|
||||
end
|
||||
""" Multiply fieldset with some vector x. """
|
||||
Base.(:*)(x::Array{Float64, 1}, fs::FieldSet) = sum(x .* fs.fields)
|
||||
""" Get length of a fieldset. """
|
||||
function Base.length(fieldset::FieldSet)
|
||||
length(fieldset.fields)
|
||||
end
|
||||
""" Return ith field from fieldset. """
|
||||
function Base.getindex(fieldset::FieldSet, i::Int64)
|
||||
fieldset.fields[i]
|
||||
end
|
||||
#""" Return last field from fieldset. """
|
||||
function Base.endof(fieldset::FieldSet)
|
||||
length(fieldset)
|
||||
end
|
||||
function Base.convert(fieldset::Type{FieldSet}, field::Field)
|
||||
FieldSet(Field[field])
|
||||
end
|
||||
|
||||
|
||||
""" Basis function. """
|
||||
type Basis
|
||||
""" Defined to dimensionless coordinate ξ∈[-1,1]^n. """
|
||||
type SpatialBasis <: AbstractBasis
|
||||
basis :: Function
|
||||
dbasisdxi :: Function
|
||||
end
|
||||
#""" Constructor of basis function. """
|
||||
#function Basis(basis)
|
||||
# Basis(basis, ForwardDiff.jacobian(basis))
|
||||
#end
|
||||
|
||||
typealias Basis SpatialBasis
|
||||
|
||||
""" Defined to to interval t∈[0, 1]. """
|
||||
type TemporalBasis <: AbstractBasis
|
||||
basis :: Function
|
||||
dbasisdt :: Function
|
||||
end
|
||||
function TemporalBasis()
|
||||
basis(t) = [1-t, t]
|
||||
dbasis(t) = [-1, 1]
|
||||
return TemporalBasis(basis, dbasis)
|
||||
end
|
||||
|
||||
function call(b::TemporalBasis, value::Number)
|
||||
b.basis(value)
|
||||
end
|
||||
|
||||
function call(b::SpatialBasis, value::Vector)
|
||||
b.basis(value)
|
||||
end
|
||||
|
||||
### INTERPOLATION IN TIME DOMAIN ###
|
||||
|
||||
function Base.call(field::Field, basis::TemporalBasis, time)
|
||||
# FieldSet -> Field -> TimeStep -> Increment -> data
|
||||
# special cases, -Inf, +Inf and ~0.0
|
||||
if time > field[end].time
|
||||
return field[end][end]
|
||||
end
|
||||
if (time < field[1].time) || abs(time-field[1].time) < 1.0e-12
|
||||
return field[1][end]
|
||||
end
|
||||
i = length(field)
|
||||
while field[i].time >= time
|
||||
i -= 1
|
||||
end
|
||||
field[i].time == time && return field[i][end]
|
||||
t1 = field[i].time
|
||||
t2 = field[i+1].time
|
||||
inc1 = field[i][end]
|
||||
inc2 = field[i+1][end]
|
||||
# TODO: may there be some reasons for "unphysical" jumps in
|
||||
# fields w.r.t time which should be taken account in some way?
|
||||
# i.e. dt between two fields → 0
|
||||
dt = t2 - t1
|
||||
b = basis.basis((time-t1)/dt)
|
||||
r = Increment[inc1, inc2]
|
||||
return dot(b, r)
|
||||
end
|
||||
function Base.call(field::DiscreteField, time)
|
||||
return Base.call(field, TemporalBasis(), time)
|
||||
end
|
||||
|
||||
function Base.call(field::Field, basis::TemporalBasis, time,
|
||||
derivative::Type{Val{:derivative}})
|
||||
# FieldSet -> Field -> TimeStep -> Increment -> data
|
||||
|
||||
if length(field) == 1
|
||||
# just one timestep, time derivative cannot be evaluated.
|
||||
error("Field length = $(length(field)), cannot evaluate time derivative")
|
||||
end
|
||||
|
||||
function eval_field(i, j)
|
||||
timesteps = TimeStep[field[i], field[j]]
|
||||
increments = Increment[timesteps[1][end], timesteps[2][end]]
|
||||
J = norm(timesteps[2].time - timesteps[1].time)
|
||||
dbasisdt = basis.dbasisdt( (time-timesteps[1].time)/J )
|
||||
return dot(dbasisdt, increments)/J
|
||||
end
|
||||
|
||||
# special cases, +Inf, -Inf, ~0.0
|
||||
if (time > field[end].time) || isapprox(time, field[end].time)
|
||||
return eval_field(endof(field)-1, endof(field))
|
||||
end
|
||||
if (time < field[1].time) || isapprox(time, field[1].time)
|
||||
return eval_field(1, 2)
|
||||
end
|
||||
|
||||
# search for a correct "bin" between time steps
|
||||
i = length(field)
|
||||
#while field[i].time >= time + 1.0e-12
|
||||
while (field[i].time > time) && !isapprox(field[i].time, time)
|
||||
i -= 1
|
||||
end
|
||||
|
||||
if isapprox(field[i].time, time)
|
||||
# This is the hard case, maybe discontinuous time
|
||||
# derivative if linear approximation.
|
||||
# we are on the "mid node" in time axis
|
||||
field1 = eval_field(i-1,i)
|
||||
field2 = eval_field(i,i+1)
|
||||
return 1/2*(field1 + field2)
|
||||
end
|
||||
|
||||
return eval_field(i, i+1)
|
||||
|
||||
end
|
||||
|
||||
### INTERPOLATION IN SPATIAL DOMAIN ###
|
||||
|
||||
function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector)
|
||||
basis = basis.basis(xi)
|
||||
sum([basis[i]*increment[i] for i=1:length(increment)])
|
||||
end
|
||||
|
||||
function Base.call(increment::Increment, basis::SpatialBasis, xi::Vector,
|
||||
geometry::Increment, gradient::Type{Val{:gradient}})
|
||||
dbasis = basis.dbasisdxi(xi)
|
||||
J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
|
||||
grad = inv(J)*dbasis
|
||||
gradf = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
|
||||
return gradf
|
||||
end
|
||||
|
||||
### INTEGRATIONPOINT ###
|
||||
|
||||
"""
|
||||
Integration point
|
||||
@@ -160,7 +282,7 @@ attributes :: Dict{Any, Any}
|
||||
material models.
|
||||
"""
|
||||
type IntegrationPoint
|
||||
xi :: Array{Float64, 1}
|
||||
xi :: Vector
|
||||
weight :: Float64
|
||||
fields :: Dict{ASCIIString, FieldSet}
|
||||
end
|
||||
@@ -168,20 +290,6 @@ function IntegrationPoint(xi, weight)
|
||||
IntegrationPoint(xi, weight, Dict())
|
||||
end
|
||||
|
||||
call(b::SpatialBasis, ip::IntegrationPoint) = b.basis(ip.xi)
|
||||
|
||||
|
||||
|
||||
# convenient functions -- maybe this is not correct place for them
|
||||
""" Evaluate basis function in point ξ. """
|
||||
call(b::Basis, xi::Vector) = b.basis(xi)
|
||||
call(b::Basis, ip::IntegrationPoint) = b.basis(ip.xi)
|
||||
Base.(:*)(basis::Basis, fs::FieldSet) = (xi, t) -> basis(xi)*fs(t)
|
||||
|
||||
#""" Interpolate field (h*f)(ξ) """
|
||||
#Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
|
||||
#""" Interpolate from set of fields with basis b, i.e. f(t) = b(t)*[f1, f2] """
|
||||
#Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
|
||||
#""" Interpolate field f using basis b. """
|
||||
#Base.(:*)(b::Basis, f::Field) = (x) -> b(x)*f
|
||||
#Base.(:*)(b::Basis, f::Array{Field}) = (t) -> b(t)*f
|
||||
|
||||
|
||||
+1
-1
@@ -155,7 +155,7 @@ function xdmf_new_field(grid, name, source, data)
|
||||
|
||||
typ = string(typeof(data))
|
||||
datatype = "unknown"
|
||||
@debug("typeof: ", typ)
|
||||
@debug("typeof: $typ")
|
||||
for j in ["Int", "Float"]
|
||||
@debug(j)
|
||||
if contains(typ, j)
|
||||
|
||||
+10
-9
@@ -1,26 +1,27 @@
|
||||
# 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 FactCheck
|
||||
using Base.Test
|
||||
|
||||
using JuliaFEM: Quad4, Field, FieldSet, CPS4, get_basis, solve!, PlaneStressElasticityProblem
|
||||
|
||||
|
||||
facts("test plane elasticity on single element, volume load") do
|
||||
function run()
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]))
|
||||
element["youngs modulus"] = FieldSet(Field(500.0))
|
||||
element["poissons ratio"] = FieldSet(Field(0.3))
|
||||
element["displacement load"] = FieldSet(Field(0.0, Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]))
|
||||
element["geometry"] = Vector[[0.0, 0.0], [10.0, 0.0], [10.0, 1.0], [0.0, 1.0]]
|
||||
element["youngs modulus"] = 500.0
|
||||
element["poissons ratio"] = 0.3
|
||||
element["displacement load"] = Vector[[0.0, -10.0], [0.0, -10.0], [0.0, -10.0], [0.0, -10.0]]
|
||||
equation = CPS4(element)
|
||||
free_dofs = [3, 4, 5, 6]
|
||||
problem = PlaneStressElasticityProblem([equation])
|
||||
solve!(problem, free_dofs; max_iterations=10)
|
||||
#solve!(equation, "displacement", free_dofs; max_iterations=10)
|
||||
disp = get_basis(element)("displacement", [1.0, 1.0])[2]
|
||||
Logging.info("displacement at tip: $disp")
|
||||
info("displacement at tip: $disp")
|
||||
# verified using Code Aster.
|
||||
@fact disp --> roughly(-8.77303119819776E+00)
|
||||
@test disp ≈ -8.77303119819776
|
||||
end
|
||||
|
||||
run()
|
||||
|
||||
@@ -0,0 +1,293 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
|
||||
module TypesTests
|
||||
|
||||
using JuliaFEM: Increment, TimeStep, AbstractField, DefaultDiscreteField, FieldSet
|
||||
using JuliaFEM: TemporalBasis, SpatialBasis, ContinuousField, DiscreteField
|
||||
using JuliaFEM: Field
|
||||
|
||||
|
||||
using Base.Test
|
||||
|
||||
function test_increment()
|
||||
info("testing Increment")
|
||||
# testing Increment
|
||||
I1 = Increment([1, 2, 3])
|
||||
I2 = Increment([2, 3, 4])
|
||||
@test dot(I1, I2) == 20
|
||||
@test dot([1,2,3], I2) == 20
|
||||
@test dot(I1, [2,3,4]) == 20
|
||||
@test 1/2*(I1+I2) == [1.5, 2.5, 3.5]
|
||||
@test I1 + 1 == [2, 3, 4]
|
||||
@test I1 - 1 == [0, 1, 2]
|
||||
@test I1*3 == [3, 6, 9]
|
||||
@test I1+I2 == [3, 5, 7]
|
||||
f = zeros(Increment, 2, 4)
|
||||
@test length(f) == 4
|
||||
|
||||
g = similar(f, ones(8))
|
||||
@test typeof(f) == typeof(g)
|
||||
@test length(f) == length(g)
|
||||
|
||||
# promotion of increment
|
||||
@test typeof(I1+1) == typeof(I1)
|
||||
@test typeof(I1-1) == typeof(I1)
|
||||
@test typeof(I1*3) == typeof(I1)
|
||||
|
||||
# FIXME
|
||||
#@test typeof(I1) == typeof(I1+I2)
|
||||
#@test typeof(I1/2) == typeof(I1)
|
||||
#@test typeof(1/2*S1) == typeof(I1)
|
||||
end
|
||||
test_increment()
|
||||
|
||||
function test_timestep()
|
||||
info("testing TimeStep")
|
||||
i1 = Increment([1, 2, 3])
|
||||
i2 = Increment([2, 3, 4])
|
||||
i3 = Increment([2, 3, 4])
|
||||
i4 = Increment([3, 4, 5])
|
||||
t1 = TimeStep(1.0, Increment[i1, i2])
|
||||
t2 = TimeStep(2.0, Increment[i3, i4])
|
||||
@test length(t1) == length(t2) == 2
|
||||
t3 = TimeStep(3.0, i1+1)
|
||||
end
|
||||
test_timestep()
|
||||
|
||||
function test_watta_fak()
|
||||
# TODO: this test will fail if Increment is typealiased to Vector
|
||||
fs = FieldSet()
|
||||
fs["discrete field"] = [1, 2, 3, 4]
|
||||
T0 = last(fs["discrete field"])
|
||||
info("last discrete field: $T0, ", typeof(T0))
|
||||
T1 = T0 + 1
|
||||
info("adding 1 to discrete field: $T1, ", typeof(T1))
|
||||
ts = TimeStep(1.0, T1)
|
||||
info("creating time step: $ts, ", typeof(ts))
|
||||
push!(fs["discrete field"], ts)
|
||||
info("last discrete field = ", last(fs["discrete field"]))
|
||||
|
||||
info("fieldset: $fs")
|
||||
|
||||
@test last(fs["discrete field"]) == [2, 3, 4, 5]
|
||||
end
|
||||
test_watta_fak()
|
||||
|
||||
function test_default_discrete_field()
|
||||
info("testing DefaultDiscreteField")
|
||||
i1 = Increment([1, 2, 3])
|
||||
i2 = Increment([2, 3, 4])
|
||||
i3 = Increment([2, 3, 4])
|
||||
i4 = Increment([3, 4, 5])
|
||||
t1 = TimeStep(1.0, Increment[i1, i2])
|
||||
t2 = TimeStep(2.0, Increment[i3, i4])
|
||||
timesteps = TimeStep[t1, t2]
|
||||
f1 = DefaultDiscreteField(timesteps)
|
||||
@test length(f1) == 2
|
||||
@test isa(f1, AbstractField) == true
|
||||
end
|
||||
test_default_discrete_field()
|
||||
|
||||
function test_fieldset()
|
||||
i1 = Increment([1, 2, 3])
|
||||
i2 = Increment([2, 3, 4])
|
||||
i3 = Increment([2, 3, 4])
|
||||
i4 = Increment([3, 4, 5])
|
||||
t1 = TimeStep(1.0, Increment[i1, i2])
|
||||
t2 = TimeStep(2.0, Increment[i3, i4])
|
||||
timesteps = TimeStep[t1, t2]
|
||||
f1 = DefaultDiscreteField(timesteps)
|
||||
info("testing adding discrete field to FieldSet")
|
||||
fs = FieldSet()
|
||||
fs["temperature"] = f1
|
||||
@test length(fs) == 1
|
||||
|
||||
info("testing adding discrete fields quickly")
|
||||
# the easy way
|
||||
fs2 = FieldSet()
|
||||
fs2["temperature"] = [1, 2, 3, 4]
|
||||
@test fs2["temperature"][end][end] == [1, 2, 3, 4]
|
||||
@test last(fs2["temperature"]) == [1, 2, 3, 4]
|
||||
|
||||
fs2 = FieldSet()
|
||||
fs2["constant scalar field"] = 1
|
||||
fs2["scalar field"] = [1, 2, 3, 4]
|
||||
fs2["vector field"] = reshape(collect(1:8), 2, 4)
|
||||
fs2["second order tensor field"] = reshape(collect(1:3*3*4), 3, 3, 4)
|
||||
fs2["fourth order tensor field"] = reshape(collect(1:3*3*3*3*4), 3, 3, 3, 3, 4)
|
||||
timestep = fs2["vector field"][end]
|
||||
@test timestep.time == 0.0
|
||||
|
||||
info("testing adding timesteps")
|
||||
# add another timestep
|
||||
fs = FieldSet()
|
||||
fs["temperature"] = [1, 2, 3, 4]
|
||||
T0 = last(fs["temperature"]) # last increment of last field
|
||||
info("last temperature = $T0")
|
||||
T1 = T0 + 1
|
||||
@test typeof(T0) == typeof(T1)
|
||||
timestep = TimeStep(1.0, Increment[T1]) # new list of increments for timestep
|
||||
push!(fs["temperature"], timestep)
|
||||
T2 = last(fs["temperature"])
|
||||
info("last temperature = $T2")
|
||||
@test last(fs["temperature"]) == [2, 3, 4, 5]
|
||||
# or more easily
|
||||
timestep = TimeStep(2.0, T1)
|
||||
push!(fs["temperature"], timestep)
|
||||
@test length(fs["temperature"].timesteps) == 3
|
||||
|
||||
info("test adding several time steps at once")
|
||||
fs3 = FieldSet()
|
||||
fs3["time series 1"] = (0.0, [1, 2, 3, 4]), (0.5, [2, 3, 4, 5]), (1.0, [1, 1, 1, 1])
|
||||
@test fs3["time series 1"][end].time == 1.0
|
||||
fs3["time series 2"] = [1, 2, 3, 4], [2, 3, 4, 5], [1, 1, 1, 1]
|
||||
@test fs3["time series 2"][end].time == 2.0
|
||||
end
|
||||
test_fieldset()
|
||||
|
||||
type MyFunnyContinuousField <: ContinuousField
|
||||
basis :: Function
|
||||
discretefield :: DiscreteField
|
||||
end
|
||||
function Base.call(field::MyFunnyContinuousField, xi::Vector, time::Number=1.0)
|
||||
data = last(field.discretefield) # get the last timestep last increment
|
||||
info("data = $data, typeof data = $(typeof(data))")
|
||||
basis = time*field.basis(xi) # evaluate basis at point ξ.
|
||||
sum([basis[i]*data[i] for i=1:length(data)]) # sum results
|
||||
end
|
||||
function test_continuous_field()
|
||||
info("testing continuous field")
|
||||
fs = FieldSet()
|
||||
fs["discrete field"] = [1, 2, 3, 4]
|
||||
basis(xi) = 1/4*[
|
||||
(1-xi[1])*(1-xi[2]),
|
||||
(1+xi[1])*(1-xi[2]),
|
||||
(1+xi[1])*(1+xi[2]),
|
||||
(1-xi[1])*(1+xi[2])]
|
||||
fs["continuous field"] = MyFunnyContinuousField(basis, fs["discrete field"])
|
||||
@test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(1+2+3+4)
|
||||
T0 = last(fs["discrete field"])
|
||||
T1 = T0 + 1.0
|
||||
ts = TimeStep(1.0, T1)
|
||||
push!(fs["discrete field"], TimeStep(1.0, T0+1.0))
|
||||
@test fs["continuous field"]([0.0, 0.0], 1.0) == 1/4*(2+3+4+5)
|
||||
end
|
||||
test_continuous_field()
|
||||
|
||||
type MyFunnyDiscreteField <: DiscreteField
|
||||
discrete_points :: Vector
|
||||
continuousfield :: ContinuousField
|
||||
end
|
||||
Base.length(field::MyFunnyDiscreteField) = length(field.discrete_points)
|
||||
Base.endof(field::MyFunnyDiscreteField) = endof(field.discrete_points)
|
||||
Base.last(field::MyFunnyDiscreteField) = Float64[field[i] for i=1:length(field)]
|
||||
function Base.getindex(field::MyFunnyDiscreteField, idx::Int64)
|
||||
field.continuousfield(field.discrete_points[idx])
|
||||
end
|
||||
|
||||
function test_discrete_field()
|
||||
info("testing discrete field")
|
||||
fs = FieldSet()
|
||||
fs["discrete field"] = [1, 2, 3, 4]
|
||||
basis(xi) = 1/4*[
|
||||
(1-xi[1])*(1-xi[2]),
|
||||
(1+xi[1])*(1-xi[2]),
|
||||
(1+xi[1])*(1+xi[2]),
|
||||
(1-xi[1])*(1+xi[2])]
|
||||
fs["continuous field"] = MyFunnyContinuousField(basis, fs["discrete field"])
|
||||
discrete_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]
|
||||
fs["discrete field 2"] = MyFunnyDiscreteField(discrete_points, fs["continuous field"])
|
||||
@test last(fs["discrete field 2"]) ≈ [
|
||||
1.7559830641437073,
|
||||
2.0893163974770410,
|
||||
2.9106836025229590,
|
||||
3.2440169358562922]
|
||||
end
|
||||
test_discrete_field()
|
||||
|
||||
|
||||
function test_interpolation_in_temporal_basis()
|
||||
info("testing interpolation on temporal basis")
|
||||
temporalbasis = TemporalBasis((t) -> [1-t, t], (t) -> [-1, 1])
|
||||
@test temporalbasis(0.2) == [0.8, 0.2]
|
||||
i1 = Increment([0.0])
|
||||
i2 = Increment([1.0])
|
||||
i3 = Increment([2.0])
|
||||
t1 = TimeStep(0.0, Increment[i1])
|
||||
t2 = TimeStep(2.0, Increment[i2])
|
||||
t3 = TimeStep(4.0, Increment[i3])
|
||||
field = Field(TimeStep[t1, t2, t3])
|
||||
@test call(field, temporalbasis, -Inf) == [0.0]
|
||||
@test call(field, temporalbasis, 0.0) == [0.0]
|
||||
@test call(field, temporalbasis, 1.0) == [0.5]
|
||||
@test call(field, temporalbasis, 2.0) == [1.0]
|
||||
@test call(field, temporalbasis, 3.0) == [1.5]
|
||||
@test call(field, temporalbasis, 4.0) == [2.0]
|
||||
@test call(field, temporalbasis, +Inf) == [2.0]
|
||||
@test call(field, temporalbasis, +Inf, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, -Inf, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, 0.0, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, 0.5, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, 1.0, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, 1.5, Val{:derivative}) == [0.5]
|
||||
@test call(field, temporalbasis, 2.0, Val{:derivative}) == [0.5]
|
||||
fs = FieldSet()
|
||||
|
||||
t = collect(linspace(0, 2, 5))
|
||||
x = 1/2*t.^2
|
||||
x2 = tuple(collect(zip(t, x))...)
|
||||
# => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
|
||||
fs["particle"] = x2
|
||||
position = call(fs["particle"], temporalbasis, 1.0)[1]
|
||||
@test position ≈ 0.50
|
||||
velocity = call(fs["particle"], temporalbasis, 2.0, Val{:derivative})[1]
|
||||
@test velocity ≈ (2.0-1.125)/0.5 # = 1.75
|
||||
velocity = call(fs["particle"], temporalbasis, 1.0, Val{:derivative})[1]
|
||||
v1 = (0.500 - 0.125)/0.5
|
||||
v2 = (1.125 - 0.500)/0.5
|
||||
info("v1 = $v1, v2 = $v2")
|
||||
info(mean([v1, v2]))
|
||||
@test velocity ≈ mean([v1, v2]) # = 1.00
|
||||
|
||||
# FIXME, returns wrong type.
|
||||
#=
|
||||
@test isa(position, Increment) == true
|
||||
@test isa(velocity, Increment) == true
|
||||
=#
|
||||
end
|
||||
test_interpolation_in_temporal_basis()
|
||||
|
||||
function test_interpolation_in_spatial_basis()
|
||||
info("testing interpolation on spatial basis")
|
||||
basis(xi) = 1/4*[
|
||||
(1-xi[1])*(1-xi[2])
|
||||
(1+xi[1])*(1-xi[2])
|
||||
(1+xi[1])*(1+xi[2])
|
||||
(1-xi[1])*(1+xi[2])]'
|
||||
dbasis(xi) = 1/4*[
|
||||
-(1-xi[2]) (1-xi[2]) (1+xi[2]) -(1+xi[2])
|
||||
-(1-xi[1]) -(1+xi[1]) (1+xi[1]) (1-xi[1])]
|
||||
spatialbasis = SpatialBasis(basis, dbasis)
|
||||
@test spatialbasis.basis([0.0, 0.0]) == 1/4*[1 1 1 1]
|
||||
|
||||
fs = FieldSet()
|
||||
fs["geometry"] = Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]]
|
||||
fs["displacement"] = (0.0, zeros(2, 4)), (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
|
||||
|
||||
X = call(last(fs["geometry"]), spatialbasis, [0.0, 0.0])
|
||||
u = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0])
|
||||
x = X+u
|
||||
@test X ≈ 1/2*[1, 1]
|
||||
@test x ≈ [9/16, 1/2]
|
||||
|
||||
gradu = call(last(fs["displacement"]), spatialbasis, [0.0, 0.0], last(fs["geometry"]), Val{:gradient})
|
||||
@test gradu ≈ [0.125 0.125; 0.0 0.0]
|
||||
end
|
||||
test_interpolation_in_spatial_basis()
|
||||
|
||||
|
||||
println("test_fields.jl: all test passing.")
|
||||
|
||||
end
|
||||
+6
-6
@@ -10,16 +10,16 @@ facts("tests on [0x1]x[0x1] domain") do
|
||||
|
||||
# volume element
|
||||
element = Quad4([1, 2, 3, 4])
|
||||
element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]))
|
||||
element["temperature thermal conductivity"] = FieldSet(Field(0.0, 6.0))
|
||||
element["temperature load"] = FieldSet(Field(0.0, [12.0, 12.0, 12.0, 12.0]))
|
||||
element["density"] = FieldSet(Field(0.0, 36.0))
|
||||
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
|
||||
|
||||
# boundary element
|
||||
boundary_element = Seg2([1, 2])
|
||||
boundary_element["geometry"] = FieldSet(Field(Vector[[0.0, 0.0], [1.0, 0.0]]))
|
||||
boundary_element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0]]
|
||||
# linear ramp from 1 to 6 in time 0 to 1
|
||||
boundary_element["temperature flux"] = FieldSet(Field[Field(0.0, 0.0), Field(1.0, 6.0)])
|
||||
boundary_element["temperature flux"] = (0.0, 0.0), (1.0, 6.0)
|
||||
|
||||
# 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)
|
||||
|
||||
Reference in New Issue
Block a user