mesh_interfaces.F90

Go to the documentation of this file.

module mesh_interfaces

  
  use choral_constants
  use choral_variables
  use real_type
  use algebra_set
  use basic_tools
  use graph_mod
  use cell_mod
  use mesh_mod

  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: define_interfaces
  public :: define_domain_boundary       !! TESTED

contains


  subroutine define_interfaces(m)
    type(mesh), intent(inout) :: m
    
    if (choral_verb>0) write(*,*) &
         & "mesh_interfaces : define_interfaces"

    select case(m%dim)
    case(1)
       call  define_1dmesh_vertexes(m) 

    case(2)
       call define_2dmesh_edges(m) 

    case(3)
       call define_3dmesh_faces(m)
          
    end select

    if (choral_verb>1) write(*,*) &
         &"  Number of interfaces           =", m%nbItf

  end subroutine define_interfaces



  subroutine define_domain_boundary(m)
    type(mesh), intent(inout)  :: m

    integer, dimension(:,:), allocatable :: tab, clTab
    integer, dimension(:)  , allocatable :: clType
    integer, dimension(CELL_MAX_NBNODES) :: cl_nodes

    integer :: cpt1, cpt2, cpt3
    integer :: ii, cl, idx, j1, j2, nNd

    if (choral_verb>1) write(*,*)&
         &"mesh_interfaces : define_domain_boundary"
    
    if (.NOT.valid(m)) call quit(&
         & "mesh_interfaces: define_domain_boundary: mesh 'm' not valid")

    if (m%nbItf<=0) call define_interfaces(m)  

    call boundary_interface_scan(tab, m)

    !! Flag the boundary cells already present in the mesh
    !! and checks wether boundary cells are missing
    !!
    cpt1 = 0
    cpt2 = 0
    do ii=1, m%nbItf

       !! internal interface
       if (tab(1,ii)==0) cycle

       !! boundary interface
       cl = tab(3,ii)
       if (cl==0) then  
          !! interface ii is a boundary cell 
          !! that does not presently belong to the mesh
          cpt2 = cpt2 + 1

       else
          !! interface ii is a boundary cell 
          !! that belongs to the mesh
          cpt1 = cpt1 +1

       end if

    end do

    !! All the boundary cells already belong to the mesh
    !!
    if (cpt2==0) then   
       call freemem(tab)

    !! Boundary cells need to be added
    !! The mesh is rebuilt including these cells.
    !!
    !! For the new cells :
    !!     cell types will be in clType
    !!     cell nodes will be in clTab
    !!
    else

       call allocmem(cltab , cell_max_nbnodes, cpt2 ) 
       call allocmem(cltype, cpt2 ) 
       !! 
       cpt3 = 0
       do ii=1, m%nbItf

          !! internal interface
          if (tab(1,ii)==0) cycle  

          !! interface ii is a boundary cell 
          !! that belongs to the mesh
          if (tab(3,ii)/=0) cycle         

          !! interface ii as a new cell
          !! whose index will be 'cpt3'
          cpt3 = cpt3 + 1

          !! the co-boundary of the interface 'ii' 
          !! is the cell 'cl' with 'hNd' nodes 
          !! and with nodes 'cl_nodes(1:nNd)'
          cl = tab(1,ii)
          j1  = m%clToNd%row(cl)
          j2  = m%clToNd%row(cl+1)
          nnd = j2-j1
          cl_nodes(1:nnd) = m%clToNd%col(j1: j2-1)
          !!
          !! The interface 'ii' is the interface 'idx'
          !! of the cell 'cl' for the local indexing 
          idx = tab(2,ii)
          !!
          select case(m%clType(cl))
          case(cell_edg, cell_edg_2)
             cltype(cpt3) = cell_vtx   
             
             cltab(1,cpt3) = cl_nodes(idx)
             
          case(cell_trg)
             cltype(cpt3) = cell_edg
             
             cltab(1:2,cpt3) = cl_nodes( cell_ed_vtx(1:2, idx, cell_trg) ) 
             
          case(cell_tet)
             cltype(cpt3) = cell_trg
             
             cltab(1:3,cpt3) = cl_nodes( cell_fc_vtx(1:3, idx, cell_tet) )
             
          case default
             call quit("mesh_interfaces: define_domain_boundary:&
                  & cell type not implemented: "&
                  & // trim(cell_name(tab(2,ii))))

          end select

       end do
       call freemem(tab)

       !! Add the new cells
       !!
       call add_cells(m, cltype, cltab)

    end if
       
    if (choral_verb>1) write(*,*)&
         &"  Mesh domain boundary           =", cpt1+cpt2, " cells"

  end subroutine define_domain_boundary

  subroutine boundary_interface_scan(tab, m)
    type(mesh)             , intent(inout)  :: m
    integer, dimension(:,:), allocatable    :: tab

    integer :: itf, itf2, itf_idx, itf_nb_vtx
    integer :: j1, j2, jj, ll
    integer :: dim, cl, cl_type, wdt

    integer, dimension(CELL_MAX_FC_NBVTX) :: itf_vtx
    integer, dimension(:), allocatable    :: t0, t1, t2

    if (choral_verb>2) write(*,*)&
         &"mesh_interfaces : boundary_interface_scan "

    if (m%nbItf<=0) call define_interfaces(m)  

    if (.NOT.valid(m)) call quit(&
         & "mesh_interfaces: boundary_interface_scan: mesh 'm' not valid")

    dim = m%dim
    call allocmem(tab, 3, m%nbItf)
    tab = 0

    !! max number of cell-neighbour for a node
    !!
    wdt = m%ndToCl%width
    allocate(t0(wdt), t1(wdt), t2(wdt))

    
    !! First determines the boundary interfaces
    !! and tab(1,itf) = co-boundary of these interfaces
    !!
    do itf=1, m%itfToCl%nl
       j1 =  m%itfToCl%row(itf)
       j2 =  m%itfToCl%row(itf+1)
       j2 = j2-j1

       !! boundary interface
       if (j2==1) then
          cl          = m%itfToCl%col(j1)
          tab(1, itf) = cl
       end if

    end do


    !! Second determines if there is a second cell
    !! with vertexes the vertexes of a boundary interface 
    !!
    do itf=1, m%itfToCl%nl
       if (tab(1, itf) == 0) cycle

       !! interface itf co-boundary is cell cl
       !!
       cl      = tab(1, itf)
       cl_type = m%clType(cl)

       !! interface itf has index itf_idx
       !! in the local numbering associated to cell cl
       !!
       j1 = m%clToItf%row(cl)
       j2 = m%clToItf%row(cl+1) - 1
       loop_1: do jj= j1, j2
          itf2 = m%clToItf%col(jj)

          if (itf == itf2) then
             itf_idx = jj-j1+1
             exit loop_1
          end if
       end do loop_1
       tab(2, itf) = itf_idx

       !! interface itf has vertexes itf_vtx(1:itf_nb_vtx)
       !!
       select case(dim)
       case(1)
          j1         = m%clToNd%row(cl) - 1
          j2         = j1 + itf_idx
          itf_vtx(1) = m%clToNd%col(j2)

       case(2)
          itf_vtx(1:2) = cell_ed_vtx(:,itf_idx, cl_type)

          j1           = m%clToNd%row(cl) - 1

          j2           = j1 + itf_vtx(1)
          itf_vtx(1)   = m%clToNd%col(j2)

          j2           = j1 + itf_vtx(2)
          itf_vtx(2)   = m%clToNd%col(j2)

       case(3)
          itf_nb_vtx = cell_fc_nbvtx(itf_idx, cl_type)

          itf_vtx(1:itf_nb_vtx) = &
               &     cell_fc_vtx(1:itf_nb_vtx,itf_idx, cl_type)

          j1 = m%clToNd%row(cl) - 1

          do jj=1, itf_nb_vtx
             j2           = j1 + itf_vtx(jj)
             itf_vtx(jj)  = m%clToNd%col(j2)
          end do

       end select

       !! determine what cell in the mesh has 
       !! the nodes itf_vtx(1:itf_nb_vtx) among its nodes
       !!
       select case(dim)
       case(1)
          call getrow(jj, t0, wdt, m%ndToCl,  itf_vtx(1))

       case(2)
          call getrow(j1, t1, wdt, m%ndToCl,  itf_vtx(1))
          call getrow(j2, t2, wdt, m%ndToCl,  itf_vtx(2))

          call cap_sorted_set(jj, t0, wdt, t1(1:j1), j1, t2(1:j2), j2)

       case(3)
          call getrow(j1, t1, wdt, m%ndToCl,  itf_vtx(1))

          loop_2: do ll=2, itf_nb_vtx
             call getrow(j2, t2, wdt, m%ndToCl,  itf_vtx(ll))

             call cap_sorted_set(jj, t0, wdt, t1(1:j1), j1, t2(1:j2), j2)
             

             if (ll==itf_nb_vtx) exit loop_2
             j1 = jj
             t1(1:j1) = t0(1:j1) 

          end do loop_2

       end select

#if DBG>0
       if ((jj==0) .OR. (jj>2)) call quit(&
            & "mesh_interfaces: boundary_interface_scan: 3")
#endif

       if (jj==2) then

          if (t0(1)==cl) then
             tab(3, itf) = t0(2)
          else
             tab(3, itf) = t0(1)
          end if
          
#if DBG>0
          j1 = m%clType(tab(3, itf) )
          j2 = cell_dim( j1 )
          if (j2 /= dim-1) call quit(&
               & "mesh_interfaces: boundary_interface_scan: 4")
#endif
       end if

    end do

  end subroutine boundary_interface_scan

  subroutine define_2dmesh_edges(m)  
    type(mesh), intent(inout) :: m

    integer, dimension(:,:), allocatable :: edTab
    integer, dimension(:)  , allocatable :: nnz
    integer :: cl, j1, j2

    integer  :: cl_t, nb_vtx, nb_ed, width, nb_itf
    real(DP) :: t1, t2

    if (choral_verb>1) write(*,*) &
         & "mesh_tools      : define_2DMesh_edges"

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    cells --> interfaces
    !!
    t1 = clock()

    !! STEP 1: compute the graph pattern 
    !!
    call allocmem(nnz, m%nbCl)
    do cl=1, m%nbCl
       cl_t = m%clType(cl)
       if ( cell_dim(cl_t) ==2 ) then
          nnz(cl) = cell_nbed(cl_t)
       else
          nnz(cl) = 0
       end if
    end do
    m%clToItf = graph(nnz, m%nbCl, 0)
    call freemem(nnz)

    !! STEP 2 : create edTab 
    !!
    call allocmem(edtab, 3, m%clToItf%nnz)
    width  = m%ndToCl%width
    do cl_t=1, cell_tot_nb
       if ( cell_dim(cl_t)     /= 2 ) cycle
       if ( m%cell_count(cl_t) == 0 ) cycle

       nb_vtx = cell_nbvtx(cl_t)
       nb_ed  = cell_nbed(cl_t)
       call cell_loop()

    end do

    !! STEP 3 : number interfaces
    !!
    call number_edges(nb_itf)
    m%nbItf = nb_itf
    m%clToItf%nc = nb_itf

    !! fill in the conectivity graph 2D-cells --> edges
    !!
    do cl=1, m%nbCl

       j1 = m%clToItf%row(cl)
       j2 = m%clToItf%row(cl+1)-1

       m%clToItf%col(j1:j2) = edtab(3,j1:j2)
       
    end do

    call freemem(edtab)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  2D-cells --> edges graph, CPU  = ", &
         & real(t2, sp)

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    interfaces --> cells
    !!
    t1 = clock() 
    call transpose(m%itfToCl, m%clToItf)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  edges --> 2D-cells graph, CPU  = ", &
         & real(t2, sp)

  contains

    !! Array of all the 2D-cell edges,
    !! internal edges are counted twice
    !!
    !! Step 1 = get the co-boundary of all 2D-cell edgees
    !!
    subroutine cell_loop()

      integer, dimension(nb_vtx)         :: vtx_nbCl
      integer, dimension(width, nb_vtx)  :: vtx_cl
      integer, dimension(nb_vtx)         :: cl_vtx
      integer, dimension(3)              :: coBord
      integer :: ed, v1, v2
      integer :: l1, l2

      !$OMP PARALLEL PRIVATE(ed, v1, v2, j1, j2, l1, l2, vtx_nbCl, vtx_cl, cl_vtx, coBord) 
      !$OMP DO   
      do cl=1, m%nbCl
         l1 = m%clType(cl) 
         if (l1 /= cl_t ) cycle

         !! cell cl vertexes
         !!
         j1 = m%clToNd%row(cl)
         j2 = j1 + nb_vtx - 1
         cl_vtx =  m%clToNd%col( j1:j2 )

         !! neighbour cells for each vertex
         !! of the cell cl
         !!
         do v1=1, nb_vtx
            l1 = cl_vtx(v1)
            j1 = m%ndToCl%row(l1)
            j2 = m%ndToCl%row(l1+1)

            l1 = j2-j1
            vtx_nbcl(v1) = l1

            !! vtx_cl is sorted !
            vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)

         end do

         !! build the array of edges 
         !!
         j1 = m%clToItf%row(cl) 
         do ed=1, nb_ed

            edtab(1, j1) = cl

            !! list the cells sharing the vertexes v1 and v2
            !! coBord is sorted !
            v1 = cell_ed_vtx(1, ed, cl_t)
            v2 = cell_ed_vtx(2, ed, cl_t)
            l1 = vtx_nbcl(v1)
            l2 = vtx_nbcl(v2)
            call cap_sorted_set(j2, cobord, 3, &
                 & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)

            !! j2 = number of cells sharing this edge
            select case(j2)

            case(1)   !! boundary edge
               edtab(2, j1) = cl 

            case(2)
               !! In this case the edge can be 
               !! either internal or not,

               if (cobord(1)==cl) then 
                  !! it is not necessary to know whether
                  !! the edge is internal or not here
                  !!
                  edtab(2, j1) = cobord(2)

               else        
                  !! it is necessary to know whether
                  !! the edge is internal or not here
                  !!
                  !! For this we check the dimension 
                  !! of the second cell
                  l1 = cell_dim( m%clType(cobord(1)) )

                  if ( l1==2 ) then           !! internal edge
                     edtab(2, j1) = cobord(1)

                  else                        !! boundary edge
                     edtab(2, j1) = cl

                  end if
               end if

            case(3) 
               !! In this case the edge is internal,
               !! but one of the three cells in the
               !! co-boundary is of dimension 1.
               !!
               l1 = cell_dim( m%clType(cobord(1)) )
               if (l1==1) then

                  if (cobord(2)==cl) then
                     edtab(2, j1) = cobord(3)
                  else
                     edtab(2, j1) = cobord(2)
                  end if

               else

                  l1 = cell_dim( m%clType(cobord(2)) )
                  if (l1==1) then

                     if (cobord(1)==cl) then
                        edtab(2, j1) = cobord(3)
                     else
                        edtab(2, j1) = cobord(1)
                     end if

                  else

                     if (cobord(1)==cl) then
                        edtab(2, j1) = cobord(2)
                     else
                        edtab(2, j1) = cobord(1)
                     end if

                  end if
               end if

            end select

            !! increment the edge
            j1 = j1 + 1
         end do
      end do
      !$OMP END DO
      !$OMP END PARALLEL

    end subroutine cell_loop


    !! Array of all the 2D-cell edges,
    !! internal edges are counted twice
    !!
    !! Step 2 = global numberingg of all the 2D-cell edgees
    !!
    subroutine number_edges(cpt)
      integer, intent(out) :: cpt

      integer :: jj, cl, cl2, ll

      cpt = 0
      do jj=1, m%clToItf%nnz
         
         cl  = edtab(1,jj)
         cl2 = edtab(2,jj)
         
         if ( cl <= cl2 ) then
            
            cpt = cpt + 1
            edtab(3,jj) = cpt
            
         else
            ll = m%clToItf%row(cl2)
            
            do while( edtab(2,ll) /= cl ) 
               ll = ll + 1
            end do
            edtab(3,jj) = edtab(3,ll)
            
         end if
       
      end do

    end subroutine number_edges

  end subroutine define_2dmesh_edges


  subroutine define_3dmesh_faces(m)  
    type(mesh), intent(inout) :: m

    integer, dimension(:,:), allocatable :: fcTab
    integer, dimension(:)  , allocatable :: nnz
    integer :: cl, j1, j2

    integer  :: cl_t, nb_vtx, nb_fc, width, nb_itf
    real(DP) :: t1, t2

    if (choral_verb>1) write(*,*) &
         & "mesh_tools      : define_3DMesh_faces"

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    cells --> interfaces
    !!
    t1 = clock()

    !! STEP 1: compute the graph pattern 
    !!
    call allocmem(nnz, m%nbCl)
    do cl=1, m%nbCl
       cl_t = m%clType(cl)
       if ( cell_dim(cl_t) ==3 ) then
          nnz(cl) = cell_nbfc(cl_t)
       else
          nnz(cl) = 0
       end if
    end do
    m%clToItf = graph(nnz, m%nbCl, 0)
    call freemem(nnz)

    !! STEP 2 : create fcTab 
    !!
    call allocmem(fctab, 3, m%clToItf%nnz)
    width  = m%ndToCl%width
    do cl_t=1, cell_tot_nb
       if ( cell_dim(cl_t)     /= 3 ) cycle
       if ( m%cell_count(cl_t) == 0 ) cycle

       nb_vtx = cell_nbvtx(cl_t)
       nb_fc  = cell_nbfc(cl_t)
       call cell_loop()

    end do

    !! STEP 3 : number faces
    !!
    call number_faces(nb_itf)
    m%nbItf = nb_itf
    m%clToItf%nc = nb_itf

    !! fill in the conectivity graph 2D-cells --> edges
    !!
    do cl=1, m%nbCl

       j1 = m%clToItf%row(cl)
       j2 = m%clToItf%row(cl+1)-1

       m%clToItf%col(j1:j2) = fctab(3,j1:j2)
       
    end do

    call freemem(fctab)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  cells --> faces graph,    CPU  = ", &
         & real(t2, sp)

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    interfaces --> cells
    !!
    t1 = clock() 
    call transpose(m%itfToCl, m%clToItf)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  faces --> cells graph,    CPU  = ", &
         & real(t2, sp)

  contains

    !! Array of all the 2D-cell edges,
    !! internal edges are counted twice
    !!
    !! Step 1 = get the co-boundary of all 2D-cell edgees
    !!
    subroutine cell_loop()

      integer, dimension(nb_vtx)         :: vtx_nbCl
      integer, dimension(width, nb_vtx)  :: vtx_cl
      integer, dimension(nb_vtx)         :: cl_vtx
      integer, dimension(width)          :: coBord, aux
      integer :: fc, v1, v2
      integer :: l1, l2, ii

      !$OMP PARALLEL PRIVATE(fc, v1, v2, j1, j2, l1, l2, vtx_nbCl, vtx_cl, cl_vtx, coBord, aux, ii) 
      !$OMP DO   
      do cl=1, m%nbCl
         l1 = m%clType(cl) 
         if (l1 /= cl_t ) cycle

         !! cell cl vertexes
         !!
         j1 = m%clToNd%row(cl)
         j2 = j1 + nb_vtx - 1
         cl_vtx =  m%clToNd%col( j1:j2 )

         !! neighbour cells for each vertex
         !! of the cell cl
         !!
         do v1=1, nb_vtx
            l1 = cl_vtx(v1)
            j1 = m%ndToCl%row(l1)
            j2 = m%ndToCl%row(l1+1)

            l1 = j2-j1
            vtx_nbcl(v1) = l1

            !! vtx_cl is sorted !
            vtx_cl(1:l1, v1) = m%ndToCl%col(j1:j2-1)

         end do

         !! build the array of faces 
         !!
         j1 = m%clToItf%row(cl) 
         do fc=1, nb_fc

            fctab(1, j1) = cl

            !! list the cells sharing the vertexes of the 
            !! face fc
            v1 = cell_fc_vtx(1, fc, cl_t)
            v2 = cell_fc_vtx(2, fc, cl_t)
            l1 = vtx_nbcl(v1)
            l2 = vtx_nbcl(v2)
            call cap_sorted_set(j2, cobord, width, &
                 & vtx_cl(1:l1, v1), l1, vtx_cl(1:l2, v2), l2)

            do ii=3, cell_fc_nbvtx(fc, cl_t)
               v2 = cell_fc_vtx(ii, fc, cl_t)
               l2 = vtx_nbcl(v2)

               l1 = j2
               aux(1:l1) = cobord(1:l1)

               call cap_sorted_set(j2, cobord(1:l1), l1, &
                 & aux(1:l1), l1, vtx_cl(1:l2, v2), l2)

            end do

            !! j2 = number of cells sharing this face
            select case(j2)
            case(1)   !! boundary face
               fctab(2, j1) = cl 

            case(2)
               !! In this case the face can be 
               !! either internal or not,

               if (cobord(1)==cl) then 
                  !! it is not necessary to know whether
                  !! the face is internal or not here
                  !!
                  fctab(2, j1) = cobord(2)

               else        
                  !! it is necessary to know whether
                  !! the face is internal or not here
                  !!
                  !! For this we check the dimension 
                  !! of the second cell
                  l1 = cell_dim( m%clType(cobord(1)) )

                  if ( l1==3 ) then           !! internal edge
                     fctab(2, j1) = cobord(1)

                  else                        !! boundary edge
                     fctab(2, j1) = cl

                  end if
               end if

            case(3) 

               !! In this case the edge is internal,
               !! but one of the three cells in the
               !! co-boundary is of dimension 2.
               !!
               l1 = cell_dim( m%clType(cobord(1)) )
               if (l1==2) then

                  if (cobord(2)==cl) then
                     fctab(2, j1) = cobord(3)
                  else
                     fctab(2, j1) = cobord(2)
                  end if

               else

                  l1 = cell_dim( m%clType(cobord(2)) )
                  if (l1==2) then

                     if (cobord(1)==cl) then
                        fctab(2, j1) = cobord(3)
                     else
                        fctab(2, j1) = cobord(1)
                     end if

                  else

                     if (cobord(1)==cl) then
                        fctab(2, j1) = cobord(2)
                     else
                        fctab(2, j1) = cobord(1)
                     end if

                  end if
               end if

            end select

            !! increment the edge
            j1 = j1 + 1
         end do
      end do
      !$OMP END DO
      !$OMP END PARALLEL

    end subroutine cell_loop


    !! Array of all the 3D-cell faces,
    !! internal faces are counted twice
    !!
    !! Step 2 = global numberingg of all the 3D-cell faces
    !!
    subroutine number_faces(cpt)
      integer, intent(out) :: cpt

      integer :: jj, cl, cl2, ll

      cpt = 0
      do jj=1, m%clToItf%nnz
         
         cl  = fctab(1,jj)
         cl2 = fctab(2,jj)
         
         if ( cl <= cl2 ) then
            
            cpt = cpt + 1
            fctab(3,jj) = cpt
            
         else
            ll = m%clToItf%row(cl2)
            
            do while( fctab(2,ll) /= cl ) 
               ll = ll + 1
            end do
            fctab(3,jj) = fctab(3,ll)
            
         end if
       
      end do

    end subroutine number_faces

  end subroutine define_3dmesh_faces



  subroutine define_1dmesh_vertexes(m)  
    type(mesh), intent(inout) :: m

    integer, dimension(:,:), allocatable :: vtxTab
    integer, dimension(:)  , allocatable :: nnz
    integer :: cl, j1, j2, cl_t

    integer  :: width, nb_itf
    real(DP) :: t1, t2

    if (choral_verb>1) write(*,*) &
         & "mesh_tools      : define_1DMesh_vertexes"

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    cells --> interfaces
    !!
    t1 = clock()

    !! STEP 1: compute the graph pattern 
    !!
    call allocmem(nnz, m%nbCl)
    do cl=1, m%nbCl
       cl_t = m%clType(cl)
       if ( cell_dim(cl_t) == 1 ) then
          nnz(cl) = 2
       else
          nnz(cl) = 0
       end if
    end do
    m%clToItf= graph(nnz, m%nbCl, 0)
    call freemem(nnz)

    !! STEP 2 : create vtxTab 
    !!
    call allocmem(vtxtab, 3, m%clToItf%nnz)
    width  = m%ndToCl%width
    call cell_loop()


    !! STEP 3 : number interfaces
    !!
    call number_vertexes(nb_itf)
    m%nbItf = nb_itf
    m%clToItf%nc = nb_itf

    !! fill in the conectivity graph 1D-cells --> vertexes
    !!
    do cl=1, m%nbCl

       j1 = m%clToItf%row(cl)
       j2 = m%clToItf%row(cl+1)-1

       m%clToItf%col(j1:j2) = vtxtab(3,j1:j2)
       
    end do

    call freemem(vtxtab)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  1D-cells --> vertex graph, CPU = ", &
         & real(t2, sp)

    !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !!
    !! Computation of the connectivity graph
    !!    interfaces --> cells
    !!
    t1 = clock() 
    call transpose(m%itfToCl, m%clToItf)
    t2 = clock() - t1
    if (choral_verb>1) write(*,*) &
         & "  vertex --> 1D-cells graph, CPU = ", &
         & real(t2, sp)

  contains

    !! Array of all the 1D-cell vertexes,
    !! internal vertexes are counted twice
    !!
    !! Step 1 = get the co-boundary of all 1D-cell vertexes
    !!
    subroutine cell_loop()

      integer, dimension(2)      :: cl_vtx
      integer, dimension(width)  :: coBord
      integer :: vtx, v, l1, l2

      do cl=1, m%nbCl
         l1 = m%clType(cl) 
         if ( cell_dim(l1) /= 1 ) cycle

         !! cell cl vertexes
         !!
         j1 = m%clToNd%row(cl)
         j2 = j1 + 1
         cl_vtx =  m%clToNd%col( j1:j2 )

         !! build the array of vertexes 
         !!
         j1 = m%clToItf%row(cl) 
         do vtx=1, 2

            vtxtab(1, j1) = cl

            v  = cl_vtx(vtx)
            l1 = m%ndToCl%row(v)
            l2 = m%ndToCl%row(v+1)

            j2 = l2-l1           
            cobord(1:j2) = m%ndToCl%col(l1:l2-1)

            !! j2 = number of cells sharing this vertex
            select case(j2)

            case(1)   !! boundary vertex
               vtxtab(2, j1) = cl 

            case(2)
               !! In this case the vertex can be 
               !! either internal or not,

               if (cobord(1)==cl) then 
                  !! it is not necessary to know whether
                  !! the vertex is internal or not here
                  !!
                  vtxtab(2, j1) = cobord(2)

               else        
                  !! it is necessary to know whether
                  !! the vertex is internal or not here
                  !!
                  !! For this we check the dimension 
                  !! of the second cell
                  l1 = cell_dim( m%clType(cobord(1)) )

                  if ( l1==1 ) then           !! internal vertex
                     vtxtab(2, j1) = cobord(1)

                  else                        !! boundary vertex
                     vtxtab(2, j1) = cl

                  end if
               end if

            case(3) 
               !! In this case the vertex is internal,
               !! but one of the three cells in the
               !! co-boundary is of dimension 1.
               !!
               l1 = cell_dim( m%clType(cobord(1)) )
               if (l1==0) then

                  if (cobord(2)==cl) then
                     vtxtab(2, j1) = cobord(3)
                  else
                     vtxtab(2, j1) = cobord(2)
                  end if

               else

                  l1 = cell_dim( m%clType(cobord(2)) )
                  if (l1==0) then

                     if (cobord(1)==cl) then
                        vtxtab(2, j1) = cobord(3)
                     else
                        vtxtab(2, j1) = cobord(1)
                     end if

                  else

                     if (cobord(1)==cl) then
                        vtxtab(2, j1) = cobord(2)
                     else
                        vtxtab(2, j1) = cobord(1)
                     end if

                  end if
               end if

            case default
               call quit("mesh_tools: define_1DMesh_vertexes: non conformal 1D mesh")

            end select

            !! increment the vertex
            j1 = j1 + 1
         end do
      end do

    end subroutine cell_loop


    !! Array of all the 2D-cell vertexes,
    !! internal vertexes are counted twice
    !!
    !! Step 2 = global numberingg of all the 2D-cell vertexes
    !!
    subroutine number_vertexes(cpt)
      integer, intent(out) :: cpt

      integer :: jj, cl, cl2, ll

      cpt = 0
      do jj=1, m%clToItf%nnz
         
         cl  = vtxtab(1,jj)
         cl2 = vtxtab(2,jj)
         
         if ( cl <= cl2 ) then
            
            cpt = cpt + 1
            vtxtab(3,jj) = cpt
            
         else
            ll = m%clToItf%row(cl2)
            
            do while( vtxtab(2,ll) /= cl ) 
               ll = ll + 1
            end do
            vtxtab(3,jj) = vtxtab(3,ll)
            
         end if
       
      end do

    end subroutine number_vertexes

  end subroutine define_1dmesh_vertexes

  subroutine add_cells(m, type, nodes)
    type(mesh), intent(inout)  :: m
    integer, dimension(:)  , intent(in) :: type
    integer, dimension(:,:), intent(in) :: nodes

    integer, dimension(:)  , allocatable :: clType
    integer, dimension(:,:), allocatable :: clTab

    integer :: nn, j1, j2, nNd, cl

    nn = size(type,1)
    if (choral_verb>1) write(*,*) &
         &"mesh_tools      : add_cells      =", nn, " new cells"

    !! For the new mesh:
    !!
    !! cell types will be in clType
    !! cell nodes will be in clTab
    !!
    call allocmem(cltab , cell_max_nbnodes, nn + m%nbCl) 
    call allocmem(cltype, nn + m%nbCl) 

    !! initialise clType and clTab
    !! with the cells already present in the mesh
    !!
    cltype(nn+1:) = m%clType
    do cl = 1, m%nbCl

       j1  = m%clToNd%row(cl)
       j2  = m%clToNd%row(cl+1)
       nnd = j2-j1

       cltab(1:nnd, cl + nn) = m%clToNd%col(j1: j2-1)

    end do


    !! Add the new boundary cells to 
    !! clType(1:cpt2) and to clTab(:, 1:cpt2)
    !! 
    cltype(1:nn)  = type
    cltab(:,1:nn) = nodes

    !! Update m%clType
    !!
    call allocmem(m%clType, nn + m%nbCl)
    m%clType = cltype
    call freemem(cltype)

    !! build the new mesh including the boundary cells
    !!
    call mesh_clear_2(m)
    call mesh_create_cltab(m, cltab)
    call freemem(cltab)
    call mesh_create_end(m)

    if (.NOT.valid(m)) call quit('mesh_interfaces: &
         & add_cells: construction not valid')
       
  end subroutine add_cells

end module mesh_interfaces