Diferencia entre revisiones de «Contaminación del lago (AMLS)»

De MateWiki
Saltar a: navegación, buscar
Línea 7: Línea 7:
  
  
 +
=Códigos de Matlab=
 +
<source lang="matlab">
 +
close all
 +
clear all
  
 +
%-------------------------------MALLA 1------------------------------------
  
 +
%%%----------PASO 1: VARIABLES-----------
  
 +
% Término fuente f(x,y)
 +
% --- Fuente física: vertido puntual normalizado (gaussiana) ---
 +
Q = 1.0;          % Caudal total de contaminante
 +
x0 = 0.5; y0 = 0.5; % Punto de vertido del contaminante
 +
sigma = 0.05;    % Anchura de la fuente
  
 +
f = @(x,y) (Q/(2*pi*sigma^2)) .* exp(-((x-x0).^2 + (y-y0).^2) / (2*sigma^2));
  
 +
% Parámetros físicos
 +
D = 0.01;  % Coeficiente de difusión
 +
k = 0.1;  % Coeficiente de reacción
 +
 +
 +
%%%----------PASO 2: CREACIÓN DE LAS MATRICES DE NODOS Y TRIANGULOS-----------
 +
 +
Vertices = load('vertices1.txt');
 +
Triangulos = load('triangulos1.txt');    %Cargamos los ficheros
 +
 +
Ne = length(Triangulos);    % número de Triangulos
 +
M = length(Vertices);      % número de nodos (vertices)
 +
 +
% Coordenadas de los nodos
 +
coord1 = zeros(M,2);
 +
coord1(:,1) = Vertices(:,2);
 +
coord1(:,2) = Vertices(:,3);
 +
 +
% Conectividad (qué nodos forman cada triángulo)
 +
nodos = zeros(Ne,3);
 +
nodos(:,1) = Triangulos(:,2);
 +
nodos(:,2) = Triangulos(:,3);
 +
nodos(:,3) = Triangulos(:,4);
 +
nodos1=nodos;
 +
%%%----------PASO 3: MATRICES LOCALES DE REFERENCIA-----------
 +
 +
% Matriz de masa local (para el término de reacción)
 +
L0 = [1/12 1/24 1/24;
 +
      1/24 1/12 1/24;
 +
      1/24 1/24 1/12];
 +
 +
% Matrices para los términos de derivadas parciales
 +
Lxx = [1/2 -1/2 0;
 +
      -1/2 1/2 0;
 +
        0    0  0];
 +
 +
Lyy = [1/2  0 -1/2;
 +
        0  0  0;
 +
      -1/2 0  1/2];
 +
 +
Lxy = [1/2  0 -1/2;
 +
      -1/2 0  1/2;
 +
        0  0  0];
 +
 +
% Vector de integración de referencia (para f)
 +
l0 = [1/6; 1/6; 1/6];
 +
 +
 +
%%%----------PASO 4: ENSAMBLADO DE LA MATRIZ GLOBAL Y EL VECTOR DE CARGAS-----------
 +
 +
A = zeros(M);    % matriz global (rigidez + reacción)
 +
b = zeros(M,1);  % vector global de cargas
 +
 +
