diff --git a/.github/workflows/build-ci-adjoint.yml b/.github/workflows/build-ci-adjoint.yml
index 4189fd3c5..494368d94 100644
--- a/.github/workflows/build-ci-adjoint.yml
+++ b/.github/workflows/build-ci-adjoint.yml
@@ -35,6 +35,12 @@ jobs:
         repository: NanoComp/libctl
         path: libctl-src
 
+    - name: Record libctl revision for the cache key
+      id: libctl
+      # Track libctl's checked-out commit so a new libctl revision invalidates
+      # the cached ~/local instead of restoring a stale (pre-mesh) build.
+      run: echo "sha=$(git -C libctl-src rev-parse HEAD)" >> "$GITHUB_OUTPUT"
+
     - name: Checkout harminv repository
       uses: actions/checkout@v6
       with:
@@ -52,7 +58,9 @@ jobs:
       id: deps-cache
       with:
         path: ~/local
-        key: deps-adjoint-${{ runner.os }}-${{ hashFiles('.github/workflows/build-ci-adjoint.yml') }}
+        # libctl's commit is part of the key so a new libctl revision rebuilds
+        # the cached deps instead of restoring a stale (pre-mesh) build.
+        key: deps-adjoint-${{ runner.os }}-libctl-${{ steps.libctl.outputs.sha }}-${{ hashFiles('.github/workflows/build-ci-adjoint.yml') }}
 
     - name: Build and install libctl
       if: steps.deps-cache.outputs.cache-hit != 'true'
diff --git a/.github/workflows/build-ci.yml b/.github/workflows/build-ci.yml
index 461024786..028093d2e 100644
--- a/.github/workflows/build-ci.yml
+++ b/.github/workflows/build-ci.yml
@@ -57,6 +57,12 @@ jobs:
         repository: NanoComp/libctl
         path: libctl-src
 
+    - name: Record libctl revision for the cache key
+      id: libctl
+      # Track libctl's checked-out commit so a new libctl revision invalidates
+      # the cached ~/local instead of restoring a stale (pre-mesh) build.
+      run: echo "sha=$(git -C libctl-src rev-parse HEAD)" >> "$GITHUB_OUTPUT"
+
     - name: Checkout harminv repository
       uses: actions/checkout@v6
       with:
@@ -80,7 +86,9 @@ jobs:
       id: deps-cache
       with:
         path: ~/local
-        key: deps-${{ runner.os }}-${{ hashFiles('.github/workflows/build-ci.yml') }}
+        # libctl's commit is part of the key so a new libctl revision rebuilds
+        # the cached deps instead of restoring a stale (pre-mesh) build.
+        key: deps-${{ runner.os }}-libctl-${{ steps.libctl.outputs.sha }}-${{ hashFiles('.github/workflows/build-ci.yml') }}
 
     - name: Build and install libctl
       if: steps.deps-cache.outputs.cache-hit != 'true'
