funcLib.f90

Go to the documentation of this file.

module funclib

  use real_type
 
  implicit none

  public :: one_r3             
  public :: emetric            
  public :: vector_field_e_x   

  public :: f_c0, f_c2, f_c3
  public :: f_c4, f_c5

contains


  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !! FE_FUNC = Finite element functions
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  function vector_field_e_x(x) result(res)
    real(RP), dimension(3)             :: res
    real(RP), dimension(3), intent(in) :: x

    res = (/1._rp, 0._rp, 0._rp/)

  end function vector_field_e_x


  function one_r3(x) result(res)
    real(RP), dimension(3), intent(in)  :: x
    real(RP)                            :: res
    
    res = 1._rp

  end function one_r3

  function emetric(x, u, v) result(res)
    real(RP), dimension(3), intent(in)  :: x, u, v
    real(RP)                            :: res
    
    res = dot_product( u, v )

  end function emetric



  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!
  !! f_ck = Ck FUNCTIONS with support [-1,1]
  !!
  !! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function f_c0(t) result(f)

    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP) :: r

    r = abs(t)
    if (r>1._rp)  then
       f = 0._rp           
    else
       f = 1._rp - r
    end if

  end function f_c0

  function f_c1(t) result(f)
    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP), parameter :: scal = 15._rp / 16._rp
    real(RP) :: r

    r = abs(t)
    if (r>1._rp) then
       f = 0._rp           
    else
       f = t**2
       f = f**2 - 2._rp*f + 1._rp
       f = f * scal
    end if

  end function f_c1

  function f_c2(t) result(f)

    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP) :: r, r3, r4, r5

    r = abs(t)
    if (r>1._rp)  then
       f = 0._rp           
    else
       r3 = r**3
       r4 = r*r3
       r5 = r*r4          
       f = 1._rp - 6._rp*r5 + 15._rp*r4 - 10._rp*r3
    end if

  end function f_c2

  function f_c3(t) result(f)

    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP) :: r, r4, r5, r6, r7

    r = abs(t)
    if (r>1._rp)    then
       f = 0._rp           
    else

       r4 = r**4  
       r5 = r4*r 
       r6 = r5*r 
       r7 = r6*r 
          
       f = 1._rp + 20._rp*r7 -70._rp*r6 + 84._rp*r5 - 35._rp*r4

    end if

  end function f_c3

  function f_c4(t) result(f)

    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP) :: r, r5, r6, r7, r8, r9

    r = abs(t)
    if (r>1._rp) then
       f = 0._rp           
    else

       r5 = r**5
       r6 = r5*r
       r7 = r6*r
       r8 = r7*r
       r9 = r8*r
       f = 1._rp - 70._rp*r9 +315._rp*r8 - 540._rp*r7 + &
            & 420._rp*r6 - 126._rp*r5

    end if

  end function f_c4


  function f_c5(t) result(f)

    real(RP) :: f
    real(RP), intent(in) :: t

    real(RP) :: r, r2, r3, r4, r5

    f = 0._rp    
    r = abs(t)
    if (r>1._rp) then
       f = 0._rp           
    else
       r2 = r*r
       r3 = r2*r
       r4 = r3*r
       r5 = r4*r
       
       f  = 2._rp * (1._rp-r)**6 * ( &
            &           1._rp    +   6._rp*r  +  21._rp*r2 + &
            &          56._rp*r3 + 126._rp*r4 + 252._rp*r5 )

    end if

  end function f_c5


end module funclib