Demo: The 2-D Laplace equation: P1 Finite Element

We present in this demo a more consistent use of the library. We aim at solving the basic Laplace equation (more precisely the Poisson problem) in 2-D using P1 finite elements.

  • As usual we start by including tha main head file. Next, since we aim at using the one of classes for the Laplace equation, we have to include the file Laplace.h.

       #include "OFELI.h"
   #include "Laplace.h"
   using namespace OFELI;
                           
   int main(int argc, char *argv[])
   {

  • The program has an argument the name of the the project file. If the argument is available, we instantiate the class IPF which manages project files.

          if (argc < 2) {
         cout << "\nUsage: " << argv[0] << " <project_file>\n";
         return 0;
      }
      IPF data("laplace_demo1 - 1.0",argv[1]);

  • As in every demo program we make use of exception handling to obtain clear error messages. Next we assign to the global variable Verbosity the value obtained from the project file. 

          Verbosity = data.getVerbose();

  • The first essential task in the program consists in constructing the mesh using the mesh file.

          Mesh ms(data.getMeshFile());

  • We have now to set some data for the equation. Namely, the Dirichlet boundary condition and the source term. They will be stored in dedicated files. To initialize them we make use of the class Prescription that helps prescribing various data. This class reads these data from a Prescription file, which can be the project file itself. The Dirichlet boundary conditions and the source are retrieved by member function get. 

          Vect<double> bc(ms), bf(ms);
      Prescription p(ms,data.getPrescriptionFile());
      p.get(BOUNDARY_CONDITION,bc);
      p.get(SOURCE,bf);

  • We declare a vector that will contain the solution once the problem solved. Then we declare an instance of class Laplace2DT3 that handles the numerical solution of the Laplace equation in 2-D using P1 finite elements (3-Node triangles). We note that the constructor used the Mesh instance and the solution vector that will contain later the solution.

          Vect u(ms);
      Laplace2DT3 eq(ms,u);

  • Now we can transmit to the equation class various prescribed data using membre function setInput: 

          eq.setInput(BOUNDARY_CONDITION,bc);
      eq.setInput(SOURCE,bf);

  • Once all data are transmitted, we can run (solve) the equation. The vector u contains now the solution and can be solved in a file. For this end, we choose here the Gmsh .pos file. 

          eq.run();
      saveField(u,"u.pos",GMSH);

  • To check the accuracy of the numerical solution, we have constructed an analytical solution of the equation. This one is also defined via the Prescription class and given as an algebraic expression in the project file. We can then compute the error norm as indicated in the source and terminate the program. 

          Vect<double> sol(ms);
      p.get(SOLUTION,sol);
      cout << "L2-Error: " << (u-sol).Norm(WNORM2) << endl;
      return 0;
   }