//-----[nacteni knihovny pro ukladani]-----
//-----[  vysledku ve >vtk< formatu  ]-----

load "iovtk";

//-----[zadavane parametry]-----

real[int] E=[2.1e5,1.1e5]; // <- Younguv modul [MPa]
real[int] nu=[0.28,0.3]; // <- Poisonova konstanta
real[int] alpha=[1.e-5,1.e-4]; // <- teplotni roztaznost [K^(-1)]
real[int] lambda(2), mu(2); // <- Lameho parametry
string nazev="prut.dat"; // <- nazev vystupniho souboru

//-----[dopocitane hodnoty]-----

for (int ii=0; ii<2; ii++) {
  lambda[ii]=E[ii]*nu[ii]/((1+nu[ii])*(1-2*nu[ii]));
  mu[ii]=E[ii]/(2*(1+nu[ii]));
 }

//-----[hranice oblasti]-----

int[int] Hrana=[0,1,2,3,4,5,6];

border Gamma1(t=0.,5.){x=t;y=0.;label=Hrana[1];};
border Gamma2(t=0.,5.){x=5.+t;y=0.;label=Hrana[1];};
border Gamma3(t=0.,1.){x=10.;y=t;label=Hrana[2];};
border Gamma4(t=0.,1.){x=10.;y=1.+t;label=Hrana[2];};
border Gamma5(t=0.,5.){x=10.-t;y=2.;label=Hrana[3];};
border Gamma6(t=0.,5.){x=5.-t;y=2.;label=Hrana[3];};
border Gamma7(t=0.,1.){x=0.;y=2.-t;label=Hrana[4];};
border Gamma8(t=0.,1.){x=0.;y=1.-t;label=Hrana[4];};
border Gamma9(t=0.,5.){x=t;y=1.;label=Hrana[5];};
border Gamma10(t=0.,5.){x=t+5.;y=1.;label=Hrana[5];};
border Gamma11(t=0.,1.){x=5.;y=t;label=Hrana[6];};
border Gamma12(t=0.,1.){x=5.;y=1.+t;label=Hrana[6];};

//-----[vykresleni hranice oblasti]-----

plot(Gamma1(10)+Gamma2(10)+Gamma3(2)+Gamma4(2)+Gamma5(10)
     +Gamma6(10)+Gamma7(2)+Gamma8(2)+Gamma9(-10)+Gamma10(-10)
     +Gamma11(2)+Gamma12(2),ps="prut_hranice.eps");

//-----[triangiangulovane oblasti dolni >Th1< a horni >Th2<]-----

mesh Th1=buildmesh(Gamma1(80)+Gamma11(120)+Gamma9(-80)+Gamma8(10));
mesh Th2=buildmesh(Gamma2(80)+Gamma3(10)+Gamma10(-80)+Gamma11(-120));
mesh Th3=buildmesh(Gamma10(80)+Gamma4(10)+Gamma5(80)+Gamma12(-120));
mesh Th4=buildmesh(Gamma9(80)+Gamma12(120)+Gamma6(80)+Gamma7(10));

//-----[precislovani hornich oblasti >Th3< a >Th4< z 0 na 1]-----

int[int] R1=[0,1];
Th3=change(Th3,region=R1);
Th4=change(Th4,region=R1);

//-----[sjednoceni obou triangulovanych oblasti]-----

mesh Th=Th1+Th2+Th3+Th4;

plot(Th,ps="prut_mesh.eps");

//-----[prostor >Vh< konecnych prvku definovanych na]-----
//-----[         triangulovane oblasti >Th<         ]-----

fespace Vh(Th, P2);
Vh u,v;
Vh uu,vv;
Vh theta=1; // <- změna teploty

//-----[definice deformaci a napeti]-----

real sqrt2=sqrt(2.);
macro epsilon(u1,u2) [dx(u1),dy(u2),(dy(u1)+dx(u2))/sqrt2] //
macro div(u,v) (dx(u)+dy(v)) //
macro sigma(u,v) [lambda[region]*(dx(u)+dy(v))+2.*mu[region]*dx(u)
                  -(3.*lambda[region]+2.*mu[region])*alpha[region]*theta,
                  lambda[region]*(dx(u)+dy(v))+2.*mu[region]*dy(v)
                  -(3.*lambda[region]+2.*mu[region])*alpha[region]*theta,
                  mu[region]*(dy(u)+dx(v))] //

//-----[reseni slabe formulace]-----

solve lame([u,v],[uu,vv])=
  int2d(Th)(2.*mu[region]*(epsilon(u,v)'*epsilon(uu,vv))
  +lambda[region]*div(u,v)*div(uu,vv))
  -int2d(Th)((3.*lambda[region]+2.*mu[region])*alpha[region]*div(uu,vv)*theta)
  +on(4,u=0,v=0);

//-----[projekce slozek napeti >sigma<]-----
//-----[   do prostoru polynomu >P1<  ]-----
//-----[na triangulovane oblasti >Th< ]-----

fespace Vh1(Th,P1);
Vh1 Shx=sigma(u,v)[0];
Vh1 Shy=sigma(u,v)[1];
Vh1 Shxy=sigma(u,v)[2];

//-----[vykrelseni posuvu do >PostScriptu<]-----

plot([u,v],ps="prut_uv.eps");

//-----[deformovana triangulovana oblast]-----

real coef=100.;
mesh Thmoved=movemesh(Th,[x+coef*u,y+coef*v]);

//-----[napasovani posuvu na deformovanou]-----
//-----[      triangulovanou oblast      ]-----
//-----[  a ulozeni do formatu ParaView  ]-----

string dataname="posuvy";
savevtk("prut_uv_moved.vtk",Thmoved,[u,v],dataname=dataname);

//-----[ulozeni posuvu ve formatu pro ParaView]-----

dataname="posuvy";
savevtk("prut_uv.vtk",Th,[u,v],dataname=dataname);

//-----[ulozeni posuvu ve formatu pro ParaView]-----

dataname="napeti";
savevtk("prut_sigma.vtk",Th,sigma(u,v),dataname=dataname);

//-----[ulozeni posuvu a napeti podel pricneho]-----
//-----[    prurezu prutu v bode >x=5.mm<     ]-----
//-----[          do souboru >nazev<          ]-----

int Nv=Th.nv;
real Eps=1.e-5;
real Xb,Yb,Ub,Vb;
real Sx,Sy,Sxy;

{
 ofstream soubor(nazev);
   soubor << "      x            y             u            v            sx           sy           sxy" 
          << endl;

   for (int ii=0; ii<Nv; ii++){
     Xb=Th(ii).x;
     Yb=Th(ii).y;
     Ub=u(Xb,Yb);
     Vb=v(Xb,Yb);
     Sx=Shx(Xb,Yb);
     Sy=Shy(Xb,Yb);
     Sxy=Shxy(Xb,Yb);
 
     if (Xb<5.+Eps && Xb>5.-Eps){
       soubor.scientific << Xb << " " << Yb 
                         << " " << Ub << " " << Vb 
                         << " " << Sx << " " << Sy << " " << Sxy 
                         << endl;
      }
    }
}