C
C NAZ11   TEMP Finite Element Method @Zitny 98
C Temperature field
C
      parameter (num=50,nem=50)
      dimension x(num),y(num),nc(3,nem),ibc(num),tbc(num)
      common am(10000),fm(100)
C Read grid data
      open(1,file='naz11.dat')
      read(1,*)
      read(1,*)nu,ne,nbc
      read(1,*)(x(i),y(i),i=1,nu)
      read(1,*)((nc(j,i),j=1,3),i=1,ne)
      read(1,*)(ibc(i),tbc(i),i=1,nbc)
C Assembly of stiffness matrix and boundary equations
      call gassem(ne,nu,nc,x,y,nbc,ibc,tbc,am,fm)
C Solution of linear equations by Gauss elimination 
      call gaussj(am,nu,nu,fm,1,1)
C Output Results
      open(2,file='naz11.out')
      write(2,1)nu,(x(i),y(i),i=1,nu)
      write(2,2)ne,((nc(j,i),j=1,3),i=1,ne)
      write(2,3)(i,fm(i),i=1,nu)
1     format(i3/(2e10.3))
2     format(i3/(3i5))
3     format(' RESULTS'/(i3,e10.3))
      end

      subroutine gassem(ne,nu,nc,x,y,nbc,ibc,tbc,am,fm)
      dimension x(nu),y(nu),nc(3,ne),am(nu,nu),fm(nu),ibc(nu),tbc(nu),
     /          al(3,3),xl(3),yl(3)
C AM - stiffness matrix of system
C FM - R.H.S. vector of system      
      do i=1,nu
      do j=1,nu
        am(i,j)=0
      enddo
        fm(i)=0
      enddo
C Element loop (NC(1..3, 1..NE)-connectivity matrix)
      do ie=1,ne
         do i=1,3
            iu=nc(i,ie)
            xl(i)=x(iu)
            yl(i)=y(iu)
         enddo
         call kloc(xl,yl,al)
C ASSEMBLY
         do i=1,3
            do j=1,3
               am(nc(i,ie),nc(j,ie))=am(nc(i,ie),nc(j,ie))+al(i,j)
            enddo
         enddo
      enddo
C Boundary conditions
      do i=1,nbc
         ib=ibc(i)
         am(ib,ib)=am(ib,ib)*1e10
         fm(ib)=am(ib,ib)*tbc(i)
      enddo
      end
      
      subroutine kloc(x,y,al)
      dimension x(3),y(3),al(3,3)
C Local STIFFNESS MATRIX AL(3 x 3)
      det=x(1)*(y(2)-y(3))+x(2)*(y(3)-y(1))+x(3)*(y(1)-y(2))
      al=0
      do i=1,3
C New trick: I1=i+1, I2=i+2 in a cyclic list 1,2,3 
         i1=mod(i,3)+1
         i2=mod(i1,3)+1
         do j=1,3
            j1=mod(j,3)+1
            j2=mod(j1,3)+1
            al(i,j)=((y(i1)-y(i2))*(y(j1)-y(j2))+
     /               (x(i2)-x(i1))*(x(j2)-x(j1)))/(2*det)
         enddo
      enddo
      end
      
      subroutine gaussj(a,n,np,b,m,mp)
C
C Numerical Recipes p. 28
C
      parameter (nmax=100)
      dimension a(np,np),b(np,mp),ipiv(nmax),indxr(nmax),indxc(nmax)
      do i=1,n
         ipiv(i)=0
      enddo
      do i=1,n
         big=0
         do j=1,n
             if(ipiv(j).ne.1)then
                   do k=1,n
                         if(ipiv(k).eq.0)then
                               if(abs(a(j,k)).ge.big)then
                                    big=abs(a(j,k))
                                    irow=j
				    icol=k
                               endif
                         elseif(ipiv(k).gt.1)then
                               pause 'Singular matrix'
                         endif
                    enddo
               endif
          enddo
          ipiv(icol)=ipiv(icol)+1
          if(irow.ne.icol)then
              do l=1,n
                 dum=a(irow,l)
                 a(irow,l)=a(icol,l)
                 a(icol,l)=dum
              enddo
              do l=1,m
                 dum=b(irow,l)
                 b(irow,l)=b(icol,l)
                 b(icol,l)=dum
              enddo
          endif
          indxr(i)=irow
          indxc(i)=icol
          if(a(icol,icol).eq.0.) pause 'Singular'
          pivinv=1/a(icol,icol)
          a(icol,icol)=1
          do l=1,n
             a(icol,l)=a(icol,l)*pivinv
          enddo
          do l=1,m
             b(icol,l)=b(icol,l)*pivinv
          enddo
          do ll=1,n
              if(ll.ne.icol)then
                  dum=a(ll,icol)   
                  a(ll,icol)=0
                  do l=1,n
                     a(ll,l)=a(ll,l)-a(icol,l)*dum
                  enddo
                  do l=1,m
                     b(ll,l)=b(ll,l)-b(icol,l)*dum
                  enddo
               endif
           enddo
      enddo
      do l=n,1,-1
         if(indxr(l).ne.indxc(l))then
            do k=1,n
               dum=a(k,indxr(l))
               a(k,indxr(l))=a(k,indxc(l))
               a(k,indxc(l))=dum
            enddo
         endif
      enddo
      end
