mirror of
https://github.com/JuliaFEM/JuliaFEM.jl.git
synced 2026-08-06 04:21:33 +00:00
symbolic fields
This commit is contained in:
Jukka Aho
@@ -25,7 +25,10 @@
|
||||
"- Complete rewrite. The main ideas proposed in earlier version didn't work.\n",
|
||||
"\n",
|
||||
"### 2015-10-29\n",
|
||||
"- Third iteration."
|
||||
"- Third iteration.\n",
|
||||
"\n",
|
||||
"### 2015-11-03\n",
|
||||
"- Starts to be ready. Introduced symbolic fields."
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -516,7 +519,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Creating continuous and discrete fields\n",
|
||||
"### Creating continuous (and discrete) fields\n",
|
||||
"\n",
|
||||
"In previous section the concept of discrete fields was demonstrated. `DefaultDiscreteField` is subtype of `DiscreteField` which is subtype of `Field`:"
|
||||
]
|
||||
@@ -548,7 +551,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"There exists another type of fields too, namely `ContinuousField`s. Like the name already suggests it stores continuous time and spatial domain and it can be used to write custom fields. It needs to be callable. In this example `ContinuousField` is created which returns 1x4 dimensional array defined in $\\boldsymbol\\xi \\in [-1,1]^2, t\\in[0,1]$:"
|
||||
"There exists another type of fields too: `ContinuousField`s. Like the name already suggests it stores continuous time and spatial domain and it can be used to write custom fields. It needs to be callable. In this example `ContinuousField` is created which returns 1x4 dimensional array defined in $\\boldsymbol\\xi \\in [-1,1]^2, t\\in[0,1]$:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -707,7 +710,7 @@
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"call (generic function with 1258 methods)"
|
||||
"call (generic function with 1264 methods)"
|
||||
]
|
||||
},
|
||||
"execution_count": 22,
|
||||
@@ -731,7 +734,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Now we create two fields, one is discrete and another is continuous which takes discrete field as parameter:"
|
||||
"Now we create two fields, one is discrete and another is continuous taking discrete field as input argument:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -796,7 +799,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Next we add another discrete field and see that continuous field which is \"connected\" to discrete fields updates accordingly:"
|
||||
"Next we add another discrete field and see that continuous field \"connecting\" to discrete field updates accordingly:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -954,13 +957,13 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"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",
|
||||
"Basically summing the above values together we have (almost) 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",
|
||||
"In earlier the concepts of `DiscreteField` and `ContinuousField` were introduced, so that now we can define discrete set of values and continuous functions. It has also been shown how fields can depend from each other such a way that we interpolate continuous field from discrete field and vice versa. \n",
|
||||
"In earlier discussion the concepts of `DiscreteField` and `ContinuousField` were introduced, so that now we can define discrete set of values and continuous functions. It has also been shown how fields can depend from each other such a way that we can interpolate continuous field from discrete field and vice versa.\n",
|
||||
"\n",
|
||||
"This motivates us to create continuous fields which are interpolated from discrete fields using some polynomial basis. By thinking this way interpolation is nothing more than just an application of the earlier results already shown.\n",
|
||||
"This motivates us to create continuous fields that are interpolated from discrete fields using some polynomial basis. By thinking this way interpolation is nothing more than just an application of the earlier results already shown.\n",
|
||||
"\n",
|
||||
"Some fields has special meaning in JuliaFEM. Typically one needs to define at least discrete field `geometry` and continuous field `basis` so that basic interpolation is working on element. \n",
|
||||
"\n",
|
||||
@@ -1132,7 +1135,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"It's also possible to take time derivatives. To do so, call `Field` with time and additional argument `Val{:derivative}`. Again, same example:"
|
||||
"It's also possible to take time derivatives using finite difference approximation. To do so, call `Field` with time and additional argument `Val{:diff}`. Again, same example:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1155,7 +1158,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"pos(1.0, Val{:derivative})"
|
||||
"velocity = pos(1.0, Val{:diff})"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1169,7 +1172,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"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",
|
||||
"To interpolate in spatial domain, call `Basis` with `Increment` and coordinate $\\boldsymbol\\xi$. Increments to interpolate are the latest ones in each time step by default. 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 $\\Omega = \\left[0,1\\right]\\times\\left[0,1\\right] \\in \\mathbb{R}^2$ domain with a displacement field\n",
|
||||
"\\begin{equation}\n",
|
||||
@@ -1248,6 +1251,7 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# evaluate geometry + displacement at midpoint, basically x = X + u\n",
|
||||
"N(X_increment, [0.0, 0.0]) + N(u_increment, [0.0, 0.0])"
|
||||
]
|
||||
},
|
||||
@@ -1255,7 +1259,7 @@
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Calculating gradient of vector field, i.e, $u_{i,j} = \\frac{\\partial u_i}{\\partial X_j}$:"
|
||||
"Calculating gradients is also possible, for scalar or vector fields, if passing additional argument `Val{:grad}`. Here's an example of gradient of vector field, i.e, $u_{i,j} = \\frac{\\partial u_i}{\\partial X_j}$:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1280,14 +1284,14 @@
|
||||
],
|
||||
"source": [
|
||||
"N = Basis(basis, dbasis)\n",
|
||||
"gradu = N(X_increment, u_increment, [0.0, 0.0], Val{:gradient})"
|
||||
"gradu = N(X_increment, u_increment, [0.0, 0.0], Val{:grad})"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Calculating time derivative of small strain tensor $\\epsilon$ is basically:"
|
||||
"One can of course combine the above results. In this case field must be first interpolated in time domain and after that in spatial domain. So calculating time derivative of small strain tensor $\\epsilon$ is basically:"
|
||||
]
|
||||
},
|
||||
{
|
||||
@@ -1311,10 +1315,211 @@
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"du = u(1.0, Val{:derivative})\n",
|
||||
"grad_du = N(X_increment, du, [0.0, 0.0], Val{:gradient})\n",
|
||||
"du = u(1.0, Val{:diff})\n",
|
||||
"grad_du = N(X_increment, du, [0.0, 0.0], Val{:grad})\n",
|
||||
"strain_rate = 1/2*(grad_du + grad_du')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Symbolic fields\n",
|
||||
"\n",
|
||||
"There's still one kind of field not yet introduced, namely *symbolic fields*:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 42,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"JuliaFEM.SymbolicField(\"displacement\")"
|
||||
]
|
||||
},
|
||||
"execution_count": 42,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"u = Field(\"displacement\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"The idea is to write equation in expression form using tiny symbolic interface and evaluate equations afterwards. Main benefits: \n",
|
||||
"\n",
|
||||
"- possibility to simplify equations symbolically\n",
|
||||
"- prettier syntax\n",
|
||||
"- \"lazy evaluation\" gives possibilities to optimize performance.\n",
|
||||
"\n",
|
||||
"Here's a preliminary example how this could work:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 43,
|
||||
"metadata": {
|
||||
"collapsed": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# in future this will be something like this:\n",
|
||||
"# u = Field(\"displacement\")\n",
|
||||
"# eps = 1/2*(grad(u) + grad(u)')\n",
|
||||
"# strain_rate = diff(eps)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 44,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stdout",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"(1 / 2) * (grad(diff(displacement)) + grad(diff(displacement))')"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"# now this is\n",
|
||||
"expr = JuliaFEM.Expression(\"1/2*(grad(diff(displacement)) + grad(diff(displacement))')\")\n",
|
||||
"print(expr.expr)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Evaluation of expression: call `Basis` with `FieldSet`, $\\boldsymbol\\xi$ and time as arguments:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 45,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"INFO: expression = (1 / 2) * (grad(diff(displacement)) + grad(diff(displacement))')\n",
|
||||
"INFO: expression = 1 / 2\n",
|
||||
"INFO: new expression: 1 / 2\n",
|
||||
"INFO: expression = grad(diff(displacement)) + grad(diff(displacement))'\n",
|
||||
"INFO: expression = grad(diff(displacement))\n",
|
||||
"INFO: inside gradient operator\n",
|
||||
"INFO: gradient args: Any[:grad,:(diff(displacement))]\n",
|
||||
"INFO: expression inside gradient\n",
|
||||
"INFO: grad: evaluate field diff(displacement)\n",
|
||||
"INFO: taking time derivative of displacement\n",
|
||||
"INFO: expression = diff(displacement)\n",
|
||||
"INFO: inside diff operator\n",
|
||||
"INFO: data = Any[[0.0,0.0],[1.0,-1.0],[2.0,3.0],[0.0,0.0]]\n",
|
||||
"INFO: evaluate geometry\n",
|
||||
"INFO: evaluate gradient\n",
|
||||
"INFO: expression = grad(diff(displacement))'\n",
|
||||
"INFO: transpose\n",
|
||||
"INFO: expression = transpose(grad(diff(displacement)))\n",
|
||||
"INFO: expression = grad(diff(displacement))\n",
|
||||
"INFO: inside gradient operator\n",
|
||||
"INFO: gradient args: Any[:grad,:(diff(displacement))]\n",
|
||||
"INFO: expression inside gradient\n",
|
||||
"INFO: evaluate gradient\n",
|
||||
"INFO: new expression: transpose([1.5 0.5\n",
|
||||
" 1.0 2.0])\n",
|
||||
"INFO: new expression: [1.5 0.5\n",
|
||||
" 1.0 2.0] + transpose([1.5 0.5\n",
|
||||
" 1.0 2.0])\n",
|
||||
"INFO: new expression: (1 / 2) * ([1.5 0.5\n",
|
||||
" 1.0 2.0] + transpose([1.5 0.5\n",
|
||||
" 1.0 2.0]))\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"X = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])\n",
|
||||
"u = Field(\n",
|
||||
" (0.5, Vector[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),\n",
|
||||
" (1.5, Vector[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]))\n",
|
||||
"fs = FieldSet()\n",
|
||||
"fs[\"geometry\"] = X\n",
|
||||
"fs[\"displacement\"] = u\n",
|
||||
"xi = [0.0, 0.0]\n",
|
||||
"time = 1.0\n",
|
||||
"result = N(expr, fs, xi, time);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 46,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
":((1 / 2) * (2x2 Array{Float64,2}:\n",
|
||||
" 1.5 0.5\n",
|
||||
" 1.0 2.0 + transpose(2x2 Array{Float64,2}:\n",
|
||||
" 1.5 0.5\n",
|
||||
" 1.0 2.0)))"
|
||||
]
|
||||
},
|
||||
"execution_count": 46,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"result"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"What happened is that fields were automatically interpolated and substituted to the equation. The result is `Expr` and can be evaluated like usual:"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 47,
|
||||
"metadata": {
|
||||
"collapsed": false
|
||||
},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"text/plain": [
|
||||
"2x2 Array{Float64,2}:\n",
|
||||
" 1.5 0.75\n",
|
||||
" 0.75 2.0 "
|
||||
]
|
||||
},
|
||||
"execution_count": 47,
|
||||
"metadata": {},
|
||||
"output_type": "execute_result"
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"strain_rate = eval(result)"
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
|
||||
+3
-2
@@ -50,9 +50,10 @@ 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("fields.jl") # fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
# fields, see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
include("fields.jl")
|
||||
include("basis.jl") # interpolation of discrete fields
|
||||
include("symbolic.jl") # a thin symbolic layer for fields
|
||||
include("types.jl") # type definitions
|
||||
|
||||
### ELEMENTS ###
|
||||
|
||||
+5
-5
@@ -13,7 +13,7 @@ end
|
||||
|
||||
""" Evaluate gradient of basis. This need geometry information to calculate Jacobian. """
|
||||
function Base.call(basis::Basis, geometry::Increment, xi::Vector,
|
||||
::Type{Val{:gradient}})
|
||||
::Type{Val{:grad}})
|
||||
dbasis = basis.dbasisdxi(xi)
|
||||
J = sum([dbasis[:,i]*geometry[i]' for i=1:length(geometry)])
|
||||
grad = inv(J)*dbasis
|
||||
@@ -30,8 +30,8 @@ end
|
||||
|
||||
""" 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})
|
||||
xi::Vector, ::Type{Val{:grad}})
|
||||
grad = basis(geometry, xi, Val{:grad})
|
||||
gradf = sum([grad[:,i]*field[i]' for i=1:length(field)])'
|
||||
return gradf
|
||||
end
|
||||
@@ -140,10 +140,10 @@ field
|
||||
time
|
||||
Time to interpolate.
|
||||
derivative
|
||||
set Val{:derivative} to activate this function
|
||||
set Val{:diff} to activate this function
|
||||
|
||||
"""
|
||||
function Base.call(field::DiscreteField, time::Number, ::Type{Val{:derivative}})
|
||||
function Base.call(field::DiscreteField, time::Number, ::Type{Val{:diff}})
|
||||
|
||||
# FieldSet -> Field -> TimeStep -> Increment -> data
|
||||
|
||||
|
||||
+116
@@ -0,0 +1,116 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
|
||||
immutable SymbolicField <: Field
|
||||
name :: ASCIIString
|
||||
end
|
||||
|
||||
immutable Expression
|
||||
expr :: Expr
|
||||
end
|
||||
|
||||
function Expression(expression::ASCIIString)
|
||||
return Expression(parse(expression))
|
||||
end
|
||||
|
||||
function Base.convert(::Type{Field}, name::ASCIIString)
|
||||
return SymbolicField(name)
|
||||
end
|
||||
|
||||
function Base.convert(::Type{Symbol}, field::SymbolicField)
|
||||
return Symbol(field.name)
|
||||
end
|
||||
|
||||
function Base.(:*)(n::Number, field::SymbolicField)
|
||||
expr = Expr(:call, :*, n, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function grad(field::SymbolicField)
|
||||
expr = Expr(:call, :grad, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function diff(field::SymbolicField)
|
||||
expr = Expr(:call, :diff, Symbol(field.name))
|
||||
return Expression(expr)
|
||||
end
|
||||
|
||||
function Base.call(basis::Basis, expr::Expression, fieldset::FieldSet,
|
||||
xi::Vector, time::Number)
|
||||
return replace(expr.expr, basis, fieldset, xi, time)
|
||||
end
|
||||
|
||||
# Unbelievable code. I don't know why or how this works.
|
||||
""" Replace symbolic fields with real arrays. """
|
||||
function Base.replace(expression::Expr, basis::Basis, fieldset::FieldSet,
|
||||
xi::Vector, time::Number, data=Dict())
|
||||
|
||||
info("expression = $expression")
|
||||
if expression.head == symbol("'")
|
||||
info("transpose")
|
||||
expr = Expr(:call, :transpose, expression.args...)
|
||||
return replace(expr, basis, fieldset, xi, time, data)
|
||||
end
|
||||
operator = expression.args[1]
|
||||
|
||||
if operator == :diff
|
||||
info("inside diff operator")
|
||||
if !haskey(data, expression)
|
||||
data[expression] = fieldset[string(expression.args[2])](time, Val{:diff})
|
||||
end
|
||||
info("data = $(data[expression])")
|
||||
#return basis(data[expression], xi)
|
||||
return data[expression]
|
||||
end
|
||||
|
||||
if operator == :grad
|
||||
info("inside gradient operator")
|
||||
info("gradient args: $(expression.args)")
|
||||
field_name = expression.args[2]
|
||||
if isa(field_name, Expr)
|
||||
info("expression inside gradient")
|
||||
end
|
||||
#if startswith(string(field_name), "diff")
|
||||
if !haskey(data, field_name)
|
||||
info("grad: evaluate field $field_name")
|
||||
if isa(field_name, Symbol)
|
||||
data[field_name] = fieldset[string(field_name)](time)
|
||||
elseif isa(field_name, Expr) && (field_name.args[1] == :diff)
|
||||
info("taking time derivative of $(field_name.args[2])")
|
||||
#data[field_name] = fieldset[string(field_name.args[2])](time, Val{:diff})
|
||||
data[field_name] = replace(field_name, basis, fieldset, xi, time, data)
|
||||
end
|
||||
end
|
||||
if !haskey(data, :geometry)
|
||||
info("evaluate geometry")
|
||||
data[:geometry] = fieldset["geometry"](time)
|
||||
end
|
||||
#return Expr(:call, basis, data[:geometry], data[field_name], xi, Val{:gradient})
|
||||
#return :(basis($(data[:geometry]), $(data[field_name]), $xi, Val{:gradient}))
|
||||
info("evaluate gradient")
|
||||
return basis(data[:geometry], data[field_name], xi, Val{:grad})
|
||||
end
|
||||
|
||||
for i in 2:length(expression.args)
|
||||
arg = expression.args[i]
|
||||
if isa(arg, Expr)
|
||||
expression.args[i] = replace(arg, basis, fieldset, xi, time, data)
|
||||
end
|
||||
if isa(arg, Symbol) && haskey(fieldset, string(arg))
|
||||
if !haskey(data, arg)
|
||||
info("evaluate field $arg")
|
||||
data[arg] = fieldset[string(arg)](time)
|
||||
end
|
||||
#expression.args[i] = Expr(:call, basis, data[arg], xi)
|
||||
#return :(basis($(data[arg]), $xi))
|
||||
expression.args[i] = basis(data[arg], xi)
|
||||
end
|
||||
end
|
||||
|
||||
info("new expression: $expression")
|
||||
return expression
|
||||
end
|
||||
|
||||
+13
-13
@@ -37,7 +37,7 @@ function test_basis_gradient_interpolation()
|
||||
# P(X) = [1.0, X[1], X[2], X[1]*X[2]]
|
||||
# basis2, dbasis2 = JuliaFEM.calculate_lagrange_basis(P, X)
|
||||
N = get_basis()
|
||||
gradN = N(X, [0.0, 0.0], Val{:gradient})
|
||||
gradN = N(X, [0.0, 0.0], Val{:grad})
|
||||
@test gradN == 1/2*[-1 1 1 -1; -1 -1 1 1]
|
||||
# @test dN([0.0, 0.0]) == dbasis2([0.5, 0.5])
|
||||
end
|
||||
@@ -56,7 +56,7 @@ function test_interpolation_of_gradient_of_scalar_increment_in_spatial_domain()
|
||||
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 = N(X, T, [0.0, 0.0], Val{:grad})
|
||||
gradT_expected(X) = [1-2*X[2] 3-2*X[1]]
|
||||
@test gradT == gradT_expected([0.5, 0.5])
|
||||
end
|
||||
@@ -86,7 +86,7 @@ function test_interpolation_of_gradient_of_vector_field()
|
||||
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(xi) = N(X, u, xi, Val{:grad})
|
||||
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
|
||||
@@ -144,13 +144,13 @@ function test_derivative_interpolation_in_temporal_basis_in_constant_velocity()
|
||||
t2 = TimeStep(2.0, Increment[i2])
|
||||
t3 = TimeStep(4.0, Increment[i3])
|
||||
field = Field(TimeStep[t1, t2, t3])
|
||||
@test field(+Inf, Val{:derivative}) == [0.5]
|
||||
@test field(-Inf, Val{:derivative}) == [0.5]
|
||||
@test field( 0.0, Val{:derivative}) == [0.5]
|
||||
@test field( 0.5, Val{:derivative}) == [0.5]
|
||||
@test field( 1.0, Val{:derivative}) == [0.5]
|
||||
@test field( 1.5, Val{:derivative}) == [0.5]
|
||||
@test field( 2.0, Val{:derivative}) == [0.5]
|
||||
@test field(+Inf, Val{:diff}) == [0.5]
|
||||
@test field(-Inf, Val{:diff}) == [0.5]
|
||||
@test field( 0.0, Val{:diff}) == [0.5]
|
||||
@test field( 0.5, Val{:diff}) == [0.5]
|
||||
@test field( 1.0, Val{:diff}) == [0.5]
|
||||
@test field( 1.5, Val{:diff}) == [0.5]
|
||||
@test field( 2.0, Val{:diff}) == [0.5]
|
||||
end
|
||||
|
||||
function test_derivative_interpolation_in_temporal_basis_in_variable_velocity()
|
||||
@@ -164,12 +164,12 @@ function test_derivative_interpolation_in_temporal_basis_in_variable_velocity()
|
||||
# => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
|
||||
pos = Field(timesteps)
|
||||
|
||||
velocity = pos(1.0, Val{:derivative})[1]
|
||||
velocity = pos(1.0, Val{:diff})[1]
|
||||
v1 = (0.500 - 0.125)/0.5
|
||||
v2 = (1.125 - 0.500)/0.5
|
||||
@test isapprox(velocity, mean([v1, v2])) # = 1.00
|
||||
|
||||
velocity = pos(2.0, Val{:derivative})[1]
|
||||
velocity = pos(2.0, Val{:diff})[1]
|
||||
@test isapprox(velocity, (2.0-1.125)/0.5) # = 1.75
|
||||
end
|
||||
|
||||
@@ -183,7 +183,7 @@ function test_derivative_interpolation_in_temporal_basis_in_variable_velocity_ch
|
||||
end
|
||||
# => ((0.0,0.0),(0.5,0.125),(1.0,0.5),(1.5,1.125),(2.0,2.0))
|
||||
pos = Field(timesteps)
|
||||
velocity = pos(1.0, Val{:derivative})
|
||||
velocity = pos(1.0, Val{:diff})
|
||||
# after interpolation, we are expecting to have same type where we started
|
||||
@test isa(velocity, Increment) == true
|
||||
end
|
||||
|
||||
@@ -0,0 +1,126 @@
|
||||
# This file is a part of JuliaFEM.
|
||||
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
|
||||
|
||||
# https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/notebooks/2015-06-14-data-structures.ipynb
|
||||
|
||||
module SymbolicFieldTests
|
||||
|
||||
using JuliaFEM
|
||||
using JuliaFEM: Basis, Field, FieldSet, Expression, diff, grad
|
||||
|
||||
using JuliaFEM.Test
|
||||
|
||||
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(basis, dbasis)
|
||||
end
|
||||
|
||||
function get_fieldset()
|
||||
X = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
|
||||
T = Field(
|
||||
(0.0, [0, 0, 0, 0]),
|
||||
(1.0, [1, 2, 3, 4]))
|
||||
u = 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]]))
|
||||
# FIXME: this is not working
|
||||
#fieldset = FieldSet("geometry" => X, "temperature" => T, "displacement" => u)
|
||||
fieldset = FieldSet()
|
||||
fieldset["geometry"] = X
|
||||
fieldset["temperature"] = T
|
||||
fieldset["displacement"] = u
|
||||
return fieldset
|
||||
end
|
||||
|
||||
function test_create_symbolic_field()
|
||||
f = Field("temperature")
|
||||
@test isa(f, Field)
|
||||
end
|
||||
|
||||
function test_evaluate_symbolic_field()
|
||||
T = Field("temperature")
|
||||
expr = Symbol(T)
|
||||
@test expr == :(temperature)
|
||||
end
|
||||
|
||||
function test_simple_math()
|
||||
T = Field("temperature")
|
||||
eq = 1/2*T
|
||||
@test eq.expr == :(0.5*temperature)
|
||||
end
|
||||
|
||||
function test_evaluate_expression()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = 1/2*T
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), 1/2*mean([1, 2, 3, 4]))
|
||||
end
|
||||
|
||||
function test_evaluate_gradient_of_scalar_field()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = grad(T)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [0.0 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_gradient_of_vector_field()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
expr = grad(u)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_strain_rate()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
#expr = grad(u)
|
||||
expr = Expression("1/2*(grad(diff(displacement)) + grad(diff(displacement))')")
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.75; 0.75 2.0])
|
||||
end
|
||||
|
||||
function test_evaluate_grad_diff()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
#expr = grad(u)
|
||||
expr = Expression("grad(diff(displacement))")
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), [1.5 0.5; 1.0 2.0])
|
||||
end
|
||||
|
||||
function test_grad_diff_simplification()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
u = Field("displacement")
|
||||
expr1 = diff(grad(u))
|
||||
expr2 = Expression("grad(diff(displacement))")
|
||||
info("expr1 = $expr1")
|
||||
info("expr2 = $expr2")
|
||||
@test expr1 == expr2
|
||||
end
|
||||
|
||||
function test_evaluate_time_derivative()
|
||||
basis = get_basis()
|
||||
fieldset = get_fieldset()
|
||||
T = Field("temperature")
|
||||
expr = diff(T)
|
||||
result = basis(expr, fieldset, [0.0, 0.0], 1.0)
|
||||
@test isapprox(eval(result), mean([1.0, 2.0, 3.0, 4.0]))
|
||||
end
|
||||
|
||||
end
|
||||
Reference in New Issue
Block a user