function MDFHelmholtzD
%
% METODO DE DIFERENCIAS FINITAS - DIRICHLET-HEMHOLTZ EN UN RECTANGULO
%
% -D(u_{xx} + u_{yy}) + au = f(x,y), (x,y) \in R
% u = g(x,y), (x,y) \in R'
%
% R ES EL RECTANGULO (0,Lx) X (0,Ly)
% R' ES SU FRONTERA
%
% RESOLUCION EXACTA DEL SISTEMA:
% - ALMACENAMIENTO SPARSE DE LA MATRIZ DE COEFICIENTES
% - FACTORIZACION DE CHOLESKY
% - RESOLUCION DE DOS SISTEMAS TRIANGULARES
%
% DATOS:
%
%   D:      COEFICIENTE DEL TERMINO DE ORDEN 2
%   a:      COEFICIENTE DEL TERMINO DE ORDEN 0
%   Lx,Ly:  SUPREMOS EN x,y
%   Nx,Ny:  PASOS EN x,y
%   idisp:  PARAMETRO DE IMPRESION
%   iplot:  PARAMETRO DE REPRESENTACION GRAFICA

%
% DATOS
%
%D = 1.;                           % COEFICIENTE DEL TERMINO DE ORDEN 2
D = 1.e-2;                        % COEFICIENTE DEL TERMINO DE ORDEN 2
a = 1.;                           % COEFICIENTE DEL TERMINO DE ORDEN 0
Lx = 1.; Ly = 1.;                 % SUPREMOS EN x,y
Nx = 10; Ny = 10;                 % NUMERO DE PASOS EN x,y
idisp = 1;                        % PARAMETRO DE IMPRESION
iplot= 1;                         % PARAMETRO DE REPRESENTACION GRAFICA

%
% CALCULOS PRELIMINARES
%
Nxx = Nx+1; Nyy = Ny+1;           % PUNTOS EN x,y
Nxint = Nx-1; Nyint = Ny-1;       % PUNTOS INTERIORES EN x,y
Nxyint = Nxint*Nyint;             % NUMERO TOTAL DE PUNTOS INTERIORES
hx = Lx/Nx; hy = Ly/Ny;           % TAMAÑOS DE LOS PASOS EN x,y
xval = linspace(0.,Lx,Nxx);       % VALORES DE x
yval = linspace(0.,Ly,Nyy);       % VALORES DE y

%
% CALCULOS AUXILIARES
%
Nxyintm1 = Nxyint-1; Nxyint2m1 = 2*Nxyint-1; Nxyint3m2 = 3*Nxyint-2;
Nxm2 = Nx-2; Nxyaux = Nxm2*Nyint;
Nxycasi = Nxyint3m2+Nxyaux; Nxyfin = Nxycasi+Nxyaux;
% Nxxyy = Nxx*Nyy;
hxp2 = hx*hx; coef0 = hx/hy; coef1 = coef0*coef0; coef1D = coef1*D; 
coef2 = 2.*D*(1.+coef1) + a*hxp2;

%
% CONSTRUCCION DE LA MATRIZ DE COEFICIENTES (HUECA, BANDA)
%
Ifil = ones(Nxyfin,1); Icol = ones(Nxyfin,1); Acoef = zeros(Nxyfin,1);

for k=1:Nxyint
    Ifil(k) = k;
    Icol(k) = k;
    Acoef(k) = coef2;
end;

for k=1:Nxyintm1
    kaux1 = Nxyint+k;
    Ifil(kaux1) = k; Icol(kaux1) = k+1; Acoef(kaux1) = -D;
    kaux2 = Nxyint2m1+k;
    Ifil(kaux2) = k+1; Icol(kaux2) = k; Acoef(kaux2) = -D;
    if (rem(k,Nxint) == 0)
        Acoef(kaux1) = 0.;
        Acoef(kaux2) = 0.;
    end;
end;

