Lagrange_t Derived Type

  1. SELF_Lagrange_t.f90
  2. SELF_Lagrange_t
  3. Lagrange_t

type, public :: Lagrange_t

A data structure for working with Lagrange Interpolating Polynomials in one, two, and three dimensions. The Lagrange data-structure stores the information necessary to interpolate between two sets of grid-points and to estimate the derivative of data at native grid points. Routines for multidimensional interpolation are based on the tensor product of 1-D interpolants. It is assumed that the polynomial degree (and the interpolation nodes) are the same in each direction. This assumption permits the storage of only one array of interpolation nodes and barycentric weights and is what allows this data structure to be flexible.


Inherited by

type~~lagrange_t~~InheritedByGraph type~lagrange_t Lagrange_t type~lagrange Lagrange type~lagrange->type~lagrange_t type~lagrange~2 Lagrange type~lagrange~2->type~lagrange_t type~lagrange~3 Lagrange type~lagrange~3->type~lagrange_t type~self_dataobj SELF_DataObj type~self_dataobj->type~lagrange interp type~twopointvector3d_t TwoPointVector3D_t type~twopointvector3d_t->type~self_dataobj type~vector3d_t Vector3D_t type~vector3d_t->type~self_dataobj type~scalar2d_t Scalar2D_t type~scalar2d_t->type~self_dataobj type~twopointvector2d_t TwoPointVector2D_t type~twopointvector2d_t->type~self_dataobj type~vector2d_t Vector2D_t type~vector2d_t->type~self_dataobj type~scalar1d_t Scalar1D_t type~scalar1d_t->type~self_dataobj type~tensor2d_t Tensor2D_t type~tensor2d_t->type~self_dataobj type~scalar3d_t Scalar3D_t type~scalar3d_t->type~self_dataobj type~tensor3d_t Tensor3D_t type~tensor3d_t->type~self_dataobj type~scalar2d~2 Scalar2D type~scalar2d~2->type~scalar2d_t type~scalar1d~2 Scalar1D type~scalar1d~2->type~scalar1d_t type~twopointvector2d TwoPointVector2D type~twopointvector2d->type~twopointvector2d_t type~vector2d~2 Vector2D type~vector2d~2->type~vector2d_t type~tensor3d~2 Tensor3D type~tensor3d~2->type~tensor3d_t type~vector3d Vector3D type~vector3d->type~vector3d_t type~vector2d Vector2D type~vector2d->type~vector2d_t type~scalar1d Scalar1D type~scalar1d->type~scalar1d_t type~scalar1d~3 Scalar1D type~scalar1d~3->type~scalar1d_t type~tensor2d~2 Tensor2D type~tensor2d~2->type~tensor2d_t type~twopointvector3d TwoPointVector3D type~twopointvector3d->type~twopointvector3d_t type~twopointvector2d~2 TwoPointVector2D type~twopointvector2d~2->type~twopointvector2d_t type~scalar3d Scalar3D type~scalar3d->type~scalar3d_t type~tensor3d Tensor3D type~tensor3d->type~tensor3d_t type~tensor2d Tensor2D type~tensor2d->type~tensor2d_t type~scalar2d Scalar2D type~scalar2d->type~scalar2d_t type~scalar3d~2 Scalar3D type~scalar3d~2->type~scalar3d_t type~vector3d~2 Vector3D type~vector3d~2->type~vector3d_t type~twopointvector3d~2 TwoPointVector3D type~twopointvector3d~2->type~twopointvector3d_t type~semhex SEMHex type~semhex->type~vector3d x, nHat type~semhex->type~scalar3d nScale, J type~semhex->type~tensor3d dxds, dsdx type~mappedvector3d_t MappedVector3D_t type~mappedvector3d_t->type~vector3d type~mappedvector3d_t->type~semhex geometry type~geometry1d Geometry1D type~geometry1d->type~scalar1d x, dxds type~mappedscalar3d_t MappedScalar3D_t type~mappedscalar3d_t->type~scalar3d type~mappedscalar3d_t->type~semhex geometry type~mappedscalar2d_t MappedScalar2D_t type~mappedscalar2d_t->type~scalar2d type~semquad SEMQuad type~mappedscalar2d_t->type~semquad geometry type~mappedtwopointvector2d_t MappedTwoPointVector2D_t type~mappedtwopointvector2d_t->type~twopointvector2d type~mappedtwopointvector2d_t->type~semquad geometry type~mappedvector2d_t MappedVector2D_t type~mappedvector2d_t->type~vector2d type~mappedvector2d_t->type~semquad geometry type~semquad->type~vector2d x, nHat type~semquad->type~tensor2d dxds, dsdx type~semquad->type~scalar2d nScale, J type~mappedscalar1d_t MappedScalar1D_t type~mappedscalar1d_t->type~scalar1d type~mappedscalar1d_t->type~geometry1d geometry type~mappedtwopointvector3d_t MappedTwoPointVector3D_t type~mappedtwopointvector3d_t->type~twopointvector3d type~mappedtwopointvector3d_t->type~semhex geometry type~dgmodel1d_t DGModel1D_t type~dgmodel1d_t->type~geometry1d geometry type~mappedscalar1d MappedScalar1D type~dgmodel1d_t->type~mappedscalar1d solution, solutionGradient, flux, source, fluxDivergence, dSdt, workSol type~mappedvector3d MappedVector3D type~mappedvector3d->type~mappedvector3d_t type~mappedscalar3d MappedScalar3D type~mappedscalar3d->type~mappedscalar3d_t type~mappedscalar3d~2 MappedScalar3D type~mappedscalar3d~2->type~mappedscalar3d_t type~mappedtwopointvector2d~2 MappedTwoPointVector2D type~mappedtwopointvector2d~2->type~mappedtwopointvector2d_t type~mappedvector2d MappedVector2D type~mappedvector2d->type~mappedvector2d_t type~mappedvector2d~2 MappedVector2D type~mappedvector2d~2->type~mappedvector2d_t type~mappedscalar1d~2 MappedScalar1D type~mappedscalar1d~2->type~mappedscalar1d_t type~dgmodel3d_t DGModel3D_t type~dgmodel3d_t->type~semhex geometry type~dgmodel3d_t->type~mappedvector3d solutionGradient, flux type~dgmodel3d_t->type~mappedscalar3d solution, source, fluxDivergence, dSdt, workSol type~mappedvector3d~2 MappedVector3D type~mappedvector3d~2->type~mappedvector3d_t type~mappedscalar2d MappedScalar2D type~mappedscalar2d->type~mappedscalar2d_t type~mappedscalar2d~2 MappedScalar2D type~mappedscalar2d~2->type~mappedscalar2d_t type~mappedtwopointvector2d MappedTwoPointVector2D type~mappedtwopointvector2d->type~mappedtwopointvector2d_t type~dgmodel2d_t DGModel2D_t type~dgmodel2d_t->type~semquad geometry type~dgmodel2d_t->type~mappedvector2d solutionGradient, flux type~dgmodel2d_t->type~mappedscalar2d solution, source, fluxDivergence, dSdt, workSol type~mappedscalar1d->type~mappedscalar1d_t type~mappedtwopointvector3d MappedTwoPointVector3D type~mappedtwopointvector3d->type~mappedtwopointvector3d_t type~mappedtwopointvector3d~2 MappedTwoPointVector3D type~mappedtwopointvector3d~2->type~mappedtwopointvector3d_t type~dgmodel1d DGModel1D type~dgmodel1d->type~dgmodel1d_t type~dgmodel1d~2 DGModel1D type~dgmodel1d~2->type~dgmodel1d_t type~ecdgmodel2d_t ECDGModel2D_t type~ecdgmodel2d_t->type~mappedtwopointvector2d twoPointFlux type~dgmodel2d DGModel2D type~ecdgmodel2d_t->type~dgmodel2d type~dgmodel2d->type~dgmodel2d_t type~dgmodel2d~2 DGModel2D type~dgmodel2d~2->type~dgmodel2d_t type~dgmodel3d DGModel3D type~dgmodel3d->type~dgmodel3d_t type~esatmo3d_t ESAtmo3D_t type~esatmo3d_t->type~mappedvector3d diffFlux type~esatmo3d_t->type~mappedscalar3d diffDiv type~ecdgmodel3d ECDGModel3D type~esatmo3d_t->type~ecdgmodel3d type~esatmo2d_t ESAtmo2D_t type~esatmo2d_t->type~mappedvector2d diffFlux type~esatmo2d_t->type~mappedscalar2d diffDiv type~ecdgmodel2d ECDGModel2D type~esatmo2d_t->type~ecdgmodel2d type~dgmodel3d~2 DGModel3D type~dgmodel3d~2->type~dgmodel3d_t type~linearshallowwater2d_t LinearShallowWater2D_t type~linearshallowwater2d_t->type~mappedscalar2d fCori type~linearshallowwater2d_t->type~dgmodel2d type~lineareuler2d_pml_t LinearEuler2D_PML_t type~lineareuler2d_pml_t->type~mappedscalar2d sigma_x, sigma_y type~lineareuler2d_t LinearEuler2D_t type~lineareuler2d_pml_t->type~lineareuler2d_t type~ecdgmodel3d_t ECDGModel3D_t type~ecdgmodel3d_t->type~mappedtwopointvector3d twoPointFlux type~ecdgmodel3d_t->type~dgmodel3d type~advection_diffusion_1d_t advection_diffusion_1d_t type~advection_diffusion_1d_t->type~dgmodel1d type~ecdgmodel2d~2 ECDGModel2D type~ecdgmodel2d~2->type~ecdgmodel2d_t type~nulldgmodel3d_t NullDGModel3D_t type~nulldgmodel3d_t->type~dgmodel3d type~esatmo2d~2 ESAtmo2D type~esatmo2d~2->type~esatmo2d_t type~lineareuler2d_pml LinearEuler2D_PML type~lineareuler2d_pml->type~lineareuler2d_pml_t type~ecdgmodel2d->type~ecdgmodel2d_t type~nulldgmodel1d_t NullDGModel1D_t type~nulldgmodel1d_t->type~dgmodel1d type~lineareuler3d_t LinearEuler3D_t type~lineareuler3d_t->type~dgmodel3d type~lineareuler2d_t->type~dgmodel2d type~lineareuler2d_pml~2 LinearEuler2D_PML type~lineareuler2d_pml~2->type~lineareuler2d_pml_t type~ecdgmodel3d~2 ECDGModel3D type~ecdgmodel3d~2->type~ecdgmodel3d_t type~esatmo3d~2 ESAtmo3D type~esatmo3d~2->type~esatmo3d_t type~linearshallowwater2d LinearShallowWater2D type~linearshallowwater2d->type~linearshallowwater2d_t type~esatmo3d ESAtmo3D type~esatmo3d->type~esatmo3d_t type~nulldgmodel2d_t NullDGModel2D_t type~nulldgmodel2d_t->type~dgmodel2d type~advection_diffusion_3d_t advection_diffusion_3d_t type~advection_diffusion_3d_t->type~dgmodel3d type~linearshallowwater2d~2 LinearShallowWater2D type~linearshallowwater2d~2->type~linearshallowwater2d_t type~burgers1d_t Burgers1D_t type~burgers1d_t->type~dgmodel1d type~advection_diffusion_2d_t advection_diffusion_2d_t type~advection_diffusion_2d_t->type~dgmodel2d type~esatmo2d ESAtmo2D type~esatmo2d->type~esatmo2d_t type~ecdgmodel3d->type~ecdgmodel3d_t type~lineareuler3d LinearEuler3D type~lineareuler3d->type~lineareuler3d_t type~lineareuler2d LinearEuler2D type~lineareuler2d->type~lineareuler2d_t type~nulldgmodel1d NullDGModel1D type~nulldgmodel1d->type~nulldgmodel1d_t type~burgers1d Burgers1D type~burgers1d->type~burgers1d_t type~lineareuler3d~2 LinearEuler3D type~lineareuler3d~2->type~lineareuler3d_t type~advection_diffusion_3d~2 advection_diffusion_3d type~advection_diffusion_3d~2->type~advection_diffusion_3d_t type~ecadvection3d_t ECAdvection3D_t type~ecadvection3d_t->type~ecdgmodel3d type~advection_diffusion_1d advection_diffusion_1d type~advection_diffusion_1d->type~advection_diffusion_1d_t type~nulldgmodel3d~2 NullDGModel3D type~nulldgmodel3d~2->type~nulldgmodel3d_t type~lineareuler2d~2 LinearEuler2D type~lineareuler2d~2->type~lineareuler2d_t type~nulldgmodel2d~2 NullDGModel2D type~nulldgmodel2d~2->type~nulldgmodel2d_t type~burgers1d~2 Burgers1D type~burgers1d~2->type~burgers1d_t type~ecadvection2d_t ECAdvection2D_t type~ecadvection2d_t->type~ecdgmodel2d type~nulldgmodel3d NullDGModel3D type~nulldgmodel3d->type~nulldgmodel3d_t type~nulldgmodel2d NullDGModel2D type~nulldgmodel2d->type~nulldgmodel2d_t type~advection_diffusion_3d advection_diffusion_3d type~advection_diffusion_3d->type~advection_diffusion_3d_t type~advection_diffusion_1d~2 advection_diffusion_1d type~advection_diffusion_1d~2->type~advection_diffusion_1d_t type~nulldgmodel1d~2 NullDGModel1D type~nulldgmodel1d~2->type~nulldgmodel1d_t type~advection_diffusion_2d advection_diffusion_2d type~advection_diffusion_2d->type~advection_diffusion_2d_t type~advection_diffusion_2d~2 advection_diffusion_2d type~advection_diffusion_2d~2->type~advection_diffusion_2d_t type~ecadvection3d ECAdvection3D type~ecadvection3d->type~ecadvection3d_t type~ecadvection2d~2 ECAdvection2D type~ecadvection2d~2->type~ecadvection2d_t type~ecadvection2d ECAdvection2D type~ecadvection2d->type~ecadvection2d_t type~ecadvection3d~2 ECAdvection3D type~ecadvection3d~2->type~ecadvection3d_t

