% =================     
% ROUTINE Stress.m
% =================     
   NNPE2=2*NNPE;
   NUMEQ=NUMNP;

   SIGXX(1:NUMNP)=0.0;
   SIGXY(1:NUMNP)=0.0;
   SIGYY(1:NUMNP)=0.0;
   SIGZZ(1:NUMNP)=0.0;
   SGEFF(1:NUMNP)=0.0;
   SKD(1:NUMNP)=0.0; 

   for I=1:NUMEL
%     ------------------------------ 
%     Begin calculation of strain
%     and stress at each quadrature 
%     point in each element.
%     ------------------------------ 
      npx=2*NP(I,:)-1;
      npy=2*NP(I,:);
      rc=NP(I,:);
      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 stress.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;
%        ---------------------------------
%        Determine R-matrix corresponding
%        to analysis being performed
%        ---------------------------------
%        IPROG = 1    PLANE STRAIN ELASTICITY
%        IPROG = 2    PLANE STRESS ELASTICITY
%        ------------------------------------
         COEF
         if IPROG == 1        
            % Plane strain
            R(1,1)=2.0*G+LAMBDA;
            R(2,2)=R(1,1); 
            R(3,3)=G;
            R(4,4)=G;
            R(1,2)=LAMBDA;
            R(2,1)=LAMBDA;
            Mu = LAMBDA/(2*(LAMBDA+G));
         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)=G;
            R(4,4)=G;
            Mu = LAMBDA/(2*(LAMBDA+G));
         end
%       --------------------------------
%       CALCULATE CONTRIBUTION TO RHS 
%       FOR STRESS
%       --------------------------------
        EPXX=DNDX*U(npx);
        EPYY=DNDY*U(npy);
        EPXY=(DNDY*U(npx)+DNDX*U(npy))/2.;

        SGXX=R(1,1)*EPXX+R(1,2)*EPYY;
        SGYY=R(2,1)*EPXX+R(2,2)*EPYY;
        SGXY=2.0*R(3,3)*EPXY;

        if IPROG == 1
           SGZZ=Mu*(SGXX+SGYY);
        else
           SGZZ=0;
        end

        for K=1:NNPE 
           NPK=NP(I,K);
           XK=XORD(NPK);
           YK=YORD(NPK);
           wght=(XK-XJ)^2 + (YK-YJ)^2;
           aa=wght;
           wght=1.0/(sqrt(wght));
           SIGXX(NPK)=SIGXX(NPK)+wght*SGXX;
           SIGYY(NPK)=SIGYY(NPK)+wght*SGYY;
           SIGXY(NPK)=SIGXY(NPK)+wght*SGXY;
           SIGZZ(NPK)=SIGZZ(NPK)+wght*SGZZ;
           SKD(NPK)=SKD(NPK)+wght;
        end
     end
%    ---------   end of volume quadrature
   end
%  ---------   end of element loop


% ----------------------------------
% Solve for nodal values of stress
% components by dividing by diagonal
% ----------------------------------
   for I=1:NUMNP
     SIGXX(I)=SIGXX(I)/SKD(I);
     SIGYY(I)=SIGYY(I)/SKD(I);
     SIGXY(I)=SIGXY(I)/SKD(I);
     SIGZZ(I)=SIGZZ(I)/SKD(I);
   end

% --------------------------------
%  Calculate nodal point value of 
%  effective stress 
% --------------------------------
   for I=1:NUMNP
      XX=SIGXX(I);
      YY=SIGYY(I);
      XY=SIGXY(I);
      ZZ=SIGZZ(I);
      XP=XORD(I);
      YP=YORD(I);
 
      SG=((XX-YY)^2 +(YY-ZZ)^2 ...
         +(ZZ-XX)^2 +6.0*XY^2);

      SGEFF(I)=sqrt(SG/2.0);
   end




