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

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

indx=[1;3;5;7];  ex_mat=e_mat(indx); 
indy=[2;4;6;8];  ey_mat=e_mat(indy);
 
% integration points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
g1=0.577350269189626; w1=1;
gp(:,1)=[-g1; g1;-g1; g1];   w(:,1)=[ w1; w1; w1; w1]; 
gp(:,2)=[-g1;-g1; g1; g1];   w(:,2)=[ w1; w1; w1; w1];
wp=w(:,1).*w(:,2);      xsi=gp(:,1);      eta=gp(:,2);

% shape functions and derivatives in isoparametric space %%%%%%%%%%%%%%%
N(:,1)=(1-xsi).*(1-eta)/4;   N(:,2)=(1+xsi).*(1-eta)/4;
N(:,3)=(1+xsi).*(1+eta)/4;   N(:,4)=(1-xsi).*(1+eta)/4;
dNr(1:2:8  ,1)=-(1-eta)/4;   dNr(1:2:8  ,2)= (1-eta)/4;
dNr(1:2:8  ,3)= (1+eta)/4;   dNr(1:2:8  ,4)=-(1+eta)/4;
dNr(2:2:8+1,1)=-(1-xsi)/4;   dNr(2:2:8+1,2)=-(1+xsi)/4;
dNr(2:2:8+1,3)= (1+xsi)/4;   dNr(2:2:8+1,4)= (1-xsi)/4;

% jacobian transformation matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
JT=dNr*[ex_mat;ey_mat]';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop over all integration points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ip=1:4
  indx=[2*ip-1; 2*ip];
  detJ=det(JT(indx,:));
  if detJ<10*eps;    disp('Jacobi determinant less than zero!');   end;
  JTinv=inv(JT(indx,:));
  dNx=JTinv*dNr(indx,:);

  F=zeros(2,2);
  for j=1:4
    jndx=[2*j-1; 2*j];
    F=F+e_spa(jndx)'*dNx(:,j)';
  end
  
  var = i_var(ip);  
  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(ip)=var;

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

  for i=1:nod
  for j=1:nod
    eni=(i-1)*2;
    enj=(j-1)*2;
  
    Ke(eni+1,enj+1)=Ke(eni+1,enj+1)+(dNx(1,i)*A(1,1,1,1)*dNx(1,j) ...
                                    +dNx(1,i)*A(1,1,1,2)*dNx(2,j) ...
                                    +dNx(2,i)*A(1,2,1,1)*dNx(1,j) ...
                                    +dNx(2,i)*A(1,2,1,2)*dNx(2,j))*detJ*wp(ip);
    Ke(eni+1,enj+2)=Ke(eni+1,enj+2)+(dNx(1,i)*A(1,1,2,1)*dNx(1,j) ...
                                    +dNx(1,i)*A(1,1,2,2)*dNx(2,j) ...
                                    +dNx(2,i)*A(1,2,2,1)*dNx(1,j) ...
                                    +dNx(2,i)*A(1,2,2,2)*dNx(2,j))*detJ*wp(ip);
    Ke(eni+2,enj+1)=Ke(eni+2,enj+1)+(dNx(1,i)*A(2,1,1,1)*dNx(1,j) ...
                                    +dNx(1,i)*A(2,1,1,2)*dNx(2,j) ...
                                    +dNx(2,i)*A(2,2,1,1)*dNx(1,j) ...
                                    +dNx(2,i)*A(2,2,1,2)*dNx(2,j))*detJ*wp(ip);
    Ke(eni+2,enj+2)=Ke(eni+2,enj+2)+(dNx(1,i)*A(2,1,2,1)*dNx(1,j) ...
                                    +dNx(1,i)*A(2,1,2,2)*dNx(2,j) ...
                                    +dNx(2,i)*A(2,2,2,1)*dNx(1,j) ...
                                    +dNx(2,i)*A(2,2,2,2)*dNx(2,j))*detJ*wp(ip);
  end
  end

end
% loop over all integration points %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%