Distribución del Potencial Eléctrico Cerebral (Modelo 2D) - CAN

Trabajo realizado por estudiantes
Título	Distribución del Potencial Eléctrico Cerebral (Modelo 2D) - CAN
Asignatura	MNEDP
Curso	2025-26
Autores	Claudia Domínguez, Ignacio Martínez, Analía Olivero
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Póster del Trabajo

También es posible ver el póster en PDF a resolución completa. Se puede acceder a los contenidos del código QR mediante este enlace.

2 Códigos de MatLab

2.1 Código: make_concentric_mesh.m

Esta función se encarga de generar el mallado.

function [points, cells, regions] = make_concentric_mesh( ...
        R, r1, ntheta_outer, ntheta_inner, nr_inner, nr_outer, exc, save_path)

    %esta función genera el mallado de la elipse añadiendo radios y elipses
    %semejantes más pequeñas interiores (anillos elípticos)

    % exc = semieje vertical / semieje horizontal
    

    % radios
    r_inner = linspace(0, r1, nr_inner + 1);
    r_outer = linspace(r1, R, nr_outer + 1);
    radii   = [r_inner, r_outer(2:end)];
    k_if    = nr_inner;  % índice de la circunferencia interfaz
    n_rings = numel(radii);

    % listas para guardar los datos
    points = [];
    %como diferentes anillos tienen diferentes tamaños (los exteriores y
    %los interiores) no nos vale con un vector/array, así que usamos celdas
    %  (más parecido a una lista de Python)
    ring_nodes = cell(1, n_rings);
    %esto genera una secuencia de celdas 1xn_rings, que contendrá los
    %índices de los nodos de nuestros anillos

    % k = 0 es el centro (origen)
    points = [0 0];
    ring_nodes{1} = 1;
    
    %ángulos
    thetas_outer = linspace(0, 2*pi, ntheta_outer+1);
    thetas_outer(end) = [];
    thetas_inner = linspace(0, 2*pi, ntheta_inner+1);
    thetas_inner(end) = [];

    % anillos interiores
    for k = 2:k_if
        r = radii(k);
        x = r * cos(thetas_inner);
        % escalado vertical (exc)
        y = exc * r * sin(thetas_inner);   
        
        %guaradamos los puntos y los nodos de este anillo
        start_idx = size(points,1) + 1;
        pts = [x(:), y(:)];
        points = [points; pts];
        ring_nodes{k} = start_idx : (start_idx + ntheta_inner - 1);
    end

    % anillos exteriores
    for k = k_if+1:n_rings
        r = radii(k);
        x = r * cos(thetas_outer);
        %escalado horizontal
        y = exc * r * sin(thetas_outer);  

        start_idx = size(points,1) + 1;
        pts = [x(:), y(:)];
        points = [points; pts];
        ring_nodes{k} = start_idx : (start_idx + ntheta_outer - 1);
    end

    % triangulación: recorremos los anillos, que tienen el mismo número de
    % puntos cada uno, tomando dos nodos del anillo interior y sus dos
    % corrrespondientes del anillo esterior, formamos un cuadrilátero y lo
    % dividimos en dos, lo que genera una triangulación.
    cells = [];
    regions = [];

    % 1. "estrella" interior (emanan del origen)
    for j = 1:ntheta_inner
        a = ring_nodes{1}(1);
        b = ring_nodes{2}(j);
        c = ring_nodes{2}(mod(j,ntheta_inner)+1);

        cells(end+1,:) = [a, b, c];
        regions(end+1) = 1;
    end

    % triangulación de Omega_1 (menos el centro) - anillos interiores,
    % Omega_1
    for k = 2:k_if
        
        inner = ring_nodes{k}; %anillo interior (de Omega_1)
        outer = ring_nodes{k+1}; %anillo exterior (de Omega_1)
        N = ntheta_inner; %número de nodos en cada anillo

        for j = 1:N
            a = inner(j); %un nodo interior
            b = inner(mod(j,N)+1); % el siguiente nodo interior
            c = outer(j); %un nodo exterior
            d = outer(mod(j,N)+1); %el siguiente nodo exterior

            cells(end+1,:) = [a b d]; %primer triángulo
            regions(end+1) = 1; %marcamos este triángulo como uno de Omega_1

            cells(end+1,:) = [a d c]; %segundo triángulo
            regions(end+1) = 1;%lo marcamos como de Omega_1
        end
    end

    % triangulación de Omega_2
    for k = k_if+1:n_rings-1
        inner = ring_nodes{k}; %anillo interior de Omega_2
        outer = ring_nodes{k+1}; %anillo exterior de Omega_2
        N = ntheta_outer; %número de nodos

        for j = 1:N
            a = inner(j); %un nodo interior
            b = inner(mod(j,N)+1); %el siguiente nodo interior
            c = outer(j); %un nodo exterior
            d = outer(mod(j,N)+1); %el siguiente nodo exterior

            cells(end+1,:) = [a b d]; %primer triángulo del cuadrilátero
            regions(end+1) = 2;%lo marcamos como de Omega_2

            cells(end+1,:) = [a d c]; %segundo triángulo
            regions(end+1) = 2;%lo marcamos como de Omega_2
        end
    end

    % lo guardamos
    save(save_path, "points", "cells", "regions");

