diff --git a/cpp/modmesh/universe/polygon.hpp b/cpp/modmesh/universe/polygon.hpp
index e144e8bc..4a3d16c7 100644
--- a/cpp/modmesh/universe/polygon.hpp
+++ b/cpp/modmesh/universe/polygon.hpp
@@ -63,8 +63,1099 @@ union Triangle3dData
     Triangle3dNamed<T> f;
 }; /* end union Triangle3dData */
 
+template <typename T>
+struct Trapezoid3dNamed
+{
+    T x0, x1, x2, x3, y0, y1, y2, y3, z0, z1, z2, z3;
+}; /* end struct Trapezoid3dNamed */
+
+template <typename T>
+union Trapezoid3dData
+{
+    T v[12];
+    Trapezoid3dNamed<T> f;
+}; /* end union Trapezoid3dData */
+
 } /* end namespace detail */
 
+/**
+ * Trapezoid in three-dimensional space.
+ *
+ * @tparam T floating-point type
+ */
+template <typename T>
+class Trapezoid3d
+    : public NumberBase<int32_t, T>
+{
+
+public:
+
+    static_assert(std::is_arithmetic_v<T>, "T in Trapezoid3d<T> must be arithmetic type");
+
+    using point_type = Point3d<T>;
+    using value_type = typename point_type::value_type;
+
+    Trapezoid3d(point_type const & p0, point_type const & p1, point_type const & p2, point_type const & p3)
+        : m_data{p0.x(), p1.x(), p2.x(), p3.x(), p0.y(), p1.y(), p2.y(), p3.y(), p0.z(), p1.z(), p2.z(), p3.z()}
+    {
+    }
+
+    Trapezoid3d() = default;
+    Trapezoid3d(Trapezoid3d const &) = default;
+    Trapezoid3d & operator=(Trapezoid3d const &) = default;
+    Trapezoid3d(Trapezoid3d &&) = default;
+    Trapezoid3d & operator=(Trapezoid3d &&) = default;
+    ~Trapezoid3d() = default;
+
+    point_type p0() const { return point_type(m_data.f.x0, m_data.f.y0, m_data.f.z0); }
+    void set_p0(point_type const & p)
+    {
+        m_data.f.x0 = p.x();
+        m_data.f.y0 = p.y();
+        m_data.f.z0 = p.z();
+    }
+    point_type p1() const { return point_type(m_data.f.x1, m_data.f.y1, m_data.f.z1); }
+    void set_p1(point_type const & p)
+    {
+        m_data.f.x1 = p.x();
+        m_data.f.y1 = p.y();
+        m_data.f.z1 = p.z();
+    }
+    point_type p2() const { return point_type(m_data.f.x2, m_data.f.y2, m_data.f.z2); }
+    void set_p2(point_type const & p)
+    {
+        m_data.f.x2 = p.x();
+        m_data.f.y2 = p.y();
+        m_data.f.z2 = p.z();
+    }
+    point_type p3() const { return point_type(m_data.f.x3, m_data.f.y3, m_data.f.z3); }
+    void set_p3(point_type const & p)
+    {
+        m_data.f.x3 = p.x();
+        m_data.f.y3 = p.y();
+        m_data.f.z3 = p.z();
+    }
+
+    value_type x0() const { return m_data.f.x0; }
+    value_type & x0() { return m_data.f.x0; }
+    void set_x0(value_type v) { m_data.f.x0 = v; }
+
+    value_type y0() const { return m_data.f.y0; }
+    value_type & y0() { return m_data.f.y0; }
+    void set_y0(value_type v) { m_data.f.y0 = v; }
+
+    value_type z0() const { return m_data.f.z0; }
+    value_type & z0() { return m_data.f.z0; }
+    void set_z0(value_type v) { m_data.f.z0 = v; }
+
+    value_type x1() const { return m_data.f.x1; }
+    value_type & x1() { return m_data.f.x1; }
+    void set_x1(value_type v) { m_data.f.x1 = v; }
+
+    value_type y1() const { return m_data.f.y1; }
+    value_type & y1() { return m_data.f.y1; }
+    void set_y1(value_type v) { m_data.f.y1 = v; }
+
+    value_type z1() const { return m_data.f.z1; }
+    value_type & z1() { return m_data.f.z1; }
+    void set_z1(value_type v) { m_data.f.z1 = v; }
+
+    value_type x2() const { return m_data.f.x2; }
+    value_type & x2() { return m_data.f.x2; }
+    void set_x2(value_type v) { m_data.f.x2 = v; }
+
+    value_type y2() const { return m_data.f.y2; }
+    value_type & y2() { return m_data.f.y2; }
+    void set_y2(value_type v) { m_data.f.y2 = v; }
+
+    value_type z2() const { return m_data.f.z2; }
+    value_type & z2() { return m_data.f.z2; }
+    void set_z2(value_type v) { m_data.f.z2 = v; }
+
+    value_type x3() const { return m_data.f.x3; }
+    value_type & x3() { return m_data.f.x3; }
+    void set_x3(value_type v) { m_data.f.x3 = v; }
+
+    value_type y3() const { return m_data.f.y3; }
+    value_type & y3() { return m_data.f.y3; }
+    void set_y3(value_type v) { m_data.f.y3 = v; }
+
+    value_type z3() const { return m_data.f.z3; }
+    value_type & z3() { return m_data.f.z3; }
+    void set_z3(value_type v) { m_data.f.z3 = v; }
+
+    point_type operator[](size_t i) const { return point_type(m_data.v[i], m_data.v[i + 4], m_data.v[i + 8]); }
+
+    point_type at(size_t i) const
+    {
+        check_size(i, 4);
+        return (*this)[i];
+    }
+
+    bool operator==(Trapezoid3d const & other) const
+    {
+        return m_data.v[0] == other.m_data.v[0] &&
+               m_data.v[1] == other.m_data.v[1] &&
+               m_data.v[2] == other.m_data.v[2] &&
+               m_data.v[3] == other.m_data.v[3] &&
+               m_data.v[4] == other.m_data.v[4] &&
+               m_data.v[5] == other.m_data.v[5] &&
+               m_data.v[6] == other.m_data.v[6] &&
+               m_data.v[7] == other.m_data.v[7] &&
+               m_data.v[8] == other.m_data.v[8] &&
+               m_data.v[9] == other.m_data.v[9] &&
+               m_data.v[10] == other.m_data.v[10] &&
+               m_data.v[11] == other.m_data.v[11];
+    }
+
+    bool operator!=(Trapezoid3d const & other) const
+    {
+        return m_data.v[0] != other.m_data.v[0] ||
+               m_data.v[1] != other.m_data.v[1] ||
+               m_data.v[2] != other.m_data.v[2] ||
+               m_data.v[3] != other.m_data.v[3] ||
+               m_data.v[4] != other.m_data.v[4] ||
+               m_data.v[5] != other.m_data.v[5] ||
+               m_data.v[6] != other.m_data.v[6] ||
+               m_data.v[7] != other.m_data.v[7] ||
+               m_data.v[8] != other.m_data.v[8] ||
+               m_data.v[9] != other.m_data.v[9] ||
+               m_data.v[10] != other.m_data.v[10] ||
+               m_data.v[11] != other.m_data.v[11];
+    }
+
+    size_t size() const { return 4; }
+
+    void mirror_x()
+    {
+        m_data.f.y0 = -m_data.f.y0;
+        m_data.f.y1 = -m_data.f.y1;
+        m_data.f.y2 = -m_data.f.y2;
+        m_data.f.y3 = -m_data.f.y3;
+        m_data.f.z0 = -m_data.f.z0;
+        m_data.f.z1 = -m_data.f.z1;
+        m_data.f.z2 = -m_data.f.z2;
+        m_data.f.z3 = -m_data.f.z3;
+    }
+
+    void mirror_y()
+    {
+        m_data.f.x0 = -m_data.f.x0;
+        m_data.f.x1 = -m_data.f.x1;
+        m_data.f.x2 = -m_data.f.x2;
+        m_data.f.x3 = -m_data.f.x3;
+        m_data.f.z0 = -m_data.f.z0;
+        m_data.f.z1 = -m_data.f.z1;
+        m_data.f.z2 = -m_data.f.z2;
+        m_data.f.z3 = -m_data.f.z3;
+    }
+
+    void mirror_z()
+    {
+        m_data.f.x0 = -m_data.f.x0;
+        m_data.f.x1 = -m_data.f.x1;
+        m_data.f.x2 = -m_data.f.x2;
+        m_data.f.x3 = -m_data.f.x3;
+        m_data.f.y0 = -m_data.f.y0;
+        m_data.f.y1 = -m_data.f.y1;
+        m_data.f.y2 = -m_data.f.y2;
+        m_data.f.y3 = -m_data.f.y3;
+    }
+
+    void mirror(Axis axis)
+    {
+        switch (axis)
+        {
+        case Axis::X: mirror_x(); break;
+        case Axis::Y: mirror_y(); break;
+        case Axis::Z: mirror_z(); break;
+        default: throw std::invalid_argument("Trapezoid3d::mirror: invalid axis"); break;
+        }
+    }
+
+private:
+
+    void check_size(size_t i, size_t s) const
+    {
+        if (i >= s)
+        {
+            throw std::out_of_range(std::format("Trapezoid3d: i {} >= size {}", i, s));
+        }
+    }
+
+    detail::Trapezoid3dData<T> m_data;
+
+}; /* end class Trapezoid3d */
+
+using Trapezoid3dFp32 = Trapezoid3d<float>;
+using Trapezoid3dFp64 = Trapezoid3d<double>;
+
+template <typename T>
+class TrapezoidPad
+    : public NumberBase<int32_t, T>
+    , public std::enable_shared_from_this<TrapezoidPad<T>>
+{
+
+private:
+
+    struct ctor_passkey
+    {
+    };
+
+public:
+
+    using real_type = T;
+    using value_type = T;
+    using point_type = Point3d<T>;
+    using trapezoid_type = Trapezoid3d<T>;
+    using point_pad_type = PointPad<T>;
+
+    template <typename... Args>
+    static std::shared_ptr<TrapezoidPad<T>> construct(Args &&... args)
+    {
+        return std::make_shared<TrapezoidPad<T>>(std::forward<Args>(args)..., ctor_passkey());
+    }
+
+    TrapezoidPad(uint8_t ndim, ctor_passkey const &)
+        : m_p0(point_pad_type::construct(ndim))
+        , m_p1(point_pad_type::construct(ndim))
+        , m_p2(point_pad_type::construct(ndim))
+        , m_p3(point_pad_type::construct(ndim))
+    {
+    }
+
+    TrapezoidPad(uint8_t ndim, size_t nelem, ctor_passkey const &)
+        : m_p0(point_pad_type::construct(ndim, nelem))
+        , m_p1(point_pad_type::construct(ndim, nelem))
+        , m_p2(point_pad_type::construct(ndim, nelem))
+        , m_p3(point_pad_type::construct(ndim, nelem))
+    {
+    }
+
+    TrapezoidPad(
+        SimpleArray<T> const & x0,
+        SimpleArray<T> const & y0,
+        SimpleArray<T> const & x1,
+        SimpleArray<T> const & y1,
+        SimpleArray<T> const & x2,
+        SimpleArray<T> const & y2,
+        SimpleArray<T> const & x3,
+        SimpleArray<T> const & y3,
+        ctor_passkey const &)
+        : m_p0(point_pad_type::construct(x0, y0))
+        , m_p1(point_pad_type::construct(x1, y1))
+        , m_p2(point_pad_type::construct(x2, y2))
+        , m_p3(point_pad_type::construct(x3, y3))
+    {
+        check_constructor_point_size(*m_p0, *m_p1, *m_p2, *m_p3);
+    }
+
+    TrapezoidPad(
+        SimpleArray<T> const & x0,
+        SimpleArray<T> const & y0,
+        SimpleArray<T> const & z0,
+        SimpleArray<T> const & x1,
+        SimpleArray<T> const & y1,
+        SimpleArray<T> const & z1,
+        SimpleArray<T> const & x2,
+        SimpleArray<T> const & y2,
+        SimpleArray<T> const & z2,
+        SimpleArray<T> const & x3,
+        SimpleArray<T> const & y3,
+        SimpleArray<T> const & z3,
+        ctor_passkey const &)
+        : m_p0(point_pad_type::construct(x0, y0, z0))
+        , m_p1(point_pad_type::construct(x1, y1, z1))
+        , m_p2(point_pad_type::construct(x2, y2, z2))
+        , m_p3(point_pad_type::construct(x3, y3, z3))
+    {
+        check_constructor_point_size(*m_p0, *m_p1, *m_p2, *m_p3);
+    }
+
+    TrapezoidPad(
+        SimpleArray<T> & x0,
+        SimpleArray<T> & y0,
+        SimpleArray<T> & x1,
+        SimpleArray<T> & y1,
+        SimpleArray<T> & x2,
+        SimpleArray<T> & y2,
+        SimpleArray<T> & x3,
+        SimpleArray<T> & y3,
+        bool clone,
+        ctor_passkey const &)
+        : m_p0(point_pad_type::construct(x0, y0, clone))
+        , m_p1(point_pad_type::construct(x1, y1, clone))
+        , m_p2(point_pad_type::construct(x2, y2, clone))
+        , m_p3(point_pad_type::construct(x3, y3, clone))
+    {
+        check_constructor_point_size(*m_p0, *m_p1, *m_p2, *m_p3);
+    }
+
+    TrapezoidPad(
+        SimpleArray<T> & x0,
+        SimpleArray<T> & y0,
+        SimpleArray<T> & z0,
+        SimpleArray<T> & x1,
+        SimpleArray<T> & y1,
+        SimpleArray<T> & z1,
+        SimpleArray<T> & x2,
+        SimpleArray<T> & y2,
+        SimpleArray<T> & z2,
+        SimpleArray<T> & x3,
+        SimpleArray<T> & y3,
+        SimpleArray<T> & z3,
+        bool clone,
+        ctor_passkey const &)
+        : m_p0(point_pad_type::construct(x0, y0, z0, clone))
+        , m_p1(point_pad_type::construct(x1, y1, z1, clone))
+        , m_p2(point_pad_type::construct(x2, y2, z2, clone))
+        , m_p3(point_pad_type::construct(x3, y3, z3, clone))
+    {
+        check_constructor_point_size(*m_p0, *m_p1, *m_p2, *m_p3);
+    }
+
+    std::shared_ptr<TrapezoidPad<T>> clone()
+    {
+        if (ndim() == 2)
+        {
+            return TrapezoidPad<T>::construct(x0(), y0(), x1(), y1(), x2(), y2(), x3(), y3());
+        }
+        else
+        {
+            return TrapezoidPad<T>::construct(x0(), y0(), z0(), x1(), y1(), z1(), x2(), y2(), z2(), x3(), y3(), z3());
+        }
+    }
+
+    TrapezoidPad() = delete;
+    TrapezoidPad(TrapezoidPad const &) = delete;
+    TrapezoidPad(TrapezoidPad &&) = delete;
+    TrapezoidPad & operator=(TrapezoidPad const &) = delete;
+    TrapezoidPad & operator=(TrapezoidPad &&) = delete;
+
+    ~TrapezoidPad() = default;
+
+    void append(trapezoid_type const & t)
+    {
+        if (ndim() == 2)
+        {
+            m_p0->append(t.x0(), t.y0());
+            m_p1->append(t.x1(), t.y1());
+            m_p2->append(t.x2(), t.y2());
+            m_p3->append(t.x3(), t.y3());
+        }
+        else
+        {
+            m_p0->append(t.x0(), t.y0(), t.z0());
+            m_p1->append(t.x1(), t.y1(), t.z1());
+            m_p2->append(t.x2(), t.y2(), t.z2());
+            m_p3->append(t.x3(), t.y3(), t.z3());
+        }
+    }
+
+    void append(point_type const & p0, point_type const & p1, point_type const & p2, point_type const & p3)
+    {
+        if (ndim() == 2)
+        {
+            m_p0->append(p0.x(), p0.y());
+            m_p1->append(p1.x(), p1.y());
+            m_p2->append(p2.x(), p2.y());
+            m_p3->append(p3.x(), p3.y());
+        }
+        else
+        {
+            m_p0->append(p0.x(), p0.y(), p0.z());
+            m_p1->append(p1.x(), p1.y(), p1.z());
+            m_p2->append(p2.x(), p2.y(), p2.z());
+            m_p3->append(p3.x(), p3.y(), p3.z());
+        }
+    }
+
+    void append(T x0, T y0, T x1, T y1, T x2, T y2, T x3, T y3)
+    {
+        m_p0->append(x0, y0);
+        m_p1->append(x1, y1);
+        m_p2->append(x2, y2);
+        m_p3->append(x3, y3);
+    }
+
+    void append(T x0, T y0, T z0, T x1, T y1, T z1, T x2, T y2, T z2, T x3, T y3, T z3)
+    {
+        m_p0->append(x0, y0, z0);
+        m_p1->append(x1, y1, z1);
+        m_p2->append(x2, y2, z2);
+        m_p3->append(x3, y3, z3);
+    }
+
+    void extend_with(TrapezoidPad<T> const & other)
+    {
+        size_t const ntrap = other.size();
+        for (size_t i = 0; i < ntrap; ++i)
+        {
+            append(other.get(i));
+        }
+    }
+
+    uint8_t ndim() const { return m_p0->ndim(); }
+
+    size_t size() const { return m_p0->size(); }
+
+    SimpleArray<T> pack_array() const
+    {
+        using shape_type = typename SimpleArray<T>::shape_type;
+        SimpleArray<T> ret(shape_type{m_p0->size(), static_cast<size_t>(ndim() * 4)});
+        if (ndim() == 3)
+        {
+            for (size_t i = 0; i < m_p0->size(); ++i)
+            {
+                ret(i, 0) = m_p0->x(i);
+                ret(i, 1) = m_p0->y(i);
+                ret(i, 2) = m_p0->z(i);
+                ret(i, 3) = m_p1->x(i);
+                ret(i, 4) = m_p1->y(i);
+                ret(i, 5) = m_p1->z(i);
+                ret(i, 6) = m_p2->x(i);
+                ret(i, 7) = m_p2->y(i);
+                ret(i, 8) = m_p2->z(i);
+                ret(i, 9) = m_p3->x(i);
+                ret(i, 10) = m_p3->y(i);
+                ret(i, 11) = m_p3->z(i);
+            }
+        }
+        else
+        {
+            for (size_t i = 0; i < m_p0->size(); ++i)
+            {
+                ret(i, 0) = m_p0->x(i);
+                ret(i, 1) = m_p0->y(i);
+                ret(i, 2) = m_p1->x(i);
+                ret(i, 3) = m_p1->y(i);
+                ret(i, 4) = m_p2->x(i);
+                ret(i, 5) = m_p2->y(i);
+                ret(i, 6) = m_p3->x(i);
+                ret(i, 7) = m_p3->y(i);
+            }
+        }
+        return ret;
+    }
+
+    void expand(size_t length)
+    {
+        m_p0->expand(length);
+        m_p1->expand(length);
+        m_p2->expand(length);
+        m_p3->expand(length);
+    }
+
+    real_type x0_at(size_t i) const { return m_p0->x_at(i); }
+    real_type y0_at(size_t i) const { return m_p0->y_at(i); }
+    real_type z0_at(size_t i) const { return m_p0->z_at(i); }
+    real_type x1_at(size_t i) const { return m_p1->x_at(i); }
+    real_type y1_at(size_t i) const { return m_p1->y_at(i); }
+    real_type z1_at(size_t i) const { return m_p1->z_at(i); }
+    real_type x2_at(size_t i) const { return m_p2->x_at(i); }
+    real_type y2_at(size_t i) const { return m_p2->y_at(i); }
+    real_type z2_at(size_t i) const { return m_p2->z_at(i); }
+    real_type x3_at(size_t i) const { return m_p3->x_at(i); }
+    real_type y3_at(size_t i) const { return m_p3->y_at(i); }
+    real_type z3_at(size_t i) const { return m_p3->z_at(i); }
+    real_type & x0_at(size_t i) { return m_p0->x_at(i); }
+    real_type & y0_at(size_t i) { return m_p0->y_at(i); }
+    real_type & z0_at(size_t i) { return m_p0->z_at(i); }
+    real_type & x1_at(size_t i) { return m_p1->x_at(i); }
+    real_type & y1_at(size_t i) { return m_p1->y_at(i); }
+    real_type & z1_at(size_t i) { return m_p1->z_at(i); }
+    real_type & x2_at(size_t i) { return m_p2->x_at(i); }
+    real_type & y2_at(size_t i) { return m_p2->y_at(i); }
+    real_type & z2_at(size_t i) { return m_p2->z_at(i); }
+    real_type & x3_at(size_t i) { return m_p3->x_at(i); }
+    real_type & y3_at(size_t i) { return m_p3->y_at(i); }
+    real_type & z3_at(size_t i) { return m_p3->z_at(i); }
+
+    real_type x0(size_t i) const { return m_p0->x(i); }
+    real_type y0(size_t i) const { return m_p0->y(i); }
+    real_type z0(size_t i) const { return m_p0->z(i); }
+    real_type x1(size_t i) const { return m_p1->x(i); }
+    real_type y1(size_t i) const { return m_p1->y(i); }
+    real_type z1(size_t i) const { return m_p1->z(i); }
+    real_type x2(size_t i) const { return m_p2->x(i); }
+    real_type y2(size_t i) const { return m_p2->y(i); }
+    real_type z2(size_t i) const { return m_p2->z(i); }
+    real_type x3(size_t i) const { return m_p3->x(i); }
+    real_type y3(size_t i) const { return m_p3->y(i); }
+    real_type z3(size_t i) const { return m_p3->z(i); }
+    real_type & x0(size_t i) { return m_p0->x(i); }
+    real_type & y0(size_t i) { return m_p0->y(i); }
+    real_type & z0(size_t i) { return m_p0->z(i); }
+    real_type & x1(size_t i) { return m_p1->x(i); }
+    real_type & y1(size_t i) { return m_p1->y(i); }
+    real_type & z1(size_t i) { return m_p1->z(i); }
+    real_type & x2(size_t i) { return m_p2->x(i); }
+    real_type & y2(size_t i) { return m_p2->y(i); }
+    real_type & z2(size_t i) { return m_p2->z(i); }
+    real_type & x3(size_t i) { return m_p3->x(i); }
+    real_type & y3(size_t i) { return m_p3->y(i); }
+    real_type & z3(size_t i) { return m_p3->z(i); }
+
+    point_type p0_at(size_t i) const { return m_p0->get_at(i); }
+    point_type p1_at(size_t i) const { return m_p1->get_at(i); }
+    point_type p2_at(size_t i) const { return m_p2->get_at(i); }
+    point_type p3_at(size_t i) const { return m_p3->get_at(i); }
+    void set_p0_at(size_t i, point_type const & p) { m_p0->set_at(i, p); }
+    void set_p1_at(size_t i, point_type const & p) { m_p1->set_at(i, p); }
+    void set_p2_at(size_t i, point_type const & p) { m_p2->set_at(i, p); }
+    void set_p3_at(size_t i, point_type const & p) { m_p3->set_at(i, p); }
+
+    point_type p0(size_t i) const { return m_p0->get(i); }
+    point_type p1(size_t i) const { return m_p1->get(i); }
+    point_type p2(size_t i) const { return m_p2->get(i); }
+    point_type p3(size_t i) const { return m_p3->get(i); }
+    void set_p0(size_t i, point_type const & p) { m_p0->set(i, p); }
+    void set_p1(size_t i, point_type const & p) { m_p1->set(i, p); }
+    void set_p2(size_t i, point_type const & p) { m_p2->set(i, p); }
+    void set_p3(size_t i, point_type const & p) { m_p3->set(i, p); }
+
+    SimpleArray<value_type> x0() { return m_p0->x(); }
+    SimpleArray<value_type> y0() { return m_p0->y(); }
+    SimpleArray<value_type> z0() { return m_p0->z(); }
+    SimpleArray<value_type> x1() { return m_p1->x(); }
+    SimpleArray<value_type> y1() { return m_p1->y(); }
+    SimpleArray<value_type> z1() { return m_p1->z(); }
+    SimpleArray<value_type> x2() { return m_p2->x(); }
+    SimpleArray<value_type> y2() { return m_p2->y(); }
+    SimpleArray<value_type> z2() { return m_p2->z(); }
+    SimpleArray<value_type> x3() { return m_p3->x(); }
+    SimpleArray<value_type> y3() { return m_p3->y(); }
+    SimpleArray<value_type> z3() { return m_p3->z(); }
+
+    std::shared_ptr<point_pad_type> p0() { return m_p0; }
+    std::shared_ptr<point_pad_type> p1() { return m_p1; }
+    std::shared_ptr<point_pad_type> p2() { return m_p2; }
+    std::shared_ptr<point_pad_type> p3() { return m_p3; }
+
+    trapezoid_type get_at(size_t i) const
+    {
+        if (ndim() == 3)
+        {
+            return trapezoid_type(point_type(x0_at(i), y0_at(i), z0_at(i)),
+                                  point_type(x1_at(i), y1_at(i), z1_at(i)),
+                                  point_type(x2_at(i), y2_at(i), z2_at(i)),
+                                  point_type(x3_at(i), y3_at(i), z3_at(i)));
+        }
+        else
+        {
+            return trapezoid_type(point_type(x0_at(i), y0_at(i), 0.0),
+                                  point_type(x1_at(i), y1_at(i), 0.0),
+                                  point_type(x2_at(i), y2_at(i), 0.0),
+                                  point_type(x3_at(i), y3_at(i), 0.0));
+        }
+    }
+    void set_at(size_t i, trapezoid_type const & t)
+    {
+        x0_at(i) = t.x0();
+        y0_at(i) = t.y0();
+        x1_at(i) = t.x1();
+        y1_at(i) = t.y1();
+        x2_at(i) = t.x2();
+        y2_at(i) = t.y2();
+        x3_at(i) = t.x3();
+        y3_at(i) = t.y3();
+        if (ndim() == 3)
+        {
+            z0_at(i) = t.z0();
+            z1_at(i) = t.z1();
+            z2_at(i) = t.z2();
+            z3_at(i) = t.z3();
+        }
+    }
+    void set_at(size_t i, point_type const & p0, point_type const & p1, point_type const & p2, point_type const & p3)
+    {
+        x0_at(i) = p0.x();
+        y0_at(i) = p0.y();
+        x1_at(i) = p1.x();
+        y1_at(i) = p1.y();
+        x2_at(i) = p2.x();
+        y2_at(i) = p2.y();
+        x3_at(i) = p3.x();
+        y3_at(i) = p3.y();
+        if (ndim() == 3)
+        {
+            z0_at(i) = p0.z();
+            z1_at(i) = p1.z();
+            z2_at(i) = p2.z();
+            z3_at(i) = p3.z();
+        }
+    }
+    // clang-format off
+    void set_at(size_t i, value_type x0, value_type y0,
+                value_type x1, value_type y1,
+                value_type x2, value_type y2,
+                value_type x3, value_type y3)
+    // clang-format on
+    {
+        x0_at(i) = x0;
+        y0_at(i) = y0;
+        x1_at(i) = x1;
+        y1_at(i) = y1;
+        x2_at(i) = x2;
+        y2_at(i) = y2;
+        x3_at(i) = x3;
+        y3_at(i) = y3;
+    }
+    // clang-format off
+    void set_at(size_t i, value_type x0, value_type y0, value_type z0,
+                value_type x1, value_type y1, value_type z1,
+                value_type x2, value_type y2, value_type z2,
+                value_type x3, value_type y3, value_type z3)
+    // clang-format on
+    {
+        x0_at(i) = x0;
+        y0_at(i) = y0;
+        x1_at(i) = x1;
+        y1_at(i) = y1;
+        x2_at(i) = x2;
+        y2_at(i) = y2;
+        x3_at(i) = x3;
+        y3_at(i) = y3;
+        if (ndim() == 3)
+        {
+            z0_at(i) = z0;
+            z1_at(i) = z1;
+            z2_at(i) = z2;
+            z3_at(i) = z3;
+        }
+    }
+
+    trapezoid_type get(size_t i) const
+    {
+        if (ndim() == 3)
+        {
+            return trapezoid_type(point_type(x0(i), y0(i), z0(i)),
+                                  point_type(x1(i), y1(i), z1(i)),
+                                  point_type(x2(i), y2(i), z2(i)),
+                                  point_type(x3(i), y3(i), z3(i)));
+        }
+        else
+        {
+            return trapezoid_type(point_type(x0(i), y0(i), 0.0),
+                                  point_type(x1(i), y1(i), 0.0),
+                                  point_type(x2(i), y2(i), 0.0),
+                                  point_type(x3(i), y3(i), 0.0));
+        }
+    }
+    void set(size_t i, trapezoid_type const & t)
+    {
+        x0(i) = t.x0();
+        y0(i) = t.y0();
+        x1(i) = t.x1();
+        y1(i) = t.y1();
+        x2(i) = t.x2();
+        y2(i) = t.y2();
+        x3(i) = t.x3();
+        y3(i) = t.y3();
+        if (ndim() == 3)
+        {
+            z0(i) = t.z0();
+            z1(i) = t.z1();
+            z2(i) = t.z2();
+            z3(i) = t.z3();
+        }
+    }
+    void set(size_t i, point_type const & p0, point_type const & p1, point_type const & p2, point_type const & p3)
+    {
+        x0(i) = p0.x();
+        y0(i) = p0.y();
+        x1(i) = p1.x();
+        y1(i) = p1.y();
+        x2(i) = p2.x();
+        y2(i) = p2.y();
+        x3(i) = p3.x();
+        y3(i) = p3.y();
+        if (ndim() == 3)
+        {
+            z0(i) = p0.z();
+            z1(i) = p1.z();
+            z2(i) = p2.z();
+            z3(i) = p3.z();
+        }
+    }
+    // clang-format off
+    void set(size_t i, value_type x0_value, value_type y0_value,
+             value_type x1_value, value_type y1_value,
+             value_type x2_value, value_type y2_value,
+             value_type x3_value, value_type y3_value)
+    // clang-format on
+    {
+        x0(i) = x0_value;
+        y0(i) = y0_value;
+        x1(i) = x1_value;
+        y1(i) = y1_value;
+        x2(i) = x2_value;
+        y2(i) = y2_value;
+        x3(i) = x3_value;
+        y3(i) = y3_value;
+    }
+    // clang-format off
+    void set(size_t i, value_type x0_value, value_type y0_value, value_type z0_value,
+             value_type x1_value, value_type y1_value, value_type z1_value,
+             value_type x2_value, value_type y2_value, value_type z2_value,
+             value_type x3_value, value_type y3_value, value_type z3_value)
+    // clang-format on
+    {
+        x0(i) = x0_value;
+        y0(i) = y0_value;
+        x1(i) = x1_value;
+        y1(i) = y1_value;
+        x2(i) = x2_value;
+        y2(i) = y2_value;
+        x3(i) = x3_value;
+        y3(i) = y3_value;
+        if (ndim() == 3)
+        {
+            z0(i) = z0_value;
+            z1(i) = z1_value;
+            z2(i) = z2_value;
+            z3(i) = z3_value;
+        }
+    }
+
+    void mirror_x()
+    {
+        size_t const ntrap = size();
+        for (size_t i = 0; i < ntrap; ++i)
+        {
+            y0(i) = -y0(i);
+            y1(i) = -y1(i);
+            y2(i) = -y2(i);
+            y3(i) = -y3(i);
+            if (ndim() == 3)
+            {
+                z0(i) = -z0(i);
+                z1(i) = -z1(i);
+                z2(i) = -z2(i);
+                z3(i) = -z3(i);
+            }
+        }
+    }
+
+    void mirror_y()
+    {
+        size_t const ntrap = size();
+        for (size_t i = 0; i < ntrap; ++i)
+        {
+            x0(i) = -x0(i);
+            x1(i) = -x1(i);
+            x2(i) = -x2(i);
+            x3(i) = -x3(i);
+            if (ndim() == 3)
+            {
+                z0(i) = -z0(i);
+                z1(i) = -z1(i);
+                z2(i) = -z2(i);
+                z3(i) = -z3(i);
+            }
+        }
+    }
+
+    void mirror_z()
+    {
+        if (ndim() != 3)
+        {
+            throw std::out_of_range(
+                std::format("TrapezoidPad::mirror_z: cannot mirror Z axis for ndim {}", int(ndim())));
+        }
+        size_t const ntrap = size();
+        for (size_t i = 0; i < ntrap; ++i)
+        {
+            x0(i) = -x0(i);
+            x1(i) = -x1(i);
+            x2(i) = -x2(i);
+            x3(i) = -x3(i);
+            y0(i) = -y0(i);
+            y1(i) = -y1(i);
+            y2(i) = -y2(i);
+            y3(i) = -y3(i);
+        }
+    }
+
+    void mirror(Axis axis)
+    {
+        switch (axis)
+        {
+        case Axis::X: mirror_x(); break;
+        case Axis::Y: mirror_y(); break;
+        case Axis::Z: mirror_z(); break;
+        default: throw std::invalid_argument("TrapezoidPad::mirror: invalid axis"); break;
+        }
+    }
+
+private:
+
+    void check_constructor_point_size(point_pad_type const & p0, point_pad_type const & p1, point_pad_type const & p2, point_pad_type const & p3)
+    {
+        if (m_p0->size() != m_p1->size() || m_p0->size() != m_p2->size() || m_p0->size() != m_p3->size())
+        {
+            throw std::invalid_argument(
+                std::format("TrapezoidPad::TrapezoidPad: p0.size() {} p1.size() {} p2.size() {} p3.size() {} are not the same",
+                            p0.size(),
+                            p1.size(),
+                            p2.size(),
+                            p3.size()));
+        }
+    }
+
+    std::shared_ptr<point_pad_type> m_p0;
+    std::shared_ptr<point_pad_type> m_p1;
+    std::shared_ptr<point_pad_type> m_p2;
+    std::shared_ptr<point_pad_type> m_p3;
+
+}; /* end class TrapezoidPad */
+
+using TrapezoidPadFp32 = TrapezoidPad<float>;
+using TrapezoidPadFp64 = TrapezoidPad<double>;
+
+/**
+ * Forward declaration of Polygon3d for use in helper classes.
+ */
+template <typename T>
+class Polygon3d;
+
+/**
+ * Forward declaration of PolygonPad for use in Polygon3d handle class.
+ */
+template <typename T>
+class PolygonPad;
+
+/**
+ * Helper class for trapezoidal decomposition of polygons.
+ *
+ * This class implements the sweep line algorithm to decompose polygons into
+ * trapezoids. The decomposition is used for polygon boolean operations.
+ *
+ * @tparam T floating-point type
+ */
+template <typename T>
+class TrapezoidalDecomposer
+{
+
+public:
+
+    using value_type = T;
+    using point_type = Point3d<T>;
+    using trapezoid_pad_type = TrapezoidPad<T>;
+    using ssize_type = ssize_t;
+
+private:
+    /**
+     * Internal structure representing a trapezoid during decomposition for the sweep line algorithm.
+     * The trap has bottom and top lines perpendicular to the Y-axis, defined by lower_y and upper_y.
+     */
+    struct YTrap
+    {
+        value_type lower_y, upper_y;
+        value_type lower_x0, lower_x1, upper_x0, upper_x1;
+        value_type lower_z0, lower_z1, upper_z0, upper_z1; // TODO: currently unused for 2D decomposition
+        size_t source_polygon;
+    }; /* end struct YTrap */
+
+public:
+
+    TrapezoidalDecomposer(uint8_t ndim)
+        : m_trapezoids(trapezoid_pad_type::construct(ndim))
+    {
+    }
+
+    TrapezoidalDecomposer() = delete;
+    TrapezoidalDecomposer(TrapezoidalDecomposer const &) = delete;
+    TrapezoidalDecomposer(TrapezoidalDecomposer &&) = delete;
+    TrapezoidalDecomposer & operator=(TrapezoidalDecomposer const &) = delete;
+    TrapezoidalDecomposer & operator=(TrapezoidalDecomposer &&) = delete;
+    ~TrapezoidalDecomposer() = default;
+
+    /**
+     * Decompose a polygon into trapezoids using vertical sweep line algorithm.
+     *
+     * @param polygon_id ID of the polygon to decompose
+     * @param points Vector of polygon vertices in order
+     * @return Pair of begin and end indices into the trapezoid pad
+     */
+    std::pair<size_t, size_t> decompose(size_t polygon_id, std::vector<point_type> const & points)
+    {
+        if (polygon_id >= m_begins.size())
+        {
+            m_begins.reserve(polygon_id + 1);
+            m_ends.reserve(polygon_id + 1);
+            while (m_begins.size() <= polygon_id)
+            {
+                m_begins.push_back(-1);
+                m_ends.push_back(-1);
+            }
+        }
+
+        if (m_begins[polygon_id] != -1)
+        {
+            return {static_cast<size_t>(m_begins[polygon_id]), static_cast<size_t>(m_ends[polygon_id])};
+        }
+
+        size_t const begin_index = m_trapezoids->size();
+
+        // TODO: Implement sweep line algorithm to generate trapezoids
+        // For now, just mark the range as empty
+
+        size_t const end_index = m_trapezoids->size();
+
+        m_begins[polygon_id] = static_cast<ssize_type>(begin_index);
+        m_ends[polygon_id] = static_cast<ssize_type>(end_index);
+
+        return {begin_index, end_index};
+    }
+
+    size_t num_trapezoids(size_t polygon_id) const
+    {
+        if (polygon_id >= m_begins.size() || m_begins[polygon_id] == -1)
+        {
+            return 0;
+        }
+        return static_cast<size_t>(m_ends[polygon_id] - m_begins[polygon_id]);
+    }
+
+    std::shared_ptr<trapezoid_pad_type> trapezoids() { return m_trapezoids; }
+    std::shared_ptr<trapezoid_pad_type const> trapezoids() const { return m_trapezoids; }
+
+    void clear()
+    {
+        m_begins.clear();
+        m_ends.clear();
+        m_trapezoids = trapezoid_pad_type::construct(m_trapezoids->ndim());
+    }
+
+private:
+
+    std::shared_ptr<trapezoid_pad_type> m_trapezoids;
+    std::vector<ssize_type> m_begins;
+    std::vector<ssize_type> m_ends;
+
+}; /* end class TrapezoidalDecomposer */
+
+using TrapezoidalDecomposerFp32 = TrapezoidalDecomposer<float>;
+using TrapezoidalDecomposerFp64 = TrapezoidalDecomposer<double>;
+
+/**
+ * This class implements the union algorithm for two polygons by decomposing them
+ * into trapezoids and merging overlapping regions.
+ *
+ * @tparam T floating-point type
+ */
+template <typename T>
+class AreaBooleanUnion
+{
+
+public:
+
+    using value_type = T;
+    using point_type = Point3d<T>;
+    using polygon_type = Polygon3d<T>;
+    using polygon_pad_type = PolygonPad<T>;
+
+    AreaBooleanUnion() = default;
+    AreaBooleanUnion(AreaBooleanUnion const &) = delete;
+    AreaBooleanUnion(AreaBooleanUnion &&) = delete;
+    AreaBooleanUnion & operator=(AreaBooleanUnion const &) = delete;
+    AreaBooleanUnion & operator=(AreaBooleanUnion &&) = delete;
+    ~AreaBooleanUnion() = default;
+
+    /**
+     * Compute union of two polygons.
+     *
+     * @param pad Polygon container holding both polygons
+     * @param polygon_id1 ID of first polygon
+     * @param polygon_id2 ID of second polygon
+     * @return polygon pad containing the union of the two polygons
+     */
+    std::shared_ptr<polygon_pad_type> compute(const std::shared_ptr<polygon_pad_type> & pad, size_t polygon_id1, size_t polygon_id2);
+
+}; /* end class AreaBooleanUnion */
+
+using AreaBooleanUnionFp32 = AreaBooleanUnion<float>;
+using AreaBooleanUnionFp64 = AreaBooleanUnion<double>;
+
+/**
+ * This class implements the intersection algorithm for two polygons by decomposing them
+ * into trapezoids and finding overlapping regions.
+ *
+ * @tparam T floating-point type
+ */
+template <typename T>
+class AreaBooleanIntersection
+{
+
+public:
+
+    using value_type = T;
+    using point_type = Point3d<T>;
+    using polygon_type = Polygon3d<T>;
+    using polygon_pad_type = PolygonPad<T>;
+
+    AreaBooleanIntersection() = default;
+    AreaBooleanIntersection(AreaBooleanIntersection const &) = delete;
+    AreaBooleanIntersection(AreaBooleanIntersection &&) = delete;
+    AreaBooleanIntersection & operator=(AreaBooleanIntersection const &) = delete;
+    AreaBooleanIntersection & operator=(AreaBooleanIntersection &&) = delete;
+    ~AreaBooleanIntersection() = default;
+
+    /**
+     * Compute intersection of two polygons.
+     *
+     * @param pad Polygon container holding both polygons
+     * @param polygon_id1 ID of first polygon
+     * @param polygon_id2 ID of second polygon
+     * @return polygon pad containing polygons forming the intersection
+     */
+    std::shared_ptr<polygon_pad_type> compute(const std::shared_ptr<polygon_pad_type> & pad, size_t polygon_id1, size_t polygon_id2);
+
+}; /* end class AreaBooleanIntersection */
+
+using AreaBooleanIntersectionFp32 = AreaBooleanIntersection<float>;
+using AreaBooleanIntersectionFp64 = AreaBooleanIntersection<double>;
+
+/**
+ * This class implements the difference algorithm for two polygons (p1 - p2) by decomposing
+ * them into trapezoids and removing overlapping regions.
+ *
+ * @tparam T floating-point type
+ */
+template <typename T>
+class AreaBooleanDifference
+{
+
+public:
+
+    using value_type = T;
+    using point_type = Point3d<T>;
+    using polygon_type = Polygon3d<T>;
+    using polygon_pad_type = PolygonPad<T>;
+
+    AreaBooleanDifference() = default;
+    AreaBooleanDifference(AreaBooleanDifference const &) = delete;
+    AreaBooleanDifference(AreaBooleanDifference &&) = delete;
+    AreaBooleanDifference & operator=(AreaBooleanDifference const &) = delete;
+    AreaBooleanDifference & operator=(AreaBooleanDifference &&) = delete;
+    ~AreaBooleanDifference() = default;
+
+    /**
+     * Compute difference of two polygons (p1 - p2).
+     *
+     * @param pad Polygon container holding both polygons
+     * @param polygon_id1 ID of first polygon
+     * @param polygon_id2 ID of second polygon to subtract
+     * @return polygon pad containing polygons forming the difference
+     */
+    std::shared_ptr<polygon_pad_type> compute(const std::shared_ptr<polygon_pad_type> & pad, size_t polygon_id1, size_t polygon_id2);
+
+}; /* end class AreaBooleanDifference */
+
+using AreaBooleanDifferenceFp32 = AreaBooleanDifference<float>;
+using AreaBooleanDifferenceFp64 = AreaBooleanDifference<double>;
+
 /**
  * Triangle in three-dimensional space.
  *
@@ -660,7 +1751,11 @@ class TrianglePad
             z2(i) = p2.z();
         }
     }
-    void set(size_t i, value_type x0_value, value_type y0_value, value_type x1_value, value_type y1_value, value_type x2_value, value_type y2_value)
+    // clang-format off
+    void set(size_t i, value_type x0_value, value_type y0_value,
+             value_type x1_value, value_type y1_value,
+             value_type x2_value, value_type y2_value)
+    // clang-format on
     {
         x0(i) = x0_value;
         y0(i) = y0_value;
@@ -669,7 +1764,11 @@ class TrianglePad
         x2(i) = x2_value;
         y2(i) = y2_value;
     }
-    void set(size_t i, value_type x0_value, value_type y0_value, value_type z0_value, value_type x1_value, value_type y1_value, value_type z1_value, value_type x2_value, value_type y2_value, value_type z2_value)
+    // clang-format off
+    void set(size_t i, value_type x0_value, value_type y0_value, value_type z0_value,
+             value_type x1_value, value_type y1_value, value_type z1_value,
+             value_type x2_value, value_type y2_value, value_type z2_value)
+    // clang-format on
     {
         x0(i) = x0_value;
         y0(i) = y0_value;
@@ -789,12 +1888,6 @@ class TrianglePad
 using TrianglePadFp32 = TrianglePad<float>;
 using TrianglePadFp64 = TrianglePad<double>;
 
-/**
- * Forward declaration of PolygonPad for use in Polygon3d handle class.
- */
-template <typename T>
-class PolygonPad;
-
 /**
  * Polygon3d handle class - lightweight view into a polygon stored in PolygonPad.
  *
@@ -934,6 +2027,9 @@ class Polygon3d
  * with each polygon defined by a range [start, end) in the node list. This avoids
  * memory duplication compared to storing line segments.
  *
+ * Note that PolygonPad assumes polygons are in a XY plane currently.
+ * TODO: Extend to 3D polygons in arbitrary planes.
+ *
  * All nodes follow the right-hand rule:
  *  - Counter-clockwise ordering = positive area polygon
  *  - Clockwise ordering = negative area polygon (holes)
@@ -965,9 +2061,11 @@ class PolygonPad
     using segment_type = Segment3d<T>;
     using point_pad_type = PointPad<T>;
     using polygon_type = Polygon3d<T>;
+    using polygon_pad_type = PolygonPad<T>;
     using segment_pad_type = SegmentPad<T>;
     using curve_pad_type = CurvePad<T>;
     using rtree_type = RTree<segment_type, BoundBox3d<T>, RTreeValueOps<segment_type, BoundBox3d<T>>>;
+    using trapezoid_decomposer_type = TrapezoidalDecomposer<T>;
 
     template <typename... Args>
     static std::shared_ptr<PolygonPad<T>> construct(Args &&... args)
@@ -978,6 +2076,7 @@ class PolygonPad
     PolygonPad(uint8_t ndim, ctor_passkey const &)
         : m_points(point_pad_type::construct(ndim))
         , m_rtree(std::make_unique<rtree_type>())
+        , m_decomposer(ndim)
     {
     }
 
@@ -1129,6 +2228,51 @@ class PolygonPad
      */
     void rebuild_rtree();
 
