Diferencia entre revisiones de «Placa Plana (J52)»

clc; clear; close all;
%Datos
h = 0.1;
u = 0:h:4;
% Para la coordenada vertical, la altura máxima de la placa es 2 por lo que 
% para mantener el paso h=0.1, necesitamos dividir esa altura en aprox 20 partes.
% y definimos un parámetro 's' que vaya de 0 a 1.
num_puntos_y=2/h; 
s=linspace(0,1,num_puntos_y + 1);
%rejilla de parámetros
[U,S]=meshgrid(u,s);

% fronteras (borde inferior y superior)
f_u=U./8;
g_u=2-U./8;

%Calcular las coordenadas físicas X e Y
% X es U
X=U;
% Y es una interpolación entre el borde inferior y superior
Y=f_u+S.*(g_u-f_u);

%Dibujar el mallado
figure(1);
clf; 
hold on;

%hacemos la rejilla y la configuramos
mesh(X,Y,zeros(size(X)),'EdgeColor',[0, 0.7, 0.7],'FaceColor','none');

%Dibujamos un contorno parecido a la imagen dada
plot(X(1,:),Y(1,:),'k','LineWidth',2)% Borde inferior
plot(X(end,:),Y(end,:),'k','LineWidth',2)% Borde superior
plot(X(:,1),Y(:,1),'k','LineWidth',2)% Borde izquierdo (pared)
plot(X(:,end),Y(:,end),'k','LineWidth',2)% Borde derecho

% Configuramos los ejes y el visionado
axis([-1 5 -1 3]);
xlabel('x');
ylabel('y');
title('Mallado de la placa plana');
view(2);
box on;
grid off;
axis equal

clc; clear; close all;
% Definición del Mallado y Dominio
h=0.1;  %Muestreo                  
x=0:h:4;                  
y=0:h:2;           
[X,Y] = meshgrid(x, y);

% Funciones que definen la placa 
f=@(x)x./8;  % Límite inferior
g=@(x)2-x./8;  % Límite superior

% Máscara para puntos dentro de la placa
EnPlaca=Y>=f(X)&Y< g(X);

% Definición de la Temperatura y el Gradiente
% Función T(x,y) 
T=(1+(Y-1).^2).*(4-X);

% Aplicar máscara
T_plot=T;
T_plot(~EnPlaca)=NaN;

% Cálculo numérico del gradiente para los vectores (o usar fórmula exacta)
% Usamos la fórmula exacta derivada arriba para mayor precisión
U=-(1+(Y-1).^2);          % Componente x del gradiente
V=2.*(Y-1).*(4-X);    % Componente y del gradiente

% Solo mostramos vectores dentro de la placa
U(~EnPlaca)=NaN;
V(~EnPlaca)=NaN;

% Generación de la Gráfica
figure(1);
clf; % Limpiar figura
hold on;
axis equal; 
axis([-0.5 4.5 -0.5 2.5]);
grid on;
title({'Isotermas y Gradiente \nabla T', 'T(x,y) = (1+(y-1)^2)(4-x)'});
xlabel('x'); ylabel('y');

% Mapa de calor de fondo
[C_fill, h_fill] = contourf(X, Y, T_plot, 20, 'LineStyle', 'none');
colormap(jet); % Mapa de colores tipo jet
cb=colorbar;
ylabel(cb,'Temperatura T');

% Curvas de nivel
[C, h_cont]=contour(X,Y,T_plot,10,'k','LineWidth',1.5);
clabel(C,h_cont,'FontSize',9,'Color','k','FontWeight','bold');

% Campo Vectorial (Gradiente)
% Submuestreamos cada 3 puntos para que no se vea una mancha de flechas
step=3; 
quiver(X(1:step:end,1:step:end),Y(1:step:end, 1:step:end), ...
U(1:step:end,1:step:end),V(1:step:end,1:step:end), ...
1.2,'k','LineWidth',1); % '1.2' es el factor de escala de las flechas

% Dibujar bordes de la placa
plot(x,f(x),'k-','LineWidth',2); % Borde inferior
plot(x,g(x),'k-','LineWidth',2); % Borde superior
plot([0 0],[f(0) g(0)],'k-','LineWidth',2); % Pared izquierda
plot([4 4],[f(4) g(4)],'k-','LineWidth',2); % Pared derecha

