feSpacexk_mod.f90

Go to the documentation of this file.

!
module fespacexk_mod

  use choral_variables
  use real_type
  use basic_tools
  use abstract_interfaces, only: r3tor, r3tor3
  use graph_mod
  use mesh_mod
  use fespace_mod
  use quadmesh_mod
  use integral

  implicit none
  private


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

  public :: clear
  public :: valid
  public :: print
  public :: interp_vect_func     !! tested
  public :: extract_component    !! tested
  public :: l2_dist
  public :: l2_dist_grad
  
  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type fespacexk

     type(fespace), pointer :: x_h=>null()

     integer :: k

     integer :: nbdof2=0

     integer , dimension(:), allocatable :: dof2todof 


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

   contains

     final :: fespacexk_clear

  end type fespacexk


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


  !! constructor
  interface fespacexk
     module procedure fespacexk_create
  end interface fespacexk

  interface valid
     module procedure fespacexk_valid
  end interface valid

  interface interp_vect_func
     module procedure fespacexk_interp_vect_func
  end interface interp_vect_func

  interface print
     module procedure fespacexk_print
  end interface print

  interface l2_dist
     module procedure l2_dist_2
  end interface l2_dist
  interface l2_dist_grad
     module procedure l2_dist_grad_2
  end interface l2_dist_grad

