subroutine VectorDataReconstruction( IBM, elements, STLNum, integralType, iter, autosave, dt )
use TessellationTypes
use MappedGeometryClass
use IBMClass
use OrientedBoundingBox
use KDClass
use MPI_Process_Info
use MPI_IBMUtilities
use omp_lib
#ifdef _HAS_MPI_
use mpi
#endif
!
! -----------------------------------------------------------------------------------------
! This function computes Vector integrals, that is, those
! in the form:
! val = \int \vec{v}ยท\vec{n}dS
! The data at the boundary point (BP) is computed through a Inverse Distance Weight
! procedure.
! -----------------------------------------------------------------------------------------
implicit none
!-arguments---------------------------------------------------------------------------------
type(IBM_type), intent(inout) :: IBM
type(element), intent(inout) :: elements(:)
integer, intent(in) :: integralType, STLNum, iter
real(kind=RP), intent(in) :: dt
!-local-variables---------------------------------------------------------------------------
real(kind=rp), allocatable :: Qsurf(:,:,:), U_xsurf(:,:,:), U_ysurf(:,:,:), U_zsurf(:,:,:)
integer :: i, j
logical :: found, autosave
if( .not. IBM% Integral(STLNum)% compute ) return
allocate( Qsurf(NCONS,IBM% NumOfInterPoints,NumOfVertices + 4), &
U_xsurf(NCONS,IBM% NumOfInterPoints,NumOfVertices + 4), &
U_ysurf(NCONS,IBM% NumOfInterPoints,NumOfVertices + 4), &
U_zsurf(NCONS,IBM% NumOfInterPoints,NumOfVertices + 4) )
call IBM% BandPoint_state( elements, STLNum, .true. )
if( .not. MPI_Process% isRoot ) return
if( IBM% stlSurfaceIntegrals(STLNum)% move ) then
if( IBM% stlSurfaceIntegrals(STLNum)% motionType .eq. ROTATION ) then
call IBM% stlSurfaceIntegrals(STLNum)% getRotationaMatrix( dt )
call OBB(STLNum)% STL_rotate( IBM% stlSurfaceIntegrals(STLNum), .true. )
elseif( IBM% stlSurfaceIntegrals(STLNum)% motionType .eq. LINEAR ) then
call IBM% stlSurfaceIntegrals(STLNum)% getDisplacement( dt )
call OBB(STLNum)% STL_translate( IBM% stlSurfaceIntegrals(STLNum), .true. )
end if
end if
!$omp parallel
!$omp do schedule(runtime) private(j,found,Qsurf,U_xsurf,U_ysurf,U_zsurf)
do i = 1, IBM% stlSurfaceIntegrals(STLNum)% NumOfObjs
do j = 1, NumOfVertices + 4
call GetSurfaceState( IBM, IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i), IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j), STLNum )
Qsurf(:,:,j) = IBM% BandRegion(STLNum)% Q (:,IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j)% nearestPoints)
U_xsurf(:,:,j) = IBM% BandRegion(STLNum)% U_x(:,IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j)% nearestPoints)
U_ysurf(:,:,j) = IBM% BandRegion(STLNum)% U_y(:,IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j)% nearestPoints)
U_zsurf(:,:,j) = IBM% BandRegion(STLNum)% U_z(:,IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j)% nearestPoints)
end do
do j = 1, NumOfVertices + 4
IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j)% VectorValue = IntegratedVectorValue( Q = Qsurf(:,:,j), &
U_x = U_xsurf(:,:,j), &
U_y = U_ysurf(:,:,j), &
U_z = U_zsurf(:,:,j), &
vertex = IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% vertices(j), &
normal = IBM% stlSurfaceIntegrals(STLNum)% ObjectsList(i)% normal, &
y = IBM% IP_Distance, &
Wallfunction = IBM% Wallfunction, &
integralType = integralType, &
InterpolationType = IBM% InterpolationType )
end do
end do
!$omp end do
!$omp end parallel
deallocate( Qsurf, U_xsurf, U_ysurf, U_zsurf )
if( IBM% stl(STLNum)% move ) then
IBM% Integral(STLNum)% ListComputed = .false.
else
IBM% Integral(STLNum)% ListComputed = .true.
end if
if( autosave ) call GenerateVectormonitorTECfile( IBM, STLNum, integralType, iter )
end subroutine VectorDataReconstruction
