Numme Lab: Finite Elements
Numme Lab: Finite Elements
Numme Lab: Finite Elements
Finite Elements
The objective of this practical work is to be able to program 2D finite elements using
python. We restrict ourselves to stationary thermal problems.
Reports : They have to be done by groups of two, and the next rules have to be followed :
• The work has to be submitted on the course website (one submission for each group).
• The report has to be a pdf file and an archive containing your scripts.
1 Introduction
The objective of this project is to write a small finite element code in python, which is able to solve
thermal analysis problems. The features of the code are the following:
• We assume 2D domains
• The elements are 3 nodes linear triangular elements formulated in the physical space.
A program skeleton is given in the file solveFE.py: you have to adapt this file for the project. This
file together with necessary python and mesh files are given in an archive named NUMMEFELab2019-2020.zip
that you can uncompress in your work directory.
• testFE.py: used to formulate your problem: mesh, physical properties and boundary con-
ditions. The formulation is given to solveFE for solving;
• solveFE.py: the core of the finite element program that reads the input, build the finite
element linear system, solve it and save the results;
• tri3ThermalDirect.py: the file where you have to code the missing elementary functions;
1
• export.py: contains functions used to export the results;
• gmshParser.py: contains a mesh class and functions to read mesh files (Gmsh file format)
and interact with the mesh. If needed, you would mainly have to interact with this class
through methods getElementConnectivityAndPhysNumber, getVertexCoordinates and
getNumVertices;
• integrationRule.py: contains Gauss quadrature rules for triangles and quadrangles (for-
mally, not needed here);
• sparseUtils.py: contains helpers that are used if the sparse version of solveFE is selected;
• Additional .msh and .geo files that you can use to test your program. .geo files contain the
definition of the geometry and .msh files the associated mesh. These files can be opened
with Gmsh. If you want to test your program with another number of elements, just edit the
.geo file (parameter nbelem), open it and click on Mesh/2D then Save mesh (format msh
version 2 only!).
3 Problem description
The problem will be specified by referring to geometric entities defined on the mesh. Any entity
in the mesh belongs to a geometrical entity than can be a point, a curve or a surface. All the
geometries given in the archive obeys the following convention:
Node 4 Line 103 Node 3
Line 104
Line 102
Surface 1000
If you want to prescribe a boundary condition on the elements lying on the bottom line of the
square, you will have to refer to line 101. Adding a source term on the domain ends-up to refer to
surface 1000.
4 Program usage
You can run the finite element code by specifying the problem you want to solve:
• The conductivity;
2
• The source term;
#Mesh file
meshName = 'square.msh'
exportName = 'temperature.pos'
The post-processing of the results is done thanks to gmsh. Open gmsh, then open your result file
file (.pos).
3
5 Work
You are asked to :
• Finish the implementation of the proposed finite element code. Look at solveFE.py: the
code is commented and some guidance is given;
• Basically, you need to fill the source file where # TODO comments are given;
– Compute the elementary stiffness matrix (implement your code in tri3ThermalDirect.py);
– Compute the elementary volume force vector (implement your code in tri3Thermal-
Direct.py);
– Compute the elementary Neumann force vector (implement your code in tri3Thermal-
Direct.py);
• You can compare your results to analytical solutions (for instance the case of a square with
zero temperature on line 101, and a prescribed flux on line 103)
• Present your results for u = 0 on line 101, r = 0, unit flux on line 103, and conductivity equal
to X X X . With X X X = groupNumber + Y Y Y and Y Y Y = 0 for group A, 100 for group B and
200 for group C. Example: X X X = 205 for group C-05.
6 Bonus questions
Only if you have finished properly the other questions. Otherwise it may lead to malus points...
1. Modify your program for a domain with two isotropic materials (see squareBimat.msh);
7 Tips
• To check your elementary stiffness matrix, here is the stiffness matrix of a unique
p triangular
element with a straight angle, unit conductivity and whose edge length are 1, 1, 2:
1.0 −0.5 −0.5
[K e ] = −0.5 0.5 0.0
−0.5 0.0 0.5
Note that the terms can be ordered in a different way in your program (depending on the
connectivity of the element). Two triangular meshes are given with the source files : tri1x1.msh
and tri2x2.msh.
• For the computation of the volume force vector, note that the sum of the interpolation func-
tions is equal to 1., and that all their integrals lead to the same result.
4
8 Heating circuit
We consider a micro-heating device whose geometry is depicted figure 1. It is composed of an
electrically resistive layer of Nichrome deposited on a glass plate. A voltage is applied on the layer,
which causes Joule heating. We are interested in quality of the device in term of homogeneity of
the temperature on the bottom layer of the glass plate. The dimensions of the plate are 4mm ×
2.1mm × 0.05mm. The domain being very thin, we model the problem as 2D, and assume that the
Joule effect can be modelled as a constant "volume" heat source r = 1 W.mm−2 . It is assumed that
the temperature on the lateral boundaries of the plate is fixed to 293K.
1. Propose a converged solution for the temperature of the problem (5% accuracy);
Hint: Set the useSparse option to True in order to decrease memory consumption.