+    /**
+     * Decompose a polygon into trapezoids using vertical sweep line algorithm.
+     *
+     * @param polygon_id ID of the polygon to decompose
+     * @return Pair of begin and end indices into the decomposer's trapezoid pad
+     * @throws std::out_of_range if polygon_id is invalid
+     */
+    std::pair<size_t, size_t> decompose_to_trapezoid(size_t polygon_id);
+
+    /**
+     * Compute union of two polygons using trapezoidal decomposition.
+     *
+     * @param p1 First polygon
+     * @param p2 Second polygon
+     * @return polygon pad forming the union
+     */
+    std::shared_ptr<polygon_pad_type> boolean_union(polygon_type const & p1, polygon_type const & p2)
+    {
+        return m_boolean_union.compute(this->shared_from_this(), p1.polygon_id(), p2.polygon_id());
+    }
+
+    /**
+     * Compute intersection of two polygons using trapezoidal decomposition.
+     *
+     * @param p1 First polygon
+     * @param p2 Second polygon
+     * @return polygon pad forming the intersection
+     */
+    std::shared_ptr<polygon_pad_type> boolean_intersection(polygon_type const & p1, polygon_type const & p2)
+    {
+        return m_boolean_intersection.compute(this->shared_from_this(), p1.polygon_id(), p2.polygon_id());
+    }
+
+    /**
+     * Compute difference of two polygons (p1 - p2).
+     *
+     * @param p1 First polygon
+     * @param p2 Second polygon to subtract
+     * @return polygon pad forming the difference
+     */
+    std::shared_ptr<polygon_pad_type> boolean_difference(polygon_type const & p1, polygon_type const & p2)
+    {
+        return m_boolean_difference.compute(this->shared_from_this(), p1.polygon_id(), p2.polygon_id());
+    }
+
 private:
 
     friend class Polygon3d<T>;
