// nacteni knihovny pro ukladani vysledku ve >vtk< formatu

load "iovtk";

// zadavane parametry

int[int] hrana=[0,1,2,3,4,5,6]; // indexy casti hranice >gamma<
real E=2.1e5; // Younguv modul [MPa]
real nu=0.3; // Poisonova konstanta
real Mk=1.0e4; // moment zatizeni [N*mm]

// dopocitane parametry

real G=E/2/(1+nu); // modul pruznosti ve smyku [MPa]

// hranice meshovane oblasti (profil L),
// hranice musi byt orientovana tak, aby 
// oblast, kterou ohranicuje, lezela nalevo

border gamma1(t=0,10){x=0;y=-t;label=hrana[1];};
border gamma2(t=0,2){x=t;y=-10;label=hrana[2];};
border gamma3(t=0,8){x=2;y=-10+t;label=hrana[3];};
border gamma4(t=0,3){x=2+t;y=-2;label=hrana[4];};
border gamma5(t=0,2){x=5;y=-2+t;label=hrana[5];};
border gamma6(t=0,5){x=5-t;y=0;label=hrana[6];};

// vykresleni hranice oblasti
// sipky naznacuji orientaci hranice

plot(gamma1(10)+gamma2(2)+gamma3(8)+gamma4(3)+gamma5(2)+gamma6(5), wait=true, ps="krut_hranice.eps");

// Triangulace >Th< oblasti uvnitr hranice >gamma<

mesh Th=buildmesh(gamma1(40)+gamma2(8)+gamma3(32)+gamma4(12)+gamma5(8)+gamma6(20));

// postor >Vh< konecnych prvku definovanych na 
// triangulovane oblasti >Th<

fespace Vh(Th, P1);

// po castech P1 spojite reseni >phi1< a testovaci funkce >v<
// na triangulovane oblasti >Th<

Vh w,Chi,v;

// nastaveni hodin v sekundach
// pro mereni delky vypoctu

real cpu=clock();

// definice slabe formulace pro posuvy >w<
// a jeji reseni

solve Laplace(w,v,solver=LU)=int2d(Th)(dx(w)*dx(v)+dy(w)*dy(v)) // bilinearni cast
                             -int1d(Th,hrana[1])(-y*v) // Neumannova podminka na hrane[1]
                             -int1d(Th,hrana[2])(x*v) // Neumannova podminka na hrane[2]
                             -int1d(Th,hrana[3])(y*v) // Neumannova podminka na hrane[3]
                             -int1d(Th,hrana[4])(x*v) // Neumannova podminka na hrane[4]
                             -int1d(Th,hrana[5])(y*v) // Neumannova podminka na hrane[5]
                             -int1d(Th,hrana[6])(-x*v); // Neumannova podminka na hrane[6]

solve Poisson(Chi,v,solver=LU)=int2d(Th)(dx(Chi)*dx(v)+dy(Chi)*dy(v)) // bilinearni cast
                               -int2d(Th)(2.*v) // prava cast
                               +on(hrana[1],Chi=0) // Dirichletova podminka na hrane[1]
                               +on(hrana[2],Chi=0) // Dirichletova podminka na hrane[2]
                               +on(hrana[3],Chi=0) // Dirichletova podminka na hrane[3]
                               +on(hrana[4],Chi=0) // Dirichletova podminka na hrane[4]
                               +on(hrana[5],Chi=0) // Dirichletova podminka na hrane[5]
                               +on(hrana[6],Chi=0); // Dirichletova podminka na hrane[6]

// korekce polarniho momentu plochy >J<

real J=int2d(Th)(x*x+y*y-y*dx(w)+x*dy(w));
cout << "J=" << J << endl;

// dalsi mozny vypocet kv.  momentu plochy >J<

J=2.*int2d(Th)(Chi);
cout << "J=" << J << endl;

// vypocet >theta<

real theta=Mk/G/J;
cout << "theta=" << theta <<endl;

// vypocet napeti >sigma_zx< a >sigma_zy<

func sigmazx=G*theta*dy(Chi);
func sigmazy=-G*theta*dx(Chi);

// ulozeni vysledku do vtk souboru

int[int] order=[1];
string dataname="w [mm]";
savevtk("krut_w.vtu",Th,theta*w,dataname=dataname,order=order);
dataname="s_zx [MPa]";
savevtk("krut_szx.vtu",Th,sigmazx,dataname=dataname,order=order);
dataname="s_zy [MPa]";
savevtk("krut_szy.vtu",Th,sigmazy,dataname=dataname,order=order);
dataname="chi";
savevtk("krut_chi.vtu",Th,Chi,dataname=dataname,order=order);

// vypis strojoveho casu

cout << "CPU time = " << (clock()-cpu) << endl;