Add real dolfinx integration tests and GitHub Actions CI

.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

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()

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

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²)"