Vectorized Matlab Codes For Linear Two-Dimensional Elasticity
Vectorized Matlab Codes For Linear Two-Dimensional Elasticity
Vectorized Matlab Codes For Linear Two-Dimensional Elasticity
Elasticity
Jonas Koko
LIMOS, Universite Blaise Pascal CNRS UMR 6158
ISIMA, Campus des Cezeaux BP 10125, 63173 Aubi`ere cedex, France
koko@isima.fr
Abstract
A vectorized Matlab implementation for the P1 finite element is provided for the twodimensional linear elasticity with mixed boundary conditions. Vectorization means
that there is no loop over triangles. Numerical experiments show that our implementation is more efficient than the standard implementation with a loop over all
triangles.
Keywords: Finite element method, elasticity, Matlab.
AMS subject classification: 65N30, 65R20, 74B05
Introduction
Matlab is nowadays a widely used tool in education, engineering and research and becomes
a standard tool in many areas. But Matlab is a matrix language, that is, Matlab is designed
for vector and matrix operations. For best performance in extensive problems, one should
take advantage of this.
We propose a Matlab implementation of the P 1 finite element for the numerical solutions of two-dimensional linear elasticity problems. Instead of mixing finite element types
(e.g. [2, 1]), we propose a P1 -triangle vectorized code able to treat medium size meshes in
acceptable time. Vectorization means that there is no loop over triangles nor nodes. Our
implementation needs only Matlab basic distribution functions and can be easily modified
and refined.
The paper is organized as follows. The model problem is described in Section 2,
followed by a finite element discretization in Section 3. The data representation used in
Matlab programs is given in Section 4. The heart of the paper is the assembling functions
of the stiffness matrix in Section 5 and the right-hand side in Section 6. Numerical
experiments are carried out in Section 7 where post-processing functions are given. Matlab
programs used for numerical experiments are given in the appendix.
MODEL PROBLEM
Model problem
i, j = 1, 2,
(2.1)
(u) = (u + uT )/2,
where and denote Lame (positive) constants. These coefficients are related (in plane
deformations) to the Young modulus E and the Poisson coefficient by
E
,
(1 + )(1 2)
E
.
2(1 + )
in ,
(2.2)
on N ,
(2.3)
u = uD ,
on D .
(2.4)
v H 1 () : v = 0 on D ,
= v H 1 () : v = uD on D .
g vds,
v V.
(2.5)
In (2.5), C = (Cijkl ) is the fourth-order elastic moduli tensor corresponding to (2.1), i.e.
ij (u) =
2
X
i, j = 1, 2,
k,l=1
where
Cijkl = ij kl + (ik jl + il jk ),
ij being the Kronecker delta.
1 i, j, k, l 2,
Let Th be a regular (in the sense of Ciarlet [3]) triangulation of . Spaces V and V D are
then replaced by their discrete approximations V h and VhD defined by
v|T P1 (T ), T Th ; vh| = 0 ,
Vh = vh C 0 ();
D
D
0
Vh = vh C (); v|T P1 (T ), T Th ; vh|D = uD .
vh Vh .
(3.1)
Let {j } be the system of piecewise global basis functions of V h , i.e. for all uh =
(u1h , u2h ) Vh
N
X
j (x)uj , = 1, 2.
uh =
j=1
We set U =
j=1
u12
u21
u22
uN
1
uN
2 ),
The stiffness matrix A = (Aij ) and the right-hand side b = bi are then given by
Z
(i ) : C(j )dx,
Aij =
Z
Z
bi =
f i dx +
g i ds.
The stiffness matrix A is sparse, symmetric and positive semi-definite before incorporating
the boundary condition u|D = uD .
In practice, integrals in (3.1) are computed as sums of integrals over all triangles, using
= T T T
the fact that
h
Z
X
X Z
X Z
(uh ) : C(vh )dx =
f vh dx +
g vh ds, vh Vh .
T Th
T Th
If we set
(T )
Aij
EN
(3.3)
f i dx,
(3.4)
g i ds,
(3.5)
(T )
bi
(E)
bi
bi =
(T )
bi
T Th
(E)
bi
EN
For the mesh, we adopt the data representation used in Matlab PDE Toolbox [4]. Nodes
coordinates and triangle vertices are stored in two arrays p(1:2,1:np) and t(1:3,1:nt),
where np is the number of nodes and nt the number of triangles. The array t contains
for each element the node numbers of the vertices numbered anti-clockwise. For the
triangulation of Figure 1, the nodes array p is (np=9)
0.0000
1.0000
1.0000
1.0000
1.0000
0.0000
0.
0.
0.5000
1.0000
1.0000
0.5000
0.5000
0.0000
0.0000
0.5000
0.5000
0.5000
6
2
9
7
3
9
8
4
9
2
5
9
3
6
9
4
7
9
1
8
9
over all edges E results in a loop over all entries of ibcneum. Dirichlet boundary conditions
are provided by a list of nodes and a list of the corresponding prescribed boundary values.
Note that Matlab supports reading data from files in ASCII format (Matlab function
load) and there exists good mesh generators written in Matlab, see e.g. [6].
u1 /x
=
(u) =
u2 /y
u1 /x + u2 /y
strain tensor
11 (u)
22 (u) .
212 (u)
11 (uh )
22 (uh ) = C(uh ),
12 (uh )
where
(5.1)
+ 2
+ 2 0 .
C=
0
0
u1
1,x
0
2,x
0
3,x
0
1
1
0
1,y
0
2,y
0
3,y ... =:
Ru,
(5.2)
(uh ) =
2|T |
2|T |
1,y 1,x 2,y 2,x 3,y 3,x
u6
where |T | is the area of the triangle T . From (5.1) and (5.2), it follows that
(vh ) : C(uh ) = (vh )t C(uh ) =
1
vt Rt CRu
4|T |2
1
t1
1 1 1
t2 = x1 x2 x3
y1 y2 y3
t3
1
Rt CR.
4|T |
(5.3)
y2 y 3 x3 x 2
0 0
1
y3 y 1 x1 x 3 .
1 0 =
2|T |
y1 y 2 x2 x 1
0 1
Matlab implementation given in [2, 5] (in a more modular form) can be summarized by
the matlab Function 5.1. This implementation, directly derived from compiled languages
as Fortran or C/C++, produces a very slow code in Matlab for large size meshes due to
the presence of the loop for. Our aim is to remove this loop by reorganizing calculations.
Let us introduce the following notations
xij = xi xj ,
so that, from (5.2)
yij = yi yj ,
i, j = 1, 2, 3,
(5.4)
(T )
(T )
A
A
12
11
A(T ) =
(T )
(T )
A22
A21
(T )
(T )
(T )
(T )
where A11 and A22 are symmetric and A21 = (A21 )t . After laborious (but elementary)
calculations, we get
e 2 + x2
y
23
32
e
(T )
A11 = y
31 y23 + x13 x32
31
13
e 12 y23 + x21 x32 y
e 12 y31 + x21 x13
y
e 2 + y 2
x
32
23
e
(T )
A22 = x
13 x32 + y31 y23
13
31
e
e
x21 x32 + y12 y23 x21 x13 + y12 y31
2
2
e
y
+ x
12
21
2
2
e
x + y
21
12
e = + 2; and
where
y23 x32 + x32 y23 y23 x13 + x32 y31 y23 x21 + x32 y12
(T )
A12 = y31 x32 + x13 y23 y31 x13 + x13 y31 y31 x21 + x13 y12 .
y12 x32 + x21 y23 y12 x13 + x21 y31 y12 x21 + x21 y12
Let us introduce the following vectors
x32
x = x13
x21
(T )
(T )
y23
y = y31 .
y12
(T )
One can verify that matrices A11 , A22 and A12 can be rewritten in a simple form, using
x and y, as
(T )
= ( + 2)yyt + xxt
(T )
= ( + 2)xxt + yyt
(T )
= yxt + xyt
A11
A22
A12
that is, for 1 i, j 3
A11ij
(T )
= ( + 2)yi yj + xi xj ,
(5.5)
(T )
= ( + 2)xi xj + yi yj ,
(5.6)
(T )
= yi xj + xi yj
(5.7)
A22ij
A12ij
With Matlab, xi and yi can be computed in a fast way for all triangles using vectorization.
Then assembling the stiffness matrix reduces to two constant loops for computing xx t , yyt
(T )
(T )
(T )
and xyt . We do not need to assemble separately A 11 , A22 and A12 . Sub matrices (5.5)(5.7) are directly assembled in the global stiffness matrix, in its standard form, since we
know their locations. Matlab vectorized implementation of the Assembly of the stiffness
matrix is presented in Function 5.2.
We assume that the volume forces f = (f 1 , f2 ) are provided at mesh nodes. The integral
(3.5) is approximated as follows
Z
x2 x 1 x3 x 1
1
f (x , y ), j = mod(i 1, 2) + 1
f i dx det
y2 y 1 y3 y 1 j c c
6
T
where (xc , yc ) is the center of mass of the triangle T . With the assumption on f ,
fj (xs , ys ) = (fj (x1 , y1 ) + fj (x2 , y2 ) + fj (x3 , y3 ))/3
(6.1)
NUMERICAL EXPERIMENTS
where {(xi , yi )}i=1,3 are vertices of the triangle T . Using the notation convention (5.4),
we have
x2 x 1 x3 x 1
(6.2)
y2 y1 y3 y1 = x21 y31 x31 y12 .
We can compute (6.1) and (6.2) over all triangles using vectorization techniques, and
assemble the result with the Matlab function sparse. Matlab implementation of the
assembly of the volume forces is presented in Function 6.1.
Integrals involving Neumann conditions are approximated using the value of g at the
center (xc , yc ) of the edge E
Z
1
g i dx |E|gj (xc , yc ), j = mod(i 1, 2) + 1
2
E
As for volume forces, |E| and gj (xc , yc ) are computed over all triangles using vectorization
techniques and assembled using Matlab function sparse, Function 6.2.
Numerical experiments
NUMERICAL EXPERIMENTS
10
The Matlab function 7.1 displays the deformed mesh with an additional function sf issued
from post-processing. It is directly derived from the Matlab function Show given in [2].
The additional function can be stresses, strains, potential energy,... evaluated at mesh
nodes.
In our numerical experiments, the additional function used in the function 7.1 is either
the Von Mises effective stress or the shear energy density |dev h |2 /(4), where
1
2
2
2
+
(h11 + h22 )2 + 2(h12
h11 h22 ).
|devh | =
2 6( + )
The stress tensor h is computed at every node as the mean value of the stress on the
corresponding patch. Matlab function 7.2 calculates approximate shear energy density
and Von mises effective stress.
In all examples, Dirichlet boundary conditions are taken into account by using large
spring constants (i.e. penalization).
NUMERICAL EXPERIMENTS
11
7.1
L-shape
The L-shape test problem is a common benchmark problem for which the exact solution
is known. The domain is described by the polygon
(1, 1), (0, 2), (2, 0), (0, 2), (1, 1), (0, 0)
The exact solution is known in polar coordinates (r, ),
ur (r, ) =
u (r, ) =
1
r [(c2 1)c1 cos(( 1)) ( + 1) cos(( + 1))] ,
2
1
r [( + 1) sin(( + 1)) + (c2 + 1)c1 sin(( 1))] .
2
(7.1)
(7.2)
12
NUMERICAL EXPERIMENTS
Mesh
Triangles/Nodes
100/66
400/231
1600/861
6400/3321
25600/13041
102400/51681
2.5
1.5
0.5
0.5
1.5
2.5
1
0.5
0.5
1.5
2.5
NUMERICAL EXPERIMENTS
13
Function 7.2 Matlab function for computing the shear energy density Sed and the Von
Mises effective stress Vms
function [Sed,Vms]=elas2dsvms(p,t,uu,E,nu)
%----------------------------------------% Two-dimensional linear elasticity
% Sed Shear energy density
% Vms Von Mises effective stress
%---------------------------------------n=size(p,2); nn=2*n;
% Lame constant
lam=E*nu/((1+nu)*(1-2*nu)); mu=E/(2*(1+nu));
% area of traingles
it1=t(1,:); it2=t(2,:); it3=t(3,:);
x21=(p(1,it2)-p(1,it1)); x32=(p(1,it3)-p(1,it2)); x31=(p(1,it3)-p(1,it1));
y21=(p(2,it2)-p(2,it1)); y32=(p(2,it3)-p(2,it2)); y31=(p(2,it3)-p(2,it1));
ar=.5*(x21.*y31-x31.*y21);
% gradient of scalar basis functions
phi1=[-y32./(2*ar) x32./(2*ar)];
phi2=[ y31./(2*ar) -x31./(2*ar)];
phi3=[-y21./(2*ar) x21./(2*ar)];
% displacements
u=uu(1:2:end); v=uu(2:2:end);
uh=[u(it1) u(it2) u(it3)]; vh=[v(it1) v(it2) v(it3)];
% strains
e11=uh(:,1).*phi1(:,1)+uh(:,2).*phi2(:,1)+uh(:,3).*phi3(:,1);
e22=vh(:,1).*phi1(:,2)+vh(:,2).*phi2(:,2)+vh(:,3).*phi3(:,2);
e12=uh(:,1).*phi1(:,2)+uh(:,2).*phi2(:,2)+uh(:,3).*phi3(:,2)...
+vh(:,1).*phi1(:,1)+vh(:,2).*phi2(:,1)+vh(:,3).*phi3(:,1);
clear uh vh
% stresses
sig11=(lam+2*mu)*e11+lam*e22; sig22=lam*e11+(lam+2*mu)*e22; sig12=mu*e12;
clear e11 e22 e12
% area of patches
arp=full(sparse(it1,1,ar,n,1)+sparse(it2,1,ar,n,1)+sparse(it3,1,ar,n,1));
% mean value of stresses on patches
sm1=ar.*sig11; sm2=ar.*sig22; sm12=ar.*sig12;
s1=full(sparse(it1,1,sm1,n,1)+sparse(it2,1,sm1,n,1)+sparse(it3,1,sm1,n,1));
s2=full(sparse(it1,1,sm2,n,1)+sparse(it2,1,sm2,n,1)+sparse(it3,1,sm2,n,1));
s12=full(sparse(it1,1,sm12,n,1)+sparse(it2,1,sm12,n,1)+sparse(it3,1,sm12,n,1));
s1=s1./arp; s2=s2./arp; s12=s12./arp;
% Shear energy density
Sed=((.5+mu*mu/(6*(mu+lam)^2))*(s1+s2).^2+2*(s12.^2-s1.*s2))/(4*mu);
% Von Mises effective stress
delta=sqrt((s1-s2).^2+4*s12.^2); sp1=s1+s2+delta; sp2=s1+s2-delta;
Vms=sqrt(sp1.^2+sp2.^2-sp1.*sp2);
7.2
A membrane problem
14
NUMERICAL EXPERIMENTS
10
15
20
25
30
35
30
35
10
15
20
25
NUMERICAL EXPERIMENTS
Appendix
Main program for the L-shape problem
%------ Two-Dimensional linear elasticity
%
Example 1: L-Shape problem
%--------------------------------------% Elastic constants
E=1e5; nu=0.3;
lam=E*nu/((1+nu)*(1-2*nu)); mu=E/(2*(1+nu));
%
Penalty=10^20;
% Mesh
load lshape231 p t ibcd
n=size(p,2); nn=2*n;
% Exact solution
ue=zeros(nn,1);
[th,r]=cart2pol(p(1,:),p(2,:));
alpha=0.544483737; omg=3*pi/4; ralpha=r.^alpha/(2*mu);
C2=2*(lam+2*mu)/(lam+mu); C1=-cos((alpha+1)*omg)/cos((alpha-1)*omg);
ur=ralpha.*(-(alpha+1)*cos((alpha+1)*th)+(C2-alpha-1)*C1*cos((alpha-1)*th));
ut=ralpha.*( (alpha+1)*sin((alpha+1)*th)+(C2+alpha-1)*C1*sin((alpha-1)*th));
ue(1:2:end)=ur.*cos(th)-ut.*sin(th);
ue(2:2:end)=ur.*sin(th)+ut.*cos(th);
% Boundary conditions
nnb=2*length(ibcd);
ibc=zeros(nnb,1); ibc(1:2:end)=2*ibcd-1; ibc(2:2:end)=2*ibcd;
ubc=zeros(nnb,1); ubc(1:2:end)=ue(2*ibcd-1); ubc(2:2:end)=ue(2*ibcd);
% Assembly of the Stiffness matrix
K=elas2dmat1(p,t,E,nu);
% Right-hand side
f=zeros(nn,1);
% Penalization of Stiffness matrix and right-hand sides
K(ibc,ibc)=K(ibc,ibc)+Penalty*speye(nnb);
f(ibc)=Penalty*ubc;
% Solution by Gaussian elimination
u=K\f;
% Shear energy density & Von Mises effective stress
[Sh,Vms]=elas2dsvms(p,t,u,E,nu);
% Show the deformed mesh and shear energy density
elas2dshow(p,t,u,3000,Sh)
15
REFERENCES
16
Penalty=10^15;
% Mesh
load exple2mesh995 p t ibcd ibcneum
n=size(p,2); nn=2*n;
% Boundary conditions
nnb=2*length(ibcd);
ibc=zeros(nnb,1); ibc(1:2:end)=2*ibcd-1; ibc(2:2:end)=2*ibcd;
% Assembly of the Stiffness matrix
K=elas2dmat2(p,t,E,nu);
% Right-hand side: body forces
f1=zeros(n,1); f2=-.75*ones(n,1);
f=elas2drhs1(p,t,f1,f2);
% Right-hand side: Neumann forces
ibcn1=ibcneum(1,:); ibcn2=ibcneum(2,:); ibcn=union(ibcn1,ibcn2);
g1=zeros(n,1); g2=zeros(n,1); g2(ibcn)=10;
g=elas2drhs2(p,ibcneum,g1,g2);
clear ibcn ibcn1 ibcn2
% Penalization of Stiffness matrix and right-hand sides
K(ibc,ibc)=K(ibc,ibc)+Penalty*speye(nnb);
b=f+g; b(ibc)=0;
% Solution by Gaussian elimination
u=K\b;
% Shear energy density & Von Mises effective stress
[Sh,Vms]=elas2dsvms(p,t,u,E,nu);
%
% Show the deformed mesh and shear energy density
elas2dshow(p,t,u,20,Vms)
References
[1] Alberty J., Carstensen C. and Funken S.A. Remarks around 50 lines of matlab:
short finite element implementation. Numer. Algorithms, 20:117137, 1999.
[2] Alberty J., Carstensen C., Funken S.A. and Klose R. Matlab implementation
of the finite element method in elasticity. Computing, 69:239263, 2002.
[3] Ciarlet P.-G. The Finite Element Method for Elliptic Problems. North-Holland,
Amsterdam, 1979.
[4] Langemyr L. et al. Partial Differential Equations Toolbox Users Guide. The Math
Works, Inc., 1995.
[5] Kwon Y.W. and Bang H. The Finite Element Method Using MATLAB. CRC Press,
New York, 2000.
[6] Persson P.-O. and Strang G. A simple mesh generator in Matlab. SIAM Rev.,
42:329345, 2004.