clear

//data
Lx=1;

K = 1;


//space discretization
N=10;
Nint = N-1;
Nbound = 2;

dx=Lx/N;

//data structure for the uniform 1D mesh of size dx of the domain (0,Lx)
// cells m = 1:N
for m=1:N
  volume(m) = dx;
  X(m) = (m-1./2)*dx;
end  

//Interior faces: isig = 1:N-1
for isig=1:Nint
 // cellint(isig,1) = 
 // cellint(isig,2) = 
 // surfaceint(isig) = 
 // Xint(isig) = 
end

//Boundary faces: isig = 1:Nbound

// left end
  isig = 1;
 // cellbound(isig) = 
 // surfacebound(isig) = 
  Xbound(isig) = 0;

// right end  
  isig = 2;
 // cellbound(isig) = 
 // surfacebound(isig) =  
    Xbound(isig) = Lx;


//transmissibilities of interior faces
for isig=1:Nint
//  m1 = cellint(isig,1);
//  m2 = cellint(isig,2);
  
//  Tint(isig) = 
  
end

//transmissibilities of boundary faces
for isig=1:Nbound
  
//  m = cellbound(isig);
  
//  Tbound(isig) = 

end


// exact solution

pi = 4*atan(1);
for m=1:N
  pe(m) = sin(pi*X(m));
end

// right hand side = -K*laplacian(pe) 
for m=1:N
  sm(m) = K*pi**2*sin(pi*X(m));
end

// Dirichlet Boundary condition
for isig=1:Nbound
 pbound(isig) = sin(pi*Xbound(isig));
end

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

//Jacobian matrix and right hand side
A=sparse([],[],[N,N]);

// source term
for m=1:N
  B(m,1) = 0;
end 
 
 
 //Loop on boundary faces
 for isig=1:Nbound
//    m = cellbound(isig);

//    A(m,m) = A(m,m) + ;
//    B(m,1) = B(m,1) + ;
 end

 //Loop on interior faces
 for isig=1:Nint
//  m1 = cellint(isig,1);
//  m2 = cellint(isig,2);
  
 
  
//  A(m1,m1) = A(m1,m1) + ;
//  A(m2,m2) = A(m2,m2) + ; 
//  A(m1,m2) = A(m1,m2) - ;
//  A(m2,m1) = A(m2,m1) - ;

end



 //solution 
 p=A\B;
 
 
 //plot
 scf(1);
 xtitle('error')
 plot (X,p-pe,'r-*')
 
 scf(2);
 xtitle('exact (red) and approximate (blue) solution') 
 plot (X,p,'b-*')
 
 scf(2);
 plot (X,pe,'r-*')
 





  
  




