///////////////////////////////////////////////////////////////////
// Mladen Jurak
// Praktikum primijenjene matematike 2, ak. god. 2004/05
// FreeFEM++
//
// Konačni elementi i varijacijska formulacija 3. 
// Dirichletov rubni uvjet.
////////////////////////////////////////////////////////////////////

// Dirichletov rubni uvjet može se implementirati na različite načine. 
// Svaka implementacija zahtijeva određenu modifikaciju matrice krutosti
// i vektora desne strane. 
//
// U FreFem++-u se Dirichletov rubni uvjet implementira na sljedeći način: 
// 1. Matrica krutosti i vektor desne strane se formiraju bez uvažavanja 
//    rubnih uvjeta, dakle kao da je zadan homogeni Neumannov rubni uvjet.
// 2. Redak matrice koji odgovara nodalnoj točki q^i u kojoj je  zadan Dirichletov
//    rubni uvjet modificira se tako da se na dijagonalu stavi neki veliki broj E,
//    po defaultu E=1.0e+30.
// 3. U vektoru desne strane na mjestu koje odgovara točki q^i doda se E g(q^i),
//    gdje je g() funkcija kojom je dan Dirichletov rubni uvjet. 
//
// ILUSTRACIJA metode:

//  Riješimo Laplaceovu jednadžbu 
//                         -u_xx -u_yy =f 
//  u jediničnom krugu, s desnom stranom 
//                           f(x,y) = x*y
//  i Dirichletovim rubnim uvjetom 
//                           g(x,y) = sin(pi*x)*cos(pi*y)

    border C(t=0,2*pi) {
                  x=cos(t);
                  y=sin(t);
                  label=1;
    }

//  Sa P1 elementima dobit ćemo 7 trokuta s jednom unutarnjom točkom.
    mesh    Th = buildmesh(C(7));
    fespace Vh(Th,P1);
    Vh u, v;
    func f = x*y;
    func g = sin(pi*x)*cos(pi*y);
//    func g=2.0;

    problem Poisson(u,v,solver=LU) =
        int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinearna forma 
       -int2d(Th)(f*v)                     // linearni funkcional
       + on(1,u=g);                        // rubni uvjet

    Poisson;
//    plot(Th, u, wait=true);

//  Isti problem možemo riješiti definirajući posebno matricu krutosti i desnu 
//  stranu. (Parametri u varf su posve formalni, što nije slučaj s problem i solve.)

    varf a(u,v) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v)) // bilinearna forma
                  +  on(1,u=g);                      // Dirichletov r.u. ugrađen

    matrix A = a(Vh,Vh);                             // Matrica krutosti

// Matrica se pamti u rijetkoj retčanoj formi (Morseov format)
// Na dijagonali se pojavila konstanta E=1.0e+30

    cout << "Matrica krutosti:"<<endl;
    cout << A<<endl;                                 // Ispiši matricu krutosti

    varf b(v,f) = int2d(Th)(f*v);                    // Isto za desnu stranu.
    matrix B = b(Vh,Vh);

    Vh F, fh=f;

    // * je operator množenja (matrica x vektor)
    // U vektoru desne strane nema modifikacije zbog rubnih uvjeta.
    F[] = B*fh[];
    cout << "Vektor desne strane prije modifikacije rubnim uvjetom:"<<endl;
    cout <<"F="<<F[]<<endl;                          // Ispiši vektor desne strane

    // Da bismo mogli modificirati desnu stranu, FreeFemm++ nudi sljedeću 
    // sintaksu: funkcija a(0,Vh) ima vrijednost 0 u  svakoj točki koja nije na 
    // Dirichletovoj granici, a E*g(q^i) u ostalim točkama.
    // NAPOMENA: Postoji greška u dokumentaciji!!!! PROVJERA:

    Vh bc; bc[] = a(0,Vh);
    cout << "Vektor bc[]:"<<endl;                       // Ispišimo bc[]
    cout <<"bc="<<bc[]<<endl;

    Vh gh=g;
    cout << "Vektor gh[] (rubne vrijednosti):"<<endl;
    cout <<"gh="<<gh[]<<endl;                           // Ispišimo rubni uvjet

    Vh mm; mm[] = bc[]./gh[];
    cout << "Vektor bc[]./gh[]:"<<endl;
    cout <<"bc[]./gh[]="<<mm[]<<endl;                   // mm mora biti 0 ili 1.0e+30

    // Korekcija desne strane. Treba samo dodati rubni uvjet množen s E:
    F[] += bc[];
    cout << "Vektor desne strane nakon modifikacije rubnim uvjetom:"<<endl;
    cout <<"F="<<F[]<<endl;

    Vh w, diff;
    w[] = A^-1*F[];                                  // Riješi sustav

    diff[]=u[]-w[];                                  // Razlika dvaju rješenja

    cout<<"Razlika izmedju dva rjesenja diff = "<<diff[].max<<endl;


//////////////////////////////////////////////////////////////////////////////////////////////
//
//  Način pamćenja matrice ovisi o izabranom "solveru".

   matrix AA=a(Vh,Vh,solver=LU);   // A = LDU  sky-line nesimetrična
   cout << "Matrica u sky-line formatu, nesimetricna:"<<endl;
   cout << "AA="<<AA<<endl;


   matrix AAAAA=a(Vh,Vh,solver=Crout);  // sky-line simetrična forma
   cout << "Matrica u sky-line formatu, simetricna:"<<endl;
   cout << "AAAAA="<<AAAAA<<endl;

   matrix AAA=a(Vh,Vh,solver=CG);  // rijetka simetrična forma
   cout << "Matrica u Morseovom formatu, simetricna:"<<endl;
   cout << "AAA="<<AAA<<endl;


   matrix AAAA=a(Vh,Vh,solver=GMRES);  // rijetka nesimetrična forma
   cout << "Matrica u Morseovom formatu, nesimetricna:"<<endl;
   cout << "AAAA="<<AAAA<<endl;

//  Puna matrica: Kao i kod vektora dimenzija mora biti zadana.
   int N=8;
   real[int,int] punaA(N,N);  // Sintaksa. Matrica nije inicijalizirana.
   punaA=0;                   // Inicijalizacija nulama

   // Incijalizacija
   for(int i=0;i<N;++i){
       if(i-1 >= 0) punaA(i-1,i)=-1.0;
       punaA(i,i)=2.0;
       if(i+1<N) punaA(i+1,i)=-1.0;
   }

   cout<<"Primjer pune matrice:"<<endl;
   for (int i=0; i<N;++i) {
      for(int j=0; j<N;j++) cout<<punaA(i,j)<<" ";
      cout<<endl;
   }

   //cout << punaA <<endl; // jednostavnija verzija

  // pretvori punu matricu u rijetku
  matrix sparseA = punaA;

  cout<<"Prethodna matrica u rijetkom formatu:"<<endl;
  cout<< sparseA  << endl;