for k=1:Nxyaux
    kaux3 = Nxyint3m2+k;
    Ifil(kaux3) = k; Icol(kaux3) = Nxint+k; Acoef(kaux3) = -coef1D;
    kaux4 = Nxyfin+k;
    Ifil(kaux4) = Nxint+k; Icol(kaux4) = k; Acoef(kaux4) = -coef1D;
end;

%
% CONSTRUCCION DEL SEGUNDO MIEMBRO
%
b = zeros(Nxyint,1);

indb = 0;
for i=1:Nyint
    yi = yval(i+1);
    for j=1:Nxint
        indb = indb+1;
        b(indb) = hxp2*fsegm(xval(j+1),yi);
    end;
end;

x0 = xval(1); xfin = xval(Nxx);
y0 = yval(1); yfin = yval(Nyy);
jfin = Nxyint-Nxm2;

for j=1:Nxint
    xj = xval(j+1);
    b(j) = b(j) + coef1D*gDir(xj,y0);
    b(jfin) = b(jfin) + coef1D*gDir(xj,yfin);
    jfin = jfin+1;
end;

j1 = 1; jfin = Nxint;
for i=1:Nyint
    yi = yval(i+1);
    b(j1) = b(j1) + D*gDir(x0,yi);
    b(jfin) = b(jfin) + D*gDir(xfin,yi);
    j1 = j1 + Nxint;
    jfin = jfin + Nxint;
end;

if (idisp > 0)
    disp(' EL SEGUNDO MIEMBRO:');
    disp(b);
end;

%
% CONSTRUCCION DE LA MATRIZ EN FORMATO SPARSE, FACTORIZACION DE CHOLESKY
% Y CALCULO DE LA SOLUCION EN NODOS INTERIORES
%
Amat = sparse(Ifil,Icol,Acoef,Nxyint,Nxyint);

if (idisp > 0)
    disp(' LA MATRIZ EN FORMATO SPARSE:');
    disp(Amat);
end;

AChol = chol(Amat);
u = AChol'\b; u = AChol\u;

%
% CONSTRUCCION DE LA SOLUCION EN EL MALLADO COMPLETO
%
uval = zeros(Nyy,Nxx);

y0 = yval(1); yfin = yval(Nyy);
for j=1:Nxx
    xj = xval(j);
    uval(1,j) = gDir(xj,y0);
    uval(Nyy,j) = gDir(xj,yfin);
end;
    
indu0 = 1;
x0 = xval(1); xfin = xval(Nxx);

for i=1:Nyint
    i1 = i+1;
    yi = yval(i1);
    uval(i1,1) = gDir(x0,yi);
    for j=1:Nxint
        j1 = j+1;
        uval(i1,j1) = u(indu0);
        indu0 = indu0 + 1;
    end;
    uval(i1,Nxx) = gDir(xfin,yi);
end;

if (idisp > 0)
    disp(' LOS VALORES DE x, y, u:');
    for i=1:Nyy
        for j=1:Nxx
            disp([xval(j),yval(i),uval(i,j)]);
        end;
    end;
end;

if (iplot > 0)
%
% REPRESENTACION GRAFICA 3D
%
    close all;
    figure;
    surf(xval,yval,uval);
%    mesh(xval,yval,uval);
    xlabel('EJE x'); ylabel('EJE y'); zlabel('SOLUCION');
    title('REPRESENTACION GRAFICA 3D');
%
% CORTES y = Cte.
%
    figure;
    Nycut = round(Nyy/2);
    plot(xval,uval(Nycut,:));
    xlabel('EJE x'); ylabel('SOLUCION');
    title('LA SOLUCION EN FUNCION DE x PARA y = Ly/2');
end;

end

function fxy = fsegm(x,y)
%
%fxy = x*x - y*y;
%fxy = 0.;
xd = x-.5; yd = y-.5;
if (xd*xd + yd*yd <= 1.e-2)
    fxy = 1.;
else
    fxy = 0.;
end;

end

function gxy = gDir(x,y)
%
%gxy = x*x - y*y;
gxy = 0.;
end




