questions about orientation of ParametricRegions

[Joe Heafner]

Hello everyone.

I teach undergraduate physics and have used Python, specifically VPython, for going on 20 years although I am by no means an expert. I am interested in incorporating computation (with Python using Jupyter Notebook) into introductory calculus-based mechanics and electromagnetic theory courses. I have an MS in physics and have taught for 30 years. I am a Python and LaTeX evangelist.

I am relatively new to SymPy and am amazed at its power. I want to bring this power into the introductory undergraduate physics curriculum.

Now that you know a bit about me, I have some questions about orienting surfaces defined with ParametricRegion. After a lot of online searching, I found this discussion

https://github.com/sympy/sympy/issues/19320

which doesn't seem to address my questions so was directed here by the SymPy Twitter feed.

Here is my specific situation. Consider the integral of a vector field over a surface to calculus flux.

from sympy import *

from sympy.vector import *

C = CoordSys3D('C')

a, b, c = symbols('a b c', positive=True)

x, y, z, F = symbols('x y z F')

F = 3*(-C.k)  # F points in the -z direction

Consider a square with vertices (0,0),(5,0),(5,5),(0,5). Treat the square as a Type 1 region, doing the y integral first.

square = ParametricRegion((x,y),(,y,0,5),(x,0,5))

vector_integrate(F, square)

I get 75 as the answer, which means that square must have an orientation in the -z direction. I can justify that by noting that the order in which I defined the parameters implies such an orientation ( yhat \times xhat = -zhat)

Now if I treat the square as a Type 2 region, doing the x integral first, I get a different answer.

square = ParametricRegion((x,y),(x,0,5),(y,0,5))

vector_integrate(F, square)

I get -75 as the answer, which means that square must have an orientation in the +z direction. I can also justify this because note that I reversed the order of the parameter definitions (xhat \times yhat = +zhat). Reversing the definitions reverses the surface's orientation. No problem.

Now here's where my trouble begins. Consider now a triangular region with vertices (0,0),(5,0),(0,5). Treat the triangle as a Type 1 region, doing the y integral first.

triangle = ParametricRegion((x,y),(x,0,5),(y,0,5-x))

vector_integrate(F, triangle)

I get 75/2, which means triangle must have an orientation in the -z direction. This does NOT seem to follow the convention from the previous example since I specified the parameter definitions in the order that I thought would give triangle an orientation in the +z direction but apparently that is wrong. Unlike with square, reversing the order of the parameter definitions doesn't change the orientation.

Now if I treat the triangle as a Type 2 region, doing the x integral first I get a different answer.

triangle = ParametricRegion((x,y),(x,0,5-y),(y,0,5))

vector_integrate(F, triangle)

I get -75/2, which means triangle must have an orientation in the +z direction, which, again, doesn't match my intuition from the square example. Also once again, reversing the order of the parameter definitions doesn't change the orientation.

With square, I can EITHER change the order of the parameters themselves (e.g. (x,y) vs. (y,x)) OR change the order of the parameter definitions (e.g. (x,0,5),(y,0,5) vs. (y,0,5),(x,0,5)). 

With triangle, it seems that changing the order of the parameters themselves (e.g. (x,y) vs. (y,x)) DOES change the surface's orientation but changing the order of the parameter definitions DOES NOT. If I keep the order of the parameters the same (e.g. always (x,y)) then the only way to change triangle's orientation is to treat it as either a Type 1 or Type 2 region.

Here is my ultimate question. What rule(s) is/are used internally to orient a surface and why does the rule(s) appear to be different for different surfaces as I have presented here? How/why does a Type 1 region's orientation differ from a Type 2 region's orientation? They seem to collapse to behavior with square when the integral is the same regardless of whether it's written for a Type 1 or Type 2 region. I have never seen that addressed in calculus UNLESS different forms are used. Am I missing something obvious? Am I missing something subtle? I can't seem to find a consistent SymPy rule for surface orientation. Does one exist?

I look forward to responses!

Joe

[Oscar Benjamin]

Hi Joe,

The vector integration functionality is relatively new. It was added
last year I think in a GSOC project. It's possible that this is a bug
or a corner case that wasn't fully considered.

Have you looked into the code at all?

Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2a688b87-2616-49ab-b62e-42b417b48cb5n%40googlegroups.com.

[Faisal Riyaz]

Hey Joe,

Thanks for pointing this out.

The sign of the result depends on the direction of the normal vector. SymPy calculates normal vector using:

where r is the position vector in terms of parameters.

Different signs for these cases are due to the determination of u and v.

Generally, Parameters u and v are assumed from the order of the bounds

