Flujo estacionario en medio poroso (RoYa)

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

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