From d70754e1e7540b5044f266dd803b67a903892e1b Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 18 Jun 2026 08:17:48 +0000
Subject: [PATCH 1/4] Initial plan


From 9e7e0daebe9afcee1d923e2214b61c221a7fc457 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 18 Jun 2026 08:27:15 +0000
Subject: [PATCH 2/4] Add tests for 12 planned coverage areas

---
 tests/conftest.py                        | 360 +++++++++++++++++++++++
 tests/test_01_fea_add_input.py           |  26 ++
 tests/test_02_fea_add_state.py           |  29 ++
 tests/test_03_fea_add_output.py          |  27 ++
 tests/test_04_fea_add_field_output.py    |  21 ++
 tests/test_05_fea_boundary_conditions.py |  12 +
 tests/test_06_utils_find_node_indices.py |  27 ++
 tests/test_07_utils_locate_dofs.py       |  43 +++
 tests/test_08_utils_import_mesh.py       |  35 +++
 tests/test_09_csdl_fea_model.py          |  25 ++
 tests/test_10_output_operation.py        |  29 ++
 tests/test_11_state_operation.py         |  72 +++++
 tests/test_12_smoke_imports.py           |   6 +
 13 files changed, 712 insertions(+)
 create mode 100644 tests/conftest.py
 create mode 100644 tests/test_01_fea_add_input.py
 create mode 100644 tests/test_02_fea_add_state.py
 create mode 100644 tests/test_03_fea_add_output.py
 create mode 100644 tests/test_04_fea_add_field_output.py
 create mode 100644 tests/test_05_fea_boundary_conditions.py
 create mode 100644 tests/test_06_utils_find_node_indices.py
 create mode 100644 tests/test_07_utils_locate_dofs.py
 create mode 100644 tests/test_08_utils_import_mesh.py
 create mode 100644 tests/test_09_csdl_fea_model.py
 create mode 100644 tests/test_10_output_operation.py
 create mode 100644 tests/test_11_state_operation.py
 create mode 100644 tests/test_12_smoke_imports.py

diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..c8c957b
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,360 @@
+import importlib
+import sys
+import types
+import numpy as np
+import pytest
+
+
+class _ParameterStore(dict):
+    def declare(self, name, default=None, types=None):
+        if name not in self:
+            self[name] = default
+
+
+class _BaseModel:
+    def __init__(self, **kwargs):
+        self.parameters = _ParameterStore(kwargs)
+        self._added = []
+
+    def add(self, model, name=None):
+        self._added.append((name, model))
+
+    def declare_variable(self, name, shape=None, val=0.0):
+        if isinstance(val, np.ndarray):
+            return val
+        if shape is None:
+            return np.array(val)
+        return np.full(shape, val)
+
+    def register_output(self, name, value):
+        return value
+
+    def print_var(self, value):
+        return None
+
+
+class _BaseOperation:
+    def __init__(self, **kwargs):
+        self.parameters = _ParameterStore(kwargs)
+        self._inputs = {}
+        self._outputs = {}
+
+    def add_input(self, name, shape=None):
+        self._inputs[name] = shape
+
+    def add_output(self, name, shape=None):
+        self._outputs[name] = shape
+
+    def declare_derivatives(self, of, wrt):
+        return None
+
+
+@pytest.fixture(scope="session", autouse=True)
+def dependency_stubs():
+    # csdl stubs
+    csdl = types.ModuleType("csdl")
+    csdl.Model = _BaseModel
+    csdl.CustomExplicitOperation = _BaseOperation
+    csdl.CustomImplicitOperation = _BaseOperation
+    csdl.custom = lambda *args, op=None: np.zeros((1,))
+    sys.modules["csdl"] = csdl
+
+    # matplotlib stub
+    matplotlib = types.ModuleType("matplotlib")
+    pyplot = types.ModuleType("matplotlib.pyplot")
+    pyplot.figure = lambda *a, **k: None
+    pyplot.plot = lambda *a, **k: None
+    pyplot.show = lambda *a, **k: None
+    matplotlib.pyplot = pyplot
+    sys.modules["matplotlib"] = matplotlib
+    sys.modules["matplotlib.pyplot"] = pyplot
+
+    # mpi4py stubs
+    mpi4py = types.ModuleType("mpi4py")
+    mpi4py.MPI = types.SimpleNamespace(
+        COMM_WORLD=types.SimpleNamespace(allreduce=lambda x, op=None: x),
+        SUM="SUM",
+    )
+    sys.modules["mpi4py"] = mpi4py
+
+    # petsc4py stubs
+    petsc4py = types.ModuleType("petsc4py")
+
+    class _DummyKSP:
+        def create(self, comm=None):
+            return self
+
+        def setOperators(self, *args, **kwargs):
+            return None
+
+        def setType(self, *args, **kwargs):
+            return None
+
+        def getPC(self):
+            return types.SimpleNamespace(
+                setType=lambda *a, **k: None,
+                setFactorSolverType=lambda *a, **k: None,
+                getASMSubKSP=lambda: [types.SimpleNamespace(
+                    setType=lambda *a, **k: None,
+                    getPC=lambda: types.SimpleNamespace(setType=lambda *a, **k: None),
+                    setTolerances=lambda *a, **k: None,
+                )],
+            )
+
+        def setFromOptions(self):
+            return None
+
+        def setUp(self):
+            return None
+
+        def solve(self, b, x):
+            return None
+
+        def setConvergenceHistory(self):
+            return None
+
+        def getConvergenceHistory(self):
+            return []
+
+        def setTolerances(self, *args, **kwargs):
+            return None
+
+    class _DummySNES:
+        def create(self):
+            return self
+
+        def setTolerances(self, **kwargs):
+            return None
+
+        def getKSP(self):
+            return _DummyKSP()
+
+        def setFunction(self, *args, **kwargs):
+            return None
+
+        def setJacobian(self, *args, **kwargs):
+            return None
+
+        def setFromOptions(self):
+            return None
+
+        def solve(self, *args, **kwargs):
+            return None
+
+        def getConvergedReason(self):
+            return 1
+
+    petsc4py.PETSc = types.SimpleNamespace(
+        InsertMode=types.SimpleNamespace(ADD_VALUES=1, INSERT=2),
+        ScatterMode=types.SimpleNamespace(REVERSE=1, FORWARD=2),
+        KSP=_DummyKSP,
+        SNES=_DummySNES,
+        Options=lambda: {},
+        Mat=lambda comm=None: object(),
+        Vec=lambda: types.SimpleNamespace(create=lambda: types.SimpleNamespace(setSizes=lambda *a: None, setUp=lambda: None, assemble=lambda: None, getArray=lambda: np.array([]))),
+    )
+    sys.modules["petsc4py"] = petsc4py
+
+    # scipy stubs (only if unavailable)
+    if "scipy" not in sys.modules:
+        scipy = types.ModuleType("scipy")
+        sparse = types.ModuleType("scipy.sparse")
+        spatial = types.ModuleType("scipy.spatial")
+
+        class _FakeKDTree:
+            def __init__(self, coords):
+                self.coords = np.asarray(coords)
+
+            def query(self, nodes):
+                nodes = np.asarray(nodes)
+                dists = []
+                idxs = []
+                for node in nodes:
+                    diff = self.coords - node
+                    ds = np.sqrt(np.sum(diff * diff, axis=1))
+                    i = int(np.argmin(ds))
+                    dists.append(ds[i])
+                    idxs.append(i)
+                return np.array(dists), np.array(idxs)
+
+        sparse.csr_matrix = lambda *args, **kwargs: types.SimpleNamespace(tocoo=lambda: object())
+        spatial.KDTree = _FakeKDTree
+        scipy.sparse = sparse
+        scipy.spatial = spatial
+        sys.modules["scipy"] = scipy
+        sys.modules["scipy.sparse"] = sparse
+        sys.modules["scipy.spatial"] = spatial
+
+    # ufl stubs
+    ufl = types.ModuleType("ufl")
+    _f = lambda *a, **k: (a, k)
+    for name in [
+        "Identity", "dot", "derivative", "TestFunction", "TrialFunction", "inner", "ds", "dS", "dx",
+        "grad", "inv", "as_vector", "sqrt", "conditional", "lt", "det", "Measure", "exp", "tr",
+        "CellDiameter", "div", "SpatialCoordinate", "FacetNormal",
+    ]:
+        setattr(ufl, name, _f)
+    sys.modules["ufl"] = ufl
+
+    # dolfinx stubs
+    dolfinx = types.ModuleType("dolfinx")
+    dolfinx.log = types.SimpleNamespace(LogLevel=types.SimpleNamespace(INFO=1), set_log_level=lambda lvl: None)
+
+    io = types.ModuleType("dolfinx.io")
+
+    class _XDMFFile:
+        def __init__(self, *args, **kwargs):
+            self.path = args[1] if len(args) > 1 else ""
+
+        def __enter__(self):
+            return self
+
+        def __exit__(self, exc_type, exc, tb):
+            return False
+
+        def read_mesh(self, name="Grid"):
+            return types.SimpleNamespace(topology=types.SimpleNamespace(dim=2, create_connectivity=lambda *a, **k: None))
+
+        def read_meshtags(self, mesh, name="Grid"):
+            return {"name": name, "path": self.path}
+
+        def write_mesh(self, mesh):
+            return None
+
+        def write_function(self, func, step):
+            return None
+
+    io.XDMFFile = _XDMFFile
+
+    mesh = types.ModuleType("dolfinx.mesh")
+    mesh.create_unit_square = lambda *a, **k: object()
+    mesh.create_rectangle = lambda *a, **k: object()
+    mesh.create_interval = lambda *a, **k: object()
+    mesh.locate_entities_boundary = lambda *a, **k: np.array([], dtype=np.int32)
+    mesh.locate_entities = lambda *a, **k: np.array([], dtype=np.int32)
+    mesh.meshtags = lambda *a, **k: object()
+
+    cpp = types.ModuleType("dolfinx.cpp")
+    cpp_mesh = types.ModuleType("dolfinx.cpp.mesh")
+    cpp_mesh.CellType = types.SimpleNamespace(quadrilateral="quadrilateral")
+    cpp_mesh.h = lambda mesh, tdim, cells: np.array([])
+    cpp.mesh = cpp_mesh
+
+    fem = types.ModuleType("dolfinx.fem")
+
+    class _DummyVector:
+        def __init__(self):
+            self.array = np.array([], dtype=float)
+
+        def getArray(self):
+            return self.array
+
+        def set(self, value):
+            self.array[:] = value
+
+        def assemble(self):
+            return None
+
+        def ghostUpdate(self, *args, **kwargs):
+            return None
+
+        def pointwiseDivide(self, b, a):
+            return None
+
+        def localForm(self):
+            class _Ctx:
+                def __enter__(self_non):
+                    return types.SimpleNamespace(set=lambda v: None)
+
+                def __exit__(self_non, exc_type, exc, tb):
+                    return False
+            return _Ctx()
+
+        def copy(self, other):
+            return None
+
+    class _DummyFunction:
+        def __init__(self, function_space=None):
+            self.function_space = function_space or types.SimpleNamespace(num_sub_spaces=1)
+            self.vector = _DummyVector()
+            self.x = types.SimpleNamespace(array=np.array([], dtype=float))
+
+        def rename(self, *args, **kwargs):
+            return None
+
+        def interpolate(self, *args, **kwargs):
+            return None
+
+        def split(self):
+            return self, None
+
+        def sub(self, idx):
+            return types.SimpleNamespace(collapse=lambda: self)
+
+    fem.form = lambda x: x
+    fem.assemble_scalar = lambda x: 0.0
+    fem.Function = _DummyFunction
+    fem.FunctionSpace = lambda mesh, desc: types.SimpleNamespace(mesh=mesh)
+    fem.VectorFunctionSpace = lambda *a, **k: types.SimpleNamespace(mesh=a[0] if a else None)
+    fem.dirichletbc = lambda *a, **k: (a, k)
+    fem.locate_dofs_geometrical = lambda *a, **k: np.array([], dtype=np.int32)
+    fem.locate_dofs_topological = lambda *a, **k: np.array([], dtype=np.int32)
+    fem.Constant = lambda mesh, v: v
+    fem.set_bc = lambda *a, **k: None
+
+    fem_petsc = types.ModuleType("dolfinx.fem.petsc")
+    fem_petsc.assemble_vector = lambda *a, **k: types.SimpleNamespace(array=np.array([]), ghostUpdate=lambda *aa, **kk: None)
+    fem_petsc.assemble_matrix = lambda *a, **k: types.SimpleNamespace(
+        assemble=lambda: None,
+        copy=lambda: types.SimpleNamespace(assemble=lambda: None, convert=lambda fmt: types.SimpleNamespace(getDenseArray=lambda: np.array([[]]))),
+        transpose=lambda *aa, **kk: types.SimpleNamespace(),
+        getValuesCSR=lambda: (np.array([0]), np.array([0]), np.array([0.0])),
+        size=(0, 0),
+        getSizes=lambda: (0, 0),
+        getComm=lambda: None,
+        multTranspose=lambda *aa, **kk: None,
+    )
+    fem_petsc.NonlinearProblem = lambda *a, **k: object()
+    fem_petsc.apply_lifting = lambda *a, **k: None
+    fem_petsc.set_bc = lambda *a, **k: None
+    fem_petsc.create_matrix = lambda *a, **k: object()
+    fem_petsc._assemble_matrix_mat = lambda *a, **k: object()
+
+    nls = types.ModuleType("dolfinx.nls")
+    nls_petsc = types.ModuleType("dolfinx.nls.petsc")
+    nls_petsc.NewtonSolver = lambda *a, **k: types.SimpleNamespace(solve=lambda *aa, **kk: None)
+
+    la = types.ModuleType("dolfinx.la")
+    la.create_petsc_vector = lambda *a, **k: types.SimpleNamespace(
+        localForm=lambda: types.SimpleNamespace(__enter__=lambda self: types.SimpleNamespace(set=lambda v: None), __exit__=lambda self, exc_type, exc, tb: False)
+    )
+
+    dolfinx.io = io
+    dolfinx.mesh = mesh
+    dolfinx.cpp = cpp
+    dolfinx.fem = fem
+    dolfinx.nls = nls
+    dolfinx.la = la
+
+    sys.modules["dolfinx"] = dolfinx
+    sys.modules["dolfinx.io"] = io
+    sys.modules["dolfinx.mesh"] = mesh
+    sys.modules["dolfinx.cpp"] = cpp
+    sys.modules["dolfinx.cpp.mesh"] = cpp_mesh
+    sys.modules["dolfinx.fem"] = fem
+    sys.modules["dolfinx.fem.petsc"] = fem_petsc
+    sys.modules["dolfinx.nls"] = nls
+    sys.modules["dolfinx.nls.petsc"] = nls_petsc
+    sys.modules["dolfinx.la"] = la
+
+    yield
+
+
+@pytest.fixture
+def fresh_import():
+    def _fresh(module_name):
+        if module_name in sys.modules:
+            del sys.modules[module_name]
+        return importlib.import_module(module_name)
+
+    return _fresh
diff --git a/tests/test_01_fea_add_input.py b/tests/test_01_fea_add_input.py
new file mode 100644
index 0000000..8e99024
--- /dev/null
+++ b/tests/test_01_fea_add_input.py
@@ -0,0 +1,26 @@
+import types
+import numpy as np
+import pytest
+
+
+def test_add_input_rejects_duplicate_name(fresh_import):
+    fea_mod = fresh_import("femo.fea.fea_dolfinx")
+    fea = fea_mod.FEA(mesh=object())
+
+    class DummyFunction:
+        def __init__(self):
+            self.function_space = object()
+            self.x = types.SimpleNamespace(array=np.zeros(3))
+
+        def rename(self, *args, **kwargs):
+            return None
+
+    fea_mod.getFuncArray = lambda f: np.zeros(3)
+
+    fn = DummyFunction()
+    fea.add_input("density", fn, init_val=2.5)
+    assert "density" in fea.inputs_dict
+    assert fea.inputs_dict["density"]["shape"] == 3
+
+    with pytest.raises(ValueError):
+        fea.add_input("density", DummyFunction())
diff --git a/tests/test_02_fea_add_state.py b/tests/test_02_fea_add_state.py
new file mode 100644
index 0000000..e4b9464
--- /dev/null
+++ b/tests/test_02_fea_add_state.py
@@ -0,0 +1,29 @@
+import types
+import numpy as np
+
+
+def test_add_state_registers_expected_keys(fresh_import):
+    fea_mod = fresh_import("femo.fea.fea_dolfinx")
+    fea = fea_mod.FEA(mesh=object())
+
+    class DummyFunction:
+        def __init__(self):
+            self.function_space = object()
+
+        def rename(self, *args, **kwargs):
+            return None
+
+    fea_mod.getFuncArray = lambda f: np.zeros(4)
+    fea_mod.Function = lambda fs: types.SimpleNamespace(function_space=fs)
+
+    residual = object()
+    fn = DummyFunction()
+    fea.add_state("u", fn, residual, arguments=["rho"], dR_du="dRdu", dR_df_list=["dRdrho"])
+
+    state = fea.states_dict["u"]
+    assert state["function"] is fn
+    assert state["residual_form"] is residual
+    assert state["shape"] == 4
+    assert state["dR_du"] == "dRdu"
+    assert state["dR_df_list"] == ["dRdrho"]
+    assert state["arguments"] == ["rho"]
diff --git a/tests/test_03_fea_add_output.py b/tests/test_03_fea_add_output.py
new file mode 100644
index 0000000..6753cd4
--- /dev/null
+++ b/tests/test_03_fea_add_output.py
@@ -0,0 +1,27 @@
+import types
+import numpy as np
+
+
+def test_add_output_handles_scalar_and_field_shapes(fresh_import):
+    fea_mod = fresh_import("femo.fea.fea_dolfinx")
+    fea = fea_mod.FEA(mesh=object())
+
+    fea_mod.getFormArray = lambda form: np.zeros(5)
+    fea_mod.derivative = lambda form, fn: ("d", form, fn)
+
+    input_fn = types.SimpleNamespace(function_space=object(), x=types.SimpleNamespace(array=np.zeros(2)), rename=lambda *a, **k: None)
+    state_fn = types.SimpleNamespace(function_space=object(), rename=lambda *a, **k: None)
+
+    fea_mod.getFuncArray = lambda f: np.zeros(2)
+    fea.add_input("rho", input_fn)
+    fea.states_dict["u"] = {"function": state_fn}
+
+    fea.add_output("compliance", "scalar", form="f1", arguments=["rho", "u"])
+    scalar = fea.outputs_dict["compliance"]
+    assert scalar["shape"] == 1
+    assert len(scalar["partials"]) == 2
+
+    fea.add_output("stress", "field", form="f2", arguments=["rho"])
+    field = fea.outputs_dict["stress"]
+    assert field["shape"] == 5
+    assert len(field["partials"]) == 1
diff --git a/tests/test_04_fea_add_field_output.py b/tests/test_04_fea_add_field_output.py
new file mode 100644
index 0000000..93f37ad
--- /dev/null
+++ b/tests/test_04_fea_add_field_output.py
@@ -0,0 +1,21 @@
+import numpy as np
+import types
+
+
+def test_add_field_output_sets_shape_and_metadata(fresh_import):
+    fea_mod = fresh_import("femo.fea.fea_dolfinx")
+    fea = fea_mod.FEA(mesh=object())
+
+    class DummyFunc:
+        def __init__(self):
+            self.vector = types.SimpleNamespace(getArray=lambda: np.zeros(6))
+
+    fea_mod.FunctionSpace = lambda mesh, desc: object()
+    fea_mod.Function = lambda V: DummyFunc()
+    fea_mod.getFuncArray = lambda f: np.zeros(6)
+
+    fea.add_field_output("disp", form="expr", arguments=["u"], record=True)
+    out = fea.outputs_field_dict["disp"]
+    assert out["shape"] == 6
+    assert out["arguments"] == ["u"]
+    assert out["record"] is True
diff --git a/tests/test_05_fea_boundary_conditions.py b/tests/test_05_fea_boundary_conditions.py
new file mode 100644
index 0000000..fd608da
--- /dev/null
+++ b/tests/test_05_fea_boundary_conditions.py
@@ -0,0 +1,12 @@
+def test_add_strong_bc_appends_with_and_without_space(fresh_import):
+    fea_mod = fresh_import("femo.fea.fea_dolfinx")
+    fea = fea_mod.FEA(mesh=object())
+
+    fea_mod.dirichletbc = lambda *args: args
+
+    fea.add_strong_bc(ubc="u0", locate_BC_list=[1, 2])
+    fea.add_strong_bc(ubc="u1", locate_BC_list=[3], function_space="V")
+
+    assert fea.bc[0] == ("u0", 1)
+    assert fea.bc[1] == ("u0", 2)
+    assert fea.bc[2] == ("u1", 3, "V")
diff --git a/tests/test_06_utils_find_node_indices.py b/tests/test_06_utils_find_node_indices.py
new file mode 100644
index 0000000..4ed8118
--- /dev/null
+++ b/tests/test_06_utils_find_node_indices.py
@@ -0,0 +1,27 @@
+import numpy as np
+
+
+def test_find_node_indices_returns_nearest(fresh_import):
+    utils = fresh_import("femo.fea.utils_dolfinx")
+
+    class SimpleKDTree:
+        def __init__(self, points):
+            self.points = np.asarray(points)
+
+        def query(self, nodes):
+            nodes = np.asarray(nodes)
+            idx = []
+            dist = []
+            for node in nodes:
+                d = np.sqrt(((self.points - node) ** 2).sum(axis=1))
+                i = int(np.argmin(d))
+                idx.append(i)
+                dist.append(d[i])
+            return np.array(dist), np.array(idx)
+
+    utils.KDTree = SimpleKDTree
+
+    coords = np.array([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0]])
+    nodes = np.array([[1.1, 0.0], [0.2, 0.0]])
+    got = utils.findNodeIndices(nodes, coords)
+    assert got.tolist() == [1, 0]
diff --git a/tests/test_07_utils_locate_dofs.py b/tests/test_07_utils_locate_dofs.py
new file mode 100644
index 0000000..7034f82
--- /dev/null
+++ b/tests/test_07_utils_locate_dofs.py
@@ -0,0 +1,43 @@
+import numpy as np
+
+
+def test_locate_dofs_converts_polar_input(fresh_import):
+    utils = fresh_import("femo.fea.utils_dolfinx")
+
+    captured = {}
+
+    def fake_find(coords, all_coords):
+        captured["coords"] = coords.copy()
+        return np.array([0, 1])
+
+    utils.findNodeIndices = fake_find
+
+    class DummyV:
+        def tabulate_dof_coordinates(self):
+            return np.array([[1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
+
+    polar = np.array([0.0, 1.0, np.pi / 2, 1.0])
+    out = utils.locateDOFs(polar, DummyV(), input="polar")
+    assert np.allclose(captured["coords"], np.array([[1.0, 0.0], [0.0, 1.0]]), atol=1e-12)
+    assert out.tolist() == [0, 1, 2, 3]
+
+
+def test_locate_dofs_keeps_cartesian_input(fresh_import):
+    utils = fresh_import("femo.fea.utils_dolfinx")
+
+    captured = {}
+
+    def fake_find(coords, all_coords):
+        captured["coords"] = coords.copy()
+        return np.array([1])
+
+    utils.findNodeIndices = fake_find
+
+    class DummyV:
+        def tabulate_dof_coordinates(self):
+            return np.array([[0.0, 0.0, 0.0], [2.0, 3.0, 0.0]])
+
+    cart = np.array([2.0, 3.0])
+    out = utils.locateDOFs(cart, DummyV(), input="cartesian")
+    assert np.allclose(captured["coords"], np.array([[2.0, 3.0]]))
+    assert out.tolist() == [2, 3]
diff --git a/tests/test_08_utils_import_mesh.py b/tests/test_08_utils_import_mesh.py
new file mode 100644
index 0000000..df4dc5c
--- /dev/null
+++ b/tests/test_08_utils_import_mesh.py
@@ -0,0 +1,35 @@
+import types
+
+
+def test_import_mesh_parses_association_table_and_returns_expected_tuple(tmp_path, fresh_import):
+    utils = fresh_import("femo.fea.utils_dolfinx")
+
+    assoc = tmp_path / "mesh_association_table.ini"
+    assoc.write_text("[ASSOCIATION TABLE]\nleft=1\nright=2\n", encoding="utf-8")
+
+    class FakeXDMF:
+        def __init__(self, comm, path, mode):
+            self.path = path
+
+        def __enter__(self):
+            return self
+
+        def __exit__(self, exc_type, exc, tb):
+            return False
+
+        def read_mesh(self, name="Grid"):
+            return types.SimpleNamespace(topology=types.SimpleNamespace(dim=2, create_connectivity=lambda *a, **k: None))
+
+        def read_meshtags(self, mesh, name="Grid"):
+            return "boundaries" if "boundaries" in self.path else "subdomains"
+
+    utils.XDMFFile = FakeXDMF
+
+    mesh, boundaries, association = utils.import_mesh(prefix="mesh", directory=str(tmp_path), subdomains=False)
+    assert boundaries == "boundaries"
+    assert association == {"left": 1, "right": 2}
+
+    mesh, boundaries, subdomains, association = utils.import_mesh(prefix="mesh", directory=str(tmp_path), subdomains=True)
+    assert boundaries == "boundaries"
+    assert subdomains == "subdomains"
+    assert association["left"] == 1
diff --git a/tests/test_09_csdl_fea_model.py b/tests/test_09_csdl_fea_model.py
new file mode 100644
index 0000000..7d3616a
--- /dev/null
+++ b/tests/test_09_csdl_fea_model.py
@@ -0,0 +1,25 @@
+import types
+
+
+def test_fea_model_define_adds_state_and_output_models(fresh_import):
+    fea_model_mod = fresh_import("femo.csdl_opt.fea_model")
+
+    made = []
+    fea_model_mod.StateModel = lambda **kwargs: ("state", kwargs)
+    fea_model_mod.OutputModel = lambda **kwargs: ("output", kwargs)
+    fea_model_mod.OutputFieldModel = lambda **kwargs: ("field", kwargs)
+
+    fake_fea = types.SimpleNamespace(
+        states_dict={"u": {"arguments": ["rho"]}},
+        outputs_dict={"mass": {"arguments": ["rho", "u"]}},
+        outputs_field_dict={"disp": {"arguments": ["u"]}},
+    )
+
+    model = fea_model_mod.FEAModel()
+    model.parameters["fea"] = [fake_fea]
+    model.define()
+
+    names = [name for name, _ in model._added]
+    assert "u_state_model" in names
+    assert "mass_output_model" in names
+    assert "disp_output_model" in names
diff --git a/tests/test_10_output_operation.py b/tests/test_10_output_operation.py
new file mode 100644
index 0000000..6c0dc2a
--- /dev/null
+++ b/tests/test_10_output_operation.py
@@ -0,0 +1,29 @@
+import numpy as np
+import types
+
+
+def test_output_operation_compute_and_derivatives(fresh_import):
+    mod = fresh_import("femo.csdl_opt.output_model")
+
+    update_calls = []
+    mod.update = lambda fn, arr: update_calls.append((fn, tuple(arr.tolist())))
+    mod.assemble = lambda form, dim=0: 7.0 if dim == 0 else np.array([7.0, 8.0])
+    mod.computePartials = lambda form, fn: ("partial", form, fn)
+
+    fn = types.SimpleNamespace()
+    fea = types.SimpleNamespace(outputs_dict={
+        "mass": {"form": "m_form", "shape": 1}
+    })
+    args_dict = {"rho": {"shape": 2, "function": fn}}
+
+    op = mod.OutputOperation(fea=fea, args_dict=args_dict, output_name="mass")
+    op.define()
+
+    outputs = {"mass": np.zeros(1)}
+    op.compute(inputs={"rho": np.array([1.0, 2.0])}, outputs=outputs)
+    assert outputs["mass"].tolist() == 7.0
+    assert len(update_calls) == 1
+
+    derivatives = {}
+    op.compute_derivatives(inputs={"rho": np.array([2.0, 3.0])}, derivatives=derivatives)
+    assert ("mass", "rho") in derivatives
diff --git a/tests/test_11_state_operation.py b/tests/test_11_state_operation.py
new file mode 100644
index 0000000..ab8fcf0
--- /dev/null
+++ b/tests/test_11_state_operation.py
@@ -0,0 +1,72 @@
+import numpy as np
+import types
+
+
+def test_state_operation_core_paths(fresh_import):
+    mod = fresh_import("femo.csdl_opt.state_model")
+
+    update_calls = []
+    mod.update = lambda fn, arr: update_calls.append((fn, tuple(arr.tolist())))
+    mod.assembleVector = lambda form: np.array([5.0])
+    mod.getFuncArray = lambda fn: np.array([9.0])
+
+    class FakeFEA:
+        def __init__(self):
+            self.states_dict = {
+                "u": {
+                    "function": object(),
+                    "residual_form": "res",
+                    "shape": 1,
+                    "d_residual": object(),
+                    "d_state": object(),
+                    "record": False,
+                    "arguments": ["rho"],
+                }
+            }
+            self.inputs_dict = {}
+            self.bc = []
+            self.linear_problem = False
+            self.record = False
+            self.opt_iter = 0
+            self.solve_calls = []
+
+        def solve(self, res, fn, bcs):
+            self.solve_calls.append((res, fn, bcs))
+
+        def solveLinearFwd(self, du, A, dR, rhs, ksp):
+            return np.array([11.0])
+
+        def solveLinearBwd(self, dR, A, du, rhs, ksp):
+            return np.array([12.0])
+
+    fea = FakeFEA()
+    arg_fn = object()
+    args_dict = {"rho": {"shape": 1, "function": arg_fn, "record": False}}
+
+    op = mod.StateOperation(debug_mode=False, fea=fea, args_dict=args_dict, state_name="u")
+    op.define()
+
+    residuals = {"u": np.zeros(1)}
+    op.evaluate_residuals(inputs={"rho": np.array([1.0])}, outputs={"u": np.array([2.0])}, residuals=residuals)
+    assert residuals["u"].tolist() == [5.0]
+
+    outputs = {"u": np.array([0.0])}
+    op.solve_residual_equations(inputs={"rho": np.array([3.0])}, outputs=outputs)
+    assert fea.opt_iter == 1
+    assert outputs["u"].tolist() == [9.0]
+    assert len(fea.solve_calls) == 1
+
+    op.du = object()
+    op.dR = object()
+    op.A = object()
+    op.ksp = None
+
+    d_outputs = {"u": np.array([0.0])}
+    d_residuals = {"u": np.array([1.0])}
+    op.apply_inverse_jacobian(d_outputs=d_outputs, d_residuals=d_residuals, mode="fwd")
+    assert d_outputs["u"].tolist() == [11.0]
+
+    d_outputs = {"u": np.array([2.0])}
+    d_residuals = {"u": np.array([0.0])}
+    op.apply_inverse_jacobian(d_outputs=d_outputs, d_residuals=d_residuals, mode="rev")
+    assert d_residuals["u"].tolist() == [12.0]
diff --git a/tests/test_12_smoke_imports.py b/tests/test_12_smoke_imports.py
new file mode 100644
index 0000000..8c81cae
--- /dev/null
+++ b/tests/test_12_smoke_imports.py
@@ -0,0 +1,6 @@
+def test_smoke_imports_core_modules(fresh_import):
+    assert fresh_import("femo.fea.fea_dolfinx") is not None
+    assert fresh_import("femo.fea.utils_dolfinx") is not None
+    assert fresh_import("femo.csdl_opt.fea_model") is not None
+    assert fresh_import("femo.csdl_opt.output_model") is not None
+    assert fresh_import("femo.csdl_opt.state_model") is not None

From f4da30471e49b5d41257b801e58d80480f16ae27 Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 18 Jun 2026 08:51:12 +0000
Subject: [PATCH 3/4] test: add docstrings across unit test files

---
 tests/test_01_fea_add_input.py           | 3 +++
 tests/test_02_fea_add_state.py           | 3 +++
 tests/test_03_fea_add_output.py          | 3 +++
 tests/test_04_fea_add_field_output.py    | 3 +++
 tests/test_05_fea_boundary_conditions.py | 3 +++
 tests/test_06_utils_find_node_indices.py | 3 +++
 tests/test_07_utils_locate_dofs.py       | 4 ++++
 tests/test_08_utils_import_mesh.py       | 3 +++
 tests/test_09_csdl_fea_model.py          | 3 +++
 tests/test_10_output_operation.py        | 3 +++
 tests/test_11_state_operation.py         | 3 +++
 tests/test_12_smoke_imports.py           | 3 +++
 12 files changed, 37 insertions(+)

diff --git a/tests/test_01_fea_add_input.py b/tests/test_01_fea_add_input.py
index 8e99024..d6291ae 100644
--- a/tests/test_01_fea_add_input.py
+++ b/tests/test_01_fea_add_input.py
@@ -1,9 +1,12 @@
+"""Unit tests for FEA input registration behavior."""
+
 import types
 import numpy as np
 import pytest
 
 
 def test_add_input_rejects_duplicate_name(fresh_import):
+    """Verify add_input stores metadata and rejects duplicate input names."""
     fea_mod = fresh_import("femo.fea.fea_dolfinx")
     fea = fea_mod.FEA(mesh=object())
 
diff --git a/tests/test_02_fea_add_state.py b/tests/test_02_fea_add_state.py
index e4b9464..5e56b6f 100644
--- a/tests/test_02_fea_add_state.py
+++ b/tests/test_02_fea_add_state.py
@@ -1,8 +1,11 @@
+"""Unit tests for FEA state registration."""
+
 import types
 import numpy as np
 
 
 def test_add_state_registers_expected_keys(fresh_import):
+    """Verify add_state records the expected dictionary fields."""
     fea_mod = fresh_import("femo.fea.fea_dolfinx")
     fea = fea_mod.FEA(mesh=object())
 
diff --git a/tests/test_03_fea_add_output.py b/tests/test_03_fea_add_output.py
index 6753cd4..1abb867 100644
--- a/tests/test_03_fea_add_output.py
+++ b/tests/test_03_fea_add_output.py
@@ -1,8 +1,11 @@
+"""Unit tests for scalar and field output registration."""
+
 import types
 import numpy as np
 
 
 def test_add_output_handles_scalar_and_field_shapes(fresh_import):
+    """Verify add_output handles both scalar and field output contracts."""
     fea_mod = fresh_import("femo.fea.fea_dolfinx")
     fea = fea_mod.FEA(mesh=object())
 
diff --git a/tests/test_04_fea_add_field_output.py b/tests/test_04_fea_add_field_output.py
index 93f37ad..7b4b91f 100644
--- a/tests/test_04_fea_add_field_output.py
+++ b/tests/test_04_fea_add_field_output.py
@@ -1,8 +1,11 @@
+"""Unit tests for field output metadata registration."""
+
 import numpy as np
 import types
 
 
 def test_add_field_output_sets_shape_and_metadata(fresh_import):
+    """Verify add_field_output stores shape, arguments, and record flag."""
     fea_mod = fresh_import("femo.fea.fea_dolfinx")
     fea = fea_mod.FEA(mesh=object())
 
diff --git a/tests/test_05_fea_boundary_conditions.py b/tests/test_05_fea_boundary_conditions.py
index fd608da..f9aeb7e 100644
--- a/tests/test_05_fea_boundary_conditions.py
+++ b/tests/test_05_fea_boundary_conditions.py
@@ -1,4 +1,7 @@
+"""Unit tests for strong boundary condition registration."""
+
 def test_add_strong_bc_appends_with_and_without_space(fresh_import):
+    """Verify add_strong_bc appends BC entries for both call signatures."""
     fea_mod = fresh_import("femo.fea.fea_dolfinx")
     fea = fea_mod.FEA(mesh=object())
 
diff --git a/tests/test_06_utils_find_node_indices.py b/tests/test_06_utils_find_node_indices.py
index 4ed8118..d998d00 100644
--- a/tests/test_06_utils_find_node_indices.py
+++ b/tests/test_06_utils_find_node_indices.py
@@ -1,7 +1,10 @@
+"""Unit tests for node index lookup utility."""
+
 import numpy as np
 
 
 def test_find_node_indices_returns_nearest(fresh_import):
+    """Verify nearest-node indices are returned for requested coordinates."""
     utils = fresh_import("femo.fea.utils_dolfinx")
 
     class SimpleKDTree:
diff --git a/tests/test_07_utils_locate_dofs.py b/tests/test_07_utils_locate_dofs.py
index 7034f82..2c50996 100644
--- a/tests/test_07_utils_locate_dofs.py
+++ b/tests/test_07_utils_locate_dofs.py
@@ -1,7 +1,10 @@
+"""Unit tests for locating DOFs from point coordinates."""
+
 import numpy as np
 
 
 def test_locate_dofs_converts_polar_input(fresh_import):
+    """Verify polar coordinate input is converted before DOF lookup."""
     utils = fresh_import("femo.fea.utils_dolfinx")
 
     captured = {}
@@ -23,6 +26,7 @@ def tabulate_dof_coordinates(self):
 
 
 def test_locate_dofs_keeps_cartesian_input(fresh_import):
+    """Verify cartesian coordinate input is used directly for DOF lookup."""
     utils = fresh_import("femo.fea.utils_dolfinx")
 
     captured = {}
diff --git a/tests/test_08_utils_import_mesh.py b/tests/test_08_utils_import_mesh.py
index df4dc5c..7173074 100644
--- a/tests/test_08_utils_import_mesh.py
+++ b/tests/test_08_utils_import_mesh.py
@@ -1,7 +1,10 @@
+"""Unit tests for mesh import helper contracts."""
+
 import types
 
 
 def test_import_mesh_parses_association_table_and_returns_expected_tuple(tmp_path, fresh_import):
+    """Verify import_mesh parses associations and return arity for both modes."""
     utils = fresh_import("femo.fea.utils_dolfinx")
 
     assoc = tmp_path / "mesh_association_table.ini"
diff --git a/tests/test_09_csdl_fea_model.py b/tests/test_09_csdl_fea_model.py
index 7d3616a..c8fbe0e 100644
--- a/tests/test_09_csdl_fea_model.py
+++ b/tests/test_09_csdl_fea_model.py
@@ -1,7 +1,10 @@
+"""Unit tests for FEAModel CSDL wiring."""
+
 import types
 
 
 def test_fea_model_define_adds_state_and_output_models(fresh_import):
+    """Verify FEAModel.define adds state, output, and field output submodels."""
     fea_model_mod = fresh_import("femo.csdl_opt.fea_model")
 
     made = []
diff --git a/tests/test_10_output_operation.py b/tests/test_10_output_operation.py
index 6c0dc2a..15b2c52 100644
--- a/tests/test_10_output_operation.py
+++ b/tests/test_10_output_operation.py
@@ -1,8 +1,11 @@
+"""Unit tests for OutputOperation compute paths."""
+
 import numpy as np
 import types
 
 
 def test_output_operation_compute_and_derivatives(fresh_import):
+    """Verify output compute and derivative paths populate expected values."""
     mod = fresh_import("femo.csdl_opt.output_model")
 
     update_calls = []
diff --git a/tests/test_11_state_operation.py b/tests/test_11_state_operation.py
index ab8fcf0..ac77ebc 100644
--- a/tests/test_11_state_operation.py
+++ b/tests/test_11_state_operation.py
@@ -1,8 +1,11 @@
+"""Unit tests for StateOperation control-flow paths."""
+
 import numpy as np
 import types
 
 
 def test_state_operation_core_paths(fresh_import):
+    """Verify residual, solve, and inverse-jacobian paths in StateOperation."""
     mod = fresh_import("femo.csdl_opt.state_model")
 
     update_calls = []
diff --git a/tests/test_12_smoke_imports.py b/tests/test_12_smoke_imports.py
index 8c81cae..d913b7a 100644
--- a/tests/test_12_smoke_imports.py
+++ b/tests/test_12_smoke_imports.py
@@ -1,4 +1,7 @@
+"""Unit smoke tests for importing core FEMO modules."""
+
 def test_smoke_imports_core_modules(fresh_import):
+    """Verify core package modules import successfully."""
     assert fresh_import("femo.fea.fea_dolfinx") is not None
     assert fresh_import("femo.fea.utils_dolfinx") is not None
     assert fresh_import("femo.csdl_opt.fea_model") is not None

From 54fdb858115f1441d2efb82fdcd635ebbbad5abd Mon Sep 17 00:00:00 2001
From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com>
Date: Thu, 18 Jun 2026 08:53:30 +0000
Subject: [PATCH 4/4] test: move unit tests into tests/unit

---
 tests/{ => unit}/test_01_fea_add_input.py           | 0
 tests/{ => unit}/test_02_fea_add_state.py           | 0
 tests/{ => unit}/test_03_fea_add_output.py          | 0
 tests/{ => unit}/test_04_fea_add_field_output.py    | 0
 tests/{ => unit}/test_05_fea_boundary_conditions.py | 0
 tests/{ => unit}/test_06_utils_find_node_indices.py | 0
 tests/{ => unit}/test_07_utils_locate_dofs.py       | 0
 tests/{ => unit}/test_08_utils_import_mesh.py       | 0
 tests/{ => unit}/test_09_csdl_fea_model.py          | 0
 tests/{ => unit}/test_10_output_operation.py        | 0
 tests/{ => unit}/test_11_state_operation.py         | 0
 tests/{ => unit}/test_12_smoke_imports.py           | 0
 12 files changed, 0 insertions(+), 0 deletions(-)
 rename tests/{ => unit}/test_01_fea_add_input.py (100%)
 rename tests/{ => unit}/test_02_fea_add_state.py (100%)
 rename tests/{ => unit}/test_03_fea_add_output.py (100%)
 rename tests/{ => unit}/test_04_fea_add_field_output.py (100%)
 rename tests/{ => unit}/test_05_fea_boundary_conditions.py (100%)
 rename tests/{ => unit}/test_06_utils_find_node_indices.py (100%)
 rename tests/{ => unit}/test_07_utils_locate_dofs.py (100%)
 rename tests/{ => unit}/test_08_utils_import_mesh.py (100%)
 rename tests/{ => unit}/test_09_csdl_fea_model.py (100%)
 rename tests/{ => unit}/test_10_output_operation.py (100%)
 rename tests/{ => unit}/test_11_state_operation.py (100%)
 rename tests/{ => unit}/test_12_smoke_imports.py (100%)

diff --git a/tests/test_01_fea_add_input.py b/tests/unit/test_01_fea_add_input.py
similarity index 100%
rename from tests/test_01_fea_add_input.py
rename to tests/unit/test_01_fea_add_input.py
diff --git a/tests/test_02_fea_add_state.py b/tests/unit/test_02_fea_add_state.py
similarity index 100%
rename from tests/test_02_fea_add_state.py
rename to tests/unit/test_02_fea_add_state.py
diff --git a/tests/test_03_fea_add_output.py b/tests/unit/test_03_fea_add_output.py
similarity index 100%
rename from tests/test_03_fea_add_output.py
rename to tests/unit/test_03_fea_add_output.py
diff --git a/tests/test_04_fea_add_field_output.py b/tests/unit/test_04_fea_add_field_output.py
similarity index 100%
rename from tests/test_04_fea_add_field_output.py
rename to tests/unit/test_04_fea_add_field_output.py
diff --git a/tests/test_05_fea_boundary_conditions.py b/tests/unit/test_05_fea_boundary_conditions.py
similarity index 100%
rename from tests/test_05_fea_boundary_conditions.py
rename to tests/unit/test_05_fea_boundary_conditions.py
diff --git a/tests/test_06_utils_find_node_indices.py b/tests/unit/test_06_utils_find_node_indices.py
similarity index 100%
rename from tests/test_06_utils_find_node_indices.py
rename to tests/unit/test_06_utils_find_node_indices.py
diff --git a/tests/test_07_utils_locate_dofs.py b/tests/unit/test_07_utils_locate_dofs.py
similarity index 100%
rename from tests/test_07_utils_locate_dofs.py
rename to tests/unit/test_07_utils_locate_dofs.py
diff --git a/tests/test_08_utils_import_mesh.py b/tests/unit/test_08_utils_import_mesh.py
similarity index 100%
rename from tests/test_08_utils_import_mesh.py
rename to tests/unit/test_08_utils_import_mesh.py
diff --git a/tests/test_09_csdl_fea_model.py b/tests/unit/test_09_csdl_fea_model.py
similarity index 100%
rename from tests/test_09_csdl_fea_model.py
rename to tests/unit/test_09_csdl_fea_model.py
diff --git a/tests/test_10_output_operation.py b/tests/unit/test_10_output_operation.py
similarity index 100%
rename from tests/test_10_output_operation.py
rename to tests/unit/test_10_output_operation.py
diff --git a/tests/test_11_state_operation.py b/tests/unit/test_11_state_operation.py
similarity index 100%
rename from tests/test_11_state_operation.py
rename to tests/unit/test_11_state_operation.py
diff --git a/tests/test_12_smoke_imports.py b/tests/unit/test_12_smoke_imports.py
similarity index 100%
rename from tests/test_12_smoke_imports.py
rename to tests/unit/test_12_smoke_imports.py