Diferencia entre revisiones de «Placa Plana (Grupo 20)»

% 1. Configuración
h = 1/10;
color_malla = [0, 0.6, 0.6]; % Color Cían Oscuro

% 2. Generar Coordenadas Paramétricas (usamos estas en vez de x e y para que quede mejor ajustado a la figura)
u = 0 : h : 4;
v = 0 : h : 2;
[U, V] = meshgrid(u, v);

% 3. Transformación Matemática
% Calculamos los bordes superior e inferior para cada punto X
y_abajo = U ./ 8;
y_arriba = 2 - (U ./ 8);

% Interpolamos: calculamos la Y real basada en la altura V
factor_altura = V ./ 2; % Va de 0 (abajo) a 1 (arriba)

X = U;
Y = y_abajo + factor_altura .* (y_arriba - y_abajo);

% 4. Visualización
figure('Color', 'w'); hold on;

% A) Dibujar las líneas VERTICALES
plot(X, Y, 'Color', color_malla, 'LineWidth', 0.5);

% B) Dibujar las líneas HORIZONTALES
plot(X', Y', 'Color', color_malla, 'LineWidth', 0.5);

% C) Dibujar el contorno 
plot(u, u./8, 'k-', 'LineWidth', 2); % Borde Abajo
plot(u, 2 - u./8, 'k-', 'LineWidth', 2); % Borde Arriba
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Borde Izquierdo
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Borde Derecho

% 5. Configuración final
axis([-1 5 -1 3]); % Zoom/Encuadre exacto
xlabel('x'); ylabel('y');
title('Mallado de la placa plana');
grid on;
box on;
hold off;

% 1. Definición del Dominio y Mallado
h = 0.2; % Paso de muestreo
x = 0 : h : 4;
y = 0 : h : 2;
[X, Y] = meshgrid(x, y);

% 2. Definir la Geometría 
y_abajo = X ./ 8;
y_arriba = 2 - (X ./ 8);

% Creamos la máscara (1 dentro, 0 fuera)
dentro = (Y >= y_abajo) & (Y <= y_arriba);

% 3. Definir la TEMPERATURA T(x,y)
% T(x,y) = (1 + (y-1)^2)(4-x)
T_func = @(x,y) (1 + (y - 1).^2) .* (4 - x);

% Calculamos la temperatura en toda la malla
T = T_func(X, Y);

% 4. Aplicar la Máscara 
T(~dentro) = NaN;
X_plot = X; Y_plot = Y;
X_plot(~dentro) = NaN;
Y_plot(~dentro) = NaN;

% 5. Cálculo del Gradiente (Nabla T)
[dTdx, dTdy] = gradient(T, h, h);

% 6. Visualización
figure('Color', 'w', 'Position', [100 100 800 600]);
hold on;

% A) Pintar Curvas de Nivel (Isotermas)
% Usamos contour para ver las líneas de temperatura constante
[C, h_cont] = contour(X, Y, T, 20, 'LineWidth', 2);
clabel(C, h_cont, 'FontSize', 9, 'Color', 'k', 'FontWeight', 'bold');
colormap('jet');

% B) Pintar Campo Vectorial (Gradiente)
% quiver dibuja las flechas sobre la superficie
quiver(X, Y, dTdx, dTdy, 'k', 'LineWidth', 1.2);

% C) Dibujar el contorno de la placa
plot([0 4], [0 0.5], 'k-', 'LineWidth', 2); % Abajo
plot([0 4], [2 1.5], 'k-', 'LineWidth', 2); % Arriba
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Izquierda
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Derecha

% 7. Configuración Final
axis equal; 
axis([-0.5 4.5 -0.5 2.5]);
grid on;
xlabel('x'); ylabel('y');
title({'Isotermas y Gradiente \nablaT', ...
'T(x,y) = (1 + (y-1)^2)(4-x)'});
colorbar;

%  CÁLCULO DEL MÁXIMO 
[max_temp, indice] = max(T(:));
[fil, col] = ind2sub(size(T), indice);
x_max = X(fil, col);
y_max = Y(fil, col);

fprintf('La temperatura máxima es %.2f en el punto (x=%.2f, y=%.2f)\n', ...
max_temp, x_max, y_max);
% Marcamos el punto máximo con una estrella roja
plot(x_max, y_max, 'p', 'MarkerFaceColor', 'r', 'MarkerEdgeColor', 'k', 'MarkerSize', 15);
hold off;

