LagrangeCoeff

Returns the coefficients for lagrange polynomial.

This function returns the coefficient of basis functions.

$$p_{n}(x,y) = \sum_{a=0}^{p}\sum_{b=0}^{q}{c_{a,b}x^{a}y^{b}}$$

This routine returns the coefficients, $c_{a,b}$.

For example, for Lagrange polynomial of order 2, on equidistance grid, we have following coefficients of Monomials basis:

coeff:

Basis	$l_1$	$l_2$	$l_3$	$l_4$	$l_5$	$l_6$	$l_7$	$l_8$	$l_9$
$1$	0	0	0	0	0	0	0	0	1
$x$	0	0	0	0	0	0.5	0	-0.5	-0
$x^2$	0	0	-0	0	0	0.5	0	0.5	-1
$y$	0	0	0	0	-0.5	0	0.5	0	-0
$xy$	0.25	-0.25	0.25	-0.25	0	0	0	0	-0
$x^2 y$	-0.25	-0.25	0.25	0.25	0.5	-0	-0.5	0	-0
$y^2$	0	0	0	0	0.5	0	0.5	0	-1
$xy^2$	-0.25	0.25	0.25	-0.25	0	-0.5	0	0.5	-0
$x^2 y^2$	0.25	0.25	0.25	0.25	-0.5	-0.5	-0.5	-0.5	1

Interface 1

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Quadrangle(order, i, xij) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial
    INTEGER(I4B), INTENT(IN) :: i
    !! ith coefficients for lagrange polynomial
    REAL(DFP), INTENT(IN) :: xij(:, :)
    !! points in xij format, size(xij,2)
    REAL(DFP) :: ans(SIZE(xij, 2))
    !! coefficients
  END FUNCTION LagrangeCoeff_Quadrangle
END INTERFACE

order

Order of Lagrange polynomial.

️܀ See example

program main

  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :), coeff(:)
  integer(i4b), allocatable:: degree(:, :)
  CHARACTER( 20 ) :: layout
  order=1
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="VEFC"

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

  degree = LagrangeDegree_Quadrangle(order=order)
  coeff = LagrangeCoeff_Quadrangle(order=order, xij=x, i = 1)

  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff(1): ")

  coeff = LagrangeCoeff_Quadrangle(order=order, xij=x, i = 2)
  CALL Display(mdencode(coeff), "coeff(2): ")

  coeff = LagrangeCoeff_Quadrangle(order=order, xij=x, i = 3)
  CALL Display(mdencode(coeff), "coeff(3): ")

end program main

degree:

0	0
1	0
0	1
1	1

coeff(1):

0.25	-0.25	-0.25	0.25

coeff(2):

0.25	0.25	-0.25	-0.25

coeff(3):

0.25	0.25	0.25	0.25

This means:

$$L_{1}(x,y) = \frac{1}{4} - \frac{1}{4} \cdot x^1 - \frac{1}{4} \cdot y^1 + \frac{1}{4} xy$$ $$L_{2}(x,y) = \frac{1}{4} + \frac{1}{4} \cdot x^1 - \frac{1}{4} \cdot y^1 - \frac{1}{4} xy$$ $$L_{3}(x,y) = \frac{1}{4}+ \frac{1}{4} \cdot x^1 + \frac{1}{4} \cdot y^1 + \frac{1}{4} xy$$

↢

Interface 2

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Quadrangle(order, i, v, isVandermonde) &
    & RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial, it should be SIZE(v,2)-1
    INTEGER(I4B), INTENT(IN) :: i
    !! coefficient for ith lagrange polynomial
    REAL(DFP), INTENT(IN) :: v(:, :)
    !! vandermonde matrix size should be (order+1,order+1)
    LOGICAL(LGT), INTENT(IN) :: isVandermonde
    !! This is just to resolve interface issue
    REAL(DFP) :: ans(SIZE(v, 1))
    !! coefficients
  END FUNCTION LagrangeCoeff_Quadrangle
END INTERFACE

️܀ See example

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

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

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

  degree = LagrangeDegree_Quadrangle(order=order)

  V = LagrangeVandermonde(order=order, xij=x, elemType=Quadrangle)

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 1, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff(1): ")

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 2, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(coeff), "coeff(2): ")

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 3, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(coeff), "coeff(3): ")
end program main

degree:

0	0
1	0
0	1
1	1

coeff(1):

0.25	-0.25	-0.25	0.25

coeff(2):