@@ -1140,13 +2284,10 @@ class PolygonPad
     SimpleCollector<ssize_type> m_ends;
     std::unique_ptr<rtree_type> m_rtree;
 
-    // TODO: Add trapezoidal map data structure for polygon boolean operations
-    // Reference: http://www0.cs.ucl.ac.uk/staff/m.slater/Teaching/CG/1997-98/Solutions/Trap/
-    // Will need:
-    //  - Trapezoid decomposition structure
-    //  - Point location query structure
-    //  - Integration with node lists for synchronization
-
+    trapezoid_decomposer_type m_decomposer;
+    AreaBooleanUnion<T> m_boolean_union;
+    AreaBooleanIntersection<T> m_boolean_intersection;
+    AreaBooleanDifference<T> m_boolean_difference;
 }; /* end class PolygonPad */
 
 using PolygonPadFp32 = PolygonPad<float>;
@@ -1201,6 +2342,12 @@ Polygon3d<T> PolygonPad<T>::add_polygon(std::vector<point_type> const & nodes)
 
     for (point_type const & node : nodes)
     {
+        // check if the point is in XY plane
+        if (node.z() != 0)
+        {
+            throw std::invalid_argument("PolygonPad::add_polygon: all nodes must lie in the XY plane (z=0)");
+        }
+
         m_points->append(node);
     }
 
