Boundary points of 2D geometry

[DeepaMahm]

Hi @marcomusy,

I want to find the coordinates of the boundary points to specify boundary conditions for 2D geometry like the below (geometry generated by specifying point cloud)

I am unsure how to get the coordinates of boundary points for the two ends of the geometry. On one end , I have to specify the Dirichlet and Neumann on the other end in deepxde library; so I need the boundary coordinates separately.

Could you please give some suggestions on how we could get the boundary coordinates in vedo?

[marcomusy]

Hi, you can try to use https://vedo.embl.es/autodocs/content/vedo/vedo/pointcloud.html#delaunay2d and then https://vedo.embl.es/autodocs/content/vedo/vedo/mesh.html#Mesh.boundaries

[DeepaMahm]

Hi @marcomusy,

Sorry for the late reply.

Thanks a lot for the suggestions. I could use delaunay and mesh boundaries; it works great!

I could get all the boundary points (colored red in the above image). I would like to get the indices/labels of the boundary points which lie within the square boxes.

Could you please suggest how we can get those labels?

I tried the following to label all boundary points

boundary_pts = Points(bpts[pids], r=10, c='red5')
  # Create a Label object for all the vertices in the mesh
  # labels = boundary_pts.labels('id', scale=1).c('green2')
  labels = boundary_pts.labels('id').z(0.01)

complete code:

def get_boundary_points(pts):
    """
    uses delaunay from vedo to generate a mesh

    :return:
    refs:
    https://github.com/marcomusy/vedo/blob/master/examples/basic/delaunay2d.py
    https://github.com/marcomusy/vedo/blob/master/examples/basic/delaunay2d.py
    """

    # Use the Delaunay triangulation algorithm to create a 2D mesh
    # from the points in gp, with the given boundary ids
    b = delaunay2d(pts, mode='xy').c('w').lc('o').lw(1)

    # Create labels for the point ids and set their z to 0.01
    labels = pts.labels('id').z(0.01)

    # vplt.show(pts, labels, dly, __doc__, bg="Mint").close()

    b.compute_normals().clean().linewidth(0.1)

    # Get the point IDs on the boundary of the mesh
    pids = b.boundaries(return_point_ids=True)

    # Get the points corresponding to the boundary IDs
    bpts = b.points()

    # Create a Points object to represent the boundary points
    boundary_pts = Points(bpts[pids], r=10, c='red5')

    # Create a Label object for all the vertices in the mesh
    # labels = boundary_pts.labels('id', scale=1).c('green2')
    labels = boundary_pts.labels('id').z(0.01)
    print(type(labels))

    print('vertices:', boundary_pts.vertices)

    # Show the mesh, boundary points, vertex labels, and docstring
    vplt.show(b, boundary_pts, labels, zoom=2).close()

But I am not sure how to print the lables / get the labels of selective boundary points

[marcomusy]

Hi Deepa, you can create a small 3D box and use box.inside_points(). if you already know their positions the ids are in the same order as the points/vertices.

Btw, If you upgrade you vedo version you will need to change b = delaunay2d(pts, mode='xy').c('w').lc('o').lw(1) to b = pts.generate_delaunay2d(mode='xy').c('w').lc('o').lw(1)

[DeepaMahm]

Thank you, I have the positions but those are the coordinates of all the boundary points. Could you please let me know how to create the 3D box after plotting the geometry (show) ? Can we mark the box interatively in the plot window that appears?

Updated : b = pts.generate_delaunay2d(mode='xy').c('w').lc('o').lw(1), thank you.

[marcomusy]

in your code pids and pts.vertices are 2 aligned numpy array if i understand correctly so you can extract/slice the pids array the way you like?

how to create the 3D box after plotting the geometry (show) ? Can we mark the box interactively in the plot window that appears?

yes, you need to create a callback function that generates a new box and place it where you click.. unless you have to do this repeatedly it's probably overkilling :) i can help you with that if needed.

[DeepaMahm]

in your code pids and pts.vertices are 2 aligned numpy array if i understand correctly so you can extract/slice the pids array the way you like?

I tried to visually readout the vertex labels of boundary points on one end of the geometry (to use pids and pts.vertices and get the vertices of selective boundary pids). The label overlaps and I couldn't read it clearly (unfortunately print(lables) didn't work, since it is a mesh object) and by changing scale (labels = boundary_pts.labels('id', scale=0.5).z(0.01)) I couldn't reduce the font size of the text that appears for labels.

yes, you need to create a callback function that generates a new box and place it where you click.. unless you have to do this repeatedly it's probably overkilling :) i can help you with that if needed.

Agree :) For now , I need to do this only once. Thank you

[marcomusy]

try labels2d() instead.

[DeepaMahm]

Thank you, the labels are readable now.

Can we pass a list of integers for creating the labels?

When I try to create the labels just for the boundary points (boundary_pts.labels2d) it labeling is from 0 to ..

labels = boundary_pts.labels2d('id', on="points", scale=0.5)

and when I try to plot the labels of all points labels = pts.labels2d('id')

the text is not readable.

[marcomusy]

Consider this option:

from vedo import *

circle = Circle(res=36)

