%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Ke,Re,Ie] =tetra_3d(e_mat,e_spa,i_var,mat)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% stiffness matrix Ke and righthand side Re for linear tets
% i_var:  [rho^0-rho_0^*]/rho_0^*    relativ change in reference density
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[nmat,ngrw] = size(mat);
 nod  = 4;
 ndim = 3;
delta = eye(ndim);

Ie = i_var;
Re = zeros(12, 1);
Ke = zeros(12,12);

indx  =[1;4;7;10];                    % indices of x coordinates
ex_mat=e_mat(indx);                   % material x coordinates of 1,2,3,4
ex_spa=e_spa(indx);                   % spatial  x coordinates of 1,2,3,4

indy  =[2;5;8;11];                    % indices of y coordinates
ey_mat=e_mat(indy);                   % material y coordinates of 1,2,3,4 
ey_spa=e_spa(indy);                   % spatial  y coordinates of 1,2,3,4

indz  =[3;6;9;12];                    % indices of z coordinates
ez_mat=e_mat(indz);                   % material z coordinates of 1,2,3,4 
ez_spa=e_spa(indz);                   % spatial  z coordinates of 1,2,3,4
 
% integration point %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
gp(1)=(ex_mat(1)+ex_mat(2)+ex_mat(3)+ex_mat(4)) / 4.0; 
gp(2)=(ey_mat(1)+ey_mat(2)+ey_mat(3)+ey_mat(4)) / 4.0; 
gp(3)=(ez_mat(1)+ez_mat(2)+ez_mat(3)+ez_mat(4)) / 4.0; 

w  = 1/6;

% shape functions in isoparametric space 
N(1)=1-gp(1)-gp(2)-gp(3);  
N(2)= +gp(1); 
N(3)=       +gp(2); 
N(4)=             +gp(3); 

% partial derivatives of shape functions wrt r and s and t
dNr(1,1)=-1;  dNr(1,2)=+1;  dNr(1,3)= 0;  dNr(1,4)= 0;
dNr(2,1)=-1;  dNr(2,2)= 0;  dNr(2,3)=+1;  dNr(2,4)= 0;
dNr(3,1)=-1;  dNr(3,2)= 0;  dNr(3,3)= 0;  dNr(3,4)=+1; 

% jacobian transformation matrix
JT=dNr*[ex_mat;ey_mat;ez_mat]';
detJ=det(JT);
if detJ<10*eps;    disp('Jacobi determinant less than zero!');   end;
JTinv=inv(JT);
dNx=JTinv*dNr;

F  = zeros(3,3);
for i=1:4
 indx=[3*i-2; 3*i-1; 3*i];
 F=F+e_spa(indx)'*dNx(:,i)';
end

var = i_var(1); 
if (ngrw==7); [A,P,var]=cnst_den(F,var,mat,ndim); end;  % density
if (ngrw==9); [A,P,var]=cnst_vol(F,var,mat,ndim); end;  % volume
Ie(1) = var;

for i=1:nod
  en=(i-1)*3;
  Re(en+ 1) =  Re(en+ 1) +(P(1,1)*dNx(1,i)' ...
						 + P(1,2)*dNx(2,i)' ...
						 + P(1,3)*dNx(3,i)') * detJ * w;
  Re(en+ 2) =  Re(en+ 2) +(P(2,1)*dNx(1,i)' ...
						 + P(2,2)*dNx(2,i)' ...	
						 + P(2,3)*dNx(3,i)') * detJ * w;
  Re(en+ 3) =  Re(en+ 3) +(P(3,1)*dNx(1,i)' ...
						 + P(3,2)*dNx(2,i)' ...	
						 + P(3,3)*dNx(3,i)') * detJ * w;
end  

for i=1:nod
for j=1:nod
  eni=(i-1)*3;
  enj=(j-1)*3;
  for r=1:ndim
  for s=1:ndim
    Ke(eni+r,enj+s)=Ke(eni+r,enj+s) ...
	+(dNx(1,i)*(A(r,1,s,1)*dNx(1,j)+A(r,1,s,2)*dNx(2,j)+A(r,1,s,3)*dNx(3,j)) ...
	 +dNx(2,i)*(A(r,2,s,1)*dNx(1,j)+A(r,2,s,2)*dNx(2,j)+A(r,2,s,3)*dNx(3,j)) ...
     +dNx(3,i)*(A(r,3,s,1)*dNx(1,j)+A(r,3,s,2)*dNx(2,j)+A(r,3,s,3)*dNx(3,j)))*detJ*w;
  end
  end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


