!
!---- Master MIMSE, 
! 
! code reformulé  par L. Weynans à partir du code F77:
!     (Universite de Savoie -- R. Hassani - Sept. 2000)
!
!     Premier programme d'initiation a la methode des elements finis
!
!     Probleme du laplacien :
!       -Laplacien(u) = f  dans W
!       u = g sur G
!
!     G etant la frontiere (ou une partie) de W
!
!     discretise avec des elements finis P1
!
!     La matrice du systeme est stockee dans un tableau bidimensionel
!     et le systeme lineaire est resolu par la methode de Cholevsky
!     
!     Details des tableaux et des variables :
!     -------------------------------------
!
!     variables entieres :
!     ------------------
!     - npmax : nbr maximum de noeuds dans le maillage
!     - nemax : nbr maximum d'elements (triangles) dans le maillage
!     - npoin : nbr effectif de noeuds ( < npmax)
!     - nelem : nbr effectif d'elements ( < nemax)
!     - ndiri : nbr de noeuds sur la frontiere
!
!
!     tableaux :
!     ---------
!     1) concernant le maillage :
!     - coord : Tableau de reels de dimension 2 x npmax (dimension
!               utile : 2 x npoin) contenant les coordonnees
!               des noeuds.
!               Pour un noeud de no. ip (ip = 1, ..., npoin) :
!                   x = coord(1,ip)
!                   y = coord(2,ip)
!     - lnods : Tableau d'entiers de dimension 3 x nemax (dimension
!               utile : 3 x nelem) contenant les connectivites des elements.
!               Pour un element de no. ie (ie = 1, ..., nelem) :
!                   ip1 = lnods(1,ie)
!                   ip2 = lnods(2,ie)
!                   ip3 = lnods(3,ie)
!               sont les no. des 3 noeuds sommets
!
!     2) concernant le systeme lineaire
!     - rigel : Tableau de reels de dimension 3x3 formant la
!               matrice de rigidite elementaire
!     - vecel : Tableau de reels de dimension 3 formant le
!               vecteur second membre elementaire
!     - rigid : Tableau de reels de dimension npmax x npmax (dimension
!               utile : npoin x npoin) formant la matrice de rigidite
!     - sdmem : Tableau de reels de dimension npmax (dimension utile : npoin)
!               formant le second membre
!     - solut : Tableau de reels de dimension npmax (dimension utile : npoin)
!               formant la solution du systeme (valeurs approchee de l'inconnue
!               aux npoin noeuds du maillage)
!     - trian : Tableau de reels de dimension npmax x npmax (dimension utile :
!               npoin x npoin) resultant de la factorisation de Cholevsky
!               de la matrice rigid 
!
!     3) concernant les conditions aux limites de Dirichlet
!     - ldiri : Tableau d'entiers de dimension npmax (dimension utile : ndiri)
!               donnant la liste des no. des noeuds de la frontiere G ou est
!               donnee la valeur de l'inconnue
!
!     Details des fonctions et routines :
!     ---------------------------------
!
!     - maillage : calcul de npoin, nelem, ndiri et creation des tables 
!                  coord, lnods et ldiri pour le carre unite 
!     - rigidel  : calcul pour un element P1 de la matrice (rigel) et du 
!                  vecteur (vecel) elementaires
!     - assemb   : assemblages d'une matrice et d'un vecteur elementaires dans
!                  la matrice globale (rigid) et le vecteur global (sdmem)
!     - dirich   : modification de la matrice et du second membre par prise
!                  en compte des conditions aux limites (methode du terme
!                  diagonal dominant)
!     - decchol  : decomposition de Cholevsky d'une matrice sym. def. pos.
!     - reschol  : resolution d'un systeme lineaire apres factorisation
!                  de Cholevsky de sa matrice 
!     - f        : fonction f donnee analytiquement
!     - g        : fonction g donnee analytiquement
!     
!
   


!*************************************************************

module fonctions

implicit none

contains

      subroutine maillage (coord,lnods,ldiri,npoin, nelem,npmax,nemax,ndiri,n)

