Dubiner

Forms Dubiner basis on reference triangle domain.

Reference triangle can be "biunit" or "unit".

The shape of ans is (M,N), where M=SIZE(xij,2) (number of points) N = 0.5*(order+1)*(order+2).

In this way, ans(j,:) denotes the values of all polynomial at jth point.

Polynomials are returned in following way:

$$P_{0,0}, P_{0,1}, \cdots , P_{0,order} \\ P_{1,0}, P_{1,1}, \cdots , P_{1,order-1} \\ P_{2,0}, P_{2,1}, \cdots , P_{2,order-2} \\ \cdots P_{order,0}$$

For example for order=3, the polynomials are arranged as:

$$P_{0,0}, P_{0,1}, P_{0,2}, P_{0,3} \\ P_{1,0}, P_{1,1}, P_{1,2} \\ P_{2,0}, P_{2,1} \\ P_{3,0}$$

Interface 1

܀ Interface

INTERFACE
  MODULE PURE FUNCTION Dubiner_Triangle(order, xij, refTriangle) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial space
    REAL(DFP), INTENT(IN) :: xij(:, :)
    !! points in reference triangle, shape functions will be evaluated
    !! at these points. SIZE(xij,1) = 2, and SIZE(xij, 2) = number of points
    CHARACTER(*), INTENT(IN) :: refTriangle
    !! "unit"
    !! "biunit"
    REAL(DFP) :: ans(SIZE(xij, 2), (order + 1)* (order + 2) / 2)
    !! shape functions
    !! ans(:, j), jth shape functions at all points
    !! ans(j, :), all shape functions at jth point
  END FUNCTION Dubiner_Triangle
END INTERFACE

️܀ order=1

program main
  use easifembase
  use easifemClasses
  implicit none
  real(dfp), allocatable :: xij(:,:), avec(:)
  integer(i4b) :: ii, jj, cnt, n
  real(dfp), allocatable :: ans(:,:)
  integer(i4b) :: order
  type( VTKPlot_ ) :: aplot
  character(len=*), parameter :: fname="./results/"

  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

  order=1
  ans = Dubiner_Triangle(&
    & order=order, xij=xij, &
    & refTriangle="UNIT")

  do ii  = 1, size(ans,2)
    call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
      & filename=fname//"Dubiner( order=" // tostring(order) // &
        & " )" // tostring(ii) // &
        & ".vtp", &
      & label="P")
  end do

end program main

Results:

️܀ order=2

program main
  use easifembase
  use easifemClasses
  implicit none
  real(dfp), allocatable :: xij(:,:), avec(:)
  integer(i4b) :: ii, jj, cnt, n
  real(dfp), allocatable :: ans(:,:)
  integer(i4b) :: order
  type( VTKPlot_ ) :: aplot
  character(len=*), parameter :: fname="./results/"

  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

  order=2
  ans = Dubiner_Triangle(&
    & order=order, xij=xij, &
    & refTriangle="UNIT")

  do ii  = 1, size(ans,2)
    call aplot%scatter3D(x=xij(1,:), y=xij(2, :), z=ans(:,ii), &
      & filename=fname//"Dubiner( order=" // tostring(order) // &
        & " )" // tostring(ii) // &
        & ".vtp", &
      & label="P")
  end do

end program main

↢

Interface 2

܀ Interface

INTERFACE
  MODULE PURE FUNCTION Dubiner_Triangle(order, x, y, refTriangle) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial space
    REAL(DFP), INTENT(IN) :: x(:), y(:)
    !! x and y coordinates, total points = SIZE(x)*SIZE(y)
    CHARACTER(*), INTENT(IN) :: refTriangle
    !! "unit"
    !! "biunit"
    REAL(DFP) :: ans(SIZE(x) * SIZE(y), (order + 1) * (order + 2) / 2)
    !! shape functions
    !! ans(:, j), jth shape functions at all points
    !! ans(j, :), all shape functions at jth point
  END FUNCTION Dubiner_Triangle
END INTERFACE

Forms Dubiner basis on reference triangle domain. Reference triangle can be biunit or unit. Here x and y are the coordinates on the line. xij is given by outerproduct of x and y.

↢