0.25	0.25	-0.25	-0.25

coeff(3):

0.25	0.25	0.25	0.25

↢

Interface 3

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Quadrangle(order, i, v, ipiv) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial, it should be SIZE(x,2)-1
    INTEGER(I4B), INTENT(IN) :: i
    !! ith coefficients for lagrange polynomial
    REAL(DFP), INTENT(INOUT) :: v(:, :)
    !! LU decomposition of vandermonde matrix
    INTEGER(I4B), INTENT(IN) :: ipiv(:)
    !! inverse pivoting mapping, compes from LU decomposition
    REAL(DFP) :: ans(SIZE(v, 1))
    !! coefficients
  END FUNCTION LagrangeCoeff_Quadrangle
END INTERFACE

️܀ See example

program main
  use easifembase
  implicit none
  integer( i4b ) :: i1, i2, order, ipType
  real( dfp ), allocatable :: x(:,:), xij(:, :), coeff(:), V(:, :)
  integer(i4b), allocatable:: degree(:, :), ipiv(:)
  CHARACTER( 20 ) :: layout
  order=1
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="VEFC"

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

  degree = LagrangeDegree_Quadrangle(order=order)
  CALL Display(mdencode(degree), "degree: " // char_lf // char_lf )

  V = LagrangeVandermonde(order=order, xij=x, elemType=Quadrangle)

  CALL reallocate(ipiv, size(x, 2) )
  CALL GetLU(A=V, IPIV=ipiv)

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 1, V=V, ipiv=ipiv)
  CALL Display(mdencode(coeff), "coeff(1): "// char_lf // char_lf)

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 2, V=V, ipiv=ipiv)
  CALL Display(mdencode(coeff), "coeff(2): "// char_lf // char_lf)

  coeff = LagrangeCoeff_Quadrangle(order=order, i = 3, V=V, ipiv=ipiv)
  CALL Display(mdencode(coeff), "coeff(3): "// char_lf // char_lf)
end program main

degree:

0	0
1	0
0	1
1	1

coeff(1):

0.25	-0.25	-0.25	0.25

coeff(2):

0.25	0.25	-0.25	-0.25

coeff(3):

0.25	0.25	0.25	0.25

↢

Interface 4

܀ Interface

INTERFACE LagrangeCoeff_Quadrangle
  MODULE FUNCTION LagrangeCoeff_Quadrangle4(order, xij, basisType, &
    & refQuadrangle, alpha, beta, lambda) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: order
    !! order of polynomial
    REAL(DFP), INTENT(IN) :: xij(:, :)
    !! points in xij format, size(xij,2)
    INTEGER(I4B), OPTIONAL, INTENT(IN) :: basisType
    !! Monomials
    !! Jacobi
    !! Legendre
    !! Chebyshev
    !! Ultraspherical
    !! Heirarchical
    CHARACTER(*), OPTIONAL, INTENT(IN) :: refQuadrangle
    !! UNIT
    !! BIUNIT
    REAL(DFP), OPTIONAL, INTENT(IN) :: alpha
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: beta
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: lambda
    !! This parameter is needed when basisType is Ultraspherical
    REAL(DFP) :: ans(SIZE(xij, 2), SIZE(xij, 2))
    !! coefficients
  END FUNCTION LagrangeCoeff_Quadrangle4
END INTERFACE LagrangeCoeff_Quadrangle

order

order of polynomial

xij

interpolation points for Lagrange polynomials. The number of rows of xij should be 2 or 3. The number of columns of xij should be equal to the total number of degrees of freedom.

basisType

• Monomials
• Jacobi
• Legendre
• Chebyshev
• Ultraspherical
• Heirarchical

refQuadrangle

Reference Quadrangle can be UNIT or BIUNIT.

alpha, beta, lambda

• alpha and beta are needed when basisType is Jacobi
• lambda is needed when basisType is Ultraspherical

️܀ Example 1

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

  order=1
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout, &
    & xij=xij)

  degree = LagrangeDegree_Quadrangle(order=order)
  coeff = LagrangeCoeff_Quadrangle(  &
  & order=order, &
  & xij=x, &
  & basisType=Monomial, &
  & refQuadrangle="BIUNIT" )

  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff: ")
end program main

degree:

0	0
1	0
0	1
1	1

coeff:

Basis	$l_1$	$l_2$	$l_3$	$l_4$
$1$	0.25	0.25	0.25	0.25
$x$	-0.25	0.25	0.25	-0.25
$y$	-0.25	-0.25	0.25	0.25
$xy$	0.25	-0.25	0.25	-0.25

܀ Example 2

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

  order=2
  xij = RefQuadrangleCoord("BIUNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Quadrangle(  &
    & order=order, &
    & ipType=ipType,  &
    & layout=layout, &
    & xij=xij)

  degree = LagrangeDegree_Quadrangle(order=order)
  coeff = LagrangeCoeff_Quadrangle(  &
  & order=order, &
  & xij=x, &
  & basisType=Monomial, &
  & refQuadrangle="BIUNIT" )

  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff: ")
end program main

degree:

0	0
1	0
2	0
0	1
1	1
2	1
0	2
1	2
2	2

coeff

Basis	$l_1$	$l_2$	$l_3$	$l_4$	$l_5$	$l_6$	$l_7$	$l_8$	$l_9$
$1$	0	0	0	0	0	0	0	0	1
$x$	0	0	0	0	0	0.5	0	-0.5	-0
$x^2$	0	0	-0	0	0	0.5	0	0.5	-1
$y$	0	0	0	0	-0.5	0	0.5	0	-0
$xy$	0.25	-0.25	0.25	-0.25	0	0	0	0	-0
$x^2 y$	-0.25	-0.25	0.25	0.25	0.5	-0	-0.5	0	-0
$y^2$	0	0	0	0	0.5	0	0.5	0	-1
$xy^2$	-0.25	0.25	0.25	-0.25	0	-0.5	0	0.5	-0
$x^2 y^2$	0.25	0.25	0.25	0.25	-0.5	-0.5	-0.5	-0.5	1

↢

Interface 5

INTERFACE LagrangeCoeff_Quadrangle
  MODULE FUNCTION LagrangeCoeff_Quadrangle5(  &
    & p,  &
    & q,  &
    & xij, &
    & basisType1, &
    & basisType2, &
    & refQuadrangle,  &
    & alpha1, &
    & beta1,  &
    & lambda1, &
    & alpha2,  &
    & beta2,  &
    & lambda2) RESULT(ans)
    INTEGER(I4B), INTENT(IN) :: p
    !! order of polynomial in x direction
    INTEGER(I4B), INTENT(IN) :: q
    !! order of polynomial in y direction
    REAL(DFP), INTENT(IN) :: xij(:, :)
    !! points in xij format, size(xij,2)
    INTEGER(I4B), INTENT(IN) :: basisType1
    !! basisType in x direction
    !! Monomials
    !! Jacobi
    !! Legendre
    !! Chebyshev
    !! Ultraspherical
    !! Heirarchical
    INTEGER(I4B), INTENT(IN) :: basisType2
    !! basisType in y direction
    !! Monomials
    !! Jacobi
    !! Legendre
    !! Chebyshev
    !! Ultraspherical
    !! Heirarchical
    CHARACTER(*), OPTIONAL, INTENT(IN) :: refQuadrangle
    !! UNIT
    !! BIUNIT
    REAL(DFP), OPTIONAL, INTENT(IN) :: alpha1
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: beta1
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: lambda1
    !! This parameter is needed when basisType is Ultraspherical
    REAL(DFP), OPTIONAL, INTENT(IN) :: alpha2
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: beta2
    !! This parameter is needed when basisType is Jacobi
    REAL(DFP), OPTIONAL, INTENT(IN) :: lambda2
    !! This parameter is needed when basisType is Ultraspherical
    REAL(DFP) :: ans(SIZE(xij, 2), SIZE(xij, 2))
    !! coefficients
  END FUNCTION LagrangeCoeff_Quadrangle5
END INTERFACE LagrangeCoeff_Quadrangle

p, q, r

order of polynomial in x, y and z direction.

xij

interpolation points for Lagrange polynomials. The number of rows of xij should be 2 or 3. The number of columns of xij should be equal to the total number of degrees of freedom.

basisType1, basisType2

Basis type in x and y direction. They can take following values:

• Monomials
• Jacobi
• Legendre
• Chebyshev
• Ultraspherical
• Heirarchical

refQuadrangle

Reference Quadrangle can be UNIT or BIUNIT.

alpha1, beta1, lambda1

• alpha1 and beta1 are needed when basisType1 is Jacobi
• lambda1 is needed when basisType1 is Ultraspherical

alpha2, beta2, lambda2

• alpha2 and beta2 are needed when basisType2 is Jacobi
• lambda2 is needed when basisType2 is Ultraspherical