//////////////////////////////////////
//                                  //
// Problem of a periodic cell       //
// MAP562     		   	    //
// Winter 2016/2017		    //
// G. Allaire, T. Wick		    //
//                                  //
//////////////////////////////////////

// In this program, we solve
// a cell problem; namely a composite
// material (two phases) with the coefficients 
// alpha and beta

// Tasks:
// 1. Compute the volume fraction of each phase
//   and check that the total volume is indeed 1
// 2. Check the properties of the homogenized 
//    tensor A* (Hint: Symmetry. Why?)
// 3. Implement the entry A21
// 4. Check the lower and upper bounds of 
//   arithmetic and harmonic type
// 5. Check that also the bounds of Hashin-Shtrikman
//   are satisfied.
// 6. Choose a non-symmetric inner cell and 
//   investigate the properties of the 
//   homogenized tensor A* (can you find mixtures 
//   of the two phases for which A* is close to 
//   the Hashin-Shtrikman bound ?)
/////////////////////////////////////
// Material parameters, "two-phases"
// The two parameters are scalar,
// thus the two phases are isotropic.
real alpha=1.;
real beta=10.;

// Coefficients of the homogenized tensor A*
real A11, A12, A21, A22 ;
string legende ;

//////////////////////////////// 			
// Definition of the domain   // 	
////////////////////////////////	
mesh Sh;

real pi=4*atan(1) ;

//  Top boundary
border a1(t=1,0) { x=t; y=1;label=1; };
// Left boundary
border a2(t=1,0)  { x=0; y=t;label=2; };
// Bottom boundary
border a3(t=0,1) { x=t; y=0;label=3; };
// Right boundary
border a4(t=0,1) { x=1; y=t;label=4; };
// Inclusion (ellipse)
border a5(t=0,2*pi) { x=0.5+0.1*cos(t); y=0.5+0.3*sin(t);label=5; };

//////////////////////////////
// Mesh			    //
//////////////////////////////
int n=8 ;
Sh= buildmesh (a1(5*n)+a2(5*n)+a3(5*n)+a4(5*n)+a5(10*n));
//mesh Th = Sh ;

plot(Sh,wait=1,ps="mesh1.eps"); 


//////////////////////////////
// Define FEM spaces
fespace Vh1(Sh,P1,periodic=[ [1,x],[3,x],[2,y],[4,y] ]);

// Solution variables: u1 and u2
// Test function: phi
Vh1 u1,u2,phi;	

// Definition of the P0 space for the coefficient
// Furthermore we define the characteristic function
fespace Vh0(Sh,P0);
Vh0 Ay , reg=region;	// Region labels

//////////////////////////////
// Definition of the coefficients for each composite
// as shown in the lecture notes
Ay = alpha*reg + beta*(1.-reg) ;			
plot(Sh,Ay,fill=1,wait=1,value=true); 



///////////////////////////
// Defining the equations//
///////////////////////////

// -div_y (A(y)(e_i + \nabla_y w_i(y))) = 0
// y -> w_i(y) is Y-periodic

problem cell1(u1,phi) =
    int2d(Sh)(
               Ay*(dx(u1)*dx(phi)+dy(u1)*dy(phi))
	)
    +int2d(Sh)( Ay*dx(phi) )
;

problem cell2(u2,phi) =
    int2d(Sh)(
               Ay*(dx(u2)*dx(phi)+dy(u2)*dy(phi))
	)
    +int2d(Sh)( Ay*dy(phi) )
;

///////////////////////////
// Solving the systems   //
///////////////////////////
cell1 ;
cell2 ;


/////////////////////////////////////////
// Evaluating quantities of interest   //
/////////////////////////////////////////
// Compute the entries of the homogenized tensor A*
A11 = int2d(Sh)( Ay*((dx(u1)+1.)*(dx(u1)+1.)+dy(u1)*dy(u1)) ) ;

A22 = int2d(Sh)( Ay*((dy(u2)+1.)*(dy(u2)+1.)+dx(u2)*dx(u2)) ) ;

A12 = int2d(Sh)( Ay*((dx(u1)+1.)*dx(u2) + dy(u1)*(dy(u2)+1.)) ) ;


/////////////////////////////////////////
// Output and visualization            //
/////////////////////////////////////////
cout<<"Homogenized coefficients, A11 = " << A11
    <<" A12 = " << A12
    <<" A21 = " << A21
    <<" A22 = "  <<A22
    <<endl;

legende="Cell 1 problem, A11 "+A11;
plot(Sh,u1,fill=1,wait=1,value=true,cmm=legende); 

legende="Cell 2 problem, A22 "+A22;
plot(Sh,u2,fill=1,wait=1,value=true,cmm=legende);