InterpolationPoint

Get the interpolation point on 1D/2D/3D elements.

Calling example

x(:,:) = InterpolationPoint(order, elemType, ipType, xij)

• order order of polynomial
• elemType element type
• ipType interpolation point
  • Equidistance
  • GaussLegendre
  • GaussLegendreLobatto
  • GaussChebyshev
  • GaussChebyshevLobatto
  • GaussJacobi
  • GaussJacobiLobatto

Interface

INTERFACE
  MODULE FUNCTION InterpolationPoint(order, elemType, ipType, &
    & xij, layout) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of interpolation
    INTEGER(I4B), INTENT(IN) :: elemType
    !! element type
    INTEGER(I4B), INTENT(IN) :: ipType
    !! interpolation point type
    !! Equidistance, GaussLegendre, GaussLegendreLobatto, GaussChebyshev,
    !! GaussChebyshevLobatto, GaussJacobi, GaussJacobiLobatto
    REAL(DFP), OPTIONAL, INTENT(IN) :: xij(:, :)
    !! domain of interpolation, default values are given by
    !! line = [-1,1]
    !! triangle = (0,0), (0,1), (1,0)
    !! quadrangle = [-1,1]x[-1,1]
    CHARACTER(*), INTENT(IN) :: layout
    !! "VEFC" Vertex, Edge, Face, Cell
    !! "INCREASING" incresing order
    !! "DECREASING" decreasing order
    !! "XYZ" First X, then Y, then Z
    !! "YXZ" First Y, then X, then Z
    REAL(DFP), ALLOCATABLE :: ans(:, :)
    !! interpolation points in xij format
  END FUNCTION InterpolationPoint
END INTERFACE

Line

️܀ Equidistance

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=Equidistance, layout=VEFC

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  ipType = Equidistance
  order=0
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)
  !!
  !!
  order=1
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)
  !!
  !!
  order=2
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)
  !!
  !!
  order=3
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)
  !!
  !!
  order=5
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)

See results

ipType=Equidistance, layout=VEFC

5.5

ipType=Equidistance, layout=VEFC

1	10

ipType=Equidistance, layout=VEFC

1	10	5.5

ipType=Equidistance, layout=VEFC

1	10	4	7

ipType=Equidistance, layout=VEFC

1	10	2.8	4.6	6.4	8.2

END PROGRAM main

܀ Equidistance INCREASING

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=Equidistance, layout=INCREASING

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  order=0
  ipType = Equidistance
  layout="INCREASING"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=INCREASING" // char_lf // char_lf)

  order=1
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=INCREASING" // char_lf // char_lf)

  order=2
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=INCREASING" // char_lf // char_lf)

  order=3
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=Equidistance, layout=INCREASING" // char_lf // char_lf)

See results

ipType=Equidistance, layout=INCREASING

5.5

ipType=Equidistance, layout=INCREASING

1	10

ipType=Equidistance, layout=INCREASING

1	5.5	10

ipType=Equidistance, layout=INCREASING

1	4	7	10

END PROGRAM main

↢

܀ GaussLegendreLobatto VEFC

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=GaussLegendreLobatto, layout=VEFC

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  ipType = GaussLegendreLobatto
  order=0
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto, layout=VEFC" // char_lf // char_lf)
  !!
  !!
  order=1
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto, layout=VEFC" // char_lf //char_lf)
  !!
  !!
  order=5
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto, layout=VEFC" // char_lf // char_lf)

See results

ipType=GaussLegendreLobatto, layout=VEFC

5.5

ipType=GaussLegendreLobatto, layout=VEFC

1	10

ipType=GaussLegendreLobatto, layout=VEFC

1	10	2.0573	4.2165	6.7835	8.9427

END PROGRAM main

܀ GaussLegendreLobatto INCREASING

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=GaussLegendreLobatto layout=INCREASING

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  order=0
  ipType = GaussLegendreLobatto
  layout="INCREASING"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto layout=INCREASING" // char_lf //char_lf)
  !
  order=1
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto layout=INCREASING" // char_lf //char_lf)
  !
  order=5
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussLegendreLobatto layout=INCREASING" // char_lf //char_lf)

See results

ipType=GaussLegendreLobatto layout=INCREASING

5.5

ipType=GaussLegendreLobatto layout=INCREASING

