Flujo estacionario en medio poroso (RoYa)
De MateWiki
| 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 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:
% 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.
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.
% 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.
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('
Por último, se muestra de nuevo el código pero modificado para el nuevo dominio.
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)'})
