Skip to content

Add orthogonal and isparallel function #195

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lijas wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Add orthogonal and isparallel function #195

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion src/Tensors.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ export tdot, dott, dotdot
export hessian, gradient, curl, divergence, laplace
export @implement_gradient
export basevec, eᵢ
export rotate, rotation_tensor
export rotate, rotation_tensor, orthogonal
export tovoigt, tovoigt!, fromvoigt, tomandel, tomandel!, frommandel
#########
# Types #
Expand Down
29 changes: 29 additions & 0 deletions src/math_ops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -378,3 +378,32 @@ function rotate(x::SymmetricTensor{4}, args...)
R = rotation_tensor(args...)
return unsafe_symmetric(otimesu(R, R) ⊡ x ⊡ otimesu(R', R'))
end

isparallel(v1::Vec{3}, v2::Vec{3}) = isapprox(Tensors.cross(v1,v2), zero(v1), atol = 1e-14)
Copy link
Member

@KnutAM KnutAM Mar 7, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we get rid of the absolute tolerance by doing this (which also generalizes to second order)

Suggested change
isparallel(v1::Vec{3}, v2::Vec{3}) = isapprox(Tensors.cross(v1,v2), zero(v1), atol = 1e-14)
isparallel(v1::Vec, v2::Vec) = abs(dot(v1,v2)) ≈ norm(v1)*norm(v2)
isparallel(a::SecondOrderTensor, b::SecondOrderTensor) = abs(dcontract(a,b)) ≈ norm(a)*norm(b)

Copy link
Member

@fredrikekre fredrikekre Mar 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

≈ also uses a tolerance though, right?

Copy link

@dkarrasch dkarrasch Mar 7, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This cries for a user-defined tolerance with some default value. Think of it conversely: any tolerance describes a cone around one vector in which all vectors are considered as parallel to the reference vector. Same with orthogonality.

Copy link
Member

@KnutAM KnutAM Mar 7, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

≈ also uses a tolerance though, right?

Yes, but a relative one, right? Which works since we are not comparing zeros. I agree with @dkarrasch though that an optional input should be used instead, perhaps isparallel(a,b; kwargs...) = isapprox(abs(dot(v1,v2)), norm(v1)*norm(v2); kwargs...)?
(I guess a pure angular tolerance could have numerical issues around zero-length vectors)

isparallel(v1::Vec{2}, v2::Vec{2}) = isapprox(v1[2]*v2[1], v1[1]*v2[2], atol = 1e-14)

"""
orthogonal(u::Vec)

Return a orthogonal vector `v` to u (v ⋅ u = 0.0)
"""
function orthogonal(u::Vec{3})
iszero(u) && error("Cannot construct a orthogonal vector to a vector with zero length")
r = rand(u)
while isparallel(r,u) || iszero(r) #Is this needed?
r = rand(u)
end
q = r - u*(r⋅u)/(u⋅u)
return q
end


function orthogonal(u::Vec{2})
iszero(u) && error("Cannot construct a orthogonal vector to a vector with zero length")
return Vec((u[2], -u[1]))
end

function orthogonal(u::Vec{1})
error("Cannot construct a orthogonal vector for a one dimensional vector")
#return Vec((0.0,))
end
19 changes: 19 additions & 0 deletions test/test_ops.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -319,5 +319,24 @@ end
@test isa(fromvoigt(SymmetricTensor{4,dim}, rand(num_components,num_components) .> 0.5), SymmetricTensor{4,dim,Bool})
end
end

@testsection "isparallel" begin
if dim == 3 || dim == 2
u = rand(Vec{dim})
v = 2u
@test Tensors.isparallel(u,v) == true
v = orthogonal(u)
@test Tensors.isparallel(u,v) == false
end
end

@testsection "orthogonal" begin
if dim == 3 || dim == 2
u = rand(Vec{dim})
v = orthogonal(u)
@test u⋅v ≈ 0.0 atol = 1e-14
end
end

end # of testsection
end # of testsection