% Señalar Máximo
plot(0,0,'rp','MarkerSize',12,'MarkerFaceColor','r');
plot(0,2,'rp','MarkerSize',12,'MarkerFaceColor','r');
text(0.1,0,'Max T','Color','white','FontWeight','bold');

hold off;

clc; clear; close all;

% Definición del Escenario
h=0.1; 
[X,Y]=meshgrid(0:h:4,0:h:2);
f=X./8; 
g=2-X./8;
EnPlaca=Y>=f&Y<=g; % Máscara

% Cálculos y Gradiente de T
Grad_x=-(1+(Y-1).^2); 
Grad_y=2.*(Y-1).*(4-X);

% Flujo de Calor Q (Ley de Fourier)
k=1; % Constante de conductividad
Qx=-k.*Grad_x;
Qy=-k.*Grad_y;

% Borrar vectores fuera de la placa
Qx(~EnPlaca)=NaN;
Qy(~EnPlaca)=NaN;

% Visualización
figure('Name','Flujo de Calor (Ley de Fourier)');
hold on;
axis equal;axis([-0.5 4.5 -0.5 2.5]);grid on;
title({'Flujo de Calor', 'El calor viaja desde lo caliente a lo frío'});
xlabel('X'); ylabel('Y');

% Temperatura de fondo
T=(1+(Y-1).^2).*(4-X);
T(~EnPlaca)=NaN;
contourf(X,Y,T,15,'LineStyle','none','FaceAlpha',0.5);
colormap(jet);c=colorbar;ylabel(c,'Temperatura');

% Dibujar el Campo Vectorial Q
% Usamos 'step' para no pintar todas las flechas y que se vea limpio
step=3; 
idx=1:step:size(X, 1);
idy=1:step:size(X, 2);

quiver(X(idx,idy),Y(idx,idy),Qx(idx,idy),Qy(idx,idy),1.5,'k','LineWidth',0.75); 
 
% Bordes
x_b=linspace(0,4, 200);
plot(x_b,f(x_b),'k','LineWidth',2);
plot(x_b,g(x_b),'k','LineWidth',2);
plot([0 0],[0 2],'k','LineWidth',2);
plot([4 4],[0.5 1.5],'k','LineWidth',2);

% Origen del flujo
plot(0,0,'ko','MarkerFaceColor','w');
plot(0,2,'ko','MarkerFaceColor','w'); 

hold off;

clc; clear; close all;

h=0.05;
[X,Y]=meshgrid(0:h:4,0:h:2);
f=X./8;g=2-X./8;
EnPlaca=Y>=f&Y<=g;

% Cálculos de Gradiente y Norma
Tx=-(1+(Y-1).^2); 
Ty=2.*(Y-1).*(4-X);

% Norma del Gradiente
NormaGrad=sqrt(Tx.^2+Ty.^2);
NormaGrad(~EnPlaca)=0; 

% Buscar el Máximo
max_val=max(NormaGrad(:));
[filas,cols]=find(NormaGrad==max_val);

% Coordenadas de los puntos máximos
x_max=X(filas,cols);
y_max=Y(filas,cols);

% Componentes del vector dirección en esos puntos
u_dir=Tx(filas,cols);
v_dir=Ty(filas,cols);

% Visualización
figure('Name','Máxima Variación de T');
hold on;axis equal;axis([-0.5 4.5 -0.5 2.5]);grid on;
title({'Puntos de Máxima Variación de Temperatura', '(Color = Magnitud del Gradiente T)'});
xlabel('x'); ylabel('y');

% Mapa de calor
NormaGrad(~EnPlaca)=NaN;
contourf(X,Y,NormaGrad,20,'LineStyle','none');
colormap(jet);colorbar;

% Contorno
plot([0 4 4 0 0],[0 0.5 1.5 2 0],'k','LineWidth',1.5);

% Puntos rojos y dirección
plot(x_max, y_max,'ro','MarkerSize',7.5,'MarkerFaceColor','r');

% Flechas de dirección
quiver(x_max,y_max,u_dir,v_dir,0.5,'k','LineWidth',2,'MaxHeadSize',0.5);

% Resultados
fprintf('La variación máxima es %.2f\n', max_val);
fprintf('Ocurre en los puntos:\n');
for i=1:length(x_max)
    fprintf('P(%d): x=%.2f, y=%.2f. Dirección del vector: [%.2f, %.2f]\n',i,x_max(i),y_max(i),u_dir(i),v_dir(i));
