block_diagonal fails with meshplex 1.16.1

[rob-rb]

When I call get_new_points I receive the traceback:

Traceback (most recent call last): File "/opt/project/lib/adms/utils/optimesh_entry.py", line 195, in _optimize new_points = get_new_points(mesh) File "/opt/conda/lib/python3.8/site-packages/optimesh/cvt/block_diagonal.py", line 73, in get_new_points rhs = -jac_uniform(mesh, mask) File "/opt/conda/lib/python3.8/site-packages/optimesh/cvt/_helpers.py", line 4, in jac_uniform cvc = mesh.get_control_volume_centroids(cell_mask=mask) File "/opt/conda/lib/python3.8/site-packages/meshplex/_mesh.py", line 235, in get_control_volume_centroids assert cell_mask.shape == (self.idx[-1].shape[-1],) AssertionError

[nschloe]

I need an MWE to reproduce.

[rob-rb]

import numpy as np
import npx

import meshplex
import optimesh

def control_mesh(mesh):
    s = mesh.idx_hierarchy.shape
    a = np.sort(mesh.idx_hierarchy.reshape(s[0], -1).T)
    _, cts = npx.unique_rows(a, return_counts=True)
    if cts.max() > 2:
        print("control_mesh found non-manifold edges", cts.max())
        for idx in np.where(cts > 2)[0]:
            print(idx, cts[idx])
    else:
        print("good_mesh")

points = [
    [0.0, 0.0, 0.0],
    [1.0, 0.5, 0.0],
    [1.0, 0.76, 1.0],
    [1.0, 0.8, 1.0],
    [1.0, 0.7, 1.0],
    [1.0, 0.7, 0.0],
]
cells = [[2, 0, 1], [2, 3, 0], [4, 0, 3], [4, 3, 5], [2, 5, 3]]

mesh = meshplex.MeshTri(points, cells)

mesh.save("out0.vtu")

control_mesh(mesh)

print(mesh.cells("points"))
print("flip 1...")
mesh.flip_until_delaunay(max_steps=3)
mesh.save("out1.vtu")

print(mesh.cells("points"))
n = mesh.remove_duplicate_cells()
print(f"{n = }")

print("flip 2...")
mesh.flip_until_delaunay(max_steps=1)
mesh.save("out2.vtu")

print(mesh.cells("points"))

print("flip 3..")
mesh.flip_until_delaunay(max_steps=1)
mesh.save("out3.vtu")

optimesh.get_new_points(mesh, "CVT diaognal")
control_mesh(mesh)

mesh = meshplex.MeshTri(mesh.points, mesh.cells("points"))
mesh.save("out4.vtu")

print(mesh.cells("points"))

[rob-rb]

it seems that rebuilding the mesh, before callins optimesh, solves the problem.

mesh = meshplex.MeshTri(points, mesh.cells("points"))

[nschloe]

I can't reproduce locally, so perhaps I fixed it this morning. Try the new meshplex 0.16.2.

[rob-rb]

yep

[nschloe]

Awesome! Closing then. Feel free to reopen if the issue reappears at some point.