% 1. Configuración del Mallado
h = 0.2;
x = 0 : h : 4;
y = 0 : h : 2;
[X, Y] = meshgrid(x, y);

% 2. Definir la Geometría 
y_abajo = X ./ 8;
y_arriba = 2 - (X ./ 8);

% Máscara lógica: determina qué puntos están dentro de la placa
dentro = (Y >= y_abajo) & (Y <= y_arriba);

% 3. Definir la Temperatura T(x,y)
T_func = @(x,y) (1 + (y - 1).^2) .* (4 - x);
T = T_func(X, Y);

% 4. Cálculos Físicos (Ley de Fourier)
% A) Calculamos el Gradiente (nabla T) mediante derivadas numéricas
[dTdx, dTdy] = gradient(T, h, h);

% B) Calculamos el Flujo de Calor Q
% Q = -k * nabla T. (Suponiendo k=1) -> Q apunta al frío
Qx = -dTdx;
Qy = -dTdy;

% 5. Limpieza de datos (Recortar la figura)
% Asignamos NaN (Not a Number) a todo lo que esté fuera de la placa.
% Esto hace que MATLAB no dibuje nada en esas zonas.
Qx(~dentro) = NaN;
Qy(~dentro) = NaN;
T(~dentro) = NaN;

% 6. Visualización
figure('Color', 'w', 'Name', 'Flujo de Calor Fourier');
hold on;

% A) FONDO: Mapa de Calor (contourf)
% 'LineStyle','none' quita las rayitas negras entre colores para que sea suave.
contourf(X, Y, T, 20, 'LineStyle', 'none');
colormap('jet'); % Paleta de colores: Azul (frío) -> Rojo (caliente)

% Barra de color explicativa
c = colorbar;
c.Label.String = 'Temperatura T(x,y)';
c.Label.FontSize = 10;

% B)El Campo Vectorial Q.
quiver(X, Y, Qx, Qy, 'k', 'LineWidth', 1.2, 'AutoScaleFactor', 1.5);

% C) Dibujo del contorno.
plot([0 4], [0 0.5], 'k-', 'LineWidth', 2); % Borde inferior
plot([0 4], [2 1.5], 'k-', 'LineWidth', 2); % Borde superior
plot([0 0], [0 2], 'k-', 'LineWidth', 2); % Borde izquierdo
plot([4 4], [0.5 1.5], 'k-', 'LineWidth', 2); % Borde derecho

% 7. Configuración Final
axis equal;
axis([-0.5 4.5 -0.5 2.5]); 
grid on;
xlabel('Eje x'); ylabel('Eje y');
title({'Vector Flujo de Calor', ...
'El calor fluye desde las capas más calientes hacia las más frías'});

hold off;

% 1. Definición del Mallado y la Región
h = 0.05; % Paso del mallado (puedes ajustarlo a h=0.1 o h=0.01)
x = 0:h:4;
y = 0:h:2; % El rango de Y va de 0 a 2 según las funciones f(x) y g(x)
[X, Y] = meshgrid(x, y);

% Definición de los límites de la placa (Enunciado J)
% f(x) = x/8;
% g(x) = 2 - x/8;
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos la máscara (Ponemos NaN fuera de la placa)
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo de las Derivadas Parciales (Gradiente)
% Usamos las fórmulas analíticas derivadas de T(x,y)
DerivadaX = -(1 + (Y - 1).^2);
DerivadaY = 2 .* (4 - X) .* (Y - 1);

% 3. Cálculo del Módulo del Gradiente
GradienteModulo = sqrt(DerivadaX.^2 + DerivadaY.^2);

% 4. Encontrar el Máximo Global
% Buscamos el valor máximo ignorando los NaNs
[MaxGradiente, IndiceLineal] = max(GradienteModulo);

% Convertimos el índice lineal a coordenadas de matriz (fila, columna)
[fila, col] = ind2sub(size(GradienteModulo), IndiceLineal);

% Obtenemos las coordenadas (x, y) y las componentes del vector en ese punto
PuntoMax = [X(fila, col), Y(fila, col)];
DirMax = [DerivadaX(fila, col), DerivadaY(fila, col)];

% Mostramos los resultados en la consola
fprintf('La variación máxima es: %.4f\n', MaxGradiente);
fprintf('Ocurre en el punto: (%.2f, %.2f)\n', PuntoMax(1), PuntoMax(2));
fprintf('Dirección del gradiente (vector): [%.2f, %.2f]\n', DirMax(1), DirMax(2));

