%======================
%  PROGRAM Steady.m
%======================

   clear
%-------------------------
%  LOAD DATA
%-------------------------
   load MESHo   -ASCII
   load NODES   -ASCII
   load NP      -ASCII

   NUMNP = MESHo(1);
   NUMEL = MESHo(2);
   NNPE  = MESHo(3);
     
   XORD=NODES(:,1)';
   YORD=NODES(:,2)';
   NPcode=NODES(:,3)';
%  ----------------------- 
%  Declare sparse matrix
%  ----------------------- 
   SK  = sparse(NUMNP,NUMNP);

%  ----------------------
%  General Initialization
%  ----------------------
   PI=4.0*atan(1.0);
   if NNPE == 3
      NSPE=3;
      NNPS=2;
   elseif NNPE == 6
      NSPE=3;
      NNPS=3;
   elseif NNPE == 4
      NSPE=4;
      NNPS=2;
   elseif NNPE == 8
      NSPE=4;
      NNPS=3;
   end

   Q(1:NUMNP,1)=0;
   QS(1:NUMNP)=0;
   HS(1:NUMNP)=0;
   PHI(1:NUMNP)=0;
   PHIa(1:NUMNP)=0;
   NPBC(1:NUMNP)=0;

% --------------------------------------
%  Ask User for type of analysis desired
%  -------------------------------------
   disp(' ')
   disp('ENTER:')
   disp('-----------------------------------------')
   disp('0 for 2D rectangular analysis')
   disp('1 for axi-symmetric analysis about x-axis')
   disp('2 for axi-symmetric analysis about y-axis')
   disp('-----------------------------------------')
   IRZ = input(' < ');

   if IRZ < 0 | IRZ > 2
      disp('--------------------------')
      disp('ERROR IN INPUT ')
      disp('IRZ has been set to 0')
      disp('--------------------------')
      IRZ = 0;
   end

%  ----------------------------------
%  Get shape function quadrature data
%  ----------------------------------
   [SF,WT,NUMQPT,NPSIDE] = SFquad(NNPE);

%  -----------------------------------
%  Include user written initialization
%  -----------------------------------
   INITIAL

%  -----------------
%  Place Q on RHS
%  -----------------
   RHS=Q;

%  -------------------------
%  Create element matrices 
%  -------------------------
   for I=1:NUMEL;
     QE(1:NNPE)=0.0;
     S(1:NNPE,1:NNPE)=0.0;
     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 coefficients
%        RXJ, RYJ, BXJ, BYJ, GVJ, HVJ
%        ---------------------------------
         COEF
%        -------------------------------
%        Create element stiffness matrix
%        -------------------------------
         DV=DETJ; 
         if IRZ == 1  
            DV=2.0*PI*XJ*DV;
         elseif IRZ == 2  
            DV=2.0*PI*YJ*DV;
         end

         DNDu = SF(2,:,J);
         DNDv = SF(3,:,J);
         QE=WT(1,J)*SF(1,:,J)'*HVJ*DV;
         S=S+((DNDX')*RXJ*(DNDX)+(DNDY')*RYJ*(DNDY) ...
            - (SFj')*BXJ*(DNDX) -(SFj')*BYJ*(DNDY) ...
            - (SFj')*GVJ*(SFj))*DV ;
      end
%     ------------ end of volume quadrature

%     ------------------------
%     Begin surface quadrature
%     ------------------------
      for J=1:NSPE
         CHKH=1;
         CHKQ=1;
         for K=1:NNPS;
            nps(K)=NPSIDE(J,K); 
            J1=NP(I,nps(K));
            npn(K)=J1;
            CHKH=CHKH*HS(J1);  
            CHKQ=CHKQ*QS(J1);  
         end

         if CHKH ~= 0 |  CHKQ ~= 0
            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)';
               QSK  = SF4nK*QS(npn)'; 
               HSK  = SF4nK*HS(npn)'; 
               PHIK = SF4nK*PHIa(npn)'; 
               dxdu = SF5nK*XORD(npn)'; 
               dydu = SF5nK*YORD(npn)'; 

               dsdu = sqrt(dxdu^2+dydu^2);
               duds = 1/dsdu;
              
               DS=dsdu;
               if IRZ == 1  
                  DS=2.0*PI*YK*DS;
               elseif IRZ == 2
                  DS=2.0*PI*XK*DS; 
               end

               %----------------------------------
               % Add quadrature values to integral
               %----------------------------------
               if CHKQ ~= 0
                 QE(nps)=QE(nps)+WT(2,K)*SF4nK'*QSK*DS;
               end
               if CHKH ~= 0
                 QE(nps)=QE(nps)+WT(2,K)*SF4nK'*HSK*PHIK*DS;
                 S(nps,nps)=S(nps,nps)+WT(2,K)*SF4nK'*HSK*SF4nK*DS;
               end

            end % of surface quadrature 

         end %  if-statement for quadrature

      end % of loop over element sides

%     ------------------------------------------------------------
%     Place completed element matrix in global SK and Q matrices
%     ------------------------------------------------------------
      RHS(rc)=RHS(rc)+QE;
      SK(rc,rc)=SK(rc,rc)+S;

   end  % Loop over elements

%  ---------------------------
%  Specify known values of PHI 
%  ---------------------------
   for I=1:NUMNP;
      if NPBC(I) == 1
         SK(I,I)=SK(I,I)*1.0E+06; 
         RHS(I)=PHI(I)*SK(I,I); 
      end
   end

%  ----------------------
%  Solve equations
%  ----------------------
   PHI = SK\RHS;
   save PHI PHI -ASCII


