subroutine BoundaryFlux_DGModel3D_t(this)
! this method uses an linear upwind solver for the
! advective flux and the bassi-rebay method for the
! diffusive fluxes
implicit none
class(DGModel3D_t),intent(inout) :: this
! Local
integer :: i,j,k,iel
real(prec) :: sL(1:this%nvar),sR(1:this%nvar)
real(prec) :: dsdx(1:this%nvar,1:3)
real(prec) :: nhat(1:3),nmag
do iEl = 1,this%solution%nElem
do k = 1,6
do j = 1,this%solution%interp%N+1
do i = 1,this%solution%interp%N+1
! Get the boundary normals on cell edges from the mesh geometry
nhat = this%geometry%nHat%boundary(i,j,k,iEl,1,1:3)
sL = this%solution%boundary(i,j,k,iel,1:this%nvar) ! interior solution
sR = this%solution%extboundary(i,j,k,iel,1:this%nvar) ! exterior solution
dsdx = this%solutiongradient%avgboundary(i,j,k,iel,1:this%nvar,1:3)
nmag = this%geometry%nScale%boundary(i,j,k,iEl,1)
this%flux%boundaryNormal(i,j,k,iEl,1:this%nvar) = this%riemannflux3d(sL,sR,dsdx,nhat)*nmag
enddo
enddo
enddo
enddo
endsubroutine BoundaryFlux_DGModel3D_t