clc 
clear all
close all

%% CONFIGURACIÓN DEL CASO 
%% ========== Caso 1 ==========
alpha = 0;
u_exact_fun = @(x,y) sin(pi*x/3).*sin(pi*y/3); 

sigma_1 = @(x,y) 1 + 0.5*sin(pi*x/3).*cos(pi*y/3);

du_dx_fun = @(x,y) (pi/3) * cos(pi*x/3) .* sin(pi*y/3);
du_dy_fun = @(x,y) (pi/3) * sin(pi*x/3) .* cos(pi*y/3); 
d2u_dx2_fun = @(x,y) -(pi/3)^2 * sin(pi*x/3) .* sin(pi*y/3);
d2u_dy2_fun = @(x,y) -(pi/3)^2 * sin(pi*x/3) .* sin(pi*y/3);

dsigma_dx_fun = @(x,y) 0.5*(pi/3)*cos(pi*x/3).*cos(pi*y/3);
dsigma_dy_fun = @(x,y) -0.5*(pi/3)*sin(pi*x/3).*sin(pi*y/3);

f_fun = @(x,y) - ( ...
    sigma_1(x,y).*(d2u_dx2_fun(x,y) + d2u_dy2_fun(x,y)) + ...
    dsigma_dx_fun(x,y).*du_dx_fun(x,y) + ...
    dsigma_dy_fun(x,y).*du_dy_fun(x,y) );

g_fun = @(x,y) - sigma_1(x,0).*du_dy_fun(x,0);  % ← ahora siempre g(x,y)

u0 = @(x,y) 0;

%% ========== Caso 2 ==========
alpha = 0;

% --- Solución exacta diseñada para ser 0 en x=0,x=3,y=3 ---
% u(x,y) = sin(pi*x/3) * y*(3-y)
u_exact_fun = @(x,y) sin(pi*x/3) .* ( y .* (3 - y) );

% --- Derivadas de u_exact ---
du_dx_fun   = @(x,y) (pi/3) * cos(pi*x/3) .* ( y .* (3 - y) );
du_dy_fun   = @(x,y) sin(pi*x/3) .* ( 3 - 2*y );
d2u_dx2_fun = @(x,y) -(pi/3)^2 * sin(pi*x/3) .* ( y .* (3 - y) );
d2u_dy2_fun = @(x,y) sin(pi*x/3) .* ( -2 );

% --- Usamos tu sigma (variable, positiva y suave) ---
sigma_1 = @(x,y) 1 + 0.5*sin(pi*x/3).*cos(pi*y/3);

% --- Derivadas de sigma (necesarias para f) ---
dsigma_dx_fun = @(x,y) 0.5*(pi/3)*cos(pi*x/3).*cos(pi*y/3);
dsigma_dy_fun = @(x,y) -0.5*(pi/3)*sin(pi*x/3).*sin(pi*y/3);

% --- Segundo miembro f(x,y) consistente con u_exact ---
f_fun = @(x,y) - ( ...
    sigma_1(x,y) .* ( d2u_dx2_fun(x,y) + d2u_dy2_fun(x,y) ) + ...
    dsigma_dx_fun(x,y) .* du_dx_fun(x,y) + ...
    dsigma_dy_fun(x,y) .* du_dy_fun(x,y) ...
);

% --- Condición de Neumann NO homogénea en y=0  ---
g_fun = @(x,y) - sigma_1(x,0) .* du_dy_fun(x,0);

% --- Condición de Dirichlet homogénea  ---
u0 = @(x,y) 0;

%% Casos extra (Generalización del problema, Diferentes a nuestro problema original)
%% ========== Caso 3 ==========
alpha = 0;
u_exact_fun = @(x,y) sin(pi*x/3).*y;

sigma_1 = @(x,y) 1 + 0.5*sin(pi*x/3).*cos(pi*y/3);

du_dx_fun = @(x,y) (pi/3)*cos(pi*x/3).*y;
du_dy_fun = @(x,y) sin(pi*x/3);
d2u_dx2_fun = @(x,y) -(pi/3)^2*sin(pi*x/3).*y;
d2u_dy2_fun = @(x,y) 0;

dsigma_dx_fun = @(x,y) 0.5*(pi/3)*cos(pi*x/3).*cos(pi*y/3);
dsigma_dy_fun = @(x,y) -0.5*(pi/3)*sin(pi*x/3).*sin(pi*y/3);

f_fun = @(x,y) - ( ...
    sigma_1(x,y).*d2u_dx2_fun(x,y) + ...
    dsigma_dx_fun(x,y).*du_dx_fun(x,y) + ...
    dsigma_dy_fun(x,y).*du_dy_fun(x,y) );

g_fun = @(x,y) - sigma_1(x,0).*du_dy_fun(x,0);
u0 = @(x,y) u_exact_fun(x,y);

%% ========== Caso 4  ==========
alpha = 0;
u_exact_fun = @(x,y) exp(-(x/3).^2).*cos(pi*y/3);

