real totaltime=clock();
include "settingANDproblem.idp"
//
// FIRST SAMPLE COMPUTED WITH k0=k0min GLOBAL SOLVE ON LARGE FEM SPACE 
//
include "ParametersAndSolve.idp";
//
// SET ISOLINES TO PLOT ACCORDING TO COMPUTED SOLUTION
//
 int mm=40;
real[int] viso(mm+1);
for (int i=0;i<viso.n;i++)
viso[i]=u[].min+i*(u[].max-u[].min)/mm;

plot(u,cmm=" SOLUTION FOR k0  = "+k0,wait=0,value=1,fill=1,viso=viso(0:viso.n-1),hsv=colorhsv);

//-------------------------------------------->
//    UNCOMMENT TO SEE THAT IS FASTER WORKING WITH MATRICES
//
//     cpu=clock();// get current cpu time
//	  Lap3; // solve no adding matrices
//	tttt=clock()-cpu;
//	
//cout << " CPU time NO ADDING MATRICES = " << clock()-cpu << endl;
////
//// SET ISOLINES TO PLOT ACCORDING TO COMPUTED SOLUTION
////
//for (int i=0;i<viso.n;i++)
//viso[i]=u[].min+i*(u[].max-u[].min)/mm;
//
//plot(u,cmm="NO ADDING MATRICES: TIME "+tttt+ " --> 
//SAMPLE FOR k0  = "+k0,wait=1,value=1,fill=1,viso=viso(0:viso.n-1),hsv=colorhsv);
//
//------------------------------------------------->

//
// OUTPUT FOR FIRST FUNCTION
//
	  cc[0]=int2d(Th3,700)(u);


cout<<" FOR k0 = "<<k0<<"; Average of U on top = "<<cc[0]<<endl;

//
// H1 norm of the solution
//
real part1=int3d(Th3)(dx(u)^2+dy(u)^2+dz(u)^2);
real part2=int3d(Th3)(u^2);
real normu=sqrt(part1+invvol*part2);
//
// FIRST REDUCED BASIS FUNCTION
//
      usample[0]=u/normu;
//              
cout <<"FUNCION DE BASE REDUCIDA  0 HECHA / REDUCED BASIS 0 DONE "<<endl;
ofstream file("sampleBASE"+0+".txt");
file<<usample[0][];
//
// KEEP PARAMETER AND INDEX ASSOCIATED
//
mui[0]=k0min;
jbest[0]=0; 
//
// LOOP COMPUTATION OF usample[j] FOR j=0,...,n-1
//
int n=0;

for (n=1;n<samples;n++)// LOOKING FOR FUNCTIONS OF REDUCED BASIS
{
//
// SO FAR WE HAVE THE  0,1,2,...,n-1 FUNCTIONS
//
cout <<"NUMBER OF BASIS FUNCTIONS SO FAR = "<<n<<" ( STARTING FROM PHI_0 ) "<<endl;
cout <<"CONSTRUCT ALL PARAMETER FREE MATRICES AND VECTOR FOR REDUCED BASIS... "<<endl;
// 
// CONSTRUCT ALL PARAMETER FREE MATRICES AND VECTOR FOR REDUCED BASIS
//
for(int i=0;i<n;i++)
{
// 
// Computational work could be saved here. Just add  terms comming from new basis function
//
for(int j=0;j<n;j++)
         {        
real aux0=int3d(Th3,0,2,4,6,8,10)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    );  
real aux1=int3d(Th3,1,5,9)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    );
real aux2=int3d(Th3,3,7,11)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    );   

 real auxtest=int3d(Th3)(usample[i]*usample[j] );   


real aux3=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)
                   (usample[i]*usample[j]);
      A0(i,j)=aux0;
      A1(i,j)=aux1;
      A2(i,j)=aux2;
      A3(i,j)=aux3;
      A4(i,j)=auxtest;       
         }
}
// 
// Ahora el lado derecho
//
for(int i=0;i<n;i++)
{
Rrhs[i]=int3d(Th3)(fh*usample[i])+int2d(Th3,100)(qh*usample[i]);
}

//
// PARAMETER SEARCH ONLY USES REDUCED BASIS AND PARAMETER FREE COMPUTATIONS
// ARE ALREADY DONE
//
include "ParameterSearch.idp";// 
//
//
cout <<"OUT OF PARAMETER SEARCH. OPTIMAL ONE = "<<bestk0<<endl;
//
// OBTENEMOS LA SOLUCION FEM PARA EL MU OPTIMO
//   	
       k0=bestk0;
      cpu=clock();
      include "ParametersAndSolve.idp";

	  cout << " CPU time = " << clock()-cpu << endl;
	  cc[n]=int2d(Th3,700)(u);
	  cout<<" Bestk0 = "<<bestk0
	                    <<"; Average of U on top = "<<cc[n]<<endl;

