#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() {
	double K = 1.0, weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, L_, 3);
   C0 delta_c, c_cache(3, (double*)0);
   do {
   	H1 x(qp), psi = x,
      	w = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
          	2, 1, qp);
   	w[0] = x; w[1] = x.pow(2);  // w
   	H0 b = INTEGRABLE_MATRIX("int, int, Quadrature", 3, 3, qp);
   	b[0][0] = E_*I_*d(psi)*d(psi)+K*((H0)psi)*((H0)psi);
   	b[0][1] = b[1][0] =  K*((H0)psi)*d(w[0]);
   	b[0][2] = b[2][0] =  K*((H0)psi)*d(w[1]);
   	b[1][1] = K*d(w[0])*d(w[0]); b[1][2] = K*d(w[0])*d(w[1]);
   	b[2][1] = K*d(w[1])*d(w[0]); b[2][2] = K*d(w[1])*d(w[1]);
   	C0 B = b |J(L_/2.0);                // LSH bilinear form
   	C0 l(3, (double*)0);          // RHS vector
   	l[0] = -M_*L_; l[1] = (((H0)w[0])*f_) | J(L_/2.0); l[2] = (((H0)w[1])*f_) | J(L_/2.0);
   	C0 c = l / B;
      delta_c &= c_cache - c;
      c_cache = c;
   	cout << c << ", " << norm(delta_c) << endl;
      //cout << l << endl;
      //cout << B << endl;
      //cout << norm(c) << endl;
      K *= 5.0;
   } while((double)norm(delta_c) > 1.e-6);
	return 0;
}