Using Multivac

This section details the way a program using the library Multivac is be written. The main idea is to gather a few objects to describe the simulation to be performed by Multivac. Of course there is an object associated with the speed of the front, but there are objects for the initial curve, for the output saver, for the numerical scheme, etc. Thanks to the following explanations is fairly easy to understand the structure and to choose the right objects for a given simulation, even for non-experts in level set methods.

The structure

A typical program using Multivac may be split into several sections:

Includes and options: includes Multivac and sets its overall options;
Time: selects the simulation length and its time step;
Domain and mesh: describes the mesh;
Speed function: provides the speed function (the speed in the normal direction to the front);
Initial curve: the initial front;
Level set function: the type of level set function to be used, which may depend on the level method used;
Initializer: the algorithm used to initialize or to reinitialize the level set function and the velocity field;
Updater: the numerical scheme used to update the level set function;
Saver: what should be saved;
Simulator: the object that handles the simulation.

A good example is the program track.cpp which is fully quoted at the end of this section. track.cpp is distributed with Multivac (in the top directory) with several other examples (all files with the "cpp" extension).

track.cpp

Includes and options

The first lines in track.cpp are:

/************************
 * INCLUDES AND OPTIONS *
 ************************/

// Seldon library and Multivac project provide error management
// using exception handling.  Exceptions that may be raised
// are selected according to debug levels.
#define MULTIVAC_DEBUG_LEVEL_2

// Define MULTIVAC_REPORT if you want to clock the time needed for
// updates and initializations (results will be displayed on screen).
#define MULTIVAC_REPORT_NO

// Multivac includes.
#include "multivac.hxx"

using namespace Multivac;

In this part, the library Multivac is included (through multivac.hxx). Before that, there are two important lines that determine the behavior of Multivac: #define MULTIVAC_DEBUG_LEVEL_2 and #define MULTIVAC_REPORT_NO.

MULTIVAC_DEBUG_LEVEL_2 selects the debug mode for Multivac. It determines how many checks (at run time) are to be performed. Of course the number of checks has an impact on the CPU time needed for the simulation. This is the reason why the user should choose the right level. There are four levels MULTIVAC_DEBUG_LEVEL_X where X ranges from 1 to 4.

During the development, it is advocated to make as many checks as possible and to choose the fourth level (i.e. to put: #define MULTIVAC_DEBUG_LEVEL_4). This way, Multivac will essentially check for the memory allocations, for the input/output operations and the indices validity.

If the simulation takes time and if Multivac is responsible for this -- not the computations for the speed function --, the final program should run in the second mode, i.e. with #define MULTIVAC_DEBUG_LEVEL_2, as in track.cpp.

The advanced user might want to know that, if MULTIVAC_DEBUG_LEVEL_X is defined, then SELDON_DEBUG_LEVEL_X is defined by Multivac. Actually all the checks, that are implied by MULTIVAC_DEBUG_LEVEL_X, are those made in the underlying C++ library Seldon (which provides the matrices and vectors used in Multivac). Therefore to get details about the debug mode, the advanced user is invited to read the documentation of the C++ library Seldon (homepage: http://spacetown.free.fr/lib/seldon/).

#define MULTIVAC_REPORT_NO may be replaced by #define MULTIVAC_REPORT to get a report of the (user) time spent in the different steps of the simulation.

The last option is the numerical accuracy of the simulation. It may be modified at the line:

  typedef double real;

In the previous line, real is set to double: Multivac will run in double-precision mode. Replace double by float to have Multivac running into single-precision mode. The mechanism is clear: most variables are declared with the type real in track.cpp and real is automatically replaced (at compile time) with its value set in the line typedef [user type] real;.

Time

The next part in track.cpp sets the time step and the length of the simulation:

  //////////
  // TIME //
  //////////

  // Final time of the simulation, the initial time being 0.
  real FinalTime = 0.1;
  // Time step.
  real Delta_t = 0.0001;

  // Number of iterations.
  int NbIterations = int (FinalTime / Delta_t);

FinalTime is of course the final time of the simulation. The initial time is assumed to be 0. Delta_t is the time step. The number of iterations is computed from the final time and the time step, but the user may directly set it. Actually, Multivac takes into account the final time and the number of iterations. The variable Delta_t may therefore be omitted.

Domain and mesh

Then the domain is defined and the mesh associated with it:

  ///////////////////
  // DOMAIN & MESH //
  ///////////////////

  // Choose the type of the mesh:
  //   1) COrthogonalMesh<real>
  //   2) --
  typedef COrthogonalMesh<real> MeshType;

  // Domain bounds (the domain is a rectangle).
  real Xmin = 0.0;
  real Xmax = 3.0;
  real Ymin = 0.0;
  real Ymax = 3.0;

  // For an orthogonal mesh, Nx and Ny are the number of
  // grid points along (x'x) and (y'y) (respectively).
  int Nx = 301;
  int Ny = Nx;

  // Choose the right constructor:
  //   1) Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny)
  //   2) --
  MeshType Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny);

Xmin, Xmax, Ymin and Ymax define the domain. The discretization is an orthogonal mesh with Nx points along x and Ny points along y. Here, Nx is set to Ny but it is not compulsory.

Notice the structure of the section: (1) the user first chooses a type for the mesh, (2) he then sets several variables and (3) he finally defines an object of the type previously chosen (MeshType = COrthogonalMesh<real>) with the relevant variables. The structure is the same in most sections.

There is no choice for the mesh: it must be COrthogonalMesh<real>. But if Multivac had managed unstructured meshes, then this option would have been proposed. The constructor is determined by the type previously chosen. In the last line, one has to choose the right constructor (Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny)), i.e. the constructor associated with the type that was chosen. In every part, the same is true: if one chooses the type i) at the beginning of the section, then one must also use the constructor i) at the end.

