BarycentricHeirarchicalBasis

Evaluate all modal basis (heirarchical polynomials) on Triangle.

Interface

܀ Interface

INTERFACE
  MODULE PURE FUNCTION BarycentricHeirarchicalBasis_Triangle1(order, &
    & pe1, pe2, pe3, lambda, refTriangle) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order in the cell of triangle, it should be greater than 2
    INTEGER(I4B), INTENT(IN) :: pe1
    !! order of interpolation on edge e1
    INTEGER(I4B), INTENT(IN) :: pe2
    !! order of interpolation on edge e2
    INTEGER(I4B), INTENT(IN) :: pe3
    !! order of interpolation on edge e3
    REAL(DFP), INTENT(IN) :: lambda(:, :)
    !! Barycenteric coordinates
    !! number of rows = 3
    !! number of cols = number of points
    CHARACTER(*), INTENT(IN) :: refTriangle
    !! reference triangle, "BIUNIT", "UNIT"
    REAL(DFP) :: ans( &
      & SIZE(xij, 2), &
      & pe1 + pe2 + pe3 + INT((order - 1) * (order - 2) / 2))
    !!
  END FUNCTION BarycentricHeirarchicalBasis_Triangle1
END INTERFACE

️܀ See example

program main
  use easifembase
  use easifemClasses
  implicit none
  real(dfp), allocatable :: xij(:,:), avec(:), lambda(:,:)
  integer(i4b) :: ii, jj, cnt, n

generate points

    n = 51
    call reallocate(avec, n)
    call reallocate(xij, 2, int((n+1)*(n+2)/2))
    avec= linspace(0.0_DFP, 1.0_DFP, n)
    cnt=0
    do ii = 1, n
      do jj = 1, n-ii+1
        cnt=cnt+1
        xij(1,cnt) = avec(ii)
        xij(2,cnt) = avec(jj)
      end do
    end do

  lambda = BarycentricCoordTriangle(xij, "UNIT")

order=1; pe1=1; pe2=1; pe3=1

  BLOCK
    real(dfp), allocatable :: ans(:,:)
    integer(i4b) :: ii, order, pe1, pe2, pe3
    type( VTKPlot_ ) :: aplot
    character(len=*), parameter :: fname="./results/"
    !!
    order=1; pe1=1; pe2=1; pe3=1
    !!
    ans = BarycentricHeirarchicalBasis_Triangle(&
      & order=order, pe1=pe1, pe2=pe2, pe3=pe3, &
      & lambda=lambda, refTriangle="UNIT")
    !!
    do ii  = 1, size(ans,2)
      call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
        & filename=fname//"Pbary( order=" // tostring(order) // &
          & ", pe1="//tostring(pe1) // &
          & ", pe2="//tostring(pe2) // &
          & ", pe3="//tostring(pe3) // &
          & " )" // tostring(ii) // &
          & ".vtp", &
        & label="P")
    end do
    !!
  END BLOCK

order=4; pe1=4; pe2=4; pe3=4

  BLOCK
    real(dfp), allocatable :: ans(:,:)
    integer(i4b) :: ii, order, pe1, pe2, pe3
    type( VTKPlot_ ) :: aplot
    character(len=*), parameter :: fname="./results/"
    !!
    order=4; pe1=4; pe2=4; pe3=4
    !!
    ans = BarycentricHeirarchicalBasis_Triangle(&
      & order=order, pe1=pe1, pe2=pe2, pe3=pe3, &
      & lambda=lambda, refTriangle="UNIT")
    !!
    do ii  = 1, size(ans,2)
      call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
        & filename=fname//"Pbary( order=" // tostring(order) // &
          & ", pe1="//tostring(pe1) // &
          & ", pe2="//tostring(pe2) // &
          & ", pe3="//tostring(pe3) // &
          & " )" // tostring(ii) // &
          & ".vtp", &
        & label="P")
    end do
    !!
  END BLOCK

