Units in `Vec` and `normal`

[pjabardo]

Vec objects have, necessarily length unit. It makes sense for Point but a vector could be something else.

As an example where this could be a problem, take the normal of a a surface, computed with method normal, It returns a vector with a unit length. If I need to compute the forces due to pressure on a surface with area A and normal n, it is kind of tricky. One common approach is to have a vector area where the direction is the normal and the norm is the area. With the current formulation this is not possible. Another common approach is to represent the the surface area as a scalar and the direction as a vector consisting of direction cosines which should necessarily be non-dimensional.

[juliohm]

Thank you for reporting the issue @pjabardo. We have some internal functions like unormalize and the StaticArrays.jl normalize that handle units differently. Can you clarify what is your suggestion?

[pjabardo]

I think Vec should be more liberal with units. It should not be restricted to length. If you take the cross product of two vectors, you get a SVector and not a Vec because of this restriction:

julia> u = Vec(1,0,0); v = Vec(0,1,0)
Vec(0.0 m, 1.0 m, 0.0 m)

julia> u × v
3-element StaticArraysCore.SVector{3, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(m^2,), 𝐋^2, nothing}}} with indices SOneTo(3):
 0.0 m^2
 0.0 m^2
 1.0 m^2

Now we have an object that can no longer be used except if manipulated somehow. And this is something useful in several operations.

I have no idea what the impact no longer restricting the units of a vector and I realize that these things are very tricky.

Take the dot product as another example

julia> u = Vec(1,0,0); v = Vec(1,1,0)
Vec(1.0 m, 1.0 m, 0.0 m)

julia> v⋅u  # This makes sense m^2
1.0 m^2

What if I just want to project v into the direction of u? I would like my result to be in m not m^2 . Of course I could use a Unitful aware vector/tensor package (or develop one of it is not available)

[juliohm]

Very useful inputs. We will think about it more deeply to see what we can do. Thanks for sharing the issue.

Em sex., 9 de ago. de 2024, 16:42, Paulo Jabardo @.***> escreveu:

I think Vec should be more liberal with units. It should not be restricted to length. If you take the cross product of two vectors, you get a SVector and not a Vec because of this restriction:

julia> u = Vec(1,0,0); v = Vec(0,1,0)Vec(0.0 m, 1.0 m, 0.0 m)

julia> u × v3-element StaticArraysCore.SVector{3, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(m^2,), 𝐋^2, nothing}}} with indices SOneTo(3): 0.0 m^2 0.0 m^2 1.0 m^2

Now we have an object that can no longer be used except if manipulated somehow. And this is something useful in several operations.

I have no idea what the impact no longer restricting the units of a vector and I realize that these things are very tricky.

Take the dot product as another example

julia> u = Vec(1,0,0); v = Vec(1,1,0)Vec(1.0 m, 1.0 m, 0.0 m)

julia> v⋅u # This makes sense m^21.0 m^2

What if I just want to project v into the direction of u? I would like my result to be in m not m^2 . Of course I could use a Unitful aware vector/tensor package (or develop one of it is not available)

— Reply to this email directly, view it on GitHub https://github.com/JuliaGeometry/Meshes.jl/issues/996#issuecomment-2278628677, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZQW3LMZDNZBFEV76ABZF3ZQULTBAVCNFSM6AAAAABMI64L5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDENZYGYZDQNRXG4 . You are receiving this because you commented.Message ID: @.***>

[juliohm]

@pjabardo we are trying to find a possible adjustment to the current design. At present, the Vec is a geometric vector, which for us means a direction in an Euclidean space that represents the physical world. In this space the Cartesian coordinates are given as lengths along the X, Y and Z axes.

Let's consider your examples...

Cross product

If you take the cross product of two vectors, you get a SVector and not a Vec because of this restriction:

julia> u = Vec(1,0,0); v = Vec(0,1,0)
Vec(0.0 m, 1.0 m, 0.0 m)

julia> u × v
3-element StaticArraysCore.SVector{3, Quantity{Float64, 𝐋^2, Unitful.FreeUnits{(m^2,), 𝐋^2, nothing}}} with indices SOneTo(3):
 0.0 m^2
 0.0 m^2
 1.0 m^2

Now we have an object that can no longer be used except if manipulated somehow. And this is something useful in several operations.

As per our current design, this result is no longer a geometric vector because lengths can't have units m^2. Say we allowed Vec with m^2, how would you add this vector to a point, to obtain a new point?

Point(1m, 2m, 3m) + Vec(1m^2, 2m^2, 3m^2)

What other operations do you have in mind?

Some random ideas come to mind:

