Ver el código de «Flujo estacionario en medio poroso (RoYa)»

{{ TrabajoED | Flujo estacionario en medio poroso (RoYa) | [[:Categoría:MNEDP|MNEDP]]|[[:Categoría:MNEDP25/26|2025-26]] | Rocío Tajuelo Díaz, Yan Wang}}

=Póster del Trabajo=

{|style="margin: 0 auto;"
| [[Archivo:Poster_RoYa2.png|900px|thumb|right|Trabajo de Métodos Numéricos en Ecuaciones Diferenciales en Derivadas Parciales sobre el flujo estacionario en medio poroso.]]
|}

=Códigos de MatLab=

A continuación, se presentan los códigos empleados para la realización del póster.

Para empezar se muestra los código utilizados para generar los documentos necesarios para las diferentes mallas. Para la malla del cuadrado de 4 triángulos se emplea el siguiente código:

{{matlab|codigo=
% Lado del cuadrado inicial.
a = 0;
b = 1;

H = [0.1 0.05 0.025]; % Pasos considerados.

% Nombres de los documentos creados

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'}; % Se emplea para los títulos de las gráficas

% 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;

% Bucle para crear las diferente mallas para cada paso.

for p = 1:3
    h = H(p);
    x = a:h:b; % División del lado en trozos de tamaño h
    xm = a+(h/2):h:b; % Coordenada x de los centroides de los cuadrados en los que se divide el cuadrado inicial.
    n = length(x);
    m = length(xm);
    apoyo = zeros(n+m,1); % Vector con las coordenadas x de los nodos que forman los cuadrados y su centroide .
    for i = 1:n
        apoyo(2*(i-1)+1) = x(i);
        if i <= m
            apoyo(2*i) = xm(i);
        end
    end
    l = length(apoyo);

    % Esta parte crea el documento con la numeración de los nodos y sus coordenada x e y.

    Verticesdoc = fopen(V{p},'w');
    formatSpec = '%d %.4f  %.4f \n';
    for i = 1:l % La numeración de los nodos se basa en el número de nodos que ya se han contado dependiendo del nodo en el que me hallo.
        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);
    
    % Este trozo crea el documento con la numeración de los elementos y los nodos que forman cada elemento.

    Elementosdoc = fopen(E{p},'w');
    formatSpec = '%d %d  %d %d \n';
    for i = 1:m
        for j = 1:m
            for k = 1:4 % El primer elemento es el triángulo de la izquierda, el segundo el de abajo, el tercero el de la derecha y el cuarto el de arriba. La numeración 
                if k == 1  % de los elementos se basa en el número de cuadrados contados dependiendo del nodo del centroide del cuadrado en el que me hallo.
                    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);

    % Esta parte crea el documento que contiene la numeración de los nodos que forman parte de la frontera Dirichlet.

    Dirichletdoc = fopen(D{p},'w');
    formatSpec = '%d \n';
    for i = 1:2 % Los nodos que forman parte de la frontera Dirichlet son aquellos cuya coordenada x es igual a 0 o 1.
        for j = 1:n % Luego, el número de nodos ya contados es proporcional.
            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);

    % Este trozo de código crea el documento que relaciona los nodos de las mallas menos finas con la malla más fija, es decir, las mallas con paso h igual a 0.1 y 0.05 con la malla con paso h igual a 0.025.

    if p == 1 % Malla con paso h igual a 0.1 con la malla más fina.
        Relaciondoc = fopen(R{p},'w');
        formatSpec = '%d %d \n';
        for i = 1:l 
            if mod(i,2) == 1 % La relación entre los primeros nodos de cada fila impar es proporcional al número de nodos que se han contado entre la 1ª y 3ª fila.
                for j = 1:n % En cada fila, la relación entre los nodos es un desplazamiento de 4 posiciones.
                    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 % La relación entre los primeros nodos de cada fila impar es proporcional al número de nodos que se han contado entre la 2ª y 4ª fila más los de la 1ª fila.
                for j = 1:m % En cada fila, la relación entre los nodos es un desplazamiento de 4 posiciones empezando por un pequeño desplazamiento.
                    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  % Malla con paso h igual a 0.05 con la malla más fina.
        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 % En cada fila, la relación entre los nodos es un desplazamiento de 2 posiciones.
                    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 % En cada fila, la relación entre los nodos es un desplazamiento de 4 posiciones empezando por un pequeño desplazamiento.
                    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
    
    % Este parte representa gráficamente los diferentes mallados para cada paso.

    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

