diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py
index 829f43e..1282ea6 100644
--- a/meanas/fdfd/operators.py
+++ b/meanas/fdfd/operators.py
@@ -236,10 +236,12 @@ def eh_full(
     else:
         pm = sparse.diags_array(numpy.where(pmc, 0, 1))    # set pm to (not PMC)
 
-    iwe = pe @ (1j * omega * sparse.diags_array(epsilon)) @ pe
-    iwm = 1j * omega
-    if mu is not None:
-        iwm *= sparse.diags_array(mu)
+    iwe = pe @ (1j * omega * sparse.diags(epsilon)) @ pe
+    if mu is None:
+        iwm = 1j * omega * sparse.eye(epsilon.size)
+    else:
+        iwm = 1j * omega * sparse.diags(mu)
+
     iwm = pm @ iwm @ pm
 
     A1 = pe @ curl_back(dxes[1]) @ pm
diff --git a/meanas/test/test_fdfd_algebra_helpers.py b/meanas/test/test_fdfd_algebra_helpers.py
index c1d1b3f..b481023 100644
--- a/meanas/test/test_fdfd_algebra_helpers.py
+++ b/meanas/test/test_fdfd_algebra_helpers.py
@@ -49,6 +49,21 @@ def _apply_matrix(op: operators.sparse.spmatrix, field: numpy.ndarray, shape: tu
     return unvec(op @ vec(field), shape)
 
 
+def _dense_e_full(mu: numpy.ndarray | None) -> numpy.ndarray:
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+    pe = numpy.where(PEC, 0.0, 1.0)
+    pm = numpy.where(PMC, 0.0, 1.0)
+
+    masked_e = pe * E_FIELD
+    curl_term = ce(masked_e)
+    if mu is not None:
+        curl_term = curl_term / mu
+    curl_term = pm * curl_term
+    curl_term = ch(curl_term)
+    return pe * (curl_term - OMEGA**2 * EPSILON * masked_e)
+
+
 def _dense_h_full(mu: numpy.ndarray | None) -> numpy.ndarray:
     ce = fd_functional.curl_forward(DXES[0])
     ch = fd_functional.curl_back(DXES[1])
@@ -78,6 +93,55 @@ def test_h_full_matches_dense_reference_with_and_without_mu() -> None:
         assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
 
 
+def test_e_full_matches_dense_reference_with_masks() -> None:
+    for mu in (None, MU):
+        matrix_result = _apply_matrix(
+            operators.e_full(OMEGA, DXES, vec(EPSILON), None if mu is None else vec(mu), vec(PEC), vec(PMC)),
+            E_FIELD,
+            SHAPE,
+        )
+        dense_result = _dense_e_full(mu)
+        assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
+
+
+def test_h_full_without_masks_matches_dense_reference() -> None:
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+    dense_result = ce(ch(H_FIELD) / EPSILON) - OMEGA**2 * MU * H_FIELD
+    matrix_result = _apply_matrix(
+        operators.h_full(OMEGA, DXES, vec(EPSILON), vec(MU)),
+        H_FIELD,
+        SHAPE,
+    )
+    assert_allclose(matrix_result, dense_result, atol=1e-10, rtol=1e-10)
+
+
+def test_eh_full_matches_manual_block_operator_with_masks() -> None:
+    pe = numpy.where(PEC, 0.0, 1.0)
+    pm = numpy.where(PMC, 0.0, 1.0)
+    ce = fd_functional.curl_forward(DXES[0])
+    ch = fd_functional.curl_back(DXES[1])
+
+    matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON), vec(MU), vec(PEC), vec(PMC)) @ numpy.concatenate(
+        [vec(E_FIELD), vec(H_FIELD)],
+    )
+    matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
+
+    dense_e = pe * ch(pm * H_FIELD) - pe * (1j * OMEGA * EPSILON * (pe * E_FIELD))
+    dense_h = pm * ce(pe * E_FIELD) + pm * (1j * OMEGA * MU * (pm * H_FIELD))
+
+    assert_allclose(matrix_e, dense_e, atol=1e-10, rtol=1e-10)
+    assert_allclose(matrix_h, dense_h, atol=1e-10, rtol=1e-10)
+
+
+def test_e2h_pmc_mask_matches_masked_unmasked_result() -> None:
+    pmc_complement = numpy.where(PMC, 0.0, 1.0)
+    unmasked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU)), E_FIELD, SHAPE)
+    masked = _apply_matrix(operators.e2h(OMEGA, DXES, vec(MU), vec(PMC)), E_FIELD, SHAPE)
+
+    assert_allclose(masked, pmc_complement * unmasked, atol=1e-10, rtol=1e-10)
+
+
 def test_poynting_h_cross_matches_negative_e_cross_relation() -> None:
     h_cross_e = _apply_matrix(operators.poynting_h_cross(vec(H_FIELD), DXES), E_FIELD, SHAPE)
     e_cross_h = _apply_matrix(operators.poynting_e_cross(vec(E_FIELD), DXES), H_FIELD, SHAPE)
