Additional boundary points added

[kayarre]

It may be expected but I was hoping that the output would return the same points along the polygon without extra points.

It appears that additional boundary points are added to the edges of a polygon. given a polygon with 58 points and a custom edge length, the output adds additional points on the edge of the mesh. converting the mesh to vtk and extracting the edge with FeatureEdges, the point count on the outer edge is 97 in the case below.

Example:

import vtk
import numpy as np
import meshio
import dmsh

poly = [[1.1207656860351562, -1.474766731262207], [0.9311693906784058, -1.6106348037719727], 
 [0.7221899032592773, -1.7151135206222534], [0.5003705620765686, -1.7881356477737427],
 [0.24269625544548035, -1.8346717357635498], [0.039923492819070816, -1.8448166847229004],
 [-0.1640581488609314, -1.8532711267471313], [-0.3358287513256073, -1.8202623128890991],
 [-0.5591791868209839, -1.757102131843567], [-0.724500298500061, -1.7027467489242554],
 [-0.8599539995193481, -1.64866304397583], [-1.0511908531188965, -1.5159660577774048],
 [-1.2156579494476318, -1.3969435691833496], [-1.3622872829437256, -1.2560036182403564],
 [-1.505843997001648, -1.0733786821365356], [-1.6370562314987183, -0.8811505436897278],
 [-1.702138900756836, -0.718871533870697], [-1.7739341259002686, -0.5288633108139038],
 [-1.8255635499954224, -0.3319735527038574], [-1.8417296409606934, -0.10016066581010818],
 [-1.848059058189392, 0.1026785671710968], [-1.8210378885269165, 0.30386221408843994],
 [-1.783565640449524, 0.5339696407318115], [-1.7262778282165527, 0.6997944712638855],
 [-1.645498275756836, 0.8552301526069641], [-1.5460537672042847, 1.0332894325256348],
 [-1.4829825162887573, 1.0954675674438477], [-1.3618178367614746, 1.258979320526123],
 [-1.2313412427902222, 1.3749150037765503], [-1.1004729270935059, 1.4900223016738892],
 [-0.9327254295349121, 1.605852484703064], [-0.7736870050430298, 1.6783183813095093], 
 [-0.5828970074653625, 1.7475013732910156], [-0.4182337522506714, 1.8048688173294067],
 [-0.1892661303281784, 1.8489644527435303], [0.014770971611142159, 1.8476674556732178],
 [0.2173810601234436, 1.833501935005188], [0.41810983419418335, 1.8017268180847168],
 [0.5875660181045532, 1.7592521905899048], [0.7782935500144958, 1.6866511106491089],
 [0.9074407815933228, 1.6185364723205566], [1.0574780702590942, 1.5281968116760254],
 [1.2203333377838135, 1.4045077562332153], [1.354766607284546, 1.2504076957702637],
 [1.4904139041900635, 1.098829746246338], [1.6049500703811646, 0.9301034808158875],
 [1.711469292640686, 0.722557544708252], [1.7588528394699097, 0.5539127588272095],
 [1.811102271080017, 0.38747504353523254], [1.8434756994247437, 0.18628628551959991],
 [1.8454340696334839, -0.01729956641793251], [1.843611478805542, -0.19182713329792023],
 [1.8114262819290161, -0.3634224832057953], [1.7626954317092896, -0.5608553290367126],
 [1.684820294380188, -0.7485829591751099], [1.6166356801986694, -0.9093696475028992],
 [1.5086168050765991, -1.0823739767074585], [1.39910089969635, -1.2188986539840698]]

def inverse_circle_generator(R, N_edges):
    # need to check that the shape of the equation is correct
    # also need to add a check that the R is passed
    R_ = R
    c = 2.0*np.pi / N_edges
    N_mult = 2.8
    k1 = N_mult * c
    k2 = (N_mult-1.0) / N_mult

    def inverse_circle(p):
        r = np.sqrt(np.sum(np.square(p), axis=0))
        r[r > R_] = R_
        r[r < 0.0] = 0.0
        h = k1 * R_*(1.0 - k2 / R_ * r)
        return h

    return inverse_circle

geo = dmsh.Polygon(poly)
edge_length = inverse_circle_generator(1.850, len(poly)*1.5)
x, cells = dmsh.generate(geo, edge_length, delta_t=0.1, tol=1.0e-6)

mesh = meshio.Mesh(x, {"triangle":cells})
meshio.write("test.vtu", mesh)

print(len(poly), x.shape)

pts = vtk.vtkPoints()
polygon = vtk.vtkPolygon()
cells = vtk.vtkCellArray()
pd = vtk.vtkPolyData()
for idx, pt in enumerate(poly):
    pts.InsertNextPoint(pt[0], pt[1], 0.0)
    polygon.GetPointIds().InsertNextId(idx)
cells.InsertNextCell(polygon)
pd.SetPoints(pts)
pd.SetPolys(cells)

writer = vtk.vtkXMLPolyDataWriter()
writer.SetFileName("polygon.vtp")
writer.SetInputConnection(pd.GetOutputPort())
writer.SetInputData(pd)
writer.Write()

The red dots are the polygon points and the green dots are the edge points from dmsh.

It looks like it keeps the polygon points but adds additional points. is there a way to constrain it to not add additional points on the boundary? does gmsh provide this kind of constraint?

Thank you so much for your time.

[nschloe]

Please always post complete code that can be copied and pasted; I'm too lazy to fiddle it together myself. :)

[kayarre]

@nschloe No problem thank you.

[nschloe]

Well, okay, this is kind of expected. What dmsh essentially does is placing a regular grid above the domain, and cut off all excess points. The remaining cells are strictly inside the domain, but not on the boundary, so then dmsh "inflates" the mesh until the outer points sit on the boundary. (Something like that.) The reason why the resulting meshes are so good is that most of the time there's no need for much inflation.

Okay, so as you can see, there's no way to strictly control the boundary points in this approach. What you can try though is to increase the initial edge length/cell size. Chosen appropriately, and with a little bit of trial and error, I think you might just get what you want.

[kayarre]

That is super helpful. Thank you so much. ~Kurt