verbosity=0;

//load "UMFPACK64"


real cputotal=clock();

macro Grad3(u) [dx(u),dy(u),dz(u)]  // EOM

//
// read the 3D mesh 
//
mesh3  Th3("Th3.mesh");

fespace Vh(Th3,P13d);

cout<<"DOFs "<<Vh.ndof<<endl;
Vh u,v,fh,qh, one=1;
real volu=int3d(Th3)(one);
real invvol=1/volu;

//
// Parameter k0 in [k0min,k0max]
// kleft, kright y Bi funciones de k0
//

real k0,k0min,k0max,kleft,kright,Bi,aa;
include "data.idp"
k0=k0min;
//
// Other parameters depend on k0
//
include "parameters.idp";

//
// variational forms free from parameters. To be computed once  
//
varf a1(u,v)= int3d(Th3,0,2,4,6,8,10)(Grad3(u)' *Grad3(v));// goes with factor k0
varf a2(u,v)= int3d(Th3,1,5,9)(Grad3(u)' *Grad3(v));// goes with factor kleft
varf a3(u,v)= int3d(Th3,3,7,11)(Grad3(u)' *Grad3(v));// goes with factor kright
varf a4(u,v)= int3d(Th3)( u*v ) ;// goes with factor aa
varf a5(u,v)=int2d(Th3,10,12,14,16,18,110,11,15,19,13,17,111
                          ,30,32,34,36,38,310,31,35,39,33,37,311,
                         40,44,48,41,45,49,20,24,28,23,27,211,199,201
                          ,299,301,399,401,499,501,599,601,699,700,701)
                          (u*v);// goes with factor Bi


varf b1(unused,v)= int3d(Th3)( fh*v )+int2d(Th3,100)( qh*v ); 

matrix ma1= a1(Vh,Vh);

matrix ma2= a2(Vh,Vh);
matrix ma3= a3(Vh,Vh);
matrix ma4= a4(Vh,Vh);
matrix ma5= a5(Vh,Vh);

Vh RHSunico;
RHSunico[]=b1(0,Vh);
//
// GLOBAL MATRIX DEPENDS ON PARAMETERS
//
matrix Aglobal,m1,m2,m3,m4,m5;
//         
// SET A COLOR SCALE FOR DRAWING PLOTS (gracias Gladys)
//
real[int] colorhsv=[
4./6.,1, 0.5,		//dark blue
4./6., 1, 1,		//blue
0.5,1,1,			//cyan
0.44, 1, 1,		//light green
0.33,1,1,		//green
0.2, 1, 1,		//light yellow
0.15,1,1,		//yellow
0.12,1,1,		//orange
0.09, 1, 1,			//light red
0.027, 1, 1			//red
];
int samples;
ifstream file("RBdim.txt");
	file>>samples;

int nn;
ifstream fileth("Meshsize.txt");
	fileth>>nn;
//
// COPY OF THE TRIANGULATION FOR COMPARISON PURPOSES
//	
mesh3  Th3c("Th3copy.mesh");
fespace Vhc(Th3c,P13d);

Vh[int] usample(samples);
//
// MATRICES FOR THE REDUCED PROBLEM
//
real[int,int] A0(samples,samples); // matriz
real[int,int] A1(samples,samples); // matriz
real[int,int] A2(samples,samples); // matriz
real[int,int] A3(samples,samples); // matriz
real[int,int] A4(samples,samples); // matriz

real[int,int] Ro(samples,samples); // matriz
real[int] Rrhsso(samples); // rhs
//
// READ DATA FROM DISK
//
cout<<"READING BASES REDUCIDAS dim = "<<samples<<endl;
for(int i=0;i<samples;i++)
{
	ifstream file("sampleBASE"+i+".txt");
	file>>usample[i][];
}

ifstream file0("A0.txt");
ifstream file1("A1.txt");
ifstream file2("A2.txt");
ifstream file3("A3.txt");
ifstream file4("A4.txt");
for(int i=0;i<samples;i++)
{
for(int j=0;j<samples;j++)
 {
   file0>>A0(i,j);
   file1>>A1(i,j);
   file2>>A2(i,j);
   file3>>A3(i,j);
   file4>>A4(i,j);
  }
}
ifstream Rrhsfile("Rrhs.txt");
for(int j=0;j<samples;j++)
 {
   Rrhsfile>>Rrhsso[j];
  }

cout<<"READING DONE  "<<endl;
cout<<"RB dim = "<<samples<<endl;