end

clc; clear; close all;

% Definición del Dominio y Mallado
h = 0.1; % Muestreo
x = 0:h:4; 
y= 0:h:2;
[X,Y]=meshgrid(x,y);

% Definición de la Placa
f=@(x)x./8;     % Límite inferior
g=@(x)2-x./8;   % Límite superior

% Máscara lógica
EnPlaca=Y>=f(X)&Y<=g(X); 

% Cálculo del Campo de Desplazamientos u
U=(1/20).*X.*Y;   % Componente horizontal
V=-(1/20).*X.^2;  % Componente vertical

% Filtrar
U(~EnPlaca)=NaN;
V(~EnPlaca)=NaN;

% Gráfica y ajuste
figure('Name', 'Campo de Desplazamientos'); 
set(gcf,'Position',[100 100 800 500]); 
hold on;
axis equal;
axis([-0.5 4.5 -1 2.5]);
grid on;
title('Campo de Desplazamientos');
xlabel('x'); ylabel('y');

% Dibujar el campo vectorial u 
quiver(X,Y,U,V,0,'r', 'LineWidth',1,'MaxHeadSize', 0.5);

% Dibujar contorno de la placa como referencia
x_borde = linspace(0,4,200);
plot(x_borde,f(x_borde),'k','LineWidth',1);
plot(x_borde,g(x_borde),'k','LineWidth',1);
plot([0 0],[f(0) g(0)],'k','LineWidth',1); % Pared izquierda
plot([4 4],[f(4) g(4)],'k','LineWidth',1); % Pared derecha

% Puntos fijos
% Buscamos puntos donde el desplazamiento sea (0,0)
% Si x=0, entonces U=0 y V=0.Por lo que toda la pared izquierda (x=0) está fija.

plot(zeros(size(y)),y,'b.','MarkerSize',15); % Puntos azules
text(0.1,1,'\leftarrow Puntos Fijos (x=0)','Color','b','FontWeight','bold');

hold off;

clc; clear; close all;

% Generación de la Malla
h=0.2; 
u=0:h:4; 
s=linspace(0,1,11);
[U,S]=meshgrid(u,s);

% Coordenadas Originales
f_x=U./8;
g_x=2-U./8;
X0=U;
Y0=f_x+S.*(g_x-f_x);
Z0=zeros(size(X0));

% Cálculo de la Deformación y Campo de desplazamientos u(x,y)
Ux=(X0.*Y0)./20;
Uy=-(X0.^2)./20;

% Coordenadas Deformadas
Xd=X0+Ux;
Yd=Y0+Uy;

% Visualización
figure('Name','Comparación Antes/Después','Color','w');

% ANTES
subplot(1,2,1);
mesh(X0,Y0,Z0,'EdgeColor',[0 0.7 0.9],'FaceColor','none');
hold on;
plot([0 0],[0 2],'k','LineWidth',3); % Pared fija
title('Antes del desplazamiento');
xlabel('X'); ylabel('Y');
axis equal; axis([-0.5 4.5 -1.5 2.5]); grid on;
view(2);

% DESPUÉS
subplot(1,2,2);
mesh(Xd,Yd,Z0,'EdgeColor',[0 0.7 0.9],'FaceColor','none');
hold on;
plot([0 0], [0 2],'k','LineWidth',3);
title('Después del desplazamiento');
xlabel('X'); ylabel('Y');
axis equal; axis([-0.5 4.5 -1.5 2.5]); grid on;
view(2);

clc; clear; close all;
% Definir Malla y Placa
h=0.05;
[X,Y]=meshgrid(0:h:4, 0:h:2);
f=X./8;g=2-X./8; % Límites
EnPlaca=Y>=f&Y<=g; % Máscara lógica

% Cálculo de la Divergencia (div=y/20)
Div=Y./20; 
Div(~EnPlaca)=NaN;

% Visualización
figure('Name','Divergencia');clf;hold on;
axis equal; axis([-0.5 4.5 -0.5 2.5]);box on;
title('Divergencia: Cambio de Volumen');
xlabel('X'); ylabel('Y');
[C, h_cont]=contourf(X, Y, Div, 20,'LineStyle','none');
colormap(jet);
c=colorbar;ylabel(c,'Valor de Divergencia');