% Este trozo dibuja gráficamente los nodos de las mallas menos finas y los de la malla más fina. Esta parte sirve para comprobar que la relación entre los nodos está bien escrita. O sea, si el nodo i de la malla menos fina se corresponde al nodo j de la malla más fina.

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
}}

El siguiente código realiza la malla del cuadrado de 2 triángulos. Es el mismo código que el anterior pero modificado ligeramente. La idea para la numeración es la misma pero modifica a esta forma de mallado.

{{matlab|codigo=

a = 0;
b = 1;

H = [0.1 0.05 0.025];

V = {'Vertices_malla1_tipo2.txt', 'Vertices_malla2_tipo2.txt', 'Vertices_malla3_tipo2.txt'};
E = {'Elementos_malla1_tipo2.txt', 'Elementos_malla2_tipo2.txt', 'Elementos_malla3_tipo2.txt'};
D = {'Dirichlet_malla1_tipo2.txt', 'Dirichlet_malla2_tipo2.txt', 'Dirichlet_malla3_tipo2.txt'};
R = {'Relacion_malla1malla3_tipo1.txt', 'Relacion_malla2malla3_tipo1.txt','Relacion_malla3malla3_tipo1.txt'};

T = {'Mallado con h=0.1', 'Mallado con h=0.05', 'Mallado con h=0.025'};

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;
    n = length(x);
    
    % Crea el documento para los nodos.

    Verticesdoc = fopen(V{p},'w');
    formatSpec = '%d %.4f  %.4f \n';
    for j = 1:n
        for i = 1:n
            nodo = [(j-1)*n + i, x(i), x(j)];
            fprintf(Verticesdoc,formatSpec,nodo);
        end
    end
    fclose(Verticesdoc);
    
    % Crea el documento para los elementos.

    Elementosdoc = fopen(E{p},'w');
    formatSpec = '%d %d  %d %d \n';
    for i = 1:n-1
        for j = 1:n-1
            triang = [(2*(j-1)+1)+(n-1)*(i-1), n*(i-1)+j, n*(i-1)+j+1, n*i+j];
            fprintf(Elementosdoc,formatSpec,triang);
            triang = [(2*(j-1)+2)+(n-1)*(i-1), n*i+j+1, n*i+j, n*(i-1)+j+1];
            fprintf(Elementosdoc,formatSpec,triang);
        end
    end
    fclose(Elementosdoc);

   % Crea el documento para los nodos de la frontera Dirichlet.

    Dirichletdoc = fopen(D{p},'w');
    formatSpec = '%d \n';
    for i = 1:2
        for j = 1:n
            if i == 1
                fprintf(Dirichletdoc,formatSpec,1+n*(j-1));
            else
                fprintf(Dirichletdoc,formatSpec,n+n*(j-1));
            end
        end

    end
    fclose(Dirichletdoc);

    % Crea el documento para la relación entre las mallas de 2 triángulos con la malla de 4 triángulos con paso h igual a 0.025.

    if p == 1
        Relaciondoc = fopen(R{p},'w');
        formatSpec = '%d %d \n';
        for i = 1:n
            for j = 1:n
                c = 8*(i-1)+1;
                nodo = [(i-1)*n + j, (n2+m2)*(c-1)/2+4*(j-1)+1];
                fprintf(Relaciondoc,formatSpec,nodo);
            end
        end
        fclose(Relaciondoc);
    elseif p == 2
        Relaciondoc = fopen(R{p},'w');
        formatSpec = '%d %d \n';
        for i = 1:n
            c = 4*(i-1)+1;
            for j = 1:n
                nodo = [(i-1)*n + j, (n2+m2)*(c-1)/2+2*(j-1)+1];
                fprintf(Relaciondoc,formatSpec,nodo);
            end
        end
        fclose(Relaciondoc);
    else
        Relaciondoc = fopen(R{p},'w');
        formatSpec = '%d %d \n';
        for i = 1:n
            c = 2*(i-1)+1;
            for j = 1:n
                nodo = [(i-1)*n + j, (n2+m2)*(c-1)/2+j];
                fprintf(Relaciondoc,formatSpec,nodo);
            end
        end
        fclose(Relaciondoc);
    end
    
   % Representa las diferentes mallas.

    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)
    plot(coord(Dir,1),coord(Dir,2),'-b','LineWidth',5);
    xlim([min(coord(:,1)),max(coord(:,1))]);
    ylim([min(coord(:,2)),max(coord(:,2))]);
    grid on
    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

% Representa los nodos de la relación para la comprobación.

for p = 1:3
    Vertices = load(V{p});
    Elememtos = load(E{p});
    Vert = load('Vertices_malla3_tipo1.txt');
    Elem = load('Elementos_malla3_tipo1.txt');
    Rel = load(R{p});
    M = size(Vertices,1);
    Nel = size(Elementos,1);
    coord = zeros(M,2);
    nodos = zeros(Nel,2);
    coord2 = zeros(size(Vert,1),2);
    nodos2 = zeros(size(Elem,1),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+4)
    xlim([min(coord(:,1)),max(coord(:,1))]);
    ylim([min(coord(:,2)),max(coord(:,2))]);
    grid on
    plot(coord(Rel(:,1),1),coord(Rel(:,1),2),'bx','MarkerSize',20);
    hold on
    plot(coord2(Rel(:,2),1),coord2(Rel(:,2),2),'r.','MarkerSize',20);
    hold off
end
}}

