function MDFHelmholtzDbis
%
% 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 APROXIMADA DEL SISTEMA:
% - ALMACENAMIENTO SPARSE DE LA MATRIZ DE COEFICIENTES
% - CONSTRUCCION DEL SEGUNDO MIEMBRO
%   (SIN CONSIDERAR CONDICIONES DE CONTORNO)
% - BLOQUEO DE LOS VALORES FRONTERA DE u
% - 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 Y AUXILIARES
%
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
%                                 % 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; Nxxyym1 = Nxxyy-1; Nxxyyaux = Nx*Nyy;
Nxxyy2m1 = 2*Nxxyy-1; Nxxyy3m2 = 3*Nxxyy-2; Nxxyycasi = Nxxyy3m2+Nxxyyaux;
Nxxyyfin = Nxxyycasi+Nxxyyaux;
hxp2 = hx*hx; coef0 = hx/hy; coef1 = coef0*coef0; coef1D = coef1*D;
coef2 = 2.*D*(1.+coef1) + a*hxp2;
Vmg = 1.e20;
%

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

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

for k=1:Nxxyym1
    kaux1 = Nxxyy+k;
    Ifil(kaux1) = k; Icol(kaux1) = k+1; Acoef(kaux1) = -D;
    kaux2 = Nxxyy2m1+k;
    Ifil(kaux2) = k+1; Icol(kaux2) = k; Acoef(kaux2) = -D;
end;

for k=1:Nxxyyaux
    kaux3 = Nxxyy3m2+k;
    Ifil(kaux3) = k; Icol(kaux3) = Nxx+k; Acoef(kaux3) = -coef1D;
    kaux4 = Nxxyyfin+k;
    Ifil(kaux4) = Nxx+k; Icol(kaux4) = k; Acoef(kaux4) = -coef1D;
end;

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

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

%
% INCORPORACION DE LA CONDICION DE CONTORNO
% (BLOQUEO DE VALORES DE u)
%

y0 = yval(1);
yfin = yval(Nyy);
jfin = Nxxyy - Nxx;
for j=1:Nxx
    jfin = jfin+1;
    Acoef(j) = Vmg; b(j) = Vmg*gDir(xval(j),y0);
    Acoef(jfin) = Vmg; b(jfin) = Vmg*gDir(xval(j),yfin);
end;

x0 = xval(1);
xfin = xval(Nxx);
ind = Nxx;
for i=2:Ny
    yi = yval(i);
    ind = ind + 1;
    Acoef(ind) = Vmg; b(ind) = Vmg*gDir(x0,yi); Iref(ind) = 1;
    ind = ind + Nx;
    Acoef(ind) = Vmg; b(ind) = Vmg*gDir(xfin,yi); Iref(ind) = 1;
end;

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

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

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

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

indu = 0;
for i=1:Nyy
    for j=1:Nxx
        indu = indu + 1;
        uval(i,j) = u(indu);
    end;
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



