interior edges

[rob-rb]

Dear Nico,

I want to write some curvature measure on interior edges.
For my mesh I have the following situation:

edges = mesh.edges["points"]
print(edges.shape)
(229941, 2)
print(mesh.is_boundary_edge.shape)
(230036,)

Where is the difference?
Is this intendend or a bug?

[nschloe]

Looks like a bug. Details?

[rob-rb]

import meshplex

mesh = meshplex._reader.read(f"./1_500_contour.obj")

def boundary_stats():
    print("b_points", mesh.points[mesh.is_boundary_point].shape)
    volumes = mesh.cell_volumes[mesh.is_boundary_cell]
    print("b_cells", volumes.shape)
    edges = mesh.edges["points"]
    print("b_edges", edges[mesh.is_boundary_edge].shape)


def callback(arg):
    return mesh.cell_volumes[mesh.is_boundary_cell] < 1e-5

boundary_stats()

while True:
    res = mesh.remove_boundary_cells(callback)
    if not res:
        break
    print("deleted", res, mesh.remove_dangling_points())

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

boundary_stats()

[rob-rb]

uncommenting line 24 -> no error

[rob-rb]

Actually I need a remove triangle method filling up the holes. Something like this:

import meshplex
import numpy as np

def join_short_edges(mesh, edge_length):
    edges = mesh.edges["points"]
    points = mesh.points
    boundary_edges = edges[mesh.is_boundary_edge]
    boundary_edge_links = points[boundary_edges]
    boundary_diffs = boundary_edge_links[:, 0] - boundary_edge_links[:,1]
    boundary_link_lengths = np.linalg.norm(boundary_diffs, axis=-1)
    short_edges = boundary_edges[boundary_link_lengths < edge_length]
    mid_points = .5*(points[short_edges[:, 0]] + points[short_edges[:, 1]])
    points = np.array(points)
    points[short_edges[:, 0]] = mid_points
    mappings = {}
    for k1, k2 in short_edges:
        mappings[k2] = k1
    boundary_cells = mesh.cells("points")[mesh.is_boundary_cell]
    dels = []
    for index, f in enumerate(boundary_cells):
        touched = False
        for j, k in enumerate(f):
            if k in mappings:
                f[j] = mappings[k]
                touched = True
        if touched:
            f = set(f)
            if len(f) < 3:
                dels.append(index)
    print("mappings", mappings)
    print("collapsing", len(dels), len(short_edges))
    inner_cells = mesh.cells("points")[np.logical_not(mesh.is_boundary_cell)]
    good_boundary_cells = [b for j, b in enumerate(boundary_cells) if not j in dels]
    faces = list(inner_cells)
    faces.extend(good_boundary_cells)
    return meshplex.MeshTri(points, faces)

#mesh = meshplex._reader.read(f"./1_500_contour.obj")
from math import *

def easy_mesh(n=4):
    coords = [(0,0)]
    n = 4
    triangles = []
    a =0
    for j in range(n):
        triangles.append((0, j+1, j+2))
        coords.append((cos(a), sin(a)))
        a += 2*pi/n


    a = -.1
    co = (cos(a), sin(a))
    coords.append(co)
    triangles.append((0, len(coords)-1, len(coords)))

    a = -.05
    co = (cos(a), sin(a))
    coords.append(co)
    triangles.append((0, len(coords)-1, 1))
    return meshplex.MeshTri(coords, triangles)

mesh = easy_mesh(4)
mesh.show()
mesh = join_short_edges(mesh, .1)
mesh.remove_dangling_points()
mesh = join_short_edges(mesh,.1)
#mesh.remove_cells([len(triangles)-1])
mesh.show()