Contents

Variables

Type-Bound Procedures

Source Code

Lagrange_t

Components

TypeVisibilityAttributesNameInitial
integer, public :: M

The number of target points.
integer, public :: N

The number of control points.
real(kind=prec), public, pointer, contiguous, dimension(:,:):: bMatrix

The boundary interpolation matrix that is used to map a grid of nodal values at the control points to the element boundaries.
real(kind=prec), public, pointer, contiguous, dimension(:):: bWeights

The barycentric weights that are calculated from the controlPoints and used for interpolation.
integer, public :: controlNodeType
real(kind=prec), public, pointer, contiguous, dimension(:):: controlPoints

The set of nodes in one dimension where data is known. To create higher dimension interpolation and differentiation operators, structured grids in two and three dimensions are created by tensor products of the controlPoints. This design decision implies that all spectral element methods supported by the Lagrange class have the same polynomial degree in each computational/spatial dimension. In practice, the controlPoints are the Legendre-Gauss, Legendre-Gauss-Lobatto, Legendre-Gauss-Radau, Chebyshev-Gauss, Chebyshev-Gauss-Lobatto, or Chebyshev-Gauss-Radau quadrature points over the domain [-1,1] (computational space). The Init routine for this class restricts controlPoints to one of these quadrature types or uniform points on [-1,1].
real(kind=prec), public, pointer, contiguous, dimension(:,:):: dMatrix

