2D Finite Elements

Time-dependent solution

Version 2022.2
Mass matrix calculation
• The calculation of the mass matrix follows a logic similar
to the one used for the stiffness matrix:
1 Dx3
For each triangle e Dy3
For each local vertex i ( 1 → 3)
2
ii=Dof Numbering(i) Dy2 e
If vertex ii is unknown Dy1
For each local vertex j ( 1 → 3)
jj=Dof Numbering(j) Dx2 3 Dx1
If vertex jj is unknown
M(ii,jj)=M(ii,jj)+rho(e)*Area(e)*(1+δ(i,j))/12
M=[Link](Me)

V =[Link]; Dof =[Link]; numDof =max(Dof);

Areas=[Link]; Mcontr=[Link];
pos = 1; row = zeros(Mcontr, 1); col = zeros(Mcontr, 1);
m = zeros(Mcontr, 1); rho = [Link];
for e = 1:size(V, 1) %main loop over each triangle
for ni = 1:3 %for each vertex of this triangle
ii = Dof(V(e, ni));
if ii > 0 %is it a degree of freedom? Yes
for nj = 1:3 %second loop
jj = Dof(V(e, nj));
if jj > 0 %add to the mass matrix
m(pos) = 1/12 * Areas(e) * rho(e) * ...
((ii == jj) + 1);
row(pos) = ii; col(pos) = jj; pos = pos + 1;
end
A
end
 se i  j
end 12
end M = 
end  A se i = j
M = sparse(row, col, m, numDof, numDof);  6
Mass matrix calculation with lumping
M=[Link](Me)

• The lumping technique allows to build a diagonal mass

matrix: we use therefore a vector to store the data
• This approach is particularly convenient when we use
Explicit Euler method
rho = [Link];
V=[Link];
Dof=[Link];
NumDof = max(Dof);
M = zeros(NumDof,1);  0 sei  j
  
Nh
for e=1:size(Tr,1)
for ni=1:3 M = diag   b  =  A
sei = j
ij
ii = Dof(V(e,ni));   
i =1
if ii > 0 3
M(ii) = M(ii) + Areas(e) / 3 * rho(e);
end
end
end
Eigenfunctions evaluation
• We need to calculate the eigenvalues of M-1A
• Specific algorithms exist, which allow to calculate
eigenvalues and eigenvectors without calculating the
inverse of M
They are implemented, e.g., in the function eigs
[V,E]=eigs(A,M,2,'SM');
Autovalore 2: 4.1813

0.5

-0.5

-1
1
0.5 1
0 0.5
0
-0.5 -0.5
-1 -1
Temporal analysis
• When we consider a first order time derivative, we can
use the well known integration methods, with M mass
matrix
𝜕𝑢 𝜕𝑢
𝜌 − ∇ ⋅ 𝜇∇𝑢 = 𝑓 𝑥, 𝑦; 𝑡 → 𝑀 = −𝐷𝑢 + 𝑏
൞ 𝜕𝑡 𝜕𝑡
𝑢 𝑡 = 0 = 𝑢0
• Explicit Euler: Mu k +1 = Mu k − Dtu k + b k t
→ convergence problems if dt too large
• Implicit Euler: (M + Dt )u = Mu + b t
k +1 k k +1

 Dt  k +1  D b k +1
b k

M + u = Mu −  u − − t
k k
• CN:
 2  2 2 2
 T
Example #1: rcv −   (kT ) = f
t
heating of a squared domain
• We study the temporal evolution of the temperature on an
rectangular aluminium plate with a hole in the middle; the
following B.C.s are applied
– Homogeneous Neumann on two opposite external sides (perfect
insulation)
– Homogenerous Dirichlet on the two other external sides
– Dirichlet, with a temperature of 20°C on the internal edges
• For the considered material we have
– density r: 2700 kg/m3
– specific heat capacity cv: 896.9 J/(kg K)
– thermal conductivity k: 237 W/(m K)  0 if t  0

• An external source is applied: f (t ) =  30t if 0  t  1000
30000 if t  1000

Example #1: definition of geometry,
B.C. and mesh
1\3

S=[Link]('mu', 237/2700/897)- 2\4

2\5
2\6
2\7
2\2
2\3 2\64
2\1 2\62
2\61
2\60
2\59
2\58
2\8 2\63 2\57

[Link]([0,0],[0.2,0.2],64);
2\9
2\10 2\56
2\55
2\11
2\12 2\54
2\53
2\13 2\52
2\14 2\51

1\2 2\15
2\16
2\17
2\50
2\49
2\48
1\4
figure;[Link]('e');
2\18
2\19 2\47
2\46
2\20
2\21 2\45
2\44
2\22
2\23 2\43
2\42
2\24
2\25 2\41
2\40
2\26
2\27 2\39
2\38
2\28
2\29
2\30
2\31 2\37
2\36
2\35
2\34
2\32
2\33

