/*********************************************************************/
/* File:   recurisve_pol.cpp                                         */
/* Author: Almedin Becirovic, JS                                     */
/* Date:   14. Apr. 2003                                             */
/*********************************************************************/


#include <fem.hpp>
namespace ngfem
{
  using namespace ngfem;



  
  template <int DIM>
  double VertexExtensionOptimal<DIM> :: coefs[SIZE][SIZE];

  template <int DIM>
  bool VertexExtensionOptimal<DIM> :: initialized = 0;


  template <int DIM>
  VertexExtensionOptimal<DIM> :: VertexExtensionOptimal ()
  {
    if (!initialized)
      {
	// ofstream tout ("test1.out");
	int i, j, k;
	Matrix<double> a(SIZE, SIZE);
	Matrix<double> b(2, SIZE);

	try
	  {
	    const IntegrationRule & ir =
	      GetIntegrationRules().SelectIntegrationRule (ET_SEGM, 2*SIZE+2);

	    Vector<double> diff_leg(SIZE); 

	    a = 0;
	    b = 0;

	    for (i = 0; i < ir.GetNIP(); i++)
	      {
		double x = ir.GetIP(i)(0);
		DerivedLegendrePolynomial (SIZE-2, 2*x-1, diff_leg);
		for (j = SIZE-1; j >= 1; j--)
		  diff_leg(j) = 2 * diff_leg(j-1);
		diff_leg(0) = 0;

		double fac = ir.GetIP(i).Weight();
		for (j = 0; j < DIM-1; j++)
		  fac *= (1-x);
	    
		for (j = 0; j < SIZE; j++)
		  for (k = 0; k < SIZE; k++)
		    a(j,k) += fac * diff_leg(j) * diff_leg(k);
	      }

	    double fac = 1;
	    for (j = 0; j < SIZE; j++)
	      {
		b(0,j) = fac;
		b(1,j) = 1;
		fac *= -1;
	      }


	    // compute coefficients
	    for (i = 1; i < SIZE; i++)
	      {
		Matrix<double> mat(i+3, i+3);
		Vector<double> rhs(i+3), sol(i+3);
		mat = 0;
		for (j = 0; j <= i; j++)
		  for (k = 0; k <= i; k++)
		    mat(j,k) = a(j,k);

		for (j = 0; j <= i; j++)
		  {
		    mat(i+1, j) = mat(j,i+1) = b(0,j);
		    mat(i+2, j) = mat(j,i+2) = b(1,j);
		  }
	    
		rhs = 0;
		rhs(i+2) = 1;

		// stabilize A
		for (j = 0; j <= i; j++)
		  {
		    for (k = 0; k <= i; k++)
		      mat(j,k) += b(0,j) * b(0,j) + b(1,j) * b(1,k);
		    rhs(j) += b(1,j);
		  }
	    
		CholeskyFactors<double> inv(mat);
		inv.Mult (rhs, sol);

		// tout << "sol = " << sol << endl;

		for (j = 0; j <= i; j++)
		  coefs[j][i] = sol(j);
	      }
	

	    /*
	    // 	    // testing
	    ofstream out ("vext.dat");
	    ofstream dout ("dvext.dat");
	    for (double x = 0; x <= 1; x += 0.01)
	      {
		out << x << " ";
		for (i = 1; i <= 10; i++)
		  out << Calc (i, x) << " ";
		out << endl;
		
		dout << x << " ";
		for (i = 1; i <= 10; i++)
		  dout << CalcDeriv (i, x) << " ";
		dout << endl;
	      }
	    */
	    initialized = 1;
	  }
	catch (Exception e)
	  {
	    cout << "Caught in VertexExtensionOptimal:" << e.What() << endl;
	  }
      }
}

