Trabajo realizado por estudiantes
Título	Ecuación del Calor Estacionaria (Grupo DIF)
Asignatura	MNEDP
Curso	2025-26
Autores	Daniel Rodríguez Calderón, Israel López Moreno, Francisco Lavao Amengual.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura
1 Póster del Trabajo

Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre la Ecuación del Calor Estacionaria en formato póster realizado por el grupo DIF.
2 Códigos de Matlab

%% MEF de la Ecuación del Calor con condiciones Dirichlet homogéneas + Neumann homogéneas
clc; clear all; close all

%% Ejemplo de Dominio Cuadrado Unidad

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30; tol = 1e-12;

[x_v, y_v] = meshgrid(linspace(0,1,n_x+1), linspace(0,1,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

triang = delaunay(nodos);
n_nodos = size(nodos,1);
n_triang = size(triang,1);

figure(1)
subplot(1,3,1)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
hold on; xlabel('x'), ylabel('y');
%shading interp; camproj perspective; axis vis3d; view(45,30)
title('\textbf{Dominio cuadrado unidad}', 'Interpreter','latex')
scatter(0.1, 0.5, 'red', '*', LineWidth = 3)
view(2); axis equal tight; 

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1;  % conductividad
c_fun = @(x,y) 0;  % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 350*exp(-10*((x-0.1).^2+(y-0.5).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4'
neumann_sides = ["gamma1","gamma2","gamma3"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
    id = triang(t,:);
    puntos = nodos(id,:);
    x1 = puntos(1,1); y1 = puntos(1,2);
    x2 = puntos(2,1); y2 = puntos(2,2);
    x3 = puntos(3,1); y3 = puntos(3,2);
    
    Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
    B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
    x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
    
    k = k_fun(x_m, y_m);
    Ke = k*Area*(B'*B);
    
    c = c_fun(x_m, y_m);
    Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
    
    rmid = r_fun(x_m, y_m);
    Re = rmid*Area/3*ones(3,1);
    
    for a = 1:3
        for b = 1:3
            I(end+1) = id(a);
            J(end+1) = id(b);
            V(end+1) = Ke(a,b) + Me(a,b);
        end
        R(id(a)) = R(id(a)) + Re(a);
    end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_gamma3 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 | is_gamma4;
neumann_nodes = zeros(length(is_gamma1), 1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
    neumann_side = neumann_sides(i);

    switch lower(neumann_side)
        case 'gamma1'
            dirichlet_nodes(is_gamma1) = false;
            neumann_nodes = neumann_nodes + is_gamma1;
        case 'gamma2'
            dirichlet_nodes(is_gamma2) = false;
            neumann_nodes = neumann_nodes + is_gamma2;
        case 'gamma3'
            dirichlet_nodes(is_gamma3) = false;
            neumann_nodes = neumann_nodes + is_gamma3;
        case 'gamma4'
            dirichlet_nodes(is_gamma4) = false;
            neumann_nodes = neumann_nodes + is_gamma4;
        otherwise
            error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
    end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2); 
subplot(3,3,1)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar; colormap(turbo);
title('\textbf{Solucion FEM (vista superior) sobre el dominio cuadrado unidad}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,2)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio cuadrado unidad}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,3)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y');
title('\textbf{Nodos de contorno sobre el dominio cuadrado unidad}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
    k_fun = @(x,y) k_vals(i);
    
    % Re-ensamblar para cada k
    I = []; J = []; V = []; R = zeros(n_nodos,1);
    for t = 1:n_triang
        id = triang(t,:);
        puntos = nodos(id,:);
        x1 = puntos(1,1); y1 = puntos(1,2);
        x2 = puntos(2,1); y2 = puntos(2,2);
        x3 = puntos(3,1); y3 = puntos(3,2);
        
        Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
        B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
        x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
        
        k = k_fun(x_m, y_m);
        Ke = k*Area*(B'*B);
        c = c_fun(x_m, y_m);
        Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
        rmid = r_fun(x_m, y_m);
        Re = rmid*Area/3*ones(3,1);
        
        for a = 1:3
            for b = 1:3
                I(end+1) = id(a);
                J(end+1) = id(b);
                V(end+1) = Ke(a,b) + Me(a,b);
            end
            R(id(a)) = R(id(a)) + Re(a);
        end
    end
    A_k = sparse(I, J, V, n_nodos, n_nodos);
    b_k = R;
    
    % Impón Dirichlet (mismo conjunto de nodos)
    A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
    A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
    b_k(ind_dir) = 0;
    
    % Resolver y almacenar max
    u_k = A_k \ b_k;
    maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,1)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio cuadrado unidad}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;


%% Ejemplo de Dominio Tipo L

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,2,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = nodos(:,1) > 1 & nodos(:,2) < 1;
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma5 | is_gamma6;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,2)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $L$}', 'Interpreter','latex')
scatter([0.5 1.5], [0.1 1.9], 'red', '*', LineWidth = 3)
view(2); axis equal tight; 

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1;  % conductividad
c_fun = @(x,y) 0;  % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 75*exp(-10*((x-0.5).^2+(y-0.1).^2)) + 75*exp(-10*((x-1.5).^2+(y-1.9).^2)); % Fuente Gaussiana

%%% --- Selección de frontera Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma6"];
g_val = 0; % flujo de Neumann en la pared seleccionada (0 según tu petición)

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
    id = triang(t,:);
    puntos = nodos(id,:);
    x1 = puntos(1,1); y1 = puntos(1,2);
    x2 = puntos(2,1); y2 = puntos(2,2);
    x3 = puntos(3,1); y3 = puntos(3,2);
    
    Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
    B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
    x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
    
    k = k_fun(x_m, y_m);
    Ke = k*Area*(B'*B);
    
    c = c_fun(x_m, y_m);
    Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
    
    rmid = r_fun(x_m, y_m);
    Re = rmid*Area/3*ones(3,1);
    
    for a = 1:3
        for b = 1:3
            I(end+1) = id(a);
            J(end+1) = id(b);
            V(end+1) = Ke(a,b) + Me(a,b);
        end
        R(id(a)) = R(id(a)) + Re(a);
    end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma2 = abs(nodos(:,1)) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma4 = abs(nodos(:,1)- 3) < tol & nodos(:,2) >= 1 & nodos(:,2) <= 2;
is_gamma5 = abs(nodos(:,2)-1) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 1;
is_beta = nodos(:,1) > 1 & nodos(:,1) <= 3 & nodos(:,2) >= 0 & nodos(:,2) < 1;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
                  is_gamma4 | is_gamma5 | is_gamma6 | is_beta;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)
    neumann_side = neumann_sides(i);

    switch lower(neumann_side)
        case 'gamma1'
            dirichlet_nodes(is_gamma1) = false;
            neumann_nodes = neumann_nodes + is_gamma1;
        case 'gamma2'
            dirichlet_nodes(is_gamma2) = false;
            neumann_nodes = neumann_nodes + is_gamma2;
        case 'gamma3'
            dirichlet_nodes(is_gamma3) = false;
            neumann_nodes = neumann_nodes + is_gamma3;
        case 'gamma4'
            dirichlet_nodes(is_gamma4) = false;
            neumann_nodes = neumann_nodes + is_gamma4;
        case 'gamma5'
            dirichlet_nodes(is_gamma5) = false;
            neumann_nodes = neumann_nodes + is_gamma5;
        case 'gamma6'
            dirichlet_nodes(is_gamma6) = false;
            neumann_nodes = neumann_nodes + is_gamma6;
        otherwise
            error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6.')
    end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2); 
subplot(3,3,4)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $L$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,5)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $L$}', 'Interpreter','latex')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,6)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside')
view(2); axis equal tight; xlabel('x'), ylabel('y')
title('\textbf{Nodos de contorno sobre el dominio tipo $L$}', 'Interpreter','latex');

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
    k_fun = @(x,y) k_vals(i);
    
    % Re-ensamblar para cada k
    I = []; J = []; V = []; R = zeros(n_nodos,1);
    for t = 1:n_triang
        id = triang(t,:);
        puntos = nodos(id,:);
        x1 = puntos(1,1); y1 = puntos(1,2);
        x2 = puntos(2,1); y2 = puntos(2,2);
        x3 = puntos(3,1); y3 = puntos(3,2);
        
        Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
        B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
        x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
        
        k = k_fun(x_m, y_m);
        Ke = k*Area*(B'*B);
        c = c_fun(x_m, y_m);
        Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
        rmid = r_fun(x_m, y_m);
        Re = rmid*Area/3*ones(3,1);
        
        for a = 1:3
            for b = 1:3
                I(end+1) = id(a);
                J(end+1) = id(b);
                V(end+1) = Ke(a,b) + Me(a,b);
            end
            R(id(a)) = R(id(a)) + Re(a);
        end
    end
    A_k = sparse(I, J, V, n_nodos, n_nodos);
    b_k = R;
    
    % Impón Dirichlet (mismo conjunto de nodos)
    A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
    A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
    b_k(ind_dir) = 0;
    
    % Resolver y almacenar max
    u_k = A_k \ b_k;
    maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,2)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $L$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;