% 5. Gráfica
figure;
hold on;
axis equal;
grid on;
title('Punto de Máxima Variación de Temperatura');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]); 

% Dibujar el contorno de la placa 
x_borde = 0:0.1:4;
plot(x_borde, x_borde/8, 'k', 'LineWidth', 1.5); % f(x)
plot(x_borde, 2 - x_borde/8, 'k', 'LineWidth', 1.5); % g(x)
line([0 0], [0 2], 'Color', 'k', 'LineWidth', 1.5); % Pared izquierda
line([4 4], [0.5 1.5], 'Color', 'k', 'LineWidth', 1.5); % Borde derecho

% A) Dibujar el punto rojo
plot3(PuntoMax(1), PuntoMax(2), 10, 'ro', 'MarkerSize', 10, 'LineWidth', 3, 'MarkerFaceColor', 'r');

% B) Dibujar el vector de dirección en ese punto (Quiver)
% Normalizamos el vector para que se vea bien en la gráfica
escala = 0.5;
quiver(PuntoMax(1), PuntoMax(2), DirMax(1), DirMax(2), ...
escala, 'r', 'LineWidth', 2, 'MaxHeadSize', 0.5);

text(PuntoMax(1)+0.2, PuntoMax(2), sprintf('Max: %.2f', MaxGradiente));

hold off;

h = 0.1; 
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la región 
Condicion = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos la máscara )
X(~Condicion) = NaN;
Y(~Condicion) = NaN;

% 2. Cálculo del Campo de Desplazamientos (u_x, u_y)

% Primero necesitamos las coordenadas polares de cada punto (rho, theta)
Rho = sqrt(X.^2 + Y.^2);
Theta = atan2(Y, X);

% Definimos la componente escalar u_theta dada por el enunciado
% u(rho, theta) = -1/20 * rho^2 * cos(theta) * e_theta
U_theta_val = -1/20 .* Rho.^2 .* cos(Theta);

% Convertimos el vector de base polar a base cartesiana
% u_x = -u_theta * sin(theta)
% u_y = u_theta * cos(theta)
Ux = -U_theta_val .* sin(Theta);
Uy = U_theta_val .* cos(Theta);


% 3. Calcular la posición deformada

% Posición final = Posición inicial + Desplazamiento
X_final = X + Ux;
Y_final = Y + Uy;


% 4. Gráfica Comparativa 

figure('Name', 'Deformación de la Placa', 'Color', 'w');

% PLACA ORIGINAL
subplot(1, 2, 1);
% Usamos mesh con Z=0 para dibujar la 'rejilla' plana
mesh(X, Y, zeros(size(X)), 'EdgeColor', 'b');
view(2); % Vista superior (2D)
axis equal; grid on;
title('Placa antes del desplazamiento');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]);

% PLACA DEFORMADA 
subplot(1, 2, 2);
mesh(X_final, Y_final, zeros(size(X)), 'EdgeColor', 'b');
view(2);
axis equal; grid on;
title('Placa después del desplazamiento');
xlabel('Eje X'); ylabel('Eje Y');
axis([-0.5 4.5 -0.5 2.5]);


 Apartado 5 

figure('Name', 'Campo Vectorial de Desplazamientos');
hold on; axis equal; grid on;
title('Campo de Desplazamientos');

% Dibujamos contorno original de referencia
plot(X(:), Y(:), '.k', 'MarkerSize', 1);

% Dibujamos los vectores (quiver)
quiver(X, Y, Ux, Uy, 0, 'r', 'LineWidth', 1.5);

xlabel('Eje X'); ylabel('Eje Y');
xlim([-0.5 4.5]); ylim([-0.5 2.5]);

