Constrained Delaunay Triangulation in 2D
Dear Gmsh developers,
many thanks for this awesome piece of software.
I am trying to conduct a constrained Delaunay Triangulation much like the one Triangle is offering.
My toy problem looks as follows:
As input I have a square with some lines all constructed by line elements. These line elements are supposed to appear as a single edges in the final mesh.
To do this I tried a couple approaches like the following suggestions:
- issue1318 -> does not work, since the "Initial Mesh Only"-Option does not insert new points and therefore gives only very poor meshes
- issue1708 -> did not work since I could not manage to embed line elements (issue 1708 deals with point constraints only)
As these issues suggest, the best course of action is to build the initial line mesh via discrete entities. However, I was not able to include / embed the line elements inside the domain (so Line 7,11,12,8,9 and 10 according to my example).
Instead I "missused" the transfinite feature by building up the model via geometric objects like lines, points and planar surfaces, whereas every used line is restricted to only hold two nodes.
My solution looks like this:
import numpy as np
import gmsh
ele_size = 0.05
nds = np.array([[0. ,0.],
[0.25,0.],
[0.75,0.],
[1. ,0.],
[1. , 1.],
[0. , 1.],
[0.5,0.3],
[0.1,0.1],
[0.2,0.2],
[0.5,0.5],
[0.7,0.5],
[0.6,0.6],
])
n_nds = nds.shape[0]
outer_cycle = np.array([0,1,2,3,4,5])
boundary_tri_lines = np.array([[1,6],
[6,2]])
single_line = np.array([[7,8]])
inner_tri_lines = np.array([[9,10],
[10,11],
[11,9]])
gmsh.initialize()
gmsh.logger.start()
# create domain / outer loop
for i,nd in enumerate(nds):
gmsh.model.geo.addPoint(nd[0], nd[1], 0,i+1 )
outer_cycle_line_tags = []
for line in np.vstack((outer_cycle,np.roll(outer_cycle,-1))).T+1:
outer_cycle_line_tags.append( gmsh.model.geo.addLine(line[0],line[1]) )
gmsh.model.geo.mesh.setTransfiniteCurve( outer_cycle_line_tags[-1],2 )
outer_curved_loop_tag = gmsh.model.geo.addCurveLoop( outer_cycle_line_tags, 1)
domain = gmsh.model.geo.addPlaneSurface([outer_curved_loop_tag], 1)
gmsh.model.geo.synchronize()
# EMBED CONSTRAINED LINES INSIDE THE DOMAIN
for line in np.vstack((single_line,inner_tri_lines,boundary_tri_lines)):
line_tag = gmsh.model.geo.addLine(line[0]+1,line[1]+1)
gmsh.model.geo.mesh.setTransfiniteCurve( line_tag,2 )
gmsh.model.geo.synchronize()
gmsh.model.mesh.embed(1, [line_tag], 2, domain)
gmsh.model.mesh.setSize(gmsh.model.getEntities(0), ele_size)
gmsh.option.setNumber("Mesh.Algorithm", 5) # 5 Delaunay # 6 Frontal-Delaunay
gmsh.option.setNumber("Mesh.MeshSizeMax", ele_size)
gmsh.model.mesh.generate(2)
gmsh.model.geo.synchronize()
gmsh.fltk.run()
gmsh.finalize()
which yields:
Now I am wondering, why there are no elements inserted in the big and small "initial" triangles? I would have suspected, that at least on additional node would have been inserted in the center of each triangle to split the triangles into three smaller elements.
Furthermore I am curious whether my method would be the intended one for a constrained Delaunay triangulation with gmsh?
Any comments are very much appreciated.