clear
syms x y x0 y0 real
syms h n positive

%    h  -->  tamaņo de los elementos (1/n)
%    n  -->  cantidad de elementos por lado del cuadrado
%    f = - 2*x*(x - 1) - 2*y*(y - 1);

h  = 1/n;
f  = - 2*x*(x - 1) - 2*y*(y - 1);

%
%  1) Calcula int (phi_i * f) dA para el elemento 
%     triangular "inferior" : 
%
%          3
%           +
%           | \
%           |   \
%           |     \
%           +-------+
%          1         2
%
%   (x0,y0) -->  coordenadas nodo 1
%

phi1 = 1 -(x-x0)/h -(y-y0)/h;
phi2 =  (x-x0)/h;
phi3 =  (y-y0)/h;


f1x    = int(phi1*f,x);
f1xint = simplify(subs(f1x,x,x0+y0+h-y) - subs(f1x,x,x0));
f1y    = int(f1xint,y);
f1inf  = simplify(subs(f1y,y,y0+h) - subs(f1y,y,y0));
pretty (f1inf)

f2x    = int(phi2*f,x);
f2xint = simplify(subs(f2x,x,x0+y0+h-y) - subs(f2x,x,x0));
f2y    = int(f2xint,y);
f2inf  = simplify(subs(f2y,y,y0+h) - subs(f2y,y,y0));
pretty (f2inf)

f3x    = int(phi3*f,x);
f3xint = simplify(subs(f3x,x,x0+y0+h-y) - subs(f3x,x,x0));
f3y    = int(f3xint,y);
f3inf  = simplify(subs(f3y,y,y0+h) - subs(f3y,y,y0));
pretty (f3inf)

%
%  2) Calcula int (phi_i * f) dA para el elemento 
%     triangular "superior" : 
%
%          2         1
%           +-------+
%             \     |
%               \   |
%                 \ |
%                   +
%                    3
%
%   (x0,y0) -->  coordenadas nodo 1
%

phi1 = 1 + (x-x0)/h + (y-y0)/h;
phi2 = -(x-x0)/h;
phi3 = -(y-y0)/h;


f1x    = int(phi1*f,x);
f1xint = simplify(subs(f1x,x,x0) - subs(f1x,x,x0+y0-h-y));
f1y    = int(f1xint,y);
f1sup  = simplify(subs(f1y,y,y0) - subs(f1y,y,y0-h));
pretty (f1sup)

f2x    = int(phi2*f,x);
f2xint = simplify(subs(f2x,x,x0) - subs(f2x,x,x0+y0-h-y));
f2y    = int(f2xint,y);
f2sup  = simplify(subs(f2y,y,y0) - subs(f2y,y,y0-h));
pretty (f2sup)

f3x    = int(phi3*f,x);
f3xint = simplify(subs(f3x,x,x0) - subs(f3x,x,x0+y0-h-y));
f3y    = int(f3xint,y);
f3sup  = simplify(subs(f3y,y,y0) - subs(f3y,y,y0-h));
pretty (f3sup)