!---- maillage regulier du carre unite par des triangles

        real*8 :: h,x,y
        integer :: npoin, nelem, npmax, nemax, ndiri,n,ie,ip
        integer:: i,j, ip1,ip2,ip3,ip4,ib,ir,it,il,id
        real*8, dimension (2,npmax) :: coord
        integer,  dimension(3,nemax)::  lnods
        integer,  dimension(npmax) :: ldiri

      print*,'nbr de points selon x et y'
      read(*,*)n

      h = 1.0 / float(n-1)
      npoin = n * n
      nelem = 2 * (n - 1) * (n - 1)

      if (npoin .gt. npmax .or. nelem .gt. nemax) then
         print*,'maillage trop fin'
         print*,'nelem =',nelem, ' nemax =',nemax
         print*,'npoin =',npoin, ' npmax =',npmax
         stop
      endif

!---- creation des coordonnees et des connectivites

      ie = 0
      ip = 0
      do j = 1, n
         y = (j-1)*h
         do i = 1, n
            x = (i-1)*h
            ip = ip + 1
            coord(1,ip) = x
            coord(2,ip) = y
            if (i .ne. n .and. j .ne. n) then
               ip1 = ip
               ip2 = ip1 + 1
               ip3 = ip1 + n
               ip4 = ip2 + n
               ie = ie + 1
               lnods(1,ie) = ip1
               lnods(2,ie) = ip2
               lnods(3,ie) = ip3
               ie = ie + 1
               lnods(1,ie) = ip4
               lnods(2,ie) = ip3
               lnods(3,ie) = ip2
            endif
         enddo
      enddo

!---- memorisation des noeuds de la frontiere a condition imposee

      print*,'indiquer par une suite de 0 ou de 1 si les faces'
      print*,'sont des faces a valeurs imposees, en commencant'
      print*,'par celle du bas et en tournant dans le sens direct'
      read(*,*)ib,ir,it,il

      id = 0
      if (ib .eq. 1) then
         do i = 1, n
            id = id + 1
            ldiri(id) = i
         enddo
      endif

      if (it .eq. 1) then
         do i = n*(n-1) + 1, n*n
            id = id + 1
            ldiri(id) = i
         enddo
      endif

      if (ir .eq. 1) then
         do i = 2, n-1
            id = id + 1
            ip = n*i
            ldiri(id) = ip
         enddo
      endif

      if (il .eq. 1) then
         do i = 2, n-1
            id = id + 1
            ip = n*(i-1) + 1
            ldiri(id) = ip
         enddo
      endif
      ndiri = id

      if (ndiri .eq. 0) then
         print*,'Attention : aucun d.d.l. impose'
         stop
      endif

     
    end subroutine maillage

!*************************************************************

      subroutine rigidel (rigel,vecel,lnods,coord, npmax,nemax,ie   )

!---- Calcul de la matrice de rigidite elementaire et du vecteur
!     elementaire pour un element P1 de numero ie

!     Les coefficients de la matrice elementaire : 
!
!             ke_ij = integrale (gradN_i . gradN_j)
!
!     sont simples a calculer puisque les fonctions d'interpolation N_i
!     (i = 1, 2, 3) sont affines sur l'element (gradN_i est constant) :
!
!             ke_ij = s * gradN_i . gradN_j
!
!     s etant l'aire du triangle
!
!     Les composantes du vecteur elementaire :
!
!            fe_i = integrale (N_i * f)
!
!     peuvent etre calculees de differentes facons suivant la fonction f
!     (integration exacte, quadratures de Gauss ou autres ..). 
!     En particulier :
!
!     Si f est supposee constante sur l'element:
!          fe_i = (s / 3) * f
!
!     Si f ne l'est pas on peut toujours l'approcher par :
!          fe_i = (s / 3) * fm  
!          ou fm est soit la valeur au centre du triangle soit la valeur
!          moyenne des 3 valeurs aux sommets : fm = (f(i) + f(j) + f(k))/3
!          Dans ce dernier cas, la formule est d'ordre h. 
!
!     Si f est affine sur l'element on peut aussi utiliser la formule :
!          fe_i = (s / 12) * (2*f(i) + f(j) + f(k))
!          qui est exacte. Si f n'est pas affine on obtient une approximation
!          d'ordre h2.
!               
!      en entree :
!       - coord : coordonnees des noeuds
!       - lnods : table de connectivite
!       - ie    : numero de l'element concerne
!
!      en sortie :
!       - rigel : matrice de rigidite elementaire (ke)
!       - vecel : vecteur elementaire (fe)
     


        real*8 :: h
        integer :: npoin, nelem, npmax, nemax, ndiri,n,i,j,ie
        real*8, dimension (2,npmax) :: coord
        integer,  dimension(3,nemax)::  lnods
        integer,  dimension(npmax) :: ldiri
        real*8, dimension (3) :: vecel
        real*8, dimension (3,3) :: rigel
        real*8 :: x1,x2,x3,y1,y2,y3,x12,x23,x31,y12,y23,y31,s,c
        real*8 ::f1,f2,f3
    
     
