#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 weight[5] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   Quadrature qp(weight, 0.0, L_, 5.0);
   J d_l(L_/4.0);
	/*double weight[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(weight, 0.0, L_, 17);
   J d_l(L_/16.0);*/
	//double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   //Quadrature qp(weight, 0.0, L_, 3);
   //J d_l(L_/2.0);
   //double weight[4] = {3.0/8.0, 9.0/8.0, 9.0/8.0, 3.0/8.0};
   //Quadrature qp(weight, 0.0, L_, 4);
   //J d_l(L_/3.0);
   // B. Define Basis Functions
   H1 x(qp),
      w = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp),
      M = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 2/*vector size*/, 1/*spatial dim.*/, qp);
   w[0] = x;     w[1] = x.pow(2);        // w = c_0 x + c_1 x^2
   M[0] = x-L_;  M[1] = (x-L_).pow(2);   // M = M_ + c_2 (x-L) + c_3 (x-L)^2

   // C. Weak Formulation
   C0 e(2), E(2);  // Cartesian basis
   C0 B =+( 0.0            *((e%e)*(E[0]%E[0])) + (d(w)*(~d(M)))                 *((e%e)*(E[0]%E[1]))+
            (d(M)*(~d(w))) *((e%e)*(E[1]%E[0])) + ( ((H0)M)*(~((H0)M))/(E_*I_) ) *((e%e)*(E[1]%E[1]))
           ) | d_l; // LSH bilinear form; use operator+() to convert to integrable_matrix
   C0 l = +( ( (-(H0)w)*f_ )           *(e*E[0])+// forced vector     : l(w, -f) = -(((H0)w), f)
             ( (-(H0)M)*(M_/(E_*I_)) ) *(e*E[1]) // essential bouncary condition of M: -B(M, M_) = -((H0)M)*M_/(E_*I_)
           ) | d_l;

   // D. Substructuring--static condensation(Note: B00 = null)
   C0 B_01 = MATRIX("int, int, C0&, int, int", 2, 2, B, 0, 2),
      B_10 = MATRIX("int, int, C0&, int, int", 2, 2, B, 2, 0),
      B_11 = MATRIX("int, int, C0&, int, int", 2, 2, B, 2, 2);
	C0 l_0 = VECTOR("int, C0&, int", 2, l, 0),
      l_1 = VECTOR("int, C0&, int", 2, l, 2);
	C0 c(4, (double*)0),
      c_0 = VECTOR("int, C0&, int", 2, c, 0),
      c_1 = VECTOR("int, C0&, int", 2, c, 2);
   cout << B_01 << endl << B_10 << endl << B_11 << endl;
   //c_1 = (l_0/B_01);                 // solve for coefficient for M
   //c_0 = (l_1 - B_11 * c_1) / B_10;  // solve for coefficient for w*/
   c = l / B;
   /*C0 M_inv = B_11.inverse(),
      BtMinvB = B_01 * M_inv * B_10;
   c_0 = BtMinvB.inverse() * (B_01*M_inv*l_1 - l_0);
   c_1 = M_inv * (l_1 - B_10 * c_0);*/
   cout << c << endl;
   return 0;
}