diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
new file mode 100644
index 0000000..16296b6
--- /dev/null
+++ b/.github/workflows/tests.yml
@@ -0,0 +1,33 @@
+name: Tests
+
+on:
+  push:
+    branches: ["main", "master"]
+  pull_request:
+
+jobs:
+  unit:
+    name: Unit tests (no FEM deps)
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: "3.11"
+      - run: pip install pytest numpy scipy
+      - run: pip install -e . --no-deps
+      - run: pytest tests/unit/ -v --tb=short
+
+  integration:
+    name: Integration tests (real dolfinx)
+    runs-on: ubuntu-latest
+    container:
+      image: dolfinx/dolfinx:v0.5.1
+    steps:
+      - uses: actions/checkout@v4
+      - name: Install dependencies
+        run: python3 -m pip install pytest scipy matplotlib
+      - name: Install femo
+        run: pip install -e . --no-deps
+      - name: Run integration tests
+        run: pytest tests/integration/ -v --tb=short
diff --git a/tests/conftest.py b/tests/conftest.py
new file mode 100644
index 0000000..58be6f7
--- /dev/null
+++ b/tests/conftest.py
@@ -0,0 +1,52 @@
+"""
+Shared pytest fixtures for femo tests.
+
+Unit tests (tests/unit/) use lightweight mock fixtures that avoid importing
+dolfinx/petsc4py so they can run without a FEM installation.
+
+Integration tests (tests/integration/) import dolfinx directly — they are
+skipped automatically when dolfinx is not installed.
+"""
+import pytest
+import numpy as np
+
+
+# ---------------------------------------------------------------------------
+# Minimal stubs used only by unit tests
+# ---------------------------------------------------------------------------
+
+class _MockFunctionSpace:
+    def __init__(self, n=10):
+        self.n = n
+
+
+class _MockFunction:
+    def __init__(self, name="u", n=10):
+        self.name = name
+        self.n = n
+        self.function_space = _MockFunctionSpace(n)
+
+        class _X:
+            array = np.zeros(n)
+        self.x = _X()
+
+    def rename(self, label, _):
+        self.name = label
+
+
+@pytest.fixture
+def mock_mesh():
+    """A minimal mesh stand-in for unit tests."""
+    class _Mesh:
+        pass
+    return _Mesh()
+
+
+@pytest.fixture
+def mock_function():
+    return _MockFunction()
+
+
+@pytest.fixture
+def mock_function_space():
+    return _MockFunctionSpace()
diff --git a/tests/integration/__init__.py b/tests/integration/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/integration/test_fea_dolfinx.py b/tests/integration/test_fea_dolfinx.py
new file mode 100644
index 0000000..de0aa02
--- /dev/null
+++ b/tests/integration/test_fea_dolfinx.py
@@ -0,0 +1,171 @@
+"""
+Integration tests for FEA using real dolfinx/ufl/petsc4py.
+
+These tests are skipped automatically if dolfinx is not installed.
+Run them inside the official container:
+
+    docker run --rm -v $(pwd):/repo -w /repo dolfinx/dolfinx:v0.9.0 \
+        pytest tests/integration/ -v
+"""
+import pytest
+import numpy as np
+
+dolfinx = pytest.importorskip("dolfinx")
+
+from mpi4py import MPI
+from dolfinx.mesh import create_unit_square, create_unit_cube
+from dolfinx.fem import (FunctionSpace, Function, dirichletbc,
+                          locate_dofs_geometrical, form, assemble_scalar,
+                          Constant)
+from petsc4py.PETSc import ScalarType
+from dolfinx.fem.petsc import assemble_vector, assemble_matrix
+from dolfinx.cpp.mesh import CellType
+import ufl
+from ufl import inner, grad, dx, TestFunction, TrialFunction
+
+from femo.fea.fea_dolfinx import FEA
+
+
+# ---------------------------------------------------------------------------
+# Fixtures
+# ---------------------------------------------------------------------------
+
+@pytest.fixture(scope="module")
+def square_mesh():
+    return create_unit_square(MPI.COMM_WORLD, 8, 8)
+
+
+@pytest.fixture(scope="module")
+def scalar_V(square_mesh):
+    return FunctionSpace(square_mesh, ("Lagrange", 1))
+
+
+# ---------------------------------------------------------------------------
+# FEA registration
+# ---------------------------------------------------------------------------
+
+def test_add_input(square_mesh, scalar_V):
+    fea = FEA(square_mesh)
+    rho = Function(scalar_V)
+    fea.add_input("rho", rho, init_val=1.0)
+
+    assert "rho" in fea.inputs_dict
+    assert np.allclose(fea.inputs_dict["rho"]["function"].x.array, 1.0)
+
+
+def test_add_input_duplicate_raises(square_mesh, scalar_V):
+    fea = FEA(square_mesh)
+    fea.add_input("rho", Function(scalar_V), init_val=0.5)
+    with pytest.raises(ValueError, match="already been used"):
+        fea.add_input("rho", Function(scalar_V), init_val=0.5)
+
+
+def test_add_state(square_mesh, scalar_V):
+    fea = FEA(square_mesh)
+    u = Function(scalar_V)
+    v = TestFunction(scalar_V)
+    res = inner(grad(u), grad(v)) * dx
+    fea.add_state("u", u, res, arguments=[])
+    assert "u" in fea.states_dict
+
+
+def test_add_scalar_output(square_mesh, scalar_V):
+    fea = FEA(square_mesh)
+    u = Function(scalar_V)
+    u.x.array[:] = 1.0
+    fea.add_input("rho", u, init_val=1.0)
+    output_form = u * dx
+    fea.add_output("vol", type="scalar", form=output_form, arguments=["rho"])
+    assert "vol" in fea.outputs_dict
+
+
+def test_add_strong_bc(square_mesh, scalar_V):
+    fea = FEA(square_mesh)
+    ubc = Function(scalar_V)
+    ubc.x.array[:] = 0.0
+
+    def left_boundary(x):
+        return np.isclose(x[0], 0.0)
+
+    dofs = locate_dofs_geometrical(scalar_V, left_boundary)
+    fea.add_strong_bc(ubc, [dofs])
+    assert len(fea.bc) == 1
+
+
+# ---------------------------------------------------------------------------
+# Poisson solve — verifies the full FEA pipeline end-to-end
+# ---------------------------------------------------------------------------
+
+def test_poisson_solve(square_mesh, scalar_V):
+    """
+    Solve -∇²u = f on the unit square with u=0 on the boundary.
+    Checks that the assembled system and Newton solve produce a non-trivial,
+    physically reasonable solution.
+    """
+    fea = FEA(square_mesh)
+    fea.PDE_SOLVER = "Newton"
+    fea.REPORT = False
+
+    u = Function(scalar_V)
+    v = TestFunction(scalar_V)
+    f_val = Constant(square_mesh, ScalarType(1.0))
+
+    # Nonlinear residual: R(u;v) = ∫ ∇u·∇v dx − ∫ f·v dx
+    res = inner(grad(u), grad(v)) * dx - inner(f_val, v) * dx
+
+    def boundary(x):
+        return (np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) |
+                np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0))
+
+    dofs = locate_dofs_geometrical(scalar_V, boundary)
+    ubc = Function(scalar_V)
+    ubc.x.array[:] = 0.0
+    fea.add_strong_bc(ubc, [dofs])
+
+    fea.add_state("u", u, res, arguments=[])
+    fea.solve(res, u, fea.bc)
+
+    u_vals = u.x.array
+    # Solution should be positive and peak near the centre (~0.073 for unit square)
+    assert np.max(u_vals) > 0.05
+    assert np.max(u_vals) < 0.15
+    # Boundary dofs should be (approximately) zero
+    assert np.allclose(u_vals[dofs], 0.0, atol=1e-10)
+
+
+# ---------------------------------------------------------------------------
+# Utility functions from utils_dolfinx
+# ---------------------------------------------------------------------------
+
+def test_getFuncArray(scalar_V):
+    from femo.fea.utils_dolfinx import getFuncArray
+    f = Function(scalar_V)
+    f.x.array[:] = 3.14
+    arr = getFuncArray(f)
+    assert isinstance(arr, np.ndarray)
+    assert np.allclose(arr, 3.14)
+
+
+def test_setFuncArray(scalar_V):
+    from femo.fea.utils_dolfinx import getFuncArray, setFuncArray
+    f = Function(scalar_V)
+    new_vals = np.ones(len(f.x.array)) * 2.71
+    setFuncArray(f, new_vals)
+    assert np.allclose(getFuncArray(f), 2.71)
+
+
+def test_create_unit_square_mesh():
+    """Smoke test: mesh creation and basic topology."""
+    mesh = create_unit_square(MPI.COMM_WORLD, 4, 4)
+    assert mesh.topology.dim == 2
+
+
+def test_assemble_mass_matrix(square_mesh, scalar_V):
+    """Verify that a mass matrix assembles without error and is non-zero."""
+    u = TrialFunction(scalar_V)
+    v = TestFunction(scalar_V)
+    a = form(inner(u, v) * dx)
+    A = assemble_matrix(a)
+    A.assemble()
+    # Frobenius norm should be positive
+    assert A.norm() > 0
diff --git a/tests/integration/test_poisson_convergence.py b/tests/integration/test_poisson_convergence.py
new file mode 100644
index 0000000..ac13205
--- /dev/null
+++ b/tests/integration/test_poisson_convergence.py
@@ -0,0 +1,78 @@
+"""
+Convergence test: verifies the Poisson solver achieves the expected
+O(h²) L2 error rate for a manufactured solution.
+
+This is the canonical correctness check for a finite element code —
+it cannot be replicated with mocks.
+"""
+import pytest
+import numpy as np
+
+dolfinx = pytest.importorskip("dolfinx")
+
+from mpi4py import MPI
+from dolfinx.mesh import create_unit_square
+from dolfinx.fem import (FunctionSpace, Function, dirichletbc,
+                          locate_dofs_geometrical, form, assemble_scalar,
+                          Constant)
+from dolfinx.cpp.mesh import CellType
+import ufl
+from ufl import inner, grad, dx, sin, pi, SpatialCoordinate
+
+from femo.fea.fea_dolfinx import FEA
+
+
+def _solve_poisson(n):
+    """
+    Solve -∇²u = 2π²·sin(πx)·sin(πy) on [0,1]², u=0 on ∂Ω.
+    Exact solution: u* = sin(πx)·sin(πy).
+    Returns (mesh_size h, L2 error).
+    """
+    mesh = create_unit_square(MPI.COMM_WORLD, n, n)
+    V = FunctionSpace(mesh, ("Lagrange", 1))
+    x = SpatialCoordinate(mesh)
+
+    u = Function(V)
+    v = ufl.TestFunction(V)
+
+    f = 2 * pi**2 * ufl.sin(pi * x[0]) * ufl.sin(pi * x[1])
+    res = inner(grad(u), grad(v)) * dx - inner(f, v) * dx
+
+    def boundary(coords):
+        return (np.isclose(coords[0], 0.0) | np.isclose(coords[0], 1.0) |
+                np.isclose(coords[1], 0.0) | np.isclose(coords[1], 1.0))
+
+    dofs = locate_dofs_geometrical(V, boundary)
+    ubc = Function(V)
+    ubc.x.array[:] = 0.0
+
+    fea = FEA(mesh)
+    fea.REPORT = False
+    fea.add_strong_bc(ubc, [dofs])
+    fea.solve(res, u, fea.bc)
+
+    # L2 error against exact solution
+    u_exact = ufl.sin(pi * x[0]) * ufl.sin(pi * x[1])
+    error = form((u - u_exact) ** 2 * dx)
+    L2 = np.sqrt(assemble_scalar(error))
+    h = 1.0 / n
+    return h, L2
+
+
+@pytest.mark.parametrize("n", [4, 8, 16])
+def test_poisson_converges(n):
+    """Each refinement should produce a non-trivial, bounded solution."""
+    h, err = _solve_poisson(n)
+    assert err > 0, "zero error suggests solve did not run"
+    assert err < 0.1, f"error {err:.4e} is too large for n={n}"
+
+
+def test_poisson_convergence_rate():
+    """
+    Check that doubling the mesh halves the error with roughly O(h²) rate.
+    Lagrange P1 on a uniform mesh → rate ≈ 2.
+    """
+    _, e1 = _solve_poisson(8)
+    _, e2 = _solve_poisson(16)
+    rate = np.log(e1 / e2) / np.log(2.0)
+    assert rate > 1.8, f"convergence rate {rate:.2f} is below expected O(h²)"
diff --git a/tests/unit/__init__.py b/tests/unit/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/unit/test_fea_registration.py b/tests/unit/test_fea_registration.py
new file mode 100644
index 0000000..a8058f6
--- /dev/null
+++ b/tests/unit/test_fea_registration.py
@@ -0,0 +1,127 @@
+"""
+Unit tests for FEA input/state/output registration logic in AbstractFEA.
+These tests use mock objects — no dolfinx installation required.
+"""
+import pytest
+import sys
+import types
+
+
+def _setup_stubs():
+    """
+    Populate sys.modules with empty stub modules and set every attribute
+    accessed via 'from X import Y' in fea_dolfinx.py and utils_dolfinx.py.
+    Safe to call multiple times — only creates stubs for modules not yet present.
+    """
+    _noop = lambda *a, **kw: None
+
+    stub_names = [
+        "dolfinx", "dolfinx.io", "dolfinx.fem", "dolfinx.fem.petsc",
+        "dolfinx.mesh", "dolfinx.nls", "dolfinx.nls.petsc",
+        "dolfinx.cpp", "dolfinx.cpp.mesh", "dolfinx.la",
+        "ufl", "petsc4py", "petsc4py.PETSc",
+        "mpi4py", "mpi4py.MPI", "matplotlib", "matplotlib.pyplot",
+        "scipy", "scipy.spatial", "scipy.sparse",
+    ]
+    for name in stub_names:
+        if name not in sys.modules:
+            sys.modules[name] = types.ModuleType(name)
+
+    # Wire parent.child references so `from parent import child` works
+    for name in stub_names:
+        parts = name.split(".")
+        if len(parts) > 1:
+            parent = sys.modules[".".join(parts[:-1])]
+            if not hasattr(parent, parts[-1]):
+                setattr(parent, parts[-1], sys.modules[name])
+
+    # dolfinx.io
+    sys.modules["dolfinx.io"].XDMFFile = _noop
+
+    # ufl  (fea_dolfinx.py and utils_dolfinx.py imports combined)
+    for sym in ["Identity", "dot", "derivative", "TestFunction", "TrialFunction",
+                "inner", "ds", "dS", "dx", "grad", "inv", "as_vector", "sqrt",
+                "conditional", "lt", "det", "Measure", "exp", "tr", "CellDiameter",
+                "SpatialCoordinate", "FacetNormal", "div"]:
+        setattr(sys.modules["ufl"], sym, _noop)
+
+    # dolfinx.fem
+    for sym in ["form", "assemble_scalar", "Function", "FunctionSpace",
+                "dirichletbc", "locate_dofs_geometrical", "locate_dofs_topological",
+                "Constant", "VectorFunctionSpace", "set_bc"]:
+        setattr(sys.modules["dolfinx.fem"], sym, _noop)
+
+    # dolfinx.fem.petsc
+    for sym in ["assemble_vector", "assemble_matrix", "NonlinearProblem",
+                "apply_lifting", "set_bc", "create_matrix", "_assemble_matrix_mat"]:
+        setattr(sys.modules["dolfinx.fem.petsc"], sym, _noop)
+
+    # dolfinx.mesh
+    for sym in ["create_unit_square", "create_rectangle", "create_interval",
+                "locate_entities_boundary", "locate_entities", "meshtags"]:
+        setattr(sys.modules["dolfinx.mesh"], sym, _noop)
+
+    # dolfinx.nls.petsc
+    sys.modules["dolfinx.nls.petsc"].NewtonSolver = _noop
+
+    # dolfinx.cpp.mesh
+    sys.modules["dolfinx.cpp.mesh"].CellType = object()
+
+    # mpi4py.MPI
+    sys.modules["mpi4py.MPI"].COMM_WORLD = object()
+
+    # scipy.spatial / scipy.sparse
+    sys.modules["scipy.spatial"].KDTree = object
+    sys.modules["scipy.sparse"].csr_matrix = _noop
+
+
+def _make_abstract_fea(mock_mesh):
+    """Return an AbstractFEA instance with all heavy deps stubbed out."""
+    _setup_stubs()
+    import importlib
+    import femo.fea.fea_dolfinx as _mod
+    importlib.reload(_mod)
+    return _mod.AbstractFEA(mock_mesh)
+
+
+class MockFunction:
+    def __init__(self, name="f"):
+        self.name = name
+    def rename(self, label, _):
+        self.name = label
+
+
+def test_add_input_stores_entry(mock_mesh):
+    fea = _make_abstract_fea(mock_mesh)
+    f = MockFunction("rho")
+    fea.add_input("rho", f)
+    assert "rho" in fea.inputs_dict
+    assert fea.inputs_dict["rho"]["function"] is f
+
+
+def test_add_input_duplicate_raises(mock_mesh):
+    fea = _make_abstract_fea(mock_mesh)
+    f = MockFunction("rho")
+    fea.add_input("rho", f)
+    with pytest.raises(ValueError, match="already been used"):
+        fea.add_input("rho", MockFunction("rho2"))
+
+
+def test_add_state_stores_entry(mock_mesh):
+    fea = _make_abstract_fea(mock_mesh)
+    u = MockFunction("u")
+    fea.add_state("u", u, residual_form=None)
+    assert "u" in fea.states_dict
+
+
+def test_add_output_stores_entry(mock_mesh):
+    fea = _make_abstract_fea(mock_mesh)
+    fea.add_output("J", form=None)
+    assert "J" in fea.outputs_dict
+
+
+def test_add_strong_bc_appends(mock_mesh):
+    fea = _make_abstract_fea(mock_mesh)
+    bc = object()
+    fea.add_strong_bc(bc)
+    assert bc in fea.bcs_list