elasticity_Neumann_2D.f90

Go to the documentation of this file.

program elasticity_neumann_2d

  use choral_constants
  use choral

  implicit none

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

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


  !!      SPACE DISCRETISATION
  !!
  !!   fe_v     = finite element method (volumic)
  !!   fe_s     = finite element method (surfacic)
  !!   qd_v     = quadrature method (volumic)  
  !!   qd_s     = quadrature method (surfacic)      
  !!   mesh_idx = index for the mesh
  !!
  integer         , parameter :: fe_v     = fe_p2_2d
  integer         , parameter :: fe_s     = fe_p2_1d
  integer         , parameter :: qd_v     = quad_gauss_trg_12
  integer         , parameter :: qd_s     = quad_gauss_edg_4
  character(LEN=1), parameter :: mesh_idx = "2"
  !!
  !! m        = mesh
  !! X_h      = finite element space
  !! Y        = [X_h]^dim
  !! qdm      = integration method
  !!
  type(mesh)      :: msh
  type(fespace)   :: X_h
  type(fespacexk) :: Y
  type(quadmesh)  :: qdm
  !!
  !! msh_file 
  character(len=100), parameter :: mesh_file= &
       & trim(GMSH_DIR)//"square/square_"//mesh_idx//".msh"

  
  !!       LINEAR SYSTEM
  !!
  !!  kry   = krylov method def.
  !!  stiff = stiffness matrix
  !!  rhs   = right hand side
  !!  u_h   = numerical solution
  !!
  type(krylov) :: kry
  type(csr)    :: stiff
  real(RP), dimension(:), allocatable :: rhs, u_h, aux
  real(RP), dimension(:), allocatable :: u_h_1, u_h_2


  !!      NUMERICAL ERRORS
  !!
  real(RP) :: err_L2, err_H1

  character(LEN=300) :: shell_com

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


  !! !!!!!!!!!!!!!!!!!!!!!  INITIALISATION
  !!
  call choral_init(verb=verb)
  write(*,*)
  write(*,*)'elasticity_Neumann_2D: start'
  write(*,*)""
  write(*,*)"    EXAMPLE OF A LINEAR ELASTICITY PROBLEM"
  write(*,*)""
  write(*,*)"            -div(A(x) e(u)) = f   on  \Omega"
  write(*,*)""
  write(*,*)"    WITH A NEUMANN BOUNDARY CONDITION "  
  write(*,*)""
  write(*,*)"            A(x)e(u)) n = g   on   \partial\Omega"
  write(*,*)""
  write(*,*)"    AND WITH NON CONSTANT LAME COEFFICIENTS:"
  write(*,*)""
  write(*,*)"    A(x) \xi = lambda(x) Tr(\xi) Id + 2 mu(x) \si"
  write(*,*)""
  write(*,*)"    The problem data 'lambda', 'mu', 'u', 'f' and 'g'"
  write(*,*)"    have been derived in "
  write(*,*)"    choral/maple/example_elasticity_Neumann.mw"
  write(*,*)""
  write(*,*)""
  write(*,*)""

  write(*,*) ""
  write(*,*)"================================== MESH ASSEMBLING"
  msh = mesh(mesh_file, 'gmsh')
  
  call print(msh)

  write(*,*) ""
  write(*,*)"================================== FINITE ELEMENT&
       & SPACE ASSEMBLING"
  x_h =fespace(msh)
  call set(x_h, fe_v)
  call set(x_h, fe_s)
  call assemble(x_h)
  call print(x_h)
  y = fespacexk(x_h, 2)
  call print(y)

  write(*,*) ""
  write(*,*)"================================== INTEGRATION METHOD&
       & ON \Omega ASSEMBLING"
  qdm = quadmesh(msh)
  call set(qdm, qd_v)
  call assemble(qdm)
  call print(qdm)

  write(*,*) ""
  write(*,*)"================================== STIFFNESS&
       & MATRIX ASSEMBLING"
  call elasticity_stiffmat(stiff, lambda, mu, &
       &                         y, qdm)

 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   1- COMPUTATION OF F"
  !!
  !!  \int_\Omega f u_i dx
  call l2_product(rhs, y, qdm, f_1, f_2)

 
  write(*,*) ""
  write(*,*)"================================== RIGHT HAND SIDE"
  write(*,*)"                                   2- COMPUTATION OF G"
  !!
  !!  \int_{Gamma, x=0} g  u_i dx
  call elasticity_neumann_rhs(aux, y, qd_s, &
       &                      g_1_x_0, g_2_x_0, f=x_le_0) 
  rhs = rhs + aux
  
  !!  \int_{Gamma, x=1} g  u_i dx
  call elasticity_neumann_rhs(aux, y, qd_s, &
       &                      g_1_x_1, g_2_x_1, f=x_ge_1) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, y=0} g  u_i dx
  call elasticity_neumann_rhs(aux, y, qd_s, &
       &                      g_1_y_0, g_2_y_0, f=y_le_0) 
  rhs = rhs + aux
  !!
  !!  \int_{Gamma, y=1} g  u_i dx
  call elasticity_neumann_rhs(aux, y, qd_s, &
       &                      g_1_y_1, g_2_y_1, f=y_ge_1) 
  rhs = rhs + aux

  write(*,*) ""
  write(*,*)"================================== SOLVE LINEAR SYSTEM&
       & Stiff U_h = RHS "
  !! 
  !! Solver setting
  kry = krylov(kry_gmres, tol=1e-6_rp, itmax=1000, restart=25, verb=2)
  call print(kry)
  !!
  !! initial guess 
  call interp_vect_func(u_h, y, u_1, u_2)
  !!
  !! linear system inversion
  call solve(u_h, kry, rhs, stiff)
  !!
  !! cleaning
  call clear(stiff)
  call freemem(rhs)

  write(*,*) ""
  write(*,*)"================================== POST-TREATMENT"
  !!
  !! extract component
  !!
  call extract_component(u_h_2, u_h, y, 2)
  call extract_component(u_h_1, u_h, y, 1)

  !!
  !! graphical display 
  !!
  !! output file def.
  shell_com = 'elasticity_Neumann_u.msh'
  !!
  !! first write  the 'finite element mesh'
  call write(x_h, trim(shell_com), 'gmsh')
  !!
  !! second save the finite element solution u_h
  !!
  call gmsh_addview(x_h, u_h_1, &
       & trim(shell_com), &
       &'Elasticity problem : numerical solution', 0.0_rp, 1)
  call gmsh_addview(x_h, u_h_2, &
       & trim(shell_com), &
       &'Elasticity problem : numerical solution', 1.0_rp, 2)
  !!
  !! shell command to visualise with gmsh
  shell_com = 'gmsh -option '//trim(gmsh_dir)//'gmsh-options-view'&
       & //' '//trim(shell_com)
  call system(shell_com)

  !!
  !! L2 and H1_0 numerical errors
  write(*,*) ""
  write(*,*) "Numerical errors"
  err_l2 = l2_dist(u_h, y, qdm, u_1, u_2)
  err_h1 = l2_dist_grad(u_h, y, qdm, grad_u1, grad_u2)
  write(*,*)"  Error | u - u_h|_L2            =", err_l2
  write(*,*)"  Error |\nabla (u - u_h)|_L2    =", err_h1

  call freemem(u_h)
  call freemem(u_h_1)
  call freemem(u_h_2)
  call freemem(aux)