1. Ignore the units of Vec for addition purposes, matching the units of the Point
2. Always remove units from Vec at construction time, and use the raw values across all algorithms

Dot product

Take the dot product as another example

julia> u = Vec(1,0,0); v = Vec(1,1,0)
Vec(1.0 m, 1.0 m, 0.0 m)

julia> v⋅u  # This makes sense m^2
1.0 m^2

What if I just want to project v into the direction of u? I would like my result to be in m not m^2 .

In this case the projection is given by

$$ \frac{v \cdot u}{u \cdot u} $$

which is unitless. It can be multiplied by $u$ to give the final result in length units again.

[pjabardo]

If you think of vectors as displacement vectors, they naturally should have length dimension. But direction vectors are very important and useful. v⃗⋅u⃗ makes perfect sense as a projection if u⃗ is a direction vector and the projection will have dimension of length.

I can see that Vec could be limited to displacement. The users of Meshes.jl that has some application in some specific area could use some other package to represent more general vectors. But, still, when doing geometric manipulations, the code is filled with ustrips and other related functions which makes the code uglier.

There is something else that I should mention, in CFD, in particular when dealing with incompressible flows, the geometry is almost always dimensionless. In this case, the units are not much of an issue since a dimensionless length is actually length whose unit is problem dependent. As an example, if I'm dealing with the flow around a sphere, all my coordinates are divided by the diameter D of the sphere and dimensionless. But what I'm really doing is using a unit of length that has length D. It would be nice to be able to specify lengths using the unit D.

[juliohm]

Thank you. We are working on this additional flexibility by first relaxing the <:Quantity restriction to <:Number and adjusting algorithms. Later on we can decide if it makes sense to always ustrip units from Vec at construction time.

Em sex., 13 de set. de 2024, 13:06, Paulo Jabardo @.***> escreveu:

If you think of vectors as displacement vectors, they naturally should have length dimension. But direction vectors are very important and useful. v⃗⋅u⃗ makes perfect sense as a projection if u⃗ is a direction vector and the projection will have dimension of length.

I can see that Vec could be limited to displacement. The users of Meshes.jl that has some application in some specific area could use some other package to represent more general vectors. But, still, when doing geometric manipulations, the code is filled with ustrips and other related functions which makes the code uglier.

There is something else that I should mention, in CFD, in particular when dealing with incompressible flows, the geometry is almost always dimensionless. In this case, the units are not much of an issue since a dimensionless length is actually length whose unit is problem dependent. As an example, if I'm dealing with the flow around a sphere, all my coordinates are divided by the diameter D of the sphere and dimensionless. But what I'm really doing is using a unit of length that has length D. It would be nice to be able to specify lengths using the unit D.

— Reply to this email directly, view it on GitHub https://github.com/JuliaGeometry/Meshes.jl/issues/996#issuecomment-2349297864, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAZQW3JFKRSTVXUOFOOLZCTZWMEQ7AVCNFSM6AAAAABMI64L5OVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDGNBZGI4TOOBWGQ . You are receiving this because you commented.Message ID: @.***>

[juliohm]

@pjabardo please consider the following 3 definitions of normal using the internal unormalize and ucross helpers:

julia> u = Vec(1, 0, 0)
Vec(1.0 m, 0.0 m, 0.0 m)

julia> v = Vec(0, 1, 0)
Vec(0.0 m, 1.0 m, 0.0 m)

julia> normal1(u::Vec, v::Vec) = unormalize(ucross(u, v))
normal1 (generic function with 1 method)

julia> normal2(u::Vec, v::Vec) = unormalize(u × v)
normal2 (generic function with 1 method)

julia> normal3(u::Vec, v::Vec) = normalize(u × v)
normal3 (generic function with 1 method)

julia> normal1(u, v)
Vec(0.0 m, 0.0 m, 1.0 m)

julia> normal2(u, v)
Vec(0.0 m^2, 0.0 m^2, 1.0 m^2)

julia> normal3(u, v)
Vec(0.0, 0.0, 1.0)

We adopted normal1 in Meshes.jl, your use case calls for normal2, and we believe that some people may prefer normal3. Is that correct? These definitions can be useful in different contexts.

We are tempted to stick to normal1 in this package, assuming that end-users can define their own variant in downstream applications. The main benefit of normal1 compared to the other variants is that it integrates well with other geometric algorithms that assume length units as discussed above. What do you think?

[pjabardo]

I think this solves most issues since Vec is no longer limited to length unit. And I believe that sticking to normal1 is reasonable since that is the current behavior and this is not something that causes problems for everyday use of the package. The possibility of easily defining other variants will certainly make things easier.

[juliohm]

Thank you @pjabardo. Closing the issue.