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

load "iovtk"

// zadavane parametry

real theta = 4.*pi/3.;
real a = 2.; //The length of the semimajor axis
real b = 1.; //The length of the semiminor axis
func z = x;

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

border Gamma1(t=0., theta){x=a*cos(t); y=b*sin(t);}
border Gamma2(t=theta, 2.*pi){x=a*cos(t); y=b*sin(t);}

// vykresleni hranice oblasti
// sipky naznacuji orientaci hranice

plot(Gamma1(100)+Gamma2(50),ps="membrana_hranice.ps");

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

mesh Th = buildmesh(Gamma1(100) + Gamma2(50));

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

fespace Vh(Th, P2); //P2 conforming triangular FEM

// reseni >phi<, testovaci funkce >w<
// a prava strana rovnice >f<
// na triangulovane oblasti >Th<

Vh phi, w, f=1;

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

real cpu=clock();

// definice slabe formulace pro >phi<
// a jeji reseni

solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w)+dy(phi)*dy(w))
                     -int2d(Th)(f*w)
                     +on(Gamma1,phi=z);

// ulozeni vysledku do vtk souboru

int[int] order=[1];
string dataname="phi";
savevtk("mambrana_phi.vtu",Th,phi,dataname=dataname,order=order);
plot(Th, wait=true, ps="membraneTh.eps"); //Plot Th

// vypis strojoveho casu

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