MortarFluxCollect_MappedVector2D_t Subroutine

  subroutine MortarFluxCollect_MappedVector2D_t(this,mesh)
    !! Replaces the big-side boundaryNormal trace on each mortar interface with the L2
    !! projection of the two small sides' boundaryNormal traces.
    !!
    !! boundaryNormal holds the Riemann-solved surface-flux integrand f* . nHat * nScale
    !! (see BoundaryFlux in the DG models). Because the small sides' nScale is half the
    !! big side's and the sub-edge coordinate Jacobian is 1/2, the projected big-side
    !! integrand is -2 * sum_k P_k g_k, where g_k are the small-side integrands and the
    !! sign accounts for the opposing outward normals. With this choice, the discrete
    !! surface integral of the big side equals minus the sum of the small sides'
    !! discrete surface integrals to roundoff, so the mortar interface is discretely
    !! conservative. Must be called after the model's BoundaryFlux and before the flux
    !! divergence is computed.
    implicit none
    class(MappedVector2D_t),intent(inout) :: this
    type(Mesh2D),intent(inout) :: mesh
    ! Local
    integer :: m,k,ivar,i,ii
    integer :: eB,sB,eS,sS,flip
    integer :: rankId,offset,N
    integer,pointer :: elemtorank(:)
    real(prec) :: fm
    real(prec) :: extBuff(1:this%interp%N+1)

    ! See https://github.com/FluidNumerics/SELF/issues/54 for the reason behind
    ! this pointer alias
    elemtorank => mesh%decomp%elemToRank(:)
    rankId = mesh%decomp%rankId
    offset = mesh%decomp%offsetElem(rankId+1)
    N = this%interp%N

    if(.not. allocated(this%mortarBuff)) then
      allocate(this%mortarBuff(1:N+1,1:4,1:mesh%nMortars,1:this%nvar,1:2))
      this%mortarBuff = 0.0_prec
    endif

    if(mesh%decomp%mpiEnabled) then
      call this%MPIMortarFluxAsync(mesh)
    endif

    ! Stage rank-local small-side integrands in the big side's edge orientation
    do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar)

      do k = 1,2
        eS = mesh%mortarInfo(2*k+1,m)
        if(elemtorank(eS) == rankId) then
          sS = mesh%mortarInfo(2*k+2,m)/10
          flip = mesh%mortarInfo(2*k+2,m)-10*sS
          if(flip == 0) then
            do i = 1,N+1
              this%mortarBuff(i,2+k,m,ivar,1) = this%boundaryNormal(i,sS,eS-offset,ivar)
            enddo
          else
            do i = 1,N+1
              this%mortarBuff(i,2+k,m,ivar,1) = this%boundaryNormal(N+2-i,sS,eS-offset,ivar)
            enddo
          endif
        endif
      enddo

    enddo

    if(mesh%decomp%mpiEnabled) then
      call mesh%decomp%FinalizeMPIExchangeAsync()

      ! Reorient small-side integrands received over MPI into the big side's orientation
      do ivar = 1,this%nvar
        do m = 1,mesh%nMortars
          eB = mesh%mortarInfo(1,m)
          if(elemtorank(eB) == rankId) then
            do k = 1,2
              eS = mesh%mortarInfo(2*k+1,m)
              sS = mesh%mortarInfo(2*k+2,m)/10
              flip = mesh%mortarInfo(2*k+2,m)-10*sS
              if(elemtorank(eS) /= rankId .and. flip == 1) then
                do i = 1,N+1
                  extBuff(i) = this%mortarBuff(N+2-i,2+k,m,ivar,1)
                enddo
                do i = 1,N+1
                  this%mortarBuff(i,2+k,m,ivar,1) = extBuff(i)
                enddo
              endif
            enddo
          endif
        enddo
      enddo
    endif

    ! Project the small-side integrands onto the big side's trace space. The factor of
    ! two converts the solution-space projection (mortarP carries the 1/2 sub-edge
    ! Jacobian) into the integrand-space projection; the sign accounts for the opposing
    ! outward normals.
    do concurrent(m=1:mesh%nMortars,ivar=1:this%nvar)

      eB = mesh%mortarInfo(1,m)
      if(elemtorank(eB) == rankId) then
        sB = mesh%mortarInfo(2,m)
        do i = 1,N+1
          fm = 0.0_prec
          do k = 1,2
            do ii = 1,N+1
              fm = fm+this%interp%mortarP(ii,i,k)*this%mortarBuff(ii,2+k,m,ivar,1)
            enddo
          enddo
          this%boundaryNormal(i,sB,eB-offset,ivar) = -2.0_prec*fm
        enddo
      endif

    enddo

  endsubroutine MortarFluxCollect_MappedVector2D_t