%==================
% ROUTINE Strain.m
%==================
   NNPE2=2*NNPE;
   NUMEQ=2*NUMNP;
   
   for I=1:NUMNP;
      IX=2*I-1;
      IY=IX+1; 
      F(IX,1)=XBC(I); 
      F(IY,1)=YBC(I); 
   end

%  ------------------------------- 
%  Define useful matrices to use
%  later.  BLNK = blank, i.e. 0.0
%  Id = Identity matrix.
%  ------------------------------- 
   BLNKv(1:NNPE)=0.0;
   BLNKs(1:NNPS)=0.0;
   Id = diag(ones(2*NNPE,1),0);

%  ------------------------------- 
%  Rotate specified traction to
%  global coordinate for surface
%  integration.
%  -------------------------------- 
   for I=1:NUMNP
     if(THETA(I)~=0)
        TH=THETA(I);
        C=cos(TH); 
        S=sin(TH); 
        tx=C*TX(I)-S*TY(I);
        ty=S*TX(I)+C*TY(I);
        TX(I)=tx;
        TY(I)=ty;
     end
   end

%  ------------------------------- 
%  FORMATION OF ELEMENT S-MATRICES 
%  -------------------------------- 
   for I=1:NUMEL

      FE(1:NNPE2,1)=0.0;
      SE(1:NNPE2,1:NNPE2)=0.0; 
      for J=1:NNPE 
        NPJ=NP(I,J); 
        NPE(J)=NPJ*2-1;
        J1=J+NNPE; 
        NPE(J1)=NPJ*2; 
      end
      rc=NP(I,:);

%     ----------------------------------------
%     Begin volume quadrature for each element
%     ----------------------------------------
      JEND=NUMQPT(1);
      for J=1:JEND; 
         XJ=0;
         YJ=0;
         RJAC=[0,0;0,0];
         SFj  = SF(1,:,J);
         DNDu = SF(2,:,J);
         DNDv = SF(3,:,J);

%        ---------------------------------
%        Determine coordinate and Jacobian
%        ---------------------------------
         XJ=SF(1,:,J)*XORD(rc)';
         YJ=SF(1,:,J)*YORD(rc)';
         RJAC(1,1)=DNDu*(XORD(rc))';
         RJAC(1,2)=DNDv*(XORD(rc))';
         RJAC(2,1)=DNDu*(YORD(rc))';
         RJAC(2,2)=DNDv*(YORD(rc))';
         DETJ=det(RJAC);
         if DETJ <=   0
           fprintf(1,'\n--------------------------')
           fprintf(1,'\n Error in steady.m        ') 
           fprintf(1,'\n DETJ =%7e',DETJ           )
           fprintf(1,'\n must be > 0.0'            )
           fprintf(1,'\n--------------------------\n')
           error
         end
%        -----------------------------
%        Determine inverse of Jacobian
%        -----------------------------
         RJACI=inv(RJAC);
%        ----------------------------------------------------
%        Determine derivative of shape functions in X-Y plane
%        ----------------------------------------------------
         DNDX=RJACI(1,1)*DNDu+RJACI(2,1)*DNDv;
         DNDY=RJACI(1,2)*DNDu+RJACI(2,2)*DNDv;
%        ---------------------------------
%        Include user written elastic
%        constants:
%        ---------------------------------
         COEF
         if IPROG == 1
            % Plane strain
            R(1,1)=2.0*G+LAMBDA;
            R(2,2)=R(1,1); 
            R(3,3)=2*G; 
            R(4,4)=2*G; 
            R(1,2)=LAMBDA; 
            R(2,1)=LAMBDA;
         elseif IPROG == 2 
            % Plane stress
            R(1,1)=4.0*G*(G+LAMBDA)/(2.0*G+LAMBDA);
            R(1,2)=2.0*G*LAMBDA/(2.0*G+LAMBDA);
            R(2,1)=R(1,2); 
            R(2,2)=R(1,1); 
            R(3,3)=2*G; 
            R(4,4)=2*G; 
         end

%        ---------------------------------
%        CALCULATE QUADRATURE CONTRIBUTION 
%        TO SE-MATRIX
%        ---------------------------------
         Np=[DNDX,      BLNKv;...
             BLNKv,     DNDY;...
             DNDY/2,    DNDX/2;...
             DNDY/2,    DNDX/2];

         SE=SE+WT(1,J)*Np'*R*Np*DETJ;

         Nf=[ SFj, BLNKv; BLNKv, SFj ];
         Bf=[GAMX;GAMY];
         FE=FE+WT(1,J)*Nf'*Bf*DETJ;
      end