% Dibujar bordes negros
x_b = linspace(0, 4, 200);
plot(x_b, x_b./8,'k','LineWidth',1);
plot(x_b, 2 - x_b./8,'k','LineWidth',1);
plot([0 0],[0 2],'k','LineWidth',1);
plot([4 4],[0.5 1.5],'k','LineWidth',1);

% Señalar puntos clave
% Máxima: Ocurre en y=2 y x=0 
plot(0,2,'kp','MarkerSize',12,'MarkerFaceColor','r');
text(0.1, 2.1,'Max Expansión','FontWeight','bold');

% Nula: Ocurre en y=0 y x=0
plot(0,0,'kp','MarkerSize',12,'MarkerFaceColor','b');
text(0.1,-0.1,'Divergencia Nula','FontWeight','bold');
hold off;

clc; clear; close all;

% Definición del Mallado
h=0.05; 
[X,Y] =meshgrid(0:h:4, 0:h:2);
Condition =(Y >= X./8) & (Y <= 2 - X./8);

% Cálculo del Rotacional
% Fórmula derivada anteriormente: Rot_z = -3x/20
% Tomamos el valor absoluto (módulo) para dibujar la magnitud
Rot_Modulo = abs(-(3 .* X) ./ 20);

% Filtrado 
Rot_Modulo(~Condition) = NaN; %(Borrar lo de fuera)

% Visualización 3D y 2D)
figure('Name', 'Rotacional', 'Position', [100, 50, 600, 800]);

% VISTA 3D
subplot(2, 1, 1); 
surf(X, Y, Rot_Modulo); 
shading interp; 
colormap(jet); 
colorbar;
axis tight; 
grid on;
view(-30, 45);
title('Vista 3D');
zlabel('Magnitud');
xlabel('x'); ylabel('y');
text(0, 1, 0.6, 'Sube como una rampa', 'Color', 'k');

% VISTA 2D
subplot(2, 1, 2);
contourf(X, Y, Rot_Modulo, 20, 'LineStyle', 'none');
hold on;
plot([0 4 4 0 0], [0 0.5 1.5 2 0], 'k', 'LineWidth', 2); % Bordes
colormap(jet); 
colorbar;
axis equal; axis([-0.5 4.5 -0.5 2.5]);
title('Vista 2D');
xlabel('x'); ylabel('y');

% Señalar el máximo en 2D
plot([4 4], [0.5 1.5], 'r', 'LineWidth', 3);
text(2.5, 2, 'Mayor giro en x=4', 'FontWeight', 'bold');

hold off;

clc; clear; close all;

% Geometría de la placa
h = 0.1;                        
x = 0:h:4;
num_puntos_y = round(2/h);
s = linspace(0,1,num_puntos_y+1);% parámetro vertical [0,1]

[U,S] = meshgrid(x,s);
f = U./8;                        
g = 2 - U./8;                   

X = U;                           
Y = f + S.*(g - f);           

% Campo de desplazamientos aplicado
% Desplazamiento dado en polares: u~(rho,theta) = -(1/20) rho^2 cos(theta) e_theta
% Convertido a cartesianas: u = [xy/20,  -x^2/20]
u = (X.*Y)/20;                   % componente en x
v = -(X.^2)/20;                  % componente en y
% w = 0 (desplazamientos planos)

% Derivadas analíticas
ux_x = Y/20;                    
ux_y = X/20;                   
vy_x = -(X/10);                 
vy_y = zeros(size(X));           

div_u = ux_x + vy_y; % divergencia = y/20

% Tensor de deformaciones
eps_xx = ux_x;                  
eps_yy = vy_y;                  
eps_xy = 0.5*(ux_y + vy_x);      

% Tensiones 
lambda = 1; mu = 1;
sigma_xx = lambda*div_u + 2*mu*eps_xx;   % = 3y/20
sigma_yy = lambda*div_u + 2*mu*eps_yy;   % = y/20
sigma_zz = lambda*div_u;                 % = y/20 (eps_zz=0)
sigma_xy = 2*mu*eps_xy;                  % = -x/20
sigma_yx = sigma_xy;                     % simétrico
sigma_xz = zeros(size(X));               % no hay variación en z
sigma_zx = zeros(size(X));               % idem

% Tensión tangencial respecto al plano ortogonal a i 
% tau_i = |sigma*i - (i·sigma·i) i| = sqrt( sigma_yx^2 + sigma_zx^2 )
tau_i = abs(sigma_yx);                   % = | -x/20 |

