function elliptic (region,nx,ny,lx,ly,ax,bx,ay,by,flagexact)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ELLIPTIC SOLVER - node centered 
%% elliptic_solver (region,nx,ny,lx,ly,ax,bx,ay,by,flagexact)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %% M. Peszynska, 04/2005
%% © Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% solve an elliptic problem in a  region 
%% of type region (generated by numgrid) with
%% size determined by nx,ny which is the number of subintervals
%% 
%% numbering of nodes is 1..nx+1, 1..ny+1
%% internal nodes are 2: nx, 2:ny
% algorithm: 
% 1. generate_ijk_grid
% 2. number grid
% 3. generate matrix 
% 4. compute rhs 
% 5. incorporate boundary conditions: both matrix and rhs
% 6. solve linear system
% 7. postprocess the results 
%% flagexact = 1 means exact solution is known and display and error
%% calculation is desired
%%%%%%%
% functions used: 
% generate_ijk_grid
% node_index (sgrid)
%%%%%%%%
% examples:
% elliptic_solver('F',100,100,40,40,0,1,0,1,0)
% elliptic_solver('Holes',100,100,10,0,0,1,0,1,0)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% input or generate data

%% sgrid is 1 for interior and 0 for everybody else 
  sgrid = generate_ijk_grid(region,nx,ny,lx,ly);