Lo siguientes códigos realizan el MEF. Primero se muestra con las mallas de 4 triángulos.

{{matlab|codigo=

% Documentos donde se encuentra la información de las mallas.

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'};

% Lado del cuadrado inicial.

a = 0;
b = 1;

% Condiciones de tipo Dirichlet.

p0 = 0.5;
p1 = 0;

% Datos del problema.

f = @(x,y) 50*exp(-100*((x-0.5).^2 + (y-0.5).^2));
K = @(x,y) 1 + 0.5.*cos(2*pi*x).*cos(2*pi*y);

% Matrices empleadas para calcular la matriz de rigidez por ensamblado.

Lxx = zeros(3,3); Lxx(1,1) = 1/2; Lxx(2,2) = 1/2; Lxx(1,2) = -1/2; Lxx(2,1) = -1/2;
Lyy = zeros(3,3); Lyy(1,1) = 1/2; Lyy(3,3) = 1/2; Lyy(3,1) = -1/2; Lyy(1,3) = -1/2;
Lxy = zeros(3,3); Lxy(1,1) = 1/2; Lxy(2,3) = 1/2; Lxy(1,3) = -1/2; Lxy(2,1) = -1/2;
l0 = (1/6)*ones(3,1);

% Variables de la tabla final sobre la convergencia.

H = [0.1 0.05 0.025];
Sol = zeros(length(a:H(end):b),3);
normL = zeros(size(H));
ordenL = zeros(size(H));
normH = zeros(size(H));
ordenH = zeros(size(H));
condA = zeros(size(H));

% Resolución del problema para las mallas de cada paso.

for p =1:3
    
    % Lectura de los documentos y creación de variables con los datos de estos.

    Dir = load(D{p});
    Vertices = load(V{p});
    Elementos = load(E{p});
    MDir = length(Dir);
    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
    
    % Matrices del problema en forma matricial.

    A = zeros(M,M);
    b = zeros(M,1);
    u = zeros(M,1);
    
    % Asignación valores a la solución de los nodos de la frontera Dirichlet.

    for i = 1:length(Dir)
        if coord(Dir(i),1) == a
            u(Dir(i)) = p0;
        else
            u(Dir(i)) = p1;
        end
    end
    
    ind = setdiff(1:M,sort(Dir)); % Nodos cuya solución hay que calcular
    
    % Ensamblado de la matriz de rigidez y el vector de cargas.

    for e = 1:Nel
        nodosK = Elementos(e, 2:end);
        Bk = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        Ck = invBk*invBk';
        fk = 0;
        Kk = 0;
        for i = 1:3 % Aproximación de la integral de f y K sobre cada elemento como un promedio.
            Kk = Kk + K(coord(nodosK(i),1),coord(nodosK(i),2));
            fk = fk + f(coord(nodosK(i),1),coord(nodosK(i),2));
        end
        Kk = Kk./3;
        fk = fk./3;
        Lk = abs(detBk)*Kk*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy);
        A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
        lk = fk*abs(detBk)*l0;
        b(nodosK) = b(nodosK) + lk;
    end
    
    for i = 1:M-MDir % Ajuste del vector de cargas.
        camb = 0;
        for j = 1:MDir
            camb = camb + A(Dir(j),ind(i))*u(Dir(j));
        end
        b(ind(i)) = b(ind(i)) - camb;
    end
    
    u(ind) = A(ind,ind)\b(ind); % Resolución del problema matricial.

    Sol(1:length(u),p) = u;
    
    % Representación de la solución calculada.

    figure(p)
    xlim([min(coord(:,1)) max(coord(:,1))]);
    ylim([min(coord(:,2)) max(coord(:,2))]);
    grid on
    trisurf(nodos,coord(:,1),coord(:,2),u);view(0,90);
    colorbar;
    title('u_h');
    
end

% Cálculo de la convergencia.

for p = 1:2
    Rel = load(R{p});
    Vertices = load(V{p});
    Elementos = load(E{p});
    Vert = load(V{end});
    M = size(Vertices,1);
    M3 = size(Vert,1);
    Nel = size(Elementos,1);
    coord = zeros(M,2);
    coord3 = zeros(M3,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:M3
        for j = 1:2
            coord3(i,j) = Vert(i,j+1);
        end
    end
    
    for i = 1:Nel
        for j = 1:3
            nodos(i,j) = Elementos(i,j+1);
        end
    end

    % Cálculo de las normas por aproximación de las soluciones de las mallas menos finas con la más fina.

    intL = 0;
    intdH = 0;
    for i = 1:Nel
        nodosK = nodos(i,:);
        Bk = zeros(2,2);
        Bk3 = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
                Bk3(j,k) = coord3(Rel(nodosK(k+1),2),j)-coord3(Rel(nodosK(1),2),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        detBk3 = det(Bk3);
        invBk3 = inv(Bk3);
        sumL = 0;
        sumdH = 0;
        uk = [Sol(nodosK(2),p)-Sol(nodosK(1),p); Sol(nodosK(3),p)-Sol(nodosK(1),p)];
        uk3 = [Sol(Rel(nodosK(2),2),3)-Sol(Rel(nodosK(1),2),3); Sol(Rel(nodosK(3),2),3)-Sol(Rel(nodosK(1),2),3)];
        for j = 1:3
            sumL = sumL + (Sol(Rel(nodosK(j),2),3) - Sol(nodosK(j),p)).^2;
            sumdH = sumdH + norm(invBk3'*uk3 - invBk'*uk).^2;
        end
        intL = intL + (abs(detBk)/3)*sumL;
        intdH = intdH + (abs(detBk)/3)*sumdH;
    end
    
    normL(p) = sqrt(intL);
    normH(p) = sqrt(intL + intdH);
    condA(p) = cond(A);
end

% Cálculo del orden en las diferentes normas.

ordenL(1) = log(normL(1)/normL(2))/log(2);
ordenH(1) = log(normH(1)/normH(2))/log(2);

% Creación de la tabla de convergencia con los datos calculados.

table(H',normL',ordenL',normH',ordenH',condA','VariableNames',{'h','norma en L2','orden L2','norma en H1','orden H1','cond(A)'})

}}

Ahora, se muestra el mismo código pero ajustado con la parte de convergencia ajustada para comparar con la malla de 4 triángulos más fina.

{{matlab|codigo=
V = {'Vertices_malla1_tipo2.txt', 'Vertices_malla2_tipo2.txt', 'Vertices_malla3_tipo2.txt', 'Vertices_malla3_tipo1.txt'};
E = {'Elementos_malla1_tipo2.txt', 'Elementos_malla2_tipo2.txt', 'Elementos_malla3_tipo2.txt', 'Elementos_malla3_tipo1.txt'};
D = {'Dirichlet_malla1_tipo2.txt', 'Dirichlet_malla2_tipo2.txt', 'Dirichlet_malla3_tipo2.txt', 'Dirichlet_malla3_tipo1.txt'};
R = {'Relacion_malla1malla3_tipo1.txt', 'Relacion_malla2malla3_tipo1.txt', 'Relacion_malla3malla3_tipo1.txt'};

a = 0;
b = 1;

p0 = 0.5;
p1 = 0;

f = @(x,y) 50*exp(-100*((x-0.5).^2 + (y-0.5).^2));
K = @(x,y) 1 + 0.5.*cos(2*pi*x).*cos(2*pi*y);


Lxx = zeros(3,3); Lxx(1,1) = 1/2; Lxx(2,2) = 1/2; Lxx(1,2) = -1/2; Lxx(2,1) = -1/2;
Lyy = zeros(3,3); Lyy(1,1) = 1/2; Lyy(3,3) = 1/2; Lyy(3,1) = -1/2; Lyy(1,3) = -1/2;
Lxy = zeros(3,3); Lxy(1,1) = 1/2; Lxy(2,3) = 1/2; Lxy(1,3) = -1/2; Lxy(2,1) = -1/2;
l0 = (1/6)*ones(3,1);

H = [0.1 0.05 0.025];
Sol = zeros(length(a:H(end):b),3);
Exacta = zeros(length(a:H(end):b),1);
coord5 = zeros(length(a:H(end):b),1);
normL = zeros(size(H));
ordenL = zeros(size(H));
normH = zeros(size(H));
ordenH = zeros(size(H));
condA = zeros(size(H));

for p =1:4

    % Lectura de documentos y creación de variables con dichos datos.

    Dir = load(D{p});
    Vertices = load(V{p});
    Elementos = load(E{p});
    MDir = length(Dir);
    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
    
    if p ==3 % Almacenamos las coordenadas de la malla de 2 triángulos más fina para uso posterior.
        coord5 = coord;
    end
    
    A = zeros(M,M);
    b = zeros(M,1);
    u = zeros(M,1);
    
    % Asignación de variables con valores conocidos

    for i = 1:length(Dir)
        if coord(Dir(i),1) == a
            u(Dir(i)) = p0;
        else
            u(Dir(i)) = p1;
        end
    end
    
    ind = setdiff(1:M,sort(Dir));
    
    % Ensamblado de la matriz de rigidez y del vector de cargas.

    for e = 1:Nel
        nodosK = Elementos(e, 2:end);
        Bk = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        Ck = invBk*invBk';
        fk = 0;
        Kk = 0;
        for i = 1:3
            Kk = Kk + K(coord(nodosK(i),1),coord(nodosK(i),2));
            fk = fk + f(coord(nodosK(i),1),coord(nodosK(i),2));
        end
        Kk = Kk./3;
        fk = fk./3;
        Lk = abs(detBk)*Kk*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy);
        A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
        lk = fk*abs(detBk)*l0;
        b(nodosK) = b(nodosK) + lk;
    end
    
    % Ajuste del vector de cargas

    for i = 1:M-MDir
        camb = 0;
        for j = 1:MDir
            camb = camb + A(Dir(j),ind(i))*u(Dir(j));
        end
        b(ind(i)) = b(ind(i)) - camb;
    end
    
    u(ind) = A(ind,ind)\b(ind);

    if p == 4  % Almacenamos la solución de la malla de 4 triángulos en una variable diferente.
        Exacta(1:length(u)) = u;
    else
        Sol(1:length(u),p) = u;
    end
    
    % Representación gráfica de la solución calculada.

    figure(p)
    xlim([min(coord(:,1)) max(coord(:,1))]);
    ylim([min(coord(:,2)) max(coord(:,2))]);
    grid on
    trisurf(nodos,coord(:,1),coord(:,2),u);view(0,90);
    colorbar;
    title('u_h');
end

% Representación gráfica de la diferencia entre la solución de la malla de 2 triángulos más fina con la solución de la malla de 4 triángulos más fina.

Rel = load(R{3});
Elementos= load(E{3});
figure(p+1)
xlim([min(coord(:,1)) max(coord(:,1))]);
ylim([min(coord(:,2)) max(coord(:,2))]);
grid on
trisurf(Elementos(:,2:end),coord5(:,1),coord5(:,2),abs(Sol(:,3)-Exacta(Rel(:,2))));view(0,90);
colorbar;
title('|u-u_h|');

% Creación de variables necesarias para la convergencia.

Vert = load(V{4});
Elem = load(E{4});
M4 = size(Vert,1);
Nel4 = size(Elem,1);
coord4 = zeros(M4,2);
nodos4 = zeros(Nel4,3);
for i = 1:M4
    for j = 1:2
        coord4(i,j) = Vert(i,j+1);
    end
end

for i = 1:Nel4
    for j = 1:3
        nodos4(i,j) = Elem(i,j+1);
    end
end

% Cálculo de la convergencia.

for p = 1:3
    Rel = load(R{p});
    Vertices = load(V{p});
    Elementos = load(E{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

    % Cálculo de las normas por aproximación de las soluciones de las mallas de 2 triángulos con la solución de la malla de 4 triángulos más fina.

    intL = 0;
    intdH = 0;
    for i = 1:Nel
        nodosK = nodos(i,:);
        Bk = zeros(2,2);
        Bk4 = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
                Bk4(j,k) = coord4(Rel(nodosK(k+1),2),j)-coord4(Rel(nodosK(1),2),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        detBk4 = det(Bk4);
        invBk4 = inv(Bk4);
        sumL = 0;
        sumdH = 0;
        uk = [Sol(nodosK(2),p)-Sol(nodosK(1),p); Sol(nodosK(3),p)-Sol(nodosK(1),p)];
        uk4 = [Exacta(Rel(nodosK(2),2))-Exacta(Rel(nodosK(1),2)); Exacta(Rel(nodosK(3),2))-Exacta(Rel(nodosK(1),2))];
        for j = 1:3
            sumL = sumL + (Exacta(Rel(nodosK(j),2)) - Sol(nodosK(j),p)).^2;
            sumdH = sumdH + norm(invBk4'*uk4 - invBk'*uk).^2;
        end
        intL = intL + (abs(detBk)/3)*sumL;
        intdH = intdH + (abs(detBk)/3)*sumdH;
    end
    
    normL(p) = sqrt(intL);
    normH(p) = sqrt(intL + intdH);
    condA(p) = cond(A);
end

% Cálculo del orden de las diferentes normas.

for i = 1:2
    ordenL(i) = log(normL(i)/normL(i+1))/log(2);
    ordenH(i) = log(normH(i)/normH(i+1))/log(2);
end

% Creación de la tabla de convergencia con los datos calculados.

table(H',normL',ordenL',normH',ordenH',condA','VariableNames',{'h','norma en L2','orden L2','norma en H1','orden H1','cond(A)'})
}}

Por último, se muestra de nuevo el código pero modificado para el nuevo dominio.

{{matlab|codigo=

V = {'Vertices_malla1_tipo3.txt', 'Vertices_malla2_tipo3.txt', 'Vertices_malla3_tipo3.txt'};
E = {'Elementos_malla1_tipo3.txt', 'Elementos_malla2_tipo3.txt', 'Elementos_malla3_tipo3.txt'};
D = {'Dirichlet_malla1_tipo3.txt', 'Dirichlet_malla2_tipo3.txt', 'Dirichlet_malla3_tipo3.txt'};
R = {'Relacion_malla1malla3_tipo3.txt', 'Relacion_malla2malla3_tipo3.txt'};

% Lectura de una de las mallas para el cálculo del centroide del dominio.

Vertices = load(V{3});
Elementos = load(E{3});
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

% Cálculo del centroide del dominio por cambio de variable al elemento de referencia.

area2 = 0;
xcmintcv = 0;
ycmintcv = 0;
for i = 1:Nel
    jac = zeros(2,2);
    triang = [nodos(i,1), nodos(i,2), nodos(i,3)];
    for j = 1:2
        for k = 1:2
            jac(j,k) = coord(triang(k+1),j)-coord(triang(1),j);
        end
    end
    detjac = abs((jac(1,1)*jac(2,2))-(jac(1,2)*jac(2,1)));
    areael = detjac*(1/2);
    area2 = area2 + areael;
    intx = @(a,b) jac(1,1)*a + jac(1,2)*b + coord(triang(1),1);
    inty = @(a,b) jac(2,1)*a + jac(2,2)*b + coord(triang(1),2);
    xcmintcv = xcmintcv + detjac*integral2(intx,0,1,0,@(z) -z+1);
    ycmintcv = ycmintcv + detjac*integral2(inty,0,1,0,@(z) -z+1);
end
xcmcv = xcmintcv/area2;
ycmcv = ycmintcv/area2;

a = 0;
b = 1;

p0 = 0.5;
p1 = 0;

f = @(x,y) 50*exp(-100*((x-xcmcv).^2 + (y-ycmcv).^2));
K = @(x,y) 1 + 0.5.*cos(2*pi*x).*cos(2*pi*y);

Lxx = zeros(3,3); Lxx(1,1) = 1/2; Lxx(2,2) = 1/2; Lxx(1,2) = -1/2; Lxx(2,1) = -1/2;
Lyy = zeros(3,3); Lyy(1,1) = 1/2; Lyy(3,3) = 1/2; Lyy(3,1) = -1/2; Lyy(1,3) = -1/2;
Lxy = zeros(3,3); Lxy(1,1) = 1/2; Lxy(2,3) = 1/2; Lxy(1,3) = -1/2; Lxy(2,1) = -1/2;
l0 = (1/6)*ones(3,1);

H = [0.1 0.05 0.025];
Sol = zeros(length(a:H(end):b),3);
normL = zeros(size(H));
ordenL = zeros(size(H));
normH = zeros(size(H));
ordenH = zeros(size(H));
condA = zeros(size(H));

for p =1:3

    % Lectura de los documentos y creación de variables con dichos datos.

    Dir = load(D{p});
    Vertices = load(V{p});
    Elementos = load(E{p});
    MDir = length(Dir);
    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
    
    A = zeros(M,M);
    b = zeros(M,1);
    u = zeros(M,1);
    
    % Asignación de valores conocidos.

    for i = 1:length(Dir)
        if coord(Dir(i),1) == a
            u(Dir(i)) = p0;
        else
            u(Dir(i)) = p1;
        end
    end
    
    ind = setdiff(1:M,sort(Dir));
    
    % Ensamblado de la matriz de rigidez y del vector de cargas.

    for e = 1:Nel
        nodosK = Elementos(e, 2:end);
        Bk = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        Ck = invBk*invBk';
        fk = 0;
        Kk = 0;
        for i = 1:3
            Kk = Kk + K(coord(nodosK(i),1),coord(nodosK(i),2));
            fk = fk + f(coord(nodosK(i),1),coord(nodosK(i),2));
        end
        Kk = Kk./3;
        fk = fk./3;
        Lk = abs(detBk)*Kk*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy);
        A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
        lk = fk*abs(detBk)*l0;
        b(nodosK) = b(nodosK) + lk;
    end
    
    % Ajuste del vector de cargas.

    for i = 1:M-MDir
        camb = 0;
        for j = 1:MDir
            camb = camb + A(Dir(j),ind(i))*u(Dir(j));
        end
        b(ind(i)) = b(ind(i)) - camb;
    end
    
    u(ind) = A(ind,ind)\b(ind);

    Sol(1:length(u),p) = u;
    
    % Representación gráfica de la solución calculada.

    figure(p)
    xlim([min(coord(:,1)) max(coord(:,1))]);
    ylim([min(coord(:,2)) max(coord(:,2))]);
    grid on
    trisurf(nodos,coord(:,1),coord(:,2),u);view(0,90);
    colorbar;
    title(strcat('u_h con h =',num2str(H(p))));
end

% Cálculo de la convergencia.

for p = 1:2
    Rel = load(R{p});
    Vertices = load(V{p});
    Elementos = load(E{p});
    Vert = load(V{end});
    M = size(Vertices,1);
    M3 = size(Vert,1);
    Nel = size(Elementos,1);
    coord = zeros(M,2);
    coord3 = zeros(M3,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:M3
        for j = 1:2
            coord3(i,j) = Vert(i,j+1);
        end
    end
    
    for i = 1:Nel
        for j = 1:3
            nodos(i,j) = Elementos(i,j+1);
        end
    end

    % Cálculo de las normas por aproximación de las soluciones de las mallas menos finas con la solución de la malla más fina.

    intL = 0;
    intdH = 0;
    for i = 1:Nel
        nodosK = nodos(i,:);
        Bk = zeros(2,2);
        Bk3 = zeros(2,2);
        for j = 1:2
            for k = 1:2
                Bk(j,k) = coord(nodosK(k+1),j)-coord(nodosK(1),j);
                Bk3(j,k) = coord3(Rel(nodosK(k+1),2),j)-coord3(Rel(nodosK(1),2),j);
            end
        end
        detBk = det(Bk);
        invBk = inv(Bk);
        detBk3 = det(Bk3);
        invBk3 = inv(Bk3);
        sumL = 0;
        sumdH = 0;
        uk = [Sol(nodosK(2),p)-Sol(nodosK(1),p); Sol(nodosK(3),p)-Sol(nodosK(1),p)];
        uk3 = [Sol(Rel(nodosK(2),2),3)-Sol(Rel(nodosK(1),2),3); Sol(Rel(nodosK(3),2),3)-Sol(Rel(nodosK(1),2),3)];
        for j = 1:3
            sumL = sumL + (Sol(Rel(nodosK(j),2),3) - Sol(nodosK(j),p)).^2;
            sumdH = sumdH + norm(invBk3'*uk3 - invBk'*uk).^2;
        end
        intL = intL + (abs(detBk)/3)*sumL;
        intdH = intdH + (abs(detBk)/3)*sumdH;
    end
    
    normL(p) = sqrt(intL);
    normH(p) = sqrt(intL + intdH);
    condA(p) = cond(A);
end

% Cálculo del orden de las diferentes normas.

ordenL(1) = log(normL(1)/normL(2))/log(2);
ordenH(1) = log(normH(1)/normH(2))/log(2);

% Creación de la tabla de convergencia con los datos calculados.

table(H',normL',ordenL',normH',ordenH',condA','VariableNames',{'h','norma en L2','orden L2','norma en H1','orden H1','cond(A)'})
}}

[[Categoría:MNEDP|MNEDP]]
[[Categoría:MNEDP25/26|2025-26]]