[Link](1).Bc([1,3]) = 1\1
[Link](0);
[Link](2).Bc(:) =
[Link](20);
figure;[Link]('bc');

Me=mesh2D(S,0.001);
figure;[Link]('d');
 Example #1: Initial solution, f=0
• We need to modify the constant terms vector at each time
step: we can “split” it in two parts: b=bConst+bVar(t)
– bBC, related to the BC.s, is constant, since the applied BC.s do not
change with respect to time
– bForce is time dependant and follows the assigned force law
D = [Link](Me);
bBC = [Link](Me);
InitialForce = 3000;
bForce =[Link](Me,InitialForce);
uStationary0 = D\(bBC+bForce);
uu = [Link](uStationary0);
figure;
[Link](uu);
ylabel('Temperature [^\circC]');
colorbar();
Example #1: Temporal evolution

M=buildMass(Me);
Tend=4000; %end time, in s
%Temporal evolution of the force
f=@(t)InitialForce+(t<1000)*(30*t)+(t>=1000)*30000;
figure;
fplot(f,[0,100]);
 Example #1: Temporal evolution
• All the discussed methods follow the following steps:
– Choice of the time step, e.g. dt = 0.1;
– Definition of the initial temperatures u = uStationary0;
– Temporal loop: for k = 1:Tend/dt,
– Save indices of the nodes which are dof: Dofs=[Link]>0;
– In each time considered step:
• calculate the new force T4=f(t);
• Solve the linear system to calculate the new solution
u=…………………
uu(Dof) = u; %update only dofs, it’s faster!
• (if required) redraw the solution, for instance with
[Link](uu,'hidemesh');%I hide the mesh
title(['t= ' num2str(t) 's']); %elapsed seconds
zlim([0 20]); %fixed vertical range
caxis([0 20]); %fixed color scale
view ([0 90]); %top view
drawnow(); %force MATLAB to immediately redraw
Explicit Euler method (EE)
k +1
Mu = Mu − Dtu + b t
k k k

