Author(s): Jukka Aho
Abstract: Description of basis data structures. In this notebook the concepts of Increment, TimeStep, DiscreteField, ContinuousField, FieldSet, SpatialBasis and TemporalBasis are intoduced. By combining these atomic structures one is able to form finite elements and interpolate it's fields in time and spatial domain.
Increment is the most atomic structure. It's a vector-like object with 1 dimension. Each element in Increment can be scalar, vector or tensor (2 or 4 order). It's easy to extend Increment to have other data types too.TimeStep is container for increments in certain time $t$.DefaultDiscreteField is container for timesteps for a single field.FieldSet is container for all fields.DefaultContinuousField can have continuous field variables.Typical element structure so far:
type MyElement <: Element
connectivity :: Array{Int, 1} # describes how dofs of this element is connected to another elements in global level
basis :: Basis # describes how to interpolate fields
fields :: ???
end
Here's the strategy in short: each non-linear iteration is Increment, what is a vector-like object containing data. Elements of Increment can be scalars, vectors or tensors. Increment belongs to TimeStep. TimeStep contains one or more increments. Then we have Field which contains all timesteps. And finally we have FieldSet which contains all fields.
FieldSet -> Field -> TimeStep -> Increment -> data
for example,
FieldSet -> Field -> TimeStep -> Increment -> data
FieldSet -> "temperature" -> 0.0 -> 1 -> [1,2,3,4]
2 -> [2,3,4,5]
FieldSet -> "temperature" -> 1.0 -> 1 -> [3,4,5,6]
2 -> [5,6,7,8]
and so on. Here we create a FieldSet containing one field temperature which contains two TimeSteps, two Increments in both of them.
In [1]:
using JuliaFEM: Increment, TimeStep, Field, FieldSet, DefaultDiscreteField
fs = FieldSet()
i1 = Increment([1, 2, 3])
i2 = Increment([2, 3, 4])
t1 = TimeStep(1.0, Increment[i1, i2])
i3 = Increment([2, 3, 4])
i4 = Increment([3, 4, 5])
t2 = TimeStep(2.0, Increment[i3, i4])
f1 = DefaultDiscreteField(TimeStep[t1, t2])
fs["temperature"] = f1
fs
Out[1]:
Accessing last increment of last timestep of temperature can be done in the following way:
In [2]:
fs["temperature"][2][2]
Out[2]:
Because the model is deep and quite heavy to type some "shortcuts" are provided, but keep in mind that everything is there in place. Here's the shortened version:
In [3]:
fs = FieldSet() # create empty fieldset
fs["temperature"] = [1, 2, 3, 4] # create temperature field with time t=0
first(fs["temperature"]), last(fs["temperature"]) # pick first and last timestep of temperature
Out[3]:
This way we can easily create discrete scalar, vector and tensor fields.
In [4]:
fs2 = FieldSet()
fs2["constant scalar field"] = 1
fs2["scalar field"] = [1, 2, 3, 4]
fs2["vector field"] = reshape(collect(1:8), 2, 4)
fs2["second order tensor field"] = reshape(collect(1:3*3*4), 3, 3, 4)
fs2["fourth order tensor field"] = reshape(collect(1:3*3*3*3*4), 3, 3, 3, 3, 4)
fs2
Out[4]:
Or even more compactly:
In [5]:
FieldSet("geometry" => [1, 2, 3, 4], "temperature" => [0, 0, 0, 0], "density" => 7850)
Out[5]:
By default using this "fast typing" fields are defined at time $t=0.0$.
In [6]:
fs2["vector field"][end].time # pick last timestep of field "vector field"
Out[6]:
Creating new field with several time steps defined can also be done compactly. Each tuple has time and increment data.
In [7]:
t0 = TimeStep(0.0, Increment([1, 2, 3, 4]))
t1 = TimeStep(0.5, Increment([2, 3, 4, 5]))
t2 = TimeStep(1.0, Increment([1, 1, 1, 1]))
fs2["time series"] = [t0, t1, t2]
info("time on last timestep: $(fs2["time series"][end].time)")
Or even without explicitly expressing time. In that case time step size is 1 second starting from 0.
In [8]:
fs2["time series 2"] = [1, 2], [2, 3], [3, 4]
fs2
Adding another field, with different time.
In [9]:
T0 = first(fs["temperature"]) # pick first timestep (or to be spesific, last increment of first timestep)
T1 = Increment(T0 + 1) # create new increment from old one
timestep = TimeStep(1.0, T1) # create new timestep at t=1.0
push!(fs["temperature"], timestep) # push to field
Out[9]:
Normal stuff like dot product works:
In [10]:
S1 = Increment([1, 2, 3])
S2 = Increment([2, 3, 4])
dot(S1, S2), dot([1, 2, 3], S2), dot(S1, [2, 3, 4])
Out[10]:
In [11]:
1/2*(S1 + S2)
Out[11]:
In [12]:
dot([1, 2], Increment[S1, S2])
Out[12]:
Creating empty Increment:
In [13]:
f = zeros(Increment, 2, 4)
Out[13]:
Create similar increment with new data:
In [14]:
g = similar(f, ones(8))
Out[14]:
In [15]:
using JuliaFEM: DefaultDiscreteField, DiscreteField, Field
DefaultDiscreteField <: DiscreteField <: Field
Out[15]:
There exists another type of fields too: ContinuousFields. 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]$:
In [16]:
using JuliaFEM: FieldSet, DefaultContinuousField, ContinuousField
In [17]:
type MyContinuousField <: ContinuousField
end
function Base.call(f::MyContinuousField, xi::Vector, time::Number)
time/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])]
end
f = MyContinuousField()
Out[17]:
In [18]:
f([0.0, 0.0], 1.0)
Out[18]:
In [19]:
fs = FieldSet()
fs["basis"] = MyContinuousField()
fs
Out[19]:
In [20]:
fs["basis"]([0.0, 0.0], 1.0)
Out[20]:
Again we have a shortcut to quickly define continuous fields:
In [21]:
fs["continuous field"] = (xi, t) -> xi[1]*xi[2]*t
fs["continuous field"]([1.0, 2.0], 3.0)
Out[21]:
Naturally we can pass another fields or even fieldsets to continuous field to make fields depend from each other. Here's another example, where ContinuousField takes another field and operates it with some function.
In [22]:
type MyContinuousField2 <: ContinuousField
basis :: Function
discrete_field :: DiscreteField
end
function Base.call(field::MyContinuousField2, xi::Vector, time::Number=1.0)
data = last(field.discrete_field) # get the last timestep last increment
basis = field.basis(xi) # evaluate basis at point ξ.
sum([basis[i]*data[i] for i=1:length(data)]) # sum results
end
Out[22]:
Now we create two fields, one is discrete and another is continuous taking discrete field as input argument:
In [23]:
fs = FieldSet()
fs["discrete field"] = [1, 2, 3, 4]
basis(xi) = 1/4*[
(1-xi[1])*(1-xi[2]),
(1+xi[1])*(1-xi[2]),
(1+xi[1])*(1+xi[2]),
(1-xi[1])*(1+xi[2])]
fs["continuous field"] = MyContinuousField2(basis, fs["discrete field"])
Out[23]:
Results:
In [24]:
fs["continuous field"]([0.0, 0.0], 1.0), last(fs["discrete field"])
Out[24]:
Next we add another discrete field and see that continuous field "connecting" to discrete field updates accordingly:
In [25]:
T0 = last(fs["discrete field"])
push!(fs["discrete field"], TimeStep(1.0, Increment(T0 + 1.0))) # push to field
fs
Out[25]:
Updated results, notice how the value of continuous field changes according to the update of discrete field.
In [26]:
fs["continuous field"]([0.0, 0.0], 1.0), last(fs["discrete field"])
Out[26]:
What we just did is that we actually interpolated discrete field using continuous functions. We evaluated discrete field using bilinear basis at midpoint of "element":
In [27]:
1/4*(1+2+3+4), 1/4*(2+3+4+5)
Out[27]:
From continuous field we can also go back to discrete fields, if needed. Here we evaluate continuous field in four discrete points, let's call them to Gauss quadrature points.
In [28]:
using JuliaFEM: DiscreteField
type MyDiscreteField <: DiscreteField
discrete_points :: Vector
continuous_field :: ContinuousField
end
Base.length(field::MyDiscreteField) = length(field.discrete_points)
Base.endof(field::MyDiscreteField) = endof(field.discrete_points)
Base.last(field::MyDiscreteField) = Float64[field[i] for i=1:length(field)]
function Base.getindex(field::MyDiscreteField, idx::Int64)
field.continuous_field(field.discrete_points[idx])
end
Out[28]:
In [29]:
discrete_points = 1.0/sqrt(3.0)*Vector[[-1, -1], [1, -1], [1, 1], [-1, 1]]
fs["discrete field 2"] = MyDiscreteField(discrete_points, fs["continuous field"])
last(fs["discrete field 2"])
Out[29]:
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.
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.
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.
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.
Interpolation of fields works of course both in time and spatial dimension. Here's a simple example showing the main concept.
In [30]:
using JuliaFEM: Basis
Yep, it's just continuous field:
In [31]:
Basis <: ContinuousField
Out[31]:
Basis + DiscreteField = Interpolation:
In [32]:
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])]
N = Basis(basis, dbasis)
Out[32]:
In [33]:
N([0.0, 0.0])
Out[33]:
In [34]:
t = linspace(0, 2, 5)
y = 1/2*t.^2
t, y
Out[34]:
In [35]:
timesteps = TimeStep[]
for (ti, yi) in zip(t, y)
increment = Increment(yi)
push!(timesteps, TimeStep(ti, increment))
end
#fs["particle position"] = t, 1/2*t.^2
pos = Field(timesteps)
pos(1.0)
Out[35]:
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:
In [36]:
velocity = pos(1.0, Val{:diff})
Out[36]:
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.
Let's have a $\Omega = \left[0,1\right]\times\left[0,1\right] \in \mathbb{R}^2$ domain with a displacement field \begin{equation} \mathbf{u}\left(X_1, X_2\right)=t\begin{bmatrix}X_{1}\left(X_{2}+1\right)\\ X_{1}\left(4X_{2}-1\right) \end{bmatrix}. \end{equation}
Gradient is \begin{equation} \mbox{grad}\left(\mathbf{u}\right)=t\begin{bmatrix}X_{2}+1 & X_{1}\\ 4X_{2}-1 & 4X_{1} \end{bmatrix} \end{equation}
We look for a displacement in center point of the domain at time $t=1.0$:
In [37]:
X = Field(Vector{Float64}[[0.0,0.0], [1.0,0.0], [1.0,1.0], [0.0,1.0]])
u = Field(
(0.5, Vector{Float64}[[0.0, 0.0], [0.5, -0.5], [1.0, 1.5], [0.0, 0.0]]),
(1.5, Vector{Float64}[[0.0, 0.0], [1.5, -1.5], [3.0, 4.5], [0.0, 0.0]]));
In [38]:
# evaluate fields in time t=1.0 -> Increments
X_increment = X(1.0)
u_increment = u(1.0)
X_increment, u_increment
Out[38]:
In [39]:
# evaluate geometry + displacement at midpoint, basically x = X + u
N(X_increment, [0.0, 0.0]) + N(u_increment, [0.0, 0.0])
Out[39]:
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}$:
In [40]:
N = Basis(basis, dbasis)
gradu = N(X_increment, u_increment, [0.0, 0.0], Val{:grad})
Out[40]:
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:
In [41]:
du = u(1.0, Val{:diff})
grad_du = N(X_increment, du, [0.0, 0.0], Val{:grad})
strain_rate = 1/2*(grad_du + grad_du')
Out[41]:
In [42]:
u = Field("displacement")
Out[42]:
The idea is to write equation in expression form using tiny symbolic interface and evaluate equations afterwards. Main benefits:
Here's a preliminary example how this could work:
In [43]:
# in future this will be something like this:
# u = Field("displacement")
# eps = 1/2*(grad(u) + grad(u)')
# strain_rate = diff(eps)
In [44]:
# now this is
expr = JuliaFEM.Expression("1/2*(grad(diff(displacement)) + grad(diff(displacement))')")
print(expr.expr)
Evaluation of expression: call Basis with FieldSet, $\boldsymbol\xi$ and time as arguments:
In [45]:
X = Field(Vector[[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
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]]))
fs = FieldSet()
fs["geometry"] = X
fs["displacement"] = u
xi = [0.0, 0.0]
time = 1.0
result = N(expr, fs, xi, time);
In [46]:
result
Out[46]:
What happened is that fields were automatically interpolated and substituted to the equation. The result is Expr and can be evaluated like usual:
In [47]:
strain_rate = eval(result)
Out[47]: