My question follows this first issue. I would like to mesh the fluidic pinball (see 2 below). It simply consists of 3 equally spaced circles, of same radius r. From left to right and up to down, we have the front cylinder, the top and the bottom cylinders. The geometry is symmetrical with respect to the y axis. I managed to mesh as I wished, except for one detail that I already wrote about in the 1.
Like the geometry, I would like the mesh to be symmetrical with respect to the y axis. So in the script, the mesh is first generated in the upper half of the geometry y>0, and then a symmetric operation with respect to the y axis is done to create the lower half mesh.
Here is the problem. As you can see in figure 2, in the upper half of the mesh, the area between front and top circles is refined by GMSH (probably due to the proximity of the circles) but this is not recreated between top and down circles. In the first issue 1, I tried to recreate that refinement using a geometry that was outside the upper half mesh (trying to create a distance field from the upper-half of bottom circle). But as stated in an answer of the issue, we cannot compute Distance/Threshold field from a geometry entity (in that case the upper half of the bottom circle) that does not lie inside the mesh (here the upper half mesh). To overcome that, I decided to create a Distance/Threshold field from the horizontal line (lines labeled 2 & 3 in the script) and transform it so it would correspond to the Distance/Threshold field from the upper half of the bottom circle. The transformation is based on MathEval field, that will transform the Distance field, applying then Threshold field. The solution I try to use is thus Distance/MathEval/Threshold. The MathEval transformation is based on the following calculus using the following geometrical illustration 3:
As I understand it, and following , a Distance field from horizontal lines (labeled 2& 3 in the script) should, for every point (green dot on the figure) of coordinate (x,y), compute the distance y (the distance from any point to the y axis is its y coordinate). This Distance field is then transformed using MathEval considering this relation (obtained applying Pythagorean theorem):
(d+r)^2=(y+c_y)^2+x^2.
Leading to ((y+c_y)^2+x^2)^{1/2}-r. In this latter equation, d would be, as I understand it, what would a Distance field from the upper half of the bottom cylinder produce.
So I am applying that in the script, and would like to observe a mesh refinement just below the top circle, same as the mesh refinement between front and top circles. But my script does not produce that kind of refinement below the top circle.
I hope I am clear in my explanations, and if you managed to understand them, would you have some solution for my problem ?
The mwe script is written below:
############################################
#------------ Mesh Creation ---------------#
############################################
import gmsh
import time
import sys
import numpy as np
import math
# cylinders radius
r = 0.5
# coordinates of the centers of the three cylinders (front, top and bottom)
c_x, c_y = [-3*math.sqrt(3)*r/2, 0, 0], [0, 3*r/2,-3*r/2]
# dimension of the mesh
gdim = 2
# Calculating the time that mesh generation will take
tic = time.perf_counter()
# Initialize and add model
gmsh.initialize()
gmsh.model.add("fluidic_pinball")
# target mesh size for any point
lc = 1
# mesh settings (used in the fields & geometry)
factors = {"cylinders": 0.0095, "cylinders_far": 1.5, "wake_start": 0.08, "wake_end": 0.3, "around_cylinders": 2}
############################################
#---------- Utility function --------------#
############################################
# Functions that are used to create the area around cylinders
def x_pm(x_c, delta_x, y_c, delta_y, R, pm):
if pm == "+":
sgn = 1
elif pm == "-":
sgn = -1
return (x_c-delta_x)/2 + sgn*(y_c-delta_y)/2*math.sqrt((4*math.pow(R, 2) - (math.pow(x_c-delta_x, 2) + math.pow(y_c-delta_y, 2)))/(math.pow(x_c-delta_x, 2)+math.pow(y_c-delta_y, 2)))+delta_x
def y_pm(x_c, delta_x, y_c, delta_y, R, pm):
if pm == "+":
sgn = -1
elif pm == "-":
sgn = 1
return (y_c-delta_y)/2 + sgn*(x_c-delta_x)/2*math.sqrt((4*math.pow(R, 2) - (math.pow(x_c-delta_x, 2) + math.pow(y_c-delta_y, 2)))/(math.pow(x_c-delta_x, 2)+math.pow(y_c-delta_y, 2)))+delta_y
def x_on_axis_pm(x_c, y_c, R, pm):
if pm == "+":
sgn = 1
elif pm == "-":
sgn = -1
return x_c + sgn*math.sqrt(math.pow(R, 2) - math.pow(y_c, 2))
############################################
#---------- Geometry creation -------------#
############################################
## 0 dim : Points definition
# Points to create the circles for the full and half cylinders
for i in range(2):
gmsh.model.geo.addPoint(c_x[i], c_y[i], 0, tag=1 + 4*i) # center
gmsh.model.geo.addPoint(c_x[i] + r, c_y[i], 0, tag=2 + 4*i)
gmsh.model.geo.addPoint(c_x[i], c_y[i] + r, 0, tag=3 + 4*i)
gmsh.model.geo.addPoint(c_x[i] - r, c_y[i], 0, tag=4 + 4*i)
gmsh.model.geo.addPoint(c_x[1], c_y[1] - r, 0, tag=9)
# Points to create the area around the cylinders
theta = 30 # angel that define the wake height level
assert 90 > theta > 0 # theta must be < 90° or the geometry will not be correct, and > 0 since gmsh can't create circle arcs >= 180°
gmsh.model.geo.addPoint(x_on_axis_pm(c_x[1], c_y[1], r*factors["around_cylinders"], "+"), 0, 0, tag=10)
gmsh.model.geo.addPoint(c_x[1] + r*factors["around_cylinders"], c_y[1], 0, tag=11)
gmsh.model.geo.addPoint(c_x[1] + factors["around_cylinders"]*r*math.cos(math.pi*theta/180), c_y[1] + factors["around_cylinders"]*r*math.sin(math.pi*theta/180), 0, tag=12)
gmsh.model.geo.addPoint(x_pm(c_x[1], c_x[0], c_y[1], c_y[0], r*factors["around_cylinders"], "-"), y_pm(c_x[1], c_x[0], c_y[1], c_y[0], r*factors["around_cylinders"], "-"), 0, tag=13)
gmsh.model.geo.addPoint(c_x[0] - factors["around_cylinders"]*r, c_y[0], 0, tag=18)
# Point to refine mesh locally
gmsh.model.geo.addPoint(0, 0, 0, tag=14)
## 1 dim: Lines and Curves
# outer boundary (without curves)
gmsh.model.geo.addLine(18, 4, tag=1)
gmsh.model.geo.addLine(2, 14, tag=2)
gmsh.model.geo.addLine(14, 10, tag=3)
# Creation of the circles for the three cylinders (4 CircleArcs for each circle)
# half front circle
gmsh.model.geo.addCircleArc(2, 1, 3, tag=4)
gmsh.model.geo.addCircleArc(3, 1, 4, tag=5)
# full top circle
gmsh.model.geo.addCircleArc(6, 5, 7, tag=6)
gmsh.model.geo.addCircleArc(7, 5, 8, tag=7)
gmsh.model.geo.addCircleArc(8, 5, 9, tag=8)
gmsh.model.geo.addCircleArc(9, 5, 6, tag=9)
# Creation of the circles arc for the area around cylinder
gmsh.model.geo.addCircleArc(10, 5, 11, tag=10)
gmsh.model.geo.addCircleArc(11, 5, 12, tag=11)
gmsh.model.geo.addCircleArc(12, 5, 13, tag=12)
gmsh.model.geo.addCircleArc(13, 1, 18, tag=13)
# # Line between cylinders
# gmsh.model.geo.synchronize()
############################################
#-------------- Plane surfaces ------------#
############################################
## 2 dim
# Creation of the plane surface: surface around cylinders
gmsh.model.geo.addCurveLoop([1, -5, -4, 2, 3, 10, 11, 12, 13], tag=1)
gmsh.model.geo.addCurveLoop([6, 7, 8, 9], tag=2)
gmsh.model.geo.addPlaneSurface([1, 2], tag=1)
gmsh.model.geo.synchronize()
############################################
#---------------- Fields ------------------#
############################################
# Fields
# Distance from the cylinders (front half and full top)
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "CurvesList", [4, 5, 6, 7, 8, 9])
gmsh.model.mesh.field.setNumber(1, "Sampling", 1000)
# Threshold field to refine the mesh around the three cylinder
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "InField", 1)
gmsh.model.mesh.field.setNumber(2, "SizeMin", lc*factors["cylinders"])
gmsh.model.mesh.field.setNumber(2, "SizeMax", lc*factors["cylinders_far"])
gmsh.model.mesh.field.setNumber(2, "DistMin", 0)
gmsh.model.mesh.field.setNumber(2, "DistMax", factors["around_cylinders"]*r)
# Distance from the cylinder (the bottom one, whose geometry is not created)
gmsh.model.mesh.field.add("Distance", 3)
gmsh.model.mesh.field.setNumbers(3, "CurvesList", [2, 3])
gmsh.model.mesh.field.setNumber(3, "Sampling", 1000)
# MathEval field
dict_arg = {"smax": lc*factors["cylinders_far"], "smin": lc*factors["cylinders"], "dmax": factors["around_cylinders"]*r, "dmin": 0, "c_y": c_y[1], "alpha": 1, "r": r}
gmsh.model.mesh.field.add("MathEval", 4)
gmsh.model.mesh.field.setString(4, "F", "(sqrt(F6+{c_y})^2+x^2)-{r}".format(**dict_arg)) # line
# Threshold field to refine the mesh around the three cylinder
gmsh.model.mesh.field.add("Threshold", 5)
gmsh.model.mesh.field.setNumber(5, "InField", 4)
gmsh.model.mesh.field.setNumber(5, "SizeMin", lc*factors["cylinders"])
gmsh.model.mesh.field.setNumber(5, "SizeMax", lc*factors["cylinders_far"])
gmsh.model.mesh.field.setNumber(5, "DistMin", 0)
gmsh.model.mesh.field.setNumber(5, "DistMax", factors["around_cylinders"]*r)
# Constant mesh size inside cylinder surroundings
gmsh.model.mesh.field.add("Constant", 6)
gmsh.model.mesh.field.setNumber(6, "VIn", lc*factors["wake_start"])
gmsh.model.mesh.field.setNumbers(6, "SurfacesList", [1])
# Final field as in t10.py, take the minimum mesh size for all Threshold fields.
gmsh.model.mesh.field.add("Min", 7)
gmsh.model.mesh.field.setNumbers(7, "FieldsList", [2, 5, 6])
gmsh.model.mesh.field.setAsBackgroundMesh(7)
############################################
#----- Mesh Saving & Post-processing ------#
############################################
# Generate and write to file
gmsh.model.mesh.generate(gdim)
############################################
#----------- Symmetry handling ------------#
############################################
# get the mesh data
m = {}
for e in gmsh.model.getEntities():
bnd = gmsh.model.getBoundary([e])
nod = gmsh.model.mesh.getNodes(e[0], e[1])
ele = gmsh.model.mesh.getElements(e[0], e[1])
m[e] = (bnd, nod, ele)
# transform the mesh and create new discrete entities to store it
def transform(m, offset_entity, offset_node, offset_element, tx, ty, tz):
for e in sorted(m):
gmsh.model.addDiscreteEntity(
e[0], e[1] + offset_entity,
[(abs(b[1]) + offset_entity) * math.copysign(1, b[1]) for b in m[e][0]])
coord = []
for i in range(0, len(m[e][1][1]), 3):
x = m[e][1][1][i] * tx
y = m[e][1][1][i + 1] * ty
z = m[e][1][1][i + 2] * tz
coord.append(x)
coord.append(y)
coord.append(z)
gmsh.model.mesh.addNodes(e[0], e[1] + offset_entity,
m[e][1][0] + offset_node, coord)
gmsh.model.mesh.addElements(e[0], e[1] + offset_entity, m[e][2][0],
[t + offset_element for t in m[e][2][1]],
[n + offset_node for n in m[e][2][2]])
if (tx * ty * tz) < 0: # reverse the orientation
gmsh.model.mesh.reverse([(e[0], e[1] + offset_entity)])
transform(m, 1000, 1000000, 1000000, 1, -1, 1)
# remove the duplicate nodes that will have been created on the internal
# boundaries
gmsh.model.mesh.removeDuplicateNodes()
gmsh.model.geo.synchronize()
############################################
#------------- Write & show ---------------#
############################################
gmsh.write("mirror_pinball_mesh.msh")
# gmsh.model.mesh.generate(gdim)
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.clear()
Also, do you know how to use LaTeX rendering in stackoverflow issues ? I would help make my post more clear.