dt=0.1;
u=uStationary0;
Ddt=D*dt;
for k=1:Tend/dt,
t=k*dt;
F=f(t-dt); %previous step!
u=u-M\(Ddt*u-dt*(bBC+bForce*(F/InitialForce));
uu(Dofs)=u;
[Link](uu,'hidemesh'); hold off;
zlim([0,20]);
caxis([0,20]);
view([0,90]);
title(['t=',num2str(t) 's']); drawnow();
end
NB! Please pay attention to the correct choice of dt: dt  2 maxwith max
largest eigenvalue, in module, of M −1 A , to be calculated as
eigs(A,M,1,'LM') and NEVER taking the inverse of M!
Implicit Euler method (IE)
(M + Dt )u k +1 = Muk + bk +1t
dt = 0.1;
u = uStationary0;
S = (M+D*dt);
for k=1:Tend/dt,
t = k*dt;
F = f(t); %current step!
u = S\(M*u+dt*(bBC+bForce*(F/InitialForce)));
uu(Dofs) = u;
[Link](uu,'hidemesh');
hold off;
zlim([0,20]);
caxis([0,20]);
view([0,90]);
title(['t=',num2str(t) 's']);
drawnow();
end
Crank-Nicholson method (CN)
 Dt  k +1  D k b k +1 + b k 
M + u = Mu −  u − t
k

 2  2 2 
dt=10;
u=uStationary0;
B=D*dt/2; CN1=M+B; CN2=M-B;
Fold=f(0);
for k=1:Tend/dt,
t=k*dt;
F=f(t);
Coeff=(F+Fold)/InitialForce/2;
u=CN1\(CN2*u+dt*(bBC+bForce*Coeff));
uu(Dofs)=u;
[Link](uu,'hidemesh');
hold off;
zlim([0,20]);
caxis([0,20]);view([0,90]);
title(['t=',num2str(t) 's']);drawnow();
Fold=F;
end
Temporal evolution
• At every time step, we need to solve a linear system that
allows to calculate the new approximation of the solution
• If the mass matrix (EE) or both the stifness and mass
matrices (IE, CN) do not change in time, it is generally
very efficient to evaluate their preconditioners (see
functions ilu and ichol).
• This approach, used in conjunction with an iterative
method starting from the solution at the previous step,
can significantly reduce the required computation time
Implicit Euler method (IE) with
preconditioning
dt = 0.1;
u = uStationary0;
S = (M+D*dt);
%u=S\(M*u+bdt*T4); since at each step I must solve the
%system Su=..., then I pre-factor S
opt = struct('type','ict','droptol',1e-4);
R = ichol(S,opt);
for k=1:Tend/dt,
t = k*dt;
F = f(t);
[u, flag] = pcg(S,M*u+dt(bBC+bVar*F/InitialForce),
1e-6,1000,R',R,u);
uu(Dofs) = u;
[Link](uu,'hidemesh');zlim([0,20]);
caxis([0,20]);view([0,90]);
title(['t=',num2str(t) 's']);drawnow();
end
NB: to estimate the time reduction, do not draw the solution!
ODExx functions
The syntax is the same for all the functions of the family
[t,y]=ode45(dydt,[T0, Tend], y0, odepar)
•t column vector of N components containg the time values the solution
was calculated
•y matrix of N rows, each one containing the solution calculated at the
corresponding time step indicated in the t vector
•dydt function to integrate. It MUST receive two input parameters, (time
and y), like dydt=@(t,y)siny)+t
•[T0, Tend] 2 components vector, indicating the initial and final integration
time
•y0 problem initial values
•odepar optional variable created by the odeset function to set
additional integration parameters such as relative or absolute tolerances,
mass matrix, maximum time step,...
odepar=odeset('Absrel',1e-6,'AbsTol',1e-3', 'NonNegative',1,
'MaxStep',1e-1, 'Mass',M);
ODExx functions
• In MATLAB, to numerically integrate (systems of)
differential equations the odexx functions are available

• All these functions have the same calling syntax, but

implement different algorithms
ODExx functions: examples
• Let’s consider some examples of Cauchy problems
integrated over the interval [0,1]
• dy dt + a sin( t ) − y 2
= 0 dy dt = − a sin( t ) + y 2
→
 y (0) = 0  y (0) = 0
>>[t,y]=ode45(@(t,y)-a*sin(t)+y.^2,[0,1],0);
• dy dt = y2
d 2 y dt 2 − exp(t + y ) = 0  1
 dy2 dt = exp(t + y1 )
 y (0) = 0 →
dy dt = 2  y1 (0) = 0
 t =0  y2 (0 ) = 2
>>[t,y]=ode45(@(t,y)[y(2);exp(t+y(1)),[0,1],[0;2]);
Example #1: M
u
= − Du + b
ODExx functions t
u0 = uStationary0;
figure;
fode = @(t,u)-D*u+(bBC+bForce*f(t)/InitialForce);
odepar = odeset('Mass',M); %Mass matrix
[t,U] = ode45(fode, [0,Tend],u0,odepar);
size(U) % 24809 x 769 -> 769 time steps
for k = 1:10:length(t) %now simply plot the solution
uu(Dofs) = U(k,:);
[Link](uu,'hidemesh');
zlim([0 TMax]);
caxis([0 TMax]);
% view([0 90]);
title(['t= ' num2str(t(k)) 's']);
drawnow();
end
 Time varying Dirichlet B.C.s
• Let’s consider the system M u t + Du = b with M and
D built on all the nodes (and not only the internal ones)
• Let nodes 1 and 3 be dof, node 2 a node on an edge with
time invariant Dirichlet B.C (u2) and node 4 on an edge with
a time-varying Dirichlet B.C. (u4(t))
 m11 m12 m13 m14  u1 t   m11 0 m13 0  u1 t   m12 u2 t + m14 u4 t 
       
 m21 m22 m23 m24  u2 t   m21 m22 m23 m24  u2 t   0 
m = +
m32 m33 m34  u3 t   m31 0 m33 0  u3 t   m42 u2 t + m44 u4 t 
 31       
m m44  u4 t   m41 m42 m44  u4 
t    
 41 m42 m43 m43 0 
• Moving to the system involving only the dof we finally have
 m11 m13  u1 t   m12   m14   m11 m13  u1 t 
   + u2 t   + u4 t   =    + u4 t mvar,4
 31
m m33  u 3 t      42 
m  44   31
m m m33  u 3 t 
=0

• We need to take into account this additional contribution

Example #2:
heating of a squared domain
• We study the same problema as Example#1 but the
following B.C.s are now applied:
– Homogeneous Neumann on two opposite external sides (perfect
insulation)
– Non homogeneous Dirichlet on the two other external sides, where
T=20°
– Dirichlet on the internal edges, with the law
0 if t  0

f D (t ) = 20 +  8t if 0  t  10
80 if t  10

• No external source is applied
Example #2: definition of geometry,
B.C. and mesh
1\3

S=[Link]('mu', 237/2700/897)- 2\4

2\5
2\6
2\7
2\2
2\3 2\64
2\1 2\62
2\61
2\60
2\59
2\58
2\8 2\63 2\57

[Link]([0,0],[0.2,0.2],64);
2\9
2\10 2\56
2\55
2\11
2\12 2\54
2\53
2\13 2\52
2\14 2\51

1\2 2\15
2\16
2\17
2\50
2\49
2\48
1\4
figure;[Link]('e');
2\18
2\19 2\47
2\46
2\20
2\21 2\45
2\44
2\22
2\23 2\43
2\42
2\24
2\25 2\41
2\40
2\26
2\27 2\39
2\38
2\28
2\29
2\30
2\31 2\37
2\36
2\35
2\34
2\32
2\33

[Link](1).Bc([1,3]) = 1\1
[Link](0);
%Assign group 1 to the edges
%of the hole
[Link](2).Bc(:) =
[Link](1,1);
figure;[Link]('bc');
Me=mesh2D(S,0.001);
figure;[Link]('d');
Example #2
• We need to split vector b into a constant term (group 0) and a time
varying term (for the nodes on the internal border, group 1))
b(t ) = bconst + bvar (t ) = bconst + bvar 1 f D (t )
• Moreover, we need to generate a proper vector mvar to take into
account the contribution from the mass matrix:
D = [Link](Me);
[bConst, bVars] = [Link](Me);
[M, mVars] = [Link](Me);
• bVars is a matrix: each column contains the vector associated to
the B.C. of the corresponding group with index > 0
– In this case, since we have just one group of edges with a
N
time dependent B.C., bVars is a vector; in general
>>bVar=bVars(:,1); D D D
• Same behavior for BuildMass: mVar=mVars(:,1);

N
Example #2: stationary solution
• We solve the linear system in order to obtain the initial distribution
of the temperature:
T0=20; %initial temperature on the internal border
uStationary0=D\(bConst+bVar*T0);
uu=zeros(size([Link].X));
uu(Dof>0)=uStationary0;
uu(DirichletNodesCircle)=T0;
uu(DirichletNodesExtBorder)=20;
• Since all the Dirichlet edges have a T=20°C
and there is no external contribution, the
temperature in each node of the domain
is exactly 20°C.
Example #2: temporal evolution
• Using the Implicit Euler method:
u
+ mvar D + Du = b(t ) →
df
M
t dt
(M + Dt )u k +1 = Mu k + t (bconst + bvar f D k +1 )− mvar ( f D k +1 − f D k )
• Therefore:
u=uStationary0; dt=1; A=(M+D*dt);
TCircleOld=fDirichlet(0);
for k=1:Tend/dt,
TCircle=fDirichlet(k*dt); DeltaT=Tcircle-TCircleOld;
u=A\(M*u+dt*(bConst+bVar*TCircle)-mVar*DeltaT);
uu(Dof)=u;
uu(DirichletNodesCircle)=TCircle;
TCircleOld=TCircle;
end
Temporal analysis of elliptic
problems (Full\06ElasticMembrane)
• Let’s consider the temporal evolution of the elastic
membrane (elliptic problem)
  2u
 2 + u +  u = f
 t

 u ( 0 ) = u0
 u ' ( 0 ) = u1



• We obtain the following system Mu(tn ) + ( D + R)u(tn ) = f (tn )
with M mass matrix, D diffusion term and R reaction term

• For simplicity, we assume constant B.C.s

Example: script RunDt
• We study the evolution of the position of a squared
membrane placed in a viscous fluid, starting from a zero
displacement, when a constant force is applied from t=0
Sh = [Link]('mu',1);
Me = mesh2D(Sh);
f = @(x,y)-4*ones(size(x));
D = [Link](Me);
b = [Link](f);
M = [Link](Me);
%calculate and plot the stationary solution
uu = [Link](pcg(A,b,1e-3,250));
figure;
[Link](uu,'hidemesh');
Dofs = [Link]>0;
Example: script RunDt
%Numerical evolution, dumped oscillation
dt = 0.01; %time steps
Tend = 10; %simulation end time
NumIter = ceil(Tend/dt); %number of steps
displ=zeros(NumIter,1); %to store the max displacement
u0 = zeros(size(b)); %null initial displacement
u1 = u0;
a = 1.5; %fluid friction coefficient
for i = 1:NumIter
[u2,flag]=pcg(M*(1+a*dt),M*((2+a*dt)*u1-u0)-
dt^2*(D*u1-b),[],30,[],[],u1);
u0=u1; u1=u2;
uu(Dofs)=u2;
[Link](uu,'hidemesh'); zlim([-1 1]);
drawnow; %to force a refresh of the figure,
%otherwise it is updated at the end of the loop
displ(i)=max(abs(u2)); %save the maximum displac.
end
Example: script RunDt
figure;
plot((1:NumIter)*dt, displ);
title('Maximum displacement');
xlabel('t [A.U.]');
ylabel('u [A.U.]'); 0.5
Maximum displacement

grid on; 0.45

0.4

0.35

0.3
u [A.U.]

0.25

0.2

0.15

0.1

0.05

0
0 1 2 3 4 5 6 7 8 9 10
t [A.U.]
Example 3: charging/discharging
heat into/from a sphere array
Continuation of the stationary analysis in Example 9