% 1. Definición de Geometría y Desplazamientos
h = 0.1; % Paso del mallado
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la Región (
% La placa existe si Y >= X/8 Y Y <= 2 - X/8
Condicion = (Y >= X./8) & (Y <= 2 - X./8);

% Aplicamos máscara 
X(~Condicion) = NaN;
Y(~Condicion) = NaN;

% Cálculo de Desplazamientos
% Coordenadas polares necesarias para la fórmula
Rho = sqrt(X.^2 + Y.^2);
Theta = atan2(Y, X);

% Valor escalar de u_theta
U_theta_val = -1/20 .* Rho.^2 .* cos(Theta);

% Proyección a cartesianas 
Ux = -U_theta_val .* sin(Theta);
Uy = U_theta_val .* cos(Theta);

% --- Posición Final Deformada ---
X_final = X + Ux;
Y_final = Y + Uy;

% 2. Gráfica Comparativa de Sólidos Mallados 
figure('Name', 'Comparación Antes/Después', 'Color', 'w', 'Position', [100, 100, 1000, 500]);

% SUBPLOT 1: Placa Original
subplot(1, 2, 1);
hold on;

mesh(X, Y, zeros(size(X)), 'EdgeColor', [0 0.7 0.9], 'FaceColor', 'none');

view(2); % Fuerza la vista superior 2D
axis equal; 
grid on;
box on; 

title('Placa antes del desplazamiento', 'FontWeight', 'bold');
xlabel('Eje X');
ylabel('Eje Y');
axis([-1.5 4.5 0 2.5]);

% SUBPLOT 2: Placa Deformada
subplot(1, 2, 2);
hold on;
% Dibujamos la malla usando las coordenadas finales (X_final, Y_final)
mesh(X_final, Y_final, zeros(size(X)), 'EdgeColor', [0 0.7 0.9], 'FaceColor', 'none');

view(2);
axis equal;
grid on;
box on;

title('Placa después del desplazamiento', 'FontWeight', 'bold');
xlabel('Eje X');
ylabel('Eje Y');
axis([-1.5 4.5 0 2.5]); 

sgtitle('Visualización del Sólido Mallado bajo Deformación'); % Título

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa 
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo Analítico de la Divergencia
Divergencia = Y ./ 20;

% 3. Encontrar Máximos, Mínimos y Nulos analíticamente en el mallado
max_val = max(Divergencia(:));
min_val = min(Divergencia(:));

% Encontramos los índices donde ocurren
[f_max, c_max] = find(Divergencia == max_val);
[f_min, c_min] = find(Divergencia == min_val);

% Puntos (x,y) correspondientes
PuntosMax = [X(f_max), Y(f_max)];
PuntosMin = [X(f_min), Y(f_min)];

% 4. Mostrar resultados en consola
fprintf('--- Resultados del Análisis de Divergencia ---\n');
fprintf('Divergencia Máxima: %.4f (Ocurre en y máximas)\n', max_val);
fprintf('Divergencia Mínima: %.4f (Ocurre en y = 0)\n', min_val);
fprintf('La divergencia es nula solo donde y = 0.\n');

% 5. Gráfica de la Divergencia
figure('Name', 'Divergencia del Desplazamiento', 'Color', 'w');
hold on;

% Dibujamos el mapa de colores
contourf(X, Y, Divergencia, 20, 'LineColor', 'none');
colormap(jet); % Mapa de colores clásico (Azul=Bajo, Rojo=Alto)
colorbar; % Barra lateral para saber qué valor es cada color

% Dibujar contorno de la placa para que quede bonito
x_borde = 0:0.1:4;
plot(x_borde, x_borde/8, 'k', 'LineWidth', 2);
plot(x_borde, 2 - x_borde/8, 'k', 'LineWidth', 2);
line([0 0], [0 2], 'Color', 'k', 'LineWidth', 2);

% Señalar Máximo y Mínimo
plot(PuntosMax(:,1), PuntosMax(:,2), 'rp', 'MarkerSize', 12, 'MarkerFaceColor', 'r'); % Estrella roja
plot(PuntosMin(:,1), PuntosMin(:,2), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b'); % Círculo azul

% Etiquetas y títulos
title('Cambio de Volumen Local');
xlabel('Eje X'); ylabel('Eje Y');
axis equal;
axis([-0.5 4.5 -0.5 2.5]);

% Leyenda manual para entender los puntos
legend('Divergencia (Color)', 'Bordes', '', '', 'Máxima Expansión', 'Mínima/Nula', ...
'Location', 'eastoutside');

hold off;

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 1. Definición del Mallado y Región
h = 0.05;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo del Rotacional (Magnitud)
%

%% Apartado 9 (Versión 2D): Tensiones Normales
clc; clear; close all;

% 1. Definición de Geometría y Material
h = 0.1;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Placa J (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN; Y(~Condition) = NaN;

lambda = 1; mu = 1;

% 2. Cálculo de Componentes
% Deformaciones y Divergencia
Div_u = Y ./ 20;
Eps_xx = Y ./ 20;
Eps_yy = zeros(size(X));
% Eps_zz es 0

% Tensiones Normales (Sigma = lambda*Div*I + 2*mu*Eps)
Sig_xx = lambda .* Div_u + 2 .* mu .* Eps_xx; % 3y/20
Sig_yy = lambda .* Div_u + 2 .* mu .* Eps_yy; % y/20
Sig_zz = lambda .* Div_u; % y/20 (Por efecto Poisson)
Sig_xy = 2 .* mu .* (-X./40); % -x/20 (Cizalladura)

% 3. Gráfica 2D (Subplot)
figure('Name', 'Apartado 9: Tensiones Normales (2D)', 'Color', 'w', 'Position', [100, 100, 1200, 400]);

% Sigma XX
subplot(1, 3, 1);
contourf(X, Y, Sig_xx, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{xx} (Eje i)'); xlabel('X'); ylabel('Y');

% Sigma YY
subplot(1, 3, 2);
contourf(X, Y, Sig_yy, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{yy} (Eje j)'); xlabel('X'); ylabel('Y');

% Sigma ZZ
subplot(1, 3, 3);
contourf(X, Y, Sig_zz, 20, 'LineColor', 'none');
axis equal; axis([-0.5 4.5 -0.5 2.5]); colorbar;
title('\sigma_{zz} (Eje k)'); xlabel('X'); ylabel('Y');

sgtitle('Componentes Normales del Tensor de Tensiones');

% 1. Definición del Mallado
h = 0.15; % Un paso un poco más grande para que las flechas no se amontonen
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% Definición de la placa 
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

% 2. Cálculo de la Tensión Tangencial (Sigma_xy)
% Recordamos del apdo 9: Sig_xy = 2*mu*eps_xy = 2*1*(-x/40) = -x/20
Sig_xy = -X ./ 20;

% 3. Definición del Vector Tangencial
% En la cara vertical (normal i), la tangencial es vertical (j).
% El vector es: (0, Sigma_xy)
U_vec = zeros(size(X)); % Componente X es 0
V_vec = Sig_xy; % Componente Y es el valor de la tensión

% 4. Gráfica
figure('Name', 'Apartado 10: Vector Tensión Tangencial', 'Color', 'w');

% Usamos quiver para dibujar las flechitas azules
% quiver(x, y, u, v)
quiver(X, Y, U_vec, V_vec, 'b', 'LineWidth', 1);

% Configuración de la vista
view(2); % Vista 2D superior
axis equal; % Ejes proporcionales
grid on; % Cuadrícula de fondo

% Títulos y Etiquetas
title('Tensiones tangenciales (Plano ortogonal a i)');
xlabel('Eje X');
ylabel('Eje Y');

% Ajuste de límites para ver bien la placa
axis([-0.5 4.5 -0.5 2.5]);

%% Apartado 11: Tensión de Von Mises (Usando autovalores 'eig')
% Inicializamos la matriz de resultados con NaNs
VonMises_Calc = NaN(size(X));
PuntoMax_VM = [-1, -1];
Max_VM_Val = -1;

% Bucle para recorrer cada punto del mallado
[filas, cols] = size(X);
for f = 1:filas
for c = 1:cols
% Si el punto está dentro de la placa (no es NaN)
if ~isnan(X(f,c))
% 1. Construir el tensor de tensiones en ese punto (3x3)
s_xx = Sig_xx(f,c);
s_yy = Sig_yy(f,c);
s_zz = Sig_zz(f,c);
s_xy = Sig_xy(f,c);

Tensor = [s_xx, s_xy, 0;
s_xy, s_yy, 0;
0, 0, s_zz];

% 2. Calcular autovalores (Tensiones Principales)
sig_principales = eig(Tensor);
s1 = sig_principales(1);
s2 = sig_principales(2);
s3 = sig_principales(3);

% 3. Fórmula de Von Mises dada en el enunciado
vm = sqrt( ((s1-s2)^2 + (s2-s3)^2 + (s3-s1)^2) / 2 );
VonMises_Calc(f,c) = vm;

% Guardar máximo para señalarlo después
if vm > Max_VM_Val
Max_VM_Val = vm;
PuntoMax_VM = [X(f,c), Y(f,c)];
end
end
end
end

% Gráfica
figure('Name', 'Apartado 11: Von Mises con Autovalores', 'Color', 'w');
contourf(X, Y, VonMises_Calc, 30, 'LineColor', 'none');
colormap(hot); % 'hot' suele usarse para mapas de "calor/tensión"
c = colorbar; c.Label.String = 'Tensión Von Mises (Pa)';

hold on;
axis equal; axis([-0.5 4.5 -0.5 2.5]);
title('Tensión de Von Mises (Cálculo mediante Autovalores)');
xlabel('Eje X'); ylabel('Eje Y');

% Señalar el punto máximo
plot(PuntoMax_VM(1), PuntoMax_VM(2), 'bp', 'MarkerSize', 14, 'MarkerFaceColor', 'b'); % Estrella azul
text(PuntoMax_VM(1)-1, PuntoMax_VM(2), sprintf('Máx: %.4f en (%.1f, %.1f)', Max_VM_Val, PuntoMax_VM), 'Color', 'b', 'FontWeight', 'bold');
hold off;

%% Apartado 12: Campo de Fuerzas Volumétricas (F = -div Sigma)
clc; clear; close all;

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

% Definición de la placa (Trapecio)
Condition = (Y >= X./8) & (Y <= 2 - X./8);
X(~Condition) = NaN;
Y(~Condition) = NaN;

lambda = 1; mu = 1;

% 2. Recalculamos las Tensiones (Necesarias para derivar)
% Usamos las fórmulas analíticas que ya conocemos
Div_u = Y ./ 20;
Eps_xx = Y ./ 20;
Eps_yy = zeros(size(X));

Sig_xx = lambda .* Div_u + 2 .* mu .* Eps_xx; % = 3y/20
Sig_yy = lambda .* Div_u + 2 .* mu .* Eps_yy; % = y/20
Sig_xy = -X ./ 20; % = -x/20 (cizalladura)

% 3. Cálculo de la Divergencia del Tensor (Numéricamente)
% F = -div(sigma)
% Fx = -( d(Sxx)/dx + d(Sxy)/dy )
% Fy = -( d(Sxy)/dx + d(Syy)/dy )

% La función gradient calcula derivadas centrales.
% Sintaxis: gradient(Matriz, paso_x, paso_y)
[dSxx_dx, dSxx_dy] = gradient(Sig_xx, h, h);
[dSxy_dx, dSxy_dy] = gradient(Sig_xy, h, h);
[dSyy_dx, dSyy_dy] = gradient(Sig_yy, h, h);

Fx = - (dSxx_dx + dSxy_dy);
Fy = - (dSxy_dx + dSyy_dy);

% 4. Comprobación de Magnitud (Para ver si es cero)
MagnitudF = sqrt(Fx.^2 + Fy.^2);
max_F = max(MagnitudF(:));

fprintf('--- Resultado Apartado 12 ---\n');
fprintf('La fuerza máxima calculada es: %e\n', max_F);

if max_F < 1e-10
fprintf('CONCLUSIÓN: La fuerza es prácticamente NULA (Equilibrio).\n');
subtitle_text = 'Fuerza \approx 0 (Equilibrio Estático)';
else
subtitle_text = 'Existen fuerzas volumétricas';
end

% 5. Gráfica
figure('Name', 'Apartado 12: Campo de Fuerzas', 'Color', 'w');
hold on;

% Dibujamos los vectores.
% Como son casi cero, MATLAB pintará puntos. Para ver "algo",
% a veces se normalizan, pero lo correcto es mostrar que no hay fuerza.
quiver(X, Y, Fx, Fy, 'r');

% Dibujar bordes
plot([0 4 4 0 0], [0 0.5 1.5 2 0], 'k--', 'LineWidth', 0.5); % Caja aproximada

title('Campo de Fuerzas Volumétricas');
subtitle(subtitle_text);
xlabel('Eje X'); ylabel('Eje Y');
axis equal; axis([-0.5 4.5 -0.5 2.5]);
grid on;

% Añadimos una nota en la gráfica explicándolo
text(1, 1, {'Fuerza Nula', '(puntos indican magnitud ~0)'}, ...
'Color', 'r', 'FontWeight', 'bold', 'HorizontalAlignment', 'center');

hold off;

%% Apartado 13: Cálculo de la Masa Total
clc; clear; close all;

% 1. Definición del Mallado
% Usamos un paso fino (h=0.01) para que la integral sea precisa
h = 0.01;
x = 0:h:4;
y = 0:h:2;
[X, Y] = meshgrid(x, y);

% 2. Definición de la Placa (Región)
Condition = (Y >= X./8) & (Y <= 2 - X./8);

% 3. Definición de la Densidad
% Función dada: d(x,y) = (4 - x) *