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

function [Integrate]=THhomogenint(Cor,mainfunction,flag)
% In this sub function is description to homogen integration for
% tetrahedral volumes. Generally application is selecting 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 #5# point integration method.
%
% Cor: Tetrahedral Element Node Global Cartesian Coordinates
% Cor=[ 1   1   1   1     
%       X1  X2  X3  X4 
%       Y1  Y2  Y3  Y4 
%       Z1  Z2  Z3  Z4 ]));
%
% Mainfunction: Parametric main function or possible a parametric matrix
% Mainfunction=f(X,Y,Z) this [X,Y,Z] depends [e1,e2,e3,e4] homogen axises
%
% flag=1:Parametric function  ,  flag=2:Homogen function


syms X Y Z e1 e2 e3 e4 t1 t2 real

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

%Tetrahedral volume
%Volume=1/6*abs(det([ 1   1   1   1     
%                     X1  X2  X3  X4 
%                     Y1  Y2  Y3  Y4 
%                     Z1  Z2  Z3  Z4 ]));

V=1/6*abs(det([ 1      1      1      1     
               px(1)  px(2)  px(3)  px(4) 
               py(1)  py(2)  py(3)  py(4) 
               pz(1)  pz(2)  pz(3)  pz(4) ]));

if flag==1
%Tetrahedral edge node coordinates

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


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

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

%e4=(1-e1-e2-e3)=>
%e4=(1-e1-(1-e1)*t2-(1-e1)*(1-t2)*t1) =>
%e4=(1-e1)*[(1-t2)*(1-t1)]

%                                       e1 = e1
%                                       e2 = (1-e1)*t2
%                                       e3 = (1-e1)*(1-t2)*t1
%                                       e4 = (1-e1)*[(1-t2)*(1-t1)]

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


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

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

%Chain-Rules
%Integra1=int(6*V*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(int(6*V*Ah*(1-e1)*(1-e1)*(1-t2),e1,0,1),t1,0,1),t2,0,1);
%               This -6*V- is obtain homogen integration constant