The derivative matrix for mapping function nodal values to a nodal values of the derivative estimate. The dMatrix is based on a strong form of the derivative.
real(kind=prec), public, pointer, contiguous, dimension(:,:):: dSplitMatrix

The split-form derivative matrix D_split = D - 0.5M^{-1}B, where B is the SBP boundary operator.

D_split is skew-symmetric under the M inner product: MD_split + D_split^TM = 0. Unlike D, it is NOT an SBP operator (the weighted symmetric part is zero, not the boundary term B). This property makes D_split ideal for the EC-DGSEM split-form volume integral:

du/dt = -(2/J) sum_n D_split[n,i] * F_EC(u_i, u_n) + (1/J) M^{-1} B^T f_Riemann

The skew-symmetry of the volume term guarantees it contributes zero to the entropy rate, so all entropy change passes through the surface term. An entropy-dissipative Riemann solver then makes the full scheme entropy-stable.

Using D_split in the volume combined with plain f_Riemann on the surface is algebraically identical to using D in the volume with the penalty (f_Riemann - f_local) on the surface (Trixi.jl convention).

In SELF index convention (dSplitMatrix(ii,i) = D_split[i-1, ii-1]): dSplitMatrix(ii,i) = dMatrix(ii,i) - 0.5(bMatrix(i,2)bMatrix(ii,2) - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i)
real(kind=prec), public, pointer, contiguous, dimension(:,:):: dgMatrix

