This function computes 3d parabola profile from treeIDs of an element.
This profile is defined by element barycentric coordinate and 3d parabola parameters. 3D parabola profile at given plane is computed in the following way:
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(tem_canonicalND_type) | :: | me |
contains plane parameter for 3d parabola |
|||
integer(kind=long_k), | intent(in) | :: | treeIds(n) |
treeIds of elements in given level |
||
type(treelmesh_type), | intent(in) | :: | tree |
global treelm mesh |
||
integer, | intent(in) | :: | n |
number of return values |
return value of a function
function tem_spatial_parabol3d_for_treeIds( me, treeIds, tree, n ) result(res)
! -------------------------------------------------------------------- !
!> contains plane parameter for 3d parabola
type( tem_canonicalND_type ) :: me
!> global treelm mesh
type( treelmesh_type ), intent(in) ::tree
!> number of return values
integer, intent(in) :: n
!> treeIds of elements in given level
integer(kind=long_k), intent(in) :: treeIds(n)
!> return value of a function
real(kind=rk) :: res(n)
! -------------------------------------------------------------------- !
real(kind=rk) :: coord(3), alpha, beta, diff(3)
real(kind=rk) :: vecAsqr, vecBsqr
real(kind=rk) :: center(3), halfvec(2,3)
!loop variables
integer :: iDir, jDir
! -------------------------------------------------------------------- !
jDir = 0
do iDir = 1, 3
if( me%active( iDir )) then
jDir = jDir + 1
halfvec( jDir, : ) = me%vec( :, iDir ) / 2._rk
end if
end do
center = me%origin + halfvec(1, :) + halfvec(2, :)
vecAsqr = dot_product( halfvec(1, :), halfvec(1, :) )
vecBsqr = dot_product( halfvec(2, :), halfvec(2, :) )
!loop over number of return values
do iDir = 1, n
!barycentric coordinate
coord = tem_BaryOfId( tree, treeIds(iDir) )
!distance between parabola center and barycentric coordinates
diff = coord - center
!projection of diff in a plane on vecA
alpha = dot_product( diff, halfvec(1, :) ) / vecAsqr
!projection of diff in a plane on vecB
beta = dot_product( diff, halfvec(2, :) ) / vecBsqr
res( iDir ) = ( 1.0_rk-alpha ) * ( 1.0_rk+alpha ) &
& * ( 1.0_rk-beta ) * ( 1.0_rk+beta )
if ( abs(alpha) .gt. 1.0_rk .or. abs(beta) .gt. 1.0_rk ) then
res( iDir ) = 0.0_rk
end if
end do
end function tem_spatial_parabol3d_for_treeIds