Speed function

This section is fundamental since it sets the velocity in the normal direction to the front.

  ////////////////////
  // SPEED FUNCTION //
  ////////////////////

  // Choose the speed function:
  //   1) CConstantSpeed<real>
  //   2) CPiecewiseConstantSpeed<real>
  //   3) CFireModel<real>
  //   4) CSimplifiedFireModel<real>
  typedef CSimplifiedFireModel<real> SpeedType;

  // Speed rate for a constant speed function.
  real SpeedRate = 0.5;

  // Second speed rate for a piecewise-constant speed function.
  real SpeedRate0 = 0.2;

  // Limit (for a piecewise-constant speed function).
  real Limit = 0.9;

  // Parameters for both simplified and full fire model.
  real U = 100.0;
  real m = 1.5;
  real c1 = 0.5;
  real epsilon0 = 0.2;

  // Last parameter for the simplified fire model.
  real alpha = 0.9;

  // Parameters for the full fire model.
  real a = 0.1;
  real b = 1.0;
  real epsilon1 = 0.003;

  // Choose the right constructor:
  //   1) F(SpeedRate)
  //   2) F(SpeedRate, SpeedRate0, Limit)
  //   3) F(U, m, c1, epsilon0, a, b, epsilon1)
  //   4) F(U, m, c1, epsilon0, alpha)
  SpeedType F(U, m, c1, epsilon0, alpha);

The user probably wants to implement his own velocity. Details about the inclusion of a new speed function are provided in the section "Adding a speed function". However the user should first understand the way it works with the speed functions that are already provided.

The same structure as for the mesh is shown. First the user selects the type of his speed function (e.g. CSimplifiedFireModel<real>); then he sets a few parameters (e.g. the wind speed U = 100.0); finally the speed-function object is declared with the right constructor. In track.cpp, the fourth type was chosen, therefore the fourth constructor must be used.

Between typedef CSimplifiedFireModel<real> SpeedType; and the constructor SpeedType F(U, m, c1, epsilon0, alpha), a few parameters are set. The only useful parameters are the variables that are arguments in the constructor, i.e. (in this case) U, m, c1, epsilon0 and alpha. Of course, one may have called the constructor directly with the numerical values: SpeedType F(100.0, 1.5, 0.5, 0.2, 0.9).

Four speed functions are provided as examples. The first one is a constant velocity. Its constructor therefore takes a single argument: F(SpeedRate). The second one is a piecewise constant speed function. The domain is split horizontally at y = Limit (= 0.9, in track.cpp), and a constant velocity is associated with each part of the domain (SpeedRate above and SpeedRate0 below). The last speed functions are more complex velocities associated with forest fire propagations. No detail is provided because the description of the fire models are out of the scope of this guide.

With regard to the management of the speed function, the declaration of the speed function (the call to its constructor) is the only requirement in the program. After this declaration, all is managed by Multivac.

Initial curve