%% Ejemplo de Dominio Tipo T

%%% --- Parámetros de la malla ---
n_x = 30; n_y = 30;
[x_v, y_v] = meshgrid(linspace(0,3,n_x+1), linspace(0,3,n_y+1));
X = x_v(:); Y = y_v(:);
nodos = [X Y];

fuera_dominio = (nodos(:,1) > 2 & nodos(:,2) < 2) | (nodos(:,1) < 1 & nodos(:,2) < 2);
nodos = nodos(~fuera_dominio,:);

tol = 1e-12;

is_gamma1 = abs(nodos(:,2)) < tol & nodos(:,1) >= 1 & nodos(:,1) <= 2;
is_gamma2 = abs(nodos(:,1)-1) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;
is_gamma3 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 1;
is_gamma4 = abs(nodos(:,1)) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma5 = abs(nodos(:,2)-3) < tol & nodos(:,1) >= 0 & nodos(:,1) <= 3;
is_gamma6 = abs(nodos(:,1)-3) < tol & nodos(:,2) >= 2 & nodos(:,2) <= 3;
is_gamma7 = abs(nodos(:,2)-2) < tol & nodos(:,1) >= 2 & nodos(:,1) <= 3;
is_gamma8 = abs(nodos(:,1)-2) < tol & nodos(:,2) >= 0 & nodos(:,2) <= 2;

