diff --git a/nhf/checks.py b/nhf/checks.py
index 468cc8e..19914ab 100644
--- a/nhf/checks.py
+++ b/nhf/checks.py
@@ -1,6 +1,30 @@
 import cadquery as Cq
 
 def binary_intersection(a: Cq.Assembly) -> Cq.Shape:
-    objs = [s.toCompound() for _, s in a.traverse() if isinstance(s, Cq.Assembly)]
+    objs = [s.toCompound() for _, s in a.traverse()
+            if isinstance(s, Cq.Assembly)]
     obj1, obj2 = objs[:2]
     return obj1.intersect(obj2)
+
+
+def pairwise_intersection(assembly: Cq.Assembly, tol: float=1e-6) -> list[(str, str, float)]:
+    """
+    Given an assembly, test the pairwise intersection volume of its components.
+    Return the pairs whose intersection volume exceeds `tol`.
+    """
+    m = {name: (i, shape.moved(loc))
+         for i, (shape, name, loc, _)
+         in enumerate(assembly)}
+    result = []
+    for name, (i1, sh1) in m.items():
+        for name2, (i2, sh2) in m.items():
+            if name == name2:
+                assert i1 == i2
+                continue
+            if i2 <= i1:
+                # Remove the upper diagonal
+                continue
+            vol = sh1.intersect(sh2).Volume()
+            if vol > tol:
+                result.append((name, name2, vol))
+    return result
diff --git a/nhf/materials.py b/nhf/materials.py
index c0adbde..0cb266d 100644
--- a/nhf/materials.py
+++ b/nhf/materials.py
@@ -31,6 +31,7 @@ class Material(Enum):
 
     WOOD_BIRCH = 0.8, _color('bisque', 0.9)
     PLASTIC_PLA = 0.5, _color('azure3', 0.6)
+    RESIN_TRANSPARENT = 1.1, _color('cadetblue2', 0.6)
     ACRYLIC_BLACK = 0.5, _color('gray50', 0.6)
 
     def __init__(self, density: float, color: Cq.Color):
diff --git a/nhf/test.py b/nhf/test.py
index c6e94a6..fc9d15b 100644
--- a/nhf/test.py
+++ b/nhf/test.py
@@ -1,9 +1,27 @@
+"""
+Unit tests for tooling
+"""
 import unittest
 import cadquery as Cq
 from nhf.build import Model, target
+import nhf.checks
+import nhf.utils
+
+# Color presets for testing purposes
+color_parent = Cq.Color(0.7, 0.7, 0.5, 0.5)
+color_child = Cq.Color(0.5, 0.7, 0.7, 0.5)
+
+def makeSphere(r: float) -> Cq.Solid:
+    """
+    Makes a full sphere. The default function makes a hemisphere
+    """
+    return Cq.Solid.makeSphere(r, angleDegrees1=-90)
 
 class BuildScaffold(Model):
 
+    def __init__(self):
+        super().__init__(name="scaffold")
+
     @target(name="obj1")
     def o1(self):
         return Cq.Solid.makeBox(10, 10, 10)
@@ -19,6 +37,115 @@ class TestBuild(unittest.TestCase):
         self.assertEqual(s.target_names, names)
         self.assertEqual(s.check_all(), len(names))
 
