AnisoMortarExchange_MappedVector3D_t Subroutine

  subroutine AnisoMortarExchange_MappedVector3D_t(this,mesh)
    !! Fills the extBoundary attribute on all faces participating in an anisotropic
    !! mortar interface; vector analogue of the scalar AnisoMortarExchange (see
    !! MappedScalar3D_t for the algorithm description).
    implicit none
    class(MappedVector3D_t),intent(inout) :: this
    type(Mesh3D),intent(inout) :: mesh
    ! Local
    integer :: m,h,ivar,idir,i,j,i2,j2
    integer :: eB,eS,sS,flip
    integer :: rankId,offset,N
    integer,pointer :: elemtorank(:)
    real(prec) :: extBuff(1:this%interp%N+1,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%anisoMortarBuff)) then
      allocate(this%anisoMortarBuff(1:N+1,1:N+1,1:4,1:mesh%nAnisoMortars,1:this%nvar,1:3))
      this%anisoMortarBuff = 0.0_prec
    endif

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

    ! Stage rank-local traces in the big face's coordinates
    do concurrent(m=1:mesh%nAnisoMortars,ivar=1:this%nvar,idir=1:3)
      block
        integer :: i,j,i2,j2,h
        integer :: eB,sB,eS,sS,flip

        eB = mesh%anisoMortarInfo(1,m)
        if(elemtorank(eB) == rankId) then
          sB = mesh%anisoMortarInfo(2,m)
          do h = 1,2
            do j = 1,N+1
              do i = 1,N+1
                this%anisoMortarBuff(i,j,h,m,ivar,idir) = &
                  this%boundary(i,j,sB,eB-offset,ivar,idir)
              enddo
            enddo
          enddo
        endif

        do h = 1,2
          eS = mesh%anisoMortarInfo(2*h+1,m)
          if(elemtorank(eS) == rankId) then
            sS = mesh%anisoMortarInfo(2*h+2,m)/10
            flip = mesh%anisoMortarInfo(2*h+2,m)-10*sS
            do j = 1,N+1
              do i = 1,N+1
                call MortarFaceMap(i,j,N,flip,i2,j2)
                this%anisoMortarBuff(i,j,2+h,m,ivar,idir) = &
                  this%boundary(i2,j2,sS,eS-offset,ivar,idir)
              enddo
            enddo
          endif
        enddo
      endblock
    enddo

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

      ! Reorient small-face traces received over MPI into the big face's coordinates
      do idir = 1,3
        do ivar = 1,this%nvar
          do m = 1,mesh%nAnisoMortars
            eB = mesh%anisoMortarInfo(1,m)
            if(elemtorank(eB) == rankId) then
              do h = 1,2
                eS = mesh%anisoMortarInfo(2*h+1,m)
                sS = mesh%anisoMortarInfo(2*h+2,m)/10
                flip = mesh%anisoMortarInfo(2*h+2,m)-10*sS
                if(elemtorank(eS) /= rankId .and. flip /= 0) then
                  do j = 1,N+1
                    do i = 1,N+1
                      call MortarFaceMap(i,j,N,flip,i2,j2)
                      extBuff(i,j) = this%anisoMortarBuff(i2,j2,2+h,m,ivar,idir)
                    enddo
                  enddo
                  do j = 1,N+1
                    do i = 1,N+1
                      this%anisoMortarBuff(i,j,2+h,m,ivar,idir) = extBuff(i,j)
                    enddo
                  enddo
                endif
              enddo
            endif
          enddo
        enddo
      enddo
    endif

    ! Compute external states :
    !  small faces get the restricted big-face trace (exact),
    !  the big face gets the L2 projection of the small-face traces
    do concurrent(m=1:mesh%nAnisoMortars,ivar=1:this%nvar,idir=1:3)
      block
        integer :: i,j,ii,i2,j2,h,split
        integer :: eB,sB,eS,sS,flip
        real(prec) :: fm
        real(prec) :: acc(1:N+1,1:N+1)

        split = mesh%anisoMortarInfo(9,m)

        do h = 1,2
          eS = mesh%anisoMortarInfo(2*h+1,m)
          if(elemtorank(eS) == rankId) then
            sS = mesh%anisoMortarInfo(2*h+2,m)/10
            flip = mesh%anisoMortarInfo(2*h+2,m)-10*sS
            do j = 1,N+1
              do i = 1,N+1
                fm = 0.0_prec
                if(split == 1) then
                  do ii = 1,N+1
                    fm = fm+this%interp%mortarR(ii,i,h)* &
                         this%anisoMortarBuff(ii,j,h,m,ivar,idir)
                  enddo
                else
                  do ii = 1,N+1
                    fm = fm+this%interp%mortarR(ii,j,h)* &
                         this%anisoMortarBuff(i,ii,h,m,ivar,idir)
                  enddo
                endif
                call MortarFaceMap(i,j,N,flip,i2,j2)
                this%extBoundary(i2,j2,sS,eS-offset,ivar,idir) = fm
              enddo
            enddo
          endif
        enddo

        eB = mesh%anisoMortarInfo(1,m)
        if(elemtorank(eB) == rankId) then
          sB = mesh%anisoMortarInfo(2,m)
          acc = 0.0_prec
          do h = 1,2
            do j = 1,N+1
              do i = 1,N+1
                fm = 0.0_prec
                if(split == 1) then
                  do ii = 1,N+1
                    fm = fm+this%interp%mortarP(ii,i,h)* &
                         this%anisoMortarBuff(ii,j,2+h,m,ivar,idir)
                  enddo
                else
                  do ii = 1,N+1
                    fm = fm+this%interp%mortarP(ii,j,h)* &
                         this%anisoMortarBuff(i,ii,2+h,m,ivar,idir)
                  enddo
                endif
                acc(i,j) = acc(i,j)+fm
              enddo
            enddo
          enddo
          do j = 1,N+1
            do i = 1,N+1
              this%extBoundary(i,j,sB,eB-offset,ivar,idir) = acc(i,j)
            enddo
          enddo
        endif
      endblock
    enddo

  endsubroutine AnisoMortarExchange_MappedVector3D_t