diff --git a/meanas/test/test_fdfd_functional.py b/meanas/test/test_fdfd_functional.py
index f4fd4bb..5f0adef 100644
--- a/meanas/test/test_fdfd_functional.py
+++ b/meanas/test/test_fdfd_functional.py
@@ -70,6 +70,15 @@ def test_eh_full_matches_sparse_operator_with_mu() -> None:
     assert_fields_match(functional_h, matrix_h)
 
 
+def test_eh_full_matches_sparse_operator_without_mu() -> None:
+    matrix_result = operators.eh_full(OMEGA, DXES, vec(EPSILON)) @ numpy.concatenate([vec(E_FIELD), vec(H_FIELD)])
+    matrix_e, matrix_h = (unvec(part, SHAPE) for part in numpy.split(matrix_result, 2))
+    functional_e, functional_h = functional.eh_full(OMEGA, DXES, EPSILON)(E_FIELD, H_FIELD)
+
+    assert_fields_match(functional_e, matrix_e)
+    assert_fields_match(functional_h, matrix_h)
+
+
 def test_e2h_matches_sparse_operator_with_mu() -> None:
     matrix_result = apply_matrix(
         operators.e2h(OMEGA, DXES, vec(MU)),
@@ -80,6 +89,16 @@ def test_e2h_matches_sparse_operator_with_mu() -> None:
     assert_fields_match(functional_result, matrix_result)
 
 
+def test_e2h_matches_sparse_operator_without_mu() -> None:
+    matrix_result = apply_matrix(
+        operators.e2h(OMEGA, DXES),
+        E_FIELD,
+    )
+    functional_result = functional.e2h(OMEGA, DXES)(E_FIELD)
+
+    assert_fields_match(functional_result, matrix_result)
+
+
 def test_m2j_matches_sparse_operator_without_mu() -> None:
     matrix_result = apply_matrix(
         operators.m2j(OMEGA, DXES),
diff --git a/meanas/test/test_fdmath_operators.py b/meanas/test/test_fdmath_operators.py
new file mode 100644
index 0000000..bb7fe31
--- /dev/null
+++ b/meanas/test/test_fdmath_operators.py
@@ -0,0 +1,89 @@
+import numpy
+import pytest
+from numpy.testing import assert_allclose, assert_array_equal
+
+from ..fdmath import operators, unvec, vec
+
+
+SHAPE = (2, 3, 2)
+SCALAR_FIELD = numpy.arange(numpy.prod(SHAPE), dtype=float).reshape(SHAPE, order='C')
+VECTOR_LEFT = (numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') + 0.5)
+VECTOR_RIGHT = (2.0 + numpy.arange(3 * numpy.prod(SHAPE), dtype=float).reshape((3, *SHAPE), order='C') / 3.0)
+
+
+def _apply_scalar_matrix(op: operators.sparse.spmatrix) -> numpy.ndarray:
+    return (op @ SCALAR_FIELD.ravel(order='C')).reshape(SHAPE, order='C')
+
+
+def _mirrored_indices(size: int, shift_distance: int) -> numpy.ndarray:
+    indices = numpy.arange(size) + shift_distance
+    indices = numpy.where(indices >= size, 2 * size - indices - 1, indices)
+    indices = numpy.where(indices < 0, -1 - indices, indices)
+    return indices
+
+
+@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
+def test_shift_circ_matches_numpy_roll(axis: int, shift_distance: int) -> None:
+    matrix_result = _apply_scalar_matrix(operators.shift_circ(axis, SHAPE, shift_distance))
+    expected = numpy.roll(SCALAR_FIELD, -shift_distance, axis=axis)
+    assert_array_equal(matrix_result, expected)
+
+
+@pytest.mark.parametrize(('axis', 'shift_distance'), [(0, 1), (1, -1), (2, 1)])
+def test_shift_with_mirror_matches_explicit_mirrored_indices(axis: int, shift_distance: int) -> None:
+    matrix_result = _apply_scalar_matrix(operators.shift_with_mirror(axis, SHAPE, shift_distance))
+    indices = [numpy.arange(length) for length in SHAPE]
+    indices[axis] = _mirrored_indices(SHAPE[axis], shift_distance)
+    expected = SCALAR_FIELD[numpy.ix_(*indices)]
+    assert_array_equal(matrix_result, expected)
+
+
+@pytest.mark.parametrize(
+    ('args', 'message'),
+    [
+        ((0, (2,), 1), 'Invalid shape'),
+        ((3, SHAPE, 1), 'Invalid direction'),
+    ],
+    )
+def test_shift_circ_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
+    with pytest.raises(Exception, match=message):
+        operators.shift_circ(*args)
+
+
+@pytest.mark.parametrize(
+    ('args', 'message'),
+    [
+        ((0, (2,), 1), 'Invalid shape'),
+        ((3, SHAPE, 1), 'Invalid direction'),
+        ((0, SHAPE, SHAPE[0]), 'too large'),
+    ],
+    )
+def test_shift_with_mirror_rejects_invalid_arguments(args: tuple[int, tuple[int, ...], int], message: str) -> None:
+    with pytest.raises(Exception, match=message):
+        operators.shift_with_mirror(*args)
+
+
+def test_vec_cross_matches_pointwise_cross_product() -> None:
+    matrix_result = unvec(operators.vec_cross(vec(VECTOR_LEFT)) @ vec(VECTOR_RIGHT), SHAPE)
+    expected = numpy.empty_like(VECTOR_LEFT)
+    expected[0] = VECTOR_LEFT[1] * VECTOR_RIGHT[2] - VECTOR_LEFT[2] * VECTOR_RIGHT[1]
+    expected[1] = VECTOR_LEFT[2] * VECTOR_RIGHT[0] - VECTOR_LEFT[0] * VECTOR_RIGHT[2]
+    expected[2] = VECTOR_LEFT[0] * VECTOR_RIGHT[1] - VECTOR_LEFT[1] * VECTOR_RIGHT[0]
+    assert_allclose(matrix_result, expected)
+
+
+def test_avg_forward_matches_half_sum_with_forward_neighbor() -> None:
+    matrix_result = _apply_scalar_matrix(operators.avg_forward(1, SHAPE))
+    expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, -1, axis=1))
+    assert_allclose(matrix_result, expected)
+
+
+def test_avg_back_matches_half_sum_with_backward_neighbor() -> None:
+    matrix_result = _apply_scalar_matrix(operators.avg_back(1, SHAPE))
+    expected = 0.5 * (SCALAR_FIELD + numpy.roll(SCALAR_FIELD, 1, axis=1))
+    assert_allclose(matrix_result, expected)
+
+
+def test_avg_forward_rejects_invalid_shape() -> None:
+    with pytest.raises(Exception, match='Invalid shape'):
+        operators.avg_forward(0, (2,))