mesh = circle.generate_mesh()
labels = mesh.labels2d('id')

# Add a pointdata array to the mesh to store the point IDs
mesh.add_ids()
mesh.color('green9').lc("white")

mboundaries = mesh.boundaries()
print(mboundaries)

pids = mboundaries.pointdata["PointID"]
print("PIDS   = ", pids)
print("POINTS = ", mesh.vertices[pids])

show(mesh, mboundaries, labels).close()

[DeepaMahm]

Hi @marcomusy, Thank you again.

Unfortunately, when I try for my use case the label text is not readable when all labels are displayed and I am not able to zoom in (screen hangs).

Is there a way to display the labels of just the boundary points?

If I can view the labels, then I could use mesh.vertices[pids_selective_bpts] to get the vertices of the boundary points that I want to select (assuming pids = labels).

[marcomusy]

Yes, try

from vedo import *

circle = Circle(res=36)
mesh = circle.generate_mesh()

# Add a pointdata array to the mesh to store the point IDs
mesh.add_ids()
mesh.color('green9').lc("white")

mboundaries = mesh.boundaries()
print(mboundaries)

pids = mboundaries.pointdata["PointID"]
print("PIDS   = ", pids)
print("POINTS = ", mesh.vertices[pids])
labels = mboundaries.labels("PointID")

show(mesh, mboundaries, labels).close()

[DeepaMahm]

Thank you so much.

mesh.add_ids() helps!

Unfortunately, I've the same issue again; the labels aren't clearly readable.

I tried to check if labels = boundary_pts.labels2d('id', on="points", scale=0.5)

works with mesh.add_ids(). But labels2d add labels for the boundary points from 0 ...

Could you please suggest if there is an alternative? Thanks so much for your kind attention

[marcomusy]

No worries, in fact there is a bug in labels2d(), add

settings.use_parallel_projection = False

then make the labels smaller so they dont overlap eg scale=0.25.

[DeepaMahm]

Hello @marcomusy ,

Instead of circle, I am now trying to create a rectangle geometry.

I tried,

from vedo import *

 rectangle = Rectangle([0.0, 0.0], [5.0, 1.0])
 mesh = rectangle.generate_mesh()

# Add a pointdata array to the mesh to store the point IDs
mesh.add_ids()
mesh.color('green9').lc("white")

mboundaries = mesh.boundaries()
print(mboundaries)

pids = mboundaries.pointdata["PointID"]
print("PIDS   = ", pids)
print("POINTS = ", mesh.vertices[pids])
labels = mboundaries.labels("PointID")

show(mesh, mboundaries, labels).close()

However, print( mesh.vertices) displays only 5 points and I see a warning (vtkDelaunay2D.cxx:1419 WARN| vtkDelaunay2D (00000175765CA0F0): Edge not recovered, polygon fill suspect)

[[2.5012121 1.014475  0.       ]
 [5.        1.        0.       ]
 [5.        0.        0.       ]
 [0.        0.        0.       ]
 [0.        1.        0.       ]]

whereas, while using Circle, I could get all the points on the surface (vtkDelaunay2D).

Could you please let me know how to create `Rectangle and generate_mesh()?

I also tried,

polygon = [
         [(82, 92, 47), (87, 88, 47),  # x,y,z of vertices
          (93, 95, 47), (88, 99, 47)],
         [[0, 1, 2, 3]],  # vertex connectivity
     ]
 mesh = Mesh(polygon)

Unfortunately, this also returns only 5 points when I try mesh.vertices.

[marcomusy]

Hi, what about triangulating a simple grid:

from vedo import *

mesh = Grid(s=[np.arange(0,5, 0.2), np.arange(0,1, 0.2)]).triangulate()

# Add a pointdata array to the mesh to store the point IDs
mesh.add_ids()
mesh.color('green9').lc("white")

mboundaries = mesh.boundaries()
print(mboundaries)

pids = mboundaries.pointdata["PointID"]
print("PIDS   = ", pids)
print("POINTS = ", mesh.vertices[pids])
labels = mboundaries.labels("PointID")

show(mesh, mboundaries, labels).close()

[DeepaMahm]

Thanks a lot; it works great!

I could get

Could you please let me know if we can create fine mesh? i.e. I would like to get more points/vertices in the domain and boundary

[marcomusy]

Hi, you can change the numpy array to any spacing np.arange(0,5, x)

[DeepaMahm]

Hi @marcomusy ,

For a geometry like the following

I would like to get all the boundaries (please refer to the lines marked with black and indicated with arrows).

Using,

b = pts.generate_delaunay2d(mode='xy').c('w').lc('o').lw(1)
b.add_ids()
#vplt.show(pts, labels, dly, __doc__, bg="Mint").close()
b.compute_normals().clean().linewidth(0.1)
boundary_pts = Points(b.points()[pids], r=10, c='red5')

I could obtain the points marked with red dots as boundaries.

Could you please help me with obtaining all the boundary points (i.e. including boundaries inside the domain)?

The original geometry is added here (image) for your kind reference and I would like to get the boundary points of all the vessel branches. The pointcloud of the geometry is included here (pointcloud), kindly see.