elasticity.F90

Go to the documentation of this file.

module elasticity

  use choral_constants
  use choral_variables
  use real_type
  use basic_tools
  use abstract_interfaces, only: r3tor
  use algebra_lin
  use cell_mod
  use graph_mod
  use mesh_mod, only: mesh, valid
  use fe_mod
  use fespace_mod
  use quad_mod
  use quadmesh_mod
  use geotsf_mod
  use csr_mod
  use csr_fromgraph
  use fespacexk_mod
  use diffusion

  implicit none
  private

  public :: elasticity_matrix_pattern    
  public :: elasticity_massmat          !! tested
  public :: elasticity_stiffmat         !! tested
  public :: l2_product                  !! tested
  public :: elasticity_neumann_rhs
  public :: elasticity_dirichlet

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

  interface l2_product
     module procedure fespacexk_l2_product
  end interface l2_product

contains


  subroutine elasticity_matrix_pattern(g, Y, qdm)
    type(graph)    , intent(inout):: g
    type(fespacexk), intent(in)   :: Y
    type(quadmesh) , intent(in)   :: qdm

    logical, dimension(:), allocatable :: b
    type(fespace), pointer :: X_h
    type(graph) :: g0
    integer     :: ii

    if (choral_verb>2) write(*,*) &
         & 'elasticity      : matrix_pattern'

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: elasticity_matrix_pattern:&
         & feSpacexk 'Y' not valid" )
    x_h => y%X_h
    allocate( b(x_h%mesh%nbCl) )

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"elasticity: elasticity_matrix_pattern:&
         & 'Y' and 'qdmh' associated to 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, y%ClToDof2, b)
    deallocate(b)

    call prod_tab(g, g0, g0)

  end subroutine elasticity_matrix_pattern



  subroutine elasticity_massmat(mass, a, Y, qdm, dofToDof)
    type(csr)      , intent(inout):: mass
    procedure(r3tor)              :: a
    type(fespacexk), intent(in)   :: Y
    type(quadmesh) , intent(in)   :: qdm
    type(graph)    , intent(in), optional :: dofToDof

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

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

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: elasticity_massMat:&
         & feSpacexk 'Y' not valid" )

    !! shortcuts
    !!
    x_h => y%X_h
    m   => x_h%mesh

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"elasticity: elasticity_massMat:&
         & 'Y' and 'qdmh' associated to different meshes" )

    !! Y = [X_h]^k, with k = dimension in space
    !!
    if(y%k /= m%dim) call quit(  &
         &"elasticity: elasticity_massMat:&
         & exponent 'Y%k' must be equal to mesh%dim" )

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

       if (valid(doftodof)<0) call quit(  &
            & "elasticity: elasticity_massMat:&
            & matrix pattern 'dofToDof' not valid " )
       
       mass = csr(doftodof)
    else
       call elasticity_matrix_pattern(gph, y, 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

          if(dim /= m%dim) call quit(  &
               &'elasticity: elasticity_massMat:&
               & integration restricted to cell&
               & with dimension mesh%dim' )

          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, j2, l2, component
      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
               do component=1, dim

                  j2 = ( p(jj) - 1) * dim + component
                  l2 = ( p(ll) - 1) * dim + component
                  call addentry_2(mass, j2, l2, mat(jj,ll))

               end do
            end do
         end do

      end do

    end subroutine cell_loop

  end subroutine elasticity_massmat



  subroutine elasticity_stiffmat(stiff, lambda, mu, &
       &                         Y, qdm, dofToDof)
    type(csr)      , intent(inout):: stiff
    procedure(r3tor)              :: lambda, mu
    type(fespacexk), intent(in)   :: Y
    type(quadmesh) , intent(in)   :: qdm
    type(graph)    , intent(in), optional :: dofToDof

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

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

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: elasticity_stiffMat:&
         & feSpacexk 'Y' not valid" )

    !! shortcuts
    !!
    x_h => y%X_h
    m   => x_h%mesh

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

    if(.NOT.associated( qdm%mesh, x_h%mesh)) call quit(  &
         &"elasticity: elasticity_stiffMat:&
         & 'Y' and 'qdmh' associated to different meshes" )

    !! Y = [X_h]^k, with k = dimension in space
    !!
    if(y%k /= m%dim) call quit(  &
         &"elasticity: elasticity_stiffMat:&
         & exponent 'Y%k' must be equal to mesh%dim" )

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

       if (valid(doftodof)<0) call quit(  &
            & "elasticity: elasticity_stiffMat:&
            & matrix pattern 'dofToDof' not valid " )
       
       stiff = csr(doftodof)
    else
       call elasticity_matrix_pattern(gph, y, qdm)
       stiff = 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

          if(dim /= m%dim) call quit(  &
               &'elasticity: elasticity_stiffMat:&
               & integration restricted to cell&
               & with dimension mesh%dim' )

          nbdof2 = nbdof*dim

          call cell_loop()

       end do
    end do

  contains

    subroutine cell_loop()

      real(RP), dimension(dim, dim, dim, nbDof):: eu, Aeu
      real(RP), dimension(dim, nbDof, nn)      :: grad_ref     
      real(RP), dimension(3  , nbDof)          :: grad    
      real(RP), dimension(nbDof2, nbDof2)      :: mat 
      real(RP), dimension(dim, dim)            :: prod    
      integer , dimension(nbDof2)              :: p       
      integer , dimension(CELL_MAX_NBNODES)    :: p_cl    
      real(RP), dimension(3, CELL_MAX_NBNODES) :: X 
      real(RP), dimension(nn)                  :: sqrt_wgt

      type(geotsf):: g
      integer     :: cl, ii, jj, ll, cc, i2, c2, cl_nd, i, j
      real(RP)    :: inv_J, lbda_Ty, mu_Tyx2, val

      !! 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(  &
            & "elasticity: elasticity_stiffMat:&
            & quadrature rule with negative weights not allowed" )
      !!
      sqrt_wgt = sqrt(sqrt_wgt)

      !! gradient of the basis functions
      !! at the quadrature nodes
      !!
      do ll=1, nn
         call scalgrad_fe(grad_ref(:,:,ll), nbdof, &
              &           quad_coord(qt)%y(:,ll), dim, ft)
         grad_ref(:,:,ll) = grad_ref(:,:,ll) * sqrt_wgt(ll)
      end do

      !! geometricat transformation K_ref --> K
      !!
      g = geotsf( dim, nn, quad_coord(qt)%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

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

            !! Lamé coefficients
            lbda_ty = lambda(g%Ty(:,ll))
            mu_tyx2 = 2._rp * mu(g%Ty(:,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

            !! With the assumption Omega \subset \R^d
            !! we have grda(:,ii) \in \R^d
            !!
            !! We now compute the symmetrised gradient
            !!
            select case(dim)
            case(3)
               do ii=1, nbdof
                  eu(:,:,1,ii) = e3(grad(1:3,ii), 1)                     
                  eu(:,:,2,ii) = e3(grad(1:3,ii), 2)                     
                  eu(:,:,3,ii) = e3(grad(1:3,ii), 3)                     
               end do

            case(2)
               do ii=1, nbdof
                  eu(:,:,1,ii) = e2(grad(1:2,ii), 1)                     
                  eu(:,:,2,ii) = e2(grad(1:2,ii), 2)                     
               end do

            case(1)
               do ii=1, nbdof
                  eu(1,1,1,ii) = grad(1,ii)
               end do
            end select

            !! Hooke tensor
            !!
            aeu = mu_tyx2 * eu
            do ii=1, nbdof
               do cc=1, dim
                  val = lbda_ty * grad(cc,ii) 

                  do jj=1, dim
                     aeu(jj,jj,cc,ii) = aeu(jj,jj,cc,ii) + val
                  end do

               end do
            end do

            !! elementary stiffness matrix
            !!
            do ii=1, nbdof
               do cc=1, dim
                  i = (ii-1)*dim + cc

                  do i2=1, nbdof
                     do c2=1, dim
                        j = (i2-1)*dim + c2

                        prod = aeu(:,:,cc,ii) * eu(:,:,c2,i2)
                        val  = sum( prod ) * inv_j

                        mat(i,j) = mat(i,j) + val
                     end do
                  end do
               end do
            end do

         end do

         !! dof associated with cell cl
         call getrow(p, nbdof2, y%clToDof2, cl)

         !! add elementary stiffness matrix to global stiffness matrix
         do i=1, nbdof2
            do j=1, nbdof2
               call addentry_2(stiff, p(i), p(j), mat(i,j))
            end do
         end do

      end do

    end subroutine cell_loop

  end subroutine elasticity_stiffmat


  function e2(v, component)
    real(RP), dimension(2,2)            :: e2
    real(RP), dimension(2) , intent(in) :: v
    integer                , intent(in) :: component

    real(RP) :: a

    select case(component)
    case(1)
       a = v(2) * 0.5_rp

       e2(1,1) = v(1)
       e2(1,2) = a
       e2(2,1) = a
       e2(2,2) = 0.0_rp       

    case(2)
       a = v(1) * 0.5_rp

       e2(1,1) = 0.0_rp       
       e2(1,2) = a
       e2(2,1) = a
       e2(2,2) = v(2)

#if DBG>0
    case default
    call quit("elasticity: e2: 1 <= component <=2" )
#endif
    end select

  end function e2

  function e3(v, component)
    real(RP), dimension(3,3)            :: e3
    real(RP), dimension(3) , intent(in) :: v
    integer                , intent(in) :: component

    real(RP) :: a, b

    select case(component)
    case(1)
       a = v(2) * 0.5_rp
       b = v(3) * 0.5_rp

       e3(1  ,1) = v(1)
       e3(2  ,1) = a
       e3(3  ,1) = b
       e3(1  ,2) = a
       e3(2:3,2) = 0.0_rp 
       e3(1  ,3) = b
       e3(2:3,3) = 0.0_rp 

    case(2)
       a = v(1) * 0.5_rp
       b = v(3) * 0.5_rp

       e3(1,1)  = 0.0_rp 
       e3(2,1)  = a
       e3(3,1)  = 0.0_rp
       e3(1,2)  = a
       e3(2,2)  = v(2)
       e3(3,2)  = b
       e3(1,3)  = 0.0_rp
       e3(2,3)  = b
       e3(3,3)  = 0.0_rp

    case(3)
       a = v(1) * 0.5_rp
       b = v(2) * 0.5_rp

       e3(1:2,1) = 0.0_rp 
       e3(3  ,1) = a
       e3(1:2,2) = 0.0_rp 
       e3(3  ,2) = b
       e3(1  ,3) = a
       e3(2  ,3) = b
       e3(3  ,3) = v(3)  

#if DBG>0
    case default
    call quit("elasticity: e3: 1 <= component <=3" )
#endif
    end select

  end function e3



  subroutine fespacexk_l2_product(FV, Y, qdm, f_1, f_2, f_3) 
    real(RP), dimension(:), allocatable :: FV
    type(fespacexk)       , intent(in)  :: Y
    type(quadmesh)        , intent(in)  :: qdm
    procedure(r3tor)                    :: f_1, f_2
    procedure(r3tor)      , optional    :: f_3

    type(fespace), pointer :: X_h
    integer :: i, dim, jj, ll
    real(RP), dimension(:), allocatable :: FV_i

    if (choral_verb>1) write(*,*) &
         & "elasticity      : feSpacexk_L2_product"

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: feSpacexk_L2_product: &
         & fe space 'Y' not valid" )

    x_h => y%X_h
    dim = x_h%mesh%dim
    if (y%k /= dim) call quit(&
         &"elasticity: feSpacexk_L2_product: &
         & fespacsexk exponent 'Y%k' incorrect" )
    if ( dim==1 ) call quit(&
         &"elasticity: feSpacexk_L2_product: &
         & dim 1 case, use diffusion::L2_product" )
    if ( (dim==3).AND.(.NOT.present(f_3)) ) call quit(&
         &"elasticity: feSpacexk_L2_product: &
         & argument 'f_3' missing" )

    call allocmem(fv, y%nbDof2)
    fv = 0._rp

    !! first component
    i = 1
    call l2_product(fv_i, f_1, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       fv(ll)  = fv_i(jj)
    end do

    !! second component
    i = 2
    call l2_product(fv_i, f_2, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       fv(ll)  = fv_i(jj)
    end do   
    if (dim==2) return

    !! third component
    i = 3
    call l2_product(fv_i, f_3, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       fv(ll)  = fv_i(jj)
    end do
    
  end subroutine fespacexk_l2_product


  subroutine elasticity_neumann_rhs(rhs, Y, quad_type, g_1, g_2, g_3,  f) 
    real(RP), dimension(:), allocatable :: rhs
    type(fespacexk)       , intent(in)  :: Y
    integer               , intent(in)  :: quad_type
    procedure(r3tor)                    :: g_1, g_2
    procedure(r3tor)      , optional    :: g_3, f

    type(quadmesh)         :: qdm
    type(fespace), pointer :: X_h
    integer :: i, dim, jj, ll
    real(RP), dimension(:), allocatable :: rhs_i

    if (choral_verb>1) write(*,*) &
         & "elasticity      : elasticity_Neumann_rhs"

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: elasticity_Neumann_rhs: &
         & fe space 'Y' not valid" )

    x_h => y%X_h
    dim = x_h%mesh%dim

    if (y%k /= dim) call quit(&
         &"elasticity: elasticity_Neumann_rhs: &
         & fespacsexk exponent 'Y%k' incorrect" )
    if ( dim==1 ) call quit(&
         &"elasticity: elasticity_Neumann_rhs: &
         & dim 1 case, use diffusion::diffusion_Neumann_rhs" )
    if ( (dim==3).AND.(.NOT.present(g_3)) ) call quit(&
         &"elasticity: elasticity_Neumann_rhs: &
         & argument 'g_3' missing" )
    if ( quad_dim(quad_type) /= dim-1 ) call quit(&
         &"elasticity: elasticity_Neumann_rhs: &
         & unproper quadrature method" )

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

    call allocmem(rhs, y%nbDof2)
    rhs = 0._rp


    !! first component
    i = 1
    call l2_product(rhs_i, g_1, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       rhs(ll)  = rhs_i(jj)
    end do

    !! second component
    i = 2
    call l2_product(rhs_i, g_2, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       rhs(ll)  = rhs_i(jj)
    end do   
    if (dim==2) return

    !! third component
    i = 3
    call l2_product(rhs_i, g_3, x_h, qdm)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       rhs(ll)  = rhs_i(jj)
    end do

  end subroutine elasticity_neumann_rhs



  subroutine elasticity_dirichlet(K, rhs, Y, g1, g2, g3,  f) 
    type(csr)             , intent(inout) :: K
    real(RP), dimension(:), intent(inout) :: rhs
    type(fespacexk)       , intent(in)    :: Y
    procedure(r3tor)                      :: g1, g2
    procedure(r3tor)      , optional      :: g3, f

    logical , dimension(:), allocatable   :: flag
    real(RP), dimension(:), allocatable   :: g1_h, g2_h, g3_h

    type(fespace), pointer :: X_h
    integer :: dim, ii, ll, j1, j2

    if (choral_verb>1) write(*,*) &
         & "elasticity      : elasticity_Dirichlet"

    if(.NOT.valid(y)) call quit(  &
         &"elasticity: elasticity_Dirichlet_rhs: &
         & fe space 'Y' not valid" )

    x_h => y%X_h
    dim = x_h%mesh%dim

    if (y%k /= dim) call quit(&
         &"elasticity: elasticity_Dirichlet_rhs: &
         & fespacsexk exponent 'Y%k' incorrect" )
    if ( dim==1 ) call quit(&
         &"elasticity: elasticity_Dirichlet_rhs: &
         & dim 1 case, use diffusion::diffusion_Dirichlet_rhs" )
    if ( (dim==3).AND.(.NOT.present(g3)) ) call quit(&
         &"elasticity: elasticity_Dirichlet_rhs: &
         & argument 'g3' missing" )

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

    if(k%nl/=y%nbDof2) call quit(  &
         &"elasticity: elasticity_Dirichlet: &
         & csr matrix 'K' has a wrong size" )

    if ( size(rhs,1) /= y%nbDof2 ) call quit(  &
         &"elasticity: elasticity_Dirichlet: &
         & 'rhs' has a wrong size" )

    !! Flag the finite element nodes in {\Gamma + f >=0}
    !!
    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(g1_h, g1, x_h)
    call interp_scal_func(g2_h, g2, x_h)
    if (dim==3)  call interp_scal_func(g3_h, g3, x_h)

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

       if ( .NOT. flag(ii) ) cycle

       !! modifies rhs(ll)
       ll = (ii-1)*dim 
       rhs(ll+1) = g1_h(ii) 
       rhs(ll+2) = g2_h(ii) 
       if (dim==3) rhs(ll+3) = g3_h(ii) 

       !! modifies K(ll,:)
       j1 = k%row(ll+1)
       if (dim==2)   j2 = k%row(ll+3)-1
       if (dim==3)   j2 = k%row(ll+4)-1
       k%a(j1:j2) = 0.0_rp

       j1 = k%diag(ll+1)
       k%a(j1) = 1.0_rp
       j1 = k%diag(ll+2)
       k%a(j1) = 1.0_rp
       if (dim==3) then
          j1 = k%diag(ll+3)
          k%a(j1) = 1.0_rp
       end if

    end do

  end subroutine elasticity_dirichlet


end module elasticity