  /*

  template <int DIM>
  VertexExtensionOptimal<DIM> :: VertexExtensionOptimal ()
  {
    if (!initialized)
      {

	ofstream tout;
	if (DIM == 3)
	  tout.open("test3d.out");
	else
	  tout.open("test2d.out");

	int i, j, k;
	Matrix<double> a(SIZE, SIZE);
	Matrix<double> b(2, SIZE);


	ofstream jacout ("jacobi.out");
	Vector<> jac(11);
	for (double x = -1; x <= 1; x += 0.01)
	  {
	    JacobiPolynomial (10, x, 1, -1, jac);
	    for (int j = 0; j <= 10; j++)
	      jacout << " " << jac(j) << " ";
	    jacout << endl;
	  }

	try
	  {
	    const IntegrationRule & ir =
	      GetIntegrationRules().SelectIntegrationRule (ET_SEGM, 2*SIZE+2);

	    Vector<double> diff_leg(SIZE); 
	    Vector<double> help(SIZE); 

	    a = 0;
	    b = 0;

	    for (i = 0; i < ir.GetNIP(); i++)
	      {
		double x = ir.GetIP(i)(0);

		JacobiPolynomial (SIZE-2, 2*x-1, 2, 0, help);
		diff_leg(0) = 0;
		for (j = 1; j < SIZE; j++)
		  diff_leg(j) = 0.5 * (j+1) * 2.0 * help(j-1);

		double fac = ir.GetIP(i).Weight();
		for (j = 0; j < DIM-1; j++)
		  fac *= (1-x);
	    
		for (j = 0; j < SIZE; j++)
		  for (k = 0; k < SIZE; k++)
		    a(j,k) += fac * diff_leg(j) * diff_leg(k);
	      }

	    JacobiPolynomial (SIZE-1, -1.0, 1, -1, help);
	    for (j = 0; j < SIZE; j++)
	      b(0,j) = help(j);
	    JacobiPolynomial (SIZE-1, 1.0, 1, -1, help);
	    for (j = 0; j < SIZE; j++)
	      b(1,j) = help(j);


	    tout << "mata = " << a << endl;
	    tout << "matb = " << b << endl;

	    tout << "diag (A), 1/diag(A),   bi*bi/aii " << endl;
	    for (j = 1; j < SIZE; j++)
	      tout << a(j,j) << ",    " << 1.0/a(j,j) << "    " << b(1,j)*b(1,j)/a(j,j) << endl;

	    // compute coefficients
	    for (i = 1; i < SIZE; i++)
	      {
		Matrix<double> mat(i+3, i+3);
		Vector<double> rhs(i+3), sol(i+3);
		mat = 0;
		for (j = 0; j <= i; j++)
		  for (k = 0; k <= i; k++)
		    mat(j,k) = a(j,k);

		for (j = 0; j <= i; j++)
		  {
		    mat(i+1, j) = mat(j,i+1) = b(0,j);
		    mat(i+2, j) = mat(j,i+2) = b(1,j);
		  }
	    
		// stabilize A
		rhs = 0;
		rhs(i+2) = 1;

		tout << "mat,1 = " << mat << endl;
		tout << "rhs,1 = " << rhs << endl;

		for (j = 0; j <= i; j++)
		  {
		    for (k = 0; k <= i; k++)
		      mat(j,k) += b(0,j) * b(0,k) + b(1,j) * b(1,k);
		    rhs(j) += b(1,j);
		  }
	    
		tout << "mat,2 = " << mat << endl;
		tout << "rhs,2 = " << rhs << endl;

		CholeskyFactors<double> inv(mat);
		inv.Mult (rhs, sol);

		tout << "sol = " << sol << endl;

		for (j = 0; j <= i; j++)
		  coefs[j][i] = sol(j);

		double schur = 0;
		for (j = 1; j <= i; j++)
		  schur += b(1,j) * b(1,j) / a(j,j);
		double lam = 1.0 / schur;

		tout << "schur = " << schur << ", lam = " << lam << endl;

		sol(0) = 0;
		for (j = 1; j <= i; j++)
		  sol(j) = (2*j+1) / ( (j+1.0)*(i*i+2*i));
		tout << "sol,2 = " << sol << endl;
	      }
	

	    // testing
	    ofstream out ("vext.dat");
	    ofstream dout ("dvext.dat");
	    for (double x = 0; x <= 1.0001; x += 0.01)
	      {
		out << x << " ";
		for (i = 1; i <= 10; i++)
		  out << Calc (i, x) << " ";
		out << endl;
		
		dout << x << " ";
		for (i = 1; i <= 10; i++)
		  dout << CalcDeriv (i, x) << " ";
		dout << endl;
	      }
	    
	    initialized = 1;
	  }
	catch (Exception e)
	  {
	    cout << "Caught in VertexExtensionOptimal:" << e.What() << endl;
	  }
      }
  }
  */



