Skip to content

Commit 80c6d65

Browse files
added continuous example to FEM and fixed parametric interface docstrings (#118)
* added continuous example to FEM and fixed parametric interface docstrings * a few changes in the FEM example (#119) * editing FEM text * temporarily disable tests * more changes * added more description on the continuum problem * reactivate tests * update title Co-authored-by: lucaferranti <49938764+lucaferranti@users.noreply.github.com> Co-authored-by: Jorge Pérez Zerpa <42485529+jorgepz@users.noreply.github.com>
1 parent 9301319 commit 80c6d65

File tree

4 files changed

+263
-36
lines changed

4 files changed

+263
-36
lines changed

docs/literate/applications/FEM_example.jl

Lines changed: 229 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -1,20 +1,20 @@
1-
# # A simple FEM application
2-
#
3-
# ## Introduction
4-
#
5-
# ### General description
6-
# The Finite Element Method is widely used to solve PDEs in Engineering applications and particularly in Structural Analysis problems [[BAT14]](@ref). The procedure consists discretizing the domain into _elements_ and assembly a system of balance equations. For linear problems, this system can be usually written as
7-
#
1+
# # Application of Interval Linear Algebra to FEM analysis
2+
# The Finite Element Method is widely used to solve PDEs in Engineering applications and particularly in Structural Analysis problems [[BAT14]](@ref). The procedure consists in discretizing the domain into _elements_ and constructing (_assembling_) a system of balance equations. For linear problems, this system can be usually written as
83
# ```math
94
# K \cdot d = f
105
# \qquad
11-
# K = \sum_e K_e
6+
# K = \sum_{e=1}^{n_e} K_e
127
# ```
13-
# where $f$ is the vector of external loads, $K_e$ is the element stiffness matrix and $d$ is the vector of unknown displacements.
8+
# where $n_e$ is the number of elements of the domain, $f$ is the vector of external loads, $K_e$ is the stiffness matrix of element $e$ in global coordinates, $K$ is the assembled stiffness matrix and $d$
9+
# is the vector of unknown displacements. This tutorial shows how IntervalLinearAlgebra can
10+
# be used to solve structural mechanics problems with uncertainty in the parameters. Particularly,
11+
# it highlights the importance of parametric interval linear systems.
1412
#
15-
# ### FEM for truss structures
16-
# A frequent and simple type of structures are _Truss structures_, which are formed by bars connected but not welded. Truss models are considered during the conceptual design of bridges or other structures.
13+
# ## Simple truss structure
1714
#
15+
# A frequent and simple type of structures are _Truss structures_, which are formed by bars connected but not welded. Truss models are usually considered during the conceptual design of bridges or other structures.
16+
#
17+
# ### Stiffness equations
1818
# The stiffness matrix of a truss element in the local coordinate system is given by
1919
# ```math
2020
# K_L = s
@@ -28,7 +28,7 @@
2828
# \right),
2929
# ```
3030
#
31-
# where $s =\frac{E A}{L}$, with $E$ being the Young modulus, $A$ the area of the cross-section and $L$ the length of that truss element.
31+
# where $s =\frac{E A}{L}$ is the stiffness, $E$ is the Young modulus, $A$ is the area of the cross-section and $L$ is the length of that truss element.
3232
#
3333
# The change-of-basis matrix is given by
3434
# ```math
@@ -40,17 +40,18 @@
4040
# 0 & 0 & \cos(\alpha) & -\sin(\alpha) \\
4141
# 0 & 0 & \sin(\alpha) & \cos(\alpha)
4242
# \end{matrix}
43-
# \right),
43+
# \right).
4444
# ```
4545
#
4646
# The system of equations for each element is written in local coordinates as
4747
# ```math
4848
# K_L d_L = f_L
4949
# ```
50-
# and using the change-of-basis we obtain
50+
# and using the change-of-basis we obtain the equations for that element in the global systems of coordinates
5151
# ```math
52-
# K_G d_G = f_G \qquad K_G = Q K_L Q^T
52+
# K_G d_G = f_G \qquad K_G = Q K_L Q^T.
5353
# ```
54+
# After the system of equations for each element is in global coordinates, the whole system is assembled.
5455
#
5556
# The unitary stiffness matrix (for $s=1$) can be computed using the following function.
5657
function unitaryStiffnessMatrix( coordFirstNode, coordSecondNode )
@@ -64,39 +65,37 @@ function unitaryStiffnessMatrix( coordFirstNode, coordSecondNode )
6465
return Kglo, length
6566
end
6667
#
67-
# ## Example problem
68+
# ### Example problem
6869
#
6970
# A problem based on Example 4.1 from [[SKA06]](@ref) is considered. The following diagram shows the truss structure considered.
7071
#
7172
# ```@raw html
7273
# <img src="../../assets/trussDiagram.svg" style="width: 100%" alt="truss diagram"/>
7374
# ```
7475
#
75-
# ### Case with fixed parameters
76+
# #### Case with fixed parameters
7677
#
7778
# The scalar parameters considered are given by
7879
E = 2e11 ; # Young modulus
7980
A = 5e-3 ; # Cross-section area
8081
# The coordinate matrix is given by
81-
## coordinates x y
8282
nodesCMatrix = [ 0.0 0.0 ;
8383
1.0 1.0 ;
8484
2.0 0.0 ;
8585
3.0 1.0 ;
8686
4.0 0.0 ];
8787
# the connectivity matrix is given by
88-
## connectivity start end
8988
connecMatrix = [ 1 2 ;
9089
1 3 ;
9190
2 3 ;
9291
2 4 ;
9392
3 4 ;
9493
3 5 ;
95-
4 5 ]
94+
4 5 ];
9695
# and the fixed degrees of freedom (supports) are defined by the vector
97-
fixedDofs = [2 9 10 ]
98-
#
99-
# calculations
96+
fixedDofs = [ 2 9 10 ];
97+
98+
# The number of elements and nodes are computed, as well as the free degrees of freedom.
10099
numNodes = size( nodesCMatrix )[1] # compute the number of nodes
101100
numElems = size( connecMatrix )[1] # compute the number of elements
102101
freeDofs = zeros(Int8, 2*numNodes-length(fixedDofs))
@@ -108,12 +107,11 @@ while indDof <= (2*numNodes)
108107
end
109108
global indDof = indDof + 1
110109
end
111-
#
112-
# assembly
110+
111+
# The global stiffness equations are computed for the unknown displacements (free dofs)
113112
KG = zeros( 2*numNodes, 2*numNodes )
114113
FG = zeros( 2*numNodes )
115114
for elem in 1:numElems
116-
print(" assembling stiffness matrix of element ", elem , "\n")
117115
indexFirstNode = connecMatrix[ elem, 1 ]
118116
indexSecondNode = connecMatrix[ elem, 2 ]
119117
dofsElem = [2*indexFirstNode-1 2*indexFirstNode 2*indexSecondNode-1 2*indexSecondNode ]
@@ -129,14 +127,14 @@ FG[4] = -1e4 ;
129127
KG = KG[ freeDofs, : ]
130128
KG = KG[ :, freeDofs ]
131129
FG = FG[ freeDofs ]
132-
130+
# and the system is solved.
133131
u = KG \ FG
134132
UG = zeros( 2*numNodes )
135133
UG[ freeDofs ] = u
136134
#
137-
# #### Deformed structure
138-
#
139-
# The reference (dashed blue line) and deformed (solid red) configurations of the structure are ploted. Since the displacements are very small, a `scaleFactor` is considered to amplify the deformation and ease the visualization.
135+
# The reference (dashed blue line) and deformed (solid red) configurations of the structure are ploted.
136+
# Since the displacements are very small, a `scaleFactor` is considered to amplify
137+
# the deformation and ease the visualization.
140138
#
141139
using Plots
142140
scaleFactor = 2e3
@@ -164,11 +162,12 @@ savefig("deformed.png") # hide
164162
#
165163
# ![](deformed.png)
166164
#
167-
# ### Problem with interval parameters
165+
# #### Problem with interval parameters
168166

169167
md"""
170168
Suppose now we have a 10% uncertainty for the stiffness $s_{23}$ associated with the third
171-
element. To model the problem, we introduce the symbolic variable `s23` using the IntervalLinearAlgebra macro `@linvars`.
169+
element. To model the problem, we introduce the symbolic variable `s23` using the
170+
IntervalLinearAlgebra macro [`@affinevars`](@ref).
172171
"""
173172

174173
using IntervalLinearAlgebra
@@ -206,12 +205,12 @@ KGp = AffineParametricArray(KGp)
206205
srange = E * A / sqrt(2) ± 0.1 * E * A / sqrt(2)
207206

208207
# To solve the system, we could of course just subsitute `srange` into the parametric matrix
209-
# `KGp` and solve the "normal" interval linear system
208+
# `KGp` and solve the "normal" interval linear system (naive approach)
210209

211210
usimple = solve(KGp(srange), Interval.(FG))
212211

213-
# This approach, however soffers from the [dependency problem](https://en.wikipedia.org/wiki/Interval_arithmetic#Dependency_problem)
214-
# and hence the compute displacement will be an overestimation of the true displacement.
212+
# This approach, however suffers from the [dependency problem](https://en.wikipedia.org/wiki/Interval_arithmetic#Dependency_problem)
213+
# and hence the computed displacements will be an overestimation of the true displacements.
215214

216215
# To mitigate this issue, algorithms to solve linear systems with parameters have been developed.
217216
# In this case we use the algorithm presented in [[SKA06]](@ref)
@@ -221,3 +220,198 @@ uparam = solve(KGp, FG, srange)
221220
# We can now compare the naive and parametric solution
222221

223222
hcat(usimple, uparam)/1e-6
223+
224+
md"""
225+
As you can see, the naive non-parametric approach significantly overestimates the displacements.
226+
It is true that for this very simple and small structure, both displacements are small, however
227+
as the number of nodes increases, the effect of the dependency problem also increases and
228+
the non-parametric approach will fail to give useful results. This is demonstrated in the
229+
next section.
230+
"""
231+
# ## A continuum mechanics problem
232+
233+
md"""
234+
In this problem a simple solid plane problem is considered. The solid is fixed on its bottom edge and loaded with a shear tension on the top edge.
235+
"""
236+
237+
# First, we set the geometry and construct a regular grid of points
238+
239+
L = [1.0, 4.0] # dimension in each direction
240+
t = 0.2 # thickness
241+
nx = 10 # number of divisions in direction x
242+
ny = 20 # number of divisions in direction y
243+
nel = [nx, ny]
244+
neltot = 2 * nx * ny; # total number of elements
245+
nnos = nel .+ 1 # number of nodes in each direction
246+
nnosx, nnosy = nnos
247+
nnostot = nnosx * nnosy ; # total number of nodes
248+
249+
# we compute the vector of indexes of the loaded nodes (bottom ones)
250+
251+
startloadnode = (nnosy - 1) * nnosx + 1 # boundary conditions
252+
endinloadnode = nnosx * nnosy
253+
LoadNodes = startloadnode:endinloadnode
254+
lins1 = range(0, L[1], length=nnosx)
255+
lins2 = range(0, L[2], length=nnosy)
256+
257+
# and construct the matrix of coordinates of the nodes
258+
259+
nodes = zeros(nnostot, 2) # nodes: first column x-coord, second column y-coord
260+
for i = 1:nnosy # first discretize along y-coord
261+
idx = (nnosx * (i-1) + 1) : (nnosx*i)
262+
nodes[idx, 1] = lins1
263+
nodes[idx, 2] = fill(lins2[i], nnosx)
264+
end
265+
266+
# The connectivity matrix Mcon is computed, considering 3-node triangular elements
267+
268+
Mcon = Matrix{Int64}(undef, neltot, 3); # connectivity matrix
269+
for j = 1:ny
270+
for i = 1:nx
271+
intri1 = 2*(i-1)+1+2*(j-1)*nx
272+
intri2 = intri1 + 1
273+
Mcon[intri1, :] = [j*nnosx+i, (j-1)*nnosx+i, j*nnosx+i+1 ]
274+
Mcon[intri2, :] = [j*nnosx+i+1, (j-1)*nnosx+i, (j-1)*nnosx+i+1]
275+
end
276+
end
277+
278+
# the undeformed mesh is plotted as follows
279+
280+
Xel = Matrix{Float64}(undef, 3, neltot); Yel = Matrix{Float64}(undef, 3, neltot)
281+
for i = 1:neltot
282+
Xel[:, i] = nodes[Mcon[i, :], 1] # the j-th column has the x value at the j-th element
283+
Yel[:, i] = nodes[Mcon[i, :], 2] # the j-th column has the y value at the j-th element
284+
end
285+
fig = plot(ratio=1, xlimits=(-1, 3), title="Undeformed mesh", xlabel="x", ylabel="y")
286+
plot!(fig, [Xel[:, 1]; Xel[1, 1]], [Yel[:, 1]; Yel[1, 1]], linecolor=:blue, linewidth=1.4, label="")
287+
for i = 2:neltot
288+
plot!(fig, [Xel[:, i]; Xel[1, i]], [Yel[:, i]; Yel[1, i]], linecolor=:blue, linewidth=1.4, label="")
289+
end
290+
savefig("undeformed2.png") # hide
291+
292+
# ![](undeformed2.png)
293+
294+
# Let us now define the material parameters. Here we assume a 10% uncertainty on the Young modulus, while Poisson ratio and the density are fixed. This can be related to a steel plate problem with an unknown composition, thus unknown exact Young modulus value.
295+
296+
ν = 0.25 # Poisson
297+
ρ = 8e3 # density
298+
299+
@affinevars E # Young modulus, defined as symbolic variable
300+
En = 200e9 # nominal value of the young modulus
301+
Erange = En ± 0.1 * En # uncertainty range of the young modulus
302+
303+
# We can now assemble the global stiffness matrix.
304+
# We set the constitutive matrix for a plane stress state.
305+
C = E / (1-ν^2) * [ 1 ν 0 ;
306+
ν 1 0 ;
307+
0 0 (1-ν)/2 ]
308+
309+
# We compute the free and fixed degrees of freedom
310+
function nodes2dofs(u)
311+
v = Vector{Int64}(undef, 2*length(u))
312+
for i in 1:length(u)
313+
v[2i-1] = 2u[i] - 1; v[2i] = 2u[i]
314+
end
315+
return v
316+
end
317+
FixNodes = 1:nnosx
318+
FixDofs = nodes2dofs(FixNodes) # first add all dofs of the nodes
319+
deleteat!(FixDofs, 3:2:(length(FixDofs)-2)) # then remove the free dofs of the nodes
320+
LibDofs = Vector(1:2*nnostot) # free degrees of fredom
321+
deleteat!(LibDofs, FixDofs)
322+
323+
# and we assemble the matrix
324+
325+
function stiffness_matrix(x, y, C, t)
326+
327+
A = det([ones(1, 3); x'; y']) / 2 # element area
328+
B = 1 / (2*A) * [y[2]-y[3] 0 y[3]-y[1] 0 y[1]-y[2] 0 ;
329+
0 x[3]-x[2] 0 x[1]-x[3] 0 x[2]-x[1] ;
330+
x[3]-x[2] y[2]-y[3] x[1]-x[3] y[3]-y[1] x[2]-x[1] y[1]-y[2] ]
331+
332+
K = B' * C * B * A * t ;
333+
return K
334+
end
335+
336+
KG = zeros(AffineExpression{Float64}, 2*nnostot, 2*nnostot);
337+
for i = 1:neltot
338+
Ke = stiffness_matrix(Xel[:, i], Yel[:, i], C, t)
339+
aux = nodes2dofs(Mcon[i, :])
340+
KG[aux, aux] .+= Ke
341+
end
342+
343+
K = AffineParametricArray(KG[LibDofs, LibDofs])
344+
345+
# Finally, we assemble the loads vector
346+
347+
areaelemsup = L[1] / nx * t
348+
349+
f = zeros(2*nnostot);
350+
351+
f[2*LoadNodes[1]] = 0.5 * areaelemsup;
352+
353+
for i in 2:(length(LoadNodes)-1)
354+
idx = 2 * LoadNodes[i]
355+
356+
f[idx-1] = 1 * areaelemsup # horizontal force
357+
end
358+
f[2*LoadNodes[end]] = 0.5 * areaelemsup
359+
360+
q = 1e9 # distributed load on the up_edge
361+
F = q * f;
362+
FLib = F[LibDofs]
363+
nothing # hide
364+
365+
# now we can solve the displacements from the parametric interval linear system and plot
366+
# minimum and maximum displacement.
367+
368+
u = solve(K, FLib, Erange) # solving
369+
370+
# plotting
371+
U = zeros(Interval, 2*nnostot)
372+
U[LibDofs] .= u
373+
374+
Ux = U[1:2:2*nnostot-1]
375+
Uy = U[2:2:2*nnostot]
376+
377+
nodesdef = hcat(nodes[:, 1] + Ux, nodes[:, 2] + Uy);
378+
379+
Xeld = Interval.(copy(Xel))
380+
Yeld = Interval.(copy(Yel))
381+
382+
# build elements coordinate vectors
383+
for i = 1:neltot
384+
Xeld[:, i] = nodesdef[Mcon[i, :], 1] # the j-th column has the x coordinate of the j-th element
385+
Yeld[:, i] = nodesdef[Mcon[i, :], 2] # the j-th column has the y coordinate of the j-th element
386+
end
387+
388+
plot!(fig, [inf.(Xeld[:, 1]); inf.(Xeld[1, 1])], [inf.(Yeld[:, 1]); inf.(Yeld[1, 1])], linecolor=:green, linewidth=1.4, label="", title="Displacements")
389+
for i = 2:neltot
390+
plot!(fig, [inf.(Xeld[:, i]); inf.(Xeld[1, i])], [inf.(Yeld[:, i]); inf.(Yeld[1, i])], linecolor=:green, linewidth=1.4, label="")
391+
end
392+
393+
plot!(fig, [sup.(Xeld[:, 1]); sup.(Xeld[1, 1])], [sup.(Yeld[:, 1]); sup.(Yeld[1, 1])], linecolor=:red, linewidth=1.4, label="", title="Displacements")
394+
for i = 2:neltot
395+
plot!(fig, [sup.(Xeld[:, i]); sup.(Xeld[1, i])], [sup.(Yeld[:, i]); sup.(Yeld[1, i])], linecolor=:red, linewidth=1.4, label="")
396+
end
397+
savefig("displacement2.png") # hide
398+
399+
# ![](displacement2.png)
400+
401+
md"""
402+
In this case, ignoring the dependency and treating the problem as a "normal" interval linear
403+
system would fail. The reason for this is that the matrix is not strongly regular, which is
404+
a necessary condition for the implemented algorithms to work.
405+
"""
406+
407+
is_strongly_regular(K(Erange))
408+
409+
# ## Conclusions
410+
411+
md"""
412+
This tutorial showed how interval methods can be useful in engineering applications dealing
413+
with uncertainty. As in most applications the elements in the matrix will depend on some
414+
common parameters, due to the dependency problem neglecting the parametric structure will
415+
result in poor results. This highlights the importance of parametric interval methods in
416+
engineering applications.
417+
"""

docs/src/api/misc.md

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,14 @@ Pages = ["misc.md"]
1111
set_multiplication_mode
1212
```
1313

14+
## Symbolic Interface
15+
16+
```@docs
17+
@affinevars
18+
AffineExpression
19+
AffineParametricArray
20+
```
21+
1422
## Others
1523
```@autodocs
1624
Modules = [IntervalLinearAlgebra]

src/IntervalLinearAlgebra.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,7 +23,8 @@ export
2323
rref,
2424
eigenbox, Rohn, Hertz, verify_eigen, bound_perron_frobenius_eigenvalue,
2525
AffineExpression, @affinevars,
26-
AffineParametricArray, AffineParametricMatrix, AffineParametricVector
26+
AffineParametricArray, AffineParametricMatrix, AffineParametricVector,
27+
Skalna06
2728

2829

2930
include("linear_systems/enclosures.jl")

0 commit comments

Comments
 (0)