clear

exec 'func_well_solution.sci'


//data
Lx = 100;
Ly = 100;

//productor well pressure 
pw = 0.;

//well radius
rw = 0.1


//productor well flow rate of the exact solution
qw = 1.


//permeability
K=1.;



//space discretization: Nx and Ny (odd to have the well at the center)
Nx = 11;
Ny = Nx;
N = Nx*Ny;
dx = Lx/Nx;
dy = Ly/Ny;

//Nint = 
//Nbound = 


//uniforme mesh of size dx dy 
// cells m = 1:N
for i=1:Nx
  Xi(i) = (i-1/2)*dx;
for j=1:Ny  
  Yj(j) = (j-1/2)*dy;
  
  m = i + (j-1)*Nx;
  
  volume(m) = dx*dy;
  X(m) = Xi(i);
  Y(m) = Yj(j);
end
end

//Interior faces: isig = 1:Nint
isig = 0;
//normal x 
for i=1:Nx-1
for j=1:Ny

  isig = isig + 1;
  
//  cellint(isig,1) = 
//  cellint(isig,2) = 
//  surfaceint(isig) = 
//  Xint(isig) = 
//  Yint(isig) = 
end
end 
//normal y
for i=1:Nx
for j=1:Ny-1

  isig = isig + 1;
  
//  cellint(isig,1) =
//  cellint(isig,2) =
//  surfaceint(isig) =
//  Xint(isig) =
//  Yint(isig) =
end
end 


//boundary faces x = 0
for j=1:Ny
  isig = j;
  i=1;
  m1 = i+(j-1)*Nx;
  
//  cellbound(isig) = 
//  surfacebound(isig) = 
//  Xbound(isig) = 
//  Ybound(isig) = 
end

//boundary faces x = Lx
for j=1:Ny
  isig = j+Ny;
  i=Nx;
  m1 = i+(j-1)*Nx;
  
//  cellbound(isig) =
//  surfacebound(isig) = 
//  Xbound(isig) = 
//  Ybound(isig) = 
end

//boundary faces y=0
for i=1:Nx
  isig = i+2*Ny;
  j=1;
  m1 = i+(j-1)*Nx;
  
//  cellbound(isig) = 
//  surfacebound(isig) = 
//  Xbound(isig) = 
//  Ybound(isig) = 
end

//boundary faces y=Ly
for i=1:Nx
  isig = i+Nx+2*Ny;
  j=Ny;
  m1 = i+(j-1)*Nx;
  
//  cellbound(isig) = 
//  surfacebound(isig) = 
//  Xbound(isig) = 
//  Ybound(isig) = 
end


for m=1:N
  //permeabilities
    perm(m) = K;
end

//transmissibilities of interior faces 
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  
  d1 = (Xint(isig)-X(m1))**2 + (Yint(isig)-Y(m1))**2;
  d1 = sqrt(d1);
  d2 = (Xint(isig)-X(m2))**2 + (Yint(isig)-Y(m2))**2;
  d2 = sqrt(d2);
  
  
//  Tint(isig) = 

end



//transmissibilities of boundary faces
for isig=1:Nbound
  m1 = cellbound(isig);
  d = (Xbound(isig)-X(m1))**2 + (Ybound(isig)-Y(m1))**2;
  d = sqrt(d);
  
//  Tbound(isig) = 

end



//well 

Nbwell = 1; //number of wells

//productor well

mw = int((Nx+1)/2) + (int((Ny+1)/2)-1)*Nx;
cellwell(1) = mw;
pi = 4*atan(1);


// Twell(1) = 




//boundary faces pressure (DIRICHLET BC)
for isig=1:Nbound
  
 r = sqrt((Xbound(isig)-X(mw))**2  + (Ybound(isig)-Y(mw))**2 );
 pressionbound(isig) = pwell_ex(pw,qw,rw,K,r);
 
end


//----------------------------------


  
// pressure computation 
//Jacobian matrix and right hand side 
A=sparse([],[],[N,N]);
 
for m=1:N
  B(m,1) = 0;
end

 //Loop on interior faces
for isig=1:Nint
  m1 = cellint(isig,1);
  m2 = cellint(isig,2);
  
  s = Tint(isig);
  
end

 //Loop on boundary faces
for isig=1:Nbound
  m = cellbound(isig);  
  s = Tbound(isig);
  
end



 //productor well: imposed pressure
 mw = cellwell(1);
 s = Twell(1);
 
 

 

 //solution
 p=A\B;
 

//
//  
//for i=1:Nx
//for j=1:Ny  
// m = i + (j-1)*Nx; 
// r = sqrt((X(m)-X(mw))**2  + (Y(m)-Y(mw))**2 );
// if (m~=mw)
//   pe(m) = pwell_ex(pw,qw,rw,K,r);
//  exy(i,j) = p(m)-pe(m);
//  pxy(i,j) = p(m);
//  pexy(i,j) = pe(m);
// else 
//  pe(m) = pwell_ex(pw,qw,rw,K,r0);
// // exy(i,j) = p(m)-pe(m);
//   exy(i,j) = 0;
//  pxy(i,j) = p(m);
// end
//end  
//end
//
//
//
//
//scf(1);
//xtitle('exact pression')
//plot3d(Xi/100,Yj/100,pexy,40,80)
//  
//scf(2);
//xtitle('approximate pressure')
//plot3d(Xi/100,Yj/100,pxy,40,80)
//
//scf(3);
//xtitle('error on pressure')
//plot3d(Xi/100,Yj/100,exy,40,80)
//
//
//
////relative l2 norm of the error 
//s = 0;
//ss = 0;
//for m=1:N
//  s = s + volume(m)*(p(m)-pe(m))**2;
//  ss = ss + volume(m)*pe(m)**2;
//end  
//
//erreur = sqrt(s)/sqrt(ss)
//
//
//print(%io(2), 'erreur ',erreur)
//
  
  
