feSpace_mod.F90

Go to the documentation of this file.

!
module fespace_mod

  use choral_constants
  use choral_variables
  use real_type
  use abstract_interfaces, only: r3tor3, r3tor
  use algebra_set
  use basic_tools
  use graph_mod
  use cell_mod
  use quad_mod
  use mesh_mod
  use fe_mod
  use geotsf_mod
  
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  public :: fespace

  public :: set, assemble             !! TESTED
  public :: clear
  public :: interp_scal_func          !! TESTED
  public :: interp_vect_func          !! TESTED
  public :: interp_scal_fe            !! TESTED
  public :: eval_scal_fe_at_points    !! TESTED
  public :: print, valid, write, gmsh_addview
  public :: flag_x_h_dof
  
  
  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type fespace

     !! DOF
     integer :: nbdof=0

     !! feType      = finite element method type per mesh cell
     !! dofCoord    = dof nodes coordinates
     !! dofType(ii) = type of dof ii
     real(RP), dimension(:,:), allocatable :: dofcoord
     integer , dimension(:)  , allocatable :: doftype 
     integer , dimension(:)  , allocatable :: fetype  

     !! associated mesh 
     type(mesh), pointer :: mesh=>null()

     !! cell dof local indexing --> dof global indexing
     !! cell to dof connectivity
     type(graph) :: cltodof

     !! dof_count(CELL_XXX) = number of dof
     !!                       having CELL_XXX as geometry
     !! fe_count(FE_XXX   ) = number of cells
     !!                       with the FE method FE_XXX
     integer, dimension(CELL_TOT_NB) :: dof_count = 0
     integer, dimension(0:FE_TOT_NB) :: fe_count  = 0

   contains

     final :: fespace_clear

  end type fespace


  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure fespace_clear
  end interface clear

  interface valid
     module procedure fespace_valid
  end interface valid

  interface fespace
     module procedure fespace_create
  end interface fespace

  interface set
     module procedure fespace_set
  end interface set

  interface print
     module procedure fespace_print
  end interface print

  interface write
     module procedure fespace_write
  end interface write

  interface assemble
     module procedure fespace_assemble
  end interface assemble

  interface interp_vect_func
     module procedure fespace_interp_vect_func
  end interface interp_vect_func

