Loading

We consider the same example as in Lesson 2 with the following modifications:

- We solve the linear system using the Conjugate Gradient method.
- In view of an iterative method we prefer to use, to handle boundary conditions, a classical substitution method rather than the penalty formulation.
- We use a defined material by giving its name in the mesh file.

The Finite Element Code

- We start like in Lesson 2.
#include "OFELI.h" #include "Therm.h" using namespace OFELI; int main(int argc, char *argv[]) {

- We expand the argument of the program:
if (argc <= 1) { cout << " Usage: lesson3 <mesh_file>" << endl; exit(1); }

- We read mesh data using an option that removes imposed degrees of
freedom for the list of unknowns. This is in contrast with the penalty technique used in
Lesson 2.
Mesh ms(argv[1],true);

- We store the matrix in a sparse format, thus using class SpMatrix<double>
SpMatrix<double> A(ms);

- Vectors b and x will store
respectively the right-hand side and the solution. Note that, since imposed degrees of freedom
are eliminated from the equations, the vectors have as sizes the actual number of equations.
The vector bc will store imposed (Dirichlet) boundary conditions,
and u is the final vector containing the
solution. Note that the vectors bc and u
are constructed using the mesh instance. This means that their size is
the total number of degrees of freedom, while the other vectors have
as size ms.getNbEq() which the number of equations,
*i.e.*the number of degrees of freedom minus the number of prescribed boundary conditions.Vect<double> b(ms.getNbEq()), x(ms.getNbEq()), bc(ms), u(ms);

- Let us present a simple method to implement boundary conditions: We use for this
the member function setNodeBC of class Vect<double>.
This function enables imposing by means of an
algebraic expression. Here, we impose on all nodes with code 1 the
value y, which is the y-coordinate of the node.
bc.setNodeBC(1,"y");

- We construct the linear system of equations just as in Lesson 2. The difference
here is that element right-hand sides need to be updated to take into account imposed boundary conditions at element
level. This is necessary when using an elimination technique for boundary conditions. The member function
updateBC is then used before assembly. Note that we have used a different member function for assembly than the
one in Lesson 2: Here this one is a member function of the equation class. Its use is simpler
since it uses as a unique argument either the matrix (for matrix assembly) or the right-hand side.
MeshElements(ms) { DC2DT3 eq(theElement); eq.Diffusion(); eq.updateBC(bc); eq.ElementAssembly(A); eq.ElementAssembly(b); }

- In order to use an iterative method to solve the linear system
use the class LinearSolver that
drives iterative solvers. Here we use the constructor with the maximal
number of iterations 1000, the tolerance for convergence criterium
1.e-8 and a verbosity parameter. We solve the
linear system using the conjugate gradient method preconditioned with the
diagonal ILU preconditioner.
LinearSolver<double> ls(1000,1.e-8,1); int nb_it = ls.solve(A,b,x,CG_SOLVER,DILU_PREC);

Note that the member function solve returns the number of performed iterations.

- We print out this number of iterations
cout << "Nb. of iterations: " << nb_it << endl;

- We incorporate boundary conditions into the solution vector:
u.insertBC(ms,x,bc);

- We finally output the solution vector and end the program
cout << "\nSolution:\n" << u; return 0; }

This is very simple: in the mesh data file, all elements have, in the present example, the code 1. If we say nothing, then a generic material (precisely called GenericMaterial with default properties is used. Otherwise, we can assign, in the mesh file a material, here Aluminium to this code, by the line

Material 1 AluminiumThis line must be given after all elements lines. Note that the file Aluminium.md must be present in the material's directory. We are now ready to test the package.

A test

We use here a finer mesh than in Lesson 2 except the line defining the material. We obtain the same solution as Lesson 2 after 15 iterations.

Copyright © 1998-2016 Rachid Touzani