used at the creation of ParametricRegion object.

But when the limits of parameters are dependent on each other such as in the case of

the triangle.

ParametricRegion((x,y),(x,0,5),(y,0,5-x)) # Bounds of y are dependent on x

SymPy uses topological sort to determine the order of integration.

This is needed because we would get an incorrect result

if we perform integration with respect to x first in this case.

You can have a look at the implementation here:

https://github.com/sympy/sympy/blob/master/sympy/vector/integrals.py#L81

Your concern is valid. The results should be consistent.

As a user, feel free to suggest how SymPy should determine the direction

of normal vector.  Personally, I think the order of bounds when defining

parametric regions should define the direction. Also, we would need to document this better.

Once decided on the behavior, I will make the necessary changes to fix this.

Best,

Faisal

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRcbYHAcncvcMtdURGQ134RD9XevcyvBDhGFumeaTP2wg%40mail.gmail.com.

[Joe Heafner]

Hello Oscar.

Thank you very much for the reply. Yes, I am indeed looking at the source. I think Faisal identified the problem in his response and I look forward to being a part of the discussion to address it. The SymPy community seems very open, and I appreciate that.

Joe Heafner
Sent from my iPad

> On Jul 19, 2021, at 05:37, Oscar Benjamin <oscar.j....@gmail.com> wrote:
>
> Hi Joe,

> You received this message because you are subscribed to a topic in the Google Groups "sympy" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/sympy/zy9gQH9gMvQ/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRcbYHAcncvcMtdURGQ134RD9XevcyvBDhGFumeaTP2wg%40mail.gmail.com.

[Joe Heafner]

Hello Faisal.

Thank you so much for the helpful reply. Your explanation finally got me to, I think, understand precisely how u and v are defined, which was previously a mystery to me. I will verify my understanding with some specific examples to make sure I understand.

I also noticed something over the past few days. The main difference between a type 1 and type 2 region is which coordinate serves as the parameter. I realized that when I use y as the parameter (i.e. treat the triangle as a type 2 region), I previously would use

ParametricRegion((x,y),(x,0,5-y),(y,0,5))

 but I realized that 

ParametricRegion((x,y),(x,0,5-y),(y,5,0)) 

made more sense given that x starts at 0, and of course this changes the sign of the resulting integral. I believe the code expects the lower limit before the upper limit, which corresponds to the first expression I was using. 

I’m not sure what convention would be “best” but I’m happy to be part of that conversation, so I will continue to give it more thought.

Thank you again!

Joe Heafner

Sent from my iPad

On Jul 19, 2021, at 08:30, Faisal Riyaz <faisalr...@gmail.com> wrote:



Hey Joe,

Thanks for pointing this out.

The sign of the result depends on the direction of the normal vector. SymPy calculates normal vector using:

[Joe Heafner]

Faisal,

After more experimentation and thought, I agree that the order of the bounds is probably best. To make sure I completely understand, let me present an example.

For ParametricRegion((x,y),(x,0,5),(y,0,5)) dr/du would point in the +x direction, dr/dv would point in the +y direction so the orientation would be in the direction of +x \times +y or the +z direction. Is that equivalent to what you’re saying? If so, then the order in which x and y appear in the definition (i.e. (x,y)) should be the same as the order in which x and y appear when their bounds are defined. In other words, the two definitions

ParametricRegion((x,y),(x,0,5),(y,0,5))

and

ParametricRegion((y,x),(y,0,5),(x,0,5))

should result in the same orientation for this rectangular region, as they currently do. Is this correct? 

However, with my triangular region example, by the new intended convention the two definitions

ParametricRegion((x,y),(x,0,5-y),(y,5,0))

and 

ParametricRegion((y,x),(y,5,5-x),(x,0,5))

should also give the same orientation, which they currently DO NOT. Is this correct?

I hope I have understood correctly. 

[Joe Heafner]

Oh dear, I think I made an error in my last example. I meant to say

ParametricRegion((x,y),(x,0,5),(y,5,5-x)) and ParametricRegion((y,x),(y,5,5-x),(x,0,5))

should give the same orientation, which they currently do not.

And 

ParametricRegion((x,y),(x,0,5-y),(y,5,0)) and ParametricRegion((y,x),(y,5,0),(x,0,5-y))

should give the same orientation, which they currently do not.

This is so easily confused. Again, I hope my reasoning is clear and correct.

[Aaron Meurer]

Was there ever an issue or pull request opened for this?

Aaron Meurer

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/3f50a9ff-35d2-417c-bd01-5edb8e32885an%40googlegroups.com.