for e = 1:Ne
 +
    % Nodos del triangulo e
 +
    nodosK = nodos(e,:);
 +
   
 +
    % Matriz del elemento (de coordenadas locales a reales)
 +
    Bk = [coord1(nodosK(2),1)-coord1(nodosK(1),1), coord1(nodosK(3),1)-coord1(nodosK(1),1);
 +
          coord1(nodosK(2),2)-coord1(nodosK(1),2), coord1(nodosK(3),2)-coord1(nodosK(1),2)];
 +
   
 +
    detBk = det(Bk);          % área*2 del triángulo
 +
    Ck = inv(Bk) * (inv(Bk))'; % matriz Ck = (Bk^-1)(Bk^-T)
 +
   
 +
    % Matriz elemental (difusión + reacción)
 +
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
 +
   
 +
    % Ensamblado en la matriz global
 +
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
 +
    A1=A;
 +
   
 +
    % Cálculo del término fuente medio en el triángulo
 +
    fk = mean([ ...
 +
        f(coord1(nodosK(1),1), coord1(nodosK(1),2)), ...
 +
        f(coord1(nodosK(2),1), coord1(nodosK(2),2)), ...
 +
        f(coord1(nodosK(3),1), coord1(nodosK(3),2)) ]);
 +
   
 +
    % Vector local de cargas
 +
    lk = fk * abs(detBk) * l0;
 +
   
 +
    % Ensamblado en el vector global
 +
    b(nodosK) = b(nodosK) + lk;
 +
end
 +
 +
%%%----------PASO 5: APLICAR CONDICIONES DE CONTORNO DIRICHLET HOMOGÉNEAS (c=0)-----------
 +
 +
% Detectamos los nodos de frontera:
 +
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
 +
edges_sorted = sort(edges,2);
 +
[uniqE, ~, ic] = unique(edges_sorted,'rows');
 +
counts = accumarray(ic,1);
 +
boundary_edges = uniqE(counts==1,:);
 +
boundary_nodes = unique(boundary_edges(:));
 +
 +
% Impongo c=0 en los nodos de la frontera
 +
for n = boundary_nodes'
 +
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
 +
    b(n) = 0;
 +
end
 +
 +
%%%----------PASO 6: RESOLUCIÓN DEL SISTEMA-----------
 +
 +
u1 = A\b;  % solución aproximada en los vertices
 +
 +
%%%----------PASO 7: GRAFICAMOS-----------
 +
 +
figure(1)
 +
trisurf(nodos, coord1(:,1), coord1(:,2), u1)
 +
title('Solución numérica c_h')
 +
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
 +
 +
 +
 +
%--------------------------------------MALLA 2--------------------------------------------------------
 +
 +
 +
Vertices = load('vertices2.txt');
 +
Triangulos = load('triangulos2.txt');   
 +
 +
Ne = length(Triangulos);   
 +
M = length(Vertices);     
 +
 +
 +
coord2 = zeros(M,2);
 +
coord2(:,1) = Vertices(:,2);
 +
coord2(:,2) = Vertices(:,3);
 +
 +
nodos = zeros(Ne,3);
 +
nodos(:,1) = Triangulos(:,2);
 +
nodos(:,2) = Triangulos(:,3);
 +
nodos(:,3) = Triangulos(:,4);
 +
nodos2=nodos;
 +
 +
L0 = [1/12 1/24 1/24;
 +
      1/24 1/12 1/24;
 +
      1/24 1/24 1/12];
 +
 +
 +
Lxx = [1/2 -1/2 0;
 +
      -1/2 1/2 0;
 +
        0    0  0];
 +
 +
Lyy = [1/2  0 -1/2;
 +
        0  0  0;
 +
      -1/2 0  1/2];
 +
 +
Lxy = [1/2  0 -1/2;
 +
      -1/2 0  1/2;
 +
        0  0  0];
 +
 +
 +
l0 = [1/6; 1/6; 1/6];
 +
 +
A = zeros(M);   
 +
b = zeros(M,1); 
 +
 +
