Diferencia entre revisiones de «Flujo estacionario en medio poroso (RoYa)»

De MateWiki
Saltar a: navegación, buscar
(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