Triangle is not generated by Point3D method

[abhishek-kuma]

tri=Triangle(Point3D(1,0,0),Point3D(0,1,0),Point3D(0,0,1)) while running this code it shows the following error

[smichr]

SymPy does not work with 3D objects other than Point, Plane and LinearEntity. Any such objects can be rotated into the x-y plane and queried there for any properties of interest.

>> from https://github.com/sympy/sympy/issues/24936 import RotationMatrix
>>> from sympy import *
>>> pts = [(0, 0, 0), (1, 2, 3), (4, 4, 5)]
>>> p = Plane(*pts)
>>> eqs = list(RotationMatrix(x=x,y=y)*Matrix(p.normal_vector))
>>> sol = solve(eqs[:2])[0]  # want normal to be parallel to Z axis so x and y components are 0
>>> r = RotationMatrix(x=sol[x],y=sol[y])
>>> t = Triangle(*[list(r*Matrix(i))[:2] for i in pts])
>>> t.area  # negative because of orientation of points
-sqrt(69)/2

put centroid into plane by doing reverse rotation

>>> _r = RotationMatrix(y=-sol[y],x=-sol[x])
>>> t.centroid
Point2D(sqrt(4485)/39, 16*sqrt(65)/39)
>>> c = Point(_r*Matrix(Point3D(_)))
>>> c in p
True
>>> c
Point3D(5/3, 2, 8/3)

See that the circumcenter satisfies requirements:

>>> cc = Point(_r*Matrix(Point3D(t.circumcenter)))
>>> len(set([cc.distance(i) for i in pts])) == 1
True

Putting this all together to return a helper function for dealing with any xy-point of interest might be like this:

def zPolygon(*pts):
    """return co-planar points in Polygon and function
    that will rotate any point in the xy-plane into the plane
    containing the pts as originally passed

    Examples
    ========

    >>> pts = [(0, 0, 0), (1, 2, 3), (4, 4, 5)]
    >>> t, f = zPolygon(*pts) # aTriangle is returned from 3 pts
    >>> f(t.centroid)
    Point3D(5/3, 2, 8/3)
    """
    assert not all(len(i) == 2 for i in pts), 'use Polygon for 2D points'
    assert len(pts) > 2, 'Point and Lines can be 3D'
    p = Plane(*pts[:3])
    assert all(i in p for i in pts[3:])
    eqs = list(RotationMatrix(x=x,y=y)*Matrix(p.normal_vector))
    # normal to be parallel to z axis: x and y components are 0
    sol = solve(eqs[:2])[0]
    r = RotationMatrix(x=sol[x],y=sol[y])
    g = Polygon(*[list(r*Matrix(i))[:2] for i in pts])
    _r = RotationMatrix(y=-sol[y],x=-sol[x])
    return g, lambda pt: Point(*_r*Matrix(list(pt)+[0]*(3-len(pt))))

cf #24936