#include "vs.h"
const double EPSILON = 1.e-12;
int main() { // three parameters approximation; Bubnov-Galerkin Method only
	// A. Bode's Integration Formula
   double weight[13] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   Quadrature qp(weight, 0.0, 1.0, 13);
   J d_l(1.0/12.0);      // per normalized length

   // B. Define Basis Functions
   H1 x(qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
            3/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = 1.0-x; phi[1] = 1.0-x.pow(2); phi[2] = 1.0-x.pow(3);
#if defined(__PETROV_GALERKIN)
#define PI 3.141592654
   H1 psi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
            3/*vector size*/, 1/*spatial dim.*/, qp);
   psi[0] = cos(PI/2.0*x); psi[1] = cos(3.0*PI/2.0*x); psi[2] = cos(5.0*PI/2.0*x);
#else
	H1 psi;
   psi &= phi;
#endif

   // C. Weighted-Residual Formulation
   C0 c(3,(double*)0), delta_c(3,(double*)0);
   do { // Newton's Method
      H1 v = c * phi;
   	C0 I = ( d(psi)(0) * ( ((H0)v)+sqrt(2.0) ) * d(v)[0] + ((H0)psi) ) | d_l,
   	   I_t = ( d(psi)(0) % (((H0)phi) * d(v)[0] +
                              d(phi)(0) * ( ((H0)v)+sqrt(2.0) )
                             )
              ) | d_l;
      delta_c = -I /I_t; // Newton's formula
      c += delta_c;                                // update increment
      cout << c << ", " << norm(delta_c) << endl;
   } while((double)norm(delta_c ) > EPSILON); // check for convergence
   cout << c << endl;
   double x_ = 0.5;
   cout << "exact solution @ 0.5: " << (sqrt(1.0+pow(x_,2))) << endl;
   cout << "approx. @ 0.5: " << ( c[0]*(1.0-x_)+c[1]*(1.0-pow(x_,2))+c[2]*(1.0-pow(x_,3))+sqrt(2.0) ) << endl;
	return 0;
}