sigma_1 = @(x,y) 1 + x.*y + sin(pi*x).*cos(pi*y);

du_dx_fun = @(x,y) (-2*x/9).*exp(-(x/3).^2).*cos(pi*y/3);
du_dy_fun = @(x,y) (-pi/3).*exp(-(x/3).^2).*sin(pi*y/3);
d2u_dx2_fun = @(x,y) (2/9*((4*x.^2)/9 - 1)).*exp(-(x/3).^2).*cos(pi*y/3);
d2u_dy2_fun = @(x,y) -(pi/3)^2.*exp(-(x/3).^2).*cos(pi*y/3);

dsigma_dx_fun = @(x,y) y + pi*cos(pi*x).*cos(pi*y);
dsigma_dy_fun = @(x,y) x - pi*sin(pi*x).*sin(pi*y);

f_fun = @(x,y) - ( ...
    sigma_1(x,y).*(d2u_dx2_fun(x,y)+d2u_dy2_fun(x,y)) + ...
    dsigma_dx_fun(x,y).*du_dx_fun(x,y) + ...
    dsigma_dy_fun(x,y).*du_dy_fun(x,y) );
g_fun = @(x,y) 0;
u0 = @(x,y) u_exact_fun(x,y);



%% Estructura Común
%% --- Configuración común ---
u_ex  = u_exact_fun;
ux_ex = du_dx_fun;
uy_ex = du_dy_fun;

hs = [0.5, 0.25, 0.125, 0.0625];
RES = zeros(4,3);
tol = 1e-10;

