%%%%%%%%% Bmatrix_2D4N  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the B matrix (derivative of shape functions), epsilon=B*u
% INPUTS
% 1.- chi=[chi(1) chi(2)]=[ji nu] where the natural derivatives will be evaluated
% 2.- xnod=[x1;y1;x2;y2;x3;y3;x4;y4], node coordinates
% OUTPUT: B matrix (size 3x8)
% Dependencies: calls 'XY_deriv_2D4N' (see below)
function B=Bmatrix_2D4N(chi,xnod)
  B=zeros(3,8);
  xyd=XY_deriv_2D4N(chi,xnod);
  for i=1:4
    B(1,2*i-1)=xyd(i,1);
    B(2,2*i)=xyd(i,2);
    B(3,2*i-1)=xyd(i,2);
    B(3,2*i)=xyd(i,1);
  end
return
end

% %%%%%% XY_deriv_2D4N  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Compute the derivative of shape functions respect real coordinates (x,y)
% % INPUT: chi=[chi(1) chi(2)]=[ji nu] where the derivatives will be evaluated
% % OUTPUT: is a matrix 2x4 containing the derivatives
% % Dependencies: calls 'natural_deriv_2D4N', 'Jacobian_2D4N' (see below)
% function xyd=XY_deriv_2D4N(chi,xnod)
% J=jacobian_2D4N(chi,xnod);
% invJ=inv(J);
% nd=natural_deriv_2D4N(chi);
% xyd=nd*invJ;
% return
% end

% %%%%%% natural_deriv_2D4N  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Compute the derivative of shape functions respect natural coordinates
% % INPUT: chi=[chi(1) chi(2)]=[ji nu] where the derivatives will be evaluated
% % OUTPUT: is a matrix 2x4 containing the derivatives
% function nd=natural_deriv_2D4N(ji)
% 
%  ji=chi(1);nu=chi(2);
%   d1_N1=-(1/4)*(1-nu);      % derivative of N1 respect ji
%   d2_N1=-(1/4)*(1-ji);      % derivative of N1 respect mu
%   d1_N2=+(1/4)*(1-nu);      % ...
%   d2_N2=-(1/4)*(1+ji);
%   d1_N3=+(1/4)*(1+nu);
%   d2_N3=+(1/4)*(1+ji);
%   d1_N4=-(1/4)*(1+nu);
%   d2_N4=+(1/4)*(1-ji);
%   nd=[d1_N1 d2_N1 ; d1_N2 d2_N2; d1_N3  d2_N3; d1_N4 d2_N4];
% return
% end

% %%%%%% jacobian_2D4N  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Compute the JACOBIAN of the transformation x,y --> ji,nu
% % QUADRILATERAL ELEMENT
% % input1: chi=[chi(1) chi(2)]=[ji nu] where the Jacobian will be evaluated
% % input2: global node coordinates [x1;y1;x2;y2;x3;y3;x4;y4]
% % output: jacobian matrix 2x2  [dx/dji dx/dnu;dy/dji dy/dnu];
% function jaco=jacobian_2D4N(chi,xnod)
% % N1=(1/4)*(1-ji)*(1-nu)
% % N2=(1/4)*(1+ji)*(1-nu)
% % N3=(1/4)*(1-ji)*(1+nu)
% % N4=(1/4)*(1+ji)*(1+nu)
%   ji=chi(1);nu=chi(2);
%   d1_N1=-(1/4)*(1-nu);      % derivative of N1 respect ji
%   d2_N1=-(1/4)*(1-ji);      % derivative of N1 respect mu
%   d1_N2=+(1/4)*(1-nu);      % ...
%   d2_N2=-(1/4)*(1+ji);
%   d1_N3=+(1/4)*(1+nu);
%   d2_N3=+(1/4)*(1+ji);
%   d1_N4=-(1/4)*(1+nu);
%   d2_N4=+(1/4)*(1-ji);
% 
%   jaco=zeros(2,2);
%   jaco(1,1)=[d1_N1 0 d1_N2 0 d1_N3 0 d1_N4 0]*xnod;  % dx/dji
%   jaco(2,1)=[0 d1_N1 0 d1_N2 0 d1_N3 0 d1_N4]*xnod;  % dy/dji
%   jaco(1,2)=[d2_N1 0 d2_N2 0 d2_N3 0 d2_N4 0]*xnod;  % dx/dnu
%   jaco(2,2)=[0 d2_N1 0 d2_N2 0 d2_N3 0 d2_N4]*xnod;  % dy/dnu
% return
% end