program forming
      parameter (nx=8,ny=4,n=nx*ny)
      dimension  p1(nx),u1(nx)
c     appel le code sructure
      call fluide(p1)
c     appel structire
      call  struc(p1,u1)

c     tracer u1
      h=1./(nx-1)
      write(*,*)'h',h
      do i=1,nx
         xi=(i-1)*h
         write(*,*)xi,u1(i)
         write(1,*)xi,u1(i)
       enddo
      stop
      end


      subroutine fluide(p1)
      parameter (nx=8,ny=4,nelems=(nx-1)*(ny-1),n=nx*ny)
      dimension a(n,n),b(n),p(n),p1(nx)
      dimension con(4,nelems)
      integer con
c     initialiser la marice

c     initialiser b


c     connectivite
      do i=1,nx-1
         do j=1,ny-1
            k=(nx-1)*(j-1)+i
            con(1,k)=k+(j-1)
            con(2,k)=con(1,k)+1
            con(3,k)=con(2,k)+nx
            con(4,k)=con(3,k)-1
         enddo
      enddo
c     imprimer con
      do k=1,nelems
         write(*,*)k,(con(j,k),j=1,4)
      enddo
c     assembmage


c     condition aux limites

      p0=1.
      do i=n-nx+1,n
         do j=1,n
            a(i,j)=0.
         enddo
         a(i,i)=1.
         p(i)=-p0
      enddo

c     resoudre A.p = b  par Jacobi

      do i=1,n
         p(i)=-2.
      enddo


c    calculer P1 pression sur les noeuds structire ( 1 a nx)
      do i=1,nx
         p1(i)=p(i)
      enddo
      return
      end

      subroutine struc(p1,u1)
      parameter (n=8,nelems=n-1)
      dimension a(n,n),b(n),ak(2,2),u0(n),u1(n),bk(2),p1(n)
      dimension con(2,nelems)
      integer con
c     initialiser a
      do i=1,n
         do j=1,n
            a(i,j)=0.
         enddo
      enddo
c     initialser b
      do i=1,n
      b(i)=0.
      enddo
c     connectivite
      do k=1,nelems
         con(1,k)=k
         con(2,k)=k+1
      enddo

      h =1./nelems
      write(*,*)'h=',h
c     assemblage
         do k=1,nelems
            ak(1,1)=+1./h
            ak(1,2)=-1./h
            ak(2,1)=-1./h
            ak(2,2)=+1./h
            bk(1)=h*p1(k)/2.
            bk(2)=h*p1(k+1)/2.
            do i=1,2
               i1=con(i,k)
               b(i1)=b(i1)+bk(i)
               do j=1,2
               j1=con(j,k)
               a(i1,j1)=a(i1,j1)+ak(i,j)
            enddo
         enddo
      enddo
c     conditions aux limites
      do j=1,n
         a(1,j)=0.
         a(n,j)=0.
      enddo
      a(1,1)=1.
      a(n,n)=1.
      b(1)=0.
      b(n)=0.
c     imprimer ma matrice
      do i=1,n
         write(*,*)(a(i,j),j=1,n)
      enddo
      do i=1,n
         write(*,*)i,b(i)
      enddo

c     resoudre Au=b
      do i=1,n
         u0(i)=0.
      enddo
      iter=0
 1    do i=1,n
         som1=0.
         do j=1,i-1
            som1=som1+a(i,j)*u0(j)
         enddo
         som2=0.
         do j=i+1,n
            som2=som2+a(i,j)*u0(j)
         enddo
         u1(i)=(b(i)-som1-som2)/a(i,i)
      enddo
c     test
      er=0.
      do i=1,n
         er=er+abs(u1(i)-u0(i))
      enddo
      tol=1.e-5
      if(er.gt.tol) then
         do i=1,n
            u0(i)=u1(i)
         enddo
         iter=iter+1
         write(*,*)iter,er
         goto 1
      endif

      return
      end