%------------------------------------* 
%          PROGRAM TransL            *
%                                    * 
%        User must supply:           *
%------------------------------------* 
% INITIAL.m  BC.m  COEF.m  MONITOR.m *
%------------------------------------* 

     clear
%    -------------------------
%    LOAD data from mesh
%    -------------------------
     load MESHo   -ASCII
     load NODES   -ASCII
     load NP      -ASCII

     NUMNP = MESHo(1);
     NUMEL = MESHo(2);
     NNPE  = MESHo(3);
     if NNPE ~= 3 & NNPE ~= 4
       error(...
       'NNPE must equal 3 or 4 in transL')
     end

%    ----------------------- 
%    Initialize variables
%    ----------------------- 
     XORD=NODES(:,1)';
     YORD=NODES(:,2)';
     NPcode=NODES(:,3)';
     NNPS=2;
     if NNPE  == 3
       NSPE=3;
     elseif NNPE == 4
       NSPE=4;
     end

%    -------------------------------------
%    Create SF array
%    -------------------------------------
     [SF,WT,NUMQPT,NPSIDE] = SFquad(NNPE);

%    ---------------------------------
%    INCLUDE user's initialization
%    ---------------------------------
     INITIAL

%    ------------------------
%    Begin transient analysis
%    ------------------------
     disp(' ')
     disp('Segment    Time      MinP        MaxP  ')
     disp('-------------------------------------- ')

     TIME=0.0;
     INCR=0;
     ISEG=0;

     while  INCR <  MNI
       Q(1:NUMNP) =0.0;
       CL(1:NUMNP)=0.0;
%      -------------------------------
%      INCLUDE user written BC.m
%      to set all boundary condtions
%      -------------------------------
       BC
%      -----------------------------
%      Formation of C and S matrices
%      -----------------------------
       for I=1:NUMEL
         QE(1:NNPE)=0.0;
         C(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,1:NNPE,J);
            DNDu = SF(2,1:NNPE,J);
            DNDv = SF(3,1:NNPE,J);

%           ---------------------------------
%           Determine coordinate and Jacobian
%           ---------------------------------
            SFj=SF(1,1:NNPE,J);
            XJ=SFj*XORD(rc)';
            YJ=SFj*YORD(rc)';
            PHIJ=SFj*PHI(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

            C= C+WT(1,J)*SFj*RC*DV;
            QE=QE+WT(1,J)*SFj*QV*DV;
            S=S+((DNDX')*RX*(DNDX)+(DNDY')*RY*(DNDY))*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

%        ------------------------------
%        Assemble into global matrices
%        ------------------------------
         Q(rc)=Q(rc)+QE;
         CL(rc)=CL(rc)+C;
         Q(rc)=Q(rc)-(S*PHI(rc)')';

     end % of loop over all elements
%    ---------------------------------
%    Increment variables to next time
%    ---------------------------------
     for I=1:NUMNP
       dPdt(I)=Q(I)/CL(I);
       if NPBC(I) == 1
         dPdt(I) = 0;
       end
     end
   
     PHI=PHI+dPdt*DTIME;
%    -------------------------
%    INCLUDE user's MONITOR.m
%    -------------------------
     MONITOR

     TIME=TIME+DTIME; 
     INCR=INCR+1; 
     if INCR >= MNI
       INCsv=INCR;
     end

     if INCR >= INCsv
       INCsv=INCsv+INTsv; 
       ISEG=ISEG+1;
       PHIp(ISEG,:)=PHI;
       minP = min(PHI);
       maxP = max(PHI);
       a = [ISEG TIME minP maxP];
       fprintf(1,'%4i  %10.2e %10.2e %10.2e \n',a)
     end

    end % of transient analysis

    save PHIp    PHIp  -ASCII
    save PHIstrt PHI   -ASCII
   
    fprintf(1,'\n NUMNP = %4i \n',NUMNP)
 
    disp(' ------------------------------------------------------ ')
    disp(' Nodal point values for the dependent variable, PHI,    ')
    disp(' corresponding to the above segments at the times shown ')
    disp(' are recorded on file PHIp in in your working directory.')
    disp(' ------------------------------------------------------ ')
    disp(' To extract any segement, e.g. 5, for plotting or other ')
    disp(' uses, type:                                            ')
    disp('                   load PHIp                            ')
    disp('                   PHI = PHIp(5,1:N)                    ')
    disp('                   save PHI PHI -ASCII                  ')
    disp(' where N is the numerical value for NUMNP shown above.  ')
    disp(' ------------------------------------------------------ ')
