Copy physical surfaces
Good morning,
In order to mesh the stator of an electrical motor, I define the geometry of a single stator slot, and then I copy and rotate it to create the full one.
Because I have many possible variation of the slot, to define my physical surface I coded a function that select all the surface inside a cylinder (by checking if each points is inside). Unfortunately it is very slow, so I tried 2 different methods:
Method 1
- Copy and rotate the initial surface of the slot
- Extrude all slots
- Select all the physical surfaces using
get_volume_in_cylinder
(very slow)
Method 2
- Extrude the first slot
- Copy the extruded volumes
- Adding each copied surface to its corresponding physical group
- Select all the physical surfaces for this slot using
get_volume_in_cylinder
- Rotate the extruded volumes (very slow)
In my opinion the fastest method would be to:
- Extrude the first slot
- Select all the physical surfaces for this slot using
get_volume_in_cylinder
- Copy and rotate the initial surface of the first slot
- Extrude all the other slots
- Adding each copied surface to its corresponding physical group
But step 5. seems really complicated because some physical surfaces are lateral and may be not created during the extrusion if they already exist for an adjacent slot. Plus there is no easy 1-1 correspondence like in method 2.
Am I taking this problem the wrong way? Do you any advice on how I could solve this?
def get_volumes_in_cylinder(origin, z_min=-inf, z_max=inf, R_min=0, R_max=inf, theta_min=-pi, theta_max=pi, vtol=1e-5, v_tags=None):
"""
Return tags of all the volumes contained in the 3D cylinder along z-axis
Warning: slow
@param origin : center of the cylinder
@param z_min : min z value of cylinder relative to origin
@param z_max : max z value of cylinder relative to origin
@param R_min : min radius of cylinder
@param R_max : max radius of cylinder
@param theta_min : min angle (if cylinder sector) relative to x-axis, between [-pi, pi]
@param theta_max : max angle (if cylinder sector) relative to x-axis, between [-pi, pi]
@param vtol : tolerance to be accepted inside the cylinder
@param v_tags : surfaces of volumes to test
"""
def is_in_cylinder(v_tag):
if R_min == 0 and R_max == inf and theta_min == -pi and theta_max == pi:
_, _, vol_z_min, _, _, vol_z_max = gmsh.model.getBoundingBox(3, v_tag)
if vol_z_min > z_min - vtol and vol_z_max < z_max + vtol:
return True
else:
for s_tag in gmsh.model.getAdjacencies(3, v_tag)[1]:
for l_tag in gmsh.model.getAdjacencies(2, s_tag)[1]:
for p_tag in gmsh.model.getAdjacencies(1, l_tag)[1]:
coord = get_point_coord(p_tag) - o
R = sqrt(coord[0] ** 2 + coord[1] ** 2)
if coord[2] > z_min - vtol and coord[2] < z_max + vtol:
if R > R_min - vtol and R < R_max + vtol:
if atan2(coord[1], coord[0]) > theta_min - R * vtol and atan2(coord[1], coord[0]) < theta_max + R * vtol:
continue # point inside cylinder
return False # point outside cylinder
return True
if is_int(origin): # if tag and not coordinates
o = get_point_coord(origin)
else:
o = np.array(origin)
gmsh.model.geo.synchronize()
res_v_tags = []
candidates = []
if v_tags is not None:
candidates = v_tags
# dim_tags = gmsh.model.getBoundary([(3, v_tag) for v_tag in v_tags], combined=True, oriented=False)
# candidates += [tag for _, tag in dim_tags]
if v_tags is None:
dim_tags = gmsh.model.getEntitiesInBoundingBox(o[0] - R_max - vtol, o[1] - R_max - vtol, o[2] + z_min - vtol, o[0] + R_max + vtol, o[1] + R_max + vtol, o[2] - z_max + vtol, dim=3)
candidates = [tag for _, tag in dim_tags]
for v_tag in candidates:
if is_in_cylinder(v_tag):
res_v_tags.append(v_tag)
return res_v_tags