MortarExchange_MappedVector2D_t Subroutine

public subroutine MortarExchange_MappedVector2D_t(this, mesh)

Fills the extBoundary attribute on all sides participating in a 2:1 nonconforming (mortar) interface; vector analogue of the scalar MortarExchange (see MappedScalar2D_t for the algorithm description).

Arguments

TypeIntentOptionalAttributesName
class(MappedVector2D_t), intent(inout) :: this
type(Mesh2D), intent(inout) :: mesh

Contents


Source Code

  subroutine MortarExchange_MappedVector2D_t(this,mesh)
    !! Fills the extBoundary attribute on all sides participating in a 2:1
    !! nonconforming (mortar) interface; vector analogue of the scalar MortarExchange
    !! (see MappedScalar2D_t for the algorithm description).
    implicit none
    class(MappedVector2D_t),intent(inout) :: this
    type(Mesh2D),intent(inout) :: mesh
    ! Local
    integer :: m,k,ivar,idir,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%MPIMortarExchangeAsync(mesh)
    endif

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

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

      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,idir)
            enddo
            if(flip == 0) then
              this%extBoundary(i,sS,eS-offset,ivar,idir) = fm
            else
              this%extBoundary(N+2-i,sS,eS-offset,ivar,idir) = 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,idir)
            enddo
          enddo
          this%extBoundary(i,sB,eB-offset,ivar,idir) = fm
        enddo
      endif

    enddo

  endsubroutine MortarExchange_MappedVector2D_t