LagrangeCoeff

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}$.

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.

️܀ 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
    !! 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_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(:, :)
  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 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_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_Triangle4(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_Triangle4
END INTERFACE

️܀ 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

↢