for e = 1:Ne
 +
   
 +
    nodosK = nodos(e,:);
 +
   
 +
   
 +
    Bk = [coord2(nodosK(2),1)-coord2(nodosK(1),1), coord2(nodosK(3),1)-coord2(nodosK(1),1);
 +
          coord2(nodosK(2),2)-coord2(nodosK(1),2), coord2(nodosK(3),2)-coord2(nodosK(1),2)];
 +
   
 +
    detBk = det(Bk);       
 +
    Ck = inv(Bk) * (inv(Bk))';
 +
   
 +
 
 +
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
 +
   
 +
   
 +
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
 +
    A2=A;
 +
   
 +
 
 +
    fk = mean([ ...
 +
        f(coord2(nodosK(1),1), coord2(nodosK(1),2)), ...
 +
        f(coord2(nodosK(2),1), coord2(nodosK(2),2)), ...
 +
        f(coord2(nodosK(3),1), coord2(nodosK(3),2)) ]);
 +
   
 +
 +
    lk = fk * abs(detBk) * l0;
 +
   
 +
 +
    b(nodosK) = b(nodosK) + lk;
 +
end
 +
 +
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
 +
edges_sorted = sort(edges,2);
 +
[uniqE, ~, ic] = unique(edges_sorted,'rows');
 +
counts = accumarray(ic,1);
 +
boundary_edges = uniqE(counts==1,:);
 +
boundary_nodes = unique(boundary_edges(:));
 +
 +
 +
for n = boundary_nodes'
 +
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
 +
    b(n) = 0;
 +
end
 +
 +
 +
u2 = A\b; 
 +
 +
figure(2)
 +
trisurf(nodos, coord2(:,1), coord2(:,2), u2)
 +
title('Solución numérica c_h')
 +
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
 +
 +
 +
%---------------------------------------MALA 3-----------------------------------------------------------
 +
 +
 +
 +
 +
Vertices = load('vertices3.txt');
 +
Triangulos = load('triangulos3.txt');   
 +
 +
Ne = length(Triangulos);   
 +
M = length(Vertices);     
 +
 +
 +
coord3 = zeros(M,2);
 +
coord3(:,1) = Vertices(:,2);
 +
coord3(:,2) = Vertices(:,3);
 +
 +
 +
nodos = zeros(Ne,3);
 +
nodos(:,1) = Triangulos(:,2);
 +
nodos(:,2) = Triangulos(:,3);
 +
nodos(:,3) = Triangulos(:,4);
 +
nodos3=nodos;
 +
 +
L0 = [1/12 1/24 1/24;
 +
      1/24 1/12 1/24;
 +
      1/24 1/24 1/12];
 +
 +
 +
Lxx = [1/2 -1/2 0;
 +
      -1/2 1/2 0;
 +
        0    0  0];
 +
 +
Lyy = [1/2  0 -1/2;
 +
        0  0  0;
 +
      -1/2 0  1/2];
 +
 +
Lxy = [1/2  0 -1/2;
 +
      -1/2 0  1/2;
 +
        0  0  0];
 +
 +
 +
l0 = [1/6; 1/6; 1/6];
 +
 +
A = zeros(M); 
 +
b = zeros(M,1); 
 +
 +
for e = 1:Ne
 +
   
 +
    nodosK = nodos(e,:);
 +
   
 +
   
 +
    Bk = [coord3(nodosK(2),1)-coord3(nodosK(1),1), coord3(nodosK(3),1)-coord3(nodosK(1),1);
 +
          coord3(nodosK(2),2)-coord3(nodosK(1),2), coord3(nodosK(3),2)-coord3(nodosK(1),2)];
 +
   
 +
    detBk = det(Bk);         
 +
    Ck = inv(Bk) * (inv(Bk))';
 +
   
 +
   
 +
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
 +
   
 +
 
 +
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
 +
    A3=A;
 +
   
 +
   
 +
    fk = mean([ ...
 +
        f(coord3(nodosK(1),1), coord3(nodosK(1),2)), ...
 +
        f(coord3(nodosK(2),1), coord3(nodosK(2),2)), ...
 +
        f(coord3(nodosK(3),1), coord3(nodosK(3),2)) ]);
 +
   
 +
   
 +
    lk = fk * abs(detBk) * l0;
 +
   
 +
 
 +
    b(nodosK) = b(nodosK) + lk;
 +
end
 +
 +
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
 +
edges_sorted = sort(edges,2);
 +
