LagrangeCoeff_Triangle

Returns the coefficients for Lagrange polynomial.

This function returns the coefficient of basis functions.

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

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

note

There are four interfaces for this function.

Interface 1

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Triangle(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_Triangle
END INTERFACE

order

Order of Lagrange polynomial.

i

ith Lagrnage polynomial.

xij

Points in xij format, size(xij,2) denotes the number of points. These are actually the interpolation points.

note

• This subroutine will use the monomial basis functions.
• Internally we compute the Vandermonde matrix and then solve it to get the coefficients.

️܀ 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 = RefTriangleCoord("UNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Triangle( order=order, &
    & ipType=ipType, layout=layout, xij=xij)
  degree = LagrangeDegree_Triangle(order=order)
  coeff = LagrangeCoeff_Triangle(order=order, xij=x, i = 1)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff(1): ")
  coeff = LagrangeCoeff_Triangle(order=order, xij=x, i = 2)
  CALL Display(mdencode(coeff), "coeff(2): ")
  coeff = LagrangeCoeff_Triangle(order=order, xij=x, i = 3)
  CALL Display(mdencode(coeff), "coeff(3): ")
end program main

degree:

x	y
0	0
1	0
0	1

coeff(1):

1	-1	-1

coeff(2):

0	1	0

coeff(3):

0	0	1

This means:

$$L_{1}(x,y) = 1\cdot x^0 y^0 - 1.0 \cdot x^1 y^0 - 1.0 \cdot x^0 y^1 = 1-x-y$$ $$L_{2}(x,y) = 0\cdot x^0 y^0 + 1.0 \cdot x^1 y^0 + 0.0 \cdot x^0 y^1 = x$$ $$L_{3}(x,y) = 0\cdot x^0 y^0 + 0.0 \cdot x^1 y^0 + 1.0 \cdot x^0 y^1 = y$$

↢

Interface 2

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Triangle(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
    !! 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 of ith Lagrange polynomial
  END FUNCTION LagrangeCoeff_Triangle
END INTERFACE

order

order of Lagrange polynomial, it should be equal to the SIZE(v,2)-1.

i

ith Lagrange polynomial.

v

Vandermonde matrix, The size should be at least (order+1,order+1).

️܀ 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 = RefTriangleCoord("UNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Triangle( order=order, &
    & ipType=ipType, layout=layout, xij=xij)
  degree = LagrangeDegree_Triangle(order=order)
  V = LagrangeVandermonde(order=order, xij=x, elemType=Triangle)
  coeff = LagrangeCoeff_Triangle(order=order, i = 1, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff(1): ")
  coeff = LagrangeCoeff_Triangle(order=order, i = 2, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(coeff), "coeff(2): ")
  coeff = LagrangeCoeff_Triangle(order=order, i = 3, V=V, isVanderMonde=.true.)
  CALL Display(mdencode(coeff), "coeff(3): ")
end program main

degree:

x	y
0	0
1	0
0	1

coeff(1):

1	-1	-1

coeff(2):

0	1	0

coeff(3):

0	0	1

This means:

$$L_{1}(x,y) = 1\cdot x^0 y^0 - 1.0 \cdot x^1 y^0 - 1.0 \cdot x^0 y^1 = 1-x-y$$ $$L_{2}(x,y) = 0\cdot x^0 y^0 + 1.0 \cdot x^1 y^0 + 0.0 \cdot x^0 y^1 = x$$ $$L_{3}(x,y) = 0\cdot x^0 y^0 + 0.0 \cdot x^1 y^0 + 1.0 \cdot x^0 y^1 = y$$

↢

Interface 3

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Triangle(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 lagrange polynomial
    REAL(DFP), INTENT(INOUT) :: v(:, :)
    !! LU decomposition of vandermonde matrix
    INTEGER(I4B), INTENT(IN) :: ipiv(:)
    !! inverse pivoting mapping, it is returned from LU decomposition
    REAL(DFP) :: ans(SIZE(v, 1))
    !! coefficients of ith Lagrange polynomial
  END FUNCTION LagrangeCoeff_Triangle
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 = RefTriangleCoord("UNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Triangle( order=order, &
    & ipType=ipType, layout=layout, xij=xij)
  degree = LagrangeDegree_Triangle(order=order)
  V = LagrangeVandermonde(order=order, xij=x, elemType=Triangle)
  call reallocate(ipiv, size(x, 2) )
  CALL GetLU(A=V, IPIV=ipiv)
  coeff = LagrangeCoeff_Triangle(order=order, i = 1, V=V, ipiv=ipiv)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff(1): ")
  coeff = LagrangeCoeff_Triangle(order=order, i = 2, V=V, ipiv=ipiv)
  CALL Display(mdencode(coeff), "coeff(2): ")
  coeff = LagrangeCoeff_Triangle(order=order, i = 3, V=V, ipiv=ipiv)
  CALL Display(mdencode(coeff), "coeff(3): ")
end program main

degree:

x	y
0	0
1	0
0	1

coeff(1):

1	-1	-1

coeff(2):

0	1	0

coeff(3):

0	0	1

This means:

$$L_{1}(x,y) = 1\cdot x^0 y^0 - 1.0 \cdot x^1 y^0 - 1.0 \cdot x^0 y^1 = 1-x-y$$ $$L_{2}(x,y) = 0\cdot x^0 y^0 + 1.0 \cdot x^1 y^0 + 0.0 \cdot x^0 y^1 = x$$ $$L_{3}(x,y) = 0\cdot x^0 y^0 + 0.0 \cdot x^1 y^0 + 1.0 \cdot x^0 y^1 = y$$

↢

Interface 4

܀ Interface

INTERFACE
  MODULE FUNCTION LagrangeCoeff_Triangle(order, xij, basisType, refTriangle) 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
    !! Heirarchical
    CHARACTER(*), OPTIONAL, INTENT(IN) :: refTriangle
    REAL(DFP) :: ans(SIZE(xij, 2), SIZE(xij, 2))
    !! coefficients
  END FUNCTION LagrangeCoeff_Triangle
END INTERFACE

This is the master function, which allows us to specify the basis type and reference triangle.

basisType

• Here, basis type is the type of basis functions we want to use for contructing the Lagrange polynomials. It can take following values:
  • Monomials (default)
  • Jacobi, Orthogonal, Legendre, Lobatto, Ultraspherical (in this case we call Dubiner functions)
  • Heirarchical

ans

It represents the coefficients of the Lagrange polynomial. ans(:, i) denotes the coefficients of ith 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 = RefTriangleCoord("UNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Triangle( order=order, &
    & ipType=ipType, layout=layout, xij=xij)
  degree = LagrangeDegree_Triangle(order=order)
  coeff = LagrangeCoeff_Triangle(order=order, xij=x)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff: ")
end program main

degree:

x	y
0	0
1	0
0	1

coeff:

L1	L2	L3
1	-0	-0
-1	1	-0
-1	0	1

This means:

$$L_{1}(x,y) = 1\cdot x^0 y^0 - 1.0 \cdot x^1 y^0 - 1.0 \cdot x^0 y^1 = 1-x-y$$ $$L_{2}(x,y) = 0\cdot x^0 y^0 + 1.0 \cdot x^1 y^0 + 0.0 \cdot x^0 y^1 = x$$ $$L_{3}(x,y) = 0\cdot x^0 y^0 + 0.0 \cdot x^1 y^0 + 1.0 \cdot x^0 y^1 = y$$

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 = RefTriangleCoord("UNIT")
  ipType=Equidistance
  layout="VEFC"
  x =InterpolationPoint_Triangle( order=order, &
    & ipType=ipType, layout=layout, xij=xij)
  degree = LagrangeDegree_Triangle(order=order)
  coeff = LagrangeCoeff_Triangle(order=order, xij=x)
  CALL Display(mdencode(degree), "degree: ")
  CALL Display(mdencode(coeff), "coeff: ")
end program main

degree:

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

coeff:

1	-0	-0	0	-0	0
-3	-1	-0	4	-0	0
2	2	0	-4	-0	0
-3	0	-1	0	-0	4
4	0	0	-4	4	-4
2	0	2	0	0	-4

↢