function [K,M] = fem(B,H,Le,E,rho,nx,ny,ncnt,NELEM,NNODE,dispVal,dispID)
%clear
%hold off
%clc
%clg
%axis('square');

Ae = B.*H;        % dot if H is an array
Je = B.*H.^3/12;  % dot if H is an array

K = zeros(2*NNODE,2*NNODE);
M = zeros(2*NNODE,2*NNODE);

for ie=1:NELEM,

  eye = ncnt(ie,1);
  jay = ncnt(ie,2);
  
% Form the local stiffness matrix, [kbar].
  L = Le(ie);
  A = Ae(ie);
  J = Je(ie);

  k = [ 12    6*L      -12   6*L;    ...

        6*L   4*L^2    -6*L  2*L^2;  ... 

       -12   -6*L       12   -6*L;   ...

        6*L   2*L^2    -6*L  4*L^2 ];


  kbar = k*(J*E/L^3);

% Form the local mass matrix
 
 m_local = [13/35+6*J/5/A/L^2  11*L/210+J/10/A/L  9/70-6*J/5/A/L^2  -13/420+J/10/A/L; ...
            11*L/210+J/10/A/L  L^2/105+2*J/15/A   13*L/420-J/10/A/L -L^2/140-J/30/A;  ...
            9/70-6*J/5/A/L^2   13*L/420-J/10/A/L   13/35+6*J/5/A/L^2 -11*L/210-J/10/A/L; ...
           -13/420+J/10/A/L   -L^2/140-J/30/A    -11*L/210-J/10/A/L  L^2/105+2*J/15/A]; 

 Mbar = rho*A*L*m_local;

% Form ID matrix to assemble kbar into the global stiffness matrix, K and
% global inertia matrix, M
  id1 = 2*(eye-1) + 1; 
  id2 = id1 + 1;
  id3 = 2*(jay-1) + 1;
  id4 = id3 + 1;
 

  K(id1,id1) = K(id1,id1) + kbar(1,1);
  K(id1,id2) = K(id1,id2) + kbar(1,2);
  K(id2,id1) = K(id2,id1) + kbar(2,1);
  K(id2,id2) = K(id2,id2) + kbar(2,2);


  K(id1,id3) = K(id1,id3) + kbar(1,3);
  K(id1,id4) = K(id1,id4) + kbar(1,4);
  K(id2,id3) = K(id2,id3) + kbar(2,3);
  K(id2,id4) = K(id2,id4) + kbar(2,4);

  K(id3,id1) = K(id3,id1) + kbar(3,1);
  K(id3,id2) = K(id3,id2) + kbar(3,2);
  K(id4,id1) = K(id4,id1) + kbar(4,1);
  K(id4,id2) = K(id4,id2) + kbar(4,2);


  K(id3,id3) = K(id3,id3) + kbar(3,3);
  K(id3,id4) = K(id3,id4) + kbar(3,4);
  K(id4,id3) = K(id4,id3) + kbar(4,3);
  K(id4,id4) = K(id4,id4) + kbar(4,4);
%-------------------------------------[K] Completed


  M(id1,id1) = M(id1,id1) + Mbar(1,1);
  M(id1,id2) = M(id1,id2) + Mbar(1,2);
  M(id2,id1) = M(id2,id1) + Mbar(2,1);
  M(id2,id2) = M(id2,id2) + Mbar(2,2);


  M(id1,id3) = M(id1,id3) + Mbar(1,3);
  M(id1,id4) = M(id1,id4) + Mbar(1,4);
  M(id2,id3) = M(id2,id3) + Mbar(2,3);
  M(id2,id4) = M(id2,id4) + Mbar(2,4);

  M(id3,id1) = M(id3,id1) + Mbar(3,1);
  M(id3,id2) = M(id3,id2) + Mbar(3,2);
  M(id4,id1) = M(id4,id1) + Mbar(4,1);
  M(id4,id2) = M(id4,id2) + Mbar(4,2);


  M(id3,id3) = M(id3,id3) + Mbar(3,3);
  M(id3,id4) = M(id3,id4) + Mbar(3,4);
  M(id4,id3) = M(id4,id3) + Mbar(4,3);
  M(id4,id4) = M(id4,id4) + Mbar(4,4);

%------------------------------------[M] completed

end

% Store unlaltered K as Ksing and M as Morig
  Ksing = K;
  Morig = M;

% Displacement boundary conditions
% --------------------------------
% dispID array contains the dof which are assigned values contained in dispVal.
[sm,sn] = size(dispID);
Ndbc = sn;

for nd=1:Ndbc,
 K = matcut(K,dispID(nd)-nd+1);
 M = matcut(M,dispID(nd)-nd+1);
end



