//  august 2015 add 
//  DEF LABEL unit CUBE before transformation 
//  6 face plan label are:
//  y = 0 ) label 1, 
//  x = 1 ) label 2, 
//  y = 1 ) label 3, 
//  x = 0 ) label 4, 
//  z = 0 ) label 5,
//  z = 1 ) label 6.
// the region numbering is 0 by default.
{
	real[int,int] mm= [ [0.5,0.0,0.5], [1.0,0.5,0.5], 
	                    [0.5,1.0,0.5], [0.0,0.5,0.5], 
						[0.5,0.5,0.0], [0.5,0.5,1.0] ];
	int[int] l6=[1,2,3,4,5,6];
	mesh3 Th=cube(2,2,2); //[x*4,y*1+0.1*sin(4*pi*x)]);
	for (int i=1; i< 7; ++i)
	{
		real s = int2d(Th,i)(1.); 
		real mx = int2d(Th,i)(x); 
		real my = int2d(Th,i)(y); 
		real mz = int2d(Th,i)(z); 	
		cout<< i << " " << s <<  " "<< mx << " " << my << " " << mz << endl;
		int j=i-1;
		real err= abs(mx-mm(j,0)) + abs(my-mm(j,1)) + abs(mz-mm(j,2)) ;
		assert(err< 0.1);
	}
}
int[int] l6=[37,42,45,40,25,57];
int r11=11;
mesh3 Th=cube(4,5,6,[x*2-1,y*2-1,z*2-1],label=l6,flags =3,region=r11); 
cout << " volume " << Th.measure << ", border area "<< Th.bordermeasure <<endl; // \index{mesh3!measure} \index{mesh3!bordermeasure}

// Check label dans  region numbering 
int err =0; 
for(int i=0; i<100; ++i)
{
    real s =int2d(Th,i)(1.);
    real sx=int2d(Th,i)(x);
    real sy=int2d(Th,i)(y);
    real sz=int2d(Th,i)(z);
    
    if( s )
    {
     int ix = (sx/s+1.5), iy=(sy/s+1.5), iz=(sz/s+1.5), ii=(ix + 4*(iy+1) + 16*(iz+1) ) ;    
      //  value of ix,iy,iz =>  face min 0 ,  face max 2  , no face 1
      cout <<" label="<< i << " s " << s << " " << ix << iy << iz << " : " << ii << endl; 
      if( i != ii ) err++;
    }
}   
real volr11 = int3d(Th,r11)(1.) ;
cout << " vol region " << 11 << ": " << volr11 << endl; 
if( (volr11 - Th.measure )>1e-8) err++;
plot(Th,fill=0); 
cout << " nb err= " << err <<endl;
assert(err==0); 
