Diferencia entre revisiones de «Flujo estacionario en medio poroso (RoYa)»
De MateWiki
(→Códigos de MatLab) |
|||
| Línea 5: | Línea 5: | ||
=Códigos de MatLab= | =Códigos de MatLab= | ||
| + | a = 0; | ||
| + | b = 1; | ||
| + | |||
| + | H = [0.1 0.05 0.025]; % Pasos considerados | ||
| + | |||
| + | V = {'Vertices_malla1_tipo1.txt', 'Vertices_malla2_tipo1.txt', 'Vertices_malla3_tipo1.txt'}; | ||
| + | E = {'Elementos_malla1_tipo1.txt', 'Elementos_malla2_tipo1.txt', 'Elementos_malla3_tipo1.txt'}; | ||
| + | D = {'Dirichlet_malla1_tipo1.txt', 'Dirichlet_malla2_tipo1.txt', 'Dirichlet_malla3_tipo1.txt'}; | ||
| + | R = {'Relacion_malla1malla3.txt', 'Relacion_malla2malla3.txt'}; | ||
| + | T = {'Mallado con h=0.1', 'Mallado con h=0.05', 'Mallado con h=0.025'}; | ||
| + | |||
| + | % Datos fijos necesarios dentro del bucle | ||
| + | h = H(3); | ||
| + | X = a:h:b; | ||
| + | Xm = a+(h/2):h:b; | ||
| + | n2 = length(X); | ||
| + | m2 = length(Xm); | ||
| + | l2 = n2+m2; | ||
| + | |||
| + | |||
| + | for p = 1:3 | ||
| + | h = H(p); | ||
| + | x = a:h:b; | ||
| + | xm = a+(h/2):h:b; | ||
| + | n = length(x); | ||
| + | m = length(xm); | ||
| + | apoyo = zeros(n+m,1); | ||
| + | for i = 1:n | ||
| + | apoyo(2*(i-1)+1) = x(i); | ||
| + | if i <= m | ||
| + | apoyo(2*i) = xm(i); | ||
| + | end | ||
| + | end | ||
| + | l = length(apoyo); | ||
| + | |||
| + | Verticesdoc = fopen(V{p},'w'); | ||
| + | formatSpec = '%d %.4f %.4f \n'; | ||
| + | for i = 1:l | ||
| + | if mod(i,2) == 1 | ||
| + | for j = 1:n | ||
| + | nodo = [(((n+m)*(floor(i/2))) + j), x(j), apoyo(i)]; | ||
| + | fprintf(Verticesdoc,formatSpec,nodo); | ||
| + | end | ||
| + | else | ||
| + | for j = 1:m | ||
| + | nodo = [(((n+m)*(i/2-1) + n) + j), xm(j) apoyo(i)]; | ||
| + | fprintf(Verticesdoc,formatSpec,nodo); | ||
| + | end | ||
| + | end | ||
| + | end | ||
| + | fclose(Verticesdoc); | ||
| + | |||
| + | Elementosdoc = fopen(E{p},'w'); | ||
| + | formatSpec = '%d %d %d %d \n'; | ||
| + | for i = 1:m | ||
| + | for j = 1:m | ||
| + | for k = 1:4 | ||
| + | if k == 1 | ||
| + | triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*i+j, (n+m)*(i-1)+j]; | ||
| + | fprintf(Elementosdoc,formatSpec,triang); | ||
| + | elseif k == 2 | ||
| + | triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*(i-1)+j, (n+m)*(i-1)+(j+1)]; | ||
| + | fprintf(Elementosdoc,formatSpec,triang); | ||
| + | elseif k == 3 | ||
| + | triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*(i-1)+(j+1), (n+m)*i+(j+1)]; | ||
| + | fprintf(Elementosdoc,formatSpec,triang); | ||
| + | else | ||
| + | triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*i+(j+1), (n+m)*i+j]; | ||
| + | fprintf(Elementosdoc,formatSpec,triang); | ||
| + | end | ||
| + | end | ||
| + | end | ||
| + | end | ||
| + | fclose(Elementosdoc); | ||
| + | |||
| + | Dirichletdoc = fopen(D{p},'w'); | ||
| + | formatSpec = '%d \n'; | ||
| + | for i = 1:2 | ||
| + | for j = 1:n | ||
| + | if i == 1 | ||
| + | fprintf(Dirichletdoc,formatSpec,(n+m)*(j-1)+1); | ||
| + | else | ||
| + | fprintf(Dirichletdoc,formatSpec,(n+m)*(j-1)+n); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | end | ||
| + | fclose(Dirichletdoc); | ||
| + | |||
| + | if p == 1 | ||
| + | Relaciondoc = fopen(R{p},'w'); | ||
| + | formatSpec = '%d %d \n'; | ||
| + | for i = 1:l | ||
| + | if mod(i,2) == 1 | ||
| + | for j = 1:n | ||
| + | c = 3*i+2*(floor(i/2)-1); | ||
| + | nodo = [(((n+m)*(floor(i/2))) + j), (n2+m2)*(c-1)/2+4*(j-1)+1]; | ||
| + | fprintf(Relaciondoc,formatSpec,nodo); | ||
| + | end | ||
| + | else | ||
| + | for j = 1:m | ||
| + | c = 4*i-3; | ||
| + | nodo = [(((n+m)*(i/2-1) + n) + j), (n2+m2)*((c-1)/2)+4*(j-1)+3]; | ||
| + | fprintf(Relaciondoc,formatSpec,nodo); | ||
| + | end | ||
| + | end | ||
| + | end | ||
| + | fclose(Relaciondoc); | ||
| + | elseif p == 2 | ||
| + | Relaciondoc = fopen(R{p},'w'); | ||
| + | formatSpec = '%d %d \n'; | ||
| + | for i = 1:l | ||
| + | c = 2*i-1; | ||
| + | if mod(i,2) == 1 | ||
| + | for j = 1:n | ||
| + | nodo = [(((n+m)*(floor(i/2))) + j), (n2+m2)*(c-1)/2+2*(j-1)+1]; | ||
| + | fprintf(Relaciondoc,formatSpec,nodo); | ||
| + | end | ||
| + | else | ||
| + | for j = 1:m | ||
| + | nodo = [(((n+m)*(i/2-1) + n) + j), (n2+m2)*((c-1)/2)+2*(j-1)+2]; | ||
| + | fprintf(Relaciondoc,formatSpec,nodo); | ||
| + | end | ||
| + | end | ||
| + | end | ||
| + | fclose(Relaciondoc); | ||
| + | end | ||
| + | |||
| + | Vertices = load(V{p}); | ||
| + | Elementos = load(E{p}); | ||
| + | Dir = load(D{p}); | ||
| + | M = size(Vertices,1); | ||
| + | Nel = size(Elementos,1); | ||
| + | coord = zeros(M,2); | ||
| + | nodos = zeros(Nel,3); | ||
| + | for i = 1:M | ||
| + | for j = 1:2 | ||
| + | coord(i,j) = Vertices(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | for i = 1:Nel | ||
| + | for j = 1:3 | ||
| + | nodos(i,j) = Elementos(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | figure(p) | ||
| + | xlim([min(coord(:,1)),max(coord(:,1))]); | ||
| + | ylim([min(coord(:,2)),max(coord(:,2))]); | ||
| + | grid on | ||
| + | plot(coord(Dir,1),coord(Dir,2),'-b','LineWidth',5); | ||
| + | hold on | ||
| + | for i=1:Nel | ||
| + | nod1 = nodos(i,1); | ||
| + | nod2 = nodos(i,2); | ||
| + | nod3 = nodos(i,3); | ||
| + | x =[coord(nod1,1), coord(nod2,1), coord(nod3,1), coord(nod1,1)]; | ||
| + | y =[coord(nod1,2), coord(nod2,2), coord(nod3,2), coord(nod1,2)]; | ||
| + | fill(x, y, 'cyan') | ||
| + | end | ||
| + | title(T{p}) | ||
| + | end | ||
| + | |||
| + | for p = 1:2 | ||
| + | Vert = load(V{p}); | ||
| + | Elem = load(E{p}); | ||
| + | Vertices = load(V{3}); | ||
| + | Elementos = load(E{3}); | ||
| + | Rel = load(R{p}); | ||
| + | M = size(Vertices,1); | ||
| + | Nel = size(Elementos,1); | ||
| + | coord = zeros(M,2); | ||
| + | nodos = zeros(Nel,3); | ||
| + | coord2 = zeros(M,2); | ||
| + | nodos2 = zeros(Nel,3); | ||
| + | for i = 1:M | ||
| + | for j = 1:2 | ||
| + | coord(i,j) = Vertices(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | for i = 1:Nel | ||
| + | for j = 1:3 | ||
| + | nodos(i,j) = Elementos(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | for i = 1:size(Vert,1) | ||
| + | for j = 1:2 | ||
| + | coord2(i,j) = Vert(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | for i = 1:size(Elem,1) | ||
| + | for j = 1:3 | ||
| + | nodos2(i,j) = Elem(i,j+1); | ||
| + | end | ||
| + | end | ||
| + | |||
| + | figure(p+3) | ||
| + | xlim([min(coord(:,1)),max(coord(:,1))]); | ||
| + | ylim([min(coord(:,2)),max(coord(:,2))]); | ||
| + | grid on | ||
| + | plot(coord(Rel(:,2),1),coord(Rel(:,2),2),'bx','MarkerSize',20); | ||
| + | hold on | ||
| + | plot(coord2(Rel(:,1),1),coord2(Rel(:,1),2),'r.','MarkerSize',20); | ||
| + | hold off | ||
| + | end | ||
[[Categoría:MNEDP|MNEDP]] | [[Categoría:MNEDP|MNEDP]] | ||
[[Categoría:MNEDP25/26|2025-26]] | [[Categoría:MNEDP25/26|2025-26]] | ||
Revisión del 21:07 16 nov 2025
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Flujo estacionario en medio poroso (RoYa) |
| Asignatura | MNEDP |
| Curso | 2025-26 |
| Autores | Rocío Tajuelo Díaz, Yan Wang |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
1 Póster del Trabajo
2 Códigos de MatLab
a = 0; b = 1;
H = [0.1 0.05 0.025]; % Pasos considerados
V = {'Vertices_malla1_tipo1.txt', 'Vertices_malla2_tipo1.txt', 'Vertices_malla3_tipo1.txt'}; E = {'Elementos_malla1_tipo1.txt', 'Elementos_malla2_tipo1.txt', 'Elementos_malla3_tipo1.txt'}; D = {'Dirichlet_malla1_tipo1.txt', 'Dirichlet_malla2_tipo1.txt', 'Dirichlet_malla3_tipo1.txt'}; R = {'Relacion_malla1malla3.txt', 'Relacion_malla2malla3.txt'}; T = {'Mallado con h=0.1', 'Mallado con h=0.05', 'Mallado con h=0.025'};
% Datos fijos necesarios dentro del bucle h = H(3); X = a:h:b; Xm = a+(h/2):h:b; n2 = length(X); m2 = length(Xm); l2 = n2+m2;
for p = 1:3
h = H(p);
x = a:h:b;
xm = a+(h/2):h:b;
n = length(x);
m = length(xm);
apoyo = zeros(n+m,1);
for i = 1:n
apoyo(2*(i-1)+1) = x(i);
if i <= m
apoyo(2*i) = xm(i);
end
end
l = length(apoyo);
Verticesdoc = fopen(V{p},'w');
formatSpec = '%d %.4f %.4f \n';
for i = 1:l
if mod(i,2) == 1
for j = 1:n
nodo = [(((n+m)*(floor(i/2))) + j), x(j), apoyo(i)];
fprintf(Verticesdoc,formatSpec,nodo);
end
else
for j = 1:m
nodo = [(((n+m)*(i/2-1) + n) + j), xm(j) apoyo(i)];
fprintf(Verticesdoc,formatSpec,nodo);
end
end
end
fclose(Verticesdoc);
Elementosdoc = fopen(E{p},'w');
formatSpec = '%d %d %d %d \n';
for i = 1:m
for j = 1:m
for k = 1:4
if k == 1
triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*i+j, (n+m)*(i-1)+j];
fprintf(Elementosdoc,formatSpec,triang);
elseif k == 2
triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*(i-1)+j, (n+m)*(i-1)+(j+1)];
fprintf(Elementosdoc,formatSpec,triang);
elseif k == 3
triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*(i-1)+(j+1), (n+m)*i+(j+1)];
fprintf(Elementosdoc,formatSpec,triang);
else
triang = [(4*(j-1)+k)+40*(i-1), ((n+m)*(i-1)+n)+j, (n+m)*i+(j+1), (n+m)*i+j];
fprintf(Elementosdoc,formatSpec,triang);
end
end
end
end
fclose(Elementosdoc);
Dirichletdoc = fopen(D{p},'w');
formatSpec = '%d \n';
for i = 1:2
for j = 1:n
if i == 1
fprintf(Dirichletdoc,formatSpec,(n+m)*(j-1)+1);
else
fprintf(Dirichletdoc,formatSpec,(n+m)*(j-1)+n);
end
end
end fclose(Dirichletdoc);
if p == 1
Relaciondoc = fopen(R{p},'w');
formatSpec = '%d %d \n';
for i = 1:l
if mod(i,2) == 1
for j = 1:n
c = 3*i+2*(floor(i/2)-1);
nodo = [(((n+m)*(floor(i/2))) + j), (n2+m2)*(c-1)/2+4*(j-1)+1];
fprintf(Relaciondoc,formatSpec,nodo);
end
else
for j = 1:m
c = 4*i-3;
nodo = [(((n+m)*(i/2-1) + n) + j), (n2+m2)*((c-1)/2)+4*(j-1)+3];
fprintf(Relaciondoc,formatSpec,nodo);
end
end
end
fclose(Relaciondoc);
elseif p == 2
Relaciondoc = fopen(R{p},'w');
formatSpec = '%d %d \n';
for i = 1:l
c = 2*i-1;
if mod(i,2) == 1
for j = 1:n
nodo = [(((n+m)*(floor(i/2))) + j), (n2+m2)*(c-1)/2+2*(j-1)+1];
fprintf(Relaciondoc,formatSpec,nodo);
end
else
for j = 1:m
nodo = [(((n+m)*(i/2-1) + n) + j), (n2+m2)*((c-1)/2)+2*(j-1)+2];
fprintf(Relaciondoc,formatSpec,nodo);
end
end
end
fclose(Relaciondoc);
end
Vertices = load(V{p});
Elementos = load(E{p});
Dir = load(D{p});
M = size(Vertices,1);
Nel = size(Elementos,1);
coord = zeros(M,2);
nodos = zeros(Nel,3);
for i = 1:M
for j = 1:2
coord(i,j) = Vertices(i,j+1);
end
end
for i = 1:Nel
for j = 1:3
nodos(i,j) = Elementos(i,j+1);
end
end
figure(p)
xlim([min(coord(:,1)),max(coord(:,1))]);
ylim([min(coord(:,2)),max(coord(:,2))]);
grid on
plot(coord(Dir,1),coord(Dir,2),'-b','LineWidth',5);
hold on
for i=1:Nel
nod1 = nodos(i,1);
nod2 = nodos(i,2);
nod3 = nodos(i,3);
x =[coord(nod1,1), coord(nod2,1), coord(nod3,1), coord(nod1,1)];
y =[coord(nod1,2), coord(nod2,2), coord(nod3,2), coord(nod1,2)];
fill(x, y, 'cyan')
end
title(T{p})
end
for p = 1:2
Vert = load(V{p});
Elem = load(E{p});
Vertices = load(V{3});
Elementos = load(E{3});
Rel = load(R{p});
M = size(Vertices,1);
Nel = size(Elementos,1);
coord = zeros(M,2);
nodos = zeros(Nel,3);
coord2 = zeros(M,2);
nodos2 = zeros(Nel,3);
for i = 1:M
for j = 1:2
coord(i,j) = Vertices(i,j+1);
end
end
for i = 1:Nel
for j = 1:3
nodos(i,j) = Elementos(i,j+1);
end
end
for i = 1:size(Vert,1)
for j = 1:2
coord2(i,j) = Vert(i,j+1);
end
end
for i = 1:size(Elem,1)
for j = 1:3
nodos2(i,j) = Elem(i,j+1);
end
end
figure(p+3)
xlim([min(coord(:,1)),max(coord(:,1))]);
ylim([min(coord(:,2)),max(coord(:,2))]);
grid on
plot(coord(Rel(:,2),1),coord(Rel(:,2),2),'bx','MarkerSize',20);
hold on
plot(coord2(Rel(:,1),1),coord2(Rel(:,1),2),'r.','MarkerSize',20);
hold off
end