Vectorized Matlab codes for linear two-dimensional elasticity

. A vectorized Matlab implementation for the linear ﬁnite element is provided for the two-dimensional linear elasticity with mixed boundary conditions. Vectorization means that there is no loop over triangles. Numerical experiments show that our implementation is more efﬁcient than the standard implementation with a loop over all triangles.


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 the best performance in large scale 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.[1,2], we propose a P 1 -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
We consider an elastic body which occupies, in its reference configuration, a bounded domain Ω in R 2 with a boundary Γ.Let {Γ D , Γ N } be a partition of Γ with Γ N possibly empty.We assume Dirichlet conditions on Γ D and Neumann conditions on Γ N .Let u = (u 1 , u 2 ) be the two-dimensional displacement field of the elastic body.Under the small deformations assumption, constitutive equations are σ ij (u) = 2µ ij (u) + λ tr ( (u))I 2 , i,j = 1, 2, (1) (u) = (∇u + ∇u T )/2, where λ and µ denote Lam é (positive) constants.These coefficients are related (in plane deformations) to the Young modulus E and the Poisson coefficient ν by and u D , the problem studied in this paper can be formulated as follows Let us introduce subspaces The variational formulation of Eqs ( 2)-( 4) is In Eq. ( 5), C = (C ijkl ) is the fourth-order elastic moduli tensor corresponding to Eq. ( 1), i.e. where δ ij being the Kronecker delta.

Finite element discretization
Let T h 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 V D h defined by is the space of polynomials of degree less or equals to 1 on the triangle T .The discrete version of Eq. ( 5) is then Let {φ j } be the system of piecewise global basis functions of V h , i.e. for all , where u j α are nodal values of u h , i.e. u αh (x j ) = u j α .Applying the standard Galerkin method to Eq. ( 6) yields The stiffness matrix A = (A ij ) and the right-hand side b = (b i ) are then given by The stiffness matrix A is sparse, symmetric and positive semi-definite before incorporating the boundary condition In practice, integrals in Eq. ( 6) are computed as sums of integrals over all triangles, using the fact that Let {ϕ i } be the linear basis functions of the triangle T or the edge E. If we set then assembling operations consist of

Data representation of the triangulation
For the mesh, we adopt the data representation used in Matlab PDE Toolbox [6].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 Fig. 1, the nodes array p is (np=9) 0.0000 1.0000 1.0000 0. 0.5000 1.0000 0.5000 0.0000 0.5000 1.0000 1.0000 0.0000 0. 1.0000 0.5000 0.0000 0.5000 0.5000 and the elements array t is (nt=8)

Assembling the stiffness matrix
As in [2,5] we adopt Voigt's representation of the linear strain tensor

Using this representation, the constitutive equation σ(u
where ) t be the vector of nodal values of u h on a triangle T .Elementary calculations provide, on a triangle T , where |T | is the area of the triangle T .From Eqs (11) and ( 12), it follows that The element stiffness matrix is therefore The element stiffness matrix Eq. ( 13) can be computed simultaneously for all indices using the fact that Matlab implementation given in [2,5] (in a more modular form) can be summarized by the Matlab Function 5.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.

Assembling the right-hand side
We assume that the volume forces f = (f 1 , f 2 ) are provided at mesh nodes.The integral Eq. ( 10) is approximated as follows where (x c , y c ) is the center of mass of the triangle T .With the assumption on f , f (x c , y c ) = (f 1 (x c , y c ), f 2 (x c , y c )) with where {(x i , y i )} i=1,3 are vertices of the triangle T .Using the notation convention Eq. ( 14), we have We can compute Eqs (18) and (19) 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 3.
In all examples, Dirichlet boundary conditions are taken into account by using large spring constants (i.e.penalization).

The exact solution is known in polar coordinates
The exponent α is the solution of the equation with ω = 3π/4 and  show informations recorded with a mesh with 13041 nodes.It appears that the assembling function elas2dmat1 takes about 98% of CPU time.Note that the time given in Table 3 for the assembling operations is the time taken by intruction K=elas2dmat1(p,t,E,nu) while the time given in Table 2 is only the time spent in the assembling function (without calling and returning operations).Table 1 shows clearly that the assembling function elas2dmat1 is the bottleneck of the program.We now compare the performances of Matlab functions elas2dmat1 and elas2dmat2 for assembling the stiffness matrix of the L-shape problem.To this end, various meshes of the L-shape are generated and we use Matlab command profile to compute the elapsed time.Table 2 shows the performances of assembling Functions 1-2.One can notice that the saving of computational time, with the vectorized function elas2dmat2, is considerable.Table 3 shows that the finite element program with elas2dmat2 is more balanced compared to elas2dmat1.
We report in Table 4 the H 1 and L 2 distances between the exact solution Eqs (20)-( 21) and the approximate solution using the assembling function elas2dmat2.The distances are computed using a 13-point Gaussian quadrature.One can notice that u − u h L 2 (Ω) → 0 and u − u h H 1 (Ω) → 0 has the mesh size h goes to zero but the convergence rates are lower than theoretical ones (2 for the L 2 -error and 1 for the H 2 -error for elliptic problems,  see e.g.[3,4]).The reason is that the exact solution u does not belong to H 2 (Ω) so that the error estimate theorems do not hold.Figure 4 shows the deformed mesh (231 nodes and 400 triangles) of the L-shape with a displacement field multiplied by a factor 3000.The grey tones visualize the approximate shear energy density and show the singularity at (0, 0).

Conclusion
We have demonstrated that, for solving isotropic linear elasticity problems in Matlab with the finite element method, the vectorized code is much more efficient than a standard implementation with a loop over triangles.
Further work is underway to derive vectorized codes for the three-dimensional linear elasticity using tetrahedral elements.The main difficulty is to invert analytically, with easy to implement formulas, the matrix  where {(x i , y i , z i )} i=1,...,4 are the vertices of a tetrahedron.Indeed, if {ϕ i } i are the linear basis functions on a tetrahedron, their gradients are given by (see e.g.[2])

Function 2 .
matrices Eqs (15)-(17) 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 2. Assembly of the stiffiness matrix with a vectorized code function K=elas2dmat2(p,t,Young,nu) %-

Fig. 3 .
Fig. 3. Matlab profile command report details with elas2dmat1 assembling function using a mesh with 13041 nodes.

Table 1
Percentage of computing effort taken by elas2dmat1

Table 3
Percentage of computing effort taken by elas2dmat2 Fig. 2. Matlab profile command report with elas2dmat1 assembling function using a mesh with 13041 nodes.

Table 4 L
2and H 1 errors and convergence rates