!---- coordonnees des 3 noeuds de l'element

      x1 = coord(1,lnods(1,ie))
      y1 = coord(2,lnods(1,ie))
      x2 = coord(1,lnods(2,ie))
      y2 = coord(2,lnods(2,ie))
      x3 = coord(1,lnods(3,ie))
      y3 = coord(2,lnods(3,ie))

!---- calcul de l'aire de l'element

      x12 = x2 - x1
      y12 = y2 - y1
      x23 = x3 - x2
      y23 = y3 - y2
      x31 = x1 - x3
      y31 = y1 - y3
      s = 0.5*(y12*x31 - x12*y31)

      if (s .le. 0.0) then
         print*,'erreur : element a surface negative ou nulle'
         print*,'element no. ',ie
         print*,'coordonnees de ses 3 sommets :'
         print*,x1,y1
         print*,x2,y2
         print*,x3,y3
         stop
      endif

!---- calcul des coefficients de la matrice elementaire

      c = 0.25 / s
      rigel(1,1) = c*(y23*y23 + x23*x23)
      rigel(1,2) = c*(y23*y31 + x23*x31)
      rigel(1,3) = c*(y23*y12 + x23*x12)
      rigel(2,2) = c*(y31*y31 + x31*x31)
      rigel(2,3) = c*(y31*y12 + x31*x12)
      rigel(3,3) = c*(y12*y12 + x12*x12)
      rigel(2,1) = rigel(1,2)
      rigel(3,1) = rigel(1,3)
      rigel(3,2) = rigel(2,3)