[uniqE, ~, ic] = unique(edges_sorted,'rows');
 +
counts = accumarray(ic,1);
 +
boundary_edges = uniqE(counts==1,:);
 +
boundary_nodes = unique(boundary_edges(:));
 +
 +
 +
for n = boundary_nodes'
 +
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
 +
    b(n) = 0;
 +
end
 +
 +
 +
 +
u3 = A\b; 
 +
 +
 +
 +
figure(3)
 +
trisurf(nodos, coord3(:,1), coord3(:,2), u3)
 +
title('Solución numérica c_h')
 +
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
 +
 +
 +
%---------------------------------------MALLA 4-------------------------------------------------------------------
 +
 +
Vertices = load('vertices4.txt');
 +
Triangulos = load('triangulos4.txt');   
 +
 +
Ne = length(Triangulos);   
 +
M = length(Vertices);     
 +
 +
 +
coord4 = zeros(M,2);
 +
coord4(:,1) = Vertices(:,2);
 +
coord4(:,2) = Vertices(:,3);
 +
 +
 +
nodos = zeros(Ne,3);
 +
nodos(:,1) = Triangulos(:,2);
 +
nodos(:,2) = Triangulos(:,3);
 +
nodos(:,3) = Triangulos(:,4);
 +
nodos4=nodos;
 +
 +
 +
 +
L0 = [1/12 1/24 1/24;
 +
      1/24 1/12 1/24;
 +
      1/24 1/24 1/12];
 +
 +
 +
Lxx = [1/2 -1/2 0;
 +
      -1/2 1/2 0;
 +
        0    0  0];
 +
 +
Lyy = [1/2  0 -1/2;
 +
        0  0  0;
 +
      -1/2 0  1/2];
 +
 +
Lxy = [1/2  0 -1/2;
 +
      -1/2 0  1/2;
 +
        0  0  0];
 +
 +
 +
l0 = [1/6; 1/6; 1/6];
 +
 +
 +
A = zeros(M);   
 +
b = zeros(M,1); 
 +
 +
for e = 1:Ne
 +
 
 +
    nodosK = nodos(e,:);
 +
   
 +
   
 +
    Bk = [coord4(nodosK(2),1)-coord4(nodosK(1),1), coord4(nodosK(3),1)-coord4(nodosK(1),1);
 +
          coord4(nodosK(2),2)-coord4(nodosK(1),2), coord4(nodosK(3),2)-coord4(nodosK(1),2)];
 +
   
 +
    detBk = det(Bk);       
 +
    Ck = inv(Bk) * (inv(Bk))';
 +
   
 +
   
 +
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
 +
   
 +
   
 +
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
 +
    A4=A;
 +
   
 +
   
 +
    fk = mean([ ...
 +
        f(coord4(nodosK(1),1), coord4(nodosK(1),2)), ...
 +
        f(coord4(nodosK(2),1), coord4(nodosK(2),2)), ...
 +
        f(coord4(nodosK(3),1), coord4(nodosK(3),2)) ]);
 +
   
 +
   
 +
    lk = fk * abs(detBk) * l0;
 +
   
 +
   
 +
    b(nodosK) = b(nodosK) + lk;
 +
end
 +
 +
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
 +
edges_sorted = sort(edges,2);
 +
[uniqE, ~, ic] = unique(edges_sorted,'rows');
 +
counts = accumarray(ic,1);
 +
boundary_edges = uniqE(counts==1,:);
 +
boundary_nodes = unique(boundary_edges(:));
 +
 +
for n = boundary_nodes'
 +
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
 +
    b(n) = 0;
 +
end
 +
 +
 +
 +
u4 = A\b; 
 +
 +
 +
figure(4)
 +
trisurf(nodos, coord4(:,1), coord4(:,2), u4)
 +
title('Solución numérica c_h')
 +
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
 +
 +
 +
 +
</source>
  
  

