%------------------------------------* 
%                                    * 
%          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

%    -----------------------------
%    General initialization of all
%    arrays and parameters
%    ------------------------------
     NUMNP = MESHo(1);
     NUMEL = MESHo(2);
     NNPE  = MESHo(3);
      
     if NNPE ~= 3 & NNPE ~= 4
       error(...
       'NNPE must equal 3 or 4 in transL')
     end

     for I=1:NUMNP;
       XORD(I)=NODES(I,1);
       YORD(I)=NODES(I,2);
       NPcode(I)=NODES(I,3);
       NPBC(I)=0.0;
     end

     NNPS=2;
     if NNPE     == 3
       NSPE=3;
     elseif NNPE == 4
       NSPE=4;
     end

     PI=3.1415926;

%    -------------------------------------
%    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
       for I=1:NUMNP
         Q(I) =0.0;
         CL(I)=0.0;
       end
%      -------------------------------
%      INCLUDE user written BC.m
%      to set all boundary condtions
%      -------------------------------
       BC
%      -----------------------------
%      Formation of C and S matrices
%      -----------------------------
       for I=1:NUMEL
         LMNT = I;
         for J=1:NNPE 
           QE(J)=0.0; 
           C(J)=0.0;
           for K=1:NNPE 
             S(J,K)=0.0;
           end
         end
%        ----------------------- 
%        Begin volume quadrature 
%        ----------------------- 
         JEND=NUMQPT(1);
         for J=1:JEND 
           RJAC(1,1)=0.0; 
           RJAC(1,2)=0.0; 
           RJAC(2,1)=0.0; 
           RJAC(2,2)=0.0; 

           XJ=0.0;
           YJ=0.0;
           PHIJ=0.0;

           for K=1:NNPE 
             NPK=NP(I,K); 
             RJAC(1,1)=RJAC(1,1)+SF(2,K,J)*XORD(NPK); 
             RJAC(1,2)=RJAC(1,2)+SF(3,K,J)*XORD(NPK); 
             RJAC(2,1)=RJAC(2,1)+SF(2,K,J)*YORD(NPK); 
             RJAC(2,2)=RJAC(2,2)+SF(3,K,J)*YORD(NPK); 
             XJ=XJ+SF(1,K,J)*XORD(NPK); 
             YJ=YJ+SF(1,K,J)*YORD(NPK); 
             PHIJ=PHIJ+SF(1,K,J)*PHI(NPK); 
           end

           DETJ=RJAC(1,1)*RJAC(2,2)...
               -RJAC(2,1)*RJAC(1,2);
           if DETJ <= 0
             error(' DETJ <= 0 in TRANSL.m')
           end

           DV=DETJ; 
           if IRZ == 1
             DV=2.0*PI*XJ*DV;
           elseif IRZ == 2
             DV=2.0*PI*YJ*DV;
           end

           RJACI(1,1)=+RJAC(2,2)/DETJ;
           RJACI(1,2)=-RJAC(1,2)/DETJ;
           RJACI(2,1)=-RJAC(2,1)/DETJ;
           RJACI(2,2)=+RJAC(1,1)/DETJ;

           for K=1:NNPE 
             DNDX(K)=RJACI(1,1)*SF(2,K,J)...
                    +RJACI(2,1)*SF(3,K,J); 
             DNDY(K)=RJACI(1,2)*SF(2,K,J)...
                    +RJACI(2,2)*SF(3,K,J); 
           end

%          -------------------------
%          INCLUDE user's COEF.m
%          -------------------------
           COEF

%          -------------------------------
%          Create element matrices
%          -------------------------------
           for K=1:NNPE 
             C(K)=C(K) + ...
                WT(1,J)*SF(1,K,J)*RC*DV;
             QE(K)=QE(K)+...
                WT(1,J)*SF(1,K,J)*QV*DV;
             for L=1:NNPE 
               S(K,L)=S(K,L)+WT(1,J)*...
                   (DNDX(K)*RX*DNDX(L) ...
                   +DNDY(K)*RY*DNDY(L))*DV;
              end
           end

         end
%        --------  end of volume quadrature

%        ------------------------ 
%        Begin surface quadrature 
%        ------------------------ 
         for J=1:NSPE
           Hchk=1;
           Qchk=1;
           for K=1:NNPS;
             J1=NP(I,NPSIDE(J,K));
             Hchk=Hchk*HS(J1);  
             Qchk=Qchk*QS(J1);  
           end
             
           if Hchk ~= 0 |  Qchk ~= 0
             KEND=NUMQPT(2); 
             for K=1:KEND
               XK   =0.0;  
               YK   =0.0;  
               DXDXI=0.0;  
               DYDXI=0.0;  
               HSK  =0.0;  
               QSK  =0.0;  
               PHIK =0.0;  
               for L=1:NNPS
                 L1=NPSIDE(J,L);
                 NPL=NP(I,L1); 
                 XL=XORD(NPL);
                 YL=YORD(NPL);
                 QSL=QS(NPL);
                 HSL=HS(NPL);
                 PHIL=PHIa(NPL);
   
                 XK   =XK   +SF(4,L,K)*XL;
                 YK   =YK   +SF(4,L,K)*YL;
                 DXDXI=DXDXI+SF(5,L,K)*XL;
                 DYDXI=DYDXI+SF(5,L,K)*YL;
                 QSK  =QSK  +SF(4,L,K)*QSL;
                 HSK  =HSK  +SF(4,L,K)*HSL;
                 PHIK =PHIK +SF(4,L,K)*PHIL;
               end

               DS=sqrt(DXDXI^2+DYDXI^2); 
               if IRZ == 1
                 DS=2.0*PI*XK*DS; 
               end
               if IRZ == 2
                 DS=2.0*PI*YK*DS; 
               end

               for L=1:NNPS
                 L1=NPSIDE(J,L);
                 if Qchk ~= 0;
                   QE(L1)=QE(L1)+WT(2,K)*SF(4,L,K)*QSK*DS; 
                 end
                 if Hchk ~= 0;
                   QE(L1)=QE(L1)+WT(2,K)*SF(4,L,K)*HSK*PHIK*DS;
                   for M=1:NNPS
                     M1=NPSIDE(J,M);
                     S(L1,M1)=S(L1,M1)...
                       +WT(2,K)*SF(4,L,K)*HSK*SF(4,M,K)*DS; 
                   end
                 end
               end

             end % surface quadrature

           end % Hchk and Qchk

         end % search through element sides

%        ----------------------
%        Create global vectors
%        ----------------------
         for J=1:NNPE 
           NPJ=NP(I,J);
           Q(NPJ)=Q(NPJ)+QE(J); 
           CL(NPJ)=CL(NPJ)+C(J); 
           for K=1:NNPE 
             NPK=NP(I,K); 
             Q(NPJ)=Q(NPJ)-S(J,K)*PHI(NPK);
           end
         end

       end % loop through elements

%      ---------------------------------
%      Increment variables to next time
%      ---------------------------------
       for I=1:NUMNP
         dPdt(I) = 0;
         if NPBC(I) ~= 1
           dPdt(I) = Q(I)/CL(I);
           PHI(I)=PHI(I)+dPdt(I)*DTIME;
         end
       end
   
%      -------------------------
%      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;
         for I=1:NUMNP
           PHIp(ISEG,I)=PHI(I);
         end

         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(' ------------------------------------------------------ ')


