This function computes 3d parabola profile from coord 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:
use projection of point on line and point on plane
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
type(tem_canonicalND_type) | :: | me |
contains parameter for 3d parabola |
|||
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 |
the number of return values |
return value of a function
function tem_spatial_parabol3d_for_coord( me, coord, n ) result(res)
! -------------------------------------------------------------------- !
!> contains parameter for 3d parabola
type( tem_canonicalND_type ) :: me
!> the number of return values
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 of a function
real(kind=rk) :: res(n)
! -------------------------------------------------------------------- !
real(kind=rk) :: alpha, beta, diff(3)
real(kind=rk) :: vecAsqr, vecBsqr
real(kind=rk) :: center(3), halfvec(2,3)
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
!distance between parabola center and barycentric coordinates
diff = coord(iDir,:) - 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.d0) .or. (abs(beta) .gt. 1.d0) ) then
res(iDir) = 0.0_rk
end if
end do
end function tem_spatial_parabol3d_for_coord