def METHOD_NAME(self):
    """
    Test that we reconstruct the exact solution on subfaces
    """
    nx = 2
    ny = 2
    nz = 2
    g = pp.StructuredTetrahedralGrid([nx, ny, nz], physdims=[1, 1, 1])
    g.compute_geometry()
    s_t = pp.fvutils.SubcellTopology(g)
    kxx = 10 * np.ones(g.num_cells)
    k = pp.SecondOrderTensor(kxx)
    bc = pp.BoundaryCondition(g)
    bc.is_dir[g.get_all_boundary_faces()] = True
    bc.is_neu[bc.is_dir] = False
    bc = pp.fvutils.boundary_to_sub_boundary(bc, s_t)
    p_b = np.sum(-g.face_centers, axis=0)
    p_b = p_b[s_t.fno_unique]
    mpfa = pp.Mpfa("flow")
    flux, bound_flux, p_t_cell, p_t_bound, *_ = mpfa._flux_discretization(
        g, k, bc, eta=0, inverter="python"
    )
    div = pp.fvutils.scalar_divergence(g)
    hf2f = pp.fvutils.map_hf_2_f(nd=1, sd=g)
    P = sps.linalg.spsolve(div * hf2f * flux, -div * hf2f * bound_flux * p_b)
    P_hf = p_t_cell * P + p_t_bound * p_b
    _, IA = np.unique(s_t.fno_unique, True)
    P_f = P_hf[IA]
    self.assertTrue(np.all(np.abs(P + np.sum(g.cell_centers, axis=0)) < 1e-10))
    self.assertTrue(np.all(np.abs(P_f + np.sum(g.face_centers, axis=0)) < 1e-10))