Loading

We consider here a 2-D steady state boundary value problem. We solve a Poisson equation (Diffusion)
with "simple" data. Concerning boundary conditions, we impose a Dirichlet (essential) boundary condition on a portion of
the domain and a homogeneous Neumann (natural) condition on the remaining boundary. Note that owing to the variational
formulation the Neumann condition is implicit (*We have nothing to do for it*).

The Finite Element Code

Here is a description of the source code. The source file can be found Example 2 in the examples directory- As usual we start by including the principal header file and the file Therm.h
that includes all classes related to heat transfer problems.
#include "OFELI.h" #include "Therm.h" using namespace OFELI;

- Our program will have as argument the mesh file name.

int main(int argc, char *argv[]) {

- We output the OFELI banner and get the program argument
banner(); if (argc <= 1) { cout << " Usage: lesson2 <mesh_file>" << endl; exit(1); }

- We construct an instance of class Mesh by using the file that contains
mesh data in file. Note that to use fist a default constructor then invoke the member function
get to read mesh data in the mesh file.

Mesh ms(argv[1]);

- The problem matrix is symmetric and will be stored in skyline format (thus using class
SkSMatrix<double>)
SkSMatrix<double> A(ms);

- Vectors b and bc (as instances
of class Vect<double>) will contain the right-hand side and the solution, and imposed boundary conditions at nodes:
Vect<double> b(ms), bc(ms);

- We assign imposed boundary conditions to vector bc by
using member function setNodeBC which allows using an interpreted function of
node coordinates: We prescribe the function y to nodes with code 1.
bc.setNodeBC(1,"y");

- We now start building the linear system. For this, we have to implement a loop over all element
meshes by using member functions topElement() and getElement(). The returned
pointer theElement enables access to current element data. The OFELI library
defines the macro MeshElements(ms) that stands for a shorthand for the line
for (ms.topElement(); (theElement=ms.getElement());)

For each element, we construct an instance of class DC2DT3 for diffusion-convection problems in 2-D using 3-node triangles. The member function Diffusion calculates the contribution to element matrix of diffusion term. We then assemble matrix and right-hand side (which is 0 here).MeshElements(ms) { DC2DT3 eq(theElement); eq.Diffusion(); eq.ElementAssembly(A); eq.ElementAssembly(b); }

Of course, in the present case, assembling the right-hand side is useless.

- Once the linear system is assembled, Dirichlet boundary conditions are imposed by a penalty
technique. This is implemented via the function Prescribe, member of all matrix classes.
A.Prescribe(b,bc);

- Solution is obtained by factorizing and backsubstituting:
A.solve(b);

Vector b contains now the solution.

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

A finite element mesh

To test this program we use a finite element mesh of a rectangle [0,3]x[0,1]. The imposed boundary conditions are

u(x,0)=0, u(x,1)=1, 0<x<3Homogeneous Neumann boundary conditions are "imposed" on the portions x=0 and x=3. The solution is then

u(x,y)=yThe mesh file is test.m. Note, in this file, that a code equal to 1 is associated to nodes with y=0 and y=1.

**You can now execute the code to obtain the exact solution.**

Copyright © 1998-2016 Rachid Touzani