1	10

ipType=GaussLegendreLobatto layout=INCREASING

1	2.0573	4.2165	6.7835	8.9427	10

END PROGRAM main

↢

GaussChebyshevLobatto VEFC

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=GaussChebyshevLobatto, layout=VEFC

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  ipType = GaussChebyshevLobatto
  order=0
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)
  !
  !
  order=1
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)
  !
  !
  order=5
  layout="VEFC"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)

See results

ipType=GaussChebyshevLobatto, layout=VEFC

0

ipType=GaussChebyshevLobatto, layout=VEFC

-1	1

ipType=GaussChebyshevLobatto, layout=VEFC

-1	1	-0.80902	-0.30902	0.30902	0.80902

END PROGRAM main

GaussChebyshevLobatto INCREASING

PROGRAM main
use easifemBase
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Line
CHARACTER(len=20) :: layout
real(dfp), parameter :: tol=1.0E-10

ipType=GaussChebyshevLobatto, layout=INCREASING

  xij = zeros(1,2, 0.0_DFP)
  xij(1,:) = [1.0, 10.0]
  ipType = GaussChebyshevLobatto
  order=0
  layout="INCREASING"
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)
  !!
  !!
  order=1
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)
  !!
  !!
  order=5
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(ip), &
    & "ipType=GaussChebyshevLobatto, layout=VEFC" // char_lf // char_lf)

See results

result

ipType=GaussChebyshevLobatto, layout=VEFC

5.5

ipType=GaussChebyshevLobatto, layout=VEFC

1	10

ipType=GaussChebyshevLobatto, layout=VEFC

1	1.8594	4.1094	6.8906	9.1406	10

END PROGRAM main

↢

More

️܀ BlythPozLegendre

This example shows the use of InterpolationPoint routine.

• In this example we test BlythPozLegendre type triangle nodes.

PROGRAM main
use easifemBase
use easifemclasses
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order, ipType
integer(i4b), parameter :: elemType=Triangle
CHARACTER(*), PARAMETER :: layout = "VEFC", fname="./results/test3"
real(dfp), parameter :: tol=1.0E-10
type(csvfile_) :: afile

Setup

  xij = zeros(2,3, 0.0_DFP)
  xij(1,:) = [0.0_DFP, 1.0_DFP, 0.0_DFP]
  xij(2,:) = [0.0_DFP, 0.0_DFP, 1.0_DFP]

ipType=Equidistance

  ipType = Equidistance
  !!
  order=5
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(TRANSPOSE(ip)), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)

See results

ipType=Equidistance, layout=VEFC

0	0
1	0
0	1
0.2	0
0.4	0
0.6	0
0.8	0
0.8	0.2
0.6	0.4
0.4	0.6
0.2	0.8
0	0.8
0	0.6
0	0.4
0	0.2
0.2	0.2
0.6	0.2
0.2	0.6
0.4	0.2
0.4	0.4
0.2	0.4

ipType=BlythPozLegendre, layout=VEFC

  ipType = BlythPozLegendre
  !!
  order=4
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(TRANSPOSE(ip)), &
    & "ipType=Equidistance, layout=VEFC" // char_lf)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

See results

-3.70074E-17	-3.70074E-17
1	-3.70074E-17
-3.70074E-17	1
0.17267	-7.40149E-17
0.5	-7.40149E-17
0.82733	-9.25186E-17
0.82733	0.17267
0.5	0.5
0.17267	0.82733
-9.25186E-17	0.82733
-7.40149E-17	0.5
-7.40149E-17	0.17267
0.22422	0.22422
0.55155	0.22422
0.22422	0.55155

ipType=BlythPozLegendre, layout=VEFC

  ipType = BlythPozLegendre
  !!
  order=6
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

ipType=BlythPozLegendre, layout=VEFC

  ipType = BlythPozLegendre
  !!
  order=9
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

END PROGRAM main

܀ BlythPozChebyshev

In this example we test BlythPozChebyshev type triangle nodes.

PROGRAM main
use easifemBase
use easifemclasses
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order
integer(i4b), parameter :: elemType=Triangle
integer(i4b), parameter :: ipType=BlythPozChebyshev
CHARACTER(*), PARAMETER :: layout = "VEFC", fname="./results/test4"
real(dfp), parameter :: tol=1.0E-10
type(csvfile_) :: afile