% Gráficos

cm = turbo; % o 'jet'
lw = 1.2;

% Mallado y contornos de la placa
figure('Name','Mallado de la placa','Color','w');
mesh(X,Y,zeros(size(X)),'EdgeColor',[0,0.6,0.7],'FaceColor','none'); hold on;
plot(X(1,:),Y(1,:),'k','LineWidth',2);       % borde inferior
plot(X(end,:),Y(end,:),'k','LineWidth',2);    % borde superior
plot(X(:,1),Y(:,1),'k','LineWidth',2);        % borde izquierdo (pared)
plot(X(:,end),Y(:,end),'k','LineWidth',2);    % borde derecho
axis([-0.5 4.5 -0.5 2.5]); axis equal; grid on; box on;
xlabel('x'); ylabel('y'); title('Mallado y contorno de la placa');

% Tensiones normales
figure('Name','Tensiones normales','Color','w');
tiledlayout(1,3,'Padding','compact','TileSpacing','compact');

nexttile;
surf(X,Y,sigma_xx); shading interp; colormap(cm); colorbar;
title('\sigma_{xx}'); xlabel('x'); ylabel('y'); zlabel('\sigma_{xx}');
axis tight; daspect([1 1 0.35]); view(45,30);

nexttile;
surf(X,Y,sigma_yy); shading interp; colormap(cm); colorbar;
title('\sigma_{yy}'); xlabel('x'); ylabel('y'); zlabel('\sigma_{yy}');
axis tight; daspect([1 1 0.35]); view(45,30);

nexttile;
surf(X,Y,sigma_zz); shading interp; colormap(cm); colorbar;
title('\sigma_{zz}'); xlabel('x'); ylabel('y'); zlabel('\sigma_{zz}');
axis tight; daspect([1 1 0.35]); view(45,30);
sgtitle('Distribución de tensiones normales (\lambda=\mu=1)','FontWeight','bold');

% Tensión tangencial respecto a i 
figure('Name','Tension tangencial respecto a i','Color','w');
surf(X,Y,tau_i); shading interp; colormap(cm); colorbar;
title('\tau_{i} = |\sigma_{xy}|'); xlabel('x'); ylabel('y'); zlabel('\tau_{i}');
axis tight; daspect([1 1 0.35]); view(45,30); grid on;

% Mmapas en planta
figure('Name','Contornos 2D de tensiones','Color','w');
tiledlayout(1,4,'Padding','compact','TileSpacing','compact');

nexttile; contourf(X,Y,sigma_xx,20,'LineColor','none'); axis equal tight; colorbar;
title('\sigma_{xx}'); xlabel('x'); ylabel('y');

nexttile; contourf(X,Y,sigma_yy,20,'LineColor','none'); axis equal tight; colorbar;
title('\sigma_{yy}'); xlabel('x'); ylabel('y');

nexttile; contourf(X,Y,sigma_zz,20,'LineColor','none'); axis equal tight; colorbar;
title('\sigma_{zz}'); xlabel('x'); ylabel('y');

nexttile; contourf(X,Y,tau_i,20,'LineColor','none'); axis equal tight; colorbar;
title('\tau_{i}'); xlabel('x'); ylabel('y');

% Definimos la tensión tangencial en el plano ortogonal al eje i (X), que es tau = [0, sigma_yx, 0].
Tau_x = zeros(size(sigma_xy)); % (debe ser nula).
Tau_y = sigma_yx;            

% La magnitud (tau_i_magnitude) se mantiene aquí, pero no se utiliza en el gráfico final.
tau_i_magnitude = abs(sigma_yx);

% Creamos la figura y configuración inicial
figure('Name','Tensión Tangencial ','Color','w');
hold on; 
axis equal; 
box on; 
grid on;
axis([-0.5 4.5 -0.5 2.5]);
title('Tensiones tangenciales respecto al plano ortogonal a i');
xlabel('x'); 
ylabel('y');

% El vector tau es vertical: (0,-x/20), según los cálculos del Ejercicio 9.
step = 3; 
idx_x = 1:step:size(X,2);
idx_y = 1:step:size(X,1);

