Loading

The present lesson demonstrates how to use an optimization problem solver. We solve the same problem as in Lesson 2 as an optimization problem where Dirichlet boundary conditions are considered as equality constraints. The optimization algorithm is the Truncated Newton algorithm described in function OptimTN. This is not the best method to solve a Laplace equation, but our purpose here is to learn how to use this class.

The Finite Element Code

- We start, as usual, by including required headers and naming the appropriate namespace.
We furthermore include the header file of the optimization definition
class called
Opt .#include "OFELI.h" #include "Opt.h" #include "User.h" using namespace OFELI;

- The program will have as argument the name of the parameter data file:
int main(int argc, char *argv[]) {

- We expand program arguments and declare an instance of class IPF for parameter
file:
if (argc <= 1) { cout << " Usage: lesson5 <parameter_file>" << endl; exit(1) } IPF data(argv[1]);

- The mesh data file is obtained from a member function of instance data.
Mesh ms(data.getMeshFile());

- We introduce boundary conditions through a user defined class (See Lesson 4).
User ud(ms);

Implementation of class User will be given later.

- n is the number of degrees of freedom.
int n = ms.getNbDOF();

- We declare some vectors. x is the solution vector, vectors low
and up will contain for each degree of freedom upper and lower bound respectively to prescribe and
pivot is a vector that will contains (after optimization) flags to indicate which constraint is
reached.
Vect<double> x(n),low(n),up(n); Vect<int> pivot(n);

- Dirichlet boundary conditions are taken into account as before :
Vect<double> bc(n); ud.setDBC(bc);

- We start now to properly define the optimization problem. For this, we declare an instance of
a class Opt that is defined separately. The constructor of this class invokes the mesh instance
and the instance ud that will be useful to transmit to the class body or boundary sources.
Opt theOpt(ms,ud);

- We also initialize the solution to zero :
x = 0;

- As said before, Dirichlet boundary conditions are considered here as equality constraints and are
then incorporated into vectors low and up via a utility function called BCAsConstraint contained in the
mentioned file OptimAux.h
BCAsConstraint(ms,bc,up,low);

- The function OptimTN can now be invoked to run the optimization algorithm.
OptimTN<Opt>(theOpt,x,low,up,pivot,100,1.e-12,1);

The reader can refer to the function OptimTN to understand the meaning of each argument of the function.

- If the algorithm has succeeded (which is the case in this example) we can output the solution
and close the program.
cout << "\nSolution :\n" << x; }

- We start by including files OFELI.h and Therm.h
since we are going to solve a heat transfer problem. We also invoke the namespace OFELI:
#include "OFELI.h" #include "Therm.h" #include "User.h" using namespace OFELI;

- This class will have a constructor that acquires the mesh and the user data instance and stores
pointers to these objects:
class Opt { public: Opt(Mesh &ms, User &ud) { _ms = &ms; _ud = &ud; }

- The other public member is the objective function.
void Objective(Vect<double> &x, double &f, Vect<double> &g) { f = 0.; g = 0.; MeshElements(*_ms) { DC2DT3 eq(theElement); Vect<double> xe(theElement,x); f += eq.Energy(xe); eq.EnerGrad(xe); eq.ElementAssembly(g); } }

Note that the arguments of the objective function are necessarily the optimization variable vector, the objective function to compute and the gradient vector to compute. Here the class DC2DT3 provides necessary material for optimization purposes. Namely, member function Energy return the energy value and function EnerGrad calculates the element energy gradient that needs be assembled into the global vector g.

- It remains to declare the mesh and user data pointers as private attributes of the class and end the class definition.
private: Mesh *_ms; User *_ud; };

A test

We use here exactly the same mesh file as in the previous lesson. Of course, we obtain the same solution as the previous lesson. The output of the program is the following:

NIT NF CG F GTG 0 1 0 6.00000000e+000 4.00000000e+001 1 2 3 1.67852990e+000 4.09059524e+001 2 3 5 1.50042385e+000 4.09106651e+001 3 4 8 1.50001179e+000 4.09109024e+001 4 5 10 1.50000062e+000 4.09109074e+001 5 6 13 1.50000000e+000 4.09109074e+001 6 7 15 1.50000000e+000 4.09109074e+001 7 8 17 1.50000000e+000 4.09109074e+001 8 9 19 1.50000000e+000 4.09109074e+001 Optimal Function Value = 1.5 Solution : 1 1.00000000e+000 2 7.49999986e-001 3 4.99999972e-001 4 2.49999994e-001 5 0.00000000e+000 6 1.00000000e+000 7 7.49999988e-001 8 4.99999982e-001 9 2.49999978e-001 10 0.00000000e+000 11 1.00000000e+000 12 7.49999989e-001 13 4.99999982e-001 14 2.49999978e-001 15 0.00000000e+000 16 1.00000000e+000 17 7.49999986e-001 18 4.99999972e-001 19 2.49999994e-001 20 0.00000000e+000

Copyright © 1998-2016 Rachid Touzani