quadMesh_integral.f90

Go to the documentation of this file.

program quadmesh_integral

  use choral
  use choral_constants

  implicit none

  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !!      VARAIBLE DEFINITION : BEGIN
  !!
  !!  

  !! verb    = verbosity level
  !!
  integer, parameter :: verb = 2

  !! msh_file  = mesh file name
  !! msh       = mesh 
  !! qdm       = integration method on the mesh
  !!
  character(len=150)      :: msh_file
  type(mesh)              :: msh
  type(quadmesh)          :: qdm

  !! integral result
  real(RP) :: int

  !!
  !!   VARAIBLE DEFINITION : END
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  call choral_init(verb=verb)

  write(*,*)
  write(*,*)'quadMesh_integral: start'
  write(*,*)""
  write(*,*)"    INTEGRAL COMPITATION ON A MESH OF [0,1]x[0,1]"
  write(*,*)"      - \Omega   = [0,1]x[0,1] "
  write(*,*)"      - \Gamma   = \partial \Omega \cap {x=1} "
  write(*,*)"      - \Gamma_r = \Gamma \cap {x=1} "
  write(*,*)""

  !!
  !! mesh file
  !!
  msh_file = trim(gmsh_dir)//"square/square_3.msh"


  write(*,*) ""
  write(*,*)"================================== ASSEMBLING THE MESH "
  msh = mesh(trim(msh_file), "gmsh")
  if (verb>1) call print(msh)

  write(*,*) ""
  write(*,*)"================================== INTEGRAL 1"
  write(*,*)"=   "
  write(*,*)"=   integration domain = \Omega "
  write(*,*)"=   "
  !!
  !! link 'qdm' to the mesh 'msh'
  !!
  qdm = quadmesh(msh)
  !!
  !! set all triangle cells to the 
  !! 6 point Gaussian quadrature rule on triangles
  !!
  call set(qdm, quad_gauss_trg_6)
  !!
  !! ends up the construction and display a short description
  !!
  call assemble(qdm)
  if (verb>1) call print(qdm)
  !!
  !! integral computation
  !!
  int = integ(func, qdm)
  print*, "   COMPUTED INTEGRAL = ", int
  print*, "   NUMERICAL ERROR   = ", &
       &   abs(int - 2._rp/3._rp)



  write(*,*) ""
  write(*,*)"================================== INTEGRAL 2"
  write(*,*)"=   "
  write(*,*)"=   integration domain = \Gamma "
  write(*,*)"=   "
  !!
  !! re-initialise 'qdm'
  !!
  qdm = quadmesh(msh)
  !!
  !! set all edge cells to the 
  !! 3 point Gaussian quadrature rule on edges
  !!
  call set(qdm, quad_gauss_edg_3)
  !!
  !! ends up the construction and display a short description
  !!
  call assemble(qdm)
  if (verb>1) call print(qdm)
  !!
  !! integral computation
  !!
  int = integ(func, qdm)
  print*, "   COMPUTED INTEGRAL = ", int
  print*, "   NUMERICAL ERROR   = ", &
       &   abs(int - 10._rp/3._rp)


  write(*,*) ""
  write(*,*)"================================== INTEGRAL 3"
  write(*,*)"=   "
  write(*,*)"=   integration domain = \Gamma_r = \Gamma \ cap {x >= 1} "
  write(*,*)"=   "
  !!
  !! re-initialise 'qdm'
  !!
  qdm = quadmesh(msh)
  !!
  !! set all edge cells where ' x >= 1' to the 
  !! 3 point Gaussian quadrature rule on edges
  !!
  call set(qdm, quad_gauss_edg_3, x_ge_1)
  !!
  !! ends up the construction and display a short description
  !!
  call assemble(qdm)
  if (verb>1) call print(qdm)
  !!
  !! integral computation
  !!
  int = integ(func, qdm)
  print*, "   COMPUTED INTEGRAL = ", int
  print*, "   NUMERICAL ERROR   = ", &
       &   abs(int - 4._rp/3._rp)


  write(*,*) ""
  write(*,*)"================================== INTEGRAL 4"
  write(*,*)"=   "
  write(*,*)"=   integration domain = \Gamma_r \cup \Omega "
  write(*,*)"=   (composite integration domain)"
  write(*,*)"=   "
  !!
  !! re-initialise 'qdm'
  !!
  qdm = quadmesh(msh)
  !!
  !! set all edge cells where ' x >= 1' to the 
  !! 3 point Gaussian quadrature rule on edges
  !!
  call set(qdm, quad_gauss_edg_3, x_ge_1)
  !!
  !! set all triangle cells to the 
  !! 6 point Gaussian quadrature rule on triangles
  !!
  call set(qdm, quad_gauss_trg_6)
  !!
  !! ends up the construction and display a short description
  !!
  call assemble(qdm)
  if (verb>1) call print(qdm)
  !!
  !! integral computation
  !!
  int = integ(func, qdm)
  print*, "   COMPUTED INTEGRAL =", int
  print*, "   NUMERICAL ERROR   = ", &
       &   abs(int - 2._rp)



contains

  !! The function to be integrated
  !!
  function func(x)
    real(RP) :: func
    real(RP), dimension(3), intent(in) :: x

    func = x(1)**2 + x(2)**2

  end function func

  function x_ge_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = x(1) - 1.0_rp

  end function x_ge_1

  
end program quadmesh_integral