/*
   An example of a program solving
     -div(k grad u) + cu = f
   in a 2D domain 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"

//Include the header file defining the different
//  data structures that we use.
#include "2d_FE.h"

// Include a file that deals with getting the data 
//   for the mesh from a few files.
#include "readmesh.cc"

// Include a file with functions describing the linear
//   finite elements in 2D. 
#include "2d_linear.cc"

// Include a file with various quadratures. Currently only
//   a 3 point rule with nodes at the midpoints is implemented.
// This rule is exact for polynomials of degree 2 and thus
//   is good enough when working with linear elements.  
#include "quadrature_2d.cc"

/************************* The Coefficients ************************/
#define PI M_PI

double k(Point p)
{
   return 1;
}

double c(Point p)
{
   return 5;
}

double f(Point p)
{
   return  (5 + 18*PI*PI)*cos(3*PI*p.x)*cos(3*PI*p.y);
}


// an expected solution (used to check our approximation)
double sol(Point p)
{
   return cos(3*PI*p.x)*cos(3*PI*p.y);
}

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

int main()
{

   std::cout << "Reading the mesh..." << std::endl;
   if(!read_mesh("data.1"))
   {
      std::cout << "Error reading the mesh..." << std::endl;
      exit(1);
   }
   
   int global_dofs = num_nodes;

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

   Create_Quadratures();
   Quadrature *quad = &MIDPOINTS_RULE;

   //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;
   }

   std::cout << "Assembling the system..." << std::endl;
   
   // 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<num_elements; 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
         Point 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*fabs(jacobian(q_point,ele))*
                  (   k(q_point) * shape_funct_x(q_point, ele, i) * shape_funct_x(q_point, ele, j)
                      +
                      k(q_point) * shape_funct_y(q_point, ele, i) * shape_funct_y(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*fabs(jacobian(q_point, ele))*
                                     f(q_point) * shape_funct(q_point, ele, i);
         }
      }
   }
   delete element_dofs;

   std::cout << "Imposing the boundary conditions..." << std::endl;

   // impose the boundary conditions
   for (int i=0; i<global_dofs; i++)
   {
      // if this node is not at the boundary we don't care about it
      if (!nodes[i].bdmarker)
         continue;
      
      // if it is - we need to impose the Dirichlet BC
      b[i] = sol(nodes[i].p);
      A[i][i]=1;
      for (int j=0; j<global_dofs; j++)
      {
         if (j!=i)
         {
            A[i][j]=0; // <-- make the i'th row the unit vector (0,0,...,0,1,0,...0)
            b[j] -= A[j][i]*b[i]; A[j][i] = 0;  // <-- and remove the known value from the system
         }
      }
   }

   std::cout << "Solving the system..." << std::endl;
   // 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 error = fabs(b[i] - sol(nodes[i].p));
      std::cout << "| b[" << std::setw(3) << i << "] - u("<< std::setw(10) << nodes[i].p.x
                << ", " << std::setw(10) << nodes[i].p.y << ")|" << std::endl 
                << "                         = |"
                << std::setw(12) << b[i] << " - " << std::setw(12) << sol(nodes[i].p)
                << " | = " << 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();
}