#include "vs.h"
#define L_ 1.0
#define E_ 1.0
#define I_ 1.0
#define f_ 1.0
#define M_ 1.0

int main() {
   for(int N = 2; N < 4; N++) {
		// A. Bode's Integration Formula--exact up to 5-th order
		double w[17] = {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, 28.0/45.0,
                              	64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   	Quadrature qp(w, 0.0, L_, 17);
   	// B. Define Basis Functions
   	H2 x((double*)0, qp),
      	phi = INTEGRABLE_VECTOR_OF_TANGENT_OF_TANGENT_BUNDLE("int, int, Quadrature",
                                      N/*vector size*/, 1/*spatial dim.*/, qp);
   	for(int i = 0; i < N; i++) phi[i] = x.pow(i+2);
   
   	// C. Variational Formulation
   	H0 d2_phi = INTEGRABLE_VECTOR("int, Quadrature", N, qp);
   	for(i = 0; i < N; i++) d2_phi[i] = dd(phi)(i)[0][0]; // degenerated hessian matrix
   	C0 B = ( E_* I_* (d2_phi%d2_phi) ) | J(L_/16.0);       // LSH bilinear form
   	C0 M_bc(N, (double*)NULL);
   	M_bc[0] = 2.0*L_*M_; M_bc[1] = 3.0*pow(L_,2.0)*M_; if(N==3) M_bc[2] = 4.0*pow(L_,3.0)*M_;
   	C0 l = ( ( ((H0)phi) * f_ ) | J(L_/16.0) ) // forced vector: (((H0)phi), f) +
           	+ M_bc;                        //                (d(phi(L))*M)

   	// D. Solutions--the Ritz Coefficients
   	C0 c = l / B;
   	for(i = 0; i < N; i++) cout << c[i] << endl;
   }
   return 0;
}