triang = delaunay(nodos);
n_nodos = size(nodos,1);

is_edge = is_gamma2 | is_gamma3 | is_gamma7 | is_gamma8;
triang_eliminar = all(is_edge(triang), 2);
triang(triang_eliminar, :) = [];

n_triang = size(triang,1);

figure(1)
subplot(1,3,3)
trimesh(triang, nodos(:,1), nodos(:,2), zeros(size(nodos,1),1))
%shading interp; camproj perspective; axis vis3d; view(45,30)
xlabel('x'), ylabel('y'); hold on;
title('\textbf{Dominio tipo $T$}', 'Interpreter','latex')
scatter(1.5, 0.1, 'red', '*', LineWidth = 3)
view(2); axis equal tight; 

%%% --- Parámetros del problema ---
k_fun = @(x,y) 1; % conductividad
c_fun = @(x,y) 0; % término de reacción
% r_fun = @(x,y) 10*exp(-50*((x-0.5).^2+(y-0.5).^2)); % fuente (como en tus códigos)
r_fun = @(x,y) 125*exp(-10*((x-1.5).^2+(y-0.1).^2)); % Fuente Gaussiana

%%% --- Selección de pared Neumann ---
% Opciones: 'gamma1','gamma2','gamma3','gamma4','gamma5','gamma6','gamma7','gamma8'
neumann_sides = ["gamma1","gamma2","gamma3","gamma5","gamma7","gamma8"];

