Conductividad eléctrica placa 2D(AHM)
De MateWiki
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Conductividad eléctrica de una placa en 2D (Grupo AHM) |
| Asignatura | MNEDP |
| Curso | 2025-26 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
.jpg)
Códigos de Matlab
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 ----
% Cálculo de EOC entre niveles
nlevels = length(hs);
EOC_L2 = nan(nlevels,1);
EOC_H1 = nan(nlevels,1);
for i = 1:nlevels-1
EOC_L2(i) = log( RES(i,1) / RES(i+1,1) ) / log(2);
EOC_H1(i) = log( RES(i,2) / RES(i+1,2) ) / log(2);
end
% Crear tabla con EOC incluidos
TRES = table( ...
hs.', ...
RES(:,1), ...
RES(:,2), ...
RES(:,3), ...
EOC_L2, ...
EOC_H1, ...
'VariableNames', {'h','L2_error','H1_error','condest_A','EOC_L2','EOC_H1'} ...
);
disp(' ');
disp('==================== TABLA FINAL (con EOC) ====================');
disp(TRES);