% Usar el comando quiver para dibujar las flechas de tensión.
quiver(X(idx_y, idx_x), Y(idx_y, idx_x),Tau_x(idx_y, idx_x),Tau_y(idx_y, idx_x),1.5,'r', 'LineWidth', 1.2);
% Cortorno de la placa 
plot(x, x/8, 'k-', 'LineWidth', 2);    
plot(x, 2 - x/8, 'k-', 'LineWidth', 2);
plot([0 0], [0 2], 'k-', 'LineWidth', 2);     
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); 

hold off;

clc; clear; close all;

% Malla
h = 0.1;
[x, y] = meshgrid(0:h:4, 0:h:2);   % ojo: x -> columnas, y -> filas
X = x; Y = y;

% Bordes de la placa
f = X./8;
g = 2 - X./8;
EnPlaca = (Y >= f) & (Y <= g);

% Definición del campo de desplazamientos en cartesianas
ux=(1/20).*X.*Y;
uy=-(1/20).*X.^2;

% Derivadas analíticas:
% ux_x=d(u_x)/dx=y/20
% ux_y=d(u_x)/dy=x/20
% uy_x=d(u_y)/dx=-x/10
% uy_y=d(u_y)/dy 0
ux_x=Y./20;
ux_y=X./20;
uy_x=-X./10;
uy_y=zeros(size(X));

% Divergencia (útil para sigma)
divU=ux_x+uy_y;

% Tensor de deformaciones simétrico 
% epsilon = 1/2(grad u + grad u^T)
% epsilon_xx = ux_x
% epsilon_yy = uy_y
% epsilon_xy = 1/2*(ux_y + uy_x)
eps_xx = ux_x;
eps_yy = uy_y;
eps_xy = 0.5*(ux_y + uy_x);

% Para el caso 3D consideramos eps_zz = 0
eps_zz = zeros(size(X));

% Coeficientes de Lamé
lambda=1;
mu=1;

% Construcción del tensor de tensiones 
%sigma = lambda*div(u)*I + 2*mu*eps
sigma_xx = lambda.*divU + 2*mu.*eps_xx;
sigma_yy = lambda.*divU + 2*mu.*eps_yy;
sigma_zz = lambda.*divU + 2*mu.*eps_zz;
sigma_xy = 2*mu.*eps_xy;    % = sigma_yx

% Von Mises y matriz para autovalores
VM = NaN(size(X));
S1 = NaN(size(X)); S2 = NaN(size(X)); S3 = NaN(size(X));

% Loop sobre puntos dentro de la placa para calcular autovalores
[idx_i,idx_j]=find(EnPlaca);
npts=numel(idx_i);

for k = 1:npts
    i = idx_i(k); j = idx_j(k);   % i: fila (y index), j: columna (x index)
    % construir sigma 3x3 en ese punto
    S = [sigma_xx(i,j), sigma_xy(i,j), 0;sigma_xy(i,j), sigma_yy(i,j),0;0,0,sigma_zz(i,j)];
    e = eig(S);          
    % ordenar autovalores (no estrictamente necesario, pero ordeno descendente)
    e = sort(e,'descend');
    s1 = e(1); s2 = e(2); s3 = e(3);
    % tensión de Von Mises (fórmula)
    sigmaVM = sqrt( ((s1-s2).^2 + (s2-s3).^2 + (s3-s1).^2) / 2 );
    VM(i,j) = sigmaVM;
    S1(i,j) = s1; S2(i,j) = s2; S3(i,j) = s3;
end

VM(~EnPlaca) = NaN;
% Buscar máximo de Von Mises dentro de la placa
[maxVM, indMax] = max(VM(:));
[rowMax, colMax] = ind2sub(size(VM), indMax);
xMax = X(rowMax, colMax);
yMax = Y(rowMax, colMax);

% Visualización
figure('Name','Tensión de Von Mises','Color','w'); clf; hold on;
axis equal; axis([-0.5 4.5 -0.5 2.5]); box on;
title('Tensión de Von Mises \sigma_{VM}');
xlabel('x'); ylabel('y');

% Contourf del campo VM
numLevels = 30;
[C, hC] = contourf(X, Y, VM, numLevels, 'LineStyle','none');
colormap(jet);
c = colorbar; ylabel(c, '\sigma_{VM}');
% Bordes de la placa
xb = linspace(0,4,200);
plot(xb, xb./8, 'k-', 'LineWidth', 1.5);         
plot(xb, 2 - xb./8, 'k-', 'LineWidth', 1.5);    
plot([0 0], [0 2], 'k-', 'LineWidth', 1.5);  
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 1.5); 

