InterpolationPoint

This routine returns the interplation points on quadrangle.

Interface

܀ Interface

INTERFACE InterpolationPoint_Quadrangle
  MODULE FUNCTION InterpolationPoint_Quadrangle1(order, ipType, xij, &
    & layout) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of element
    INTEGER(I4B), INTENT(IN) :: ipType
    !! interpolation point type
    REAL(DFP), OPTIONAL, INTENT(IN) :: xij(:, :)
    !! four vertices of quadrangle in xij format
    CHARACTER(*), INTENT(IN) :: layout
    !! VEFC, INCREASING
    REAL(DFP), ALLOCATABLE :: ans(:, :)
    !! interpolation points in xij format
  END FUNCTION InterpolationPoint_Quadrangle1
END INTERFACE InterpolationPoint_Quadrangle

xij

• xij contains nodal coordinates of quadrangle in xij format.
• The number of rows in xij can be 2 or 3 (for 2D or 3D)
• The number of columns in xij is 4
• If xij is absent then biunit quadrangle is assumed.

ipType

• ipType is interpolation point type, it can take following values

• Equidistance, uniformly/evenly distributed points

• GaussLegendreLobatto

• GaussChebyshev1Lobatto

• GaussJacobiLobatto

• GaussUltrasphericalLobatto

layout

layout specifies the arrangement of interpolation points.

It has two options:

VEFC

In "VEFC" format (vertex, edge, face, cell), the first four entries denote the vertex points, then we have edges, and then internal nodes.

The internal nodes also follow the same convention. Please read Gmsh manual on this topic.

INCREASING

In "INCREASING" format the interpolation points are stored column-wise, that is first points along a y-axis are stored first.

For example for order = 3:

VEFC

$x$	-1	1	1	-1	-0.33333	0.33333	1	1	0.33333	-0.33333	-1	-1	-0.33333	-0.33333	0.333 33	0.33333
$y$	-1	-1	1	1	-1	-1	-0.33333	0.33333	1	1	0.33333	-0.33333	-0.33333	0.33333	0.3333 3	-0.33333

INCREASING

$x$	-1	-1	-1	-1	-0.33333	-0.33333	-0.33333	-0.33333	0.33333	0.33333	0.33333	0.33333	1	1	1	1
$y$		-1	-0.33333	0.33333	1	-1	-0.33333	0.33333	1	-1	-0.33333	0.33333	1	-1	-0.33333	0.33333