// plot(u,cmm="SAMPLE PARA mu optimo = "+mu+ ". Average of U on top =  "+
//int2d(Th3,200)(u),wait=0,nbiso=nbiso,value=1,fill=1);             
//
//  USE OF GRAM-SCHMID TO ORTOGONALIZE u WITH RESPECT TO 
//  THE REST OF ELEMENTS IN THE REDUCED BASIS 
//  AND OBTAIN A ORTHOGONAL BASIS  
//
Vh ugs=u;
for (int j=0;j<n;j++)
{
real p1=int3d(Th3)(dx(u)*dx(usample[j])+dy(u)*dy(usample[j])+dz(u)*dz(usample[j])
                  );
real p2=int3d(Th3)(u*usample[j]);
real prod=p1+invvol*p2;

ugs=ugs-prod*usample[j];
}
//
// COMPUTE NORM
//
real part1= int3d(Th3)(dx(ugs)^2+dy(ugs)^2+dz(ugs)^2 );
real part2= int3d(Th3)(ugs^2);
real normugs=sqrt(part1+part2*invvol);
//
// TESTING ORTHOGONALITY wrt 0,1,2,...,n-1
//
real orto=0;
for (int j=0;j<n;j++)
{ // SCALAR PRODUCT ugs vs USAMPLE[j] for j=0,1,2,..,n-1
real aux0=int3d(Th3)( dx(ugs)*dx(usample[j])
                       +dy(ugs)*dy(usample[j])
                       +dz(ugs)*dz(usample[j])                      
                      );                          
real auxnot=int3d(Th3)(ugs*usample[j]);

      real ortoij=(aux0+invvol*auxnot)/normugs;
      //
 cout<<"REDUCED BASIS: Orthogonality n-th basis = "<<n<<"  versus j-th = "<<j<<" is "<<ortoij<<endl; 
      orto=orto+ abs(ortoij); // ADD ALL ORTHOGONALITIES   
}
cout<<"REDUCED BASIS FUNCION  =  "<<n<<
" DONE (starting from Psi_0). TOTAL ORTONORMALITY = "<<orto<<endl;
//
// STORE NEW FUNCTION
//
     usample[n]=ugs/normugs; 
//
// STORE NEW PARAMETER  
//          
      mui[n]=bestk0;

      cout <<"ALL CHOSEN PARAMETERS "<<mui<<endl;
      cout <<"INDEX FOR THIS PARAMETERS "<<jbest<<endl;

  if(abs(orto)>1e-4)// ORTHOGONALITY LOSS NO MORE COMPUTING
{
	cout<<"ORTOGONALITY LOSS WITH BASIS FUNCTION = "<<n<<". GIVEN BY "<<orto<<endl;

samples=n;// DO NOT CONSIDER THE NEWLY COMPUTED FUNCTION

	break;
} 

}
//
// END SEARCH FOR REDUCED FUNCTIONS
//
// CHOSEN PARAMETERS
//
cout<<"Parameters = "<<mui<<endl;      
//
// PARAMETER FREE MATRICES AND VECTORS WITH THE REDUCED BASIS
//
cout<<"COMPUTING MATRIX FOR REDUCED BASIS wait...."<<endl; 
for(int i=0;i<samples;i++)
{
for(int j=0;j<samples;j++)
         {  
real aux0=int3d(Th3,0,2,4,6,8,10)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    );  
real aux1=int3d(Th3,1,5,9)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    );
real aux2=int3d(Th3,3,7,11)(
                        dx(usample[i])*dx(usample[j])
                       +dy(usample[i])*dy(usample[j])
                       +dz(usample[i])*dz(usample[j])
                    ); 
real auxtest=int3d(Th3)(usample[i]*usample[j] );   


real   aux3=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)(usample[i]*usample[j]);

      A0(i,j)=aux0;
      A1(i,j)=aux1;
      A2(i,j)=aux2;
      A3(i,j)=aux3;
      A4(i,j)=auxtest;

         }
}
//
// SAVE ON DISK ALL FUNCTIONS AND DATA
//
cout<<"SIZE OF REDUCED BASIS  = "<<samples<<endl;
ofstream fileNB("RBdim.txt");
fileNB<<samples;
//
// SAVE ON DISK BASIS FUNCTIONS
//
for(int i=0;i<samples;i++)
{
	ofstream file("sampleBASE"+i+".txt");
	file<<usample[i][];
}
//
// SAVE ON DISK PARAMETER FREE MATRICES AND VECTORS WITH 
// THE REDUCED BASIS
//

ofstream file0("A0.txt");
ofstream file1("A1.txt");
ofstream file2("A2.txt");
ofstream file3("A3.txt");
ofstream file4("A4.txt"); 
for(int i=0;i<samples;i++)
{
for(int j=0;j<samples;j++)
 {
   file0<<A0(i,j)<<endl;
   file1<<A1(i,j)<<endl;
   file2<<A2(i,j)<<endl;
   file3<<A3(i,j)<<endl;
   file4<<A4(i,j)<<endl;
  }
}

// 
// SAVE ON DISK RHS
//

for(int i=0;i<samples;i++)
{
Rrhs[i]=int3d(Th3)(fh*usample[i])+int2d(Th3,100)(qh*usample[i]);
}
ofstream Rrhsfile("Rrhs.txt");
for(int j=0;j<samples;j++)
Rrhsfile<<Rrhs[j]<<endl;

cout << " CPU TOTAL TIME = " << clock()-totaltime << endl;