contains

  subroutine fespace_clear(X_h)
    type(fespace), intent(inout) :: X_h

    x_h%nbDof = 0
    x_h%dof_count = 0
    x_h%fe_count  = 0

    x_h%mesh => null()

    call freemem(x_h%feType)
    call freemem(x_h%dofType)
    call freemem(x_h%dofCoord)
    call clear(x_h%clToDof)

  end subroutine fespace_clear

  function fespace_create(m) result(X_h)
    type(fespace)         :: X_h
    type(mesh)   , target :: m

    if (choral_verb>0) write(*,*) 'feSpace_mod     : create' 
    call clear(x_h)
    if(.NOT.valid(m)) call quit( &
         & "feSpace_mod: feSpace_create:  mesh 'm' not valid" )
    x_h%mesh => m
    call allocmem( x_h%feType, m%nbCl )
    x_h%feType = fe_none
    if (choral_verb>1) write(*,*) &
         & '  Number of cells                =', x_h%mesh%nbCl

  end function fespace_create

  subroutine fespace_set(X_h, ft)
    type(fespace), intent(inout) :: X_h
    integer      , intent(in)    :: ft

    integer :: ii, jj, geo, cl_geo, cpt

    geo = fe_geo(ft)
    cpt = 0
    
    do ii=1, x_h%mesh%nbCl

       jj     = x_h%mesh%clType(ii)
       cl_geo = cell_geo(jj)

       if (geo==cl_geo) then
          x_h%feType(ii) = ft
          cpt = cpt + 1
       end if

    end do

    if (choral_verb>0) write(*,*) &
         &'feSpace_mod     : set            =', &
         &     cpt, ' cells ', fe_name(ft)
    
  end subroutine fespace_set



  subroutine fespace_print(X_h)
    type(fespace), intent(in) :: X_h

    integer :: ii
    
    write(*,*)"feSpace_mod     : print"
    write(*,*)"  Number of DOF                  =", x_h%nbDof
    
    write(*,*) '  DOF per geometrical type'
    do ii=1, cell_tot_nb
       if ( x_h%dof_count(ii)>0) then
          write(*,*) '                        ', &
               & cell_name(ii), ' DOF =', &
               &     x_h%dof_count(ii)
       end if
    end do
    write(*,*) '  Number of cells per fe method'
    do ii=0, fe_tot_nb
       if ( x_h%fe_count(ii)>0) then
          write(*,*) '               ', &
               &fe_name(ii), '  =', &
               &     x_h%fe_count(ii)
       end if
    end do

  end subroutine fespace_print


  function fespace_valid(X_h) result(b)
    type(fespace), intent(in) :: x_h
    logical :: b

    b = .false.
    if(.NOT. ( associated(x_h%mesh) )) return
    if(.NOT. ( valid(x_h%mesh)      )) return
    if( valid(x_h%clToDof) < 0       ) return
    if(.NOT. ( allocated(x_h%dofCoord) )) return
    if(.NOT. ( allocated(x_h%dofType)   )) return
    if(.NOT. ( allocated(x_h%feType)    )) return

    b =  x_h%nbdof>0 
    b =  b .AND. ( all( shape(x_h%dofCoord)  == (/3,x_h%nbdof/)) )
    b =  b .AND. ( size(x_h%dofType,1) == x_h%nbdof )
    b =  b .AND. ( size(x_h%feType,1)  == x_h%mesh%nbCl )
    b =  b .AND. ( x_h%clToDof%nl      == x_h%mesh%nbCl )

  end function fespace_valid


  subroutine fespace_write (X_h, fileName, fileFormat)
    type(fespace)   , intent(in) :: X_h
    character(len=*), intent(in) :: fileName, fileFormat

    if (choral_verb>1) write(*,*)&
         & "feSpace_mod     : write format   = ",&
         & trim(fileformat), " to ", trim(filename)

    select case(trim(fileformat))
    case("gmsh")
       call gmsh_write(x_h, filename)

    case default
       call quit( "feSpace_mod: feSpace_write: uncorrect file format" )

    end select

  end subroutine fespace_write

  subroutine gmsh_write (X_h, fileName)
    type(fespace)   , intent(in) :: X_h
    character(len=*), intent(in) :: fileName

    integer :: ii, cpt, ft, wdt, s1

    integer, dimension(:), allocatable :: p1    

    open(unit = 999,file = filename)

    write(999,fmt="(A)") '$MeshFormat'
    write(999,fmt="(A)") '2.2 0 8'
    write(999,fmt="(A)") '$EndMeshFormat'
    write(999,fmt="(A)") '$Nodes'
    write(999,*) x_h%nbdof
    do ii=1, x_h%nbdof
       write(999,*) ii, x_h%dofCoord(:, ii)
    end do
    write(999,fmt="(A)") '$EndNodes'
    write(999,fmt="(A)") '$Elements'

    !! count cells with non void fe method
    !!
    cpt = size(x_h%feType,1)
    do ii=1, size(x_h%feType,1)
       if (x_h%fetype(ii)==fe_none) cpt = cpt-1
    end do
    write(999,*) cpt

    !! write cells with non void fe method
    !!
    cpt = 0
    wdt = x_h%clToDof%width
    call allocmem(p1, wdt)
    do ii=1, size(x_h%feType,1)

       ft = x_h%fetype(ii)
       
       if (ft==fe_none) cycle

       cpt = cpt + 1
       call getrow(s1, p1, wdt, x_h%clToDof, ii)

       select case(ft)

       case(fe_p1_1d)  ! line, order 1
          write(999,*) cpt, 1, 0, p1(1:s1)

       case(fe_p2_1d)  ! line, order 2
          write(999,*) cpt, 8, 0, p1(1:s1)

       case(fe_p3_1d)  ! line, order 3
          write(999,*) cpt, 26, 0, p1(1:s1)

       case(fe_p1_2d, fe_p1_2d_disc_ortho&
            & )        ! triangle, order 1
          write(999,*) cpt, 2, 0, p1(1:s1)

       case(fe_p2_2d)  ! triangle, order 2
          write(999,*) cpt, 9, 0, p1(1:s1)

       case(fe_p3_2d)  ! triangle, order 3
          write(999,*) cpt, 21, 0, p1(1:s1)

       case(fe_p1_3d)  ! tetra, order 1
          write(999,*) cpt, 4, 0, p1(1:4)

       case(fe_p2_3d)  ! tetra, order 1
          write(999,*) cpt, 11, 0, p1(1:8), p1(10), p1(9)

       case default
          call quit('feSpace_mod: gmsh_write: cell type not supported')

       end select

    end do
    write(999,fmt="(A)") '$EndElements'

    close(999)

  end subroutine gmsh_write

  subroutine gmsh_addview(X_h, v, fileName, vname, time, idx)
    type(fespace)         , intent(in) :: X_h
    real(RP), dimension(:), intent(in) :: v
    character(len=*)      , intent(in) :: fileName, vname
    real(RP)              , intent(in) :: time
    integer               , intent(in) :: idx

    integer :: ii

    if (choral_verb>1) write(*,*)&
         &"feSpace_mod     : gmsh_addView   = "&
         &//trim(filename) 

    if (size(v,1)/=x_h%nbDof) call quit( &
         & "feSpace_mod: gmsh_addView: fe function 'v' uncorrect" )

    open(unit = 999,file = filename, position='append', status='old')

    write(999,fmt="(A)") '$NodeData'
    write(999,fmt="(A)") '1'
    write(999,fmt="(A)") '"'//trim(vname)//'"'
    write(999,fmt="(A)") '1'
    write(999,*) time
    write(999,fmt="(A)") '3'
    write(999,*) idx
    write(999,fmt="(A)") '1'
    write(999,*) x_h%nbDof
    do ii=1, x_h%nbDof
       write(999,*) ii, v(ii)    
    end do
    write(999,fmt="(A)") '$EndNodeData'
    close(999)

  end subroutine gmsh_addview


  
  subroutine interp_scal_fe(u2, X_h2, u1, X_h1)
    real(RP), dimension(:), intent(out) :: u2
    type(fespace)         , intent(in)  :: X_h1, X_h2
    real(RP), dimension(:), intent(in)  :: u1

    integer, dimension(:), allocatable :: T
    type(mesh), pointer :: m1    =>null()
    type(graph)         :: ndToNd
    real(RP)            :: dist
    integer :: w1, nbDof, dim, ft

    if (choral_verb>0) write(*,*)&
         & "feSpace_mod     : interp_scal_fe"

    if(.NOT.valid(x_h1)) call quit( &
         & "feSpace_mod: interp_scal_fe:&
         & fe space 'X_h1' not valid" )

    if(.NOT.valid(x_h2)) call quit( &
         & "feSpace_mod: interp_scal_fe:&
         & fe space 'X_h2' not valid" )

    if (size(u1,1)/=x_h1%nbDof) call quit( &
         & "feSpace_mod: interp_scal_fe:&
         & fe function 'u1' uncorrect")

    if (size(u2,1)/=x_h2%nbDof) call quit( &
         & "feSpace_mod: interp_scal_fe:&
         & fe function 'u2' uncorrect")

    !! Srortcut to m1
    !!
    m1 => x_h1%mesh

    !! Node to node connectivity for m1
    !!
    call graph_prod(ndtond, x_h1%mesh%NdToCl, x_h1%mesh%clToNd)

    !! Search the closest node of m1
    !! for all dof
    !!
    call closest_node(t, dist, x_h2%dofCoord, m1, ndtond)


    !! Locate each dof of X_h2
    !! inside a cell of m1
    !!
    call locate_nodes(t, x_h2%dofCoord, x_h1, ndtond)
    call clear(ndtond)

    !! Conversion
    !!
    w1 = m1%clToNd%width
    do ft=1, fe_tot_nb
       if (x_h1%fe_count(ft)==0) cycle

       nbdof = fe_nbdof(ft)
       dim   = fe_dim(ft)
       
       call dof_loop()

    end do
    m1    =>null()
    
  contains

    subroutine dof_loop()

      integer , dimension(w1)    :: p1
      integer , dimension(nbDof) :: p2      
      real(RP), dimension(nbDof) :: u1_loc, val
      real(RP), dimension(3,w1)  :: Y
      real(RP), dimension(3)     :: x,z
      integer :: ii, cl1, s1
      
      do ii=1, size(u2,1)

         !! dof is in cell cl1 of the mesh m1
         cl1 = t(ii)
         if ( x_h1%feType(cl1) /= ft ) cycle

         !! dof coordinates
         x = x_h2%dofCoord(:,ii)

         !! p1(1:s1) = nodes of cell cl1
         !!
         call getrow(s1, p1, w1, m1%clToNd, cl1)
            
         !! Cell node coordinates
         !!
         call getndcoord(y, w1, p1(1:s1), s1, m1)

         !! Local coordinates
         !!
         call geotsf_t_inv(z, x, y(1:3,1:s1), s1, m1%clType(cl1) )

         !! values of u1 at the nodes in cl1
         !!
         call getrow(p2, nbdof, x_h1%clToDof, cl1)
         u1_loc = u1(p2)

         select case(x_h2%dofType(ii))
         case(fe_dof_lag)
         
            ! values of the base functions at z 
            call scal_fe(val, nbdof, z(1:dim), dim, ft)
            
            u2(ii) = sum( val*u1_loc )
            
         case default
            call quit( "feSpace_mod: interp_scal_fe:&
                 & DOF type not supported")
            
         end select
      end do

    end subroutine dof_loop
       
  end subroutine interp_scal_fe



  subroutine eval_scal_fe_at_points(val, P, u, X_h)
    real(RP), dimension(:)  , allocatable :: val
    type(fespace)           , intent(in)  :: X_h
    real(RP), dimension(:)  , intent(in)  :: u
    real(RP), dimension(:,:), intent(in)  :: P

    integer, dimension(:), allocatable :: T
    type(mesh), pointer :: m =>null()
    type(graph)         :: ndToNd
    real(RP) :: dist
    integer  :: w1, nbDof, dim, ft

    if (choral_verb>0) write(*,*)&
         & "feSpace_mod     : eval_scal_fe_at_points"

    if(.NOT.valid(x_h))         call quit( &
         & "feSpace_mod: eval_scal_fe_at_points:&
         & fe space 'X_h' not valid" )

    if (size(u,1)/=x_h%nbDof)   call quit( &
         & "feSpace_mod: eval_scal_fe_at_points:&
         & fe function 'u' uncorrect")

    if (size(p,1)/=3 )           call quit( &
         & "feSpace_mod: eval_scal_fe_at_points:&
         & argument 'P' must have shape 3 x nb_points")

    !! allocation for val
    !!
    call allocmem(val, size(p,2))

    !! alias for X_h%mesh
    !!
    m => x_h%mesh

    !! Node to node connectivity for m
    !!
    call graph_prod(ndtond, m%NdToCl, m%clToNd)

    !! Search the cosest node of m
    !! for all points P(:,i)
    !!
    call closest_node(t, dist, p, m, ndtond)


    !! Locate each point P(:,i)
    !! inside a cell of X_h%mesh with non void
    !! feType
    !!
    call locate_nodes(t, p, x_h, ndtond)
    call clear(ndtond)

    !! Conversion
    !!
    w1 = m%clToNd%width
    do ft=1, fe_tot_nb
       if (x_h%fe_count(ft)==0) cycle

       nbdof = fe_nbdof(ft)
       dim   = fe_dim(ft)
       
       call dof_loop()

    end do
    m    =>null()
    
    contains

    subroutine dof_loop()

      integer , dimension(w1)    :: p1
      integer , dimension(nbDof) :: p2      
      real(RP), dimension(nbDof) :: u_loc, u_base_val
      real(RP), dimension(3,w1)  :: Y
      real(RP), dimension(3)     :: x,z
      integer :: ii, cl, s1
      
      do ii=1, size(val,1)

         !! point coordinates
         x = p(:,ii)

         !! the point P(:,i) is in the cell cl of the mesh m
         cl = t(ii)
         if ( x_h%feType(cl) /= ft ) cycle

         !! p1(1:s1) = nodes of cell cl
         !!
         call getrow(s1, p1, w1, m%clToNd, cl)
            
         !! Cell node coordinates
         !!
         call getndcoord(y, w1, p1(1:s1), s1, m)

         !! z = local coordinates of P(:,i) in the ref. cell
         !!
         call geotsf_t_inv(z, x, y(1:3,1:s1), s1, m%clType(cl) )

         !! values of u at the finite element nodes 
         !! in the cell cl1
         !!
         call getrow(p2, nbdof, x_h%clToDof, cl)
         u_loc = u(p2)
      
         ! values of the base functions at z 
         call scal_fe(u_base_val, nbdof, z(1:dim), dim, ft)
         
         val(ii) = sum( u_base_val*u_loc )
      end do

    end subroutine dof_loop
       
  end subroutine eval_scal_fe_at_points



  subroutine locate_nodes(T, P, X_h, ndToNd)
    integer , dimension(:)  , intent(inout):: T
    real(RP), dimension(:,:), intent(in)   :: P
    type(fespace)           , intent(in)   :: X_h
    type(graph)             , intent(in)   :: ndToNd   

    type(mesh), pointer :: m=>null()
    integer :: nbPts, w1, w2, w3, w4, cpt1, cpt2, cpt3
    real(RP) :: a
    
    if (choral_verb>0) write(*,*) &
         & 'feSpace_mod     : locate_nodes' 

    if(.NOT.valid(x_h)) call quit( &
         & "feSpace_mod: locate_nodes:&
         & fe space 'X_h' not valid" )

    if (size(p,1)/= 3) call quit( &
         & "feSpace_mod: locate_nodes:&
         & array 'P' first dimension must be 3" )

    !! Number of points
    !!
    nbpts=size( p, 2)

    if (size(t,1)/=nbpts) call quit( &
         & "feSpace_mod: locate_nodes:&
         & node tab 'T' uncorrect" )

    m => x_h%mesh

    cpt1 = 0
    cpt2 = 0
    cpt3 = 0
    
    w1 =  m%ndToCl%width
    w2 =  m%clToNd%width
    w3 =  ndtond%width
    w4 =  w1*w3

    call loop()

    a = re(cpt1, nbpts)*100._rp
    if (choral_verb>1) write(*,*) &
         & "  locate_nodes, 1st search       =",&
         & real(a, SP), '  % '

    if (cpt2>0) then
       a = re(cpt2, nbpts)*100._rp   
       if (choral_verb>1) write(*,*) &
            &"  locate_nodes, 2nd search       =",&
            & real(a, SP), '  % '
    end if

    if (cpt3>0) then
       a = re(cpt3, nbpts)*100._rp   
       if (choral_verb>1) write(*,*) &
            & "  locate_nodes, 3rd search       =",&
            & real(a, SP), '  % '
    end if

  contains

    subroutine loop()
      integer , dimension(  w1) :: p1
      integer , dimension(  w2) :: p2
      integer , dimension(  w3) :: p3
      integer , dimension(  w4) :: p4, bf1
      real(RP), dimension(3,w1) :: Y
 
      integer :: s1, s2, s3, s4
    
      integer  :: ii, jj, cl
      logical  :: bool
      real(RP), dimension(3) :: x
    
      do ii=1, nbpts
         x = p(:,ii)
         bool = .false.
         !!
         !! 1st search among the neighbour cells of node T(ii)
         !!
         !! p1(1:s1) = neighbour cells of the node T(ii)
         !!
         call getrow(s1, p1, w1, m%ndToCl, t(ii) )

         !! Loop on the neighbour cells
         !!
         do jj=1, s1
            cl = p1(jj)
            if (x_h%feType(cl) == fe_none) cycle
            
            !! p2(1:s2) = nodes of cell cl
            !!
            call getrow(s2, p2, w2, m%clToNd, cl )
            
            !! Cell node coordinates
            !!
            call getndcoord(y, w1, p2(1:s2), s2, m)
            
            bool = belongstocell(x, y(1:3, 1:s2), s2, m%clType(cl))
            
            if (bool) then
               cpt1 = cpt1 + 1
               exit
            end if
         end do
             
         if (.NOT.bool) then
            
            !!
            !! 2nd search among the neighbour of the neighbour cells
            !!
            !! p3(1:s3) = nodes connected with  the node T(ii)
            !!
            call getrow(s3, p3, w3, ndtond, t(ii) )
            !!
            !! p4(1:s4) = neighbour cells to nodes in p3(1:s3)
            !!
            call union(s4, p4, bf1, w4, p1, w1, &
                 &     m%ndToCl, p3(1:s3), s3)

            !!
            !! Loop on the neighbour cells
            !!
            do jj=1, s4
               cl = p4(jj)
               if (x_h%feType(cl) == fe_none) cycle

               !! p2(1:s2) = nodes of cell cl
               !!
               call getrow(s2, p2, w2, m%clToNd, cl )
               
               !! Cell node coordinates
               !!
               call getndcoord(y, w1, p2(1:s2), s2, m)

               bool = belongstocell(x, y(1:3, 1:s2), s2, m%clType(cl))
          
               if (bool) then
                  cpt2 = cpt2 + 1
                  exit
               end if
            end do

         end if

         if (.NOT.bool) then
            cpt3 = cpt3 + 1
            !!
            !! 3rd search : choose a close cell
            !!
            cl = closest_cell(x, x_h)
         end if

         t(ii) = cl

      end do

    end subroutine loop

  end subroutine locate_nodes
  

  function closest_cell(x, X_h) result(cl)
    real(RP), dimension(3), intent(in) :: x
    type(fespace)         , intent(in) :: X_h
    integer                            :: cl

    integer , dimension(  CELL_MAX_NBNODES) :: p1
    real(RP), dimension(3,CELL_MAX_NBNODES) :: Y
    integer  :: ii, jj, nbVtx, s1
    real(RP) :: a, dist
    type(mesh), pointer :: m=>null()

    m    => x_h%mesh
    cl   =  0
    dist =  1e127_rp

    
    do ii=1, m%nbCl
          
       if (x_h%feType(ii) == fe_none) cycle

       !! p1(1:s1) = nodes of cell cl1
       !!
       call getrow(s1, p1, cell_max_nbnodes, m%clToNd, ii)
            
       !! Cell node coordinates
       !!
       call getndcoord(y, cell_max_nbnodes, p1(1:s1), s1, m)

       a  = 0._rp
       jj = m%clType(ii)
       nbvtx = cell_nbvtx(jj)
       do jj=1, nbvtx
          a = a + sum(abs( y(:,jj) - x ))
       end do
       a = a / re(nbvtx)
       if (a<dist) then
          dist = a
          cl = ii
       end if
        
    end do

    m=>null()

  end function closest_cell



  subroutine fespace_assemble(X_h)
    type(fespace), intent(inout) :: X_h

    integer, dimension(:,:), allocatable :: dofTab
    integer, dimension(X_h%mesh%nbCl)    :: cl_row

    integer :: wdt

    if (choral_verb>0) write(*,*)&
         & "feSpace_mod     : assemble"

    if ( sum(x_h%feType(:))==0 ) call quit( &
         'feSpace_mod: assemble: feType not defined' )

    call freemem(x_h%dofType)
    call freemem(x_h%dofCoord)
    x_h%dof_count = 0
    x_h%fe_count  = 0

    call initialise_doftab(doftab, cl_row, x_h)

    wdt = x_h%mesh%ndToCl%width
    call scan_dof(doftab, x_h, wdt)

    call build_cltodof(doftab, cl_row, x_h)
    call freemem( doftab )    

    call x_h_dofcoord(x_h)

    if (choral_verb>1) write(*,*) &
         & "  Number of DOF                  =", x_h%nbDof

    if (.NOT.valid(x_h)) call quit(&
         & 'feSpace_mod: feSpace_assemble: assembling not valid' )

  end subroutine fespace_assemble
  

  subroutine x_h_dofcoord(X_h)
    type(fespace), intent(inout) :: X_h

    type(mesh), pointer :: m=>null()
    type(geotsf)      :: g
    integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
    integer , dimension(FE_MAX_NBDOF)       :: p_dof    
    real(RP), dimension(3, CELL_MAX_NBNODES):: X 
    integer :: ft, fe_typ, cl, sz, ll, dof
    integer :: dim, nbDof
    
    call allocmem( x_h%dofCoord, 3, x_h%nbDof)

    m=>x_h%mesh
    do ft =1, fe_tot_nb

       if( x_h%fe_count(ft) == 0 ) cycle

       dim   = fe_dim(ft)
       nbdof = fe_nbdof(ft)
       g = geotsf(dim, nbdof, fe_dof_coord(ft)%y )

       do cl=1, x_h%mesh%nbCl
          fe_typ = x_h%feType(cl)
          
          if (ft /= fe_typ) cycle

          ! cell node coordinates
          ! sz = nbNodes
          call getrow(sz, p_cl, cell_max_nbnodes, m%clToNd, cl)
          do ll=1, sz
             x(1:3,ll) = m%nd( 1:3, p_cl(ll) )          
          end do

          ! dof coordinates
          call assemble(g, m%clType(cl), x(1:3, 1:sz), sz)

          ! dof indexes
          call getrow(sz, p_dof, fe_max_nbdof, x_h%clTodof, cl)

          do ll=1, sz
             
             dof = p_dof(ll)
             x_h%dofCoord(:,dof) = g%Ty(:,ll)
             
          end do
                    
       end do
    end do
    m=>null()

  end subroutine x_h_dofcoord


  subroutine initialise_doftab(dofTab, cl_row, X_h)
    integer, dimension(:,:), allocatable, intent(inout) :: dofTab
    integer, dimension(:), intent(out) :: cl_row
    type(fespace), intent(inout) :: X_h

    integer, dimension(X_h%mesh%nbCl) :: nnz
    integer :: ii, cpt, fe_typ, nn

    !! The array dofTab contains all local dof
    !!
    !! jj =  number of dof of cell ii
    !!
    !! nnz(ii) = jj
    !!
    !! the rows in dofTab associated with the cell ii
    !! are : cl_row(ii) + 1 ... cl_row(ii) + jj
    !!
    cpt   = 0
    do ii = 1, x_h%mesh%nbCl
       fe_typ = x_h%feType(ii)

       x_h%fe_count(fe_typ) = x_h%fe_count(fe_typ) + 1

       if (fe_typ==0) then
          nnz( ii ) = 0

       else
          cl_row( ii ) = cpt
          nn  = fe_nbdof(fe_typ)
          cpt = cpt + nn
          nnz( ii ) = nn

       end if
    end do

    x_h%clToDof = graph(nnz, size(nnz,1), 0)
    call allocmem(doftab, 3, cpt  )
    doftab = 0

  end subroutine initialise_doftab


  subroutine scan_dof(dofTab, X_h, wdt)
    integer, dimension(:,:), intent(inout) :: dofTab
    type(fespace), intent(inout) :: X_h
    integer      , intent(in)    :: wdt

    type(mesh), pointer :: m=>null()
    integer , dimension(wdt, CELL_MAX_NBNODES):: vtx_toCl
    integer , dimension(wdt) :: dof1_cl, temp
    integer , dimension(CELL_MAX_NBNODES):: cl1_vtx, cl2_nd
    integer , dimension(CELL_MAX_NBNODES):: vtx_nbCl 
    integer , dimension(CELL_MAX_NBVTX)  :: dof1_vtx, dof2_vtx 
    integer , dimension(CELL_MAX_NBVTX)  :: dof1_vtx_g, dof2_vtx_g 
    real(RP), dimension(CELL_MAX_NBVTX)  :: dof1_bary, dof2_bary
    integer , dimension(CELL_MAX_NBVTX)  :: sort1, sort2
    integer , dimension(2) :: p_3
    integer :: p_2

    integer :: dof_cpt, ii, jj, ll, nn, vtx
    integer :: cl1, cl1_nbVtx
    integer :: cl2, cl2_nbNd, cl2_ind
    integer :: dof1, dof1_type, dof1_geo
    integer :: dof1_nbCl
    integer :: dof2, dof2_type, dof2_geo
    integer :: geo1, ft1, ft2
    logical :: bool

    m=>x_h%mesh
    !!    
    !! Loop on the mesh cells
    !!
    dof_cpt   = 0
    do cl1 = 1, m%nbCl
       ft1 = x_h%feType( cl1 )
       if ( ft1 == fe_none ) cycle

       !! cell cl1 geometrical type
       geo1 = fe_geo( ft1 )

       !! cl1_nbVtx  = number of vertexes of cell cl1
       !! cl1_vtx(:) = vertexes of cl1 (global indexzs)
       !!
       call getrow(cl1_nbvtx, cl1_vtx, cell_max_nbnodes, &
            &         m%clToNd, cl1)
       !!
       !! We only consider the vertexes
       !!
       cl1_nbvtx = cell_nbvtx( geo1 )


       !! For the vertex ii of cell cl1 :
       !!   vtx_nbCl(ii)   = number of neighbour cells 
       !!   vtx_toCl(:,ii) = neighbour cells
       !!
       do ii=1, cl1_nbvtx
          call getrow(vtx_nbcl(ii), vtx_tocl(:,ii), wdt, &
               & m%ndToCl, cl1_vtx(ii))
       end do

       !! For the vertex ii of cell cl1 :
       !!   Only consider the neighbour cells with index < cl1
       !!
       !!   vtx_nbCl(ii)   = number of neighbour cells
       !!                     with index < cl1
       !!
       !! IMPORTANT : the graph m%ndCl is sorted ! (see transpose) 
       !!
       do ii=1, cl1_nbvtx
          inner_loop: do jj=1, vtx_nbcl(ii)
             ll = vtx_tocl(jj,ii)
             if ( ll == cl1 ) then
                vtx_nbcl(ii) = jj-1
                exit inner_loop
             end if
          end do inner_loop
       end do

       !! Loop on the dof of cell cl1
       !!
       dof_loop: do dof1=1, fe_nbdof(ft1)
          dof_cpt = dof_cpt + 1
          dof1_geo  = fe_dof_geo(dof1, ft1)
          dof1_type = fe_dof_type(dof1, ft1)


          select case(dof1_geo)

          case(cell_tet)
             ! new dof
             x_h%dof_count(cell_tet) = x_h%dof_count(cell_tet) + 1
             doftab(3, dof_cpt) = dof1_type

          case(cell_vtx)   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

             !! vertex local index in cl1
             !!
             dof1_vtx(1) = fe_dof_vtx( 1, dof1, ft1)

             !!   dof1_nbCl  = number of cells associated with dof1
             !!   dof1_cl(:) = cells associated with dof1
             !!                with index < cl1
             !!
             vtx = dof1_vtx(1)
             dof1_nbcl = vtx_nbcl(vtx)
             if (dof1_nbcl==0) goto 10   ! new dof
             dof1_cl( 1: dof1_nbcl) = vtx_tocl( 1: dof1_nbcl, vtx)

             !!   dof1_vtx_g(:) = vtx associated with dof1         
             !!            global indexes in the mesh
             !!
             dof1_vtx_g(1) = cl1_vtx( dof1_vtx(1) )

             !! Search if dof1 also is a dof of a cell
             !! in dof1_cl( 1: dof1_nbCl)
             !!
             ngb_cl_a: do cl2_ind=1, dof1_nbcl

                cl2 = dof1_cl(cl2_ind)
                ft2 = x_h%feType(cl2)
                if ( ft2 == fe_none ) cycle ngb_cl_a

                !! cl2_nbNd  = number of nodes of cell cl2
                !! cl2_nd(:) = nodes of cell cl2 (global indexzs)
                !!
                call getrow(cl2_nbnd, cl2_nd, &
                     & cell_max_nbnodes, m%clToNd, cl2)

                !! Loop on the dof of cl2
                !!
                dof2_a: do dof2=1, fe_nbdof(ft2)

                   dof2_geo = fe_dof_geo(dof2, ft2)
                   if ( dof2_geo /= dof1_geo) cycle dof2_a 

                   !!  dof2_vtx(1) = vtx associated with dof2
                   !!                (local indexes in cl2)
                   !!
                   dof2_vtx( 1 ) = fe_dof_vtx( 1, dof2, ft2)

                   !!   dof2_vtx_g(1) = vtx associated with dof1
                   !!            global indexes in mesh
                   dof2_vtx_g(1) = cl2_nd( dof2_vtx(1) )

                   !! We verify if dof1 and dof2 are associated with
                   !! the same vertex in the mesh
                   !!
                   bool = ( dof1_vtx_g(1) == dof2_vtx_g(1) )
                   if ( .NOT. bool ) cycle dof2_a  

                   !! dof type comparison
                   !!
                   dof2_type = fe_dof_type(dof2, ft2)
                   if ( dof2_type /= dof1_type) cycle dof2_a

                   ! dof1 == dof2 : dof1 is not a new dof        
                   doftab(1:2, dof_cpt) = (/ cl2, dof2 /)
                   cycle dof_loop

                end do dof2_a

             end do ngb_cl_a
10           continue
             ! new dof
             x_h%dof_count(cell_vtx) = x_h%dof_count(cell_vtx) + 1
             doftab(3, dof_cpt) = dof1_type
             
          case(cell_edg)   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

             !! vertex local index in cl1
             !!
             dof1_vtx(1:2) = fe_dof_vtx( 1:2, dof1, ft1)

             !! barycentric coord
             !!
             dof1_bary(1:2) = fe_dof_bary(1:2, dof1, ft1)

             !!   dof1_cl(:) = cells associated with dof1
             !!                with index < cl1
             !!
             !! vtx 1
             ii = dof1_vtx(1)
             jj = dof1_vtx(2)

             ll = vtx_nbcl(ii)
             if (ll==0) goto 20   ! new dof

             nn = vtx_nbcl(jj)
             if (nn==0) goto 20   ! new dof

             call cap_sorted_set(dof1_nbcl, dof1_cl, wdt, &
                  & vtx_tocl( 1: ll, ii), ll, &
                  & vtx_tocl( 1: nn, jj), nn)

             if (dof1_nbcl==0) goto 20       ! new dof


             !!   dof1_vtx_g(1:2) = vtx associated with dof1         
             !!            global indexes in the mesh
             !!            these indexes are sorted
             !!
             dof1_vtx_g(1:2) = cl1_vtx( dof1_vtx(1:2) )
             sort1(1:2) = dof1_vtx_g(1:2)
             call order_2(sort1(1:2))

             !! Search if dof1 also is a dof of a cell
             !! in dof1_cl( 1: dof1_nbCl)
             !!
             dof1_type    = fe_dof_type(dof1, ft1)

             ngb_cl_b: do cl2_ind=1, dof1_nbcl

                cl2 = dof1_cl(cl2_ind)
                ft2 = x_h%feType(cl2)
                if ( ft2 == fe_none ) cycle ngb_cl_b

                !! cl2_nbNd  = number of nodes of cell cl2
                !! cl2_nd(:) = nodes of cell cl2 (global indexzs)
                !!
                call getrow(cl2_nbnd, cl2_nd, &
                     & cell_max_nbnodes, m%clToNd, cl2)

                !! Loop on the dof of cl2
                !!
                dof2_b: do dof2=1, fe_nbdof(ft2)

                   dof2_geo = fe_dof_geo(dof2, ft2)
                   if ( dof2_geo /= dof1_geo) cycle dof2_b 

                   !!  dof2_vtx(1:2) = vtx associated with dof2
                   !!                  (local indexes in cl2)
                   !!
                   dof2_vtx(1:2) = fe_dof_vtx( 1:2, dof2, ft2)

                   !!   dof2_vtx_g(1) = vtx associated with dof1
                   !!            global indexes in mesh
                   dof2_vtx_g(1:2) = cl2_nd( dof2_vtx(1:2) )
                   sort2(1:2) = dof2_vtx_g(1:2)
                   call order_2(sort2(1:2) )

                   !! We verify if dof1 and dof2 are associated with
                   !! the same vertexes in the mesh
                   !!
                   bool = all( sort1(1:2) == sort2(1:2) )

                   if ( .NOT. bool ) cycle dof2_b  

                   !! Verify barycentric coordinates equality
                   !!
                   dof2_bary(1:2) = fe_dof_bary(1:2, dof2, ft2)

                   !! permutation between the two vertex
                   !! triplet
                   !!
                   call get_perm_2(p_2, dof1_vtx_g(1:2)&
                        &             , dof2_vtx_g(1:2) )

                   !! barycentric coord. equality
                   !!
                   bool = equal( dof1_bary(1), dof2_bary( p_2 ) )
                   if ( .NOT. bool ) cycle dof2_b  

                   !! dof type comparison
                   !!
                   dof2_type = fe_dof_type(dof2, ft2)
                   if ( dof2_type /= dof1_type) cycle dof2_b

                   ! dof1 == dof2 : dof1 is not a new dof        
                   doftab(1:2, dof_cpt) = (/ cl2, dof2 /)
                   cycle dof_loop

                end do dof2_b

             end do ngb_cl_b

20           continue

             ! new dof
             x_h%dof_count(cell_edg) = x_h%dof_count(cell_edg) + 1
             doftab(3, dof_cpt) = dof1_type

          case(cell_trg)   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

             !! vertex local index in cl1
             !!
             dof1_vtx(1:3) = fe_dof_vtx( 1:3, dof1, ft1)

             !! barycentric coord
             !!
             dof1_bary(1:3) = fe_dof_bary(1:3, dof1, ft1)

             !!   dof1_cl(:) = cells associated with dof1
             !!                with index < cl1
             !!
             !! vtx 1
             vtx = dof1_vtx(1)
             dof1_nbcl = vtx_nbcl(vtx)
             if (dof1_nbcl==0) goto 30   ! new dof

             dof1_cl( 1: dof1_nbcl) = vtx_tocl( 1: dof1_nbcl, vtx)

             !! vtx 2, 3
             do ii=2, 3     
                vtx = dof1_vtx(ii)
                jj  = vtx_nbcl(vtx)
                if (jj==0) goto 30   ! new dof

                call cap_sorted_set(ll, temp, wdt, &
                     & dof1_cl( 1: dof1_nbcl), dof1_nbcl, &
                     & vtx_tocl( 1: jj, vtx), jj)

                if (ll==0) goto 30       ! new dof

                dof1_nbcl = ll
                dof1_cl( 1: ll) = temp( 1: ll)

             end do

             !!   dof1_vtx_g(1:3) = vtx associated with dof1         
             !!            global indexes in the mesh
             !!            these indexes are sorted
             !!
             dof1_vtx_g(1:3) = cl1_vtx( dof1_vtx(1:3) )
             sort1(1:3) = dof1_vtx_g(1:3)
             call sort(sort1(1:3), 3)

             !! Search if dof1 also is a dof of a cell
             !! in dof1_cl( 1: dof1_nbCl)
             !!
             dof1_type    = fe_dof_type(dof1, ft1)

             ngb_cl_c: do cl2_ind=1, dof1_nbcl

                cl2 = dof1_cl(cl2_ind)
                ft2 = x_h%feType(cl2)
                if ( ft2 == fe_none ) cycle ngb_cl_c

                !! cl2_nbNd  = number of nodes of cell cl2
                !! cl2_nd(:) = nodes of cell cl2 (global indexzs)
                !!
                call getrow(cl2_nbnd, cl2_nd, &
                     & cell_max_nbnodes, m%clToNd, cl2)

                !! Loop on the dof of cl2
                !!
                dof2_c: do dof2=1, fe_nbdof(ft2)

                   dof2_geo = fe_dof_geo(dof2, ft2)
                   if ( dof2_geo /= dof1_geo) cycle dof2_c 

                   !!  dof2_vtx(1:3) = vtx associated with dof2
                   !!                  (local indexes in cl2)
                   !!
                   dof2_vtx(1:3) = fe_dof_vtx( 1:3, dof2, ft2)

                   !!   dof2_vtx_g(1) = vtx associated with dof1
                   !!            global indexes in mesh
                   !!            these indexes are sorted
                   dof2_vtx_g(1:3) = cl2_nd( dof2_vtx(1:3) )
                   sort2(1:3) = dof2_vtx_g(1:3)
                   call sort(sort2(1:3), 3)

                   !! We verify if dof1 and dof2 are associated with
                   !! the same vertexes in the mesh
                   !!
                   bool = all( sort1(1:3) == sort2(1:3) )
                   if ( .NOT. bool ) cycle dof2_c  

                   !! Verify barycentric coordinates equality
                   !!
                   dof2_bary(1:3) = fe_dof_bary(1:3, dof2, ft2)

                   !! permutation between the two vertex
                   !! triplet
                   !!
                   call get_perm_3(p_3, dof1_vtx_g(1:3)  &
                        &             ,  dof2_vtx_g(1:3) )

                   !! barycentric coord. equality
                   !!
                   bool = equal( dof1_bary(1), dof2_bary( p_3(1) ) )
                   bool = bool .AND. &
                        & equal( dof1_bary(2), dof2_bary( p_3(2) ) )

                   if ( .NOT. bool ) cycle dof2_c  

                   !! dof type comparison
                   !!
                   dof2_type = fe_dof_type(dof2, ft2)
                   if ( dof2_type /= dof1_type) cycle dof2_c

                   ! dof1 == dof2 : dof1 is not a new dof        
                   doftab(1:2, dof_cpt) = (/ cl2, dof2 /)
                   cycle dof_loop

                end do dof2_c

             end do ngb_cl_c

30           continue
             ! new dof
             x_h%dof_count(cell_trg) = x_h%dof_count(cell_trg) + 1
             doftab(3, dof_cpt) = dof1_type

          case default

             call quit( 'feSpace_mod: scan_dof:&
                  & cell type not supported')
          end select


       end do dof_loop

    end do
    m=>null()
  end subroutine scan_dof


  subroutine build_cltodof(dofTab, cl_row, X_h)
    integer, dimension(:,:), intent(inout) :: dofTab
    integer, dimension(:)  , intent(in)    :: cl_row
    type(fespace)          , intent(inout) :: X_h

    integer :: cpt, ii, jj, cl, dof, fe_typ, j1, j2 
    
    !! counts all new dof
    !! give to all new dof an index
    !!
    cpt   = 0
    do ii = 1, size( doftab, 2)
       if ( doftab(1, ii) == 0 ) then
          cpt = cpt + 1
          doftab(1, ii) = cpt
       end if
    end do

    !! Number of DOF
    if (cpt/=sum( x_h%dof_count)) call quit( &
         & 'feSpace_mod: build_clToDof: 1')
    if (cpt == 0 )                call quit( &
         & 'feSpace_mod: build_clToDof: 2')
    x_h%nbDof      = cpt
    x_h%clToDof%nc = cpt

    !! DOF types
    call allocmem(x_h%dofType, x_h%nbDof)
    do ii = 1, size( doftab, 2)
       if ( doftab(2, ii) == 0 ) then
          dof = doftab(1, ii)
          x_h%dofType(dof) = doftab(3, ii)
       end if
    end do

    !!  finds the new dof equal to a local dof
    !!  when it is not a new dof
    !!
    do ii = 1, size( doftab, 2)
       
       if ( doftab(2, ii) == 0 ) cycle

       jj = ii
       do while( doftab(2, jj) /= 0 )

          cl  = doftab(1, jj)
          dof = doftab(2, jj)
          
          jj = cl_row( cl ) + dof

       end do
       doftab(1, ii) = doftab(1, jj)

    end do

    !!  Fills in the cell to dof graph
    !!
    do cl = 1, x_h%mesh%nbCl
       fe_typ = x_h%feType(cl)

       if (fe_typ==0) cycle

       jj = fe_nbdof( fe_typ )

       j1 = cl_row(cl) + 1
       j2 = cl_row(cl) + jj

       call setrow(x_h%clToDof, doftab(1, j1:j2), jj, cl)
       
    end do

  end subroutine build_cltodof

  subroutine get_perm_3(p, t1, t2)
    integer, dimension(2), intent(out) :: p
    integer, dimension(3), intent(in)  :: t1, t2

    if (t1(1)==t2(1)) then
       p(1) = 1
       if (t1(2)==t2(2)) then
          p(2) = 2
       else
          p(2) = 3
       end if

    else if (t1(1)==t2(2)) then
       p(1) = 2
       if (t1(2)==t2(1)) then
          p(2) = 1
       else
          p(2) = 3
       end if

    else
       p(1) = 3
       if (t1(2)==t2(2)) then
          p(2) = 2
       else
          p(2) = 1
       end if

    end if

  end subroutine get_perm_3

  subroutine get_perm_2(p, t1, t2)
    integer, intent(out) :: p
    integer, dimension(2), intent(in)  :: t1, t2

    if (t1(1)==t2(1)) then
       p = 1

    else
       p = 2

    end if

  end subroutine get_perm_2

  subroutine interp_scal_func(uh, u, X_h)
    real(RP), dimension(:), allocatable :: uh
    type(fespace)   , intent(in) :: X_h
    procedure(r3tor)             :: u

    integer :: ii

    if (choral_verb>1) write(*,*) &
         & "feSpace_mod     : interp_scal_func"

    if (.NOT.valid(x_h)) call quit(&
         &  "feSpace_mod: interp_scal_func:&
         & fe space 'X_h' not valid" )

    call allocmem(uh, x_h%nbDof)

    do ii=1, x_h%nbDof
       select case (x_h%dofType(ii))
       case(fe_dof_lag)
          uh(ii) = u(x_h%dofCoord(:,ii))
       case default
          call quit("feSpace_mod: interp_scal_func:&
               & DOF type not supported" )
       end select
    end do

  end subroutine interp_scal_func



  subroutine fespace_interp_vect_func(phi_h, phi, X_h)
    real(RP), dimension(:), allocatable :: phi_h
    type(fespace)         , intent(in)  :: X_h
    procedure(r3tor3)                   :: phi

    type(mesh), pointer :: m
    integer  :: dim, nbDof, ft, nbItf

    if (choral_verb>1) write(*,*) &
         & "feSpace_mod     : interp_vect_func"

    m=>x_h%mesh

    if (m%nbItf<=0) call quit("feSpace_mod:&
         & feSpace_interp_vect_func: mesh interfaces not defined" )

    call allocmem(phi_h, x_h%nbDof)
    phi_h = 0._rp

    do ft=1, fe_tot_nb
       if (x_h%fe_count(ft)==0) cycle

       dim   = fe_dim(ft)          ! dimension of K_ref
       nbdof = fe_nbdof(ft)        ! fe number of DOF
       nbitf = cell_nbitf( fe_geo( ft) )  ! interface number

       call cell_loop()

    end do

  contains

    subroutine cell_loop()

      integer , dimension(nbDof)               :: p       
      integer , dimension(CELL_MAX_NBNODES)    :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES) :: X 
      real(RP), dimension(nbDof)               :: val
      integer , dimension(CELL_MAX_NBITF)      :: eps 

      real(RP), dimension(3) :: vec
      type(geotsf):: g
      integer     :: cl, ii, jj, cl_nd, nb_itf
      real(RP)    :: msK

      real(RP), dimension(3,nbItf) :: nf 
      real(RP), dimension(nbItf)   :: msf


      !! geometricat transformation K_ref --> K
      !
      g = geotsf(dim, nbdof, fe_dof_coord(ft)%y )
      
      !! CELL LOOP
      !!
      do cl=1, m%nbCl

         ii = x_h%feType(cl)
         if (ii/=ft) cycle

         ! cell node coordinates
         call getrow(cl_nd, p_cl, cell_max_nbnodes, m%clToNd, cl)
         do ii=1, cl_nd
            x(1:3,ii) = m%nd( 1:3, p_cl(ii) )          
         end do

         ! transformation T : K_ref --> K = T(K_ref)
         call assemble(g, m%clType(cl), x(1:3, 1:cl_nd), cl_nd)

         ! unit normals and measures of the cell interfaces 
         call cell_ms_itf_n(msk, nf, msf, nbitf, g) 

#if DBG>0
         !!
         !! checks the cell orientation
         if (.NOT.cell_orientation(m, cl)) call quit(&
         &  "feSpace_mod: feSpace_interp_vect_func:&
         & cell orientation incorrect" )
#endif
         !!
         !! orientation of the interfaces of cell cl
         call interface_orientation(nb_itf, &
              & eps, cell_max_nbitf, m, cl)

         select case(ft)

         case(fe_rt0_1d, fe_rt0_2d)

            ! local dof : 
            do ii=1, nbdof
               val(ii) = sum( nf(:,ii)  * phi(g%Ty(:,ii)) )
               val(ii) = val(ii) * msf(ii)

               if (eps(ii)==1) cycle
               val(ii) = - val(ii)
            end do

         case(fe_rt1_1d)

            ! local dof : 
            do ii=1, 2
               val(ii) = sum( nf(:,ii)  * phi(g%Ty(:,ii)) )
               val(ii) = val(ii) * msf(ii)

               if (eps(ii)==1) cycle
               val(ii) = - val(ii)
            end do

            vec = phi(g%Ty(:,3)) 
            val(3) = vec(1)

         case(fe_rt1_2d_2)

            ! local dof :
            do ii=1, 3
               jj = 2*ii - 1

               val(jj) = sum( nf(:,ii)  * phi(g%Ty(:,jj)) )
               val(jj) = val(jj) * msf(ii)

               jj = jj + 1
               val(jj) = sum( nf(:,ii)  * phi(g%Ty(:,jj)) )
               val(jj) = val(jj) * msf(ii)

               if ( eps(ii) == 1 ) cycle 
               val(jj-1:jj) = -val(jj-1:jj)

            end do

            vec = phi(g%Ty(:,7)) * 2.0_rp * msk
            val(7) = sum( vec * nf(:,3) )
            val(8) = sum( vec * nf(:,1) )

            vec    = x(:,2) - x(:,1)  ! AB
            val(7) = val(7) / sum( vec * nf(:,3) )

            vec    = x(:,3) - x(:,1)  ! AC
            val(8) = val(8) / sum( vec * nf(:,1) )

         case default
            call quit("feSpace_mod: feSpace_interp_vect_func:&
                 & fe type not supported")

         end select

         ! global dof indexes
         call getrow(p, nbdof, x_h%clTodof, cl)
         phi_h(p) = val

      end do

    end subroutine cell_loop

  end subroutine fespace_interp_vect_func


  subroutine flag_x_h_dof(dof_flag, X_h, dim, f)
    logical, dimension(:), allocatable :: dof_flag
    type(fespace)        , intent(in)  :: X_h
    integer              , intent(in)  :: dim
    procedure(r3tor)     , optional    :: f

    logical , dimension(:), allocatable   :: flag

    integer, dimension(FE_MAX_NBDOF) :: dof_tab
    integer :: ii, jj, nb_dof

    if (choral_verb>1) write(*,*) &
         & "feSpace_mod     : flag_X_h_dof"

    if (.NOT.valid(x_h)) call quit('feSpace_mod:&
         & flag_X_h_dof: X_h not valid')

    call flag_mesh_cells(flag, x_h%mesh, dim, f)

    call allocmem(dof_flag, x_h%nbDof)
    dof_flag = .false.
    do ii=1, x_h%mesh%nbCl

       !! cell ii is not a boundary cell so that f >=0 
       if (.NOT.flag(ii)) cycle

       !! flag all finite element nodes linked to cell cl
       call getrow(nb_dof, dof_tab, fe_max_nbdof, x_h%clToDof, ii)
       do jj=1, nb_dof
          dof_flag( dof_tab(jj) ) = .true.
       end do

    end do
    
  end subroutine flag_x_h_dof

end module fespace_mod