% ===== BUCLE EN MALLAS =====
for k = 1:4
    fprintf('--- Procesando Malla %d ---\n', k);

    Vertices  = load(sprintf('verticesPoligono4_malla%d.txt', k));
    Elementos = load(sprintf('ElementosPoligono4_malla%d.txt', k));
    Dirichlet = load(sprintf('DirichletPoligono4_malla%d.txt', k));
    
    Dir = Dirichlet(:,1);
    X = Vertices(:,2);
    Y = Vertices(:,3);
    TRI = Elementos(:,2:4);
    M = length(X); Nel = size(TRI,1);

    Neumann = find(abs(Y-0)<tol);
    Neumann = setdiff(Neumann, Dir);

    % ---- Ensamblado ----
    A = sparse(M,M);
    b = zeros(M,1);

    Lxx = [1/2 -1/2 0; -1/2 1/2 0; 0 0 0];
    Lyy = [1/2 0 -1/2; 0 0 0; -1/2 0 1/2];
    Lxy = [1/2 0 -1/2; -1/2 0 1/2; 0 0 0];
    l0  = [1/6;1/6;1/6];

    for e = 1:Nel
        nodes = TRI(e,:);
        xK = X(nodes); yK = Y(nodes);
        BK = [xK(2)-xK(1), xK(3)-xK(1); yK(2)-yK(1), yK(3)-yK(1)];
        absdetB = abs(det(BK));
        CK = inv(BK)*inv(BK)';

        sigmaK = mean(sigma_1(xK,yK));

        LK = absdetB * sigmaK * ( CK(1,1)*Lxx + CK(1,2)*(Lxy+Lxy') + CK(2,2)*Lyy );
        fK = mean( f_fun(xK,yK) );
        lK = fK * absdetB * l0;

        A(nodes,nodes)=A(nodes,nodes)+LK;
        b(nodes)=b(nodes)+lK;
    end

    % ---- Construcción Neu ----
    edges = [TRI(:,[1 2]); TRI(:,[2 3]); TRI(:,[3 1])];
    edges_sorted = sort(edges,2);
    [u_edges,~,ic] = unique(edges_sorted,'rows');
    counts = accumarray(ic,1);
    boundary_edges = u_edges(counts==1,:);

    Neu = [];
    for e2 = 1:size(boundary_edges,1)
        n1 = boundary_edges(e2,1);
        n2 = boundary_edges(e2,2);
        if abs(Y(n1))<tol && abs(Y(n2))<tol && ~ismember(n1,Dir) && ~ismember(n2,Dir)
            Neu = [Neu; n1 n2];
        end
    end

    % ---- Integral Parte Neumann ----
    for e2 = 1:size(Neu,1)
        n1 = Neu(e2,1);
        n2 = Neu(e2,2);

        Ledge = sqrt( (X(n2)-X(n1))^2 + (Y(n2)-Y(n1))^2 );

        g1 = g_fun(X(n1),Y(n1));
        g2 = g_fun(X(n2),Y(n2));

        b_edge = (Ledge/6)*[2*g1+g2; g1+2*g2];   % ← fórmula CORRECTA

        b([n1 n2]) = b([n1 n2]) + b_edge;
    end

    % ---- Imposición Condición Dirichlet ----
    u = zeros(M,1);
    u(Dir) = arrayfun(@(i)u0(X(i),Y(i)), Dir);

    ind = setdiff(1:M, Dir);
    b(ind) = b(ind) - A(ind,Dir)*u(Dir);
    u(ind) = A(ind,ind)\b(ind);

    % ---- Cálculo de errores ----
    errL2_sq = 0; errH1_sq = 0;

    for e = 1:Nel
        nodes = TRI(e,:);
        xloc = X(nodes); yloc = Y(nodes);

        BK = [xloc(2)-xloc(1), xloc(3)-xloc(1); ...
              yloc(2)-yloc(1), yloc(3)-yloc(1)];
        detBK = abs(det(BK));
        BinvT = inv(BK)';

        ue = u_ex(xloc,yloc);
        uhe = u(nodes);

        errL2_sq = errL2_sq + detBK/3 * sum( (ue(:)-uhe(:)).^2 );

        grad_uh = BinvT * [uhe(2)-uhe(1); uhe(3)-uhe(1)];

        g1 = [ux_ex(xloc(1),yloc(1)); uy_ex(xloc(1),yloc(1))];
        g2 = [ux_ex(xloc(2),yloc(2)); uy_ex(xloc(2),yloc(2))];
        g3 = [ux_ex(xloc(3),yloc(3)); uy_ex(xloc(3),yloc(3))];
        grad_u_avg = (g1+g2+g3)/3;

        errH1_sq = errH1_sq + detBK * sum( (grad_u_avg-grad_uh).^2 );
    end

    RES(k,1) = sqrt(errL2_sq);
    RES(k,2) = sqrt(errH1_sq);
    RES(k,3) = condest(A(ind,ind));

    % --- Visualización ---
        
        % ---- Recalcular nodos Neumann para representación ----
        Neumann = find(abs(Y - 0) < tol);
        Neumann = setdiff(Neumann, Dir);

        % ---- Calcular solución exacta para representar ----
        u_exact = u_ex(X, Y);

        % ---- Representación de la malla y los bordes ----
        figure(1); 
        sgtitle('Evolución de Mallados y Condiciones de Contorno');
        % Seleccionamos el subplot actual para la malla
        subplot(2, 2, k); 
        
        % Trazado de la malla base (en z=0)
        trisurf(TRI, X, Y, zeros(size(X)), ...
            'FaceColor','cyan', 'EdgeColor','k', 'FaceAlpha',0.7);
        hold on;
        
        % Dirichlet en rojo (z=0)
        plot3(X(Dir), Y(Dir), zeros(size(Dir)), ...
            'r*', 'MarkerSize',8, 'LineWidth',1.5);
            
        % Neumann en verde (z=0)
        plot3(X(Neumann), Y(Neumann), zeros(size(Neumann)), ...
            'g*', 'MarkerSize',8, 'LineWidth',1.5);
            
        title(sprintf('Mallado %d (h=%.4f)', k, hs(k)));
        xlabel('x'); ylabel('y'); zlabel('z=0');
        axis equal; grid on; view(3);
        
        % Añadimos una unica leyenda para claridad
        if k == 1
            legend('Malla','Dirichlet','Neumann', 'Location','best');
        end
        hold off;


        % -------------------------------------------------------------------------
        % VISUALIZACIÓN: SUBPLOTS (2x2)
        % -------------------------------------------------------------------------
        
         
        figure(2); 
        sgtitle('Comparación Solución Exacta vs. Solución FEM por Mallado');
        
        % Dentro del bucle, seleccionamos el subplot actual
        subplot(2, 2, k); 
        
        hold on; grid on; view(3);
        
        % Sol exacta (Superficie oscura)
        trisurf(TRI, X, Y, u_exact, ...
            'FaceColor',[0.2 0.2 0.2], ...
            'EdgeColor','none', 'FaceAlpha',0.8);
        
        % FEM encima (Mallado claro)
        trisurf(TRI, X, Y, u, ...
            'FaceColor','none', ...
            'EdgeColor',[0.5 1.0 1.0], ...
            'LineWidth',1.5);
        
        title(sprintf('Malla %d (h=%.4f)', k, hs(k))); % Usamos hs(k) para mostrar el tamaño de malla
        xlabel('x'); ylabel('y'); zlabel('u');
        axis equal; 
        colorbar;
        legend('Sol. Exacta','Sol. FEM','Location','best');
        hold off;
end
% ---- Tabla final de errores ----
TRES = table(hs.', RES(:,1), RES(:,2), RES(:,3), ...
    'VariableNames', {'h','L2_error','H1_error','condest_A'});
disp(TRES)

pL2 = polyfit(log(hs(:)), log(RES(:,1)), 1);
pH1 = polyfit(log(hs(:)), log(RES(:,2)), 1);

fprintf('Orden L2 ≈ %.3f\n', pL2(1));
fprintf('Orden H1 ≈ %.3f\n', pH1(1));
fprintf('Escalado esperado cond(A) ~ h^{-2}\n');