% Computes the Element stiffness matrix
% INPUTS:
% 1.- Material properties (structure propMat: propMat.E, propMat.nu)
% 2.- Nodal coordinates xnod=[x1;y1;x2;y2;x3;y3;x4;y4] 
% 3.- Number of Gauss points used for the integration, npoints=1, 2 or 3
% OUTPUT: Element stiffness matrix of a 2D4N element (matrix 8x8)
% Dependencies: Calls 1) Bmatrix_2D4N, 2) jacobian_2D4N 

function stiff_elem=element_stiffness_2D4N(propMat,xnod,npoints)
%%% Definition number of Gauss points for integration
  if npoints==1
    w(1)=2;
    x(1)=0.;
  elseif npoints==2
    w(1)=1.;
    w(2)=1.;
    x(1)=-0.577350269189626;
    x(2)=-x(1);
  elseif npoints==3
    w(1)= 0.888888888888889;
    w(2)=0.555555555555556;
    w(3)=w(2);
    x(1)=0.;
    x(2)=-0.774596669241483;
    x(3)=-x(2);
  end

% Definition of constitutive matrix: hook law
  c1=propMat.E/(1-propMat.nu*propMat.nu);
  c2=propMat.nu*c1;
  G=propMat.E/(2*(1+propMat.nu));
  C=[c1 c2 0;c2 c1 0; 0 0 G];
 
  integral=zeros(8,8);
% Gauss sumation
  for i=1:npoints
    for j=1:npoints
      B=Bmatrix_2D4N([x(i) x(j)],xnod);
      fun=B'*C*B;
      detJac=det(jacobian_2D4N([x(i) x(j)],xnod));
      integral=integral+fun*detJac*w(i)*w(j);
    end
  end

  stiff_elem=integral;
  return
  
end