[Integrate]=THhomogenint(Cor,mainfunction,flag)

function [Integrate]=THhomogenint(Cor,mainfunction,flag)
% In  this  sub  function  is  description  to  homogen   integration  for
% triangular  areas.  Generally  application  is s electing Gauss-Legendre
% homogen integration method. Also, Gauss-Legendre  constants are derivate
% this homogen integrate technique. This function is using the direct
% homogen integration method and this applications are collate
% Gauss-Legendre #7# point integration method.
%
% Cor: Triangular Element Nodes Global Cartesian Coordinates
% Cor=[ 1   1   1        
%       X1  X2  X3   
%       Y1  Y2  Y3]));
%
% Mainfunction: Parametric main function or possible a parametric matrix
% Mainfunction=f(X,Y) this [X,Y] depends [e1,e2,e3] homogen axises
%
% flag=1:Parametric function  ,  flag=2:Homogen function


syms X Y e1 e2 e3 t1 real

px=Cor(:,1) ;
py=Cor(:,2) ;


%Triangular area
%Areal = 1/2*abs(det([ 1   1   1        
%                     X1  X2  X3   
%                     Y1  Y2  Y3 ]));

A=1/2*abs(det([ 1      1      1        
               px(1)  px(2)  px(3) 
               py(1)  py(2)  py(3) ]));

if flag==1
%Triangular edge node coordinates

%Cartesian axises are connecting homogen axises            
%X = e1*px(1) + e2*px(2) + e3*px(3)
%Y = e1*py(1) + e2*py(2) + e3*py(3)


%Homogen axises depends
%e2=(1-e1)*t1

%e3=(1-e1-e2) => 
%e3=(1-e1-(1-e1)*t1) =>
%e3=(1-e1)*(1-t1)



% Homogen axis(1)                       e1 = e1
% Homogen axis(2)                       e2 = (1-e1)*t1
% Homogen axis(3)                       e3 = (1-e1)*(1-t1)


Acon=subs(mainfunction,{X,Y}, ...
                       {e1*px(1) + e2*px(2) + e3*px(3) , ...
                        e1*py(1) + e2*py(2) + e3*py(3) }); 


Ah =subs(Acon,{e2,e3},...
              {(1-e1)*t1,...
               (1-e1)*(1-t1)});

end
if flag==2
    Ah =subs(mainfunction,{e2,e3},...
                          {(1-e1)*t1,...
                           (1-e1)*(1-t1)});
end

%Chain-Rules
%Integra1=int(2*A*Ah*(1-e1)*(1-e1)*(1-t2),e1,0,1);
%Integra2=int(Integra1,t1,0,1);
%Integra3=int(Integra2,t2,0,1);

Integrate=int(int(2*A*Ah*(1-e1),e1,0,1),t1,0,1);