MortarExchange_MappedScalar2D_t Subroutine

  subroutine MortarExchange_MappedScalar2D_t(this,mesh)
    !! Fills the extBoundary attribute on all sides that participate in a 2:1
    !! nonconforming (mortar) interface.
    !!
    !! Each small side receives the big side's trace restricted (exactly) to its
    !! sub-edge with the interp's mortarR matrices. The big side receives the L2
    !! projection of the two small-side traces onto its trace space with the interp's
    !! mortarP matrices. Traces are staged in mortarBuff in the big side's edge
    !! orientation; the flip recorded in mesh%mortarInfo translates between the small
    !! sides' orientations and the big side's orientation. Must be called after
    !! BoundaryInterp and may be called before, after, or without SideExchange (the
    !! side sets touched are disjoint from conforming sides).
    implicit none
    class(MappedScalar2D_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))
      this%mortarBuff = 0.0_prec
    endif

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

    ! Stage rank-local traces in the big side's edge orientation
    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
          this%mortarBuff(i,1,m,ivar) = this%boundary(i,sB,eB-offset,ivar)
          this%mortarBuff(i,2,m,ivar) = this%boundary(i,sB,eB-offset,ivar)
        enddo
      endif

      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) = this%boundary(i,sS,eS-offset,ivar)
            enddo
          else
            do i = 1,N+1
              this%mortarBuff(i,2+k,m,ivar) = this%boundary(N+2-i,sS,eS-offset,ivar)
            enddo
          endif
        endif
      enddo

    enddo

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

      ! Reorient small-side traces 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)
                enddo
                do i = 1,N+1
                  this%mortarBuff(i,2+k,m,ivar) = extBuff(i)
                enddo
              endif
            enddo
          endif
        enddo
      enddo
    endif

    ! Compute external states :
    !  small sides get the restricted big-side trace (exact),
    !  the big side gets the L2 projection of the small-side traces
    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
          do i = 1,N+1
            fm = 0.0_prec
            do ii = 1,N+1
              fm = fm+this%interp%mortarR(ii,i,k)*this%mortarBuff(ii,k,m,ivar)
            enddo
            if(flip == 0) then
              this%extBoundary(i,sS,eS-offset,ivar) = fm
            else
              this%extBoundary(N+2-i,sS,eS-offset,ivar) = fm
            endif
          enddo
        endif
      enddo

      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)
            enddo
          enddo
          this%extBoundary(i,sB,eB-offset,ivar) = fm
        enddo
      endif

    enddo

  endsubroutine MortarExchange_MappedScalar2D_t