#include "include\vs.h"
#define PI 3.141597
int main() {
   double f_0 = 1.0,
	// A. 2-D 5x5-point Bode's Integration Formula
          bode[5] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0},
	       weight[5][5];
   for(int i = 0; i < 5; i++)
   for(int j = 0; j < 5; j++) weight[i][j] = bode[i] * bode[j];
   Quadrature qp(weight[0], 0.0, 1.0, 5, 0.0, 1.0, 5);
   J d_a(pow( (1.0/4.0), 2.0));      // per normalized area

   // B. Define Basis Functions
   H1 x(2, (double*)0, qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
                                      9/*vector size*/, 2/*spatial dim.*/, qp);
   for(int m = 0; m < 3; m++)
   for(int n = 0; n < 3; n++)
   	phi[m*3+n] = sin((m+1.0)*PI*x[0])*sin((n+1.0)*PI*x[1]);

   // C. Variational Formulation
   H0 b_diag(9, (double*)0, qp); // LHS bilinear form is a diagonal matrix
   for(int i = 0; i < 9; i++)
   	b_diag[i] = d(phi[i]).pow(2);
   C0 B = b_diag | d_a;
   C0 l = ( ((H0)phi) * f_0 ) | d_a; // RHS forced vector

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