  /*
  template <int DIM>
  VertexExtensionOptimal<DIM> :: VertexExtensionOptimal ()
  {
    if (!initialized)
      {

	ofstream tout;
	if (DIM == 3)
	  tout.open("test3d.out");
	else
	  tout.open("test2d.out");

	int i, j, k;
	Matrix<double> a(SIZE, SIZE);
	Matrix<double> b(2, SIZE);


	ofstream jacout ("jacobi.out");
	Vector<> jac(11);
	for (double x = -1; x <= 1; x += 0.01)
	  {
	    JacobiPolynomial (10, x, 1, -1, jac);
	    for (int j = 0; j <= 10; j++)
	      jacout << " " << jac(j) << " ";
	    jacout << endl;
	  }

	try
	  {
	    const IntegrationRule & ir =
	      GetIntegrationRules().SelectIntegrationRule (ET_SEGM, 2*SIZE+4);

	    Vector<double> diff_leg(SIZE); 
	    Vector<double> help(SIZE); 
	    Vector<double> help2(SIZE); 

	    a = 0;
	    b = 0;

	    for (i = 0; i < ir.GetNIP(); i++)
	      {
		double x = ir.GetIP(i)(0);
		double hx = 2*x-1;

		JacobiPolynomial (SIZE-1, hx, 1, 1, help);
		JacobiPolynomial (SIZE-1, hx, 2, 2, help2);

		diff_leg(0) = 1;
		for (j = 1; j < SIZE; j++)
		  diff_leg(j) = help(j) + (1+hx)*0.5*(j+3)*help2(j-1);

		
		double fac = ir.GetIP(i).Weight();
		for (j = 0; j < DIM-1; j++)
		  fac *= (1-x);
	    
		for (j = 0; j < SIZE; j++)
		  for (k = 0; k < SIZE; k++)
		    a(j,k) += fac * diff_leg(j) * diff_leg(k);
	      }

	    JacobiPolynomial (SIZE-1, -1.0, 1, -1, help);
	    for (j = 0; j < SIZE; j++)
	      b(0,j) = help(j);
	    JacobiPolynomial (SIZE-1, 1.0, 1, -1, help);
	    for (j = 0; j < SIZE; j++)
	      b(1,j) = help(j);


	    tout << "mata = " << a << endl;
	    tout << "matb = " << b << endl;

	    tout << "diag (A), 1/diag(A),   bi*bi/aii " << endl;
	    for (j = 1; j < SIZE; j++)
	      tout << a(j,j) << ",    " << 1.0/a(j,j) << "    " << b(1,j)*b(1,j)/a(j,j) << endl;

	    // compute coefficients
	    for (i = 1; i < SIZE; i++)
	      {
		Matrix<double> mat(i+3, i+3);
		Vector<double> rhs(i+3), sol(i+3);
		mat = 0;
		for (j = 0; j <= i; j++)
		  for (k = 0; k <= i; k++)
		    mat(j,k) = a(j,k);

		for (j = 0; j <= i; j++)
		  {
		    mat(i+1, j) = mat(j,i+1) = b(0,j);
		    mat(i+2, j) = mat(j,i+2) = b(1,j);
		  }
	    
		// stabilize A
		rhs = 0;
		rhs(i+2) = 1;

		tout << "mat,1 = " << mat << endl;
		tout << "rhs,1 = " << rhs << endl;

		for (j = 0; j <= i; j++)
		  {
		    for (k = 0; k <= i; k++)
		      mat(j,k) += b(0,j) * b(0,k) + b(1,j) * b(1,k);
		    rhs(j) += b(1,j);
		  }
	    
		tout << "mat,2 = " << mat << endl;
		tout << "rhs,2 = " << rhs << endl;

		CholeskyFactors<double> inv(mat);
		inv.Mult (rhs, sol);

		tout << "sol = " << sol << endl;

		for (j = 0; j <= i; j++)
		  coefs[j][i] = sol(j);

		double schur = 0;
		for (j = 1; j <= i; j++)
		  schur += b(1,j) * b(1,j) / a(j,j);
		double lam = 1.0 / schur;

		tout << "schur = " << schur << ", lam = " << lam << endl;

		sol(0) = 0;
		for (j = 1; j <= i; j++)
		  sol(j) = (2*j+1) / ( (j+1.0)*(i*i+2*i));
		tout << "sol,2 = " << sol << endl;
	      }
	

	    // testing
	    ofstream out ("vext.dat");
	    ofstream dout ("dvext.dat");
	    for (double x = 0; x <= 1.0001; x += 0.01)
	      {
		out << x << " ";
		for (i = 1; i <= 10; i++)
		  out << Calc (i, x) << " ";
		out << endl;
		
		dout << x << " ";
		for (i = 1; i <= 10; i++)
		  dout << CalcDeriv (i, x) << " ";
		dout << endl;
	      }
	    
	    initialized = 1;
	  }
	catch (Exception e)
	  {
	    cout << "Caught in VertexExtensionOptimal:" << e.What() << endl;
	  }
      }
  }
  */




  VertexExtensionOptimal<2> dummy_vextopt2;
  VertexExtensionOptimal<3> dummy_vextopt3;






  double IntegratedLegendreMonomialExt :: coefs[SIZE][2];



  namespace
  {
    class Init
    {
      public:
      Init ()
      {
	IntegratedLegendreMonomialExt :: CalcCoeffs();
      }
    };

    Init init;
  }


}