contains 

  function lambda(x) result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = 1.0_rp

  end function lambda

  function mu(x) result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = 1.0_rp

  end function mu

  function u_1(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = cos(pi*x(1) ) * sin(pi*x(2))

  end function u_1
  function u_2(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = sin(pi*x(1) ) * cos(pi*x(2))

  end function u_2



  function grad_u1(x) result(grad)
    real(RP), dimension(3)             :: grad
    real(RP), dimension(3), intent(in) :: x
    
    grad(1) = -pi * sin(pi*x(1)) * sin(pi*x(2))
    grad(2) =  pi * cos(pi*x(1)) * cos(pi*x(2))
    grad(3) = 0._rp

  end function grad_u1
  function grad_u2(x) result(grad)
    real(RP), dimension(3)             :: grad
    real(RP), dimension(3), intent(in) :: x
    
    grad(1) =  pi * cos(pi*x(1)) * cos(pi*x(2))
    grad(2) = -pi * sin(pi*x(1)) * sin(pi*x(2))
    grad(3) = 0._rp

  end function grad_u2


  function f_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = 6.0_rp * pi**2 * u_1(x)

  end function f_1
  function f_2(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = 6.0_rp * pi**2 * u_2(x)

  end function f_2


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

  end function x_le_0


  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

  function y_le_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = -x(2)

  end function y_le_0


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

  end function y_ge_1


  function g_1_x_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  0.0_rp

  end function g_1_x_0

  function g_2_x_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r = -2.0_rp * cos(pi*x(2)) * pi

  end function g_2_x_0

  function g_1_x_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  0.0_rp

  end function g_1_x_1

  function g_2_x_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  -2.0_rp * cos(pi*x(2)) * pi

  end function g_2_x_1

  function g_1_y_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  -2.0_rp * cos(pi*x(1)) * pi

  end function g_1_y_0

  function g_2_y_0(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  0.0_rp

  end function g_2_y_0

  function g_1_y_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  -2.0_rp * cos(pi*x(1)) * pi

  end function g_1_y_1

  function g_2_y_1(x)  result(r)
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x
  
    r =  0.0_rp

  end function g_2_y_1


end program elasticity_neumann_2d