%%% --- Ensamblaje global ---
I = []; J = []; V = []; R = zeros(n_nodos,1);

for t = 1:n_triang
    id = triang(t,:);
    puntos = nodos(id,:);
    x1 = puntos(1,1); y1 = puntos(1,2);
    x2 = puntos(2,1); y2 = puntos(2,2);
    x3 = puntos(3,1); y3 = puntos(3,2);
    
    Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
    B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
    x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
    
    k = k_fun(x_m, y_m);
    Ke = k*Area*(B'*B);
    
    c = c_fun(x_m, y_m);
    Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
    
    rmid = r_fun(x_m, y_m);
    Re = rmid*Area/3*ones(3,1);
    
    for a = 1:3
        for b = 1:3
            I(end+1) = id(a);
            J(end+1) = id(b);
            V(end+1) = Ke(a,b) + Me(a,b);
        end
        R(id(a)) = R(id(a)) + Re(a);
    end
end

A = sparse(I, J, V, n_nodos, n_nodos);
b = R;

%%% --- Identificación de nodos de contorno y clasificación Dirichlet/Neumann ---
tol = 1e-12;

% Por defecto: Dirichlet en todas las paredes
dirichlet_nodes = is_gamma1 | is_gamma2 | is_gamma3 |...
                  is_gamma4 | is_gamma5 | is_gamma6 |...
                  is_gamma7 | is_gamma8;
neumann_nodes = zeros(length(is_gamma1),1);

% Reemplazar la pared seleccionada por Neumann
for i = 1:numel(neumann_sides)

    neumann_side = neumann_sides(i);

    switch lower(neumann_side)
        case 'gamma1'
            dirichlet_nodes(is_gamma1) = false;
            neumann_nodes = neumann_nodes + is_gamma1;
        case 'gamma2'
            dirichlet_nodes(is_gamma2) = false;
            neumann_nodes = neumann_nodes + is_gamma2;
        case 'gamma3'
            dirichlet_nodes(is_gamma3) = false;
            neumann_nodes = neumann_nodes + is_gamma3;
        case 'gamma4'
            dirichlet_nodes(is_gamma4) = false;
            neumann_nodes = neumann_nodes + is_gamma4;
        case 'gamma5'
            dirichlet_nodes(is_gamma5) = false;
            neumann_nodes = neumann_nodes + is_gamma5;
        case 'gamma6'
            dirichlet_nodes(is_gamma6) = false;
            neumann_nodes = neumann_nodes + is_gamma6;
        case 'gamma7'
            dirichlet_nodes(is_gamma7) = false;
            neumann_nodes = neumann_nodes + is_gamma7;
        case 'gamma8'
            dirichlet_nodes(is_gamma8) = false;
            neumann_nodes = neumann_nodes + is_gamma8;
        otherwise
            error('neumann_side no reconocido. Usa gamma1,gamma2,gamma3,gamma4,gamma5,gamma6,gamma7,gamma8.')
    end
end

ind_dir = find(dirichlet_nodes);
ind_neu = find(neumann_nodes);

%%% --- Aplicación de condiciones Dirichlet homogéneas (en nodos dirichlet_nodes) ---
% Igual que en tu primer código: anular filas/columnas y poner 1 en diagonal
A(ind_dir, :) = 0; A(:, ind_dir) = 0;
A(sub2ind(size(A), ind_dir, ind_dir)) = 1;
b(ind_dir) = 0;

%%% --- Resolución del sistema ---
u = A\b;
u = u(:);

%%% --- Visualización ---
figure(2)
subplot(3,3,7)
trisurf(triang, nodos(:,1), nodos(:,2), u, 'EdgeColor', 'none');
view(2); axis equal tight; colorbar;
title('\textbf{Solucion FEM (vista superior) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

subplot(3,3,8)
trisurf(triang, nodos(:,1), nodos(:,2), u)
shading interp; camproj perspective; axis vis3d; view(125,30)
xlabel('x'), ylabel('y'), zlabel('u')
title('\textbf{Solucion FEM (superficie) sobre el dominio tipo $T$}', 'Interpreter','latex')
xlabel('x'), ylabel('y')

% Mostrar en la malla qué nodos son Dirichlet/Neumann
subplot(3,3,9)
plot(nodos(:,1), nodos(:,2), '.k'); hold on
plot(nodos(ind_dir,1), nodos(ind_dir,2), 'or', 'MarkerFaceColor','r');
plot(nodos(ind_neu,1), nodos(ind_neu,2), 'sb', 'MarkerFaceColor','b');
legend('Interior/otros','Dirichlet','Neumann','Location','bestoutside');
xlabel('x'), ylabel('y'); 
title('\textbf{Nodos de contorno sobre el dominio tipo $T$}', 'Interpreter','latex');
view(2); axis equal tight; xlabel('x'), ylabel('y')

%%% --- Barrido de parámetros k (igual que en tus scripts originales) ---
k_vals = [0.1 1 5 10 15 20];
maxU = zeros(size(k_vals));

for i = 1:length(k_vals)
    k_fun = @(x,y) k_vals(i);
    
    % Re-ensamblar para cada k
    I = []; J = []; V = []; R = zeros(n_nodos,1);
    for t = 1:n_triang
        id = triang(t,:);
        puntos = nodos(id,:);
        x1 = puntos(1,1); y1 = puntos(1,2);
        x2 = puntos(2,1); y2 = puntos(2,2);
        x3 = puntos(3,1); y3 = puntos(3,2);
        
        Area = abs(det([x2-x1 x3-x1; y2-y1 y3-y1]))/2;
        B = [y2-y3 y3-y1 y1-y2; x3-x2 x1-x3 x2-x1]/(2*Area);
        x_m = mean(puntos(:,1)); y_m = mean(puntos(:,2));
        
        k = k_fun(x_m, y_m);
        Ke = k*Area*(B'*B);
        c = c_fun(x_m, y_m);
        Me = (c*Area/12)*[2 1 1; 1 2 1; 1 1 2];
        rmid = r_fun(x_m, y_m);
        Re = rmid*Area/3*ones(3,1);
        
        for a = 1:3
            for b = 1:3
                I(end+1) = id(a);
                J(end+1) = id(b);
                V(end+1) = Ke(a,b) + Me(a,b);
            end
            R(id(a)) = R(id(a)) + Re(a);
        end
    end
    A_k = sparse(I, J, V, n_nodos, n_nodos);
    b_k = R;
    
    % Impón Dirichlet (mismo conjunto de nodos)
    A_k(ind_dir, :) = 0; A_k(:, ind_dir) = 0;
    A_k(sub2ind(size(A_k), ind_dir, ind_dir)) = 1;
    b_k(ind_dir) = 0;
    
    % Resolver y almacenar max
    u_k = A_k \ b_k;
    maxU(i) = max(u_k);
end

figure(3)
subplot(1,3,3)
plot(k_vals, maxU, '-o', 'LineWidth', 1.5);
xlabel('$k$', 'Interpreter','latex'); ylabel('\textbf{max \{u\}}', 'Interpreter','latex');
title('\textbf{Efecto de $k$ en la solucion sobre el dominio tipo $T$}', 'Interpreter','latex');
grid on; hold on; view(2); axis square;
Ecuación del Calor Estacionaria (Grupo DIF)

1 Póster del Trabajo

2 Códigos de Matlab

Menú de navegación

Herramientas personales

Espacios de nombres

Variantes

Vistas

Más

Buscar

Navegación

Herramientas