contains

  subroutine fespacexk_clear(Y)
    type(fespacexk), intent(inout) :: Y

    y%X_h       => null()
    y%nbDof2    =  0
    call freemem(y%dof2ToDof)
    call clear(  y%clToDof2)

  end subroutine fespacexk_clear

  function fespacexk_create(X_h, k) result(Y)
    type(fespacexk)             :: Y
    type(fespace)  , target     :: X_h
    integer        , intent(in) :: k

    integer :: cl, ii, jj, ll, wdt, cpt, nn
    integer, dimension(:), allocatable :: nnz, row, row2

    if (choral_verb>0) write(*,*) &
         & 'feSpacexk_mod   : create       k =', int(k,1)
    call clear(y)
    if (k<1) call quit( &
         & "feSpacexk_mod: feSpacexk_create:&
         & wrong argument 'k'" )
    if (.NOT.valid(x_h)) call quit( &
         &"feSpacexk_mod: feSpacexk_create:&
         & fe space 'X_h' not valid" )

    !! parameters
    y%X_h => x_h
    y%k      = k
    y%nbDof2 = x_h%nbDof * k

    !! correspondance dof of Y --> dof of X_h
    call allocmem( y%dof2ToDof, y%nbDof2)
    do ii=1, y%nbDof2
       y%dof2ToDof(ii) = (ii-1)/k + 1
    end do

    !! connectivity cell --> dof of Y
    !!
    call allocmem( nnz, y%X_h%clToDof%nl )
    do cl=1, x_h%clToDof%nl

       jj = x_h%clToDof%row(cl)
       ll = x_h%clToDof%row(cl+1)
       nnz(cl) = ll - jj

    end do
    nnz = nnz * k
    y%clToDof2 = graph(nnz, size(nnz,1), y%nbDof2)
    call freemem(nnz)
    !!
    !! fill in the graph
    !!
    wdt = x_h%clToDof%width
    call allocmem(row , wdt)
    call allocmem(row2, wdt*k)
    do cl=1, x_h%clToDof%nl
       
       call getrow(ll, row, wdt, x_h%clToDof, cl)
       cpt = 0
       do ii=1, ll  
          nn = (row(ii)-1) * k
          do jj=1, k
             cpt = cpt + 1             
             row2(cpt) = nn + jj
          end do          
       end do
       nn = ll * k
       call setrow(y%clToDof2, row2(1:nn), nn, cl)
       
    end do

    if (.NOT.valid(y)) call quit(&
         & 'feSpacexk_mod : feSpacexk_create:&
         & assembling not valid' )

  end function fespacexk_create
  

  function fespacexk_valid(Y) result(b)
    type(fespacexk), intent(in) :: y
    logical :: b

    b = .false.
    if(.NOT. valid(y%X_h)           ) return
    if( valid(y%clToDof2) < 0       ) return
    if(.NOT. allocated(y%dof2ToDof) ) return

    b =             y%nbdof2>0 
    b =  b .AND. (  y%k     >0  )
    b =  b .AND. ( (y%k * y%X_h%nbDof  ) == y%nbdof2 )
    b =  b .AND. ( size( y%dof2ToDof, 1) == y%nbDof2 )

  end function fespacexk_valid



  subroutine fespacexk_interp_vect_func(u_h, Y, u_1, u_2, u_3)
    real(RP), dimension(:), allocatable :: u_h
    type(fespacexk)       , intent(in)  :: Y
    procedure(r3tor)                    :: u_1, u_2
    procedure(r3tor)      , optional    :: u_3

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

    if (choral_verb>1) write(*,*) &
         & "feSpacexk_mod   : interp_vetc_func"

    if (.NOT.valid(y)) call quit(&
         &"feSpacexk_mod: feSpacexk_interp_vect_func:&
         & fespacsexk 'Y' not valid" )

    x_h => y%X_h
    dim = x_h%mesh%dim
    if (y%k /= dim) call quit(&
         &"feSpacexk_mod: feSpacexk_interp_vect_func:&
         & fespacsexk exponent 'Y%k' incorrect" )
    if ( dim==1 ) call quit(&
         &"feSpacexk_mod: feSpacexk_interp_vect_func:&
         & dim 1 case, use feSpace_mod::interp_scal_func" )
    if ( (dim==3).AND.(.NOT.present(u_3)) ) call quit(&
         &"feSpacexk_mod: feSpacexk_interp_vect_func:&
         & argument 'u3' missing" )

    call allocmem(u_h   , y%nbDof2)

    !! first component
    i = 1
    call interp_scal_func(u_ih, u_1, x_h)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       u_h(ll) = u_ih(jj)
    end do

    !! second component
    i = 2
    call interp_scal_func(u_ih, u_2, x_h)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       u_h(ll) = u_ih(jj)
    end do   
    if (dim==2) return

    !! third component
    i = 3
    call interp_scal_func(u_ih, u_3, x_h)
    do jj=1, x_h%nbDof
       ll      = (jj-1) * y%k + i
       u_h(ll) = u_ih(jj)
    end do

  end subroutine fespacexk_interp_vect_func

  subroutine fespacexk_print(Y)
    type(fespacexk), intent(in) :: Y

    write(*,*)"feSpacexk_mod   : print"
    write(*,*)"  Number of DOF                  =", y%nbDof2
    write(*,*)"  Exponent                       =", y%k
    write(*,*)"  Valid                          =", valid(y)

  end subroutine fespacexk_print


  subroutine extract_component(u_c, u, Y, c)
    real(RP), dimension(:), allocatable :: u_c
    real(RP), dimension(:), intent(in)  :: u
    type(fespacexk)       , intent(in)  :: Y
    integer               , intent(in)  :: c

    integer :: k, jj, ll

    if (choral_verb>1) write(*,*) &
         & "feSpacexk_mod   : extract_component"

    if (.NOT.valid(y)) call quit(&
         &"feSpacexk_mod: extract_component:&
         & fespacsexk 'Y' not valid" )

    if ( (c <1) .OR. (c> y%k) ) call quit(&
         &"feSpacexk_mod: extract_component:&
         & component 'c' incorrect" )

    if ( size(u,1) /= y%nbDof2 ) call quit(&
         &"feSpacexk_mod: extract_component:&
         & wrong size for the finite element function 'u'" )

    call allocmem(u_c, y%X_h%nbDof)

    k = y%k
    ll = c
    do jj=1, y%X_h%nbDof
       u_c(jj) = u(ll)
       ll      = ll + k
    end do

  end subroutine extract_component


  function l2_dist_2(uh, Y, qdm, u_1, u_2, u_3) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: uh
    type(fespacexk)       , intent(in) :: Y
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor)                   :: u_1, u_2
    procedure(r3tor)      , optional   :: u_3

    real(RP), dimension(:), allocatable :: uh_i
    real(RP) :: itg

    if (.NOT.valid(y)) call quit(&
         &"feSpacexk_mod: L2_dist_2:&
         & fespacsexk 'Y' not valid" )
    if (.NOT.valid(qdm)) call quit(&
         &"feSpacexk_mod: L2_dist_2:&
         & quadMesh 'qdm' not valid" )
    if(.NOT.associated( qdm%mesh, y%X_h%mesh)) call quit(  &
         &"feSpacexk_mod: L2_dist_2:&
         & 'Y' and 'qdmh' associated to different meshes" )
    if (size(uh,1) /= y%nbDof2) call quit(  &
         &"feSpacexk_mod: L2_dist_2:&
         & wrong size for the finite element function 'uh'")
    if ( y%k==1 ) call quit(&
         &"feSpacexk_mod: L2_dist_2:&
         & dim 1 case, use integral::L2_dist" )
    if ( y%k>3 ) call quit(&
         &"feSpacexk_mod: L2_dist_2: dim > 3 " )
    if (( y%k==3 ) .AND. (.NOT.present(u_3)))  call quit(&
         &"feSpacexk_mod: L2_dist_2:&
         & argument 'u_3' not provided " )

    if (choral_verb>1) write(*,*) &
         & 'feSpacexk_mod   : L2_dist_2'

    !! first component
    call extract_component(uh_i, uh, y, 1)
    itg = l2_dist(u_1, uh_i, y%X_h, qdm)
    dist = itg**2

    !! second component
    call extract_component(uh_i, uh, y, 2)
    itg = l2_dist(u_2, uh_i, y%X_h, qdm)
    dist = dist + itg**2

   !! third component
    if (y%k==3) then
       call extract_component(uh_i, uh, y, 3)
       itg = l2_dist(u_3, uh_i, y%X_h, qdm)
       dist = dist + itg**2
    end if

    dist = sqrt(dist)

  end function l2_dist_2


  function l2_dist_grad_2(uh, Y, qdm, grad_u1,&
       & grad_u2, grad_u3) result(dist)
    real(RP)                           :: dist
    real(RP), dimension(:), intent(in) :: uh
    type(fespacexk)       , intent(in) :: Y
    type(quadmesh)        , intent(in) :: qdm
    procedure(r3tor3)                  :: grad_u1, grad_u2
    procedure(r3tor3)     , optional   :: grad_u3

    real(RP), dimension(:), allocatable :: uh_i
    real(RP) :: itg

    if (.NOT.valid(y)) call quit(&
         &"feSpacexk_mod: L2_dist_grad_2:&
         & fespacsexk 'Y' not valid" )
    if (.NOT.valid(qdm)) call quit(&
         &"feSpacexk_mod: L2_dist_grad_2:&
         & quadMesh 'qdm' not valid" )
    if(.NOT.associated( qdm%mesh, y%X_h%mesh)) call quit(  &
         &"feSpacexk_mod: L2_dist_grad_2:&
         & 'Y' and 'qdmh' associated to different meshes" )
    if (size(uh,1) /= y%nbDof2) call quit(  &
         &"feSpacexk_mod: L2_dist_grad_2:&
         & wrong size for the finite element function 'uh'")
    if ( y%k==1 ) call quit(&
         &"feSpacexk_mod: L2_dist_grad_2:&
         & dim 1 case, use integral::L2_dist_grad" )
    if ( y%k>3 ) call quit(&
         &"feSpacexk_mod: L2_dist_grad_2: dim > 3 " )
    if (( y%k==3 ) .AND. (.NOT.present(grad_u3)))  call quit(&
         &"feSpacexk_mod: L2_dist_grad_2:&
         & argument 'grad_u3' not provided " )

    if (choral_verb>1) write(*,*) &
         & 'feSpacexk_mod   : L2_dist_grad_2'

    !! first component
    call extract_component(uh_i, uh, y, 1)
    itg = l2_dist_grad(grad_u1, uh_i, y%X_h, qdm)
    dist = itg**2

    !! second component
    call extract_component(uh_i, uh, y, 2)
    itg = l2_dist_grad(grad_u2, uh_i, y%X_h, qdm)
    dist = dist + itg**2

   !! third component
    if (y%k==3) then
       call extract_component(uh_i, uh, y, 3)
       itg = l2_dist_grad(grad_u3, uh_i, y%X_h, qdm)
       dist = dist + itg**2
    end if

    dist = sqrt(dist)

  end function l2_dist_grad_2



end module fespacexk_mod