This function calculates the sigma for the box viscosity spongelayer and multiply with targetState 'viscosity'. This function is currectly used to define viscosity sponge in musubi.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| type(tem_spongeLayer_box_type) | :: | me |
Spacetime function to evaluate |
|||
| real(kind=rk), | intent(in) | :: | coord(n,3) |
barycentric Ids of an elements. 1st index goes over number of elements and 2nd index goes over x,y,z coordinates |
||
| integer, | intent(in) | :: | n |
Number of arrays to return |
return value
function viscSpongelayer_box_for_coord(me, coord, n) &
& result(res)
! --------------------------------------------------------------------------
!> Spacetime function to evaluate
type(tem_spongeLayer_box_type) :: me
!> Number of arrays to return
integer, intent(in) :: n
!> barycentric Ids of an elements.
!! 1st index goes over number of elements and
!! 2nd index goes over x,y,z coordinates
real(kind=rk), intent( in ) :: coord(n,3)
!> return value
real(kind=rk) :: res(n)
! --------------------------------------------------------------------------
integer :: i
real(kind=rk) :: sigma, origin(3), extent(3), box_max(3), proj_len
real(kind=rk) :: coordLoc(3), normal, vec_min(3), vec_max(3)
real(kind=rk) :: rad, vec_minSqr(3), vec_maxSqr(3)
! --------------------------------------------------------------------------
origin(:) = me%origin
extent(:) = me%extent
box_max(:) = origin(:) + extent(:)
do i = 1,n
coordLoc = coord(i,:)
vec_min(:) = coordLoc(:) - origin(:)
vec_max(:) = coordLoc(:) - box_max(:)
vec_minSqr(:) = vec_min(:)**2
vec_maxSqr(:) = vec_max(:)**2
normal = 1.0_rk
rad = 0.0_rk
! Bottom-South-West: -x,-y,-z
if (all(coordLoc(:) < origin(:))) then
rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_minSqr(3)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
& .and. coordLoc(3) < origin(3)) then ! Bottom-South-East: x, -y, -z
rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_minSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
& .and. coordLoc(3) < origin(3)) then ! Bottom-North-West: -x, y, -z
rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_minSqr(3)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2) &
& .and. coordLoc(3) < origin(3)) then ! Bottom-North-East: x, y, -z
rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_minSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2) &
& .and. coordLoc(3) > box_max(3)) then ! Top-South-West: -x, -y, z
rad = sqrt( max(vec_minSqr(1), vec_minSqr(2), vec_maxSqr(3)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2) &
& .and. coordLoc(3) > box_max(3)) then ! Top-South-East: x, -y, z
rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2), vec_maxSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2) &
& .and. coordLoc(3) > box_max(3)) then ! Top-North-West: -x, y, z
rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )
else if (all(coordLoc(:) > box_max(:))) then ! Bottom-North-East: x, y, z
rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2), vec_maxSqr(3)) )
else if (coordLoc(2) < origin(2) .and. coordLoc(3) < origin(3)) then
! Botton-South: -y,-z
rad = sqrt( max(vec_minSqr(2), vec_minSqr(3)) )
else if (coordLoc(2) < origin(2) .and. coordLoc(3) > box_max(3)) then
! Top-South: -y, z
rad = sqrt( max(vec_minSqr(2), vec_maxSqr(3)) )
else if (coordLoc(2) > box_max(2) .and. coordLoc(3) < origin(3)) then
! Botom-North: -y, z
rad = sqrt( max(vec_maxSqr(2), vec_minSqr(3)) )
else if (coordLoc(2) > box_max(2) .and. coordLoc(3) > box_max(3)) then
! Top-North: y, z
rad = sqrt( max(vec_maxSqr(2), vec_maxSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(3) < origin(3)) then
! Botton-West: -x,-z
rad = sqrt( max(vec_minSqr(1), vec_minSqr(3)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(3) < origin(3)) then
! Bottom-East: x,-z
rad = sqrt( max(vec_maxSqr(1), vec_minSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(3) > box_max(3)) then
! Top-West: -x, z
rad = sqrt( max(vec_minSqr(1), vec_maxSqr(3)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(3) > box_max(3)) then
! Top-East: x, z
rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(3)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(2) < origin(2)) then
! South-West: -x,-y
rad = sqrt( max(vec_minSqr(1), vec_minSqr(2)) )
else if (coordLoc(1) < origin(1) .and. coordLoc(2) > box_max(2)) then
! North-West: -x, y
rad = sqrt( max(vec_minSqr(1), vec_maxSqr(2)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(2) < origin(2)) then
! South-East: x, -y
rad = sqrt( max(vec_maxSqr(1), vec_minSqr(2)) )
else if (coordLoc(1) > box_max(1) .and. coordLoc(2) > box_max(2)) then
! North-East: x, y
rad = sqrt( max(vec_maxSqr(1), vec_maxSqr(2)) )
else if (coordLoc(1) < origin(1)) then ! West: -x
normal = -1_rk
rad = vec_min(1)
else if (coordLoc(2) < origin(2)) then ! South: -y
normal = -1_rk
rad = vec_min(2)
else if (coordLoc(3) < origin(3)) then ! Bottom: -z
normal = -1_rk
rad = vec_min(3)
else if (coordLoc(1) > box_max(1)) then ! East: x
normal = 1_rk
rad = vec_max(1)
else if (coordLoc(2) > box_max(2)) then ! North: y
normal = 1_rk
rad = vec_max(2)
else if (coordLoc(3) > box_max(3)) then ! Top: z
normal = 1_rk
rad = vec_max(3)
end if
proj_len = rad*normal/me%thickness
sigma = 1.0_rk + (me%dampFactor-1.0_rk)*((proj_len)**me%dampExponent)
if (proj_len > 0 .and. proj_len < 1) then
res(i) = sigma * me%targetState(1)
else if (proj_len > 1) then
res(i) = me%dampFactor * me%targetState(1)
else
res(i) = me%targetState(1)
end if
end do
end function viscSpongelayer_box_for_coord