AnisoMortarFluxCollect_MappedVector3D_t Subroutine

  subroutine AnisoMortarFluxCollect_MappedVector3D_t(this,mesh)
    !! Replaces the big-face boundaryNormal trace on each anisotropic mortar interface
    !! with the L2 projection of the two small faces' boundaryNormal traces.
    !!
    !! As in MortarFluxCollect, boundaryNormal holds f* . nHat * nScale. The small faces'
    !! nScale is one half of the big face's and the half's coordinate Jacobian is 1/2,
    !! so the projected big-face integrand is -2 * sum_h P_h g_h, with mortarP(:,:,h)
    !! along the split coordinate and the identity along the other, where g_h are the
    !! small-face integrands and the sign accounts for the opposing outward normals. The
    !! discrete surface integral of the big face then equals minus the sum of the small
    !! faces' to roundoff. Must be called after the model's BoundaryFlux and before the
    !! flux divergence is computed.
    implicit none
    class(MappedVector3D_t),intent(inout) :: this
    type(Mesh3D),intent(inout) :: mesh
    ! Local
    integer :: m,h,ivar,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%MPIAnisoMortarFluxAsync(mesh)
    endif

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

        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,1) = &
                  this%boundaryNormal(i2,j2,sS,eS-offset,ivar)
              enddo
            enddo
          endif
        enddo
      endblock
    enddo

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

      ! Reorient small-face integrands received over MPI into the big face's coordinates
      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,1)
                  enddo
                enddo
                do j = 1,N+1
                  do i = 1,N+1
                    this%anisoMortarBuff(i,j,2+h,m,ivar,1) = extBuff(i,j)
                  enddo
                enddo
              endif
            enddo
          endif
        enddo
      enddo
    endif

    ! Project the small-face integrands onto the big face's trace space. The factor of
    ! two converts the solution-space projection (mortarP carries the 1/2 half
    ! Jacobian) into the integrand-space projection; the sign accounts for the opposing
    ! outward normals.
    do concurrent(m=1:mesh%nAnisoMortars,ivar=1:this%nvar)
      block
        integer :: i,j,ii,h,split
        integer :: eB,sB
        real(prec) :: fm
        real(prec) :: acc(1:N+1,1:N+1)

        eB = mesh%anisoMortarInfo(1,m)
        if(elemtorank(eB) == rankId) then
          sB = mesh%anisoMortarInfo(2,m)
          split = mesh%anisoMortarInfo(9,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,1)
                  enddo
                else
                  do ii = 1,N+1
                    fm = fm+this%interp%mortarP(ii,j,h)* &
                         this%anisoMortarBuff(i,ii,2+h,m,ivar,1)
                  enddo
                endif
                acc(i,j) = acc(i,j)+fm
              enddo
            enddo
          enddo
          do j = 1,N+1
            do i = 1,N+1
              this%boundaryNormal(i,j,sB,eB-offset,ivar) = -2.0_prec*acc(i,j)
            enddo
          enddo
        endif
      endblock
    enddo

  endsubroutine AnisoMortarFluxCollect_MappedVector3D_t