This section describes how the initial curve (the front at the initial time) is defined.

  ///////////////////
  // INITIAL CURVE //
  ///////////////////

  // Choose the initial curve:
  //   1) CCircle<real>
  //   2) CTwoCircles<real>
  //   3) CThreeCircles<real>
  //   4) CIsland<real>
  //   5) CIsland0<real>
  //   6) CSetOfPoints<real>
  typedef CCircle<real> InitialCurveType;

  // For a circle, center coordinates and radius.
  real CircleCenterX = 1.0;
  real CircleCenterY = 1.5;
  real CircleRadius = 0.5;

  // For a second circle, center coordinates and radius.
  real CircleCenterX0 = 1.65;
  real CircleCenterY0 = 1.6;
  real CircleRadius0 = 0.3;

  // For a third circle, center coordinates and radius.
  real CircleCenterX1 = 0.50;
  real CircleCenterY1 = 1.0;
  real CircleRadius1 = 0.25;

  // File containing a front defined by a set of points.
  string InitialFrontFile = "[full path]/[set].pts";
  // Is the (Xmin, Ymin) outside the front?
  bool origin_out = false;

  // If set to 'true', outside and inside are swapped.
  bool reversed = false;

  // Choose the right constructor:
  //   1) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                   reversed)
  //   2, 4) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                      CircleCenterX0, CircleCenterY0, CircleRadius0,
  //                      reversed)
  //   3, 5) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                      CircleCenterX0, CircleCenterY0, CircleRadius0,
  //                      CircleCenterX1, CircleCenterY1, CircleRadius1,
  //                      reversed)
  //   6) InitialCurve(InitialFrontFile, 0, origin_out ^ reversed)
  InitialCurveType InitialCurve(CircleCenterX, CircleCenterY, CircleRadius);

The initial front is a closed curve or, if not closed, a curve whose ends are outside the simulation domain. The normal to the front is directed to the outside of the front. Hence, if the speed function is positive, the front expands in the outside (of the front itself), and it shrinks otherwise.

The following initial fronts may be defined:

CCircle: a single circle.
CTwoCircles: two circles (that may overlap).
CThreeCircles: three circles (that may overlap).
CIsland: creates two circles. This is the same as CTwoCircles, but the center of the second circle is assumed to be outside the front. The second circle should therefore be inside the first one and it is referred as an "island".
CIsland0: three circles. The third circle is an island.
CSetOfPoints: the initial curve is defined by a set of points. One puts a sorted list of points in a text file and Multivac assumes that those points are linked by segments (the first point in the list is linked to the second point which in turn is linked to the third point, etc.). If the curve is closed, the last point of the list must be the same as the first point. Furthermore one chooses whether the origin of the domain (for the simulation) is inside or outside the front. Notice that only a single curve may be described: it is not possible to describe two disconnected curves.

M.cpp, provided with Multivac, provides an example for CSetOfPoints. It simulates the letter M that shrinks under a constant speed function. The file M.pts describes the shape of the letter M:

0.5 0.2
0.5 2.8
0.9 2.8
1.5 2.0
2.1 2.8
2.5 2.8
2.5 0.2
2.1 0.2
2.1 1.6
1.5 1.2
0.9 1.6
0.9 0.2
0.5 0.2

Notice the parameter reversed that enables to swap the inside and the outside of the initial front.

Level set function

  ////////////////////////
  // LEVEL SET FUNCTION //
  ////////////////////////

  // Choose the level set function type:
  //   1) COrthogonalLevelSet<real>
  //   2) --
  typedef COrthogonalLevelSet<real> LevelSetType;

  // Choose the right constructor:
  //   1) Phi
  //   2) --
  LevelSetType Phi;

There is no choice. Easy!

Initializer

The initializer is an important component which is in charge of the initializations and the re-initializations of the level set function and the velocity field.

  /////////////////
  // INITIALIZER //
  /////////////////

  // Choose an initializer:
  //   1) CNarrowBandNeverInit<real> *
  //   2) CNarrowBandExtension<real> *
  //   3) CFastMarchingNeverInit<real> **
  //         * Narrow band level set method
  //        ** Fast marching method
  typedef CNarrowBandNeverInit<real> InitializerType;

  InitializerType Initializer;