Revisión del 20:42 14 nov 2025

Trabajo realizado por estudiantes
Título Contaminación del lago (AMLS)
Asignatura MNEDP
Curso 2025-26
Autores
  • Antonio Lozano Fernández,
  • Miguel Cazorla Pedraza,
  • Lara lancelloti,
  • Sara Berti.
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Hola a todos esto e suna prueba


Códigos de Matlab

close all
clear all

%-------------------------------MALLA 1------------------------------------

%%%----------PASO 1: VARIABLES-----------

% Término fuente f(x,y)
% --- Fuente física: vertido puntual normalizado (gaussiana) ---
Q = 1.0;          % Caudal total de contaminante
x0 = 0.5; y0 = 0.5; % Punto de vertido del contaminante
sigma = 0.05;     % Anchura de la fuente

f = @(x,y) (Q/(2*pi*sigma^2)) .* exp(-((x-x0).^2 + (y-y0).^2) / (2*sigma^2));

% Parámetros físicos
D = 0.01;   % Coeficiente de difusión
k = 0.1;   % Coeficiente de reacción


%%%----------PASO 2: CREACIÓN DE LAS MATRICES DE NODOS Y TRIANGULOS-----------

Vertices = load('vertices1.txt');
Triangulos = load('triangulos1.txt');    %Cargamos los ficheros

Ne = length(Triangulos);     % número de Triangulos
M = length(Vertices);       % número de nodos (vertices)

% Coordenadas de los nodos
coord1 = zeros(M,2);
coord1(:,1) = Vertices(:,2);
coord1(:,2) = Vertices(:,3);

% Conectividad (qué nodos forman cada triángulo)
nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos1=nodos;
%%%----------PASO 3: MATRICES LOCALES DE REFERENCIA-----------

% Matriz de masa local (para el término de reacción)
L0 = [1/12 1/24 1/24;
      1/24 1/12 1/24;
      1/24 1/24 1/12];

% Matrices para los términos de derivadas parciales
Lxx = [1/2 -1/2 0;
       -1/2 1/2 0;
        0    0  0];

Lyy = [1/2  0 -1/2;
        0   0   0;
       -1/2 0  1/2];

Lxy = [1/2  0 -1/2;
       -1/2 0  1/2;
        0   0   0];

% Vector de integración de referencia (para f)
l0 = [1/6; 1/6; 1/6];


%%%----------PASO 4: ENSAMBLADO DE LA MATRIZ GLOBAL Y EL VECTOR DE CARGAS-----------

A = zeros(M);    % matriz global (rigidez + reacción)
b = zeros(M,1);  % vector global de cargas

for e = 1:Ne
    % Nodos del triangulo e
    nodosK = nodos(e,:);
    
    % Matriz del elemento (de coordenadas locales a reales)
    Bk = [coord1(nodosK(2),1)-coord1(nodosK(1),1), coord1(nodosK(3),1)-coord1(nodosK(1),1);
          coord1(nodosK(2),2)-coord1(nodosK(1),2), coord1(nodosK(3),2)-coord1(nodosK(1),2)];
    
    detBk = det(Bk);          % área*2 del triángulo
    Ck = inv(Bk) * (inv(Bk))'; % matriz Ck = (Bk^-1)(Bk^-T)
    
    % Matriz elemental (difusión + reacción)
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
    
    % Ensamblado en la matriz global
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
    A1=A;
    
    % Cálculo del término fuente medio en el triángulo
    fk = mean([ ...
        f(coord1(nodosK(1),1), coord1(nodosK(1),2)), ...
        f(coord1(nodosK(2),1), coord1(nodosK(2),2)), ...
        f(coord1(nodosK(3),1), coord1(nodosK(3),2)) ]);
    
    % Vector local de cargas
    lk = fk * abs(detBk) * l0;
    
    % Ensamblado en el vector global
    b(nodosK) = b(nodosK) + lk;