diff --git a/python/geom.py b/python/geom.py
index 26f7b5270..116afc5f1 100755
--- a/python/geom.py
+++ b/python/geom.py
@@ -1343,6 +1343,171 @@ def __init__(
         super().__init__(center=center, **kwargs)
 
 
+class Mesh(GeometricObject):
+    """
+    Triangulated 3D mesh object.
+
+    The mesh must be a closed (watertight) manifold with consistently
+    oriented triangles. Open meshes will produce a warning and may
+    give incorrect results for subpixel smoothing.
+    """
+
+    def __init__(self, vertices, triangles, center=None, **kwargs):
+        """
+        Construct a `Mesh`.
+
+        + **`vertices` [list of `Vector3`]** — Vertex positions.
+
+        + **`triangles` [list of tuples of 3 ints]** — Triangle vertex indices
+          (0-based). Each tuple defines a face with outward normal determined
+          by the right-hand rule.
+
+        + **`center` [`Vector3`, optional]** — If specified, the mesh is
+          translated so its centroid is at this position. If omitted, the
+          centroid of the vertices is used.
+        """
+        self.vertices = [Vector3(*v) for v in vertices]
+        self.triangles = [tuple(t) for t in triangles]
+
+        centroid = sum(self.vertices, Vector3(0)) * (1.0 / len(self.vertices))
+        if center is not None:
+            center = Vector3(*center)
+            shift = center - centroid
+            self.vertices = [v + shift for v in self.vertices]
+        else:
+            center = centroid
+
+        super().__init__(center=center, **kwargs)
+
+    @classmethod
+    def from_stl(cls, filename, material=None, center=None, scale=1.0):
+        """Load mesh from STL file (binary or ASCII)."""
+        vertices, triangles = _parse_stl(filename, scale)
+        kwargs = {}
+        if material is not None:
+            kwargs["material"] = material
+        return cls(vertices=vertices, triangles=triangles, center=center, **kwargs)
+
+    @classmethod
+    def from_obj(cls, filename, material=None, center=None, scale=1.0):
+        """Load mesh from Wavefront OBJ file."""
+        vertices, triangles = _parse_obj(filename, scale)
+        kwargs = {}
+        if material is not None:
+            kwargs["material"] = material
+        return cls(vertices=vertices, triangles=triangles, center=center, **kwargs)
+
+
+def _parse_stl(filename, scale=1.0):
+    """Parse an STL file (binary or ASCII) and return (vertices, triangles).
+
+    STL stores per-face vertices without sharing, so vertex welding is
+    performed to build a shared vertex array and triangle index array.
+    """
+    import struct
+
+    with open(filename, "rb") as f:
+        header = f.read(80)
+        data = f.read()
+
+    # Try binary format first
+    is_binary = True
+    try:
+        num_tris = struct.unpack("<I", data[:4])[0]
+        expected_size = num_tris * 50 + 4
+        if len(data) != expected_size:
+            is_binary = False
+    except struct.error:
+        is_binary = False
+
+    raw_verts = []
+    if is_binary:
+        offset = 4
+        for _ in range(num_tris):
+            # Skip normal (12 bytes), read 3 vertices (36 bytes), skip attr (2 bytes)
+            offset += 12  # normal
+            for _ in range(3):
+                x, y, z = struct.unpack("<fff", data[offset : offset + 12])
+                raw_verts.append((x * scale, y * scale, z * scale))
+                offset += 12
+            offset += 2  # attribute byte count
+    else:
+        # ASCII format
+        text = (header + data).decode("ascii", errors="ignore")
+        import re
+
+        for m in re.finditer(
+            r"vertex\s+([-+eE.\d]+)\s+([-+eE.\d]+)\s+([-+eE.\d]+)", text
+        ):
+            x, y, z = float(m.group(1)), float(m.group(2)), float(m.group(3))
+            raw_verts.append((x * scale, y * scale, z * scale))
+
+    # Vertex welding: merge vertices within tolerance using a spatial hash
+    tol = 1e-8
+    inv_tol = 1.0 / tol if tol > 0 else 1e8
+    vertex_map = {}
+    vertices = []
+    triangles = []
+
+    def _quantize(v):
+        return (
+            int(round(v[0] * inv_tol)),
+            int(round(v[1] * inv_tol)),
+            int(round(v[2] * inv_tol)),
+        )
+
+    for i in range(0, len(raw_verts), 3):
+        tri = []
+        for j in range(3):
+            v = raw_verts[i + j]
+            key = _quantize(v)
+            if key not in vertex_map:
+                vertex_map[key] = len(vertices)
+                vertices.append(Vector3(v[0], v[1], v[2]))
+            tri.append(vertex_map[key])
+        triangles.append(tuple(tri))
+
+    return vertices, triangles
+
+
+def _parse_obj(filename, scale=1.0):
+    """Parse a Wavefront OBJ file and return (vertices, triangles).
+
+    Quad faces are automatically triangulated. OBJ is 1-indexed;
+    indices are converted to 0-indexed.
+    """
+    vertices = []
+    triangles = []
+
+    with open(filename) as f:
+        for line in f:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            parts = line.split()
+            if parts[0] == "v" and len(parts) >= 4:
+                x, y, z = float(parts[1]), float(parts[2]), float(parts[3])
+                vertices.append(Vector3(x * scale, y * scale, z * scale))
+            elif parts[0] == "f":
+                # OBJ face indices may include texture/normal refs: v/vt/vn
+                face_indices = []
+                for p in parts[1:]:
+                    idx = int(p.split("/")[0])
+                    # OBJ is 1-indexed; handle negative indices
+                    if idx < 0:
+                        idx = len(vertices) + idx
+                    else:
+                        idx -= 1
+                    face_indices.append(idx)
+                # Triangulate (fan triangulation for convex polygons)
+                for i in range(1, len(face_indices) - 1):
+                    triangles.append(
+                        (face_indices[0], face_indices[i], face_indices[i + 1])
+                    )
+
+    return vertices, triangles
+
+
 class Matrix:
     """
     The `Matrix` class represents a 3x3 matrix with c1, c2, and c3 as its columns.
diff --git a/python/meep.i b/python/meep.i
index 47fa6df8b..17e06c956 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -1631,6 +1631,7 @@ PyObject *_get_array_slice_dimensions(meep::fields *f, const meep::volume &where
         MaterialGrid,
         Matrix,
         Medium,
+        Mesh,
         MultilevelAtom,
         NoisyDrudeSusceptibility,
         NoisyLorentzianSusceptibility,
diff --git a/python/tests/test_mesh.py b/python/tests/test_mesh.py
new file mode 100644
index 000000000..c8ebc3fb4
--- /dev/null
+++ b/python/tests/test_mesh.py
@@ -0,0 +1,207 @@
+"""Integration tests for the Mesh geometric object type."""
+
+import meep as mp
+import numpy as np
+import struct
+import tempfile
+import os
+import unittest
+
+
+def make_cube_mesh(material, center=mp.Vector3(0, 0, 0)):
+    """Create a unit cube mesh centered at the given point."""
+    verts = [
+        (-0.5, -0.5, -0.5),
+        (0.5, -0.5, -0.5),
+        (0.5, 0.5, -0.5),
+        (-0.5, 0.5, -0.5),
+        (-0.5, -0.5, 0.5),
+        (0.5, -0.5, 0.5),
+        (0.5, 0.5, 0.5),
+        (-0.5, 0.5, 0.5),
+    ]
+    tris = [
+        (0, 2, 1),
+        (0, 3, 2),  # -z
+        (4, 5, 6),
+        (4, 6, 7),  # +z
+        (0, 1, 5),
+        (0, 5, 4),  # -y
+        (2, 3, 7),
+        (2, 7, 6),  # +y
+        (0, 4, 7),
+        (0, 7, 3),  # -x
+        (1, 2, 6),
+        (1, 6, 5),  # +x
+    ]
+    return mp.Mesh(vertices=verts, triangles=tris, center=center, material=material)
+
+
+def write_cube_stl(path):
+    """Write a unit cube as a binary STL file."""
+    faces = [
+        ((-0.5, -0.5, -0.5), (0.5, 0.5, -0.5), (0.5, -0.5, -0.5)),
+        ((-0.5, -0.5, -0.5), (-0.5, 0.5, -0.5), (0.5, 0.5, -0.5)),
+        ((-0.5, -0.5, 0.5), (0.5, -0.5, 0.5), (0.5, 0.5, 0.5)),
+        ((-0.5, -0.5, 0.5), (0.5, 0.5, 0.5), (-0.5, 0.5, 0.5)),
+        ((-0.5, -0.5, -0.5), (0.5, -0.5, -0.5), (0.5, -0.5, 0.5)),
+        ((-0.5, -0.5, -0.5), (0.5, -0.5, 0.5), (-0.5, -0.5, 0.5)),
+        ((0.5, 0.5, -0.5), (-0.5, 0.5, -0.5), (-0.5, 0.5, 0.5)),
+        ((0.5, 0.5, -0.5), (-0.5, 0.5, 0.5), (0.5, 0.5, 0.5)),
+        ((-0.5, -0.5, -0.5), (-0.5, -0.5, 0.5), (-0.5, 0.5, 0.5)),
+        ((-0.5, -0.5, -0.5), (-0.5, 0.5, 0.5), (-0.5, 0.5, -0.5)),
+        ((0.5, -0.5, -0.5), (0.5, 0.5, -0.5), (0.5, 0.5, 0.5)),
+        ((0.5, -0.5, -0.5), (0.5, 0.5, 0.5), (0.5, -0.5, 0.5)),
+    ]
+    with open(path, "wb") as f:
+        f.write(b"\x00" * 80)
+        f.write(struct.pack("<I", len(faces)))
+        for v0, v1, v2 in faces:
+            f.write(struct.pack("<fff", 0, 0, 0))  # normal (ignored)
+            for v in (v0, v1, v2):
+                f.write(struct.pack("<fff", *v))
+            f.write(struct.pack("<H", 0))
+
+
+class TestMesh(unittest.TestCase):
+    def test_stl_import(self):
+        """STL import: correct vertex/face counts and closure."""
+        with tempfile.NamedTemporaryFile(suffix=".stl", delete=False) as f:
+            stl_path = f.name
+            write_cube_stl(stl_path)
+        try:
+            mesh = mp.Mesh.from_stl(stl_path, material=mp.Medium(epsilon=12))
+            self.assertEqual(len(mesh.vertices), 8)
+            self.assertEqual(len(mesh.triangles), 12)
+        finally:
+            os.unlink(stl_path)
+
+    def test_point_in_mesh(self):
+        """Points inside and outside a cube mesh."""
+        mat = mp.Medium(epsilon=12)
+        cell = mp.Vector3(2, 2, 2)
+        cube = make_cube_mesh(mat)
+        sim = mp.Simulation(
+            cell_size=cell,
+            geometry=[cube],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=False,
+        )
+        sim.init_sim()
+
+        # Sample single points via get_epsilon_grid
+        def eps_at(x, y, z):
+            return np.real(
+                sim.get_epsilon_grid(np.array([x]), np.array([y]), np.array([z]))
+            ).item()
+
+        self.assertAlmostEqual(eps_at(0, 0, 0), 12.0, places=1)  # inside
+        self.assertAlmostEqual(eps_at(0.3, 0.3, 0.3), 12.0, places=1)
+        self.assertAlmostEqual(eps_at(0.8, 0, 0), 1.0, places=1)  # outside
+        self.assertAlmostEqual(eps_at(0, 0.8, 0), 1.0, places=1)
+
+    def test_volume(self):
+        """Mesh cube volume matches analytic value via grid sampling."""
+        mat = mp.Medium(epsilon=12)
+        cube = make_cube_mesh(mat)
+        sim = mp.Simulation(
+            cell_size=mp.Vector3(2, 2, 2),
+            geometry=[cube],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=False,
+        )
+        sim.init_sim()
+
+        # Sample a uniform grid and count points inside the cube
+        n = 40
+        xtics = np.linspace(-1, 1, n)
+        ytics = np.linspace(-1, 1, n)
+        ztics = np.linspace(-1, 1, n)
+        eps = np.real(sim.get_epsilon_grid(xtics, ytics, ztics))
+        inside = np.sum(eps > 6)
+        total = eps.size
+        fill_fraction = inside / total
+        expected = 1.0 / 8.0  # unit cube in 2x2x2 cell
+        self.assertAlmostEqual(fill_fraction, expected, delta=0.02)
+
+    def test_eps_averaging(self):
+        """Subpixel smoothing produces intermediate epsilon at interfaces."""
+        mat = mp.Medium(epsilon=12)
+        cell = mp.Vector3(2, 2, 2)
+        cube = make_cube_mesh(mat)
+
+        sim_off = mp.Simulation(
+            cell_size=cell,
+            geometry=[cube],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=False,
+        )
+        sim_off.init_sim()
+        eps_off = np.real(
+            sim_off.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)
+        )
+
+        sim_on = mp.Simulation(
+            cell_size=cell,
+            geometry=[cube],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=True,
+        )
+        sim_on.init_sim()
+        eps_on = np.real(
+            sim_on.get_array(center=mp.Vector3(), size=cell, component=mp.Dielectric)
+        )
+
+        unique_off = len(np.unique(np.round(eps_off, 2)))
+        unique_on = len(np.unique(np.round(eps_on, 2)))
+        self.assertGreater(
+            unique_on,
+            unique_off,
+            "eps_averaging should produce more unique epsilon values",
+        )
+
+    def test_mesh_vs_block(self):
+        """Cube mesh epsilon grid matches native Block."""
+        mat = mp.Medium(epsilon=12)
+        cell = mp.Vector3(2, 2, 2)
+        n = 30
+        xtics = np.linspace(-1, 1, n)
+        ytics = np.linspace(-1, 1, n)
+        ztics = np.array([0.0])
+
+        cube_mesh = make_cube_mesh(mat)
+        sim_mesh = mp.Simulation(
+            cell_size=cell,
+            geometry=[cube_mesh],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=False,
+        )
+        sim_mesh.init_sim()
+        eps_mesh = np.real(sim_mesh.get_epsilon_grid(xtics, ytics, ztics)).flatten()
+
+        block = mp.Block(size=mp.Vector3(1, 1, 1), material=mat)
+        sim_block = mp.Simulation(
+            cell_size=cell,
+            geometry=[block],
+            resolution=4,
+            dimensions=3,
+            eps_averaging=False,
+        )
+        sim_block.init_sim()
+        eps_block = np.real(sim_block.get_epsilon_grid(xtics, ytics, ztics)).flatten()
+
+        mismatches = np.sum(np.abs(eps_mesh - eps_block) > 0.5)
+        self.assertEqual(
+            mismatches,
+            0,
+            f"Cube mesh should match Block exactly, got {mismatches} mismatches",
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/python/typemap_utils.cpp b/python/typemap_utils.cpp
index 5352bf1eb..d4b827c51 100644
--- a/python/typemap_utils.cpp
+++ b/python/typemap_utils.cpp
@@ -1038,6 +1038,63 @@ static int pyprism_to_prism(PyObject *py_prism, geometric_object *p) {
   ;
 }
 
+static int pymesh_to_mesh(PyObject *py_mesh, geometric_object *o) {
+  SWIG_PYTHON_THREAD_SCOPED_BLOCK;
+  material_type material;
+  vector3 center;
+
+  if (!get_attr_material(py_mesh, &material) || !get_attr_v3(py_mesh, &center, "center")) {
+    return 0;
+  }
+
+  PyObject *py_verts = PyObject_GetAttrString(py_mesh, "vertices");
+  if (!py_verts || !PyList_Check(py_verts)) {
+    meep::abort("Expected Mesh.vertices to be a list\n");
+    return 0;
+  }
+
+  int num_vertices = PyList_Size(py_verts);
+  vector3 *vertices = new vector3[num_vertices];
+  for (Py_ssize_t i = 0; i < num_vertices; ++i) {
+    if (!pyv3_to_v3(PyList_GetItem(py_verts, i), &vertices[i])) {
+      delete[] vertices;
+      Py_DECREF(py_verts);
+      return 0;
+    }
+  }
+  Py_DECREF(py_verts);
+
+  PyObject *py_tris = PyObject_GetAttrString(py_mesh, "triangles");
+  if (!py_tris || !PyList_Check(py_tris)) {
+    delete[] vertices;
+    meep::abort("Expected Mesh.triangles to be a list\n");
+    return 0;
+  }
+
+  int num_triangles = PyList_Size(py_tris);
+  int *triangles = new int[3 * num_triangles];
+  for (Py_ssize_t i = 0; i < num_triangles; ++i) {
+    PyObject *tri = PyList_GetItem(py_tris, i);
+    if (!PyTuple_Check(tri) || PyTuple_Size(tri) != 3) {
+      delete[] vertices;
+      delete[] triangles;
+      Py_DECREF(py_tris);
+      meep::abort("Each triangle must be a tuple of 3 integers\n");
+      return 0;
+    }
+    for (int j = 0; j < 3; ++j)
+      triangles[3 * i + j] = (int)PyLong_AsLong(PyTuple_GetItem(tri, j));
+  }
+  Py_DECREF(py_tris);
+
+  *o = make_mesh_with_center(material, center, vertices, num_vertices, triangles, num_triangles);
+
+  delete[] vertices;
+  delete[] triangles;
+
+  return 1;
+}
+
 static int py_gobj_to_gobj(PyObject *po, geometric_object *o) {
   SWIG_PYTHON_THREAD_SCOPED_BLOCK;
   int success = 0;
@@ -1050,6 +1107,7 @@ static int py_gobj_to_gobj(PyObject *po, geometric_object *o) {
   else if (go_type == "Block") { success = pyblock_to_block(po, o); }
   else if (go_type == "Ellipsoid") { success = pyellipsoid_to_ellipsoid(po, o); }
   else if (go_type == "Prism") { success = pyprism_to_prism(po, o); }
+  else if (go_type == "Mesh") { success = pymesh_to_mesh(po, o); }
   else {
     PyErr_Format(PyExc_TypeError, "Error: %s is not a valid GeometricObject type", go_type.c_str());
     return 0;