!---- calcul des composantes du vecteur elementaire
!     (formule de quadrature obtenue en supposant que 
!      f est affine sur l'element)

      c = s / 12.0

!---- calcul de f aux 3 sommets

      f1 = f(x1,y1)
      f2 = f(x2,y2)
      f3 = f(x3,y3)

!---- integration

      vecel(1) = c * (2.*f1 + f2 + f3)
      vecel(2) = c * (2.*f2 + f3 + f1)
      vecel(3) = c * (2.*f3 + f1 + f2)

      
    end subroutine rigidel

!*************************************************************

      subroutine assemb (rigel,vecel,rigid,sdmem,lnods,npmax,nemax,ie   )

!---- Assemblage de la matrice elementaire de l'element ie, rigel, 
!     dans la matrice de rigidite globale, rigid :
!
!           K(ip,jp) = K(ip,jp) + ke(i,j)
!
!     avec ip = lnods(i,ie) et jp = lnods(j,ie)
!
!     Assemblage du vecteur elementaire de l'element ie, vecel, 
!     dans le vecteur second membre global, sdmem
!
!          F(ip) = fe(i)
!
!     avec ip = lnods(i,ie)
!
!
!     en entree :
!       - rigel : matrice de rigidite elementaire (ke)
!       - vecel : vecteur elementaire (fe)
!       - rigid : matrice de rigidite globale (K)
!       - sdmem : vecteur second membre global (F)
!       - lnods : table de connectivite
!       - ie    : numero de l'element concerne
!
!     en sortie :
!       - rigid : matrice de rigidite globale actualisee
!       - sdmem : vecteur second membre global actualise

        

        real*8 :: h
        integer :: npoin, nelem, npmax, nemax, ndiri,n,i,j,in,ip,jn,ie,jp
        real*8, dimension (2,npmax) :: coord
        integer,  dimension(3,nemax)::  lnods
        integer,  dimension(npmax) :: ldiri
        real*8, dimension (3) :: vecel
        real*8, dimension (3,3) :: rigel
        real*8 :: x1,x2,x3,y1,y2,y3,x12,x23,x31,y12,y23,y31,s,c
        real*8 ::f1,f2,f3
        real*8, dimension (npmax,npmax) :: rigid
        real*8, dimension (npmax) :: sdmem
    
    
!---- boucle sur les 3 sommets de l'element

      do in = 1, 3

!----    numero du noeud du sommet in

         ip = lnods(in,ie)

!----    sommation de la ligne in de la matrice elementaire
!        dans la matrice globale

         do jn = 1, 3
            jp = lnods(jn,ie)
            rigid(ip,jp) = rigid(ip,jp) + rigel(in,jn)
         enddo

!----    sommation de la composante in du vecteur elementaire
!        dans le vecteur global

         sdmem(ip) = sdmem(ip) + vecel(in)
      enddo

   
    end subroutine assemb

!*************************************************************

      subroutine dirich (rigid,sdmem,coord,ldiri, npmax,ndiri)

!---- Prise en compte des conditions de Dirichlet u = g sur G.
!     La liste des no. des noeuds formant G est donnee par ldiri.
!     
!     Methode du "terme diagonale dominant" :
!
!     pour un noeud ip de G on modifie : 
!     - le coefficient diagonal correspondant de la matrice de rigidite : 
!                 K(ip,ip) = K(ip,ip) + C
!     - la composante correspondante du second membre :
!                 F(ip) = C * g
!     
!     ou C est une constante "tres grande"
!
!
    
        integer ::  npmax,  ndiri,i,ip
        real*8, dimension (2,npmax) :: coord
        integer,  dimension(npmax) :: ldiri
        real*8 :: grand
        real*8, dimension (npmax,npmax) :: rigid
        real*8, dimension (npmax) :: sdmem

   
      grand = 1.e20

!---- boucle sur les noeuds de G

      do i = 1, ndiri

!----    numero du noeud

         ip = ldiri(i)

!----    modification de la matrice de rigidite

         rigid(ip,ip) = rigid(ip,ip) + grand

!----    modification du second membre

         sdmem(ip) = grand * g(coord(1,ip),coord(2,ip))
      enddo

    end subroutine dirich

!*************************************************************

     function f (x,y)

       real*8 :: f,a,x,y

!---- Fonction utilisateur definissant analytiquement la fonction f
!     donnant la valeur du laplacien de la fonction inconnue

  
      a = 6.28                        ! la solution de -Laplacien(u)=f
      f = 2.0*a*a*sin(a*x)*cos(a*y)   ! est alors u = sin(ax)*cos(ay)
      f = 4.0

    end function f

!*************************************************************

      function g (x,y)
!---- Fonction utilisateur definissant analytiquement la fonction g
!     donnant la valeur de l'inconnue sur la frontiere
     

        real*8 :: a,g,x,y

      a = 6.28
      g = sin(a*x)*cos(a*y)
      g = (x-0.5)**2 + (y-0.5)**2
     g = 0.

    end function g

!*****************************************************************

      subroutine decchol (A,R,n,nmax)

!---- decomposition de Cholevsky 
!
!                    t
!             A = R.R
!
!     d'une matrice A(n,n) symetrique definie positive

!---- dimensionnement de A et R avec les dimensions reelles
!     de ces tableaux (tels qu'ils ont ete declare dans le
!     programme appelant)

      real*8, dimension(nmax,nmax):: A,R
      integer :: n, nmax,j,k,i
      real*8 :: c,rjj

      do j = 1, n
         c = 0.0
         do k = 1, j
            c = c + R(j,k)*R(j,k)
         enddo
         R(j,j) = sqrt(A(j,j) - c)
         rjj = 1.0 / R(j,j)
         do i = j+1, n
            c = 0.0
            do k = 1, j
               c = c + R(i,k)*R(j,k)
            enddo
            R(i,j) = (A(i,j) - c)*rjj
         enddo
      enddo

    end subroutine decchol

!*****************************************************************

      subroutine reschol (R,x,b,n,nmax)

!                                             t
!---- Resolution du systeme lineaire nxn : R.R x = b
!     ou R est une matrice n x n triangulaire inferieure
!

!---- dimensionnement de R avec les dimensions reelles de ce tableaux

      real*8, dimension(nmax,nmax):: R
      real*8 , dimension(nmax) :: x, b
      integer :: i, j,n,nmax
      real*8 :: c

!---- 1) descente : resolution de R z = b (z est range dans b)

      do i = 1, n
         c = 0.0
         do j = 1, i-1
            c = c + R(i,j)*b(j)
         enddo
         b(i) = (b(i) - c) / R(i,i)
      enddo

