diffusion.F90

Go to the documentation of this file.

module diffusion

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use abstract_interfaces, only: r3tor, r3xr3xr3tor
  use algebra_lin
  use cell_mod
  use graph_mod
  use mesh_mod
  use mesh_interfaces
  use fe_mod
  use fespace_mod
  use quad_mod
  use quadmesh_mod
  use geotsf_mod
  use csr_mod
  use csr_fromgraph

  implicit none
  private

  public :: diffusion_matrix_pattern
  public :: diffusion_massmat            ! tested
  public :: diffusion_stiffmat           ! tested
  public :: diffusion_massmat_vect       ! tested
  public :: diffusion_mixed_divmat       ! tested
  public :: l2_product                   ! tested 
  public :: diffusion_neumann_rhs        ! tested 
  public :: diffusion_dirichlet

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

  interface l2_product
     module procedure fespace_l2_product
  end interface l2_product


contains


  subroutine diffusion_matrix_pattern(g, X_h, qdm)
    type(graph)   , intent(inout):: g
    type(fespace) , intent(in)   :: X_h
    type(quadmesh), intent(in)   :: qdm

    logical, dimension(X_h%mesh%nbCl) :: b
    type(graph) :: g0
    integer     :: ii

    if (choral_verb>2) write(*,*) &
         & 'diffusion       : diffusion_matrix_pattern'

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

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"diffusion: diffusion_matrix_pattern:&
         & 'X_h' and 'qdm' associated with different meshes" )

    do ii=1, x_h%mesh%nbCl
       if (   (qdm%quadType(ii)==quad_none) .OR.&
            & (x_h%feType(ii)==fe_none)  ) then      
          b(ii) = .false.
       else
          b(ii) = .true.
       end if
    end do

    call graph_extract(g0, x_h%ClToDof, b)
    call prod_tab(g, g0, g0)

  end subroutine diffusion_matrix_pattern


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

  subroutine diffusion_massmat(mass, a, X_h, qdm, dofToDof)
    type(csr)     , intent(inout):: mass
    procedure(r3tor)             :: a
    type(fespace) , intent(in)   :: X_h
    type(quadmesh), intent(in)   :: qdm
    type(graph)   , intent(in), optional :: dofToDof

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

    if (choral_verb>1) write(*,*) &
         & "diffusion       : massMat"
    if (choral_verb>2) write(*,*)&
         & "  pattern provided               = ",&
         & present(doftodof)

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

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"diffusion: diffusion_massMat:&
         & 'X_h' and 'qdm' associated with different meshes" )


    !! connectivity graph 
    !!   DOF  -->  DOF
    if (present(doftodof)) then

       if ( valid(doftodof) < 0 ) call quit(  &
            & "diffusion: diffusion_massMat:&
            & matrix pattern 'dofToDof' not valid " )
       
       mass = csr(doftodof)
    else
       call diffusion_matrix_pattern(gph, x_h, qdm)
       mass = csr(gph)
       call clear(gph)
    end if
    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, nbDof)    :: mat     
      real(RP), dimension(nbDof, nbDof, nn):: mat0     
      real(RP), dimension(nbDof)           :: u       
      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

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

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

      !! u = basis functions evaluated at quad node ll
      !!
      do ll=1, nn             
         call scal_fe(u, nbdof, y(:,ll), dim, ft)
         do ii=1, nbdof
            do jj=1, nbdof
               mat0(ii,jj,ll) = u(ii)*u(jj)*wgt(ll)
            end do
         end do
      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(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_j(g)

         ! elementary mass matrix on cell cl
         mat = 0._rp

         ! loop on quad nodes
         do ll=1, nn
            s = a(g%Ty(:,ll))*g%Jy(ll)             

            do ii=1, nbdof
               do jj=1, ii
                  mat(ii,jj) = mat(ii,jj) + &
                       & s* mat0(ii,jj,ll)
               end do
            end do
         end do

         ! completion by symmetry
         do ii=2, nbdof
            do jj=1, ii-1
               mat(jj,ii) = mat(ii,jj) 
            end do
         end do

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

         ! add elementary mass matrix to global mass matrix
         do jj=1, nbdof
            do ll=1, nbdof
               call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
            end do
         end do

      end do

    end subroutine cell_loop

  end subroutine diffusion_massmat


  subroutine diffusion_stiffmat(stiff, b, X_h, qdm, dofToDof)
    type(csr)     , intent(inout) :: stiff
    procedure(r3xr3xr3tor)        :: b
    type(fespace) , intent(in)    :: X_h
    type(quadmesh), intent(in)    :: qdm
    type(graph)   , intent(in), optional :: dofToDof

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

    if (choral_verb>1) write(*,*) &
         & "diffusion       : stiffMat"
    if (choral_verb>2) write(*,*) &
         & "  pattern provided               = ",&
         & present(doftodof)

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

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"diffusion: diffusion_stiffMat:&
         & 'X_h' and 'qdm' associated with different meshes" )

    !! connectivity graph 
    !!   DOF  -->  DOF
    if (present(doftodof)) then

       if ( valid(doftodof) < 0 ) call quit(  &
            & "diffusion: diffusion_stiffMat:&
            & matrix pattern 'dofToDof' not valid " )
       
       stiff = csr(doftodof)
    else
       call diffusion_matrix_pattern(gph, x_h, qdm)
       stiff = csr(gph)
       call clear(gph)
    end if

    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)     :: grad_ref     
      real(RP), dimension(3  , nbDof)         :: grad    
      real(RP), dimension(nbDof, nbDof)       :: mat     
      integer , dimension(nbDof)              :: p       
      integer , dimension(CELL_MAX_NBNODES)   :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES):: X 
      real(RP), dimension(nn)                 :: sqrt_wgt
      real(RP), dimension(dim, nn)            :: y

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cl_nd
      real(RP)    :: s, inv_J
      real(RP), dimension(3) :: Ty

      !! square roots of the quad weights
      !!
      sqrt_wgt = quad_wgt(qt)%y
      !!
      !! negative weights not allowed
      if ( minval(sqrt_wgt) < 0.0_rp ) call quit(  &
            & "diffusion: diffusion_stiffMat:&
            & quadrature rule with negative weights not allowed" )
      !!
      sqrt_wgt = sqrt(sqrt_wgt)
      y        = quad_coord(qt)%y


      !! gradient of the basis functions
      !! at the quadrature nodes
      !!
      do ll=1, nn
         call scalgrad_fe(grad_ref(:,:,ll), nbdof, &
              &           y(:,ll), dim, ft)
         grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
      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(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_jc(g)

         !! elementary stiffness matrix on cell cl
         !!
         mat = 0._rp

         !! loop on quad nodes
         do ll=1, nn

            !! node coordinates
            ty = g%Ty(:,ll)

            !! Jacobian inverse
            inv_j = 1.0_rp / g%Jy(ll)

            !! grad(:,ii) = Cy(:,:,ll) * grad_ref(:,ii,ll)
            !!
            select case(dim)
            case(3)
               do ii=1, nbdof
                  call matvecprod_3x3( grad(:,ii), &
                       & g%Cy(:,:,ll), grad_ref(:,ii,ll))
               end do
               
            case(2)
               do ii=1, nbdof
                  call matvecprod_3x2( grad(:,ii), &
                       & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
               end do
               
            case(1)
               do ii=1, nbdof
                  call matvecprod_3x1( grad(:,ii), &
                       & g%Cy(:,:,ll), grad_ref(:,ii,ll) )
               end do
               
            end select

            do ii=1, nbdof
               do jj=1, ii
                  s =  b( ty, grad(:,ii), grad(:,jj) )
                  mat(ii,jj) = mat(ii,jj) + s * inv_j                   
               end do
            end do

         end do
         !!
         !! completion by symmetry
         do ii=2, nbdof
            do jj=1, ii-1
               mat(jj,ii) = mat(ii,jj) 
            end do
         end do

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

         ! add elementary stiffness matrix to global stiffness matrix
         do jj=1, nbdof
            do ll=1, nbdof
               call addentry_2(stiff, p(jj), p(ll), mat(jj,ll))
            end do
         end do
      end do

    end subroutine cell_loop

  end subroutine diffusion_stiffmat



  subroutine diffusion_neumann_rhs(rhs, g, X_h, quad_type, f) 
    real(RP), dimension(:), allocatable :: rhs
    procedure(r3tor)                    :: g
    type(fespace)         , intent(in)  :: X_h
    integer               , intent(in)  :: quad_type
    procedure(r3tor)      , optional    :: f

    type(quadmesh)      :: qdm

    if (choral_verb>1) write(*,*) &
         & "diffusion       : diffusion_Neumann_rhs"

    if(.NOT.valid(x_h)) call quit(  &
         &"diffusion: diffusion_Neumann_rhs: &
         & fe space 'X_h' not valid" )
    if ( quad_dim(quad_type) /= x_h%mesh%dim-1 ) call quit(&
         &"diffusion: diffusion_Neumann_rhs: &
         & unproper quadrature method" )

    qdm = quadmesh(x_h%mesh)
    call set(qdm, quad_type, f)
    call assemble(qdm)
    call l2_product(rhs, g, x_h, qdm)

  end subroutine diffusion_neumann_rhs



  subroutine diffusion_dirichlet(K, rhs, g, X_h, rho, f) 
    type(csr)             , intent(inout) :: K
    real(RP), dimension(:), intent(inout) :: rhs
    procedure(r3tor)                      :: g
    type(fespace)         , intent(in)    :: X_h
    real(RP)              , intent(in)    :: rho
    procedure(r3tor)      , optional      :: f

    logical , dimension(:), allocatable   :: flag
    real(RP), dimension(:), allocatable   :: g_h

    integer :: ii, jj

    if (choral_verb>1) write(*,*) &
         & "diffusion       : diffusion_Dirichlet"

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

    if(.NOT.valid(k)) call quit(  &
         &"diffusion: diffusion_Dirichlet: &
         & csr matrix 'K' not valid" )

    if(k%nl/=x_h%nbDof) call quit(  &
         &"diffusion: diffusion_Dirichlet: &
         & csr matrix 'K' has a wrong size" )

    if ( size(rhs,1) /= x_h%nbDof ) call quit(  &
         &"diffusion: diffusion_Dirichlet: &
         & 'rhs' has a wrong size" )

    call flag_x_h_dof(flag, x_h, x_h%mesh%dim-1, f)

    !! interpolate the boundary source 'g'
    !! on all finite element nodes
    call interp_scal_func(g_h, g, x_h)

    !! penalise rhs(ii) and K(ii,ii)
    !! on index ii such that
    !! such that flag(ii) = .TRUE.
    call getdiag(k)
    do ii=1, x_h%nbDof

       if ( .NOT. flag(ii) ) cycle

       !! penalise rhs(ii)
       rhs(ii) = rhs(ii) + g_h(ii) * rho

       !! penalise K(ii,ii)
       jj = k%diag(ii)
       k%a(jj) = k%a(jj) + rho

    end do

  end subroutine diffusion_dirichlet



  subroutine fespace_l2_product(FV, f, X_h, qdm) 
    real(RP), dimension(:), allocatable :: FV
    type(fespace)         , intent(in)  :: X_h
    type(quadmesh)        , intent(in)  :: qdm
    procedure(r3tor)                    :: f

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

    if (choral_verb>1) write(*,*) &
         & "diffusion       : feSpace_L2_product"

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

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"diffusion: feSpace_L2_product:&
         & 'X_h' and 'qdm' associated with different meshes" )

    call allocmem(fv, x_h%nbDof)
    fv = 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     
      real(RP), dimension(nbDof)         :: xx, u    
      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

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

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

      !! u = basis functions evalated at quad node ll
      !!
      do ll=1, nn             
         call scal_fe(u, nbdof,  y(:,ll), dim, ft)
         val(:,ll) = u * wgt(ll)
      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(1:3,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)

         xx = 0._rp
         ! loop on quadrature nodes
         do ll=1, nn
            fx = f(g%Ty(:,ll)) 
            fx = fx * g%Jy(ll) 
            xx   = xx + val(:,ll)*fx
         end do

         fv(p) = fv(p) + xx

      end do

    end subroutine cell_loop

  end subroutine fespace_l2_product

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!
  !!    VECTOR FINITE ELEMENTS
  !!
  !!
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


  subroutine diffusion_massmat_vect(mass, b, X_h, qdm, dofToDof)
    type(csr)     , intent(inout) :: mass
    procedure(r3xr3xr3tor)        :: b
    type(fespace)  , intent(in)   :: X_h
    type(quadmesh), intent(in)    :: qdm
    type(graph)   , intent(in), optional :: dofToDof

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

    if (choral_verb>1) write(*,*) &
         & "diffusion       : massMat_vect"
    if (choral_verb>2) write(*,*) &
         & "  pattern provided               = ",&
         & present(doftodof)

    m => x_h%mesh

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

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"diffusion: diffusion_massMat_vect:&
         & 'X_h' and 'qdm' associated with different meshes" )

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

    !! connectivity graph 
    !!   DOF  -->  DOF
    if (present(doftodof)) then

       if ( valid(doftodof) < 0 ) call quit(  &
            & "diffusion: diffusion_massMat_vect:&
            & matrix pattern 'dofToDof' not valid " )
       
       mass = csr(doftodof)
    else
       call diffusion_matrix_pattern(gph, x_h, qdm)
       mass = csr(gph)
       call clear(gph)
    end if

    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) :: phi
      real(RP), dimension(3  , nbDof, nn) :: DT_phi    
      real(RP), dimension(nbDof, nbDof)   :: mat     
      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 

      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

      !! negative weights not allowed
      if ( minval(wgt) < 0.0_rp ) call quit(  &
            & "diffusion: diffusion_massMat_vect:&
            & quadrature rule with negative weights not allowed" )

      do ll=1, nn
         call vect_fe(phi(:,:,ll), nbdof, y(:,ll), dim, ft)
         phi(:,:,ll) = phi(:,:,ll) * sqrt(wgt(ll))
      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(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_phi(:,ii) \in \R^3 =
         !!   DT(:,:,ll).phi(:,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_phi(:,ii,ll), &
                       & g%DTy(:,:,ll), phi(:,ii,ll) )
               end do
            end do

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

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

         end select

         !! elementary mass matrix on cell cl
         !!
         mat = 0._rp
         !!
         !! Jacobian inverse
         g%Jy = 1._rp / g%Jy
         !!
         !! loop on quad nodes
         do ll=1, nn

            do ii=1, nbdof
               do jj=1, ii
                  s =  b( g%Ty(:,ll), dt_phi(:,ii,ll), &
                       &               dt_phi(:,jj,ll)  )
                  mat(ii,jj) = mat(ii,jj) + s * g%Jy(ll)                
               end do
            end do

         end do
         !!
         !! completion by symmetry
         do ii=2, nbdof
            do jj=1, ii-1
               mat(jj,ii) = mat(ii,jj) 
            end do
         end do

         !! Modification of the elementary mass matrix
         !! to take into account the interface orientation
#if DBG>0
         !!
         !! checks the cell orientation
         if (.NOT.cell_orientation(m, cl)) call quit(  &
              & "diffusion: diffusion_massMat_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 ii = 1, nbdof
               if (eps(ii)==1) cycle
               mat(ii, :) = -mat(ii, :)
               mat( :,ii) = -mat( :,ii)
            end do

         case(fe_rt1_1d)
            do ii = 1, 2
               if (eps(ii)==1) cycle
               mat(ii, :) = -mat(ii, :)
               mat( :,ii) = -mat( :,ii)
            end do

         case(fe_rt1_2d_2)
            do ii = 1, 3
               if (eps(ii)==1) cycle

               mat( 2*ii-1, :) = -mat( 2*ii-1, :)
               mat( :, 2*ii-1) = -mat( :, 2*ii-1)
               mat( 2*ii  , :) = -mat( 2*ii  , :)
               mat( :, 2*ii  ) = -mat( :, 2*ii  )         
            end do

         case default
            call quit("diffusion: diffusion_massMat_vect:&
                 & fe type not supported" )
         end select

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

         ! add elementary mass matrix to global mass matrix
         do jj=1, nbdof
            do ll=1, nbdof
               call addentry_2(mass, p(jj), p(ll), mat(jj,ll))
            end do
         end do

      end do

    end subroutine cell_loop

  end subroutine diffusion_massmat_vect


  subroutine diffusion_mixed_divmat(divMat, X_s, X_v, qdm)
    type(csr)     , intent(inout) :: divMat
    type(fespace) , intent(in)    :: X_s, X_v
    type(quadmesh), intent(in)    :: qdm

    type(graph)         :: dof_sToDof_v

    type(mesh), pointer :: m
    integer  :: dim, nbDof_s, nbDof_v, nn, ft_s, ft_v, qt

    if (choral_verb>1) write(*,*) &
         & "diffusion       : mixed_divMat"

    m => x_s%mesh

    if(.NOT.valid(x_s)) call quit(  &
         &"diffusion: diffusion_mixed_divMat:&
         & scalar fe space 'X_s' not valid" )

    if(.NOT.valid(x_v)) call quit(  &
         &"diffusion: diffusion_mixed_divMat:&
         & vector fe space 'X_v' not valid" )

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

    if(.NOT.associated( qdm%mesh, x_s%mesh)) call quit(  &
         &"diffusion: diffusion_mixed_divMat:&
         & 'X_s' and 'qdm' associated with different meshes" )

    if(.NOT.associated( qdm%mesh, x_v%mesh)) call quit(  &
         &"diffusion: diffusion_mixed_divMat:&
         & 'X_v' and 'qdm' associated with different meshes" )

    !! connectivity graph 
    !!   scalar DOF  -->  vector DOF
    call prod_tab(dof_stodof_v, x_s%clToDof, x_v%clToDof)
    divmat = csr(dof_stodof_v)
    call clear(dof_stodof_v)

    do ft_s=1, fe_tot_nb
       if (x_s%fe_count(ft_s)==0) cycle

       do ft_v=1, fe_tot_nb
          if (x_v%fe_count(ft_v)==0) cycle

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

             dim     = quad_dim(qt)      ! dimension of K_ref
             nbdof_s = fe_nbdof(ft_s)    ! scalar fe number of DOF
             nbdof_v = fe_nbdof(ft_v)    ! vector fe number of DOF
             nn      = quad_nbnodes(qt)  ! quad number of nodes

             call cell_loop()

          end do
       end do
    end do

  contains

    subroutine cell_loop()

      real(RP), dimension(nbDof_v, nn)        :: div
      real(RP), dimension(nbDof_s, nn)        :: u
      real(RP), dimension(nbDof_s, nbDof_v)   :: mat_0, mat   
      integer , dimension(nbDof_s)            :: p_s
      integer , dimension(nbDof_v)            :: p_v
      real(RP), dimension(     nn)            :: wgt
      real(RP), dimension(dim, nn)            :: y
      integer , dimension(CELL_MAX_NBITF)     :: eps 

      integer     :: cl, ii, jj, ll, nb_itf
      real(RP)    :: s

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

      !! Scalar basis functions and
      !! divergence of the vector basis functions
      !! evaluated at the quad nodes
      !!
      do ll=1, nn
         call scal_fe(u(:,ll), nbdof_s, y(:,ll), dim, ft_s)
         call vectdiv_fe(div(:,ll), nbdof_v, y(:,ll), dim, ft_v)
      end do

      !! elementary matrix 
      !!
      mat_0 = 0._rp
      !!
      !! loop on quad nodes
      do ll=1, nn

         do ii=1, nbdof_s
            do jj=1, nbdof_v
               s = u(ii,ll) * div(jj,ll) * wgt(ll) 
               mat_0(ii,jj) = mat_0(ii,jj) + s             
            end do
         end do

      end do


      !! CELL LOOP
      !!
      do cl=1, m%nbCl

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

         ii = x_s%feType(cl)
         if (ii/=ft_s) cycle

         ii = x_v%feType(cl)
         if (ii/=ft_v) cycle

         !! elementary matrix on cell cl
         !!
         mat = mat_0
         !!
         !! Modification of the elementary matrix
         !! to take into account the interface orientation
#if DBG>0
         !!
         !! checks the cell orientation
         if (.NOT.cell_orientation(m, cl)) call quit(  &
              &"diffusion: diffusion_mixed_divMat:&
              & 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_v)
         case(fe_rt0_1d, fe_rt0_1d_dual, fe_rt0_2d)
            do ii = 1, nbdof_v
               if (eps(ii)==1) cycle
               
               mat(:, ii) = -mat(:, ii)
            end do

         case(fe_rt1_1d, fe_rt1_1d_dual,  &
              &fe_rt1_1d_2, fe_rt1_1d_2_dual)
            do ii = 1, 2
               if (eps(ii)==1) cycle
               
               mat(:, ii) = -mat(:, ii)
            end do

         case(fe_rt1_2d_2)
            do ii = 1, 3
               if (eps(ii)==1) cycle

               mat( :, 2*ii-1) = -mat( :, 2*ii-1)
               mat( :, 2*ii  ) = -mat( :, 2*ii  )         
            end do

         case default
            call quit("diffusion: diffusion_mixed_divMat:&
                 & fe type not supported" )
         end select

         ! dof associated with cell cl indexes
         call getrow(p_s, nbdof_s, x_s%clToDof, cl)
         call getrow(p_v, nbdof_v, x_v%clToDof, cl)

         ! add elementary matrix to global divMat matrix
         do ii=1, nbdof_s
            do jj=1, nbdof_v
               call addentry_2(divmat, p_s(ii), p_v(jj), mat(ii,jj))
            end do
         end do

      end do

    end subroutine cell_loop

  end subroutine diffusion_mixed_divmat


end module diffusion