end

2.2 Código: calculate_con_mesh_params.m

Este código se encarga de escoger unos parámetros concretos para el mallado a partir de un valor concreto de [math]h[/math].

function [ntheta_inner, ntheta_outer, nr_inner, nr_outer] = calculate_con_mesh_params(h, R, r1)
    %calcula los parámetros de generate_mesh necesarios para conseguir una
    %cierta "h" (es un poco heurístico e impreciso, pero funciona!)
    circum_inner = 2*pi*r1;
    circum_outer = 2*pi*R;

    %divisiones angulares
    ntheta_inner = max(3, ceil(circum_inner / h));
    ntheta_outer = max(3, ceil(circum_outer / h));

    %divisiones de tipo anillo
    nr_inner = max(1, ceil(r1 / h));
    nr_outer = max(1, ceil((R - r1) / h));
end

2.3 Código: generate_mesh.m

Este código se encarga de generar el mallado para una [math]h[/math] concreta y una "excentricidad" dada.

function generate_mesh(h, name, exc)
    %calculamos los valores que nos dan esta h mínima
    [ntheta_inner, ntheta_outer, nr_inner, nr_outer] = ...
        calculate_con_mesh_params(h, 1.0, 0.95);
    
    ntheta_outer = max([3, ntheta_inner, ntheta_outer]);
    ntheta_inner = ntheta_outer;
    
    %generamos el mallado (y se guarda)
    make_concentric_mesh( ...
        1.0, 0.95, ntheta_outer, ntheta_inner, nr_inner, nr_outer, ...
        exc, name + ".mat");
end

2.4 Código: plot_mesh.m

Este código representa el mallado que hemos generado con las funciones anteriores, diferenciando en las dos regiones [math]\Omega_1[/math] y [math]\Omega_2[/math].

function plot_mesh(mesh, h)
    %representación del mallado en dos partes Omega_1 y Omega_2

    %cargamos cada parte del mallado
    T = mesh.cells;
    P = mesh.points;
    R = mesh.regions;
    
    % Índices de triángulos de cada región
    idx1 = (R == 1);
    idx2 = (R == 2);
    
    figure('Position',[100 100 900 700]); hold on;
    
    %hacemos un triplot de los triángulos de Omega_1 (a un color)
    triplot(T(idx1,:), P(:,1), P(:,2), color = '#80B3FF');
    %hacemos otro triplot de los triángulos de Omega_2 (a otro color)
    triplot(T(idx2,:), P(:,1), P(:,2), 'blue');
    
    title("Representación del mallado por regiones");
    xlabel('X');
    ylabel('Y');
    axis equal;
    saveas(gcf, "mesh_representation h="+num2str(h)+".pdf");
    hold off;
end

2.5 Código: plot_solution.m

Este código representa la solución que calculan los códigos posteriores.

function plot_solution(P, T, u, elev, azim, levels, prefijo)
    
    %Representamos la superficie
    figure('Name','Rerepsetnación como superficie');
    %hacemos un trisurf sin marcar los triángulos (queda feo si hay muchos)
    trisurf(T, P(:,1), P(:,2), u, 'EdgeColor','none');
    xlabel('x'); ylabel('y'); zlabel('u');
    %los parámetros de vista iniciales son parámetros
    view(elev, azim);
    axis equal; axis tight;
    colorbar;
    title('Solución FEM (superficie 3D)');

    %lo guardamos
    saveas(gcf, prefijo + "_surface.png");
    
    %otra figura para el mapa de calor de los niveles
    figure('Name','Representación como contorno');

    %como no podemos hacerlo circular (no sabemos), lo hacemos cuadrado
    xg = linspace(min(P(:,1)), max(P(:,1)), 300);
    yg = linspace(min(P(:,2)), max(P(:,2)), 300);

    [X,Y] = meshgrid(xg, yg);
    
    %función para interpolar los datos de los nodos al resto del mallado
    Fint = scatteredInterpolant(P(:,1), P(:,2), u, 'natural');
    Ugrid = Fint(X,Y);
    
    %contornos (sin líneas en los bordes)
    contourf(X, Y, Ugrid, levels, 'LineColor','none');

    axis equal; axis tight;
    xlabel('x'); ylabel('y');
    colormap('jet');
    colorbar;
    title('Solución FEM (contorno)');
    %lo guardamos
    saveas(gcf, prefijo + "_contour.png");

end

2.6 Código: assemble_system.m

Este código ensambla el sistema lineal que resulta del método de elementos finitos.

