diff --git a/src/classy_blocks/optimize/cell.py b/src/classy_blocks/optimize/cell.py
index 45e6330..bba3565 100644
--- a/src/classy_blocks/optimize/cell.py
+++ b/src/classy_blocks/optimize/cell.py
@@ -14,7 +14,7 @@ class NoCommonSidesError(Exception):
 
 
 class CellBase(abc.ABC):
-    side_names: ClassVar[List[str]]
+    side_names: ClassVar[List[OrientType]]
     side_indexes: ClassVar[List[IndexType]]
     edge_pairs: ClassVar[List[Tuple[int, int]]]
 
@@ -75,6 +75,19 @@ def add_neighbour(self, candidate: "CellBase") -> bool:
         except NoCommonSidesError:
             return False
 
+    @property
+    def boundary(self) -> Set[int]:
+        """Returns a list of indexes that define sides on boundary"""
+        boundary = set()
+
+        for i, side_name in enumerate(self.side_names):
+            side_indexes = self.side_indexes[i]
+
+            if self.neighbours[side_name] is None:
+                boundary.update({self.indexes[si] for si in side_indexes})
+
+        return boundary
+
     @property
     def points(self) -> NPPointListType:
         """A list of points defining this cell, as a numpy array"""
diff --git a/src/classy_blocks/optimize/junction.py b/src/classy_blocks/optimize/junction.py
index 43d569d..b20bcec 100644
--- a/src/classy_blocks/optimize/junction.py
+++ b/src/classy_blocks/optimize/junction.py
@@ -19,7 +19,7 @@ class IndexedLink:
 
 
 class Junction:
-    """A class that collects Cells/Blocks that
+    """A class that collects Cells that
     share the same Vertex"""
 
     def __init__(self, points: NPPointListType, index: int):
@@ -28,7 +28,7 @@ def __init__(self, points: NPPointListType, index: int):
 
         self.cells: Set[CellBase] = set()
 
-        self.connections: List[Junction] = []
+        self.neighbours: List[Junction] = []
 
         self.clamp: Optional[ClampBase] = None
         self.links: List[IndexedLink] = []
@@ -58,8 +58,8 @@ def add_connection(self, to: "Junction") -> bool:
         for cell in self.cells:
             for connection in cell.connections:
                 if connection.indexes == junction_indexes:
-                    if to not in self.connections:
-                        self.connections.append(to)
+                    if to not in self.neighbours:
+                        self.neighbours.append(to)
                         return True
 
         return False
@@ -73,6 +73,15 @@ def add_clamp(self, clamp: ClampBase) -> None:
     def add_link(self, link: LinkBase, follower_index: int) -> None:
         self.links.append(IndexedLink(link, follower_index))
 
+    @property
+    def is_boundary(self) -> bool:
+        """Returns True if this junction lies on boundary"""
+        for cell in self.cells:
+            if self.index in cell.boundary:
+                return True
+
+        return False
+
     @property
     def quality(self) -> float:
         """Returns average quality of all cells at this junction;
diff --git a/tests/test_modify/test_grid.py b/tests/test_modify/test_grid.py
index 4ff8624..66b73a2 100644
--- a/tests/test_modify/test_grid.py
+++ b/tests/test_modify/test_grid.py
@@ -1,17 +1,24 @@
 import numpy as np
 from parameterized import parameterized
 
+from classy_blocks.construct.flat.sketches.grid import Grid as GridSketch
+from classy_blocks.construct.stack import ExtrudedStack
+from classy_blocks.mesh import Mesh
 from classy_blocks.optimize.grid import Grid
+from classy_blocks.util import functions as f
 from tests.fixtures.mesh import MeshTestCase
 
 
 class GridTests(MeshTestCase):
+    def get_grid(self, mesh: Mesh) -> Grid:
+        points = np.array([vertex.position for vertex in mesh.vertices])
+        addresses = [block.indexes for block in mesh.blocks]
+
+        return Grid(points, addresses)
+
     def setUp(self):
         super().setUp()
-
-        points = np.array([vertex.position for vertex in self.mesh.vertices])
-        addresses = [block.indexes for block in self.mesh.blocks]
-        self.grid = Grid(points, addresses)
+        self.grid = self.get_grid(self.mesh)
 
     def test_cells_quantity(self):
         self.assertEqual(len(self.grid.cells), len(self.mesh.blocks))
@@ -19,6 +26,32 @@ def test_cells_quantity(self):
     def test_junctions_quantity(self):
         self.assertEqual(len(self.grid.junctions), len(self.mesh.vertices))
 
+    def test_junction_boundary(self):
+        # In this case, ALL junctions are on boundary
+        for junction in self.grid.junctions:
+            self.assertTrue(junction.is_boundary)
+
+    def test_junction_internal(self):
+        sketch = GridSketch([0, 0, 0], [1, 1, 0], 2, 2)
+        stack = ExtrudedStack(sketch, 1, 2)
+
+        mesh = Mesh()
+        mesh.add(stack)
+
+        try:
+            mesh.assemble()  # will fail because there are no chops
+        except:
+            pass
+
+        grid = self.get_grid(mesh)
+
+        for junction in grid.junctions:
+            if f.norm(junction.point - f.vector(0.5, 0.5, 0.5)) < 0.01:
+                self.assertFalse(junction.is_boundary)
+                continue
+
+            self.assertTrue(junction.is_boundary)
+
     @parameterized.expand(((0, 1), (1, 2), (2, 3), (3, 1), (4, 1), (5, 2), (6, 3), (7, 1)))
     def test_junction_cells(self, index, count):
         """Each junction contains cells that include that vertex"""