+class TestChecks(unittest.TestCase):
+
+    def intersect_test_case(self, offset):
+        assembly = (
+            Cq.Assembly()
+            .add(Cq.Solid.makeBox(10, 10, 10),
+                 name="c1",
+                 loc=Cq.Location((0, 0, 0)))
+            .add(Cq.Solid.makeBox(10, 10, 10),
+                 name="c2",
+                 loc=Cq.Location((0, 0, offset)))
+        )
+        coll = nhf.checks.pairwise_intersection(assembly)
+        if -10 < offset and offset < 10:
+            self.assertEqual(len(coll), 1)
+        else:
+            self.assertEqual(coll, [])
+
+    def test_intersect(self):
+        for offset in [9, 10, 11, -10]:
+            with self.subTest(offset=offset):
+                self.intersect_test_case(offset)
+
+class TestUtils(unittest.TestCase):
+
+    def test_tag_point(self):
+        """
+        A board with 3 holes of unequal sizes. Each hole is marked
+        """
+        p4x, p4y = 5, 5
+        p3x, p3y = 0, 0
+        p2x, p2y = -5, 0
+        board = (
+            Cq.Workplane('XY')
+            .box(15, 15, 5)
+            .faces("<Z")
+            .workplane()
+            .pushPoints([(p4x, p4y)])
+            .hole(4, depth=2)
+            .pushPoints([(p3x, p3y)])
+            .hole(3, depth=1.5)
+            .pushPoints([(p2x, p2y)])
+            .hole(2, depth=1)
+        )
+        board.moveTo(p4x, p4y).tagPoint("h4")
+        board.moveTo(p3x, p3y).tagPoint("h3")
+        board.moveTo(p2x, p2y).tagPoint("h2")
+
+        assembly = (
+            Cq.Assembly()
+            .add(board, name="board", color=color_parent)
+            .add(makeSphere(2), name="s4", color=color_child)
+            .add(makeSphere(1.5), name="s3", color=color_child)
+            .add(makeSphere(1), name="s2", color=color_child)
+            .constrain("board?h4", "s4", "Point")
+            .constrain("board?h3", "s3", "Point")
+            .constrain("board?h2", "s2", "Point")
+            .solve()
+        )
+        self.assertEqual(nhf.checks.pairwise_intersection(assembly), [])
+        bbox = assembly.toCompound().BoundingBox()
+        self.assertAlmostEqual(bbox.xlen, 15)
+        self.assertAlmostEqual(bbox.ylen, 15)
+        self.assertAlmostEqual(bbox.zlen, 7)
+
+    def test_tag_plane(self):
+        p4x, p4y = 5, 5
+        p3x, p3y = 0, 0
+        p2x, p2y = -5, 0
+        board = (
+            Cq.Workplane('XY')
+            .box(15, 15, 5)
+            .faces("<Z")
+            .workplane()
+            .pushPoints([(p4x, p4y)])
+            .hole(4, depth=2)
+            .pushPoints([(p3x, p3y)])
+            .hole(3, depth=1.5)
+            .pushPoints([(p2x, p2y)])
+            .hole(2, depth=1)
+        )
+        board.moveTo(p4x, p4y).tagPlane("h4")
+        board.moveTo(p3x, p3y).tagPlane("h3")
+        board.moveTo(p2x, p2y).tagPlane("h2")
+
+        def markedCylOf(r):
+            cyl = (
+                Cq.Workplane('XY')
+                .cylinder(radius=r, height=r)
+            )
+            cyl.faces("<Z").tag("mate")
+            return cyl
+
+        assembly = (
+            Cq.Assembly()
+            .add(board, name="board", color=color_parent)
+            .add(markedCylOf(2), name="c4", color=color_child)
+            .add(markedCylOf(1.5), name="c3", color=color_child)
+            .add(markedCylOf(1), name="c2", color=color_child)
+            .constrain("board?h4", "c4?mate", "Plane", param=0)
+            .constrain("board?h3", "c3?mate", "Plane", param=0)
+            .constrain("board?h2", "c2?mate", "Plane", param=0)
+            .solve()
+        )
+        self.assertEqual(nhf.checks.pairwise_intersection(assembly), [])
+        bbox = assembly.toCompound().BoundingBox()
+        self.assertAlmostEqual(bbox.xlen, 15)
+        self.assertAlmostEqual(bbox.ylen, 15)
+        self.assertAlmostEqual(bbox.zlen, 5)
 
 if __name__ == '__main__':
     unittest.main()
diff --git a/nhf/touhou/houjuu_nue/__init__.py b/nhf/touhou/houjuu_nue/__init__.py
index fe26116..9b5af33 100644
--- a/nhf/touhou/houjuu_nue/__init__.py
+++ b/nhf/touhou/houjuu_nue/__init__.py
@@ -34,10 +34,11 @@ import unittest
 import cadquery as Cq
 from nhf import Material, Role
 from nhf.build import Model, TargetKind, target
