This function calculates the sigma for the 2d box shape sponge layer from coord for polynomial n6. Sponge profile: where, \sigma_A - sponge strength, L - thickness, x0 - start of sponge.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(tem_spongeLayer_box_type) | :: | me |
Spatial sponge layer 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 spongeLayer_box2d_roundCornerPolyn6_for_coord(me, coord, n) &
& result(res)
! --------------------------------------------------------------------------
!> Spatial sponge layer 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(2), extent(2), coordLoc(2), normal
real(kind=rk) :: proj_len1, proj_len2
real(kind=rk) :: vec_corMin(2), vec_corMax(2)
real(kind=rk) :: vec_corMinSqr(2), vec_corMaxSqr(2)
real(kind=rk) :: rad, const_fac, box_max(2)
real(kind=rk) :: corInRad, corOutRad, corMin(2), corMax(2)
logical :: isCorner
! --------------------------------------------------------------------------
origin(:) = me%origin(1:2)
extent(:) = me%extent(1:2)
box_max(:) = origin(:) + extent(:)
corInRad = me%corner_radius
corOutRad = corInRad + me%thickness
corMin(:) = origin(:) + corInRad
corMax(:) = box_max(:) - corInRad
const_fac = 729_rk/(16_rk*me%thickness**6)
do i = 1,n
coordLoc = coord(i,1:2)
vec_corMin(:) = coordLoc(:) - corMin(:)
vec_corMax(:) = coordLoc(:) - corMax(:)
vec_corMinSqr(:) = vec_corMin(:)**2
vec_corMaxSqr(:) = vec_corMax(:)**2
proj_len1 = 0_rk
proj_len2 = 0_rk
normal = 1_rk
rad = 0_rk
isCorner = .true.
if (coordLoc(1) < corMin(1) .and. coordLoc(2) < corMin(2)) then
! South-West: -x,-y
rad = sqrt(vec_corMinSqr(1) + vec_corMinSqr(2))
else if (coordLoc(1) < corMin(1) .and. coordLoc(2) > corMax(2)) then
! North-West: -x, y
rad = sqrt(vec_corMinSqr(1) + vec_corMaxSqr(2))
else if (coordLoc(1) > corMax(1) .and. coordLoc(2) < corMin(2)) then
! South-East: x, -y
rad = sqrt(vec_corMaxSqr(1) + vec_corMinSqr(2))
else if (coordLoc(1) > corMax(1) .and. coordLoc(2) > corMax(2)) then
! North-East: x, y
rad = sqrt(vec_corMaxSqr(1) + vec_corMaxSqr(2))
else
isCorner = .false.
if (coordLoc(1) < origin(1)) then ! West: -x
normal = -1_rk
rad = coordLoc(1) - origin(1)
else if (coordLoc(2) < origin(2)) then ! South: -y
normal = -1_rk
rad = coordLoc(2) - origin(2)
else if (coordLoc(1) > box_max(1)) then ! East: x
normal = 1_rk
rad = coordLoc(1) - box_max(1)
else if (coordLoc(2) > box_max(2)) then ! North: y
normal = 1_rk
rad = coordLoc(2) - box_max(2)
end if
end if
if (isCorner) then
proj_len1 = rad - corInRad
proj_len2 = corOutRad - rad
else
proj_len1 = rad*normal
proj_len2 = (me%thickness*normal - rad)*normal
end if
if (proj_len1 > 0) then
sigma = const_fac * proj_len2**2 * (proj_len1**4)
res(i) = min(1.0_rk, sigma) * me%dampFactor
else
res(i) = 0.0_rk
end if
end do
end function spongeLayer_box2d_roundCornerPolyn6_for_coord