elasticity_Dirichlet_2D.f90

Go to the documentation of this file.

program elasticity_dirichlet_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)  
  !!   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
  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
  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=1)
  write(*,*)
  write(*,*)'elasticity_Dirichlet_2D: start'
  write(*,*)""
  write(*,*)"    EXAMPLE OF A LINEAR ELASTICITY PROBLEM"
  write(*,*)""
  write(*,*)"            -div(A(x) e(u)) = f   on  \Omega"
  write(*,*)""
  write(*,*)"    WITH A DIRICHLET 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' and 'f' "
  write(*,*)"    have been derived in "
  write(*,*)"    choral/maple/example_elasticity_Dirichlet.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(*,*)"================================== DIRICHLET CONDITION"
  call elasticity_dirichlet(stiff, rhs, y, u_1, u_2, gam) 


  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 
  !!
  if (verb>0) then
     !! output file def.
     shell_com = 'elasticity_Dirichlet_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)

  end if

  !!
  !! 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_u_1, grad_u_2)
  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)

contains 

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

    r = 1.0_rp + cc4(x) * 0.2_rp

  end function lambda

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

    r = 1.0_rp + ss4(x) * 0.2_rp

  end function mu

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

    r = cs(x)

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

    r = sc(x)

  end function u_2



  function grad_u_1(x) result(grad)
    real(RP), dimension(3)             :: grad
    real(RP), dimension(3), intent(in) :: x
    
    grad(1) = -pi * ss(x)
    grad(2) =  pi * cc(x)
    grad(3) = 0._rp

  end function grad_u_1
  function grad_u_2(x) result(grad)
    real(RP), dimension(3)             :: grad
    real(RP), dimension(3), intent(in) :: x
    
    grad(1) =  pi * cc(x)
    grad(2) = -pi * ss(x)
    grad(3) = 0._rp

  end function grad_u_2


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

    r = 6.0_rp

    r = r + 0.2_rp*cc4(x) 
    r = r + 0.4_rp*ss4(x) 
    r = r * cs(x)  

    r = r + cs4(x)*ss(x) * re(8,5)
    r = r - sc4(x)*ss(x) * re(8,5)
    r = r - sc4(x)*cc(x) * re(8,5)

    r = r *  pi**2  

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

    r = 6.0_rp

    r = r + 0.2_rp*cc4(x) 
    r = r + 0.4_rp*ss4(x) 
    r = r * sc(x)  

    r = r + sc4(x)*ss(x) * re(8,5)
    r = r - cs4(x)*cc(x) * re(8,5)
    r = r - cs4(x)*ss(x) * re(8,5)

    r = r *  pi**2  

  end function f_2

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

  end function gam


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

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

  end function cc
  function ss(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

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

  end function ss
  function cs(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

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

  end function cs
  function sc(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

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

  end function sc


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

    r = cos(4.0_rp*pi*x(1) ) * cos(4.0_rp*pi*x(2))

  end function cc4
  function ss4(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

    r = sin(4.0_rp*pi*x(1) ) * sin(4.0_rp*pi*x(2))

  end function ss4
  function cs4(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

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

  end function cs4
  function sc4(x) result(r) 
    real(RP)                           :: r
    real(RP), dimension(3), intent(in) :: x

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

  end function sc4



end program elasticity_dirichlet_2d