integral.F90

Go to the documentation of this file.

module integral

  use choral_constants
  use choral_variables
  use real_type, only: rp
  use basic_tools
  use abstract_interfaces, only: r3xrxrtor, r3tor, r3xr3xr3tor, r3tor3
  use algebra_lin
  use cell_mod
  use fe_mod
  use graph_mod
  use mesh_mod
  use fespace_mod
  use quad_mod
  use quadmesh_mod
  use geotsf_mod

  implicit none
  private

  public :: integ                      ! partially tested
  public :: integ_scal_fe_grad_proj
  public :: l2_dist, l2_dist_grad, l2_dist_grad_proj
  public :: l2_dist_vect

  !       %----------------------------------------%
  !       |                                        |
  !       |          GENERIc SUBROUTINES           |
  !       |                                        |
  !       %----------------------------------------%

  interface integ
     module procedure integ_scal_func         ! tested
     module procedure integ_scal_fe
     module procedure integ_scal_fe_grad      ! tested
     module procedure integ_scal_fe_grad_2
  end interface integ

  interface l2_dist
     module procedure l2_dist_scalfe
  end interface l2_dist
  interface l2_dist_grad
     module procedure l2_dist_grad_scalfe
  end interface l2_dist_grad

contains


  function integ_scal_func(u, qdm) result(int)
    real(RP)                   :: int
    procedure(r3tor)           :: u
    type(quadmesh), intent(in) :: qdm

    type(mesh), pointer :: m
    integer  :: qt, dim, nn

    if (choral_verb>1) write(*,*) &
         & 'integral        : integ_scal_func'

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: integ_scal_func:&
         & quad mesh 'qdmh' not valid" )

    int =  0.0_rp
    m   => qdm%mesh

    do qt=1, quad_tot_nb
       if (qdm%quad_count(qt)==0) cycle

       dim   = quad_dim(qt)        ! dimension of K_ref
       nn    = quad_nbnodes(qt)    ! quad number of nodes

       call cell_loop()

    end do

  contains

    subroutine cell_loop()

      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(nn)      :: wgt, u_x
      real(RP), dimension(dim, nn) :: y

      type(geotsf):: g
      integer     :: cl, jj, ll, cl_nd
      real(RP)    :: int2

      int2 = 0._rp

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y

      g = geotsf(dim, nn, y)
      
      ! Loop on mesh cells
      do cl=1, m%nbCl

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

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

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

         ! loop on quadrature nodes
         do ll=1, nn
            u_x(ll) = u(g%Ty(:,ll))
         end do
         u_x  = u_x  * g%Jy
         u_x  = u_x  * wgt
         int2 = int2 + sum(u_x)

      end do

      int = int + int2

    end subroutine cell_loop

  end function integ_scal_func





  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!
  !!    SCALAR FINITE ELEMENTS
  !!
  !!
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  function integ_scal_fe(E, u, uh, X_h, qdm) result(int)
    real(RP)                           :: int
    procedure(r3xrxrtor)               :: E
    procedure(r3tor)                   :: u
    real(RP), dimension(:), intent(in) :: uh
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm

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

    if (choral_verb>1) write(*,*) &
         & 'integral        : integ_scal_fe'

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

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: integ_scal_fe:&
         & quad mesh 'qdmh' not valid" )

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"integral: integ_scal_fe:&
         & 'X_h' and 'qdmh' associated to different meshes" )

    if(size(uh,1) /= x_h%nbDof) call quit(  &
         &"integral: integ_scal_fe:&
         & fe function 'uh' has not the expected size" )

    int = 0._rp

    m => x_h%mesh

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

       do qt=1, quad_tot_nb
          if (qdm%quad_count(qt)==0) cycle
          if (quad_geo(qt) /= fe_geo(ft) ) cycle

          dim   = quad_dim(qt)        ! dimension of K_ref
          nbdof = fe_nbdof(ft)        ! fe number of DOF
          nn    = quad_nbnodes(qt)    ! quad number of nodes

          call cell_loop()

       end do
    end do

  contains

    subroutine cell_loop()

      real(RP), dimension(nbDof, nn)          :: val
      integer , dimension(nbDof)              :: p       
      real(RP), dimension(nbDof)              :: v
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(     nn) :: wgt, E_x
      real(RP), dimension(dim, nn) :: y

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd
      real(RP)    :: int2, uh_x, u_x

      int2 = 0._rp

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y

      !! fe basis functions at quad nodes
      !!
      do ll=1, nn
         call scal_fe(val(:,ll), nbdof, y(:,ll), dim, ft)
      end do

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

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

         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 ll=1, cl_nd
            x(:,ll) = m%nd( :, p_cl(ll) )          
         end do

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

         ! dof associated with cell cl indexes
         call getrow(p, nbdof, x_h%clToDof, cl)

         ! local dof of the fe function uh
         v = uh(p)

         ! loop on quad nodes
         do ll=1, nn
            ! x = T(node ll)

            uh_x   = dot_product(val(:,ll), v)  ! u(x)
            u_x    = u(g%Ty(:,ll))              ! uh(x)
            e_x(ll)= e(g%Ty(:,ll), u_x, uh_x)   ! E(x, u(x), uh(x))

         end do

         ! update integral
         e_x = e_x * g%Jy
         e_x = e_x * wgt
         int2 = int2 + sum(e_x)

      end do

      int = int + int2

    end subroutine cell_loop

  end function integ_scal_fe





  function integ_scal_fe_grad(E, phi, uh, X_h, qdm) result(int)
    real(RP)                             :: int
    procedure(r3xr3xr3tor)               :: E
    real(RP)  , dimension(:), intent(in) :: uh
    type(fespace)           , intent(in) :: X_h
    type(quadmesh)          , intent(in) :: qdm
    procedure(r3tor3)                    :: phi

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

    if (choral_verb>1) write(*,*) &
         & 'integral        : Integ_scal_fe_grad'

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

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: integ_scal_fe_grad:&
         & quad mesh 'qdmh' not valid" )

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"integral: integ_scal_fe_grad:&
         & 'X_h' and 'qdmh' associated to different meshes" )

    if(size(uh,1) /= x_h%nbDof) call quit(  &
         &"integral: integ_scal_fe_grad:&
         & fe function 'uh' has not the expected size" )

    int = 0._rp

    m => x_h%mesh

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

       do qt=1, quad_tot_nb
          if (qdm%quad_count(qt)==0) cycle
          if (quad_geo(qt) /= fe_geo(ft) ) cycle

          dim   = quad_dim(qt)        ! dimension of K_ref
          nbdof = fe_nbdof(ft)        ! fe number of DOF
          nn    = quad_nbnodes(qt)    ! quad number of nodes

          call cell_loop()

       end do
    end do

  contains

    subroutine cell_loop()      
      real(RP), dimension(dim, nbDof, nn)  :: val
      integer , dimension(nbDof)           :: p       
      real(RP), dimension(nbDof)           :: v
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(3):: grad_uh_x, w, phi_x
      real(RP), dimension(     nn) :: wgt
      real(RP), dimension(dim, nn) :: y

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd
      real(RP)    :: int2, E_x

      int2 = 0._rp

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y

      !! gradient of the fe basis functions
      !! evaluated at quad nodes
      !!
      do ll=1, nn
         call scalgrad_fe(val(:,:,ll), nbdof,  y(:,ll), dim, ft)
      end do

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

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

         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 ll=1, cl_nd
            x(:,ll) = m%nd( :, p_cl(ll) )          
         end do

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

         ! dof associated with cell cl indexes
         call getrow(p, nbdof, x_h%clToDof, cl)

         ! local dof of the fe function uh
         v = uh(p)

         ! loop on quad nodes
         do ll=1, nn
            ! x = T(node ll)

            ! compute the gradient
            grad_uh_x(:) = 0._rp
            select case(dim)
            case(3)
               do jj=1, nbdof
                  call matvecprod_3x3( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll))
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            case(2)
               do jj=1, nbdof
                  call matvecprod_3x2( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            case(1)
               do jj=1, nbdof
                  call matvecprod_3x1( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            end select
            grad_uh_x = grad_uh_x / g%Jy(ll)

            ! phi(x)
            phi_x = phi(g%Ty(:,ll))

            ! update integral
            e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
            int2 = int2 + e_x * wgt(ll) * g%Jy(ll)

         end do
      end do

      int = int + int2

    end subroutine cell_loop

  end function integ_scal_fe_grad


  function integ_scal_fe_grad_proj(E, phi, uh, X_h, qdm) result(int)
    real(RP)                             :: int
    procedure(r3xr3xr3tor)               :: E
    real(RP)  , dimension(:), intent(in) :: uh
    type(fespace)           , intent(in) :: X_h
    type(quadmesh)          , intent(in) :: qdm
    procedure(r3tor3)                    :: phi

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

    if (choral_verb>1) write(*,*) &
         & 'integral        : Integ_scal_fe_grad_proj'

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

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: integ_scal_fe_grad_proj:&
         & quad mesh 'qdmh' not valid" )

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"integral: integ_scal_fe_grad_proj:&
         & 'X_h' and 'qdmh' associated to different meshes" )

    if(size(uh,1) /= x_h%nbDof) call quit(  &
         &"integral: integ_scal_fe_grad_proj:&
         & fe function 'uh' has not the expected size" )

    int = 0._rp

    m => x_h%mesh

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

       do qt=1, quad_tot_nb
          if (qdm%quad_count(qt)==0) cycle
          if (quad_geo(qt) /= fe_geo(ft) ) cycle

          dim   = quad_dim(qt)        ! dimension of K_ref
          nbdof = fe_nbdof(ft)        ! fe number of DOF
          nn    = quad_nbnodes(qt)    ! quad number of nodes

          call cell_loop()

       end do
    end do

  contains

    subroutine cell_loop()

      real(RP), dimension(dim, nbDof, nn)  :: val
      integer , dimension(nbDof)           :: p       
      real(RP), dimension(nbDof)           :: v
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(3):: grad_uh_x, w, phi_x
      real(RP), dimension(     nn) :: wgt
      real(RP), dimension(dim, nn) :: y

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd
      real(RP)    :: int2, E_x

      int2 = 0._rp

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y


      !! gradient of the fe basis functions
      !! evaluated at quad nodes
      !!
      do ll=1, nn
         call scalgrad_fe(val(:,:,ll), nbdof,  y(:,ll), dim, ft)
      end do

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

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

         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 ll=1, cl_nd
            x(:,ll) = m%nd( :, p_cl(ll) )          
         end do

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

         ! dof associated with cell cl indexes
         call getrow(p, nbdof, x_h%clToDof, cl)

         ! local dof of the fe function uh
         v = uh(p)

         ! loop on quad nodes
         do ll=1, nn
            ! x = T(node ll)

            ! compute the gradient
            grad_uh_x(:) = 0._rp
            select case(dim)
            case(3)
               do jj=1, nbdof
                  call matvecprod_3x3( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll))
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            case(2)
               do jj=1, nbdof
                  call matvecprod_3x2( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            case(1)
               do jj=1, nbdof
                  call matvecprod_3x1( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh_x = grad_uh_x + v(jj) * w
               end do

            end select
            grad_uh_x = grad_uh_x / g%Jy(ll)

            ! phi(x)
            phi_x = phi(g%Ty(:,ll))

            ! tangential component phi_T(x)
            if (dim==3) then
               ! nothing

            else if (dim==2) then
               w = crossprod_3d( g%DTy(:,1,ll), g%DTy(:,2,ll))
               w = w / norm2(w)
               phi_x = phi_x - sum(phi_x * w)*w

            else if (dim==1) then
               phi_x = sum( g%Cy(:,1,ll)*phi_x ) * g%Cy(:,1,ll)

            end if

            ! update integral
            e_x = e(g%Ty(:,ll), phi_x, grad_uh_x)
            int2 = int2 + e_x * wgt(ll) * g%Jy(ll)

         end do
      end do

      int = int + int2

    end subroutine cell_loop

  end function integ_scal_fe_grad_proj





  function integ_scal_fe_grad_2(E, uh1, uh2, X_h, qdm) result(int)
    real(RP)                             :: int
    procedure(r3xr3xr3tor)               :: E
    real(RP)  , dimension(:), intent(in) :: uh1, uh2
    type(fespace)           , intent(in) :: X_h
    type(quadmesh)          , intent(in) :: qdm

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

    if (choral_verb>1) write(*,*) &
         & 'integral        : Integ_scal_fe_grad_2'

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

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: integ_scal_fe_grad_2:&
         & quad mesh 'qdmh' not valid" )

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"integral: integ_scal_fe_grad_2:&
         & 'X_h' and 'qdmh' associated to different meshes" )

    if(size(uh1,1) /= x_h%nbDof) call quit(  &
         &"integral: integ_scal_fe_grad_2:&
         & fe function 'uh1' has not the expected size" )

    if(size(uh2,1) /= x_h%nbDof) call quit(  &
         &"integral: integ_scal_fe_grad_2:&
         & fe function 'uh2' has not the expected size" )


    int = 0._rp

    m => x_h%mesh

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

       do qt=1, quad_tot_nb
          if (qdm%quad_count(qt)==0) cycle
          if (quad_geo(qt) /= fe_geo(ft) ) cycle

          dim   = quad_dim(qt)        ! dimension of K_ref
          nbdof = fe_nbdof(ft)        ! fe number of DOF
          nn    = quad_nbnodes(qt)    ! quad number of nodes

          call cell_loop()

       end do
    end do

  contains

    subroutine cell_loop()

      real(RP), dimension(dim, nbDof, nn)  :: val
      integer , dimension(nbDof)           :: p       
      real(RP), dimension(nbDof)           :: v1, v2
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(3):: grad_uh1_x, grad_uh2_x, w
      real(RP), dimension(     nn) :: wgt
      real(RP), dimension(dim, nn) :: y

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd
      real(RP)    :: int2, E_x

      int2 = 0._rp

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y


      !! gradient of the fe basis functions
      !! evaluated at quad nodes
      !!
      do ll=1, nn
         call scalgrad_fe(val(:,:,ll), nbdof,  y(:,ll), dim, ft)
      end do

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

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

         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 ll=1, cl_nd
            x(:,ll) = m%nd( :, p_cl(ll) )          
         end do

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

         ! dof associated with cell cl indexes
         call getrow(p, nbdof, x_h%clToDof, cl)

         ! local dof of the fe function uh
         v1 = uh1(p)
         v2 = uh2(p)

         ! loop on quad nodes
         do ll=1, nn
            ! x = T(node ll)

            ! compute the gradient
            grad_uh1_x = 0._rp
            grad_uh2_x = 0._rp
            select case(dim)
            case(3)
               do jj=1, nbdof
                  call matvecprod_3x3( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll))
                  grad_uh1_x = grad_uh1_x + v1(jj) * w
                  grad_uh2_x = grad_uh2_x + v2(jj) * w
               end do

            case(2)
               do jj=1, nbdof
                  call matvecprod_3x2( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh1_x = grad_uh1_x + v1(jj) * w
                  grad_uh2_x = grad_uh2_x + v2(jj) * w

               end do

            case(1)
               do jj=1, nbdof
                  call matvecprod_3x1( w, &
                       & g%Cy(:,:,ll), val(:,jj,ll) )
                  grad_uh1_x = grad_uh1_x + v1(jj) * w
                  grad_uh2_x = grad_uh2_x + v2(jj) * w

               end do

            end select
            grad_uh1_x = grad_uh1_x / g%Jy(ll)
            grad_uh2_x = grad_uh2_x / g%Jy(ll)

            ! update integral
            e_x = e(g%Ty(:,ll), grad_uh1_x, grad_uh2_x)
            int2 = int2 + e_x * wgt(ll) * g%Jy(ll)

         end do
      end do

      int = int + int2

    end subroutine cell_loop

  end function integ_scal_fe_grad_2




  function e_l2(x, u1, u2) result(r)  
    real(RP) :: r
    real(RP), dimension(3), intent(in) :: x
    real(RP),               intent(in) :: u1, u2
    r = (u1 - u2)**2
  end function e_l2

  function e_l2_vect(x, p1, p2) result(r)  
    real(RP) :: r
    real(RP), dimension(3), intent(in) :: x, p1, p2
    r = sum( (p1-p2)**2 ) 
  end function e_l2_vect


  function l2_dist_scalfe(u, uh, X_h, qdm) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: uh
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor)                   :: u

    if (choral_verb>1) write(*,*) &
         & 'integral        : L2_dist_scalFE'
    dist = integ(e_l2, u, uh, x_h, qdm)
    dist = sqrt(dist)

  end function l2_dist_scalfe


  function l2_dist_grad_scalfe(phi, uh, X_h, qdm) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: uh
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor3)                  :: phi

    if (choral_verb>1) write(*,*) &
         & 'integral        : L2_dist_grad_scalFE'
    dist = integ(e_l2_vect, phi, uh, x_h, qdm)
    dist = sqrt(dist)

  end function l2_dist_grad_scalfe

  function l2_dist_grad_proj(phi, uh, X_h, qdm) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: uh
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor3)                  :: phi

    if (choral_verb>1) write(*,*) &
         & 'integral        : L2_dist_grad_proj'
    dist = integ_scal_fe_grad_proj(e_l2_vect, phi, uh, x_h, qdm)
    dist = sqrt(dist)

  end function l2_dist_grad_proj



  function l2_dist_vect(phi, phi_h, X_h, qdm) result(dist)
    real(RP)                           :: dist
    procedure(r3tor3)                  :: phi
    real(RP), dimension(:), intent(in) :: phi_h
    type(fespace)         , intent(in) :: X_h
    type(quadmesh)        , intent(in) :: qdm

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

    if (choral_verb>1) write(*,*) &
         & "integral        : L2_dist_vect"

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

    if(.NOT.valid(qdm)) call quit(  &
         &"integral: L2_dist_vect:&
         & quad mesh 'qdmh' not valid" )

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"integral: L2_dist_vect:&
         & 'X_h' and 'qdmh' associated to different meshes" )

    if(size(phi_h,1) /= x_h%nbDof) call quit(  &
         &"integral: L2_dist_vect:&
         & fe function 'phi_h' has not the expected size" )

    m => x_h%mesh
    if (m%nbItf<=0) call quit( "integral: L2_dist_vect:&
         & mesh interfaces not defined")

    dist = 0.0_rp

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

       do qt=1, quad_tot_nb
          if (qdm%quad_count(qt)==0) cycle
          if (quad_geo(qt) /= fe_geo(ft) ) cycle

          dim   = quad_dim(qt)        ! dimension of K_ref
          nbdof = fe_nbdof(ft)        ! fe number of DOF
          nn    = quad_nbnodes(qt)    ! quad number of nodes

          call cell_loop()

       end do
    end do

    dist = sqrt( dist )

  contains

    subroutine cell_loop()

      real(RP), dimension(dim, nbDof, nn) :: base
      real(RP), dimension(3  , nbDof, nn) :: DT_base    
      real(RP), dimension(nbDof)          :: v
      integer , dimension(nbDof)          :: p       
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(     nn)        :: wgt
      real(RP), dimension(dim, nn)        :: y
      integer , dimension(CELL_MAX_NBITF) :: eps 
      real(RP), dimension(3)              :: phi_x, phi_h_x

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd, nb_itf
      real(RP)    :: s

      !! quad nodes and weights
      !!
      wgt = quad_wgt(qt)%y
      y   = quad_coord(qt)%y

      do ll=1, nn
         call vect_fe(base(:,:,ll), nbdof, y(:,ll), dim, ft)
      end do

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

         jj = qdm%quadType(cl)
         if (jj/=qt) cycle

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

         !! p = dof indexes associated with cell cl
         call getrow(p, nbdof, x_h%clToDof, cl)

         !! v = local dof
         v = phi_h( p )

         !! Modification of the local dof
         !! to take into account the interface orientation
#if DBG>0
         !!
         !! checks the cell orientation
         if (.NOT.cell_orientation(m, cl)) call quit(  &
              &   "integral : L2_dist_vect:&
              & cell orientation error" )
#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)
            do jj=1, nbdof
               if (eps(jj)==1) cycle
               v(jj) = -v(jj)
            end do

         case(fe_rt1_1d, fe_rt1_1d_2)
            do jj=1, 2
               if (eps(jj)==1) cycle
               v(jj) = -v(jj)
            end do

         case(fe_rt1_2d_2)
            do jj=1, 3
               if (eps(jj)==1) cycle
               ll = 2*jj-1
               v(ll:ll+1) = -v(ll:ll+1)
            end do

         case default
            call quit( "integral: L2_dist_vect:&
                 & fe type not supported" )

         end select


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

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

         !! DT_base(:,ii) \in R^3 =
         !!   DT(:,:,ll).base(:,ii,ll)
         !!   with DT(:,:,ll) = DT at node y(:,ll)
         !!
         select case(dim)

         case(3)
            do ll=1, nn
               do ii=1, nbdof
                  call matvecprod_3x3( dt_base(:,ii,ll), &
                       &   g%DTy(:,:,ll), base(:,ii,ll) )
               end do
            end do

         case(2)
            do ll=1, nn
               do ii=1, nbdof
                  call matvecprod_3x2( dt_base(:,ii,ll), &
                       &   g%DTy(:,:,ll), base(:,ii,ll) )
               end do
            end do

         case(1)
            do ll=1, nn
               do ii=1, nbdof
                  call matvecprod_3x1( dt_base(:,ii,ll), &
                       &   g%DTy(:,:,ll), base(:,ii,ll) )
               end do
            end do

         end select


         !! Update integral
         !!
         do ll=1, nn
            ! x = Ty(:,ll)
            ! phi_h( x ) 
            phi_h_x = dt_base(:,1,ll)*v(1)
            do jj=2, nbdof
               phi_h_x = phi_h_x + dt_base(:,jj,ll)*v(jj)
            end do
            s = 1.0_rp / g%Jy(ll)
            phi_h_x = phi_h_x * s


            !! phi_x = phi(x)
            phi_x = phi(g%Ty(:,ll)) 
            phi_x = (phi_h_x - phi_x)**2  
            dist = dist + sum(phi_x) * wgt(ll) * g%Jy(ll)

         end do

      end do

    end subroutine cell_loop

  end function l2_dist_vect


end module integral