diff --git a/src/PhysicalMeshes.jl b/src/PhysicalMeshes.jl
index a2ea9ac..5f136b6 100644
--- a/src/PhysicalMeshes.jl
+++ b/src/PhysicalMeshes.jl
@@ -78,6 +78,15 @@ export
     AbstractLine2D, AbstractLine3D,
     Line, Line2D,
 
+    # Plane
+    Plane,
+    coplanar,
+
+    # Polygon
+    Polygon2D, Polygon3D,
+    polygon_rect, polygon_regular,
+    isconvex,
+
     # Ray
     Ray2D, Ray3D,
     intersect,
@@ -123,6 +132,9 @@ abstract type AbstractLine{T} <: AbstractGeometryType{T} end
 abstract type AbstractLine2D{T} <: AbstractLine{T} end
 abstract type AbstractLine3D{T} <: AbstractLine{T} end
 
+abstract type AbstractPlane{T} <: AbstractGeometryType{T} end
+abstract type AbstractPlane3D{T} <: AbstractPlane{T} end
+
 abstract type AbstractPolygon{T} <: AbstractGeometryType{T} end
 abstract type AbstractPolygon2D{T} <: AbstractPolygon{T} end
 abstract type AbstractPolygon3D{T} <: AbstractPolygon{T} end
@@ -156,6 +168,7 @@ include("core/predicates.jl")
 include("Sphere.jl")
 
 include("Line.jl")
+include("Plane.jl")
 include("Polygon.jl")
 include("Ray.jl")
 include("Cube.jl")
diff --git a/src/Plane.jl b/src/Plane.jl
new file mode 100644
index 0000000..57848c0
--- /dev/null
+++ b/src/Plane.jl
@@ -0,0 +1,16 @@
+"""
+$(TYPEDEF)
+$(TYPEDFIELDS)
+"""
+struct Plane{T} <: AbstractPlane{T}
+    a::PVector{T}
+    b::PVector{T}
+    c::PVector{T}
+end
+
+normal(a,b,c) = normalize(cross(b-a, ustrip(c-a)))
+normal(plane::Plane) = normal(plane.a, plane.b, plane.c)
+
+distance(p::AbstractPoint3D, plane::Plane) = abs(dot(normal(plane), p-plane.a))
+
+coplanar(p::AbstractPoint3D, plane::Plane, threshold::Number) = distance(p, plane) <= threshold
\ No newline at end of file
diff --git a/src/Polygon.jl b/src/Polygon.jl
index bc57a04..6bb0a66 100644
--- a/src/Polygon.jl
+++ b/src/Polygon.jl
@@ -3,7 +3,15 @@ $(TYPEDEF)
 $(TYPEDFIELDS)
 """
 struct Polygon2D{T} <: AbstractPolygon2D{T}
-    points
+    vertices
+end
+
+function Polygon2D(vertices)
+    if length(vertices) < 3
+        error("Number of vertices less than 3, a concrete polygon should have more than 2 vertices")
+    end
+
+    Polygon2D(vertices)
 end
 
 """