The derivative matrix for mapping function nodal values to a nodal values of the derivative estimate. The dgMatrix is based on a weak form of the derivative. It must be used with bMatrix to account for boundary contributions in the weak form.
real(kind=prec), public, pointer, contiguous, dimension(:,:):: iMatrix

The interpolation matrix (transpose) for mapping data from the control grid to the target grid.
real(kind=prec), public, pointer, contiguous, dimension(:,:,:):: mortarP

Mortar L2-projection matrices, the adjoints of mortarR under the exact L2 inner product on [-1,1]: given traces g_k(s) on the two sub-edges, the projected big-side trace is u(i) = sum_k sum_ii mortarP(ii,i,k)*g_k(ii).

The matrices are built as P_k = M^{-1} B_k with the exact (dense) 1-D mass matrix M(m,i) = int l_m l_i dxi and B_k(m,j) = (1/2) int l_m(xi_k(s)) l_j(s) ds, evaluated with an internal Gauss rule that is exact for the degree 2N integrands. Consequently sum_k matmul(P_k,R_k) = I (a polynomial split across both sub-edges is recovered) sum_m w_m (P_k g)_m = (1/2) sum_j w_j g_k(j) (discrete conservation) for any control node type. mortarP carries the 1/2 sub-edge Jacobian appropriate for projecting solution traces; when projecting surface-flux integrands (which carry the small side's nScale = nScale_big/2), scale by 2 (see MortarFluxCollect in the 2-D mapped data classes).
real(kind=prec), public, pointer, contiguous, dimension(:,:,:):: mortarR

Mortar restriction matrices for 2:1 nonconforming (mortar) interfaces.

An element edge trace u(xi) = sum_i u_i l_i(xi) on [-1,1] ("big" side) is restricted to the two half-edges ("small" sides) of its 2:1 neighbors. Sub-edge k=1 occupies xi in [-1,0] and k=2 occupies xi in [0,1], with the sub-edge coordinate s in [-1,1] related to the big-edge coordinate through xi_1(s) = (s-1)/2 and xi_2(s) = (s+1)/2.

mortarR(ii,i,k) = l_ii( xi_k(s_i) ), where s_i are the control points; following the SELF matrix convention, the restricted trace is u_k(i) = sum_ii mortarR(ii,i,k)*u(ii). Restriction of a degree N polynomial is exact.
real(kind=prec), public, pointer, contiguous, dimension(:):: qWeights

The quadrature weights for discrete integration. The quadradture weights depend on the type of controlPoints provided; one of Legendre-Gauss, Legendre-Gauss-Lobatto, Legendre-Gauss-Radau, Chebyshev-Gauss, Chebyshev-Gauss-Lobatto, Chebyshev-Gauss Radau, or Uniform. If Uniform, the quadrature weights are constant .
integer, public :: targetNodeType
real(kind=prec), public, pointer, contiguous, dimension(:):: targetPoints

The set of nodes in one dimension where data is to be interpolated to. To create higher dimension interpolation and differentiation operators, structured grids in two and three dimensions are created by tensor products of the targetPoints. In practice, the targetPoints are set to a uniformly distributed set of points between [-1,1] (computational space) to allow for interpolation from unevenly spaced quadrature points to a plotting grid.

Type-Bound Procedures

procedure, public :: CalculateBarycentricWeights

procedure, public :: CalculateDerivativeMatrix

procedure, public :: CalculateInterpolationMatrix

procedure, public :: CalculateLagrangePolynomials

  • public function CalculateLagrangePolynomials(this, sE) result(lAtS)

    Arguments

    TypeIntentOptionalAttributesName
    class(Lagrange_t) :: this
    real(kind=prec) :: sE

    Return Value real(kind=prec)(0:this%N)

procedure, public :: CalculateMortarMatrices

procedure, public :: Free => Free_Lagrange_t

  • public subroutine Free_Lagrange_t(this)

    Frees all memory (host and device) associated with an instance of the Lagrange_t class

    Arguments

    TypeIntentOptionalAttributesName
    class(Lagrange_t), intent(inout) :: this

    Lagrange_t class instance

procedure, public :: Init => Init_Lagrange_t

  • public subroutine Init_Lagrange_t(this, N, controlNodeType, M, targetNodeType)

    Initialize an instance of the Lagrange_t class On output, all of the attributes for the Lagrange_t class are allocated and values are initialized according to the number of control points, number of target points, and the types for the control and target nodes. If a GPU is available, device pointers for the Lagrange_t attributes are allocated and initialized.

    Arguments

    TypeIntentOptionalAttributesName
    class(Lagrange_t), intent(out) :: this

    Lagrange_t class instance
    integer, intent(in) :: N

    The number of control points for interpolant
    integer, intent(in) :: controlNodeType

    The integer code specifying the type of control points. Parameters are defined in SELF_Constants.f90. One of GAUSS(=1), GAUSS_LOBATTO(=2), or UNIFORM(=3)
    integer, intent(in) :: M

    The number of target points for the interpolant
    integer, intent(in) :: targetNodeType

    The integer code specifying the type of target points. Parameters are defined in SELF_Constants.f90. One of GAUSS(=1), GAUSS_LOBATTO(=2), or UNIFORM(=3)

procedure, public :: WriteHDF5 => WriteHDF5_Lagrange_t

  • public subroutine WriteHDF5_Lagrange_t(this, fileId)

    Arguments

    TypeIntentOptionalAttributesName
    class(Lagrange_t), intent(in) :: this
    integer(kind=HID_T), intent(in) :: fileId

Source Code

  type,public :: Lagrange_t
    !! A data structure for working with Lagrange Interpolating Polynomials in one, two, and three dimensions.
    !! The Lagrange data-structure stores the information necessary to interpolate between two
    !! sets of grid-points and to estimate the derivative of data at native grid points. Routines for
    !! multidimensional interpolation are based on the tensor product of 1-D interpolants. It is
    !! assumed that the polynomial degree (and the interpolation nodes) are the same in each direction.
    !! This assumption permits the storage of only one array of interpolation nodes and barycentric
    !! weights and is what allows this data structure to be flexible.

    integer :: N
      !! The number of control points.

    integer :: controlNodeType

    integer :: M
      !! The number of target points.

    integer :: targetNodeType

    real(prec),pointer,contiguous,dimension(:) :: controlPoints
      !! The set of nodes in one dimension where data is known.
      !! To create higher dimension interpolation and differentiation operators, structured grids in two and three
      !! dimensions are created by tensor products of the controlPoints. This design decision implies that all
      !! spectral element methods supported by the Lagrange class have the same polynomial degree in each
      !! computational/spatial dimension. In practice, the controlPoints are the Legendre-Gauss, Legendre-Gauss-Lobatto,
      !! Legendre-Gauss-Radau, Chebyshev-Gauss, Chebyshev-Gauss-Lobatto, or Chebyshev-Gauss-Radau quadrature points over
      !! the domain [-1,1] (computational space). The Init routine for this class restricts controlPoints to one of
      !! these quadrature types or uniform points on [-1,1].

    real(prec),pointer,contiguous,dimension(:) :: targetPoints
      !! The set of nodes in one dimension where data is to be interpolated to. To create higher dimension interpolation
      !! and differentiation operators, structured grids in two and three dimensions are created by tensor products of
      !! the targetPoints. In practice, the targetPoints are set to a uniformly distributed set of points between [-1,1]
      !! (computational space) to allow for interpolation from unevenly spaced quadrature points to a plotting grid.

    real(prec),pointer,contiguous,dimension(:) :: bWeights
      !! The barycentric weights that are calculated from the controlPoints and used for interpolation.

    real(prec),pointer,contiguous,dimension(:) :: qWeights
      !! The quadrature weights for discrete integration. The quadradture weights depend on the type of controlPoints
      !! provided; one of Legendre-Gauss, Legendre-Gauss-Lobatto, Legendre-Gauss-Radau, Chebyshev-Gauss,
      !! Chebyshev-Gauss-Lobatto, Chebyshev-Gauss Radau, or Uniform. If Uniform, the quadrature weights are constant
      !! $$dx = \frac{2.0}{N+1}$$.

    real(prec),pointer,contiguous,dimension(:,:) :: iMatrix
      !! The interpolation matrix (transpose) for mapping data from the control grid to the target grid.

    real(prec),pointer,contiguous,dimension(:,:) :: dMatrix
      !! The derivative matrix for mapping function nodal values to a nodal values of the derivative estimate. The
      !! dMatrix is based on a strong form of the derivative.

    real(prec),pointer,contiguous,dimension(:,:) :: dgMatrix
      !! The derivative matrix for mapping function nodal values to a nodal values of the derivative estimate. The dgMatrix is based
      !! on a weak form of the derivative. It must be used with bMatrix to account for boundary contributions in the weak form.

    real(prec),pointer,contiguous,dimension(:,:) :: dSplitMatrix
      !! The split-form derivative matrix D_split = D - 0.5*M^{-1}*B, where B is the SBP boundary operator.
      !!
      !! D_split is skew-symmetric under the M inner product: M*D_split + D_split^T*M = 0. Unlike D, it is
      !! NOT an SBP operator (the weighted symmetric part is zero, not the boundary term B). This property
      !! makes D_split ideal for the EC-DGSEM split-form volume integral:
      !!
      !!   du/dt = -(2/J) sum_n D_split[n,i] * F_EC(u_i, u_n)  +  (1/J) M^{-1} B^T f_Riemann
      !!
      !! The skew-symmetry of the volume term guarantees it contributes zero to the entropy rate, so all
      !! entropy change passes through the surface term. An entropy-dissipative Riemann solver then makes
      !! the full scheme entropy-stable.
      !!
      !! Using D_split in the volume combined with plain f_Riemann on the surface is algebraically identical
      !! to using D in the volume with the penalty (f_Riemann - f_local) on the surface (Trixi.jl convention).
      !!
      !! In SELF index convention (dSplitMatrix(ii,i) = D_split[i-1, ii-1]):
      !!   dSplitMatrix(ii,i) = dMatrix(ii,i)
      !!                       - 0.5*(bMatrix(i,2)*bMatrix(ii,2)
      !!                              - bMatrix(i,1)*bMatrix(ii,1)) / qWeights(i)

    real(prec),pointer,contiguous,dimension(:,:) :: bMatrix
      !! The boundary interpolation matrix that is used to map a grid of nodal values at the control points to the element boundaries.

    real(prec),pointer,contiguous,dimension(:,:,:) :: mortarR
      !! Mortar restriction matrices for 2:1 nonconforming (mortar) interfaces.
      !!
      !! An element edge trace u(xi) = sum_i u_i l_i(xi) on [-1,1] ("big" side) is restricted to
      !! the two half-edges ("small" sides) of its 2:1 neighbors. Sub-edge k=1 occupies
      !! xi in [-1,0] and k=2 occupies xi in [0,1], with the sub-edge coordinate s in [-1,1]
      !! related to the big-edge coordinate through xi_1(s) = (s-1)/2 and xi_2(s) = (s+1)/2.
      !!
      !! mortarR(ii,i,k) = l_ii( xi_k(s_i) ), where s_i are the control points; following the
      !! SELF matrix convention, the restricted trace is u_k(i) = sum_ii mortarR(ii,i,k)*u(ii).
      !! Restriction of a degree N polynomial is exact.

    real(prec),pointer,contiguous,dimension(:,:,:) :: mortarP
      !! Mortar L2-projection matrices, the adjoints of mortarR under the exact L2 inner
      !! product on [-1,1]: given traces g_k(s) on the two sub-edges, the projected big-side
      !! trace is u(i) = sum_k sum_ii mortarP(ii,i,k)*g_k(ii).
      !!
      !! The matrices are built as P_k = M^{-1} B_k with the exact (dense) 1-D mass matrix
      !! M(m,i) = int l_m l_i dxi and B_k(m,j) = (1/2) int l_m(xi_k(s)) l_j(s) ds, evaluated
      !! with an internal Gauss rule that is exact for the degree 2N integrands. Consequently
      !!   sum_k matmul(P_k,R_k) = I   (a polynomial split across both sub-edges is recovered)
      !!   sum_m w_m (P_k g)_m = (1/2) sum_j w_j g_k(j)   (discrete conservation)
      !! for any control node type. mortarP carries the 1/2 sub-edge Jacobian appropriate for
      !! projecting *solution* traces; when projecting surface-flux integrands (which carry the
      !! small side's nScale = nScale_big/2), scale by 2 (see MortarFluxCollect in the 2-D
      !! mapped data classes).

  contains

    procedure,public :: Init => Init_Lagrange_t
    procedure,public :: Free => Free_Lagrange_t

    procedure,public :: WriteHDF5 => WriteHDF5_Lagrange_t

    procedure,public :: CalculateBarycentricWeights
    procedure,public :: CalculateInterpolationMatrix
    procedure,public :: CalculateDerivativeMatrix
    procedure,public :: CalculateLagrangePolynomials
    procedure,public :: CalculateMortarMatrices

  endtype Lagrange_t