end

%%%----------PASO 5: APLICAR CONDICIONES DE CONTORNO DIRICHLET HOMOGÉNEAS (c=0)-----------

% Detectamos los nodos de frontera:
edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));

% Impongo c=0 en los nodos de la frontera
for n = boundary_nodes'
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
    b(n) = 0;
end

%%%----------PASO 6: RESOLUCIÓN DEL SISTEMA-----------

u1 = A\b;   % solución aproximada en los vertices

%%%----------PASO 7: GRAFICAMOS-----------

figure(1)
trisurf(nodos, coord1(:,1), coord1(:,2), u1)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight



%--------------------------------------MALLA 2--------------------------------------------------------


Vertices = load('vertices2.txt');
Triangulos = load('triangulos2.txt');    

Ne = length(Triangulos);     
M = length(Vertices);       


coord2 = zeros(M,2);
coord2(:,1) = Vertices(:,2);
coord2(:,2) = Vertices(:,3);

nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos2=nodos;

L0 = [1/12 1/24 1/24;
      1/24 1/12 1/24;
      1/24 1/24 1/12];


Lxx = [1/2 -1/2 0;
       -1/2 1/2 0;
        0    0  0];

Lyy = [1/2  0 -1/2;
        0   0   0;
       -1/2 0  1/2];

Lxy = [1/2  0 -1/2;
       -1/2 0  1/2;
        0   0   0];


l0 = [1/6; 1/6; 1/6];

A = zeros(M);    
b = zeros(M,1);  

for e = 1:Ne
    
    nodosK = nodos(e,:);
    
    
    Bk = [coord2(nodosK(2),1)-coord2(nodosK(1),1), coord2(nodosK(3),1)-coord2(nodosK(1),1);
          coord2(nodosK(2),2)-coord2(nodosK(1),2), coord2(nodosK(3),2)-coord2(nodosK(1),2)];
    
    detBk = det(Bk);         
    Ck = inv(Bk) * (inv(Bk))'; 
    
   
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
    
    
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
    A2=A;
    
   
    fk = mean([ ...
        f(coord2(nodosK(1),1), coord2(nodosK(1),2)), ...
        f(coord2(nodosK(2),1), coord2(nodosK(2),2)), ...
        f(coord2(nodosK(3),1), coord2(nodosK(3),2)) ]);
    
 
    lk = fk * abs(detBk) * l0;
    

    b(nodosK) = b(nodosK) + lk;
end

edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));


for n = boundary_nodes'
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
    b(n) = 0;
end


u2 = A\b;  

figure(2)
trisurf(nodos, coord2(:,1), coord2(:,2), u2)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight


%---------------------------------------MALA 3-----------------------------------------------------------




Vertices = load('vertices3.txt');
Triangulos = load('triangulos3.txt');    

Ne = length(Triangulos);     
M = length(Vertices);       


coord3 = zeros(M,2);
coord3(:,1) = Vertices(:,2);
coord3(:,2) = Vertices(:,3);


nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos3=nodos;

L0 = [1/12 1/24 1/24;
      1/24 1/12 1/24;
      1/24 1/24 1/12];


Lxx = [1/2 -1/2 0;
       -1/2 1/2 0;
        0    0  0];

Lyy = [1/2  0 -1/2;
        0   0   0;
       -1/2 0  1/2];

Lxy = [1/2  0 -1/2;
       -1/2 0  1/2;
        0   0   0];


l0 = [1/6; 1/6; 1/6];

A = zeros(M);   
b = zeros(M,1);  