%     -----------    end of volume quadrature

%     ------------------------------ 
%     BEGIN SURFACE QUADRATURE ON 
%     EACH SIDE OF ELEMENT I
%     ------------------------------ 
      for J=1:NSPE
%        --------------------------------
%        CHECK IF QUADRATURE IS NECESSARY 
%        --------------------------------
         CHKTX=1.0;
         CHKTY=1.0;
         clear nps;
         for K=1:NNPS
            J1=NP(I,NPSIDE(J,K));
            CHKTX=CHKTX*TX(J1);  
            CHKTY=CHKTY*TY(J1);  
            nps(K)=NPSIDE(J,K); 
            npn(K)=J1;
         end
         
         if CHKTX ~= 0 |  CHKTY ~= 0 
%           ---------------------------------------
%           BEGIN SURFACE QUADRATURE ON SIDE J
%           ---------------------------------------
            nps=[nps,nps+NNPE];
            KEND=NUMQPT(2);
            for K=1:KEND; 
               %----------------------------
               % Calculate quadrature point
               % values of parameters and 
               % mapping dsdu value.
               %----------------------------
               SF4nK=SF(4,1:NNPS,K);
               SF5nK=SF(5,1:NNPS,K);
               XK   = SF4nK*XORD(npn)';
               YK   = SF4nK*YORD(npn)';
               DXDXI= SF5nK*XORD(npn)'; 
               DYDXI= SF5nK*YORD(npn)'; 

               TXK=0.0;
               TYK=0.0;
               if CHKTX ~= 0 
                 TXK  = SF4nK*(TX(npn)'-TY(npn)'); 
               end
               if CHKTY ~= 0 
                 TYK  = SF4nK*(TX(npn)'+TY(npn)');
               end
               Trct=[TXK;TYK];

               %----------------------------------
               % Add quadrature values to integral
               %----------------------------------
               DETJS=sqrt(DXDXI^2+DYDXI^2); 
               Nf=[SF4nK',BLNKs';BLNKs',SF4nK'];
               FE(nps)=FE(nps)+Nf*Trct*WT(2,K)*DETJS;

            end % of surface quadrature 
         end %  if-statement for quadrature
      end % of loop over element sides

%     -------------------------------
%     ROTATE SE and FE TO X'-Y' AXES 
%         (Local Coordinates) 
%     -------------------------------
      for J=1:NNPE 
         TH=THETA(NP(I,J)); 
         if TH ~= 0
            C=cos(TH); 
            S=sin(TH); 
            t = [C,S;-S,C];
            p=[J,J+NNPE]; 
            T=Id;
            T(p,p)=t;
            FE=T*FE;
            SE=T*SE*T';
         end
      end

%     -----------------------------------
%     PLACE SE and FE IN GLOBAL SK MATRIX
%     -----------------------------------
      F(NPE)=F(NPE)+FE; 
      SK(NPE,NPE)=SK(NPE,NPE)+SE;

   end
%  --------  all elements are now assembled 


%  ----------------------------
%  SPECIFY BOUNDARY CONDITIONS
%  ----------------------------
   for I=1:NUMNP
      NBCI=NPBC(I);
      if NBCI == 3 | NBCI == 1
         NPX=2*I-1; 
         SK(NPX,NPX)=SK(NPX,NPX)*1.0E+06; 
         F(NPX) =SK(NPX,NPX)*XBC(I);
      end

     if NBCI == 3 | NBCI == 2
        NPY=2*I; 
        SK(NPY,NPY)=SK(NPY,NPY)*1.0E+06; 
        F(NPY) =SK(NPY,NPY)*YBC(I);
      end
   end


%  ----------------------
%  Solve equations
%  ----------------------
   U = SK\F;

%  ---------------------------- 
%  ROTATE SOLUTION VECTOR TO
%  GLOBAL COORDINATES
%  ---------------------------- 
   for I=1:NUMEQ
      F(I)=U(I); 
   end
   for I=1:NUMNP
      IY =2*I; 
      IX =IY-1;
      TH =THETA(I);
      C  =cos(TH);
      S  =sin(TH);
      U(IX)=F(IX)*C-F(IY)*S;
      U(IY)=F(IX)*S+F(IY)*C;
   end