Multivac is able to solve the unstationary and the stationary problems. Solving the unstationary problem means evolving in time the front. Solving the stationary problem means computing the arrival time of the front at any mesh point. The underlying equations are Hamilton-Jacobi equations, one with initial conditions and the other with boundary conditions. Multivac solves the unstationary problem with the so-called narrow band level set method and the stationary problem with the fast marching method. The stationary problem should only be used in very particular problems, but it is solved much faster (an example is provided with Multivac: M_fast.cpp). Usually one solves the unstationary problem because of the dependencies of the speed function. To get details, one may refer, for example, to Level Set Methods and Fast Marching Methods, J. A. Sethian, Cambridge University Press.

If one solves the stationary problem, one must choose CFastMarchingNeverInit<real>. For the unstationary problem, there are two possibilities.

CNarrowBandNeverInit<real> is a simple initializer. It demands that the speed function is defined everywhere in the mesh, not only on the front itself. Indeed, during the initializations, it computes the speed function in the whole narrow band, not only on the front. If the speed function only depends on the position and the time, this initializer is fine. If the velocity depends on the properties of the front such as the curvature, this initializer should be discarded. It may still be used if the velocity depends on the normal to the front because it would use the normal to the local level set, which is a fair approximation. Last but not the least, if computing the velocity has a high computational cost, this initializer should be discarded (because it computes the velocity in the whole narrow band).

In practice, one should switch to CNarrowBandExtension<real> if the velocity depends on the curvature, if the velocity can only be computed on the front, or if computing the velocity implies high computational costs.

Updater

The updater handles the numerical scheme, i.e. it updates the level set function at every time step. It also handles the properties of the narrow band (for unstationary problems).

  /////////////
  // UPDATER //
  /////////////

  // Choose an initializer:
  //   1) CNarrowBandFirstOrderEngquistOsher<real> *
  //   2) CNarrowBandFirstOrderLaxFriedrichs<real> *
  //   3) CNarrowBandEno2EngquistOsher<real> *
  //   4) CFastMarchingFirstOrderEngquistOsher<real> **
  //         * Narrow band level set method
  //        ** Fast marching method
  typedef CNarrowBandFirstOrderEngquistOsher<real> UpdaterType;

  // Choose the right constructor and choose its parameters:
  //   1, 2, 3) Updater(TubeSemiWidth, BarrierWidth, OutSpaceWidth)
  //              | TubeSemiWidth: number of grid points on each side of the front.
  //              | BarrierWidth: number of grid points (on each side of the front)
  //                        that imply tube reconstruction when reached.
  //              | OutSpaceWidth: number of grid points (on each side of the front)
  //                         that must not be reached.
  //              | Examples: Updater(6, 3, 1) or Updater(12, 5, 1).
  //   4) Updater(TMax)
  //        | TMax: time greater than all arrival times that will be computed.
  UpdaterType Updater(6, 3, 1);

For stationary problems, one must choose CFastMarchingFirstOrderEngquistOsher<real> which uses a first-order numerical scheme (Engquist-Osher). One has to provide a bound TMax for the arrival times at mesh points. The simulation will stop as soon as a computed arrival-time is greater than TMax or the arrival times at every mesh point are computed.

For unstationary problems, three numerical schemes are available including a second-order scheme (CNarrowBandEno2EngquistOsher). Details about those numerical schemes may be easily found in the literature dedicated to the level set methods. The user also sets the widths associated with the narrow band. The band -- also called tube -- may be split into three zones (the front is in bold, in the middle):

The widths of the three zones are:

TubeSemiWidth: the number of grid points on each side of the front (zones 1 + 2 + 3).
BarrierWidth: the number of grid points on each side of the front in the zones 2 and 3.
OutSpaceWidth: the number of grid points on each side of the front in the zone 3.

Below is a summary (where F stands for the front, T stands for TubeSemiWidth, B for BarrierWidth and O for OutSpaceWidth):

F
F       (zone 1)           |   (zone 2)    |  (zone 3)  |   (outside the band)
F       (zone 1)           |   (zone 2)    <-----O----->|
F       (zone 1)           <-------------B------------->|
F<--------------------------T-------------------------->|
F

The user has to choose those values carefully. The widths have an impact on the accuracy and on the computational cost. If the tube is too large, the computations take place in wide area, which increases the computational cost. If the tube is too narrow, the computational cost is also high because the tube will be often rebuilt (to keep the front inside). To give a rough idea, the tube semi-width may be set in [6, 12]. BarrierWidth may range from 3 to 6. OutSpaceWidth should be set to 1.