for e = 1:Ne
    
    nodosK = nodos(e,:);
    
    
    Bk = [coord3(nodosK(2),1)-coord3(nodosK(1),1), coord3(nodosK(3),1)-coord3(nodosK(1),1);
          coord3(nodosK(2),2)-coord3(nodosK(1),2), coord3(nodosK(3),2)-coord3(nodosK(1),2)];
    
    detBk = det(Bk);          
    Ck = inv(Bk) * (inv(Bk))'; 
    
    
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
    
  
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
    A3=A;
    
    
    fk = mean([ ...
        f(coord3(nodosK(1),1), coord3(nodosK(1),2)), ...
        f(coord3(nodosK(2),1), coord3(nodosK(2),2)), ...
        f(coord3(nodosK(3),1), coord3(nodosK(3),2)) ]);
    
    
    lk = fk * abs(detBk) * l0;
    
   
    b(nodosK) = b(nodosK) + lk;
end

edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));


for n = boundary_nodes'
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
    b(n) = 0;
end



u3 = A\b;   



figure(3)
trisurf(nodos, coord3(:,1), coord3(:,2), u3)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight


%---------------------------------------MALLA 4-------------------------------------------------------------------

Vertices = load('vertices4.txt');
Triangulos = load('triangulos4.txt');    

Ne = length(Triangulos);    
M = length(Vertices);       


coord4 = zeros(M,2);
coord4(:,1) = Vertices(:,2);
coord4(:,2) = Vertices(:,3);


nodos = zeros(Ne,3);
nodos(:,1) = Triangulos(:,2);
nodos(:,2) = Triangulos(:,3);
nodos(:,3) = Triangulos(:,4);
nodos4=nodos;



L0 = [1/12 1/24 1/24;
      1/24 1/12 1/24;
      1/24 1/24 1/12];


Lxx = [1/2 -1/2 0;
       -1/2 1/2 0;
        0    0  0];

Lyy = [1/2  0 -1/2;
        0   0   0;
       -1/2 0  1/2];

Lxy = [1/2  0 -1/2;
       -1/2 0  1/2;
        0   0   0];


l0 = [1/6; 1/6; 1/6];


A = zeros(M);    
b = zeros(M,1);  

for e = 1:Ne
   
    nodosK = nodos(e,:);
    
    
    Bk = [coord4(nodosK(2),1)-coord4(nodosK(1),1), coord4(nodosK(3),1)-coord4(nodosK(1),1);
          coord4(nodosK(2),2)-coord4(nodosK(1),2), coord4(nodosK(3),2)-coord4(nodosK(1),2)];
    
    detBk = det(Bk);         
    Ck = inv(Bk) * (inv(Bk))'; 
    
    
    Lk = abs(detBk) * ( D*(Ck(1,1)*Lxx + Ck(1,2)*(Lxy+Lxy') + Ck(2,2)*Lyy) + k*L0 );
    
    
    A(nodosK,nodosK) = A(nodosK,nodosK) + Lk;
    A4=A;
    
    
    fk = mean([ ...
        f(coord4(nodosK(1),1), coord4(nodosK(1),2)), ...
        f(coord4(nodosK(2),1), coord4(nodosK(2),2)), ...
        f(coord4(nodosK(3),1), coord4(nodosK(3),2)) ]);
    
    
    lk = fk * abs(detBk) * l0;
    
    
    b(nodosK) = b(nodosK) + lk;
end

edges = [nodos(:,[1,2]); nodos(:,[2,3]); nodos(:,[3,1])];
edges_sorted = sort(edges,2);
[uniqE, ~, ic] = unique(edges_sorted,'rows');
counts = accumarray(ic,1);
boundary_edges = uniqE(counts==1,:);
boundary_nodes = unique(boundary_edges(:));

for n = boundary_nodes'
    A(n,:) = 0; A(:,n) = 0; A(n,n) = 1;
    b(n) = 0;
end



u4 = A\b;   


figure(4)
trisurf(nodos, coord4(:,1), coord4(:,2), u4)
title('Solución numérica c_h')
view(0,90); axis equal; colormap(turbo); colorbar; axis tight
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Contaminación_del_lago_(AMLS)&oldid=86228»