cout<<" \n\n Creating boundaries \n";
/*
polygon layout:
       _3_
    4/    \2
   5|      |1
    6\ ___/8
         7
sides identification
1-5
2-6
3-7
4-8
*/

border C1 (t = 5*pi/4, 3*pi/4) { x = sqrt((sqrt(2)-1)/2)*cos(t) + sqrt((sqrt(2)+1)/2)  ;
				                 y = sqrt((sqrt(2)-1)/2)*sin(t)  ; label=1    ;}

border C2 (t = 3*pi/2, pi)     { x = sqrt((sqrt(2)-1)/2)*cos(t) +  sqrt(2)/2*sqrt((sqrt(2)+1)/2)  ; 	 
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) +  sqrt(2)/2*sqrt((sqrt(2)+1)/2) ; label=2;}

border C3 (t = 7*pi/4, 5*pi/4) { x = sqrt((sqrt(2)-1)/2)*cos(t)   ;  
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) + sqrt((sqrt(2)+1)/2); label=3  ;}

border C4 (t = 2*pi, 3*pi/2)   { x = sqrt((sqrt(2)-1)/2)*cos(t) - sqrt(2)/2*sqrt((sqrt(2)+1)/2)  ;  
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) + sqrt(2)/2*sqrt((sqrt(2)+1)/2) ; label=4 ;}

border C5 (t = 9*pi/4, 7*pi/4) { x = sqrt((sqrt(2)-1)/2)*cos(t) - sqrt((sqrt(2)+1)/2)  ; 
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) ; label=5  ;}

border C6 (t = pi/2, 0)        { x = sqrt((sqrt(2)-1)/2)*cos(t) - sqrt(2)/2*sqrt((sqrt(2)+1)/2)  ;  
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) -sqrt(2)/2*sqrt((sqrt(2)+1)/2) ; label=6 ;}

border C7 (t = 3*pi/4, pi/4)   { x = sqrt((sqrt(2)-1)/2)*cos(t)   ;  
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) -  sqrt((sqrt(2)+1)/2) ; label=7 ;}

border C8 (t = pi, pi/2)       { x = sqrt((sqrt(2)-1)/2)*cos(t) + sqrt(2)/2*sqrt((sqrt(2)+1)/2)  ;	 
				                 y = sqrt((sqrt(2)-1)/2)*sin(t) - sqrt(2)/2*sqrt((sqrt(2)+1)/2) ; label=8 ;}

cout<<"\n\nBuilding mesh\n";

int vertnum = 30;
mesh Th = buildmesh(C1(vertnum)+C2(vertnum)+C3(vertnum)+C4(vertnum)+C5(vertnum)+C6(vertnum)+C7(vertnum)+C8(vertnum),fixeborder=true);

//Adapting the mesh messes up the periodic boundary conditions
//func metric = 4/(1-x^2-y^2)^2;
//mesh ATh = adaptmesh(Th,metric);

cout<<"\n\n Creating function space \n";
real side = sqrt(3/4*2^(7/8)+2^(3/8)+1/2);
//dist15 = (x,y).(0,side)
func dist15 = y * side;
//dist26 = (x,y).(-sqrt(2)*side,sqrt(2)*side)
func dist26 = (y-x)*sqrt(2)*side;
//dist37 = (x,y).(side,0)
func dist37 = x * side;
//dist48 = (x,y).(sqrt(2)*side,sqrt(2)*side)
func dist48 = (y+x)*sqrt(2)*side;
func perio = [[1,dist15],[5,dist15],[2,dist26],[6,dist26],[3,dist37],[7,dist37],[4,dist48],[8,dist48]];

//function space with periodic boundary conditions and the solution of \nabla ^2 u=-xy
fespace Vhp(Th,P1, periodic=perio);
Vhp up,vp;
func g = x*y;
solve Poissonp(up,vp,solver=LU) = int2d(Th)(dx(up)*dx(vp)+dy(up)*dy(vp))-int2d(Th)(g*vp);
plot(up,wait=1);

//function space without periodic boundary conditions and the solution of \nabla^2 u=-xy with u=0 at the boundaries:
fespace Vh(Th,P1);
Vh u,v;
func f = x*y;
solve Poisson(u,v,solver=LU) = int2d(Th)(dx(u)*dx(v)+dy(u)*dy(v))-int2d(Th)(f*v)+on(C1,C2,C3,C4,C5,C6,C7,C8,u=0);
plot(u,wait=1);