/*
   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 = 20;
double l = 1;
double h = l/n;

// Include a file with functions describing the Lagrangian
//   finite elements. 
#include "1d_lagrangian.cc"

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

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

double S = 100;
double q = 200;
double D = 8.8*pow(10,7);

double k(double x)
{
   return 1;
}

double c(double x)
{
   return 4;
}

double f(double x)
{
   return -120*x*x+40*x*x*x*x;
}


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

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

int main()
{
   int global_dofs = n+1+n*(DEGREE-1);

   double** A = new double*[global_dofs];
   double*  b = new double[global_dofs];

   Create_Quadratures();
   Quadrature *quad = &FIVE_POINT_GAUSS;

   //create the empty matrix and rhs
   for (int i=0; i<global_dofs; i++)
   {
      A[i] = new double[global_dofs];
      for (int k=0; k<global_dofs; 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<global_dofs; i++)
   {
      for (int j=0; j<global_dofs; 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<global_dofs; 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[global_dofs-1] = sol(l);
      A[global_dofs-1][global_dofs-1]=1;
      for (int i=0; i<global_dofs-1; i++)
      {
         A[global_dofs-1][i]=0; //<-- make the last row the unit vector (0,0,...0,1)
         b[i] -= A[i][global_dofs-1]*b[global_dofs-1]; A[i][global_dofs-1]=0; //<-- remove the known value from the system
      }
   }

   // solve the system.
   gaussian_elimination(A,b,global_dofs);

   // Print out the numerical solution, the exact solution and the error
   //  at the nodes.
   double max = 0;
   for (int i=0; i<global_dofs; i++)
   {
      double step = 1.0/(global_dofs-1);
      double error = fabs(b[i] - sol(i*step));
      std::cout << "| b[" << std::setw(6) << i << "] - u("<< std::setw(6) << i*step <<")| = |"
                << std::setw(12) << b[i] << " - " << std::setw(12) << sol(i*step)
                << " | = " << 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<global_dofs; i++)
      delete A[i];
   delete A;
   delete b;

   Destroy_Quadratures();
}