%clf;pcolor(sgrid'),shading interp;colorbar;caxis([-1 1]);title('Numerical grid');axis equal;pause;
%% grid is numbered for interior and 0 for everybody else 
[grid,ind] = node_index(sgrid,'natural');

%% physical domain: generate (uniform) grid
if nargin <=5
     ax = 0; bx = 1; ay = 0; by =1;
end;

dx = (bx-ax)/(nx); dy = (by-ay)/ (ny);
x = ax:dx:bx; y = ay:dy:by;

%% read, set, or generate perms and velocities
perm1 = ones(nx+1,ny+1);perm2 = perm1;


%% Reynolds number
scale = 0;
%scale = 1e+2;
vely = rand(nx,ny)*scale; velx =rand(nx,ny)*scale; % interesting
%vely = ones(nx,ny)*scale; velx =ones(nx,ny)*scale;
%velx = zeros(nx,ny); vely = zeros(nx,ny);
fprintf('Peclet number*h = %g\n',max(max(velx(:)),max(vely(:)))*dx);

%% ! assign boundary conditions: assign -1 to Neumann nodes and 0 to
%% Dirichlet nodes 
fprintf('generating boundary conditions ..\n');
sgrid = assign_boundary(sgrid,x,y);
[ngrid,nind,list_neumann,list_dirichlet] = ...
    index_boundary(sgrid,grid,ind);
%ngrid',pause,

%% mark noninterior and non-Dirichlet nodes and non-Neumann nodes as strays
stray =  1 - abs(sgrid);
for k=1: size(list_dirichlet,1) stray (list_dirichlet(k,1),list_dirichlet(k,2)) = 0; end;
%sgrid',pause,stray',pause;

%% show grid with boundary type
pcolor(sgrid'),shading interp;colorbar;caxis([-1 1]);title('Grid and boundary conditions');axis equal;pause;
clf;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  transmissibilities and matrix of system
tran1 = perm1/ (dx*dx); tran2 = perm2/ (dy*dy);
velx = velx / (2*dx); vely = vely / (2*dy);

fprintf('creating matrix  ...\n');

tic,
[A] = form_matrix (ngrid.*sgrid,nind,dx,dy,x,y,tran1,tran2,velx,vely);
toc,

fprintf('creating rhs and boundary data ...\n');

%% generate values for rhs and neumann data and dirichlet data 
Rhs_data = generate_rhs (ngrid.*sgrid,nind,x,y);
[Neumann_data, Dirichlet_data] = generate_boundary_data(list_neumann,list_dirichlet, x,y);

elliptic_solve (ngrid,sgrid,nind,dx,dy,x,y,tran1,tran2,velx,vely,...
		 A, Rhs_data, Neumann_data, Dirichlet_data, stray, flagexact);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function elliptic_solve (ngrid,sgrid,nind,dx,dy,x,y,tran1,tran2,velx,vely,...
			  A, rhs, Neumann_data, Dirichlet_data, stray, flagexact)
    
tic,
[f] = form_rhs (ngrid.*sgrid,nind,dx,dy,x,y,tran1,tran2,velx,vely,...
		   rhs, Neumann_data, Dirichlet_data);
toc,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% solve the system 

fprintf('solving linear system ... \n');
tic,u = A \ f; toc,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% postprocess results 
fprintf('visualization ...\n');

[gx,gy] = meshgrid(x,y);
[nx,ny] = size(ngrid);
nx = nx -1; ny = ny -1;

spread_nu = zeros(size(ngrid));
spread_ex = zeros(size(ngrid));

%ngrid',pause

for i=1:nx+1
    for j=1:ny+1
        ind = ngrid(i,j);
        if ind ~= 0
            spread_nu (i,j) = u(abs(ind));
            if flagexact == 1
              spread_ex (i,j) = exfun(x(i),y(j));    
            else
              spread_ex (i,j) = nan;    
            end;
        else 
            if stray (i,j) == 1 
                spread_nu(i,j) = nan;
                spread_ex (i,j) = nan;
            else
                spread_nu(i,j) = Dirichlet_data(i,j);
            end;
            if flagexact == 1
              spread_ex (i,j) = exfun(x(i),y(j));
            end;
        end;
    end
end

%spread_nu',pause,

%%%

if flagexact == 1 

  subplot(1,2,1);surfc ( gx, gy, spread_ex' ); shading interp;title('Exact solution');axis equal;
  subplot(1,2,2);surfc ( gx, gy, spread_nu');shading interp;title('Numerical solution'),axis equal;
  fprintf('dx=%g dy=%g err_inf =%g\n',dx,dy,norm(spread_ex(:)-spread_nu(:),inf));

else

%  surfc ( gx, gy, spread_nu');shading interp;title('Numerical solution'),view(2);axis equal;
  contourf ( gx, gy, spread_nu');title('Numerical solution'),axis equal;
  
   colorbar;
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = rhsfun(x,y)
%v = 2*pi*pi*sin(pi*x)*sin(pi*y);
v = pi*pi*sin(pi*x);
%v = 0;

function v = exfun (x,y)
%v = sin(pi*x)*sin(pi*y);
%v = sin(pi*x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [bgrid,bind,neumann_ij,dirichlet_ij] = index_boundary(sgrid,grid,ind)
[nx,ny] = size (grid);

ninterior = size(ind,2);
agrid = abs(sgrid) - sign(grid);

[agrid,aind] = node_index(agrid,'natural');
bgrid = (agrid + ninterior * sign(agrid)) + grid;

%fprintf('bgrid ...'); bgrid',pause 

nneumann = size(aind,2);
bind = ind;
neumann = [];
neumann_ij = [];

for j=1: nneumann 
    bind ( j + ninterior ).i = aind (j).i; 
    bind ( j + ninterior ).j = aind (j).j;
    neumann = [neumann; bgrid(aind (j).i,aind (j).j)];
    neumann_ij = [neumann_ij; (aind (j).i),(aind (j).j)];
end;
%neumann_ij,

%%%%%%%%%% create a list of all dirichlet nodes: those are next to interior or Neumann nodes

ndirichlet = 0;
dgrid = sparse(nx,ny);
dirichlet_ij = [];

cgrid = sparse(1 - abs(sgrid));
[indi,indj] = find(cgrid);
%fprintf('cgrid...');
%full(cgrid)',pause,

for k=1:size(indi)
    i = indi(k);j=indj(k);
    
    if i>1
        kl = cgrid(i-1,j); 
        if kl==0 && dgrid(i,j) ==0,
            ndirichlet = ndirichlet + 1;
            dgrid(i,j)= ndirichlet; 
            dirichlet_ij = [dirichlet_ij; i,j];
        end;
    end;
    if i<nx
        kr = cgrid(i+1,j); 
        if kr==0 && dgrid(i,j) ==0,
            ndirichlet = ndirichlet + 1;
            dgrid(i,j)= ndirichlet; 
            dirichlet_ij = [dirichlet_ij; i,j];
        end;
    end;
    if j>1
        kb = cgrid(i,j-1); 
        if kb==0 && dgrid(i,j) ==0,
            ndirichlet = ndirichlet + 1;
            dgrid(i,j)= ndirichlet; 
            dirichlet_ij = [dirichlet_ij; i,j];
        end;
    end;
    
    if j<ny 
        kt = cgrid(i,j+1); 
        if kt==0 && dgrid(i,j) ==0,
            ndirichlet = ndirichlet + 1;
            dgrid(i,j)= ndirichlet; 
            dirichlet_ij = [dirichlet_ij; i,j];
        end;
    end;
%       i,j,ndirichlet,pause,
end
%       dirichlet_ij,pause,

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = generate_rhs (grid,ind,x,y)
n = size(ind,2);
f = zeros(n,1);

for k = 1:n
    i = ind(k).i; j=ind(k).j;
    if grid(i,j) > 0  %% internal node
         
        f(k) = rhsfun(x(i),y(j)) ;
        
    elseif grid(i,j) <0 %% Neumann node with higher order approximation
        
     %   f(k) = rhsfun(x(i),y(j)) / 2;   to implement this, lapl must be
     %   changed - computed so that full laplacian in neumann nodes is
     %   computed
    end;        
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Neumann_data,Dirichlet_data] = generate_boundary_data (listn, listd, x, y)
neum_data = zeros(size(listn,1),1);
dir_data = zeros(size(listd,1),1);

for k=1:size(listd,1)
    dir_data(k) = dirfun (x(listd(k,1)),y(listd(k,2)));
end;

for k=1:size(listn,1)
    neum_data(k) = neumfun (x(listn(k,1)),y(listn(k,2)));
end;

%% convert ddat to useful representation
Dirichlet_data = sparse(size(x,2),size(y,2));
for k = 1:size(listd,1) 
    Dirichlet_data(listd(k,1),listd(k,2)) = ...
        dir_data(k,1); 
end;

%% convert ddat to useful representation
Neumann_data = sparse(size(x,2),size(y,2));
for k = 1:size(listn,1) 
    Neumann_data(listn(k,1),listn(k,2)) = ...
        neum_data(k,1) ; 
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = dirfun (xx,yy)
%v = exfun(xx,yy);
%return;
% in (0,1)^2
%if xx == 0 v = 1; else v =0; end;
if xx == -1 v = 1 - yy*yy; else v =0; end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function v = neumfun (xx,yy)
%if yy == 0 v = -pi*sin(pi*xx)*cos(pi*yy); end;
%if yy == 1 v = pi*sin(pi*xx)*cos(pi*yy); end;
v = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function grid = assign_boundary(sgrid,x,y)
%%%%%%% depending on x,y assign type of boundary nodes 
[nx,ny] = size (sgrid);
grid = sgrid;

%%% these are face numbers: use only -1 to denote Neumann edge 
% left: grid(1,2:ny-1) = -1;   
% right: grid(nx,2:ny-1) = -2;
% bottom: grid(2:nx-1,1) = -3;
% top: grid(2:nx-1,ny) = -4;

% for N-S project: 
grid(nx,2:ny-1) = -1;

%% top and bottom and interior inclusions have Neumann condition
%grid (2:nx-1,1) = -1;
%grid (2:nx-1,ny) = -1;

%%% mark interiors of Holes as zeros
for j=2:nx-1 for k=2:ny-1 

      if grid(j,k) == 0 && ...
         (max(grid(j+1,k),0) + max(grid(j-1,k),0) + max(grid(j,k+1),0) + max(grid(j,k-1),0))>0, 
             grid(j,k) = 0; 
      end;
end; end;

