function ScalarSurfaceIntegral(mesh, zoneID, integralType, iter) result(val)
!
! -----------------------------------------------------------
! This function computes scalar integrals, that is, those
! in the form:
! val = \int \vec{v}ยท\vec{n}dS
! Implemented integrals are:
! * Surface: computes the zone surface.
! * Mass flow: computes the mass flow across the zone.
! * Flow: computes the volumetric flow across the zone.
! -----------------------------------------------------------
!
implicit none
class(HexMesh), intent(inout), target :: mesh
integer, intent(in) :: zoneID
integer, intent(in) :: integralType, iter
real(kind=RP) :: val, localval
!
! ---------------
! Local variables
! ---------------
!
integer :: zonefID, fID, eID, fIDs(6), ierr
class(Element), pointer :: elements(:)
!
! Initialization
! --------------
val = 0.0_RP
!
! Loop the zone to get faces and elements
! ---------------------------------------
elements => mesh % elements
!$omp parallel private(fID, eID, fIDs) shared(elements,mesh,NodalStorage,zoneID,integralType,val,&
!$omp& computeGradients)
!$omp single
do zonefID = 1, mesh % zones(zoneID) % no_of_faces
fID = mesh % zones(zoneID) % faces(zonefID)
eID = mesh % faces(fID) % elementIDs(1)
fIDs = mesh % elements(eID) % faceIDs
!$omp task depend(inout:elements(eID))
call elements(eID) % ProlongSolutionToFaces(NCONS, mesh % faces(fIDs(1)),&
mesh % faces(fIDs(2)),&
mesh % faces(fIDs(3)),&
mesh % faces(fIDs(4)),&
mesh % faces(fIDs(5)),&
mesh % faces(fIDs(6)) )
if ( computeGradients ) then
call elements(eID) % ProlongGradientsToFaces(NGRAD, mesh % faces(fIDs(1)),&
mesh % faces(fIDs(2)),&
mesh % faces(fIDs(3)),&
mesh % faces(fIDs(4)),&
mesh % faces(fIDs(5)),&
mesh % faces(fIDs(6)) )
end if
!$omp end task
end do
!$omp end single
!
! Loop the zone to get faces and elements
! ---------------------------------------
!$omp do private(fID) reduction(+:val) schedule(runtime)
do zonefID = 1, mesh % zones(zoneID) % no_of_faces
!
! Face global ID
! --------------
fID = mesh % zones(zoneID) % faces(zonefID)
!
! Compute the integral
! --------------------
val = val + ScalarSurfaceIntegral_Face(mesh % faces(fID), integralType)
end do
!$omp end do
!$omp end parallel
#ifdef _HAS_MPI_
localval = val
call mpi_allreduce(localval, val, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD, ierr)
#endif
end function ScalarSurfaceIntegral
