 SUBROUTINE init_mesh (msh)
!=========================== DOCUMENTATION ============================= 
!DOCDEBh*    UTI/PRE/init_mesh  *
!
!    NOM
! init_mesh (msh)
!    DESCRIPTION
!*
!     Calcul les donnees geometriques (surfaces, volumes, normales)
!     a partir de donnees elementaires (liste des points et connectivite faces-cellules)  
!     Calculate geometric data (surfaces, volumes, normals)
!     from basic data (list of points and face-cell connectivity)  
!*
! 
! 
!    ENTREES / INPUT
!    msh: le maillage / the mesh
!
!    SORTIES / OUTPUT
!    msh: le maillage / the mesh
!
!    SOURCE
! Ce logiciel est regi par la licence [CeCILL-B]
! This software is governed by the [CeCILL-B] license
!=========================== DEBUT DES DECLARATIONS ==================== 
!.1-----  Implicit, Use
  use geo_typ
  use geo, only : triangle
  IMPLICIT NONE
!.2-----  Declaration
  type(mesh), intent(inout) :: msh !! Maillage / Mesh
  real(kind=kind(0.d0))     :: x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4
  real(kind=kind(0.d0))     :: xn1,yn1,zn1,surf1,xm1,ym1,zm1,xn2,yn2,zn2,surf2,xm2,ym2,zm2
  real(kind=kind(0.d0))     :: xn,yn,zn,surf,xm,ym,zm,d13,d24
  integer                   :: i1,i2,i3,i4,nf,iv,iv1,iv2
!=========================== DEBUT DU CODE EXECUTABLE ==================

  DO nf=1,msh%nb_face
  ! on decoupe chaque face en 2 triangles
  ! coplanaires !!!!!
      i1=msh%list_face(nf)%vertex(1) ! n° du sommet 1
      i2=msh%list_face(nf)%vertex(2)
      i3=msh%list_face(nf)%vertex(3)
      
      x1=msh%list_vertex(i1)%x
      y1=msh%list_vertex(i1)%y
      z1=msh%list_vertex(i1)%z
      
      x2=msh%list_vertex(i2)%x
      y2=msh%list_vertex(i2)%y
      z2=msh%list_vertex(i2)%z
      
      x3=msh%list_vertex(i3)%x
      y3=msh%list_vertex(i3)%y
      z3=msh%list_vertex(i3)%z
      
      IF (msh%list_face(nf)%nbvertex.eq.4) THEN
          i4=msh%list_face(nf)%vertex(4)
          
          x4=msh%list_vertex(i4)%x
          y4=msh%list_vertex(i4)%y
          z4=msh%list_vertex(i4)%z
          
          d13=(x3-x1)**2+(y3-y1)**2+(z3-z1)**2
          d24=(x4-x2)**2+(y4-y2)**2+(z4-z2)**2
          
          IF (d13.le.d24) THEN
              CALL triangle(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn1,yn1,zn1,surf1,xm1,ym1,zm1)
              CALL triangle(x1,y1,z1,x3,y3,z3,x4,y4,z4,xn2,yn2,zn2,surf2,xm2,ym2,zm2)
          ELSE
              CALL triangle(x1,y1,z1,x2,y2,z2,x4,y4,z4,xn1,yn1,zn1,surf1,xm1,ym1,zm1)
              CALL triangle(x2,y2,z2,x3,y3,z3,x4,y4,z4,xn2,yn2,zn2,surf2,xm2,ym2,zm2)
          END IF
          
          surf=surf1+surf2
          
          IF (surf.ne.0) THEN
              xn=(surf1*xn1+surf2*xn2)/surf
              yn=(surf1*yn1+surf2*yn2)/surf
              zn=(surf1*zn1+surf2*zn2)/surf
              xm=(surf1*xm1+surf2*xm2)/surf
              ym=(surf1*ym1+surf2*ym2)/surf
              zm=(surf1*zm1+surf2*zm2)/surf
          ELSE
              xn=-1
              yn=0
              zn=0
              xm=0
              ym=0
              zm=0
          END IF

      ELSE
          CALL triangle(x1,y1,z1,x2,y2,z2,x3,y3,z3,xn,yn,zn,surf,xm,ym,zm)
      END IF

      msh%list_face(nf)%center%x=xm
      msh%list_face(nf)%center%y=ym
      msh%list_face(nf)%center%z=zm
      
      msh%list_face(nf)%surf=surf
      
      msh%list_face(nf)%vnorm%v(1)=xn
      msh%list_face(nf)%vnorm%v(2)=yn
      msh%list_face(nf)%vnorm%v(3)=zn

  END DO

! Calcul des rapports Volume/surface
   DO iv=1,msh%nb_cell
     IF(msh%list_cell(iv)%idom.eq.msh%numdom)msh%list_cell(iv)%h=0.d0
    end DO

  DO nf=1,msh%nb_face
      iv1=msh%list_face(nf)%ic1
      iv2=msh%list_face(nf)%ic2
      IF (iv1.gt.0.and.iv2.ne.-999) THEN   !!!!!!!!!!!!! ON ESSAYE !!!!!!!!!!!
          IF(msh%list_cell(iv1)%idom.eq.msh%numdom)msh%list_cell(iv1)%h=msh%list_cell(iv1)%h+msh%list_face(nf)%surf
      END IF
      IF (iv2.gt.0) THEN
          IF(msh%list_cell(iv2)%idom.eq.msh%numdom)msh%list_cell(iv2)%h=msh%list_cell(iv2)%h+msh%list_face(nf)%surf
      END IF
  END DO

  msh%small=1.d20
  msh%big  =0.d0
  DO iv=1,msh%nb_cell
    IF(msh%list_cell(iv)%idom.eq.msh%numdom)then
        msh%list_cell(iv)%h=msh%list_cell(iv)%vol/msh%list_cell(iv)%h
        if(msh%list_cell(iv)%h.gt.msh%big)  msh%big  =msh%list_cell(iv)%h
        if(msh%list_cell(iv)%h.lt.msh%small)msh%small=msh%list_cell(iv)%h
     endif
  END DO
  write( *,1000)msh%numdom,msh%small
  write( *,1010)msh%numdom,msh%small
1000    format('INIT_MESH-INFO Domain ',i3,' size of smallest cell  :',e13.5)
1010    format('INIT_MESH-INFO Domain ',i3,' size of largest cell   :',e13.5)
  msh%list_face(:)%flevel=0
  msh%list_cell(:)%clevel=0
  msh%nb_level_cell=int(dlog(msh%big*1.05/msh%small)/dlog(2.0d0))
  msh%nb_level_face=msh%nb_level_cell

!===========================   FIN DE LA ROUTINE    ====================
END SUBROUTINE INIT_MESH
!DOCFIN 