@@ -1209,7 +2356,8 @@ Polygon3d<T> PolygonPad<T>::add_polygon(std::vector<point_type> const & nodes)
     m_begins.push_back(begin_index);
     m_ends.push_back(end_index);
 
-    polygon_type polygon(this->shared_from_this(), polygon_id, typename polygon_type::ctor_passkey());
+    std::shared_ptr<PolygonPad<T> const> const_this = this->shared_from_this();
+    polygon_type polygon(const_this, polygon_id, typename polygon_type::ctor_passkey());
     rebuild_polygon_rtree(polygon);
 
     return polygon;
@@ -1469,6 +2617,52 @@ void PolygonPad<T>::rebuild_polygon_rtree(polygon_type const & polygon)
     }
 }
 
+template <typename T>
+std::pair<size_t, size_t> PolygonPad<T>::decompose_to_trapezoid(size_t polygon_id)
+{
+    if (polygon_id >= num_polygons())
+    {
+        throw std::out_of_range(
+            std::format("PolygonPad::decompose_to_trapezoid: polygon_id {} >= num_polygons {}",
+                        polygon_id,
+                        num_polygons()));
+    }
+
+    polygon_type polygon = get_polygon(polygon_id);
+    std::vector<point_type> points;
+    points.reserve(polygon.nnode());
+    for (size_t i = 0; i < polygon.nnode(); ++i)
+    {
+        points.push_back(polygon.node(i));
+    }
+
+    return m_decomposer.decompose(polygon_id, points);
+}
+
+template <typename T>
+std::shared_ptr<PolygonPad<T>> AreaBooleanUnion<T>::compute(const std::shared_ptr<PolygonPad<T>> & pad, size_t polygon_id1, size_t polygon_id2)
+{
+    // TODO: A proper implementation would merge overlapping regions using trapezoidal decomposition
+    auto empty_pad = PolygonPad<T>::construct(pad->ndim());
+    return empty_pad;
+}
+
+template <typename T>
+std::shared_ptr<PolygonPad<T>> AreaBooleanIntersection<T>::compute(const std::shared_ptr<PolygonPad<T>> & pad, size_t polygon_id1, size_t polygon_id2)
+{
+    // TODO:  A proper implementation would find overlapping regions using trapezoidal decomposition
+    auto empty_pad = PolygonPad<T>::construct(pad->ndim());
+    return empty_pad;
+}
+
+template <typename T>
+std::shared_ptr<PolygonPad<T>> AreaBooleanDifference<T>::compute(const std::shared_ptr<PolygonPad<T>> & pad, size_t polygon_id1, size_t polygon_id2)
+{
+    // TODO:  A proper implementation would subtract overlapping regions using trapezoidal decomposition
+    auto empty_pad = PolygonPad<T>::construct(pad->ndim());
+    return empty_pad;
+}
+
 } /* end namespace modmesh */
 
 // vim: set ff=unix fenc=utf8 et sw=4 ts=4 sts=4:
diff --git a/cpp/modmesh/universe/pymod/wrap_polygon.cpp b/cpp/modmesh/universe/pymod/wrap_polygon.cpp
index cb7e9f96..32d5ae9a 100644
--- a/cpp/modmesh/universe/pymod/wrap_polygon.cpp
+++ b/cpp/modmesh/universe/pymod/wrap_polygon.cpp
@@ -143,7 +143,29 @@ WrapPolygonPad<T>::WrapPolygonPad(pybind11::module & mod, const char * pyname, c
                 return results;
             },
             py::arg("box"))
-        .def("rebuild_rtree", &wrapped_type::rebuild_rtree);
+        .def("rebuild_rtree", &wrapped_type::rebuild_rtree)
+        .def(
+            "decompose_to_trapezoid",
+            [](wrapped_type & self, size_t polygon_id)
+            {
+                return self.decompose_to_trapezoid(polygon_id);
+            },
+            py::arg("polygon_id"))
+        .def(
+            "boolean_union",
+            &wrapped_type::boolean_union,
+            py::arg("p1"),
+            py::arg("p2"))
+        .def(
+            "boolean_intersection",
+            &wrapped_type::boolean_intersection,
+            py::arg("p1"),
+            py::arg("p2"))
+        .def(
+            "boolean_difference",
+            &wrapped_type::boolean_difference,
+            py::arg("p1"),
+            py::arg("p2"));
 }
 
 void wrap_polygon(pybind11::module & mod)