% Punto de máxima tensión de Von Mises
plot(xMax, yMax, 'rp', 'MarkerSize', 14, 'MarkerFaceColor', 'r');
text(xMax+0.08, yMax+0.08, sprintf('Máx VM = %.3g', maxVM), 'FontWeight','bold','Color','r');

% Mejora estética
set(gcf,'Color','w');
shading interp;

hold off;

% Información
fprintf('Máximo Von Mises = %.6g en (x,y) = (%.4g, %.4g)\n', maxVM, xMax, yMax);

 
% Definición del paso de muestreo (h) y el vector de coordenadas X.
h = 0.1;                         
x = 0:h:4;                       
num_puntos_y = round(2/h);
s = linspace(0,1,num_puntos_y+1);

% Generación del mallado físico (X, Y) dentro de la geometría trapezoidal.
[U,S] = meshgrid(x,s);
f = U./8;
g = 2 - U./8;
X = U;
Y = f + S.*(g - f);

% Creación de la máscara lógica para identificar los puntos interiores de la placa.
mask = (Y >= X/8) & (Y <= 2 - X./8);

% CAMPO DE FUERZA F
% Cálculo analítico: F = -div(sigma) resulta en F = [0, 0].
% Inicializa la componente horizontal Fx a cero.
Fx = zeros(size(X));

% Inicializa la componente vertical Fy a cero.
Fy = zeros(size(X));

% Aplica la máscara: los puntos fuera de la placa se establecen a NaN (no se dibujan).
Fx(~mask) = NaN;
Fy(~mask) = NaN;

% Configuración inicial de la figura y los ejes para la visualización.
figure('Name','Campo de Fuerzas F','Color','w');
hold on; 
axis equal; 
box on; 
grid on;
axis([-0.5 4.5 -0.5 2.5]);
title('Campo de Fuerzas Volumétricas ');
xlabel('x'); 
ylabel('y');

% DIBUJO DEL CAMPO VECTORIAL
step = 3; 
idx_x = 1:step:size(X,2);
idx_y = 1:step:size(X,1);
quiver(X(idx_y, idx_x), Y(idx_y, idx_x),Fx(idx_y, idx_x),Fy(idx_y, idx_x),0.5,'b', 'LineWidth', 1.2);

%Resultado analítico de F=0.
text(2, 1, '¡El campo de fuerzas nulo!','HorizontalAlignment', 'center', 'FontSize', 14, 'Color', [0.5 0 0], 'FontWeight', 'bold');
text(2, 0.7, 'El equilibrio se satisface en el interior de la placa.','HorizontalAlignment', 'center', 'FontSize', 10, 'Color', [0.5 0 0]);

% Dibuja las líneas que delimitan el contorno de la placa.
plot(x, x/8, 'k-', 'LineWidth', 2);
plot(x, 2 - x/8, 'k-', 'LineWidth', 2);
plot([0 0], [0 2], 'k-', 'LineWidth', 2);
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2);
hold off;

clc; clear; close all;
% Mallado
h=0.05; 
[X,Y]=meshgrid(0:h:4,0:h:2);

% Forma de la placa
Condition=(Y>=X./8)&(Y<=2-X./8);

% Definición de la Densidad : d(x,y)=(4-x)*|y|
Densidad=(4-X).*abs(Y);

% Filtrado para la Gráfica
Densidad(~Condition)=NaN;

% Visualización 3D
figure('Name','Densidad', 'Color','w');

% Superficie 3D
surf(X,Y,Densidad);

% Configuración
shading interp; % Elimina las líneas negras para que se vea suave
colormap(jet);  % Mapa de colores (De azul a rojo)
colorbar;       % Barra de referencia
caxis([0 8]);   % Fija la escala de colores

% Etiquetas y Título
title({'Densidad'},'FontSize',13);
xlabel('Eje X','FontWeight','bold');
ylabel('Eje Y','FontWeight','bold');
zlabel('Densidad (kg/m^2)','FontWeight','bold');

% Ajustes
axis tight;      % Ajusta los ejes al contenido
view(-45, 30);   % Ángulo de vista 
grid on;
box on;

% Proyección de sombra en el suelo
hold on;
% Dibujamos el contorno plano abajo del todo (en z=0) para referencia
contour(X,Y,Densidad,20,'z','Offset',0); 
hold off;