order=5; pe1=5; pe2=5; pe3=5

  BLOCK
    real(dfp), allocatable :: ans(:,:)
    integer(i4b) :: ii, order, pe1, pe2, pe3
    type( VTKPlot_ ) :: aplot
    character(len=*), parameter :: fname="./results/"
    !!
    order=5; pe1=5; pe2=5; pe3=5
    !!
    ans = BarycentricHeirarchicalBasis_Triangle(&
      & order=order, pe1=pe1, pe2=pe2, pe3=pe3, &
      & lambda=lambda, refTriangle="UNIT")
    !!
    do ii  = 1, size(ans,2)
      call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
        & filename=fname//"Pbary( order=" // tostring(order) // &
          & ", pe1="//tostring(pe1) // &
          & ", pe2="//tostring(pe2) // &
          & ", pe3="//tostring(pe3) // &
          & " )" // tostring(ii) // &
          & ".vtp", &
        & label="P")
    end do
    !!
  END BLOCK

order=5; pe1=1; pe2=5; pe3=3

  BLOCK
    real(dfp), allocatable :: ans(:,:)
    integer(i4b) :: ii, order, pe1, pe2, pe3
    type( VTKPlot_ ) :: aplot
    character(len=*), parameter :: fname="./results/"
    !!
    order=5; pe1=1; pe2=5; pe3=3
    !!
    ans = BarycentricHeirarchicalBasis_Triangle(&
      & order=order, pe1=pe1, pe2=pe2, pe3=pe3, &
      & lambda=lambda, refTriangle="UNIT")
    !!
    do ii  = 1, size(ans,2)
      call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
        & filename=fname//"Pbary( order=" // tostring(order) // &
          & ", pe1="//tostring(pe1) // &
          & ", pe2="//tostring(pe2) // &
          & ", pe3="//tostring(pe3) // &
          & " )" // tostring(ii) // &
          & ".vtp", &
        & label="P")
    end do
    !!
  END BLOCK

end program main

↢

Interface 2

܀ Interface

INTERFACE
  MODULE PURE FUNCTION BarycentricHeirarchicalBasis_Triangle(order, lambda, &
    & refTriangle) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of approximation on triangle
    REAL(DFP), INTENT(IN) :: lambda(:, :)
    !! Barycenteric coordinates
    !! number of rows = 3
    !! number of cols = number of points
    CHARACTER(*), INTENT(IN) :: refTriangle
    !! reference triangle, "BIUNIT", "UNIT"
    REAL(DFP) :: ans( &
      & SIZE(xij, 2), &
      & INT((order + 1) * (order + 2) / 2))
    !!
  END FUNCTION BarycentricHeirarchicalBasis_Triangle
END INTERFACE

️܀ See example

program main
  use easifembase
  use easifemClasses
  implicit none
  real(dfp), allocatable :: xij(:,:), avec(:), lambda(:,:)
  integer(i4b) :: ii, jj, cnt, n

  n = 51
  call reallocate(avec, n)
  call reallocate(xij, 2, int((n+1)*(n+2)/2))
  avec= linspace(0.0_DFP, 1.0_DFP, n)
  cnt=0
  do ii = 1, n
    do jj = 1, n-ii+1
      cnt=cnt+1
      xij(1,cnt) = avec(ii)
      xij(2,cnt) = avec(jj)
    end do
  end do

  lambda = BarycentricCoordTriangle(xij, "UNIT")

  BLOCK
    real(dfp), allocatable :: ans(:,:)
    integer(i4b) :: ii, order
    type( VTKPlot_ ) :: aplot
    character(len=*), parameter :: fname="./results/"
    !!
    order=5
    !!
    ans = BarycentricHeirarchicalBasis_Triangle(order=order, &
      & lambda=lambda, refTriangle="UNIT")
    !!
    do ii  = 1, size(ans,2)
      call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
        & filename=fname//"Pbary( order=" // tostring(order) // &
          & " )" // tostring(ii) // &
          & ".vtp", &
        & label="P")
    end do
    !!
  END BLOCK
end program main

↢