diff --git a/tests/test_polygon3d.py b/tests/test_polygon3d.py
index d69d57eb..a937c16a 100644
--- a/tests/test_polygon3d.py
+++ b/tests/test_polygon3d.py
@@ -332,6 +332,214 @@ def test_polygon_identity_same_pad_different_id(self):
         self.assertFalse(polygon1.is_same(polygon2))
         self.assertTrue(polygon1.is_not_same(polygon2))
 
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_union_simple(self):
+        """Test polygon boolean union with two overlapping squares."""
+        pad = self.PolygonPad(ndim=2)
+
+        # First square: (0,0) to (2,2)
+        square1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(2.0, 0.0, 0.0),
+            self.Point(2.0, 2.0, 0.0),
+            self.Point(0.0, 2.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(square1_nodes)
+
+        # Second square: (1,1) to (3,3) - overlaps with first
+        square2_nodes = [
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(3.0, 1.0, 0.0),
+            self.Point(3.0, 3.0, 0.0),
+            self.Point(1.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(square2_nodes)
+
+        # Call boolean_union - should not crash even if implementation is TODO
+        try:
+            result = pad.boolean_union(polygon1, polygon2)
+            # Implementation is not complete yet, but should return a list
+            self.assertIsInstance(result, self.PolygonPad)
+        except NotImplementedError:
+            # Expected if method raises NotImplementedError
+            pass
+
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_intersection_simple(self):
+        """Test polygon boolean intersection with two overlapping squares."""
+        pad = self.PolygonPad(ndim=2)
+
+        # First square: (0,0) to (2,2)
+        square1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(2.0, 0.0, 0.0),
+            self.Point(2.0, 2.0, 0.0),
+            self.Point(0.0, 2.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(square1_nodes)
+
+        # Second square: (1,1) to (3,3) - overlaps with first
+        # Intersection should be (1,1) to (2,2)
+        square2_nodes = [
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(3.0, 1.0, 0.0),
+            self.Point(3.0, 3.0, 0.0),
+            self.Point(1.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(square2_nodes)
+
+        # Call boolean_intersection - should not crash even
+        # if implementation is TODO
+        try:
+            result = pad.boolean_intersection(polygon1, polygon2)
+            # Implementation is not complete yet, but should return a list
+            self.assertIsInstance(result, self.PolygonPad)
+        except NotImplementedError:
+            # Expected if method raises NotImplementedError
+            pass
+
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_difference_simple(self):
+        """Test polygon boolean difference with two overlapping squares."""
+        pad = self.PolygonPad(ndim=2)
+
+        # First square: (0,0) to (2,2)
+        square1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(2.0, 0.0, 0.0),
+            self.Point(2.0, 2.0, 0.0),
+            self.Point(0.0, 2.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(square1_nodes)
+
+        # Second square: (1,1) to (3,3) - overlaps with first
+        # Difference (polygon1 - polygon2) should be L-shaped region
+        square2_nodes = [
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(3.0, 1.0, 0.0),
+            self.Point(3.0, 3.0, 0.0),
+            self.Point(1.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(square2_nodes)
+
+        # Call boolean_difference - should not crash even
+        # if implementation is TODO
+        try:
+            result = pad.boolean_difference(polygon1, polygon2)
+            # Implementation is not complete yet, but should return a list
+            self.assertIsInstance(result, self.PolygonPad)
+        except NotImplementedError:
+            # Expected if method raises NotImplementedError
+            pass
+
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_union_non_overlapping(self):
+        """Test polygon boolean union with two non-overlapping squares."""
+        pad = self.PolygonPad(ndim=2)
+
+        # First square: (0,0) to (1,1)
+        square1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(1.0, 0.0, 0.0),
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(0.0, 1.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(square1_nodes)
+
+        # Second square: (2,2) to (3,3) - no overlap
+        square2_nodes = [
+            self.Point(2.0, 2.0, 0.0),
+            self.Point(3.0, 2.0, 0.0),
+            self.Point(3.0, 3.0, 0.0),
+            self.Point(2.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(square2_nodes)
+
+        # Call boolean_union - should not crash
+        try:
+            result = pad.boolean_union(polygon1, polygon2)
+            self.assertIsInstance(result, self.PolygonPad)
+        except NotImplementedError:
+            pass
+
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_intersection_non_overlapping(self):
+        """Test polygon boolean intersection with two non-overlapping squares.
+        """
+        pad = self.PolygonPad(ndim=2)
+
+        # First square: (0,0) to (1,1)
+        square1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(1.0, 0.0, 0.0),
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(0.0, 1.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(square1_nodes)
+
+        # Second square: (2,2) to (3,3) - no overlap
+        # Intersection should be empty
+        square2_nodes = [
+            self.Point(2.0, 2.0, 0.0),
+            self.Point(3.0, 2.0, 0.0),
+            self.Point(3.0, 3.0, 0.0),
+            self.Point(2.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(square2_nodes)
+
+        # Call boolean_intersection - should not crash
+        try:
+            result = pad.boolean_intersection(polygon1, polygon2)
+            self.assertIsInstance(result, self.PolygonPad)
+        except NotImplementedError:
+            pass
+
+    # TODO: Verify and implement boolean operations
+    # once the C++ side is complete.
+    def test_boolean_operations_triangle(self):
+        """Test polygon boolean operations with triangular polygons."""
+        pad = self.PolygonPad(ndim=2)
+
+        # Triangle 1: Right triangle at origin
+        triangle1_nodes = [
+            self.Point(0.0, 0.0, 0.0),
+            self.Point(2.0, 0.0, 0.0),
+            self.Point(0.0, 2.0, 0.0)
+        ]
+        polygon1 = pad.add_polygon(triangle1_nodes)
+
+        # Triangle 2: Right triangle shifted
+        triangle2_nodes = [
+            self.Point(1.0, 1.0, 0.0),
+            self.Point(3.0, 1.0, 0.0),
+            self.Point(1.0, 3.0, 0.0)
+        ]
+        polygon2 = pad.add_polygon(triangle2_nodes)
+
+        # Test all three operations - should not crash
+        try:
+            result_union = pad.boolean_union(polygon1, polygon2)
+            self.assertIsInstance(result_union, self.PolygonPad)
+        except NotImplementedError:
+            pass
+
+        try:
+            result_intersection = pad.boolean_intersection(polygon1, polygon2)
+            self.assertIsInstance(result_intersection, self.PolygonPad)
+        except NotImplementedError:
+            pass
+
+        try:
+            result_difference = pad.boolean_difference(polygon1, polygon2)
+            self.assertIsInstance(result_difference, self.PolygonPad)
+        except NotImplementedError:
+            pass
+
 
 class Polygon3dFp32TC(Polygon3dTB, unittest.TestCase):