quad_mod.f90

Go to the documentation of this file.

module quad_mod

  use choral_constants
  use choral_env
  use real_type
  use io, only: read
  use cell_mod
  
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%

  public :: quad_init          ! TESTED

  public :: quad_name
  public :: quad_nbnodes
  public :: quad_geo
  public :: quad_order
  public :: quad_dim
  public :: quad_coord
  public :: quad_wgt

  !       %----------------------------------------%
  !       |                                        |
  !       |               QUAD  ARRAYS             |
  !       |                                        |
  !       %----------------------------------------%


  character(LEN=13), dimension(0:QUAD_TOT_NB) :: quad_name

  integer , dimension(QUAD_TOT_NB):: quad_nbnodes

  integer , dimension(QUAD_TOT_NB):: quad_geo

  integer , dimension(QUAD_TOT_NB):: quad_dim

  integer , dimension(QUAD_TOT_NB):: quad_order

  type(r_2d), dimension(QUAD_TOT_NB):: quad_coord

  type(r_1d), dimension(QUAD_TOT_NB):: quad_wgt


contains
  

  subroutine quad_init(b)
    logical, intent(in) :: b

    integer :: qt, geo, dim, nn

    if (b) write(*,*) "quad_mod        : quad_init"

    !! set quad name, geo, order, nbNodes
    !!
    call def_quad_xxx()

    !! set quad dimension
    !!
    do qt=1, quad_tot_nb
       geo = quad_geo(qt)
       dim = cell_dim(geo)
       quad_dim(qt) = dim
    end do

    !! allocate memory to store node coordinates
    !!                          and weights 
    !!
    do qt=1, quad_tot_nb

       nn  = quad_nbnodes(qt)
       dim = quad_dim(qt)

       quad_coord(qt) = r_2d( dim, nn)
       quad_wgt(qt)   = r_1d( nn)

    end do

    !! set GUAD_COORD AND QUAD_WGT
    !!
    call def_quad_data()

  end subroutine quad_init

  subroutine def_quad_xxx()

    quad_name(quad_none   ) = "VOID"
    quad_name(quad_gauss_edg_1 ) = "Gauss EDG 1  " 
    quad_name(quad_gauss_edg_2 ) = "Gauss EDG 2  " 
    quad_name(quad_gauss_edg_3 ) = "Gauss EDG 3  " 
    quad_name(quad_gauss_edg_4 ) = "Gauss EDG 4  " 
    quad_name(quad_gauss_trg_1 ) = "Gauss TRG 1  " 
    quad_name(quad_gauss_trg_3 ) = "Gauss TRG 3  " 
    quad_name(quad_gauss_trg_6 ) = "Gauss TRG 6  " ! STRANG_5 
    quad_name(quad_gauss_trg_12) = "Gauss TRG 12 " ! STRANG_9
    quad_name(quad_gauss_trg_13) = "Gauss TRG 13 " ! STRANG_10
    quad_name(quad_gauss_trg_19) = "Gauss TRG 19 " ! TOMS584_19
    quad_name(quad_gauss_trg_28) = "Gauss TRG 28 " ! TOMS612_28
    quad_name(quad_gauss_trg_37) = "Gauss TRG 37 " ! TOMS706_37
    quad_name(quad_gauss_tet_1 ) = "Gauss TET 1  " 
    quad_name(quad_gauss_tet_4 ) = "Gauss TET 4  " ! keast1 1
    quad_name(quad_gauss_tet_15) = "Gauss TET 15 " ! keast1 6 
    quad_name(quad_gauss_tet_31) = "Gauss TET 31 " ! keast1 8
    quad_name(quad_gauss_tet_45) = "Gauss TET 45 " ! keast1 9 

    quad_geo(quad_gauss_edg_1 ) = cell_edg
    quad_geo(quad_gauss_edg_2 ) = cell_edg
    quad_geo(quad_gauss_edg_3 ) = cell_edg
    quad_geo(quad_gauss_edg_4 ) = cell_edg
    quad_geo(quad_gauss_trg_1 ) = cell_trg
    quad_geo(quad_gauss_trg_3 ) = cell_trg
    quad_geo(quad_gauss_trg_6 ) = cell_trg
    quad_geo(quad_gauss_trg_12) = cell_trg
    quad_geo(quad_gauss_trg_13) = cell_trg
    quad_geo(quad_gauss_trg_19) = cell_trg
    quad_geo(quad_gauss_trg_28) = cell_trg
    quad_geo(quad_gauss_trg_37) = cell_trg
    quad_geo(quad_gauss_tet_1 ) = cell_tet
    quad_geo(quad_gauss_tet_4 ) = cell_tet
    quad_geo(quad_gauss_tet_15) = cell_tet
    quad_geo(quad_gauss_tet_31) = cell_tet
    quad_geo(quad_gauss_tet_45) = cell_tet

    quad_nbnodes(quad_gauss_edg_1 ) = 1
    quad_nbnodes(quad_gauss_edg_2 ) = 2
    quad_nbnodes(quad_gauss_edg_3 ) = 3
    quad_nbnodes(quad_gauss_edg_4 ) = 4
    quad_nbnodes(quad_gauss_trg_1 ) = 1
    quad_nbnodes(quad_gauss_trg_3 ) = 3
    quad_nbnodes(quad_gauss_trg_6 ) = 6
    quad_nbnodes(quad_gauss_trg_12) = 12
    quad_nbnodes(quad_gauss_trg_13) = 13
    quad_nbnodes(quad_gauss_trg_19) = 19
    quad_nbnodes(quad_gauss_trg_28) = 28
    quad_nbnodes(quad_gauss_trg_37) = 37
    quad_nbnodes(quad_gauss_tet_1 ) = 1 
    quad_nbnodes(quad_gauss_tet_4 ) = 4
    quad_nbnodes(quad_gauss_tet_15) = 15  
    quad_nbnodes(quad_gauss_tet_31) = 31  
    quad_nbnodes(quad_gauss_tet_45) = 45  

    quad_order(quad_gauss_edg_1 ) = 1 
    quad_order(quad_gauss_edg_2 ) = 3 
    quad_order(quad_gauss_edg_3 ) = 5 
    quad_order(quad_gauss_edg_4 ) = 7 
    quad_order(quad_gauss_trg_1 ) = 1 
    quad_order(quad_gauss_trg_3 ) = 2 
    quad_order(quad_gauss_trg_6 ) = 4 
    quad_order(quad_gauss_trg_12) = 6 
    quad_order(quad_gauss_trg_13) = 7
    quad_order(quad_gauss_trg_19) = 8
    quad_order(quad_gauss_trg_28) = 10  ! under converging /= 11
    quad_order(quad_gauss_trg_37) = 10  ! under converging /= 13
    quad_order(quad_gauss_tet_1 ) = 1
    quad_order(quad_gauss_tet_4 ) = 2
    quad_order(quad_gauss_tet_15) = 5  
    quad_order(quad_gauss_tet_31) = 7  
    quad_order(quad_gauss_tet_45) = 8  

  end subroutine def_quad_xxx


  subroutine def_quad_data()

    real(RP) :: a
    integer  :: qt

    qt = quad_gauss_edg_1
    quad_wgt(qt)%y = 1.0_rp
    quad_coord(qt)%y = 0.5_rp

    qt = quad_gauss_edg_2
    quad_wgt(qt)%y = 0.5_rp
    quad_coord(qt)%y(1,1) = 0.21132486540518702_rp
    quad_coord(qt)%y(1,2) = 0.788675134594813_rp

    qt = quad_gauss_edg_3
    quad_wgt(qt)%y(1) = 0.44444444444444444444444444444444444_rp
    quad_wgt(qt)%y(2) = 0.27777777777777777777777777777777778_rp
    quad_wgt(qt)%y(3) = 0.27777777777777777777777777777777778_rp
    quad_coord(qt)%y(1,1) = 0.5_rp
    quad_coord(qt)%y(1,2) = 0.11270166537925852_rp
    quad_coord(qt)%y(1,3) = 0.8872983346207415_rp

    qt = quad_gauss_edg_4
    quad_wgt(qt)%y(1) = 0.17392742256872692868653197461099970_rp
    quad_wgt(qt)%y(2) = 0.17392742256872692868653197461099970_rp
    quad_wgt(qt)%y(3) = 0.32607257743127307131346802538900030_rp
    quad_wgt(qt)%y(4) = 0.32607257743127307131346802538900030_rp
    a = 0.06943184420297371238802675555359525_rp
    quad_coord(qt)%y(1,1) = a
    quad_coord(qt)%y(1,2) = 1._rp - a
    a = 0.33000947820757186759866712044837766_rp
    quad_coord(qt)%y(1,3) = a
    quad_coord(qt)%y(1,4) = 1._rp-a

    qt = quad_gauss_trg_1
    quad_wgt(qt)%y = 0.5_rp 
    quad_coord(qt)%y(:,1) = re(1, 3)

    qt = quad_gauss_trg_3
    quad_wgt(qt)%y = re(1, 6)
    quad_coord(qt)%y(:,1) = re(1, 6)
    quad_coord(qt)%y(:,2) = (/ re(2, 3), re(1, 6)/)
    quad_coord(qt)%y(:,3) = (/ re(1, 6), re(2, 3)/)

    qt = quad_gauss_trg_6
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang5_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y * 0.5_rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang5_x.txt')


    qt = quad_gauss_trg_12
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang9_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y*0.5_rp    
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang9_x.txt')

    qt = quad_gauss_trg_13
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang10_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y*0.5_rp    
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/strang10_x.txt')

    qt = quad_gauss_trg_19
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/TOMS584_19_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y*0.5_rp    
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/TOMS584_19_x.txt')


    qt = quad_gauss_trg_28
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/toms612_28_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y*0.5_rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/toms612_28_x.txt')

    qt = quad_gauss_trg_37
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/toms706_37_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y*0.5_rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/toms706_37_x.txt')


    qt = quad_gauss_tet_1
    quad_wgt(qt)%y = re(1, 6)
    quad_coord(qt)%y = re(1, 4)

    qt = quad_gauss_tet_4
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast1_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y / 6._rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast1_x.txt')

    qt = quad_gauss_tet_15
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast6_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y / 6._rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast6_x.txt')

    qt = quad_gauss_tet_31
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast8_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y / 6._rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast8_x.txt')

    qt = quad_gauss_tet_45
    call read(quad_wgt(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast9_w.txt')
    quad_wgt(qt)%y = quad_wgt(qt)%y / 6._rp
    call read(quad_coord(qt)%y, &
         & trim(choral_dir)//'SRC/quad_data/keast9_x.txt')

  end subroutine def_quad_data


end module quad_mod