import numpy as np
import pytest

from cymsim.mesh import disk_mesh


@pytest.mark.parametrize("h", [0.03, 0.012])
def test_disk_mesh_is_a_valid_manifold(h):
    R = 0.2286
    m = disk_mesh(R, h)
    assert (m.area > 0).all()
    # polygon area approaches pi R^2 from below
    assert 0.98 * np.pi * R**2 < m.area.sum() <= np.pi * R**2 + 1e-12
    edges = np.sort(np.concatenate([m.tri[:, [0, 1]], m.tri[:, [1, 2]], m.tri[:, [2, 0]]]), axis=1)
    _, count = np.unique(edges, axis=0, return_counts=True)
    assert set(count) == {1, 2}                      # interior edges shared by two triangles
    assert (count == 1).sum() == len(m.boundary)    # boundary edges = rim nodes
    assert np.allclose(np.hypot(*m.xy[m.boundary].T), R)
    # element quality: inradius / circumradius, 1 for equilateral
    L = m.edge_lengths()
    a, b, c = L.T
    s = (a + b + c) / 2
    q = 8 * (s - a) * (s - b) * (s - c) / (a * b * c)
    assert q.min() > 0.5
    assert 0.5 * h < L.mean() < 1.5 * h


def test_cst_gradient_is_exact_for_linear_fields():
    m = disk_mesh(0.2, 0.03)
    f = 2.0 * m.xy[:, 0] - 3.0 * m.xy[:, 1] + 1.0
    fe = f[m.tri]
    assert np.allclose((m.Gx * fe).sum(1), 2.0)
    assert np.allclose((m.Gy * fe).sum(1), -3.0)


def test_locate_and_interpolate():
    m = disk_mesh(0.2, 0.03)
    f = m.xy[:, 0] ** 1 + 0.5 * m.xy[:, 1]
    nodes, w = m.locate(0.05, -0.07)
    assert np.isclose(w.sum(), 1.0) and (w >= -1e-9).all()
    assert np.isclose(m.interpolate(f, 0.05, -0.07), 0.05 - 0.035)
    with pytest.raises(ValueError):
        m.locate(0.3, 0.0)