viscSpongelayer_box_for_coord Function

private function viscSpongelayer_box_for_coord(me, coord, n) result(res)

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.

Arguments

Type IntentOptional 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 real(kind=rk), (n)

return value


Called by

proc~~viscspongelayer_box_for_coord~~CalledByGraph proc~viscspongelayer_box_for_coord viscSpongelayer_box_for_coord interface~tem_viscspongelayer_box_for tem_viscSpongeLayer_box_for interface~tem_viscspongelayer_box_for->proc~viscspongelayer_box_for_coord proc~viscspongelayer_box_for_treeids viscSpongelayer_box_for_treeIDs interface~tem_viscspongelayer_box_for->proc~viscspongelayer_box_for_treeids proc~viscspongelayer_box_for_treeids->proc~viscspongelayer_box_for_coord proc~tem_spatial_for_coord tem_spatial_for_coord proc~tem_spatial_for_coord->interface~tem_viscspongelayer_box_for proc~tem_spatial_for_treeids tem_spatial_for_treeIDs proc~tem_spatial_for_treeids->interface~tem_viscspongelayer_box_for interface~tem_spatial_for tem_spatial_for interface~tem_spatial_for->proc~tem_spatial_for_coord interface~tem_spatial_for->proc~tem_spatial_for_treeids proc~tem_spatial_scalar_for_index tem_spatial_scalar_for_index interface~tem_spatial_for->proc~tem_spatial_scalar_for_index proc~tem_spatial_scalar_for_index->proc~tem_spatial_for_coord proc~tem_spacetime_for_coord tem_spacetime_for_coord proc~tem_spacetime_for_coord->interface~tem_spatial_for proc~tem_spacetime_for_treeids tem_spacetime_for_treeIDs proc~tem_spacetime_for_treeids->interface~tem_spatial_for proc~tem_spacetime_scalar_for_index tem_spacetime_scalar_for_index proc~tem_spacetime_scalar_for_index->interface~tem_spatial_for proc~tem_spacetime_vector_for_coord tem_spacetime_vector_for_coord proc~tem_spacetime_vector_for_coord->interface~tem_spatial_for proc~tem_spacetime_vector_for_index tem_spacetime_vector_for_index proc~tem_spacetime_vector_for_index->interface~tem_spatial_for proc~tem_spacetime_vector_for_treeids tem_spacetime_vector_for_treeIDs proc~tem_spacetime_vector_for_treeids->interface~tem_spatial_for proc~tem_spatial_scalar_storeval tem_spatial_scalar_storeVal proc~tem_spatial_scalar_storeval->interface~tem_spatial_for proc~tem_spatial_vector_storeval tem_spatial_vector_storeVal proc~tem_spatial_vector_storeval->interface~tem_spatial_for

Source Code

  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