Setup

  xij = zeros(2,3, 0.0_DFP)
  xij(1,:) = [0.0_DFP, 1.0_DFP, 0.0_DFP]
  xij(2,:) = [0.0_DFP, 0.0_DFP, 1.0_DFP]

  order=4
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(TRANSPOSE(ip)), &
    & "ipType=BlythPozChebyshev, layout=VEFC" // char_lf)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

See results

ipType=BlythPozChebyshev, layout=VEFC

0	0
1	0
0	1
0.14645	-2.43879E-14
0.5	-3.44909E-14
0.85355	-2.43879E-14
0.85355	0.14645
0.5	0.5
0.14645	0.85355
-2.43879E-14	0.85355
-3.44909E-14	0.5
-2.43879E-14	0.14645
0.21548	0.21548
0.56904	0.21548
0.21548	0.56904

  order=6
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

  order=9
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_BlythPozChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

END PROGRAM main

܀ IsaacLegendre

In this example we test IsaacLegendre type triangle nodes.

PROGRAM main
use easifemBase
use easifemclasses
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order
integer(i4b), parameter :: elemType=Triangle
integer(i4b), parameter :: ipType=IsaacLegendre
CHARACTER(*), PARAMETER :: layout = "VEFC", fname="./results/test5"
real(dfp), parameter :: tol=1.0E-10
type(csvfile_) :: afile

Setup

  xij = zeros(2,3, 0.0_DFP)
  xij(1,:) = [0.0_DFP, 1.0_DFP, 0.0_DFP]
  xij(2,:) = [0.0_DFP, 0.0_DFP, 1.0_DFP]

  order=4
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(TRANSPOSE(ip)), &
    & "ipType=IsaacLegendre, layout=VEFC" // char_lf)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

See results

ipType=IsaacLegendre, layout=VEFC

-8.32667E-17	-8.32667E-17
1	-8.32667E-17
-8.32667E-17	1
0.17267	-4.59259E-17
0.5	5.55112E-17
0.82733	-4.59259E-17
0.82733	0.17267
0.5	0.5
0.17267	0.82733
-4.59259E-17	0.82733
5.55112E-17	0.5
-4.59259E-17	0.17267
0.22216	0.22216
0.55569	0.22216
0.22216	0.55569

  !!
  order=6
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  !!
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

  !!
  order=9
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacLegendre_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

END PROGRAM main

܀ IsaacChebyshev

• In this example we test IsaacChebyshev type triangle nodes.

PROGRAM main
use easifemBase
use easifemclasses
implicit none
real(dfp), allocatable :: ip(:,:), xij(:, :)
integer(i4b) :: order
integer(i4b), parameter :: elemType=Triangle
integer(i4b), parameter :: ipType=IsaacChebyshev
CHARACTER(*), PARAMETER :: layout = "VEFC", fname="./results/test6"
real(dfp), parameter :: tol=1.0E-10
type(csvfile_) :: afile

setup

  xij = zeros(2,3, 0.0_DFP)
  xij(1,:) = [0.0_DFP, 1.0_DFP, 0.0_DFP]
  xij(2,:) = [0.0_DFP, 0.0_DFP, 1.0_DFP]

  order=4
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  call display(MdEncode(TRANSPOSE(ip)), &
    & "ipType=IsaacChebyshev, layout=VEFC" // char_lf)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

See results

ipType=IsaacChebyshev, layout=VEFC

0	0
1	0
0	1
0.14645	0
0.5	0
0.85355	0
0.85355	0.14645
0.5	0.5
0.14645	0.85355
0	0.85355
0	0.5
0	0.14645
0.20995	0.20995
0.58009	0.20995
0.20995	0.58009

  order=6
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

  order=9
  ip = InterpolationPoint(order=order, ipType=ipType, elemType=elemType, &
    & xij=xij, layout=layout)
  !!
  call afile%initiate( &
    & filename=fname//"_IsaacChebyshev_m="//tostring(order)//".csv", &
    & status="NEW", &
    & action="WRITE")
  call afile%open()
  call afile%write("x, y")
  call afile%write(ip, orient="TRANSPOSE")
  call afile%deallocate()

END PROGRAM main

↢