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 types defined

Browse Source
This commit is contained in:
Jukka Aho
2015-10-09 01:44:08 +03:00
parent 90c2a6f969
commit 91665f6de4
9 changed files with 787 additions and 454 deletions
Download Patch File Download Diff File Expand all files Collapse all files
notebooks/2015-08-29-developing-juliafem.ipynb
+539 -202
View File
File diff suppressed because it is too large Load Diff
src/JuliaFEM.jl
+3 -1
View File
@@ -7,7 +7,9 @@ using Lexicon
using Logging
@Logging.configure(level=DEBUG)
include("types.jl") # type definitions
include("types.jl") # type definitions
include("interpolate.jl") # interpolation routines
### ELEMENTS ###
include("elements.jl")
src/elements.jl
+44 -130
View File
@@ -8,9 +8,11 @@ Related notebooks
2015-08-29-developing-juliafem.ipynb
=#
using JuliaFEM: interpolate
using FactCheck
using ForwardDiff
abstract Element
#= ELEMENT DEFINITIONS
@@ -143,9 +145,10 @@ function test_element(eltype)
fld = Field(0.0, collect(1:n))
Logging.info("Creating new scalar field $fld")
Logging.info("Pushing field to element.")
new_field!(el, :field1)
push_field!(el, :field1, fld)
@fact el[:field1][1] --> fld
new_fieldset!(el, "field1")
add_field!(el, "field1", fld)
fieldset = get_fieldset(el, "field1")
@fact fieldset[1] --> fld
mid = zeros(dim)
try
@@ -164,8 +167,9 @@ function test_element(eltype)
end
Logging.info("Interpolating scalar field at $mid")
f(field, xi, t) = el(xi)*el[field](t)
i = f(:field1, mid, 0.0)
#f(field, xi, t) = el(xi)*el[field](t)
#i = f(:field1, mid, 0.0)
i = interpolate(el, "field1", mid, 0.0)
Logging.info("Value: $i")
Logging.info("Element $eltype passed tests.")
end
@@ -188,21 +192,22 @@ get_dbasisdxi(el::Element, xi::Vector) = el.basis.dbasisdxi(xi)
"""
Interpolate field on element.
"""
function interpolate(el::Element, field::Symbol, xi::Vector, t::Number)
get_basis(el, xi)*el[field](t)
end
function interpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
interpolate(el, Symbol(field), xi, t)
function interpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number)
fieldset = get_fieldset(el, symbol(field_name))
field = interpolate(fieldset, t)
basis = get_basis(el)
interpolate(basis, field, xi)
end
"""
Interpolate derivative of field on element.
"""
function dinterpolate(el::Element, field::Symbol, xi::Vector, t::Number)
get_dbasisdxi(el, xi)*el[field](t)
end
function dinterpolate(el::Element, field::ASCIIString, xi::Vector, t::Number)
dinterpolate(el, Symbol(field), xi, t)
function dinterpolate(el::Element, field_name::Union{Symbol, ASCIIString}, xi::Vector, t::Number)
#get_dbasisdxi(el, xi)*el[field](t)
fieldset = get_fieldset(el, symbol(field_name))
field = interpolate(fieldset, t)
basis = get_basis(el)
dinterpolate(basis, field, xi)
end
"""
@@ -210,155 +215,64 @@ Get jacobian of element evaluated at point ξ on element in reference configurat
Parameters
----------
el::Element
xi::Vector
geometry_field::Any, optional
time::Number, optional, default=0.0
el :: Element
xi :: Vector
geometry_field :: Any, optional
time :: Number
Returns
-------
Vector or Matrix
depending on element type
Notes
-----
Big "J" comes from reference (undeformed) configuration.
"""
function get_Jacobian(el::Element, xi, t, geometry_field=:Geometry)
function get_jacobian(el::Element, xi, t, geometry_field=symbol("geometry"))
dinterpolate(el, geometry_field, xi, t)
end
"""
Get jacobian of element evaluated at point ξ on element in current configuration.
Notes
-----
Small "j" comes from current (deformed) configuration.
"""
function get_jacobian(el::Element, xi, t, geometry_field=:Geometry, displacement_field=:displacement)
dbasisdxi = get_dbasisdxi(el, xi)
X = get_field(el, geometry_field)(t)
u = get_field(el, displacement_field)(t)
j = dbasisdxi*(X+u)
return j
end
"""
Evaluate partial derivatives of basis, dbasis/dX
"""
function get_dbasisdX(el::Element, xi, t)
dbasisdxi = get_dbasisdxi(el, xi)
J = get_Jacobian(el, xi, t)
J = get_jacobian(el, xi, t)
dbasisdxi*inv(J)
end
"""
Evaluate partial derivatives of basis, dbasis/dx
"""
function get_dbasisdx(el::Element, xi, t)
dbasisdxi = get_dbasisdxi(el, xi)
j = get_jacobian(el, xi, t)
dbasisdxi*inv(j)
""" Create new empty set of fields for element. """
function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString})
el.fields[symbol(field_name)] = FieldSet()
end
function new_fieldset!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field)
new_fieldset!(el, symbol(field_name))
add_field!(el, symbol(field_name), field)
end
""" Create new empty field of some type. """
function new_field!(el::Element, field_name::Symbol)
el.fields[field_name] = Field[]
end
function new_field!(el::Element, field_name::Symbol, field::Field)
new_field!(el, field_name)
push_field!(el, field_name, field)
end
function new_field!(el::Element, field_name::ASCIIString, field::Field)
new_field!(el, Symbol(field_name), field)
end
function new_field!(el::Element, field_name::ASCIIString)
new_field!(el, Symbol(field_name))
""" Add new field to fieldset of element. """
function add_field!(el::Element, field_name::Union{Symbol, ASCIIString}, field::Field)
push!(el.fields[symbol(field_name)], field)
end
""" Push to existing set field of fields. """
function push_field!(el::Element, field_name::Symbol, field::Field)
push!(el.fields[field_name], field)
""" Get fieldset. """
function get_fieldset(el::Element, field_name::Union{Symbol, ASCIIString})
el.fields[symbol(field_name)]
end
function push_field!(el::Element, field_name::ASCIIString, field::Field)
push_field!(el, Symbol(field_name), field)
""" Get fieldset, convenient function. """
function Base.getindex(el::Element, field_name::Union{Symbol, ASCIIString})
get_fieldset(el, field_name)
end
""" Get field variable. """
function get_field(el::Element, field_name::Symbol)
el.fields[field_name]
end
function get_field(el::Element, field_name::ASCIIString)
el.fields[Symbol(field_name)]
end
function Base.getindex(el::Element, field_name::Union{ASCIIString, Symbol})
get_field(el, field_name)
end
#=
"""
Evaluate some field in point ξ on element using basis functions.
Parameters
----------
el :: Element
field :: Any
xi :: Vector
Returns
-------
Scalar, Vector, Tensor, depending on what is type of field to interpolate.
Notes
-----
This has another version which returns multiple values for set of coordinates {ξᵢ}.
dinterpolate returns derivatives.
Examples
--------
>>> field = [1.0, 2.0, 3.0, 4.0]
>>> set_field(el, :temperature, field)
>>> interpolate(el, :temperature, [0.0, 0.0])
15.0
"""
function interpolate(el::Element, field, xi::Number)
interpolate(el, field, [xi])
end
function interpolate(el::Element, field, xi::Vector)
field = get_field(el, field)
sum(get_basis(el, xi) .* field)
end
function interpolate(el::Element, field, xis::Array{Vector, 1})
field = get_field(el, field)
interpolate_(xi) = sum(get_basis(el, xi) .* field)
map(interpolate_, xis)
end
function dinterpolate(el::Element, field, xi::Number)
dinterpolate(el, field, [xi])
end
function dinterpolate(el::Element, field, xi::Vector)
fld = get_field(el, field)
dbasis = get_dbasisdxi(el, xi)
if isa(dbasis, Vector)
return sum(dbasis .* fld)
end
return sum([fld[i]*dbasis[i,:] for i in 1:length(fld)])
end
=#
"""
calculate "local" normals in elements, in a way that
n = Nᵢnᵢ gives some reasonable results for ξ ∈ [-1, 1]
"""
function calculate_normals!(el::Element, t, field_name=:Normals)
function calculate_normals!(el::Element, t, field_name=symbol("normals"))
new_field!(el, field_name, Vector)
for xi in Vector[[-1.0], [1.0]]
t = dinterpolate(el, :Geometry, xi)
@@ -371,7 +285,7 @@ end
"""
Alter normal field such that normals of adjacent elements are averaged.
"""
function average_normals!(elements, normal_field=:Normals)
function average_normals!(elements, normal_field=symbol("normals"))
d = Dict()
for el in elements
c = get_connectivity(el)
src/equations.jl
+3 -3
View File
@@ -35,6 +35,8 @@ get_dbasisdx(eq::Equation, ip::IntegrationPoint) = get_dbasisdx(get_element(eq),
interpolate(eq::Equation, field::Union{ASCIIString, Symbol}, ip::IntegrationPoint) = interpolate(get_element(el), field, ip.xi)
integrate_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
integrate_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
get_lhs(eq::Equation, t::Number) = has_lhs(eq) ? integrate(eq, get_lhs, t) : nothing
get_rhs(eq::Equation, t::Number) = has_rhs(eq) ? integrate(eq, get_rhs, t) : nothing
"""
@@ -48,7 +50,7 @@ function get_detJ(el::Element, ip::IntegrationPoint, t::Float64)
get_detJ(el, ip.xi, t)
end
function get_detJ(el::Element, xi::Vector, t::Float64)
J = get_Jacobian(el, xi, t)
J = get_jacobian(el, xi, t)
s = size(J)
return s[1] == s[2] ? det(J) : norm(J)
end
@@ -79,5 +81,3 @@ function set_global_dofs!(eq::Equation, dofs)
eq.global_dofs = dofs
end
# Equations for heat problems
#include("heat_equations.jl")
src/interpolate.jl
+57
View File
@@ -0,0 +1,57 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM: Basis, Field, FieldSet, diff
"""
Interpolate field u using basis N in point xi.
"""
function interpolate{T}(N::Basis, u::Field{Vector{T}}, xi::Array{Float64,1})
N(xi)*u
end
"""
Interpolate field u using basis N in set of points xi. Convenient function.
"""
function interpolate{T}(N::Basis, u::Field{Vector{T}}, xis::Array{Array{Float64,1},1})
T[N(xi)*u for xi in xis]
end
function interpolate{T}(N::Basis, u::Field{T}, xi::Array{Float64,1})
u.values
end
"""
Interpolate a field from fieldset for some time t.
"""
function interpolate(fields::FieldSet, t::Number)
if length(fields) == 0
throw("Empty set of fields.")
end
if t <= fields[1].time
return Field(t, fields[1].values)
end
if t >= fields[end].time
return fields[end]
end
i = length(fields)
while fields[i].time >= t
i -= 1
end
if fields[i].time == t
return fields[i]
end
#Logging.debug("doing linear interpolation between fields $i and $(i+1)")
f1 = fields[i]
t1 = f1.time
f2 = fields[i+1]
t2 = f2.time
dt = t2 - t1
nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
f = Field(t, nw)
return f
end
function dinterpolate(N::Basis, u::Field, xi::Array{Float64, 1})
dN = diff(N)
dN(xi)*u
end
src/problems.jl
+3 -3
View File
@@ -10,9 +10,9 @@ get_equation(pr::Type{Problem}, el::Type{Element}) = nothing
"""
Add new element to problem
"""
function add_element!(pr::Problem, el::Element)
eq = get_equation(typeof(pr), typeof(el))
push!(pr.equations, eq(el))
function add_element!(problem::Problem, element::Element)
equation = get_equation(typeof(problem), typeof(element))
push!(problem.equations, equation(element))
end
"""
src/types.jl
+37 -69
View File
@@ -5,89 +5,55 @@
using ForwardDiff
""" Field. """
""" Field is a fundamental type which holds some values in some time t """
type Field{T}
time :: Float64
increment :: Int64
values :: T
end
""" Initialize field. """
function Field(time, values)
Field(time, 1, values)
end
""" Get length of a field (number of basis functions in practice). """
Base.length(f::Field) = length(f.values)
function Base.length(f::Field)
length(f.values)
end
""" Get field discrete value at point i. """
Base.getindex(f::Field, i::Int64) = f.values[i]
""" Interpolate field h(ξ)*f = x*f """
function interpolate{T}(x::Vector, f::Field{Vector{T}})
function Base.getindex(f::Field, i::Int64)
f.values[i]
end
""" Multiply field with some constant k. """
function Base.(:*)(k::Number, f::Field)
Field(f.time, k*f.values)
end
""" Multiply field with some vector x. """
function Base.(:*)(x::Vector, f::Field)
@assert length(x) == length(f)
sum([f[i]*x[i] for i in 1:length(f)])
end
function interpolate{T}(x::Matrix, f::Field{Vector{T}})
""" 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)])
end
function interpolate(x::Vector, f::Field)
f.values*x
end
Base.(:*)(x::Union{Vector, Matrix}, f::Field) = interpolate(x, f)
""" Interpolate field (h*f)(ξ) """
Base.(:*)(f::Function, fld::Field) = (x) -> f(x)*fld
""" Multiply field with some constant k. """
Base.(:*)(k::Float64, f::Field) = Field(f.time, k*f.values)
""" 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
""" Interpolate from set of fields x*[f1, f2] where x is evaluated basis. """
"""
FieldSet is array of fields, each field maybe having different time and/or increment.
"""
typealias FieldSet Array{Field, 1}
""" Multiply fieldset with some vector x. """
Base.(:*)(x::Array{Float64, 1}, f::Array{Field}) = sum(x .* f)
""" 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 a field from finite set of fields some time t ∈ R.
"""
function call(fields :: Array{Field, 1}, t::Number)
if length(fields) == 0
throw("Empty set of fields.")
end
if t <= fields[1].time
return Field(t, fields[1].values)
end
if t >= fields[end].time
return fields[end]
end
i = length(fields)
while fields[i].time >= t
i -= 1
end
if fields[i].time == t
return fields[i]
end
#Logging.debug("doing linear interpolation between fields $i and $(i+1)")
f1 = fields[i]
t1 = f1.time
f2 = fields[i+1]
t2 = f2.time
dt = t2 - t1
nw = (t2-t)/dt*f1.values + (t-t1)/dt*f2.values
f = Field(t, nw)
return f
""" Add new field to fieldset. """
function add_field!(fs::FieldSet, field::Field)
push!(fs, field)
end
function call(field::Field, t::Float64)
Field(t, field.increment, field.values)
end
""" Basis function. """
@@ -95,21 +61,23 @@ type Basis
basis :: Function
dbasisdxi :: Function
end
""" Constructor of basis function. """
function Basis(basis)
Basis(basis, ForwardDiff.jacobian(basis))
end
""" 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
""" Evaluate basis function in point ξ. """
call(b::Basis, xi) = b.basis(xi)
""" Get partial derivative of basis function. """
∂(h::Basis) = h.dbasisdxi
diff(h::Basis) = h.dbasisdxi
derivative(h::Basis) = h.dbasisdxi
# convenient functions
""" Evaluate basis function in point ξ. """
call(b::Basis, xi) = b.basis(xi)
#""" 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
test/test_elements.jl
+45 -3
View File
@@ -2,7 +2,49 @@
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using FactCheck
using JuliaFEM: test_element
using JuliaFEM: Element, Basis, FieldSet
# prototype element
type MockElement <: Element
connectivity :: Array{Int, 1}
basis :: Basis
fields :: Dict{Symbol, FieldSet}
end
function MockElement(connectivity)
h(xi) = [
(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4]
dh(xi) = [
-(1-xi[2])/4.0 -(1-xi[1])/4.0
(1-xi[2])/4.0 -(1+xi[1])/4.0
(1+xi[2])/4.0 (1+xi[1])/4.0
-(1+xi[2])/4.0 (1-xi[1])/4.0]
basis = Basis(h, dh)
MockElement(connectivity, basis, Dict())
end
JuliaFEM.get_number_of_basis_functions(el::Type{MockElement}) = 4
JuliaFEM.get_element_dimension(el::Type{MockElement}) = 2
using JuliaFEM: test_element
facts("test test_element against mock element") do
test_element(MockElement)
end
using JuliaFEM: new_fieldset!, add_field!, Field, get_fieldset
facts("test adding fieldsets and fields to element") do
el = MockElement([1, 2, 3, 4])
fieldset = new_fieldset!(el, "geometry")
field1 = Field(0.0, [0.0, 0.0, 0.0, 0.0])
add_field!(el, "geometry", field1)
field2 = Field(1.0, [1.0, 1.0, 1.0, 1.0])
add_field!(fieldset, field2)
fields = get_fieldset(el, "geometry")
@fact length(fields) --> 2
@fact fields[1] --> field1
@fact fields[2] --> field2
end
using JuliaFEM: Quad4
test_element(Quad4)
test/test_types.jl
+56 -43
View File
@@ -1,18 +1,11 @@
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM: Basis, Field, get_field, diff
using JuliaFEM: Basis, Field, FieldSet, interpolate, dinterpolate
using FactCheck
facts("test fields and interpolation") do
# simple interpolation in domain [-1, 1]
N = Basis((ξ) -> [0.5*(1.0-ξ[1]), 0.5*(1.0+ξ[1])])
u = Field(0.0, [0.0, 1.0])
@fact N([0.0])*u --> 0.5
@fact (N*u)([0.0]) --> 0.5
# multiply of field with constant
facts("test fields") do
# multiple field with some constant
u1 = Field(0.0, [0.0, 1.0])
u2 = 3.0*u1
@fact u1.time --> 0.0
@@ -24,48 +17,68 @@ facts("test fields and interpolation") do
u2 = Field(0.0, [1.0, 2.0])
u3 = u1 + u2
@fact u3.values --> [1.0, 3.0]
end
# interpolation between two fields in time domain
facts("test interpolation of fields") do
# interpolation of field in spatial domain
N = Basis((xi) -> [0.5*(1.0-xi[1]), 0.5*(1.0+xi[1])])
u = Field(0.0, [0.0, 1.0])
@fact interpolate(N, u, [0.0]) --> 0.5
# interpolation of fieldset in time domain
u1 = Field(0.0, [0.0, 1.0])
u2 = Field(0.0, [1.0, 2.0])
t = Basis((t) -> [1-t, t])
u = Field[u1, u2]
u2 = (t*u)(0.5)
@fact u2.values --> [0.5, 1.5]
u2 = Field(1.0, [1.0, 2.0])
u = FieldSet([u1, u2])
@fact interpolate(u, 0.5).values --> [0.5, 1.5]
@fact interpolate(u, 0.5).time --> 0.5
# interpolation in set of fields is defined for every time value
# interpolation of fieldset is defined for every time value:
u1 = Field(0.0, [0.0, 0.0])
u2 = Field(1.0, [1.0, 2.0])
u3 = Field(2.0, [0.5, 1.5])
u = Field[u1, u2, u3]
@fact u(-1.0).values --> [0.0, 0.0] # "out of range -" -> first known value
@fact u(0.0).values --> [0.0, 0.0]
@fact u(1.0).values --> [1.0, 2.0]
@fact u(2.0).values --> [0.5, 1.5]
@fact u(3.0).values --> [0.5, 1.5] # "out of range +" -> last known value
@fact u(0.5).values --> [0.5, 1.0]
@fact u(1.5).values --> [0.75, 1.75]
u = FieldSet([u1, u2, u3])
@fact interpolate(u, -1.0).values --> [0.0, 0.0] # "out of range -" -> first known value
@fact interpolate(u, 0.0).values --> [0.0, 0.0]
@fact interpolate(u, 1.0).values --> [1.0, 2.0]
@fact interpolate(u, 2.0).values --> [0.5, 1.5]
@fact interpolate(u, 3.0).values --> [0.5, 1.5] # "out of range +" -> last known value
@fact interpolate(u, 0.5).values --> [0.5, 1.0]
@fact interpolate(u, 1.5).values --> [0.75, 1.75]
# use Inf to get very first or last value of field
@fact u(-Inf).values --> [0.0, 0.0]
@fact u(+Inf).values --> [0.75, 1.75]
@fact interpolate(u, -Inf).values --> [0.0, 0.0]
@fact interpolate(u, +Inf).values --> [0.5, 1.5]
# multidimensional interpolation with and without derivatives
X = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
h = Basis((xi) ->
[(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4])
# midpoint of field
@fact (h*X)([0.0, 0.0]) --> [0.5, 0.5]
@fact h([0.0, 0.0])*X --> [0.5, 0.5]
# derivatives of field at midpoint
@fact diff(h)([0.0, 0.0])*X --> [0.5 0.0; 0.0 0.5]
@fact (diff(h)*X)([0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
h(xi) = [
(1-xi[1])*(1-xi[2])/4
(1+xi[1])*(1-xi[2])/4
(1+xi[1])*(1+xi[2])/4
(1-xi[1])*(1+xi[2])/4]
dh(xi) = [
-(1-xi[2])/4.0 -(1-xi[1])/4.0
(1-xi[2])/4.0 -(1+xi[1])/4.0
(1+xi[2])/4.0 (1+xi[1])/4.0
-(1+xi[2])/4.0 (1-xi[1])/4.0]
N = Basis(h, dh)
# multiplying scalar field with a vector -> vector
b = Basis((xi) -> [1/2*(1-xi[1]), 1/2*(1+xi[1])])
f = Field(0.0, 100.0)
@fact b(0.0) * f --> [50.0, 50.0]
X = Field(0.0, Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
# get midpoint of field in spatial domain
@fact interpolate(N, X, [0.0, 0.0]) --> [0.5, 0.5]
# derivatives of field at midpoint
@fact dinterpolate(N, X, [0.0, 0.0]) --> [0.5 0.0; 0.0 0.5]
# interpolate of scalar field -> scalar
H = Field(0.0, 6.0)
@fact interpolate(N, H, [0.0, 0.0]) --> 6.0
# multiplying scalar field with a vector -> vector
# this is actually not so good idea...
#h(xi) = [1/2*(1-xi[1]), 1/2*(1+xi[1])]
#dh(xi) = [-1/2 1/2]'
#N = Basis(h, dh)
#f = Field(0.0, 100.0)
#@fact interpolate(N, f, [0.0]) --> [50.0, 50.0]
end