grid 3d from polytopal domain

[alessiofumagalli]

Here a minimal example of a cube with a side divided in two parts

    side_1 = np.array([[0, 0, 1, 1], [0, 1, 1, 0], [0, 0, 0, 0]])
    side_2 = np.array([[0, 1, 1, 0], [0, 0, 1, 1], [1, 1, 1, 1]])

    side_3a = np.array([[0, 0, 0, 0], [0, 0, 1, 1], [0, 0.5, 0.5, 0]])
    side_3b = np.array([[0, 0, 0, 0], [0, 0, 1, 1], [0.5, 1, 1, 0.5]])

    side_4 = np.array([[1, 1, 1, 1], [0, 1, 1, 0], [0, 0, 1, 1]])

    side_5 = np.array([[0, 1, 1, 0, 0], [0, 0, 0, 0, 0], [0, 0, 1, 1, 0.5]])
    side_6 = np.array([[0, 1, 1, 0, 0], [1, 1, 1, 1, 1], [0, 0, 1, 1, 0.5]])
    sides = [side_1, side_2, side_3a, side_3b, side_4, side_5, side_6]

    domain = pp.Domain(polytope=sides)
    network = pp.create_fracture_network(domain=domain)

    mesh_size = 0.5
    mesh_kwargs = {"mesh_size_frac": mesh_size, "mesh_size_min": mesh_size / 10}
    mdg = network.mesh(mesh_kwargs)

The creation of the grid is blocked by

porepy/src/porepy/geometry/intersections.py

Line 868 in 12e01f8

assert np.sum(dot_prod_from_main[:-1] == 0) == 2

Without this line the code works and produce a pp compatible grid.
Do you have any idea if it possible to relax that assert?

Thanks for you kind help,

Alessio

[keileg]

The point is, the domain boundary has a hanging node (or rather edge) of sorts? I would not drop the relevant line, it serves as a sanity check in a computation which is far too complex compared to what we really would like (this is a case of an implementation which has snowballed beyond recognition). So while everything works out in your case, and in standard use, I really want to keep the check.

Have you tried to run the code ignoring assertions?