/*==============================================================================

                                    O  F  E  L  I

                            Object  Finite  Element  Library

  ==============================================================================

   Copyright (C) 1998 - 2004 Rachid Touzani

   This program is free software; you can redistribute it and/or modify it under
   the terms of the GNU General Public License as published by the Free 
   Software Foundation; Version 2 of the License.

   This program is distributed in the hope that it will be useful, but WITHOUT
   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 
   details.

   You should have received a copy of the GNU General Public License 
   along with this program; if not, write to the :

   Free Software Foundation
   Inc., 59 Temple Place - Suite 330
   Boston, MA  02111-1307, USA

  ==============================================================================

              An example of a Finite Element Code using OFELI

            Solution of a 1-D Elliptic problem using P1 Finite elements  

  ==============================================================================*/

#include "OFELI.h"
using namespace OFELI;

int main(int argc, char *argv[])
{
   double   Lmin=0, Lmax=1;
   int      N = 10;
   double   f(double x);
   double   error(const Mesh &ms, const Vect<double> &u);

// Read and output mesh data
   banner();
   if (argc > 1)
     N = atoi(argv[1]);
   Mesh ms(Lmin,Lmax,N);
   int NbN = N+1;
   ms.Verbose(10);
   cout << ms;

// Declare problem data (matrix, rhs, boundary conditions, body forces)
   TrMatrix<double> a(NbN);
   Vect<double> b(NbN);

// Build matrix and R.H.S.
   double h = (Lmax-Lmin)/double(N);
   for (int i=2; i<NbN; i++) {
      double x = ms.PtrNode(i)->Coord(1);
      a.Set(i,i,2./h);
      a.Set(i,i+1,-1./h);
      a.Set(i,i-1,-1./h);
      b(i) = f(x)*h;
   }

// Impose boundary conditions
   a.Set(1,1,1.); a.Set(1,2,0.);
   b(1) = 0;
   a.Set(NbN,NbN,1.); a.Set(NbN-1,NbN,0.);
   b(NbN) = 0;

// Solve the linear system of equations
   a.Solve(b);

// Output solution and error
   cout << "\nSolution :\n" << b;
   cout << "Error = " << error(ms,b) << endl;

#ifdef WITH_PAUSE
   system("PAUSE");
#endif
   return 0;
}


double f(double x)
{
   double a = 10;
   return 2*a*(1-2*a*x*x)*exp(-a*x*x);
}


double error(const Mesh &ms, const Vect<double> &u)
{
   double err=0;
   double a = 10;
   for (int i=1; i<=u.Size(); i++) {
      double x = ms.PtrNode(i)->Coord(1);
      err = max(err,fabs(u(i) - exp(-a*x*x) - x*(1-exp(-a)) + 1));
   }
   return err;
}