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