️܀ Equidistance

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :)
  CHARACTER( 20 ) :: layout

  order=3
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="VEFC"

  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout,  &
    & xij=xij)

  call display( Mdencode(x), "xij (order="//tostring(order)//")=" )


  order=3
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="INCREASING"

  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout,  &
    & xij=xij)

  call display( Mdencode(x), "xij (order="//tostring(order)//")=" )

end program main

See results

xij (order=3), VEFC:

$x$	-1	1	1	-1	-0.33333	0.33333	1	1	0.33333	-0.33333	-1	-1	-0.33333	-0.33333	0.333 33	0.33333
$y$	-1	-1	1	1	-1	-1	-0.33333	0.33333	1	1	0.33333	-0.33333	-0.33333	0.33333	0.3333 3	-0.33333

xij (order=3)=, INCREASING

$x$	-1	-1	-1	-1	-0.33333	-0.33333	-0.33333	-0.33333	0.33333	0.33333	0.33333	0.33333	1	1	1	1
$y$		-1	-0.33333	0.33333	1	-1	-0.33333	0.33333	1	-1	-0.33333	0.33333	1	-1	-0.33333	0.33333

GaussLegendreLobatto

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :)
  CHARACTER( 20 ) :: layout

  order=3
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=GaussLegendreLobatto
  layout="VEFC"

  x =InterpolationPoint_Quadrangle( &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout,  &
    & xij=xij)

  call display( Mdencode(TRANSPOSE(x)), "xij (order="//tostring(order)//")=" )
  call blanklines(nol=2)
end program main

See results

xij (order=5)=

-1	-1
1	-1
1	1
-1	1
-0.44721	-1
0.44721	-1
1	-0.44721
1	0.44721
0.44721	1
-0.44721	1
-1	0.44721
-1	-0.44721
-0.44721	-0.44721
-0.44721	0.44721
0.44721	0.44721
0.44721	-0.44721

GaussChebyshevLobatto

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :)
  CHARACTER( 20 ) :: layout

  order=3
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=GaussChebyshevLobatto
  layout="VEFC"

  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout,  &
    & xij=xij)

  call display( Mdencode(TRANSPOSE(x)), "xij (order="//tostring(order)//")=" )
  call blanklines(nol=2)
end program main

See results

xij (order=3) =

0	0
1	0
0	1
0.11747	0
0.35738	0
0.64262	0
0.88253	0
0.88253	0.11747
0.64262	0.35738
0.35738	0.64262
0.11747	0.88253
0	0.88253
0	0.64262
0	0.35738
0	0.11747
0.15829	0.15829
0.68343	0.15829
0.15829	0.68343
0.4133	0.17339
0.4133	0.4133
0.17339	0.4133

↢

GaussJacobiLobatto

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :)
  CHARACTER( 20 ) :: layout
  REAL( DFP )  :: alpha, beta

  order=3
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=GaussJacobiLobatto
  layout="VEFC"
  alpha = 0.0_DFP
  beta = 0.0_DFP

  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType, &
    & layout=layout, &
    & xij=xij, &
    & alpha=alpha, &
    & beta = beta)

  call display( Mdencode(TRANSPOSE(x)), "xij (order="//tostring(order)//")=" )
  call blanklines(nol=2)
end program main

See results

xij (order=3)=

-1	-1
1	-1
1	1
-1	1
-0.44721	-1
0.44721	-1
1	-0.44721
1	0.44721
0.44721	1
-0.44721	1
-1	0.44721
-1	-0.44721
-0.44721	-0.44721
-0.44721	0.44721
0.44721	0.44721
0.44721	-0.44721

GaussUltrasphericalLobatto

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :)
  CHARACTER( 20 ) :: layout
  real( dfp ) :: lambda

  order=3
  lambda = 0.5_DFP
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=GaussUltrasphericalLobatto
  layout="VEFC"
  x =InterpolationPoint_Quadrangle( &
    & order=order, &
    & ipType=ipType, &
    & layout=layout, &
    & xij=xij, &
    & lambda=lambda)

  call display( Mdencode(TRANSPOSE(x)), "xij (order="//tostring(order)//")=" )
  call blanklines(nol=2)
end program main

See results

xij (order=3)=

-1	-1
1	-1
1	1
-1	1
-0.44721	-1
0.44721	-1
1	-0.44721
1	0.44721
0.44721	1
-0.44721	1
-1	0.44721
-1	-0.44721
-0.44721	-0.44721
-0.44721	0.44721
0.44721	0.44721
0.44721	-0.44721

↢

Interface 2

Interface

INTERFACE InterpolationPoint_Quadrangle
  MODULE FUNCTION InterpolationPoint_Quadrangle2(  &
    & p, q, ipType1, ipType2, layout, xij, alpha1, beta1, &
    & lambda1, alpha2, beta2, lambda2) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: p
    !! order of element in x direction
    INTEGER(I4B), INTENT(IN) :: q
    !! order of element in y direction
    INTEGER(I4B), INTENT(IN) :: ipType1
    !! interpolation point type in x direction
    INTEGER(I4B), INTENT(IN) :: ipType2
    !! interpolation point type in y direction
    CHARACTER(*), INTENT(IN) :: layout
    !! VEFC, INCREASING
    REAL(DFP), OPTIONAL, INTENT(IN) :: xij(:, :)
    !! four vertices of quadrangle in xij format
    REAL(DFP), OPTIONAL, INTENT(IN) :: alpha1
    !! Jacobi parameter
    REAL(DFP), OPTIONAL, INTENT(IN) :: beta1
    !! Jacobi parameter
    REAL(DFP), OPTIONAL, INTENT(IN) :: lambda1
    !! Ultraspherical parameter
    REAL(DFP), OPTIONAL, INTENT(IN) :: alpha2
    !! Jacobi parameter
    REAL(DFP), OPTIONAL, INTENT(IN) :: beta2
    !! Jacobi parameter
    REAL(DFP), OPTIONAL, INTENT(IN) :: lambda2
    !! Ultraspherical parameter
    REAL(DFP), ALLOCATABLE :: ans(:, :)
    !! interpolation points in xij format
  END FUNCTION InterpolationPoint_Quadrangle2
END INTERFACE InterpolationPoint_Quadrangle