-from nhf.parts.joints import HirthJoint
+from nhf.parts.joints import HirthJoint, TorsionJoint
 from nhf.parts.handle import Handle
 import nhf.touhou.houjuu_nue.wing as MW
 import nhf.touhou.houjuu_nue.trident as MT
+import nhf.utils
 
 @dataclass
 class Parameters(Model):
@@ -89,6 +90,16 @@ class Parameters(Model):
 
     # Exterior radius of the wing root assembly
     wing_root_radius: float = 40
+    wing_root_wall_thickness: float = 8
+
+    shoulder_joint: TorsionJoint = field(default_factory=lambda: TorsionJoint(
+        radius_axle=8
+    ))
+
+    # Two holes on each side (top and bottom) are used to attach the shoulder
+    # joint. This governs the distance between these two holes
+    shoulder_attach_dist: float = 25
+    shoulder_attach_diam: float = 8
 
     """
     Heights for various wing joints, where the numbers start from the first joint.
@@ -164,7 +175,8 @@ class Parameters(Model):
         result.faces(">Y").tag("mount")
         plane = result.faces(">Y").workplane()
         for tag, x, y in self.harness_wing_base_pos:
-            conn = [(px + x, py + y) for px, py in self.hs_joint_harness_conn()]
+            conn = [(px + x, py + y) for px, py
+                    in self.hs_joint_harness_conn()]
             for i, (px, py) in enumerate(conn):
                 (
                     plane
@@ -250,7 +262,30 @@ class Parameters(Model):
         Generate the wing root which contains a Hirth joint at its base and a
         rectangular opening on its side, with the necessary interfaces.
         """