Saver

  ///////////
  // SAVER //
  ///////////

  // Choose the saver type:
  //   1) CNeverSave<real>
  //        | Nothing is saved.
  //   2) CCurvesSaver<real>
  //        | The front is constructed and saved at each time step.
  //        | Not relevant for the fast marching method.
  //   3) CSaveLastCurve<real>
  //        | Saves the last curve, after all calculations.
  //        | Not relevant for the fast marching method.
  //   4) CSaveAtTheEnd<real>
  //        | Saves the level set, after all calculations.
  //        | Relevant for the fast marching method.
  typedef CCurvesSaver<real> SaverType;

  // Number of curves that will be saved.
  // Set NbCurves to 0 in order to save all curves.
  int NbCurves = 10;

  // Output directory.
  string Directory = "/home/vivien/Computations/Fronts/";

  // Stores time at each time step.
  string TimeFile = Directory + "Time";

  // Stores front points.
  string CurvesFile = Directory + "Curves";
  // Stores number of points on fronts at each time step.
  string CurveLengthsFile = Directory + "CurveLengths";
  // Stores the level set function.
  string PhiFile = Directory + "Phi";
  // Stores the speed function.
  string FFile = Directory + "F";

  // Stores grid coordinates along the (x'x) axe.
  string XFile = Directory + "X";
  // Stores grid coordinates along the (y'y) axe.
  string YFile = Directory + "Y";
  // Stores mesh points.
  string PointsFile = Directory + "Points";
  // Stores mesh edges.
  string EdgesFile = Directory + "Edges";
  // Stores mesh triangles.
  string TrianglesFile = Directory + "Triangles";

  // Indicates the number of iterations between saves.
  int Period = NbIterations / (NbCurves==0?NbIterations:NbCurves);

  SaverType Saver(TimeFile, CurvesFile,
		  CurveLengthsFile, PhiFile, FFile, XFile, YFile,
		  PointsFile, EdgesFile, TrianglesFile, Period);

For stationary problems, one should to choose CSaveAtTheEnd<real>, which simply saves the solution, namely the level set function.

For unstationary problems, there are three classes: CNeverSave<real>, CCurvesSaver<real> and CSaveLastCurve<real>. Their descriptions in the code above are clear. Notice that the constructor is the same for all types: SaverType Saver(TimeFile, CurvesFile, ...).

The output files are always text files. For unstationary problems, the fronts are stored in CurvesFile as a sorted list of points (which are the intersections between the edges of the mesh and the front), one point per line. There is no delimiter between two curves, may they be associated with the same time step or not. For each time step, the number of saved points is written in CurveLengthsFile, which enables to identify at which time step a given point was saved.

Simulation

  /////////////////////
  /////////////////////
  ////  SIMULATOR  ////
  /////////////////////
  /////////////////////

  CSimulator<real, MeshType, SpeedType, InitialCurveType,
    LevelSetType, InitializerType, UpdaterType, SaverType>
    Simulator(Mesh, F, InitialCurve, Phi,
	      Initializer, Updater, Saver,
	      NbIterations, FinalTime);


  ////////////////////////////////
  // SIMULATION INITIALIZATIONS //
  ////////////////////////////////

  Simulator.Init();


  ////////////////
  // SIMULATION //
  ////////////////

  Simulator.Run();

Nothing should be modified by the user. The simulation is initialized and is then launched.

Other examples

The user is welcome to browse the other examples, i.e. the files with the extension "cpp" in the top directory of Multivac. He should quickly be able to have his own simulations running. The next sections "Displaying the output" and "Adding a speed function" should also help.

track.cpp: full content

track.cpp is distributed with Multivac (in its top directory), with many other examples (files with the extension "cpp"). It is also provided in this documentation as a file and it is quoted below:

/************************
 * INCLUDES AND OPTIONS *
 ************************/

// Seldon library and Multivac project provide error management
// using exception handling.  Exceptions that may be raised
// are selected according to debug levels.
#define MULTIVAC_DEBUG_LEVEL_2

// Define MULTIVAC_REPORT if you want to clock the time needed for
// updates and initializations (results will be displayed on screen).
#define MULTIVAC_REPORT_NO

// Multivac includes.
#include "multivac.hxx"

using namespace Multivac;


/*****************
 * MAIN FUNCTION *
 *****************/

