mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
defining basis. still needs some rethinking...
This commit is contained in:
Jukka Aho
@@ -4,3 +4,9 @@
|
||||
docs/build/html
|
||||
*.swp
|
||||
*.lnk
|
||||
*.mess
|
||||
*.log
|
||||
*.aux
|
||||
*.out
|
||||
*.tex
|
||||
*.export
|
||||
|
||||
+175
-60
@@ -1,88 +1,219 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
abstract AbstractBasis
|
||||
abstract Basis <: ContinuousField
|
||||
|
||||
""" Defined to dimensionless coordinate ξ∈[-1,1]^n. """
|
||||
type SpatialBasis <: AbstractBasis
|
||||
### ELEMENT BASIS
|
||||
|
||||
""" This is the normal "user defined" basis functions familiar from school books. """
|
||||
type ElementBasis <: Basis
|
||||
basis :: Function
|
||||
dbasisdxi :: Function
|
||||
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)
|
||||
function Basis(basis::Function, dbasisdxi::Function)
|
||||
return ElementBasis(basis, dbasisdxi)
|
||||
end
|
||||
|
||||
function call(b::TemporalBasis, value::Number)
|
||||
b.basis(value)
|
||||
function Base.call(basis::ElementBasis, xi::Vector, time::Number=0.0)
|
||||
basis.basis(xi) # passing time does not make much sense actually for this...
|
||||
end
|
||||
|
||||
function call(b::SpatialBasis, value::Vector)
|
||||
b.basis(value)
|
||||
""" Interpolate increment in spatial domain using ElementBasis. """
|
||||
function Base.call(basis::ElementBasis, increment::Increment, xi::Vector)
|
||||
basis = basis.basis(xi)
|
||||
sum([basis[i]*increment[i] for i=1:length(increment)])
|
||||
end
|
||||
|
||||
### INTERPOLATION IN TIME DOMAIN ###
|
||||
### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD
|
||||
|
||||
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
|
||||
""" Here we add field we are wanting to interpolate with ElementBasis. """
|
||||
type ElementFieldBasis <: Basis
|
||||
element_basis :: ElementBasis
|
||||
field :: DiscreteField
|
||||
time_extrapolation :: Symbol
|
||||
time_interpolation :: Symbol
|
||||
end
|
||||
|
||||
function Basis(basis::Function, dbasisdxi::Function, field::DiscreteField,
|
||||
time_extrapolation=:linear, time_interpolation=:linear)
|
||||
element_basis = ElementBasis(basis, dbasisdxi)
|
||||
return ElementFieldBasis(element_basis, field, time_extrapolation,
|
||||
time_interpolation)
|
||||
end
|
||||
|
||||
function Base.call(basis::ElementFieldBasis, xi::Vector, time::Number)
|
||||
increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation)
|
||||
return basis.element_basis(increment, xi)
|
||||
end
|
||||
|
||||
""" Interpolate discrete field in time domain. """
|
||||
function Base.call(field::DiscreteField, time::Number,
|
||||
time_extrapolation::Symbol=:linear,
|
||||
time_interpolation::Symbol=:linear)
|
||||
|
||||
# special cases, only 1 timestep defined or time = -Inf -> return first ts
|
||||
if (length(field) == 1) || (time == -Inf)
|
||||
return field[1][end]
|
||||
end
|
||||
|
||||
# special case, time = +Inf -> return last ts
|
||||
if time == +Inf
|
||||
return field[end][end]
|
||||
end
|
||||
|
||||
# very likely we are always near some defined timestep, usually field
|
||||
# defined only on t = 0.0, test neighbourhood for timesteps
|
||||
for i=1:length(field)
|
||||
if isapprox(field[i].time, time)
|
||||
return field[i][end]
|
||||
end
|
||||
end
|
||||
|
||||
# special case: out of time domain in positive direction, very likely
|
||||
# to happen in incremental constitutive models
|
||||
if time > field[end].time
|
||||
if time_extrapolation == :constant
|
||||
# constant time extrapolation, return last field
|
||||
return field[end][end]
|
||||
else
|
||||
# multiple fields, pick last and second last and do linear interpolation
|
||||
f1 = field[end-1]
|
||||
f2 = field[end]
|
||||
dt = abs(f2.time - f1.time)
|
||||
i1 = f1[end]
|
||||
i2 = f2[end]
|
||||
di = i2 - i1
|
||||
increment = Increment(i2 + di./dt * (time-f2.time))
|
||||
return increment
|
||||
end
|
||||
end
|
||||
|
||||
# special case: out of time domain in negative direction
|
||||
if time < field[1].time
|
||||
if time_extrapolation == :constant
|
||||
# constant time extrapolation, return first field
|
||||
return field[1][end]
|
||||
else
|
||||
# multiple fields, pick first and second and do linear interpolation
|
||||
f1 = field[1]
|
||||
f2 = field[2]
|
||||
dt = abs(f2.time - f1.time)
|
||||
i1 = f1[end]
|
||||
i2 = f2[end]
|
||||
di = i2 - i1
|
||||
increment = Increment(i1 - di./dt * (f1.time - time))
|
||||
return increment
|
||||
end
|
||||
end
|
||||
|
||||
# find correct bin and perform interpolation
|
||||
i = length(field)
|
||||
while field[i].time >= time
|
||||
i -= 1
|
||||
end
|
||||
field[i].time == time && return field[i][end]
|
||||
|
||||
if time_interpolation == :linear
|
||||
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)
|
||||
dt = abs(t2 - t1)
|
||||
t = (time-t1)/dt
|
||||
increment = Increment((1-t)*inc1 + t*inc2)
|
||||
return increment
|
||||
end
|
||||
|
||||
if time_interpolation == :constant
|
||||
# nearest neightbour interpolation, i.e. pick nearest defined field
|
||||
t1 = field[i].time
|
||||
t2 = field[i+1].time
|
||||
dt1 = abs(t1-time)
|
||||
dt2 = abs(t2-time)
|
||||
if dt1 < dt2
|
||||
return field[i][end]
|
||||
else
|
||||
return field[i+1][end]
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
|
||||
function Base.call(field::Field, basis::TemporalBasis, time,
|
||||
derivative::Type{Val{:derivative}})
|
||||
### ELEMENT GRADIENT BASIS = ELEMENT BASIS + GEOMETRY
|
||||
|
||||
""" Gradient of ElementBasis, needs geometry information. """
|
||||
type ElementGradientBasis <: Basis
|
||||
element_basis :: ElementBasis
|
||||
geometry :: DiscreteField
|
||||
time_extrapolation :: Symbol
|
||||
time_interpolation :: Symbol
|
||||
end
|
||||
|
||||
function ElementGradientBasis(element_basis::ElementBasis, geometry::DiscreteField)
|
||||
return ElementGradientBasis(element_basis, geometry, :linear, :linear)
|
||||
end
|
||||
|
||||
function Base.call(basis::ElementGradientBasis, xi::Vector, time::Number=0.0)
|
||||
dbasis = basis.element_basis.dbasisdxi(xi)
|
||||
geometry = basis.geometry(time, basis.time_extrapolation, basis.time_interpolation)
|
||||
J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
|
||||
grad = inv(J)*dbasis
|
||||
return grad
|
||||
end
|
||||
|
||||
### ELEMENT FIELD GRADIENT BASIS = ELEMENT GRADIENT BASIS + FIELD
|
||||
|
||||
""" Gradient of ElementFieldBasis, needs field to interpolate. """
|
||||
type ElementFieldGradientBasis <: Basis
|
||||
element_gradient_basis :: ElementGradientBasis
|
||||
field :: DiscreteField
|
||||
time_extrapolation :: Symbol
|
||||
time_interpolation :: Symbol
|
||||
end
|
||||
|
||||
function ElementFieldGradientBasis(element_gradient_basis::ElementGradientBasis,
|
||||
field::DiscreteField)
|
||||
return ElementFieldGradientBasis(element_gradient_basis, field, :linear, :linear)
|
||||
end
|
||||
|
||||
function Base.call(basis::ElementFieldGradientBasis, xi::Vector, time::Number=0.0)
|
||||
grad = basis.element_gradient_basis(xi, time)
|
||||
increment = basis.field(time, basis.time_extrapolation, basis.time_interpolation)
|
||||
gradf = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
|
||||
return gradf
|
||||
end
|
||||
|
||||
### INTERPOLATION IN TIME DOMAIN ###
|
||||
|
||||
function Base.call(field::DiscreteField, time::Number,
|
||||
derivative::Type{Val{:derivative}},
|
||||
time_extrapolation::Symbol=:linear,
|
||||
time_interpolation::Symbol=:linear)
|
||||
|
||||
# FieldSet -> Field -> TimeStep -> Increment -> data
|
||||
|
||||
time_extrapolation == :linear || error("$time_extrapolation not implemented")
|
||||
time_interpolation == :linear || error("$time_interpolation not implemented")
|
||||
|
||||
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
|
||||
t1 = field[i]
|
||||
t2 = field[j]
|
||||
J = abs(t2.time - t1.time)
|
||||
t = (time-t1.time)/J
|
||||
result = 1/J*((1-t)*t1[end] + t*t2[end])
|
||||
return Increment(result)
|
||||
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
|
||||
@@ -106,19 +237,3 @@ function Base.call(field::Field, basis::TemporalBasis, time,
|
||||
|
||||
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
|
||||
|
||||
|
||||
+21
-5
@@ -131,7 +131,7 @@ end
|
||||
# FIXME: having some serious problems here to get tuple form working.
|
||||
|
||||
# 3. DefaultDiscreteField
|
||||
type DefaultDiscreteField <: DiscreteField
|
||||
immutable DefaultDiscreteField <: DiscreteField
|
||||
timesteps :: Vector{TimeStep}
|
||||
#=
|
||||
function DefaultDiscreteField(data::Array)
|
||||
@@ -170,6 +170,26 @@ function Base.size(field::DefaultDiscreteField)
|
||||
return size(field.timesteps)
|
||||
end
|
||||
|
||||
function Base.length(field::DefaultDiscreteField)
|
||||
return length(field.timesteps)
|
||||
end
|
||||
|
||||
function Base.start(::DefaultDiscreteField)
|
||||
return 1
|
||||
end
|
||||
|
||||
function Base.next(field::DefaultDiscreteField, state)
|
||||
return (field[state+1], state+1)
|
||||
end
|
||||
|
||||
function Base.done(field::DefaultDiscreteField, state)
|
||||
return state > length(field)
|
||||
end
|
||||
|
||||
function eltype(::Type{DefaultDiscreteField})
|
||||
return TimeStep
|
||||
end
|
||||
|
||||
function Base.linearindexing(::Type{DefaultDiscreteField})
|
||||
return LinearFast()
|
||||
end
|
||||
@@ -178,10 +198,6 @@ function Base.getindex(field::DefaultDiscreteField, i::Int)
|
||||
return field.timesteps[i]
|
||||
end
|
||||
|
||||
function Base.length(field::DefaultDiscreteField)
|
||||
return length(field.timesteps)
|
||||
end
|
||||
|
||||
function Base.endof(field::DefaultDiscreteField)
|
||||
return endof(field.timesteps)
|
||||
end
|
||||
|
||||
+1
-1
@@ -23,6 +23,6 @@ function IntegrationPoint(xi, weight)
|
||||
IntegrationPoint(xi, weight, Dict())
|
||||
end
|
||||
|
||||
call(b::SpatialBasis, ip::IntegrationPoint) = b.basis(ip.xi)
|
||||
call(N::ElementBasis, ip::IntegrationPoint) = N(ip.xi)
|
||||
|
||||
|
||||
|
||||
+157
-31
@@ -5,9 +5,163 @@ module BasisTests
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
using JuliaFEM: get_basis, grad, FieldSet, Field, Quad4
|
||||
using JuliaFEM
|
||||
using JuliaFEM: Basis, ElementGradientBasis, ElementFieldGradientBasis, Field
|
||||
|
||||
function get_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])]
|
||||
|
||||
return basis, dbasis
|
||||
end
|
||||
|
||||
|
||||
function test_basic_interpolation()
|
||||
basis, dbasis = get_basis()
|
||||
b = Basis(basis, dbasis)
|
||||
@test b([0.0, 0.0]) == 1/4*[1 1 1 1]
|
||||
@test b([0.0, 0.0], 1.0) == 1/4*[1 1 1 1]
|
||||
end
|
||||
|
||||
function test_basic_interpolation_of_field()
|
||||
# in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
|
||||
temperature = Field(
|
||||
(0.0, [0.0, 0.0, 0.0, 0.0]),
|
||||
(1.0, [1.0, 2.0, 3.0, 4.0]))
|
||||
basis, dbasis = get_basis()
|
||||
b = Basis(basis, dbasis, temperature)
|
||||
T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
|
||||
@test b([0.0, 0.0], 0.0) == T([0.5, 0.5], 0.0)
|
||||
@test b([0.0, 0.0], 0.6) == T([0.5, 0.5], 0.6)
|
||||
@test b([0.0, 0.0], 1.0) == T([0.5, 0.5], 1.0)
|
||||
end
|
||||
|
||||
function test_linear_time_extrapolation_of_field()
|
||||
temperature = Field(
|
||||
(0.0, [0.0, 0.0, 0.0, 0.0]),
|
||||
(1.0, [1.0, 2.0, 3.0, 4.0]))
|
||||
basis, dbasis = get_basis()
|
||||
b = Basis(basis, dbasis, temperature, :linear)
|
||||
T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
|
||||
@test b([0.0, 0.0], -1.0) == T([0.5, 0.5], -1.0)
|
||||
@test b([0.0, 0.0], 3.0) == T([0.5, 0.5], 3.0)
|
||||
end
|
||||
|
||||
function test_constant_time_extrapolation_of_field()
|
||||
temperature = Field(
|
||||
(0.0, [0.0, 0.0, 0.0, 0.0]),
|
||||
(1.0, [1.0, 2.0, 3.0, 4.0]))
|
||||
basis, dbasis = get_basis()
|
||||
b = Basis(basis, dbasis, temperature, :constant)
|
||||
T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
|
||||
@test b([0.0, 0.0], -1.0) == T([0.5, 0.5], 0.0)
|
||||
@test b([0.0, 0.0], 3.0) == T([0.5, 0.5], 1.0)
|
||||
end
|
||||
|
||||
function test_time_extrapolation_of_field_with_single_timestep()
|
||||
temperature = Field([1.0, 2.0, 3.0, 4.0])
|
||||
basis, dbasis = get_basis()
|
||||
b = Basis(basis, dbasis, temperature)
|
||||
@test b([0.0, 0.0], 1.0) == mean([1.0, 2.0, 3.0, 4.0])
|
||||
end
|
||||
|
||||
function test_gradient_interpolation_empty_gradient()
|
||||
X = [0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]'
|
||||
geometry = Field(X)
|
||||
# P(X) = [1.0, X[1], X[2], X[1]*X[2]]
|
||||
# basis2, dbasis2 = JuliaFEM.calculate_lagrange_basis(P, X)
|
||||
basis, dbasis = get_basis()
|
||||
N = Basis(basis, dbasis)
|
||||
dN = ElementGradientBasis(N, geometry)
|
||||
@test dN([0.0, 0.0]) == 1/2*[-1 1 1 -1; -1 -1 1 1]
|
||||
# @test dN([0.0, 0.0]) == dbasis2([0.5, 0.5])
|
||||
end
|
||||
|
||||
function test_gradient_interpolation_of_scalar_field()
|
||||
# in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
|
||||
geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
|
||||
|
||||
temperature = Field([1, 2, 3, 4])
|
||||
basis, dbasis = get_basis()
|
||||
N = Basis(basis, dbasis)
|
||||
|
||||
dN = ElementGradientBasis(N, geometry)
|
||||
dT = ElementFieldGradientBasis(dN, temperature)
|
||||
dT_expected(X) = [1-2*X[2] 3-2*X[1]]
|
||||
|
||||
@test dT([0.0, 0.0]) == dT_expected([0.5, 0.5])
|
||||
end
|
||||
|
||||
function test_interpolation_of_vector_field()
|
||||
# in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
|
||||
geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
|
||||
displacement = Field(
|
||||
(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]),
|
||||
(1.0, Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]))
|
||||
|
||||
basis, dbasis = get_basis()
|
||||
X = Basis(basis, dbasis, geometry)
|
||||
u = Basis(basis, dbasis, displacement)
|
||||
u_expected(X,t) = [1/4*t*X[1]*X[2], 0]
|
||||
# x = X + u
|
||||
x = X([0.0, 0.0], 1.0) + u([0.0, 0.0], 1.0)
|
||||
@test isapprox(x, [9/16, 1/2])
|
||||
@test isapprox(u([0.0, 0.0], 1.0), u_expected([0.5, 0.5], 1.0))
|
||||
end
|
||||
|
||||
function test_interpolation_of_gradient_of_vector_field()
|
||||
# in unit square, u(X) = t*[X[1]*X[2]/4, X[1]*(X[1]+X[2])/2]
|
||||
# => u_i,j = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
|
||||
geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
|
||||
displacement = Field(
|
||||
(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]]),
|
||||
(1.0, Vector[[0.0, 0.0], [0.0, 0.5], [0.25, 1.0], [0.0, 0.0]]))
|
||||
|
||||
basis, dbasis = get_basis()
|
||||
N = Basis(basis, dbasis)
|
||||
dN = ElementGradientBasis(N, geometry)
|
||||
dU = ElementFieldGradientBasis(dN, displacement)
|
||||
dU_expected(X, t) = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
|
||||
|
||||
@test isapprox(dU([0.0, 0.0], 1.0), dU_expected([0.5, 0.5], 1.0))
|
||||
end
|
||||
|
||||
# TODO: how on earth make this work without some serious spaghetti code
|
||||
function test_time_derivative_gradient_interpolation_of_field()
|
||||
# in unit square, u(X) = t*[X[1]*X[2]/4, X[1]*(X[1]+X[2])/2]
|
||||
# => u_i,j = t*[X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
|
||||
# => d(u_i,j)/dt = [X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
|
||||
geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
|
||||
displacement = Field(
|
||||
(0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]]),
|
||||
(1.0, Vector[[0.0, 0.0], [0.0, 0.5], [0.25, 1.0], [0.0, 0.0]]))
|
||||
|
||||
basis, dbasis = get_basis()
|
||||
N = Basis(basis, dbasis)
|
||||
# wanted
|
||||
#u = Basis(basis, dbasis, displacement)
|
||||
#L = grad(diff(u))
|
||||
#D = 1/2*(L + L')
|
||||
#@text isapprox(D([0.0, 0.0], 1.0), ...)
|
||||
|
||||
xi = [0.0, 0.0]
|
||||
time = 1.0
|
||||
grad = ElementGradientBasis(N, geometry)(xi, time)
|
||||
increment = displacement(time, Val{:derivative}, :linear, :linear)
|
||||
diffgradu = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
|
||||
diffgradu_expected(X, t) = [X[2]/4 X[1]/4; X[1]/2+(X[1]+X[2])/2 X[1]/2]
|
||||
@test diffgradu == diffgradu_expected([0.5, 0.5], 1.0)
|
||||
end
|
||||
|
||||
#=
|
||||
"""basic continuum interpolations"""
|
||||
function test_basic_interpolations()
|
||||
|
||||
@@ -42,10 +196,9 @@ function test_basic_interpolations()
|
||||
@test isapprox(U, [1.24235 0.13804; 0.13804 1.02149])
|
||||
end
|
||||
|
||||
=#
|
||||
|
||||
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])
|
||||
@@ -90,31 +243,4 @@ function test_interpolation_in_temporal_basis()
|
||||
@test isa(velocity, Increment) == true
|
||||
end
|
||||
|
||||
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
|
||||
|
||||
end
|
||||
|
||||
+18
-6
@@ -6,8 +6,7 @@ module FieldTests
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM: Increment, TimeStep, Field, DefaultDiscreteField, FieldSet
|
||||
using JuliaFEM: TemporalBasis, SpatialBasis, ContinuousField, DiscreteField
|
||||
using JuliaFEM: DefaultContinuousField
|
||||
using JuliaFEM: ContinuousField, DiscreteField, DefaultContinuousField
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
@@ -189,7 +188,9 @@ function test_default_discrete_field_quick_way_two_timesteps_with_vector_value()
|
||||
end
|
||||
|
||||
function test_default_discrete_field_quick_way_set_time_vector_also()
|
||||
f1 = DefaultDiscreteField( (0.5, [1, 2, 3]), (1.0, [3, 4, 5]) )
|
||||
f1 = DefaultDiscreteField(
|
||||
(0.5, [1, 2, 3]),
|
||||
(1.0, [3, 4, 5]))
|
||||
@test isa(f1[1], TimeStep)
|
||||
@test isa(f1[2], TimeStep)
|
||||
@test isa(f1[1][1], Increment)
|
||||
@@ -200,6 +201,20 @@ function test_default_discrete_field_quick_way_set_time_vector_also()
|
||||
@test f1[2].time == 1.0
|
||||
end
|
||||
|
||||
function test_default_discrete_field_for_loop()
|
||||
field = DefaultDiscreteField(
|
||||
(0.5, [1, 2, 3]),
|
||||
(1.0, [3, 4, 5]),
|
||||
(1.5, [4, 5, 6]))
|
||||
timesteps = [ts for ts in field]
|
||||
@test timesteps[1].time == 0.5
|
||||
@test timesteps[2].time == 1.0
|
||||
@test timesteps[3].time == 1.5
|
||||
@test timesteps[1][end] == [1, 2, 3]
|
||||
@test timesteps[2][end] == [3, 4, 5]
|
||||
@test timesteps[3][end] == [4, 5, 6]
|
||||
end
|
||||
|
||||
function test_default_continuous_field()
|
||||
|
||||
function myfield(xi::Vector, time::Float64)
|
||||
@@ -244,13 +259,10 @@ function test_adding_timesteps()
|
||||
fs = FieldSet()
|
||||
fs["temperature"] = [1, 2, 3, 4]
|
||||
T0 = last(fs["temperature"]) # last increment of last field
|
||||
@debug("last temperature T0 = $T0")
|
||||
T1 = Increment(T0 + 1)
|
||||
@debug("typeof T1 = $(typeof(T1))")
|
||||
timestep = TimeStep(1.0, Increment[T1]) # new list of increments for timestep
|
||||
push!(fs["temperature"], timestep)
|
||||
T2 = last(fs["temperature"])
|
||||
@debug("last temperature T2 = $T2")
|
||||
@test length(fs["temperature"]) == 2
|
||||
@test last(fs["temperature"]) == [2, 3, 4, 5]
|
||||
@test fs["temperature"][end].time == 1.0
|
||||
|
||||
Reference in New Issue
Block a user