/*
   An example of a program solving
     -(ku')' + cu = f
   in the 1D interval [0,l] using linear
   finite elements.
*/  
   
#include <iostream>
#include <iomanip>
#include "math.h"

//Include our own previously written code...
//  in general we will be using iterative solvers,
//  but for this particular example (since we're
//  assembling the full global matrix) we'll use
//  the Gaussian Elimination direct solver that
//  we have written for Math609.
#include "gaussian.cc"


int n = 25;
double l = 10;
double h = l/n;

// Include a file with functions describing the linear
//   finite elements. If we wish to change our program to
//   work with quadratic elements, here's where we'll make
//   the changes.
#include "1d_linear.cc"

// Include a file with various quadratures (currently only
//   the Simpson's rule for the unit interval is implemented).
#include "quadrature.cc"

/************************* The Coefficients ************************/

double k(double x)
{
   return 2;
}

double c(double x)
{
   return 3;
}

double f(double x)
{
   return -12*x+3*x*x*x;
}


// an expected solution (used to check our approximation)
double sol(double x)
{
   return x*x*x;
}

/************************* The Main Code ************************/

int main()
{

   double** A = new double*[n+1];
   double*  b = new double[n+1];

   Create_Quadratures();
   Quadrature *quad = &SIMPSONS_RULE;

   //create the empty matrix and rhs
   for (int i=0; i<=n; i++)
   {
      A[i] = new double[n+1];
      for (int k=0; k<=n; k++)
         A[i][k]=0;
      b[i] = 0;
   }
   
   // cycle over each of the elements
   //   and fill in the local contributions
   int num_dofs = dofs_per_element();
   int* element_dofs = new int[num_dofs];
   for (int ele=0; ele<n; ele++)
   {
      fill_element_dofs(element_dofs, ele);

      // for all the points in our quadrature
      for (int l = 0; l < quad->size; l++)
      {
         // get the values of the point and the weight
         double q_point = map_to_FE(quad->points[l], ele);
         double q_weight = quad->weights[l];
         for (int i=0; i<num_dofs; i++)
         {
            for (int j=0; j<num_dofs; j++)
            {
               // and calculate the value at a particualar point that is included in the
               //   quadrature used in computing the integral (k*phi_i',phi_j')+(c*phi_i,phi_j)
               A[element_dofs[i]][element_dofs[j]] += q_weight*jacobian(q_point,ele)*
                  (   k(q_point) * shape_funct_x(q_point, ele, i) * shape_funct_x(q_point, ele, j)
                      +
                      c(q_point) * shape_funct(q_point, ele, i) * shape_funct(q_point, ele, j)
                  );
            }

            // ... and in a simillar manner - in computing the integral (f,phi_i)
            //  for the right hand side
            b[element_dofs[i]] += q_weight*jacobian(q_point, ele)*
                                     f(q_point) * shape_funct(q_point, ele, i);
         }
      }
   }
   delete element_dofs;

/*
   // Print out the system we're solving for debug purposes
   for (int i=0; i<=n; i++)
   {
      for (int j=0; j<=n; j++)
         std::cout << A[i][j] << " ";
      std::cout << "b=" << b[i] << std::endl;
   }
*/

   // impose the boundary conditions
   {
      // force the boundary condition at 0
      b[0] = sol(0);
      A[0][0]=1;
      for (int i=1; i<=n; i++)
      {
         A[0][i]=0; // <-- make the first row the unit vector (1,0,0...0)
         b[i] -= A[i][0]*b[0]; A[i][0] = 0;  // <-- and remove the known value from the system
      }
      // force the boundary condition at L
      b[n] = sol(l);
      A[n][n]=1;
      for (int i=0; i<n; i++)
      {
         A[n][i]=0; b[n]=sol(l); //<-- make the last row the unit vector (0,0,...0,1)
         b[i] -= A[i][n]*b[n]; A[i][n]=0; //<-- remove the known value from the system
      }
   }

   // solve the system.
   gaussian_elimination(A,b,n+1);

   // Print out the numerical solution, the exact solution and the error
   //  at the nodes.
   double max = 0;
   for (int i=0; i<=n; i++)
   {
      double error = fabs(b[i] - sol(i*h));
      std::cout << "| b[" << std::setw(6) << i << "] - u("<< std::setw(6) << i*h <<")| = |"
                << std::setw(12) << b[i] << " - " << std::setw(12) << sol(i*h)
                << " | = " << std::setw(12) << error << std::endl;
      if (max < error)
         max = error;
   }
   std::cout << std::endl << "Discrete l-infiniy norm : " << max << std::endl;

   for (int i=0; i<=n; i++)
      delete A[i];
   delete A;
   delete b;

   Destroy_Quadratures();
}