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

data structures ready

Browse Source
This commit is contained in:
Jukka Aho
2015-11-03 12:14:19 +02:00
parent 1e9d33c0f0
commit 7a36a0ac3d
4 changed files with 348 additions and 278 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/tutorials/2015-06-14-data-structures.ipynb
+103 -37
View File
@@ -707,7 +707,7 @@
{
"data": {
"text/plain": [
"call (generic function with 1275 methods)"
"call (generic function with 1258 methods)"
]
},
"execution_count": 22,
@@ -954,7 +954,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Basically summing the above values together we have just done numerical integration over element area. By using these two simple concepts we are able to construct very interesting results.\n",
"Basically summing the above values together we have (almost) just done numerical integration over element area. By using these two simple concepts we are able to construct rest of the results.\n",
"\n",
"## Interpolation\n",
"\n",
@@ -1024,7 +1024,7 @@
{
"data": {
"text/plain": [
"JuliaFEM.ElementBasis(basis,dbasis)"
"JuliaFEM.Basis(basis,dbasis)"
]
},
"execution_count": 32,
@@ -1070,7 +1070,7 @@
"source": [
"### Interpolation in time domain\n",
"\n",
"To interpolate in time domain, call `DiscreteField` given time. Result is a `Increment` interpolated to that time. Here we interpolate the position of particle moving $x = \\frac{1}{2}t^2$ at time $t=1.0$."
"To interpolate in time domain, call `DiscreteField` given time. Result is an `Increment` interpolated to that time. Here we interpolate the position of particle moving $y = \\frac{1}{2}t^2$ at time $t=1.0$."
]
},
{
@@ -1093,8 +1093,8 @@
],
"source": [
"t = linspace(0, 2, 5)\n",
"x = 1/2*t.^2\n",
"t, x"
"y = 1/2*t.^2\n",
"t, y"
]
},
{
@@ -1118,8 +1118,8 @@
],
"source": [
"timesteps = TimeStep[]\n",
"for (ti, xi) in zip(t, x)\n",
" increment = Increment(xi)\n",
"for (ti, yi) in zip(t, y)\n",
" increment = Increment(yi)\n",
" push!(timesteps, TimeStep(ti, increment))\n",
"end\n",
"\n",
@@ -1132,7 +1132,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"It's also possible to take time derivatives. To do so, call `Field` with `TemporalBasis`, time, and additional argument `Val{:derivative}`. Again, same example:"
"It's also possible to take time derivatives. To do so, call `Field` with time and additional argument `Val{:derivative}`. Again, same example:"
]
},
{
@@ -1169,9 +1169,23 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"To interpolate in spatial domain, call `Increment` with `SpatialBasis` and coordinate $\\boldsymbol\\xi$. Increments to interpolate are the latest ones in each time step. Result depends from the content of the field. If it is scalar field, result will be scalar, if it's vector the result will be vector and so on.\n",
"To interpolate in spatial domain, call `Basis` with `Increment` and coordinate $\\boldsymbol\\xi$. Increments to interpolate are the latest ones in each time step. Result depends from the content of the field. If it is scalar field, result will be scalar, if it's vector the result will be vector and so on.\n",
"\n",
"Let's have a $\\left[0,1\\right]\\times\\left[0,1\\right] \\in \\mathbb{R}^2$ domain and $u_5 = 0.25$ displacement in upper right corner pointint to the $x_1$ direction. We seek for a center point of this at time $t=1.0$."
"Let's have a $\\Omega = \\left[0,1\\right]\\times\\left[0,1\\right] \\in \\mathbb{R}^2$ domain with a displacement field\n",
"\\begin{equation}\n",
"\\mathbf{u}\\left(X_1, X_2\\right)=t\\begin{bmatrix}X_{1}\\left(X_{2}+1\\right)\\\\\n",
"X_{1}\\left(4X_{2}-1\\right)\n",
"\\end{bmatrix}.\n",
"\\end{equation}\n",
"\n",
"Gradient is\n",
"\\begin{equation}\n",
"\\mbox{grad}\\left(\\mathbf{u}\\right)=t\\begin{bmatrix}X_{2}+1 & X_{1}\\\\\n",
"4X_{2}-1 & 4X_{1}\n",
"\\end{bmatrix}\n",
"\\end{equation}\n",
"\n",
"We look for a displacement in center point of the domain at time $t=1.0$:"
]
},
{
@@ -1182,12 +1196,10 @@
},
"outputs": [],
"source": [
"geometry = Field(Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
"displacement = Field(\n",
" (0.0, Vector[[0.0, 0.0], [0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]),\n",
" (1.0, Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]))\n",
"X = Basis(basis, dbasis, geometry)\n",
"u = Basis(basis, dbasis, displacement);"
"X = Field(Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])\n",
"u = Field(\n",
" (0.5, Vector{Float64}[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),\n",
" (1.5, Vector{Float64}[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]));"
]
},
{
@@ -1200,9 +1212,7 @@
{
"data": {
"text/plain": [
"2-element Array{Float64,1}:\n",
" 0.5625\n",
" 0.5 "
"([[0.0,0.0],[1.0,0.0],[1.0,1.0],[0.0,1.0]],Any[[0.0,0.0],[1.0,-1.0],[2.0,3.0],[0.0,0.0]])"
]
},
"execution_count": 38,
@@ -1211,15 +1221,10 @@
}
],
"source": [
"x2(xi, t) = X(xi, t) + u(xi, t)\n",
"x2([0.0, 0.0], 1.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating gradient of vector field, i.e, $u_{i,j}$:"
"# evaluate fields in time t=1.0 -> Increments\n",
"X_increment = X(1.0)\n",
"u_increment = u(1.0)\n",
"X_increment, u_increment"
]
},
{
@@ -1232,9 +1237,9 @@
{
"data": {
"text/plain": [
"2x2 Array{Float64,2}:\n",
" 0.125 0.125\n",
" 0.0 0.0 "
"2-element Array{Float64,1}:\n",
" 1.25\n",
" 1.0 "
]
},
"execution_count": 39,
@@ -1243,11 +1248,72 @@
}
],
"source": [
"# this needs some redesign.\n",
"N = FEM.Basis(basis, dbasis)\n",
"dN = FEM.ElementGradientBasis(N, geometry)\n",
"gradu = FEM.ElementFieldGradientBasis(dN, displacement)\n",
"gradu([0.0, 0.0], 1.0)"
"N(X_increment, [0.0, 0.0]) + N(u_increment, [0.0, 0.0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating gradient of vector field, i.e, $u_{i,j} = \\frac{\\partial u_i}{\\partial X_j}$:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x2 Array{Float64,2}:\n",
" 1.5 0.5\n",
" 1.0 2.0"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"N = Basis(basis, dbasis)\n",
"gradu = N(X_increment, u_increment, [0.0, 0.0], Val{:gradient})"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating time derivative of small strain tensor $\\epsilon$ is basically:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"2x2 Array{Float64,2}:\n",
" 1.5 0.75\n",
" 0.75 2.0 "
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"du = u(1.0, Val{:derivative})\n",
"grad_du = N(X_increment, du, [0.0, 0.0], Val{:gradient})\n",
"strain_rate = 1/2*(grad_du + grad_du')"
]
}
],
src/basis.jl
+109 -84
View File
@@ -1,53 +1,44 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
abstract Basis <: ContinuousField
### ELEMENT BASIS
""" This is the normal "user defined" basis functions familiar from school books. """
type ElementBasis <: Basis
type Basis <: ContinuousField
basis :: Function
dbasisdxi :: Function
end
function Basis(basis::Function, dbasisdxi::Function)
return ElementBasis(basis, dbasisdxi)
end
function Base.call(basis::ElementBasis, xi::Vector, time::Number=0.0)
""" Evaluate basis. """
function Base.call(basis::Basis, xi::Vector, time::Number=0.0)
basis.basis(xi) # passing time does not make much sense actually for this...
end
""" Interpolate increment in spatial domain using ElementBasis. """
function Base.call(basis::ElementBasis, increment::Increment, xi::Vector)
""" Evaluate gradient of basis. This need geometry information to calculate Jacobian. """
function Base.call(basis::Basis, geometry::Increment, xi::Vector,
::Type{Val{:gradient}})
dbasis = basis.dbasisdxi(xi)
J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
grad = inv(J)*dbasis
return grad
end
### INTERPOLATION IN SPATIAL DOMAIN ###
""" Interpolate increment in spatial domain using Basis. """
function Base.call(basis::Basis, increment::Increment, xi::Vector)
basis = basis.basis(xi)
sum([basis[i]*increment[i] for i=1:length(increment)])
end
### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD
""" 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
""" Return gradient of increment in spatial domain using Basis.. """
function Base.call(basis::Basis, geometry::Increment, field::Increment,
xi::Vector, ::Type{Val{:gradient}})
grad = basis(geometry, xi, Val{:gradient})
gradf = sum([grad[:,i]*field[i]' for i=1:length(field)])'
return gradf
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
### INTERPOLATION IN TIME DOMAIN ###
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. """
""" Interpolate discrete field in time domain. Return Increment. """
function Base.call(field::DiscreteField, time::Number,
time_extrapolation::Symbol=:linear,
time_interpolation::Symbol=:linear)
@@ -139,62 +130,23 @@ function Base.call(field::DiscreteField, time::Number,
end
### ELEMENT GRADIENT BASIS = ELEMENT BASIS + GEOMETRY
""" Interpolate time derivative of field in some time t. This assumes linear
interpolation in time which is then differentiated.
""" Gradient of ElementBasis, needs geometry information. """
type ElementGradientBasis <: Basis
element_basis :: ElementBasis
geometry :: DiscreteField
time_extrapolation :: Symbol
time_interpolation :: Symbol
end
Parameters
----------
field
Discrete field to interpolate. Must have timesteps and increments defined
time
Time to interpolate.
derivative
set Val{:derivative} to activate this function
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)
"""
function Base.call(field::DiscreteField, time::Number, ::Type{Val{:derivative}})
# 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")
@@ -237,3 +189,76 @@ function Base.call(field::DiscreteField, time::Number,
end
### ELEMENT FIELD BASIS = ELEMENT BASIS + FIELD
#=
""" 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
=#
### 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 grad(N::ElementBasis, f::ElementFieldBasis, X::ElementFieldBasis)
f.time_extrapolation == X.time_extrapolation || error("interpolation mismatch")
f.time_interpolation == X.time_interpolation || error("interpolation mismatch")
dN = ElementGradientBasis(N, X.field, f.time_extrapolation, f.time_interpolation)
dfdX = ElementFieldGradientBasis(dN, f.field, f.time_extrapolation, f.time_interpolation)
return dfdX
end
function grad(N::ElementBasis, f::DiscreteField, X::DiscreteField)
dN = ElementGradientBasis(N, X)
dfdX = ElementFieldGradientBasis(dN, f)
end
function ElementGradientBasis(element_basis::ElementBasis, geometry::DiscreteField)
return ElementGradientBasis(element_basis, geometry, :linear, :linear)
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
=#
### INTERPOLATION IN TIME DOMAIN ###
src/types.jl
+1 -1
View File
@@ -23,6 +23,6 @@ function IntegrationPoint(xi, weight)
IntegrationPoint(xi, weight, Dict())
end
call(N::ElementBasis, ip::IntegrationPoint) = N(ip.xi)
call(N::Basis, ip::IntegrationPoint) = N(ip.xi)
test/test_basis.jl
+135 -156
View File
@@ -6,7 +6,7 @@ module BasisTests
using JuliaFEM.Test
using JuliaFEM
using JuliaFEM: Basis, ElementGradientBasis, ElementFieldGradientBasis, Field
using JuliaFEM: Basis, Field
using JuliaFEM: Increment, TimeStep
function get_basis()
@@ -21,186 +21,103 @@ function get_basis()
-(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
return Basis(basis, dbasis)
end
### Test interpolation in spatial domain
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]
function test_basis_interpolation()
N = get_basis()
@test N([0.0, 0.0]) == 1/4*[1 1 1 1]
@test N([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)
# when going to \pm infinity, return the last one.
@test b([0.0, 0.0], -Inf) == T([0.5, 0.5], 0.0)
@test b([0.0, 0.0], +Inf) == T([0.5, 0.5], 1.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)
function test_basis_gradient_interpolation()
X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
# 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]
N = get_basis()
gradN = N(X, [0.0, 0.0], Val{:gradient})
@test gradN == 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()
function test_interpolation_of_scalar_increment_in_spatial_domain()
# in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
T = Increment([1.0, 2.0, 3.0, 4.0])
N = get_basis()
T_interpolated = N(T, [0.0, 0.0])
@test T_interpolated == T_known([0.5, 0.5])
end
function test_interpolation_of_gradient_of_scalar_increment_in_spatial_domain()
# 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])
X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
T = Increment([1.0, 2.0, 3.0, 4.0])
N = get_basis()
gradT = N(X, T, [0.0, 0.0], Val{:gradient})
gradT_expected(X) = [1-2*X[2] 3-2*X[1]]
@test gradT == gradT_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)
geometry = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
displacement = Increment(Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]])
N = get_basis()
X = N(geometry, [0.0, 0.0])
u = N(displacement, [0.0, 0.0])
x = X+u
u_expected(X) = [1/4*X[1]*X[2], 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))
@test isapprox(u, u_expected([0.5, 0.5]))
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]]))
# in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
# => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
X = Increment([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.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]
# displacement = Field(
# (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
# (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
@test isapprox(dU([0.0, 0.0], 1.0), dU_expected([0.5, 0.5], 1.0))
u = Increment([0.0 0.0; 1.0 -1.0; 2.0 3.0; 0.0 0.0]')
N = get_basis()
gradu(xi) = N(X, u, xi, Val{:gradient})
gradu_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
@test isapprox(gradu([0.0, 0.0]), gradu_expected([0.5, 0.5]))
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]]))
### Test interpolation in time domain
# wanted
#u = get_basis(element, "displacement")
#L = grad(diff(u))
#D = 1/2*(L + L')
#@test isapprox(D([0.0, 0.0], 1.0), ...)
basis, dbasis = get_basis()
N = Basis(basis, dbasis)
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)
function test_linear_time_extrapolation_of_field()
#T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
T = Field(
(0.0, [0.0, 0.0, 0.0, 0.0]),
(1.0, [1.0, 2.0, 3.0, 4.0]))
@test T(-1.0) == -1.0*[1.0, 2.0, 3.0, 4.0]
@test T( 3.0) == 3.0*[1.0, 2.0, 3.0, 4.0]
# when going to \pm infinity, return the last one.
@test T(-Inf) == 0.0*[1.0, 2.0, 3.0, 4.0]
@test T(+Inf) == 1.0*[1.0, 2.0, 3.0, 4.0]
end
#=
"""basic continuum interpolations"""
function test_basic_interpolations()
element = Quad4([1, 2, 3, 4])
element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0])
element["displacement"] = (
Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]],
Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
# from my old home works
basis = get_basis(element)
dbasis = grad(basis)
@test isapprox(basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0), [9/16, 1/2])
gradu = dbasis("displacement", [0.0, 0.0], 1.0)
epsilon = 1/2*(gradu + gradu')
rotation = 1/2*(gradu - gradu')
X = basis("geometry", [0.0, 0.0], 1.0)
k = 0.25
epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
@test isapprox(epsilon, epsilon_wanted)
@test isapprox(rotation, rotation_wanted)
F = I + gradu
@test isapprox(F, [X[2]*k+1 X[1]*k; 0 1])
C = F'*F
@test isapprox(C, [(X[2]*k+1)^2 (X[2]*k+1)*X[1]*k; (X[2]*k+1)*X[1]*k X[1]^2*k^2+1])
E = 1/2*(F'*F - I)
@test isapprox(E, [1/2*(X[2]*k + 1)^2-1/2 1/2*(X[2]*k+1)*X[1]*k; 1/2*(X[2]*k + 1)*X[1]*k 1/2*X[1]^2*k^2])
U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I)
@test isapprox(U, [1.24235 0.13804; 0.13804 1.02149])
function test_constant_time_extrapolation_of_field()
#T_known(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
T = Field(
(0.0, [0.0, 0.0, 0.0, 0.0]),
(1.0, [1.0, 2.0, 3.0, 4.0]))
@test T(-1.0, :constant) == [0.0, 0.0, 0.0, 0.0]
@test T( 3.0, :constant) == [1.0, 2.0, 3.0, 4.0]
end
=#
function test_time_extrapolation_of_field_with_single_timestep()
T = Field([1.0, 2.0, 3.0, 4.0])
@test T(1.0) == [1.0, 2.0, 3.0, 4.0]
end
function test_interpolation_in_temporal_basis()
i1 = Increment(0.0)
@@ -210,9 +127,6 @@ function test_interpolation_in_temporal_basis()
t2 = TimeStep(2.0, Increment[i2])
t3 = TimeStep(4.0, Increment[i3])
field = Field(TimeStep[t1, t2, t3])
info("field(1.0) = $(field(1.0))")
@test field(-Inf) == [0.0]
@test field( 0.0) == [0.0]
@test field( 1.0) == [0.5]
@@ -274,4 +188,69 @@ function test_derivative_interpolation_in_temporal_basis_in_variable_velocity_ch
@test isa(velocity, Increment) == true
end
#=
function test_time_derivative_gradient_interpolation_of_field()
# in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
# => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
# => d(u_i,j)/dt = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
geometry = Field([0.0 0.0; 1.0 0.0; 1.0 1.0; 0.0 1.0]')
displacement = Field(
(0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
(1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))
# wanted
#u = get_basis(element, "displacement")
#L = grad(diff(u))
#D = 1/2*(L + L')
#@test isapprox(D([0.0, 0.0], 1.0), ...)
basis, dbasis = get_basis()
N = Basis(basis, dbasis)
xi = [0.0, 0.0]
time = 1.2
grad = ElementGradientBasis(N, geometry)(xi, time)
increment = displacement(time, Val{:derivative})
diffgradu = sum([grad[:,i]*increment[i]' for i=1:length(increment)])'
diffgradu_expected(X, t) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
@test diffgradu == diffgradu_expected([0.5, 0.5], 1.2)
end
"""basic continuum interpolations"""
function test_basic_interpolations()
element = Quad4([1, 2, 3, 4])
element["geometry"] = Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]]
element["temperature"] = ([0.0, 0.0, 0.0, 0.0], [1.0, 2.0, 3.0, 4.0])
element["displacement"] = (
Vector[[0.0, 0.0], [0.0, 0.0], [0.00, 0.0], [0.0, 0.0]],
Vector[[0.0, 0.0], [0.0, 0.0], [0.25, 0.0], [0.0, 0.0]])
# from my old home works
basis = get_basis(element)
dbasis = grad(basis)
@test isapprox(basis("geometry", [0.0, 0.0], 1.0) + basis("displacement", [0.0, 0.0], 1.0), [9/16, 1/2])
gradu = dbasis("displacement", [0.0, 0.0], 1.0)
epsilon = 1/2*(gradu + gradu')
rotation = 1/2*(gradu - gradu')
X = basis("geometry", [0.0, 0.0], 1.0)
k = 0.25
epsilon_wanted = [X[2]*k 1/2*X[1]*k; 1/2*X[1]*k 0]
rotation_wanted = [0 k/2*X[1]; -k/2*X[1] 0]
@test isapprox(epsilon, epsilon_wanted)
@test isapprox(rotation, rotation_wanted)
F = I + gradu
@test isapprox(F, [X[2]*k+1 X[1]*k; 0 1])
C = F'*F
@test isapprox(C, [(X[2]*k+1)^2 (X[2]*k+1)*X[1]*k; (X[2]*k+1)*X[1]*k X[1]^2*k^2+1])
E = 1/2*(F'*F - I)
@test isapprox(E, [1/2*(X[2]*k + 1)^2-1/2 1/2*(X[2]*k+1)*X[1]*k; 1/2*(X[2]*k + 1)*X[1]*k 1/2*X[1]^2*k^2])
U = 1/sqrt(trace(C) + 2*sqrt(det(C)))*(C + sqrt(det(C))*I)
@test isapprox(U, [1.24235 0.13804; 0.13804 1.02149])
end
=#
end