-        return MW.wing_root(joint=self.hs_hirth_joint)
+        return MW.wing_root(
+            joint=self.hs_hirth_joint,
+            shoulder_attach_dist=self.shoulder_attach_dist,
+            shoulder_attach_diam=self.shoulder_attach_diam,
+            wall_thickness=self.wing_root_wall_thickness,
+        )
+
+    @target(name="shoulder")
+    def shoulder_parent_joint(self) -> Cq.Assembly:
+        result = (
+            self.shoulder_joint.rider()
+            .copyWorkplane(Cq.Workplane(
+                'YZ', origin=Cq.Vector((100, 0, self.wing_root_wall_thickness))))
+            .rect(30, 10, centered=(True, False))
+            .extrude("next")
+            .copyWorkplane(Cq.Workplane(
+                'YX', origin=Cq.Vector((60, 0, self.wing_root_wall_thickness))))
+            .hole(self.shoulder_attach_diam)
+            .moveTo(0, self.shoulder_attach_dist)
+            .hole(self.shoulder_attach_diam)
+        )
+        result.moveTo(0, 0).tagPlane('conn0')
+        result.moveTo(0, self.shoulder_attach_dist).tagPlane('conn1')
+        return result
 
     def wing_r1_profile(self) -> Cq.Sketch:
         """
@@ -295,10 +330,10 @@ class Parameters(Model):
     #    Assemblies      #
     ######################
 
-    def trident_assembly(self):
+    def trident_assembly(self) -> Cq.Assembly:
         return MT.trident_assembly(self.trident_handle)
 
-    def harness_assembly(self):
+    def harness_assembly(self) -> Cq.Assembly:
         harness = self.harness()
         result = (
             Cq.Assembly()
@@ -317,7 +352,19 @@ class Parameters(Model):
         result.solve()
         return result
 
-    def wings_assembly(self):
+    def wing_r1_assembly(self) -> Cq.Assembly:
+        result = (
+            Cq.Assembly()
+            .add(self.wing_root(), name="r1")
+            .add(self.shoulder_parent_joint(), name="shoulder_parent_top",
+                 color=Material.RESIN_TRANSPARENT.color)
+            .constrain("r1/scaffold?conn_top0", "shoulder_parent_top?conn0", "Plane")
+            .constrain("r1/scaffold?conn_top1", "shoulder_parent_top?conn1", "Plane")
+            .solve()
+        )
+        return result
+
+    def wings_assembly(self) -> Cq.Assembly:
         """
         Assembly of harness with all the wings
         """
diff --git a/nhf/touhou/houjuu_nue/wing.py b/nhf/touhou/houjuu_nue/wing.py
index ea3a5c3..d9ab898 100644
--- a/nhf/touhou/houjuu_nue/wing.py
+++ b/nhf/touhou/houjuu_nue/wing.py
@@ -6,6 +6,7 @@ import math
 import cadquery as Cq
 from nhf import Material, Role
 from nhf.parts.joints import HirthJoint
+import nhf.utils
 
 def wing_root_profiles(
         base_sweep=150,
@@ -129,7 +130,8 @@ def wing_root_profiles(
 def wing_root(joint: HirthJoint,
               bolt_diam: int = 12,
               union_tol=1e-4,
-              attach_diam=8,
+              shoulder_attach_diam=8,
+              shoulder_attach_dist=25,
               conn_width=40,
               conn_height=100,
               wall_thickness=8) -> Cq.Assembly:
@@ -139,7 +141,7 @@ def wing_root(joint: HirthJoint,
     tip_centre = Cq.Vector((-150, 0, -80))
     attach_points = [
         (15, 0),
-        (40, 0),
+        (15 + shoulder_attach_dist, 0),
     ]
     root_profile, middle_profile, tip_profile = wing_root_profiles(
         conn_width=conn_width,
@@ -183,7 +185,7 @@ def wing_root(joint: HirthJoint,
             Cq.Workplane('bottom', origin=tip_centre + Cq.Vector((0, -50, 0)))
         )
         .pushPoints(attach_points)
-        .hole(attach_diam)
+        .hole(shoulder_attach_diam)
     )
     # Generate attach point tags
 
@@ -200,12 +202,8 @@ def wing_root(joint: HirthJoint,
             )
         )
         for i, (px, py) in enumerate(attach_points):
-            (
-                plane
-                .moveTo(px, py)
-                .eachpoint(Cq.Vertex.makeVertex(0, 0, 0))
-                .tag(f"conn_{side}{i}")
-            )
+            tag = f"conn_{side}{i}"
+            plane.moveTo(px, py).tagPlane(tag)
 
     result.faces("<Z").tag("base")
     result.faces(">X").tag("conn")
diff --git a/nhf/utils.py b/nhf/utils.py
new file mode 100644
index 0000000..7cc020d
--- /dev/null
+++ b/nhf/utils.py
@@ -0,0 +1,38 @@
+"""
+Marking utilities for `Cq.Workplane`
+
+Adds the functions to `Cq.Workplane`:
+1. `tagPoint`
+2. `tagPlane`
+"""
+import cadquery as Cq
+
+
+def tagPoint(self, tag: str):
+    """
+    Adds a vertex that can be used in `Point` constraints.
+    """
+    vertex = Cq.Vertex.makeVertex(0, 0, 0)
+    self.eachpoint(vertex.moved, useLocalCoordinates=True).tag(tag)
+
+Cq.Workplane.tagPoint = tagPoint
+
+
+def tagPlane(self, tag: str, axis='Z'):
+    """
+    Adds a phantom `Cq.Edge` in the given location which can be referenced in a
+    `Axis`, `Point`, or `Plane` constraint.
+    """
+    x, y, z = 0, 0, 0
+    if axis in ('z', 'Z'):
+        z = 1
+    elif axis in ('y', 'Y'):
+        y = 1
+    elif axis in ('x', 'X'):
+        x = 1
+    else:
+        assert False, "Axis must be one of x,y,z"
+    edge = Cq.Edge.makeLine((-x, -y, -z), (x, y, z))
+    self.eachpoint(edge.moved, useLocalCoordinates=True).tag(tag)
+
+Cq.Workplane.tagPlane = tagPlane