!                                t
!---- 2) montee : resolution de R x = z      

      do i = n, 1, -1
         c = 0.0
         do j = i+1, n
            c = c + R(j,i)*x(j)
         enddo
         x(i) = (b(i) - c) / R(i,i)
       enddo

     end subroutine reschol

   end module fonctions


   program Poisson 
     
     use fonctions
     
     implicit none
     
     !---- nbr max de noeuds et d'elements
     
     
     !---- dimensionnement des divers tableaux
     
     integer, parameter :: npmax=1000, nemax= 2000
     integer :: npoin, nelem,  ndiri,n,i,j,in,ip,jn,ie,jp,k
     real*8, dimension (2,npmax) :: coord
     integer,  dimension(3,nemax)::  lnods
     integer,  dimension(npmax) :: ldiri
     real*8, dimension (3) :: vecel
     real*8, dimension (3,3) :: rigel
     real*8 :: x1,x2,x3,y1,y2,y3,x12,x23,x31,y12,y23,y31,s,c
     real*8 ::f1,f2,f3
     real*8, dimension (npmax,npmax) :: rigid,trian
     real*8, dimension (npmax) :: sdmem,solut
    

  
!---- creation d'un maillage regulier pour le carre unite

      call maillage (coord,lnods,ldiri,npoin,nelem,npmax,nemax,ndiri,n)

!---- mise a zero de la matrice de rigidite (rigid), 
!     du vecteur second membre (sdmem) et du vecteur solution (solut)

      do ip = 1, npoin
         do jp = 1, npoin
            rigid(ip,jp) = 0.0
         enddo
         sdmem(ip) = 0.0
         solut(ip) = 0.0
      enddo

!---- calcul de la matrice de rigidite et du second membre

      do ie = 1, nelem

!----    matrice de rigidite et vecteur elementaires

         call rigidel (rigel,vecel,lnods,coord, npmax,nemax,ie)

!----    assemblage de rigid et de sdmem

         call assemb (rigel,vecel,rigid,sdmem,lnods,  npmax,nemax,ie)
      enddo

!---- prise en compte des conditions aux limites 
!     (modification de rigid et sdmem)

      call dirich (rigid,sdmem,coord,ldiri,npmax,ndiri)

!---- resolution par factorisation de cholevsky

      call decchol(rigid,trian,npoin,npmax)
      call reschol(trian,solut,sdmem,npoin,npmax)

!---- ecriture de la solution dans un fichier

      open(1,file='solut')
    !  do ip = 1, npoin
    !     write(1,*)(coord(i,ip),i=1,2),solut(ip)
    !  enddo 

      ip = 0
      do j = 1,n
         do i = 1,n
            ip = ip + 1
            write(1,*)(coord(k,ip),k=1,2),solut(ip)   
            print*, i,j, ip         
         enddo
         write(1,*)
      enddo

          

    end program Poisson
