subroutine CalculateMortarMatrices(this)
implicit none
class(Lagrange_t),intent(inout) :: this
! Local
integer :: i,ii,k,q,nq
real(prec) :: gpts(0:this%N+1)
real(prec) :: gwts(0:this%N+1)
real(prec) :: lbig(0:this%N)
real(prec) :: lsub(0:this%N)
real(prec) :: xik
real(real64) :: mmat(1:this%N+1,1:this%N+1)
real(real64) :: mwork(1:this%N+1,1:this%N+1)
real(real64) :: bmat(1:this%N+1,1:this%N+1,1:2)
real(real64) :: pmat(1:this%N+1,1:this%N+1)
nq = this%N+1 ! Internal quadrature rule has nq+1 = N+2 points; exact for degree 2N+3
call LegendreQuadrature(nq,gpts,gwts,GAUSS)
mmat = 0.0_real64
bmat = 0.0_real64
do q = 0,nq
! Values of the control-point Lagrange basis at the internal quadrature point; used
! both for the sub-edge (small side) basis in B_k and for the big-side basis in M.
lsub = this%CalculateLagrangePolynomials(gpts(q))
do i = 1,this%N+1
do ii = 1,this%N+1
mmat(ii,i) = mmat(ii,i)+ &
real(gwts(q),real64)*real(lsub(ii-1),real64)*real(lsub(i-1),real64)
enddo
enddo
do k = 1,2
! Map the sub-edge coordinate s = gpts(q) to the big-edge coordinate
if(k == 1) then
xik = 0.5_prec*(gpts(q)-1.0_prec) ! xi in [-1,0]
else
xik = 0.5_prec*(gpts(q)+1.0_prec) ! xi in [0,1]
endif
lbig = this%CalculateLagrangePolynomials(xik)
! B_k(m,j) = (1/2) int l_m(xi_k(s)) l_j(s) ds ; the 1/2 is the sub-edge Jacobian
do i = 1,this%N+1 ! j index (sub-edge basis)
do ii = 1,this%N+1 ! m index (big-edge basis)
bmat(ii,i,k) = bmat(ii,i,k)+ &
0.5_real64*real(gwts(q),real64)* &
real(lbig(ii-1),real64)*real(lsub(i-1),real64)
enddo
enddo
enddo
enddo
do k = 1,2
! Solve M P_k = B_k for the projection matrix; the solver overwrites its
! matrix argument with the factorization, so factor a fresh copy of M
mwork = mmat
pmat = bmat(1:this%N+1,1:this%N+1,k)
call SolveLinearSystems(this%N+1,mwork,pmat)
! Store the transpose, following the SELF operator convention
! result(i) = sum_ii matrix(ii,i)*input(ii)
do i = 1,this%N+1
do ii = 1,this%N+1
this%mortarP(ii,i,k) = real(pmat(i,ii),prec)
enddo
enddo
! Restriction: evaluate the big-side basis at the sub-edge images of the control points
do i = 1,this%N+1
if(k == 1) then
xik = 0.5_prec*(this%controlPoints(i)-1.0_prec)
else
xik = 0.5_prec*(this%controlPoints(i)+1.0_prec)
endif
lbig = this%CalculateLagrangePolynomials(xik)
do ii = 1,this%N+1
this%mortarR(ii,i,k) = lbig(ii-1)
enddo
enddo
enddo
endsubroutine CalculateMortarMatrices