int main()
{

  // To catch exceptions.
  TRY;
  
  // real type: double, float...
  typedef double real;


  //////////
  // TIME //
  //////////

  // Final time of the simulation, the initial time being 0.
  real FinalTime = 0.1;
  // Time step.
  real Delta_t = 0.0001;

  // Number of iterations.
  int NbIterations = int (FinalTime / Delta_t);


  ///////////////////
  // DOMAIN & MESH //
  ///////////////////

  // Choose the type of the mesh:
  //   1) COrthogonalMesh<real>
  //   2) --
  typedef COrthogonalMesh<real> MeshType;

  // Domain bounds (the domain is a rectangle).
  real Xmin = 0.0;
  real Xmax = 3.0;
  real Ymin = 0.0;
  real Ymax = 3.0;

  // For an orthogonal mesh, Nx and Ny are the number of
  // grid points along (x'x) and (y'y) (respectively).
  int Nx = 301;
  int Ny = Nx;

  // Choose the right constructor:
  //   1) Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny)
  //   2) --
  MeshType Mesh(Xmin, Xmax, Ymin, Ymax, Nx, Ny);


  ////////////////////
  // SPEED FUNCTION //
  ////////////////////

  // Choose the speed function:
  //   1) CConstantSpeed<real>
  //   2) CPiecewiseConstantSpeed<real>
  //   3) CFireModel<real>
  //   4) CSimplifiedFireModel<real>
  typedef CSimplifiedFireModel<real> SpeedType;

  // Speed rate for a constant speed function.
  real SpeedRate = 0.5;

  // Second speed rate for a piecewise-constant speed function.
  real SpeedRate0 = 0.2;

  // Limit (for a piecewise-constant speed function).
  real Limit = 0.9;

  // Parameters for both simplified and full fire model.
  real U = 100.0;
  real m = 1.5;
  real c1 = 0.5;
  real epsilon0 = 0.2;

  // Last parameter for the simplified fire model.
  real alpha = 0.9;

  // Parameters for the full fire model.
  real a = 0.1;
  real b = 1.0;
  real epsilon1 = 0.003;

  // Choose the right constructor:
  //   1) F(SpeedRate)
  //   2) F(SpeedRate, SpeedRate0, Limit)
  //   3) F(U, m, c1, epsilon0, a, b, epsilon1)
  //   4) F(U, m, c1, epsilon0, alpha)
  SpeedType F(U, m, c1, epsilon0, alpha);


  ///////////////////
  // INITIAL CURVE //
  ///////////////////

  // Choose the initial curve:
  //   1) CCircle<real>
  //   2) CTwoCircles<real>
  //   3) CThreeCircles<real>
  //   4) CIsland<real>
  //   5) CIsland0<real>
  //   6) CSetOfPoints<real>
  typedef CCircle<real> InitialCurveType;

  // For a circle, center coordinates and radius.
  real CircleCenterX = 1.0;
  real CircleCenterY = 1.5;
  real CircleRadius = 0.5;

  // For a second circle, center coordinates and radius.
  real CircleCenterX0 = 1.65;
  real CircleCenterY0 = 1.6;
  real CircleRadius0 = 0.3;

  // For a third circle, center coordinates and radius.
  real CircleCenterX1 = 0.50;
  real CircleCenterY1 = 1.0;
  real CircleRadius1 = 0.25;

  // File containing a front defined by a set of points.
  string InitialFrontFile = "[full path]/[set].pts";
  // Is the (Xmin, Ymin) outside the front?
  bool origin_out = false;

  // If set to 'true', outside and inside are swapped.
  bool reversed = false;

  // Choose the right constructor:
  //   1) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                   reversed)
  //   2, 4) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                      CircleCenterX0, CircleCenterY0, CircleRadius0,
  //                      reversed)
  //   3, 5) InitialCurve(CircleCenterX, CircleCenterY, CircleRadius,
  //                      CircleCenterX0, CircleCenterY0, CircleRadius0,
  //                      CircleCenterX1, CircleCenterY1, CircleRadius1,
  //                      reversed)
  //   6) InitialCurve(InitialFrontFile, 0, origin_out ^ reversed)
  InitialCurveType InitialCurve(CircleCenterX, CircleCenterY, CircleRadius);


  ////////////////////////
  // LEVEL SET FUNCTION //
  ////////////////////////

  // Choose the level set function type:
  //   1) COrthogonalLevelSet<real>
  //   2) --
  typedef COrthogonalLevelSet<real> LevelSetType;

  // Choose the right constructor:
  //   1) Phi
  //   2) --
  LevelSetType Phi;


  /////////////////
  // INITIALIZER //
  /////////////////

  // Choose an initializer:
  //   1) CNarrowBandNeverInit<real> *
  //   2) CNarrowBandExtension<real> *
  //   3) CFastMarchingNeverInit<real> **
  //         * Narrow band level set method
  //        ** Fast marching method
  typedef CNarrowBandNeverInit<real> InitializerType;

  InitializerType Initializer;


  /////////////
  // UPDATER //
  /////////////

  // Choose an initializer:
  //   1) CNarrowBandFirstOrderEngquistOsher<real> *
  //   2) CNarrowBandFirstOrderLaxFriedrichs<real> *
  //   3) CNarrowBandEno2EngquistOsher<real> *
  //   4) CFastMarchingFirstOrderEngquistOsher<real> **
  //         * Narrow band level set method
  //        ** Fast marching method
  typedef CNarrowBandFirstOrderEngquistOsher<real> UpdaterType;

  // Choose the right constructor and choose its parameters:
  //   1, 2, 3) Updater(TubeSemiWidth, BarrierWidth, OutSpaceWidth)
  //              | TubeSemiWidth: number of grid points on each side of the front.
  //              | BarrierWidth: number of grid points (on each side of the front)
  //                        that imply tube reconstruction when reached.
  //              | OutSpaceWidth: number of grid points (on each side of the front)
  //                         that must not be reached.
  //              | Examples: Updater(6, 3, 1) or Updater(12, 5, 1).
  //   4) Updater(TMax)
  //        | TMax: time greater than all arrival times that will be computed.
  UpdaterType Updater(6, 3, 1);


  ///////////
  // SAVER //
  ///////////

  // Choose the saver type:
  //   1) CNeverSave<real>
  //        | Nothing is saved.
  //   2) CCurvesSaver<real>
  //        | The front is constructed and saved at each time step.
  //        | Not relevant for the fast marching method.
  //   3) CSaveLastCurve<real>
  //        | Saves the last curve, after all calculations.
  //        | Not relevant for the fast marching method.
  //   4) CSaveAtTheEnd<real>
  //        | Saves the level set, after all calculations.
  //        | Relevant for the fast marching method.
  typedef CCurvesSaver<real> SaverType;

  // Number of curves that will be saved.
  // Set NbCurves to 0 in order to save all curves.
  int NbCurves = 10;

  // Output directory.
  string Directory = "/home/vivien/Computations/Fronts/";

  // Stores time at each time step.
  string TimeFile = Directory + "Time";

  // Stores front points.
  string CurvesFile = Directory + "Curves";
  // Stores number of points on fronts at each time step.
  string CurveLengthsFile = Directory + "CurveLengths";
  // Stores the level set function.
  string PhiFile = Directory + "Phi";
  // Stores the speed function.
  string FFile = Directory + "F";

  // Stores grid coordinates along the (x'x) axe.
  string XFile = Directory + "X";
  // Stores grid coordinates along the (y'y) axe.
  string YFile = Directory + "Y";
  // Stores mesh points.
  string PointsFile = Directory + "Points";
  // Stores mesh edges.
  string EdgesFile = Directory + "Edges";
  // Stores mesh triangles.
  string TrianglesFile = Directory + "Triangles";

  // Indicates the number of iterations between saves.
  int Period = NbIterations / (NbCurves==0?NbIterations:NbCurves);

  SaverType Saver(TimeFile, CurvesFile,
		  CurveLengthsFile, PhiFile, FFile, XFile, YFile,
		  PointsFile, EdgesFile, TrianglesFile, Period);



  /////////////////////
  /////////////////////
  ////  SIMULATOR  ////
  /////////////////////
  /////////////////////

  CSimulator<real, MeshType, SpeedType, InitialCurveType,
    LevelSetType, InitializerType, UpdaterType, SaverType>
    Simulator(Mesh, F, InitialCurve, Phi,
	      Initializer, Updater, Saver,
	      NbIterations, FinalTime);


  ////////////////////////////////
  // SIMULATION INITIALIZATIONS //
  ////////////////////////////////

  Simulator.Init();


  ////////////////
  // SIMULATION //
  ////////////////

  Simulator.Run();

  // To catch exceptions.
  END;

  return 0;

}