function [A, F, P, T] = assemble_system(mesh_path, f_rhs)

    data = load(mesh_path, "points", "cells", "regions");
    P = data.points;    % Nx2
    T = data.cells;     % Mx3
    R = data.regions;   % Mx1 (región de cada elemento)

    N = size(P, 1);
    M = size(T, 1);

    %conductividades
    sigma1 = 0.33;   % cerebro
    sigma2 = 0.16;   % cráneo

    %matriz global y vector de cargas
    A = sparse(N, N);
    F = zeros(N, 1);

    %bucle sobre elementos
    for k = 1:M
        nodes  = T(k, :);        %índices de los nodos del triángulo, 1x3
        coords = P(nodes, :);    %coordenada, 3x2

        x1 = coords(1,1); 
        y1 = coords(1,2);
    
        x2 = coords(2,1); 
        y2 = coords(2,2);
    
        x3 = coords(3,1); 
        y3 = coords(3,2);
    
        area = 0.5 * abs( (x2-x1)*(y3-y1) - (y2-y1)*(x3-x1) );

        % Escogemos sigma según región
        if R(k) == 1
            sigma = sigma1;
        else
            sigma = sigma2;
        end
        
        %matriz BK (igual que en clase)
        BK = [coords(2,1) - coords(1,1), coords(3,1) - coords(1,1);
                       coords(2,2) - coords(1,2), coords(3,2) - coords(1,2)];
       
        %matriz de rigidez local (igual que en clase, pero en otra función
        %para que sea más limpio)
        Ke = local_stiffness_L(2*area, BK, sigma);

        %calculamos la aportación al vector de cargas diferenciando entre
        %la región interior y la región exterior
        if R(k) == 1
            %hacemos la media de la f, igual que en clase
            fk = 1/3 * (f_rhs(x1,y1) + f_rhs(x2,y2)+ f_rhs(x3,y3));
            fe = local_load_L(BK, fk);
        else
            fe = zeros(3, 1);
        end

        % Añadimos las aportaciones al vector de carga y a la matriz de
        % rigidez de este triángulo
        F(nodes) = F(nodes) + fe;
        A(nodes, nodes) = A(nodes, nodes) + Ke;
    end
end




function Ke = local_stiffness_L(alpha, BK, sigma)
    % Matriz de rigidez local
    % Matriz de deformación Ck = inv(BK) @ inv(BK).T
    invBk = inv(BK);

    Ck = invBk*invBk';

    % Matrices L
    Lxx = 0.5 * [1, -1, 0; -1, 1, 0; 0, 0, 0];
    Lyy = 0.5 * [1, 0, -1; 0, 0, 0; -1, 0, 1];
    Lxy = 0.5 * [1, -1, 0; 0, 0, 0; -1, 1, 0];

    % Matriz local de rigidez
    LK = Ck(1,1)*Lxx + Ck(1,2)*(Lxy + Lxy') + Ck(2,2)*Lyy;
    Ke = sigma * 0.5 * abs(alpha) * LK;
end


function fe = local_load_L(BK, fK)
    %vector de carga local de 3x1 usando:
    % fe = fK * |det(BK)| * ell0
    %donde ell0 = [1/6; 1/6; 1/6] (apuntes)

    ell0 = [1/6; 1/6; 1/6];
    area2 = abs(det(BK));  % |det(BK)| = 2 * area
    fe = fK * area2 * ell0;
end

2.7 Código: main.m

La función principal, aplica las funciones anteriores a un problema concreto, generando el mallado y usándolo con una [math]f[/math] concreta, definida en el mismo código.

%escogemos una h y un valor de "excentricidad"
h = 0.025;
%esta "excentricidad" no es rigurosa, es simplemente un parámetro para
%regular cómo de achatada es la elipse de nuestra región, ¡no todas las
%cabezas son iguales!
excentricidad = 1; 

% generamos el mallado
mesh_name = "head_mesh";
generate_mesh(h, mesh_name, excentricidad)
mesh_path = mesh_name + ".mat";

%cargamos el mallado desde el archivo
mesh = load("head_mesh.mat", "points", "cells", "regions");

plot_mesh(mesh, h)

%parámetros para el doble dipolo
%intensidad del primer dipolo suave
A1 = 2; 
%intensidad del segundo dipolo suave
A2 = -2;
%"radio" inverso del dipolo
TD = 0.02;

%definimos la función f
function val = f_rhslack(x, y, A1, A2, TD)
    % % Fuente intracraneal f(x,y).
    % val = x.^2 - y.^2 ;

    x0 = 0; y0 = 0;
    val = A1 * (x./ TD).*exp(-((x-x0).^2 + (y-y0).^2) / TD);
    % x0 = -0.1; y0 = -0.1;
    % val = val + A2 * (y./ TD).*exp(-((x-x0).^2 + (y-y0).^2) / TD);
end
%handle de función con los valores de A1, A2 y TD definidos arriba
f_rhs = @(x,y) f_rhslack(x,y, A1, A2, TD);

%ensamblamos el sistema
[A, F, P, T] = assemble_system(mesh_path, f_rhs);

% Resolución del sistema
u = A \ F;

name = "SolucionFEMcabeza";

%Representación de la solución
plot_solution(P, T, u, 35, -60, 30, name)