@@ -11,5 +19,88 @@ $(TYPEDEF)
 $(TYPEDFIELDS)
 """
 struct Polygon3D{T} <: AbstractPolygon3D{T}
-    points
+    vertices
+end
+
+function Polygon3D(vertices)
+    if length(vertices) < 3
+        error("Number of vertices less than 3, a concrete polygon should have more than 2 vertices")
+    end
+
+    if !coplanar(vertices, one(vertices[1].x) * 1e-6)
+        @warn "The polygon is not coplanar with 1e-6 threshold!"
+    end
+
+    Polygon3D(vertices)
+end
+
+function isconvex(polygon::Polygon2D)
+    
+end
+
+"""
+$(TYPEDSIGNATURES)
+The internal angles of a convex polygon are all smaller than π,
+or equivalently, the cross products of any two adjacent edges have the same direction.
+"""
+function isconvex(polygon::Polygon3D)
+    n = length(polygon.vertices)
+    if n < 4
+        return true
+    end
+
+    N = normal(polygon)
+    for i in 1:n
+        v1 = polygon.vertices[i]
+        v2 = polygon.vertices[mod1(i+1, n)]
+        v3 = polygon.vertices[mod1(i+2, n)]
+        if dot(cross(v2 - v1, v3 - v2), N) < 0
+            return false
+        end
+    end
+    return true
+end
+
+"""
+$(TYPEDSIGNATURES)
+Compute the unit normal vector from the first three vertices
+"""
+normal(polygon::Polygon3D) = normal(polygon.vertices[1:3]...)
+
+"""
+$(TYPEDSIGNATURES)
+
+First construct a plane from the first three vertices, then iteratively check whether all other vertices are coplanar with this plane or not.
+"""
+function coplanar(vertices, threshold)
+    n = length(vertices)
+    
+    # A concrete polygon have more than 2 vertices
+    if n == 3
+        return true
+    else
+        plane = Plane(vertices[1:3]...)
+        for i = 4:n
+            if !coplanar(vertices[i], plane, threshold)
+                return false
+            end
+        end
+        return true
+    end
+end
+
+coplanar(polygon::Polygon3D, threshold) = coplanar(polygon.vertices, threshold)
+
+function is_inbound(pos, polygon)
+
+end
+
+### Some common polygons
+
+function polygon_rect()
+    
+end
+
+function polygon_regular()
+    
 end
\ No newline at end of file
diff --git a/src/Ray.jl b/src/Ray.jl
index b692c2c..854c6d4 100644
--- a/src/Ray.jl
+++ b/src/Ray.jl
@@ -29,6 +29,8 @@ struct Ray3D{T<:Number} <: AbstractRay3D{T}
     n::PVector{T}
 end
 
+### Ray2D reflect from Line2D
+
 function intersect(ray::Ray2D, line::Line2D)
     AB = line.b - line.a
     ray_norm = PVector2D(-ray.n.y, ray.n.x) # ⟂ ray
@@ -39,15 +41,15 @@ function intersect(ray::Ray2D, line::Line2D)
     AX = ray.x - line.a
     t1 = dot(PVector2D(-AB.y, AB.x), AX) / dot_ray_norm_AB
     t2 = dot(AX, ray_norm) / dot_ray_norm_AB
-    if t1 >= 0 && t2 >= 0 && t2 <= 1
-        intersection_point = ray.x + ray.n * t1
+    if t1 >= 0 && t2 >= 0 && t2 <= one(t2)
+        intersection_point = ray.x + ray.n * t1 # rescaling
         return true, intersection_point
     else
         return false, PVector2D()  # no intersection
     end
 end
 
-reflect(vec::PVector2D, n::PVector2D) = vec - 2*ustrip(n)*dot(vec, ustrip(n))
+reflect(vec::AbstractPoint, n::AbstractPoint) = vec - 2*ustrip(n)*dot(vec, ustrip(n))
 
 function reflect(ray::Ray2D, line::Line2D)
     hit, intersection = intersect(ray, line)
@@ -56,4 +58,29 @@ function reflect(ray::Ray2D, line::Line2D)
     else
         return nothing
     end
+end
+
+### Ray3D reflect from Plane
+
+function intersect(ray::Ray3D, plane::Plane)
+    N = normal(plane)
+    denom = dot(ray.n, N)
+    if norm(denom) < 1e-6 * unit(denom)
+        return false, PVector(0, 0, 0)  # ∥, no intersection
+    end
+    t = dot(plane.a - ray.x, N) / denom
+    if t < 0
+        return false, PVector(0, 0, 0)  # no intersection
+    end
+    intersection_point = ray.x + ray.n * t # rescaling
+    return true, intersection_point
+end
+
+function reflect(ray::Ray3D, plane::Plane)
+    hit, intersection = intersect(ray, plane)
+    if hit
+        return Ray3D(intersection, reflect(ray.n, normal(plane)))
+    else
+        return nothing
+    end
 end
\ No newline at end of file
diff --git a/test/testRay.jl b/test/testRay.jl
index 8264e51..ba31a92 100644
--- a/test/testRay.jl
+++ b/test/testRay.jl
@@ -1,4 +1,4 @@
-@testset "Ray Unitless" begin
+@testset "Ray2D Unitless" begin
     ray = Ray2D(PVector2D(), PVector2D(1.0, 1.0))
     @test ray.theta == π/4
 
@@ -16,7 +16,7 @@
     @test isnothing(reflect(ray, Line2D(PVector(2, 1), PVector(2, 0))))
 end
 
-@testset "Ray Unitful" begin
+@testset "Ray2D Unitful" begin
     ray = Ray2D(PVector2D(u"m"), PVector2D(1.0u"m", 1.0u"m"))
     @test ray.theta == π/4
 
@@ -32,4 +32,42 @@ end
     @test ray_reflect.theta == 3/4 * π
 
     @test isnothing(reflect(ray, Line2D(PVector(2u"m", 1u"m"), PVector(2u"m", 0u"m"))))
+end
+
+@testset "Ray3D Unitless" begin
+    ray = Ray3D(PVector(), PVector(0.0, 1.0, 1.0))
+    plane = Plane(PVector(1, 1, 1), PVector(1, 1, -1), PVector(-1, 1, 1))
+    @test normal(plane) == PVector(0.0, 1.0, 0.0)
+    hit, intersection_point = intersect(ray, plane)
+    @test intersection_point == PVector(0.0, 1.0, 1.0)
+    ray_reflect = reflect(ray, plane)
+    @test ray_reflect.n == PVector(0.0, -1.0, 1.0)
+    
+    #: zero incident angle
+    ray = Ray3D(PVector(), PVector(0.0, 1.0, 0.0))
+    plane = Plane(PVector(1, 1, 1), PVector(1, 1, -1), PVector(-1, 1, 1))
+    @test normal(plane) == PVector(0.0, 1.0, 0.0)
+    hit, intersection_point = intersect(ray, plane)
+    @test intersection_point == PVector(0.0, 1.0, 0.0)
+    ray_reflect = reflect(ray, plane)
+    @test ray_reflect.n == PVector(0.0, -1.0, 0.0)
+end
+
+@testset "Ray3D Unitful" begin
+    ray = Ray3D(PVector(u"m"), PVector(0.0, 1.0, 1.0, u"m"))
+    plane = Plane(PVector(1, 1, 1, u"m"), PVector(1, 1, -1, u"m"), PVector(-1, 1, 1, u"m"))
+    @test normal(plane) == PVector(0.0, 1.0, 0.0, u"m")
+    hit, intersection_point = intersect(ray, plane)
+    @test intersection_point == PVector(0.0, 1.0, 1.0, u"m")
+    ray_reflect = reflect(ray, plane)
+    @test ray_reflect.n == PVector(0.0, -1.0, 1.0, u"m")
+    
+    #: zero incident angle
+    ray = Ray3D(PVector(u"m"), PVector(0.0, 1.0, 0.0, u"m"))
+    plane = Plane(PVector(1, 1, 1, u"m"), PVector(1, 1, -1, u"m"), PVector(-1, 1, 1, u"m"))
+    @test normal(plane) == PVector(0.0, 1.0, 0.0, u"m")
+    hit, intersection_point = intersect(ray, plane)
+    @test intersection_point == PVector(0.0, 1.0, 0.0, u"m")
+    ray_reflect = reflect(ray, plane)
+    @test ray_reflect.n == PVector(0.0, -1.0, 0.0, u"m")
 end
\ No newline at end of file