Skip to content
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.

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

Added support to integrate scalar/vector fields over objects of geometry module #19650

Merged
merged 6 commits into from Jul 3, 2020
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
@@ -14,7 +14,7 @@
from sympy.vector.orienters import (AxisOrienter, BodyOrienter,
SpaceOrienter, QuaternionOrienter)
from sympy.vector.operators import Gradient, Divergence, Curl, Laplacian, gradient, curl, divergence
from sympy.vector.parametricregion import ParametricRegion
from sympy.vector.parametricregion import (ParametricRegion, parametric_region_list)
from sympy.vector.integrals import (ParametricIntegral, vector_integrate)

__all__ = [
@@ -40,5 +40,5 @@
'Gradient', 'Divergence', 'Curl', 'Laplacian', 'gradient', 'curl',
'divergence',

'ParametricRegion', 'ParametricIntegral', 'vector_integrate',
'ParametricRegion', 'parametric_region_list', 'ParametricIntegral', 'vector_integrate',
]
@@ -1,11 +1,12 @@
from sympy.core.basic import Basic
from sympy import simplify
from sympy import S, simplify
from sympy.core import Basic, diff
from sympy.matrices import Matrix
from sympy.vector import CoordSys3D, Vector, ParametricRegion
from sympy.vector import CoordSys3D, Vector, ParametricRegion, parametric_region_list
from sympy.vector.operators import _get_coord_sys_from_expr
from sympy.integrals import Integral, integrate
from sympy.core.function import diff
from sympy.utilities.iterables import topological_sort, default_sort_key
from sympy.geometry.entity import GeometryEntity


class ParametricIntegral(Basic):
"""
@@ -32,7 +33,7 @@ class ParametricIntegral(Basic):
8*pi

>>> ParametricIntegral(C.j + C.k, ParametricRegion((r*cos(theta), r*sin(theta)), r, theta))
ParametricIntegral(C.j + C.k, ParametricRegion((r*cos(theta), r*sin(theta)), r, theta))
0

"""

@@ -48,7 +49,7 @@ def __new__(cls, field, parametricregion):
coord_sys = next(iter(coord_set))

if parametricregion.dimensions == 0:
return super().__new__(cls, field, parametricregion)
return S.Zero

base_vectors = coord_sys.base_vectors()
base_scalars = coord_sys.base_scalars()
@@ -135,6 +136,7 @@ def field(self):
def parametricregion(self):
return self.args[1]


def vector_integrate(field, *region):
"""
Compute the integral of a vector/scalar field
@@ -150,10 +152,31 @@ def vector_integrate(field, *region):
>>> vector_integrate(C.x*C.i, region)
12

Integrals over special regions can also be calculated using geometry module.
>>> from sympy.geometry import Point, Circle, Triangle
>>> c = Circle(Point(0, 2), 5)
>>> vector_integrate(C.x**2 + C.y**2, c)
290*pi
>>> triangle = Triangle(Point(-2, 3), Point(2, 3), Point(0, 5))
>>> vector_integrate(3*C.x**2*C.y*C.i + C.j, triangle)
-8

>>> vector_integrate(12*C.y**3, (C.y, 1, 3))
240
>>> vector_integrate(C.x**2*C.z, C.x)
C.x**3*C.z/3

"""
if len(region) == 1:
if isinstance(region[0], ParametricRegion):
return ParametricIntegral(field, region[0])

if isinstance(region[0], GeometryEntity):
regions_list = parametric_region_list(region[0])

result = 0
friyaz marked this conversation as resolved.
Show resolved Hide resolved
for reg in regions_list:
result += vector_integrate(field, reg)
return result

return integrate(field, *region)
@@ -1,5 +1,10 @@
from sympy.core.basic import Basic
from sympy.core.containers import Tuple
from functools import singledispatch
from sympy import pi
from sympy.core import Basic, Tuple
from sympy.core.symbol import _symbol
from sympy.solvers import solve
from sympy.geometry import Point, Segment, Curve, Ellipse, Polygon


class ParametricRegion(Basic):
"""
@@ -82,3 +87,84 @@ def parameters(self):
@property
def dimensions(self):
return len(self.limits)


@singledispatch
def parametric_region_list(reg):
"""
Returns a list of ParametricRegion objects representing the geometric region.

Examples
========

>>> from sympy.abc import t
>>> from sympy.vector import parametric_region_list
>>> from sympy.geometry import Point, Curve, Ellipse, Segment, Polygon

>>> p = Point(2, 5)
>>> parametric_region_list(p)
[ParametricRegion((2, 5))]

>>> c = Curve((t**3, 4*t), (t, -3, 4))
>>> parametric_region_list(c)
[ParametricRegion((t**3, 4*t), (t, -3, 4))]

>>> e = Ellipse(Point(1, 3), 2, 3)
>>> parametric_region_list(e)
[ParametricRegion((2*cos(t) + 1, 3*sin(t) + 3), (t, 0, 2*pi))]

