function [u] = ejemplobidimensional(n)
% function [u] = ejemplobidimensional(n)
%
% Resuelve
%    -d2u/dx2 - d2u/dy2 = 2*x*(x - 1) + 2*y*(y - 1),  0 <= x,y <= 1
%     u(x,0) = u(x,1) = u(0,y) = u(1,y) = 0
%
% Solucion exacta:
%     u  =  x*(x-1)*y*(y-1)
%
% Elementos finitos lineales
% Malla cuadrada n*n
%        n = cantidad de elementos por lado
%
%            14    15
%    13+-----+-----+-----+16
%      | \   | \   | \   |
%      |   \ |   \ |   \ |
%    9 +-----+-----+-----+12
%      | \   | \   | \   |
%      |   \ |   \ |   \ |
%    5 +-----+-----+-----+8
%      | \   | \   | \   |
%      |   \ |   \ |   \ |
%      +-----+-----+-----+
%      1     2     3     4

% cantidad de gdl total
nn = (n+1)*(n+1);
% cantidad de elementos total
nel = 2*n*n;

% inicializacion de matriz de rigidez y vector de cargas globales
A = zeros(nn,nn);
b = zeros(nn,1);

%  lazo sobre los elementos
%  trato primero los elementos inferiores, luego los superiores
for kx = 1:n
    for ky = 1:n
        %  1) Tratamiento elemento triangular "inferior"
        %
        %          3
        %           +
        %           | \
        %           |   \
        %           |     \
        %           +-------+
        %          1         2
        %
        %coordenadas nodo 1 del elemento k inferior
        x0 = (kx-1)*1/n;      y0 = (ky-1)*1/n;
        
        % matriz de rigidez elemental inferior
        Ainf = ...
            [2  -1  -1;
            -1   1   0;
            -1   0   1] / 2;
        
        %vector de cargas elemental inferior
        binf = [...
             -((x0^2-x0+y0^2-y0)/3*n^2 + (x0+y0-1)/6*n      + 1/15)/n^4;
             -((x0^2-x0+y0^2-y0)/3*n^2 + (4*x0+2*y0-3)/12*n + 2/15)/n^4;
             -((x0^2-x0+y0^2-y0)/3*n^2 + (2*x0+4*y0-3)/12*n + 2/15)/n^4];
        
        %vector de localizacion
        locel = [(ky-1)*(n+1)+kx  (ky-1)*(n+1)+kx+1  (ky)*(n+1)+kx];
        
        %ensamble de matriz de rigidez
        for i=1:3
            for j=1:3
                A(locel(i),locel(j)) = A(locel(i),locel(j)) + Ainf(i,j);
            end
        end
        
        %ensamble de vector de cargas
        for i=1:3
            b(locel(i)) = b(locel(i)) + binf(i);
        end
        %  2) Tratamiento elemento triangular "superior"
        %
        %          2         1
        %           +-------+
        %             \     |
        %               \   |
        %                 \ |
        %                   +
        %                    3
        
        %coordenadas nodo 1 del elemento superior
        x0 = kx*1/n;      y0 = ky*1/n;
        
        % matriz de rigidez elemental superior
        Asup = ...
            [2  -1  -1;
            -1   1   0;
            -1   0   1] / 2;
        
        %vector de cargas elemental superior
        bsup = [...
             -((x0^2-x0+y0^2-y0)/3*n^2 - (x0+y0-1)/6*n      + 1/15)/n^4;
             -((x0^2-x0+y0^2-y0)/3*n^2 - (4*x0+2*y0-3)/12*n + 2/15)/n^4;
             -((x0^2-x0+y0^2-y0)/3*n^2 - (2*x0+4*y0-3)/12*n + 2/15)/n^4];
        
        %vector de localizacion
        locel = [(ky)*(n+1)+kx+1  (ky)*(n+1)+kx  (ky-1)*(n+1)+kx+1];
        
        %ensamble de matriz de rigidez
        for i=1:3
            for j=1:3
                A(locel(i),locel(j)) = A(locel(i),locel(j)) + Asup(i,j);
            end
        end
        
        %ensamble de vector de cargas
        for i=1:3
            b(locel(i)) = b(locel(i)) + bsup(i);
        end
    end
end

% Lista de nodos fijos (nodos donde u=0)
fix  = [ 1:n+1 ...
    n+2:n+1:(n+1)*(n-1)+1  ...
    2*(n+1):n+1:(n+1)*(n)  ...
    n*(n+1)+1:(n+1)*(n+1)];

% solucion sistema de ecuaciones
%
% Modificacion sistema ecuaciones para
%    introducir condiciones de borde homogeneas

b(fix)     = zeros(size(fix,1));
A(:,fix)   = zeros(size(A,1),length(fix));
A(fix,:)   = zeros(length(fix),size(A,1));
A(fix,fix) = diag(ones(length(fix),1));

u = A \ b;