>>> s = Segment(Point(1, 3), Point(2, 6))
>>> parametric_region_list(s)
[ParametricRegion((t + 1, 3*t + 3), (t, 0, 1))]

>>> p1, p2, p3, p4 = [(0, 1), (2, -3), (5, 3), (-2, 3)]
>>> poly = Polygon(p1, p2, p3, p4)
>>> parametric_region_list(poly)
[ParametricRegion((2*t, 1 - 4*t), (t, 0, 1)), ParametricRegion((3*t + 2, 6*t - 3), (t, 0, 1)),\
ParametricRegion((5 - 7*t, 3), (t, 0, 1)), ParametricRegion((2*t - 2, 3 - 2*t), (t, 0, 1))]

"""
raise ValueError("SymPy cannot determine parametric representation of the region.")


@parametric_region_list.register(Point)
def _(obj):
return [ParametricRegion(obj.args)]


@parametric_region_list.register(Curve)
def _(obj):
definition = obj.arbitrary_point(obj.parameter).args
bounds = obj.limits
return [ParametricRegion(definition, bounds)]


@parametric_region_list.register(Ellipse)
def _(obj, parameter='t'):
definition = obj.arbitrary_point(parameter).args
t = _symbol(parameter, real=True)
bounds = (t, 0, 2*pi)
return [ParametricRegion(definition, bounds)]


@parametric_region_list.register(Segment)
def _(obj, parameter='t'):
t = _symbol(parameter, real=True)
definition = obj.arbitrary_point(t).args

for i in range(0, 3):
lower_bound = solve(definition[i] - obj.points[0].args[i], t)
upper_bound = solve(definition[i] - obj.points[1].args[i], t)

if len(lower_bound) == 1 and len(upper_bound) == 1:
bounds = t, lower_bound[0], upper_bound[0]
break

definition_tuple = obj.arbitrary_point(parameter).args
return [ParametricRegion(definition_tuple, bounds)]


@parametric_region_list.register(Polygon)
def _(obj, parameter='t'):
l = [parametric_region_list(side, parameter)[0] for side in obj.sides]
return l
@@ -1,12 +1,14 @@
from sympy import sin, cos, pi, S, sqrt
from sympy.testing.pytest import raises
from sympy.vector.coordsysrect import CoordSys3D
from sympy.vector.integrals import ParametricIntegral
from sympy.vector.integrals import ParametricIntegral, vector_integrate
from sympy.vector.parametricregion import ParametricRegion
from sympy.abc import x, y, z, u, v, r, t, theta, phi
from sympy.geometry import Point, Segment, Curve, Circle, Polygon, Plane

C = CoordSys3D('C')

def test_lineintegrals():
def test_parametric_lineintegrals():
halfcircle = ParametricRegion((4*cos(theta), 4*sin(theta)), (theta, -pi/2, pi/2))
assert ParametricIntegral(C.x*C.y**4, halfcircle) == S(8192)/5

@@ -24,7 +26,7 @@ def test_lineintegrals():
field2 = C.y*C.i + C.z*C.j + C.z*C.k
assert ParametricIntegral(field2, ParametricRegion((cos(t), sin(t), t**2), (t, 0, pi))) == -5*pi/2 + pi**4/2

def test_surfaceintegrals():
def test_parametric_surfaceintegrals():

semisphere = ParametricRegion((2*sin(phi)*cos(theta), 2*sin(phi)*sin(theta), 2*cos(phi)),\
(theta, 0, 2*pi), (phi, 0, pi/2))
@@ -41,7 +43,7 @@ def test_surfaceintegrals():
assert ParametricIntegral(-15.6*C.y*C.k, triangle1) == ParametricIntegral(-15.6*C.y*C.k, triangle2)
assert ParametricIntegral(C.z, triangle1) == 10*C.z

def test_volumeintegrals():
def test_parametric_volumeintegrals():

cube = ParametricRegion((x, y, z), (x, 0, 1), (y, 0, 1), (z, 0, 1))
assert ParametricIntegral(1, cube) == 1
@@ -58,3 +60,35 @@ def test_volumeintegrals():
assert ParametricIntegral(C.x*C.i + C.j - 100*C.k, region_under_plane1) == \
ParametricIntegral(C.x*C.i + C.j - 100*C.k, region_under_plane2)
assert ParametricIntegral(2*C.x, region_under_plane2) == -9

def test_vector_integrate():
halfdisc = ParametricRegion((r*cos(theta), r* sin(theta)), (r, -2, 2), (theta, 0, pi))
assert vector_integrate(C.x**2, halfdisc) == 4*pi
vector_integrate(C.x, ParametricRegion((t, t**2), (t, 2, 3))) == -17*sqrt(17)/12 + 37*sqrt(37)/12

assert vector_integrate(C.y**3*C.z, (C.x, 0, 3), (C.y, -1, 4)) == 765*C.z/4

s1 = Segment(Point(0, 0), Point(0, 1))
assert vector_integrate(-15*C.y, s1) == S(-15)/2
s2 = Segment(Point(4, 3, 9), Point(1, 1, 7))
assert vector_integrate(C.y*C.i, s2) == -6

curve = Curve((sin(t), cos(t)), (t, 0, 2))
assert vector_integrate(5*C.z, curve) == 10*C.z

c1 = Circle(Point(2, 3), 6)
assert vector_integrate(C.x*C.y, c1) == 72*pi
c2 = Circle(Point(0, 0), Point(1, 1), Point(1, 0))
assert vector_integrate(1, c2) == c2.circumference

triangle = Polygon((0, 0), (1, 0), (1, 1))
assert vector_integrate(C.x*C.i - 14*C.y*C.j, triangle) == 0
p1, p2, p3, p4 = [(0, 0), (1, 0), (5, 1), (0, 1)]
poly = Polygon(p1, p2, p3, p4)
assert vector_integrate(-23*C.z, poly) == -161*C.z - 23*sqrt(17)*C.z

point = Point(2, 3)
assert vector_integrate(C.i*C.y - C.z, point) == ParametricIntegral(C.y*C.i, ParametricRegion((2, 3)))

pl = Plane(Point(1, 1, 1), Point(2, 3, 4), Point(2, 2, 2))
raises(ValueError, lambda: vector_integrate(C.x*C.z*C.i + C.k, pl))
@@ -1,12 +1,14 @@
from sympy import sin, cos, pi
from sympy.vector.coordsysrect import CoordSys3D
from sympy.vector.parametricregion import ParametricRegion
from sympy.vector.parametricregion import ParametricRegion, parametric_region_list
from sympy.geometry import Point, Segment, Curve, Ellipse, Line, Parabola, Polygon
from sympy.testing.pytest import raises
from sympy.abc import a, b, r, t, x, y, z, theta, phi


C = CoordSys3D('C')

def test_parametricregion():
def test_ParametricRegion():

point = ParametricRegion((3, 4))
assert point.definition == (3, 4)
@@ -42,7 +44,7 @@ def test_parametricregion():
assert circle.definition == (r*cos(theta), r*sin(theta))
assert circle.dimensions == 1

halfdisc = ParametricRegion((r*cos(theta), r* sin(theta)), (r, -2, 2), (theta, 0, pi))
halfdisc = ParametricRegion((r*cos(theta), r*sin(theta)), (r, -2, 2), (theta, 0, pi))
assert halfdisc.definition == (r*cos(theta), r*sin(theta))
assert halfdisc.parameters == (r, theta)
assert halfdisc.limits == {r: (-2, 2), theta: (0, pi)}
@@ -65,3 +67,30 @@ def test_parametricregion():

raises(ValueError, lambda: ParametricRegion((a*t**2, 2*a*t), (a, -2)))
raises(ValueError, lambda: ParametricRegion((a, b), (a**2, sin(b)), (a, 2, 4, 6)))


def test_parametric_region_list():

point = Point(-5, 12)
assert parametric_region_list(point) == [ParametricRegion((-5, 12))]

e = Ellipse(Point(2, 8), 2, 6)
assert parametric_region_list(e, t) == [ParametricRegion((2*cos(t) + 2, 6*sin(t) + 8), (t, 0, 2*pi))]

c = Curve((t, t**3), (t, 5, 3))
assert parametric_region_list(c) == [ParametricRegion((t, t**3), (t, 5, 3))]

s = Segment(Point(2, 11, -6), Point(0, 2, 5))
assert parametric_region_list(s, t) == [ParametricRegion((2 - 2*t, 11 - 9*t, 11*t - 6), (t, 0, 1))]
s1 = Segment(Point(0, 0), (1, 0))
assert parametric_region_list(s1, t) == [ParametricRegion((t, 0), (t, 0, 1))]
s2 = Segment(Point(1, 2, 3), Point(1, 2, 5))
assert parametric_region_list(s2, t) == [ParametricRegion((1, 2, 2*t + 3), (t, 0, 1))]
s3 = Segment(Point(12, 56), Point(12, 56))
assert parametric_region_list(s3) == [ParametricRegion((12, 56))]

poly = Polygon((1,3), (-3, 8), (2, 4))
assert parametric_region_list(poly, t) == [ParametricRegion((1 - 4*t, 5*t + 3), (t, 0, 1)), ParametricRegion((5*t - 3, 8 - 4*t), (t, 0, 1)), ParametricRegion((2 - t, 4 - t), (t, 0, 1))]

p1 = Parabola(Point(0, 0), Line(Point(5, 8), Point(7,8)))
raises(ValueError, lambda: parametric_region_list(p1))