class Geometry:
def __init__(self, data: dF.FDM() = None, geom: FiQuSGeometry() = None,
geom_folder: str = None, verbose: bool = False):
"""
Class to generate geometry
:param data: FiQuS data model
:param geom: ROXIE geometry data
:param verbose: If True more information is printed in python console.
"""
self.data: dF.FDM() = data
self.geom: FiQuSGeometry() = geom.Roxie_Data
# move cooling holes to a desired position
if self.data.magnet.solve.thermal.collar_cooling.move_cooling_holes:
self.geom.iron.key_points = self.move_keypoints(self.geom.iron.key_points, self.data.magnet.solve.thermal.collar_cooling.move_cooling_holes)
self.geom_folder = geom_folder
self.verbose: bool = verbose
self.md = dM.MultipoleData()
self.gu = GmshUtils(self.geom_folder, self.verbose)
self.gu.initialize(verbosity_Gmsh=self.data.run.verbosity_Gmsh)
self.occ = gmsh.model.occ
self.model_file = os.path.join(self.geom_folder, self.data.general.magnet_name)
self.blk_ins_lines = {} # for meshed insulation
self.ins_wire_lines = {} # for meshed insulation
self.block_coil_mid_pole_blks = {}
self.nc = {'collar': 'c', 'iron_yoke': 'i', 'poles': 'p'}
self.inv_nc = {v: k for k, v in self.nc.items()} #invert naming convention
if self.data.magnet.geometry.electromagnetics.symmetry != 'none':
self.symmetric_loop_lines = {'x': [], 'y': []}
self.symmetric_bnds = {'x_p': {'pnts': [], 'line_pnts': []}, 'y_p': {'pnts': [], 'line_pnts': []},
'x_n': {'pnts': [], 'line_pnts': []}, 'y_n': {'pnts': [], 'line_pnts': []}}
def clear(self):
self.md = dM.MultipoleData()
self.block_coil_mid_pole_blks = {}
gmsh.clear()
def ending_step(self, gui: bool = False):
if gui:
self.gu.launch_interactive_GUI()
else:
if gmsh.isInitialized():
gmsh.clear()
gmsh.finalize()
def saveHalfTurnCornerPositions(self):
self.occ.synchronize()
iH, iL, oH, oL, iHr, iLr, oHr, oLr = [], [], [], [], [], [], [], []
for po in self.geom.coil.physical_order:
block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[
po.winding].blocks[po.block]
for halfTurn_nr, halfTurn in block.half_turns.items():
ht = halfTurn.corners.insulated
ht_b = halfTurn.corners.bare
iHr.append([ht_b.iH.x, ht_b.iH.y])
iLr.append([ht_b.iL.x, ht_b.iL.y])
oHr.append([ht_b.oH.x, ht_b.oH.y])
oLr.append([ht_b.oL.x, ht_b.oL.y])
iH.append([ht.iH.x, ht.iH.y])
iL.append([ht.iL.x, ht.iL.y])
oH.append([ht.oH.x, ht.oH.y])
oL.append([ht.oL.x, ht.oL.y])
with open(f"{self.model_file}.crns", 'w') as f:
json.dump({'iH': iH, 'iL': iL, 'oH': oH, 'oL': oL,
'iHr': iHr, 'iLr': iLr, 'oHr': oHr, 'oLr': oLr}, f)
def saveStrandPositions(self, run_type):
symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'
ht_nr = 0
std_nr = 0
parser_x, parser_y, blocks, ht, std, pole_blocks = [], [], [], [], [], []
for po in self.geom.coil.physical_order:
block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[
po.winding].blocks[po.block]
if po.pole == 1: pole_blocks.append(po.block)
for halfTurn_nr, halfTurn in block.half_turns.items():
ht_nr += 1
for strand_group_nr, strand_group in halfTurn.strand_groups.items():
for strand_nr, strand in strand_group.strand_positions.items():
std_nr += 1
blocks.append(po.block)
ht.append(ht_nr)
std.append(std_nr)
parser_x.append(strand.x)
parser_y.append(strand.y)
mirrored = {}
condition = {2: [1, -1], 3: [1, 1], 4: [-1, 1]}
if symmetry == 'xy': mirroring = {2: [-1, 1], 3: [-1, -1], 4: [1, -1]}
elif symmetry == 'x': mirroring = {3: [1, -1], 4: [1, -1]}
elif symmetry == 'y': mirroring = {2: [-1, 1], 3: [-1, 1]}
else: mirroring = {}
if mirroring:
df = pd.DataFrame({'parser_x': parser_x, 'parser_y': parser_y}, index=std)
for qdr, mrr in mirroring.items():
subdf = df[(condition[qdr][0] * df['parser_x'] < 0) & (condition[qdr][1] * df['parser_y'] < 0)]
for strand, x, y in zip(subdf.index, subdf['parser_x'], subdf['parser_y']):
mirrored[strand] = df[(df['parser_x'] == mrr[0] * x) & (df['parser_y'] == mrr[1] * y)].index.item()
with open(f"{self.model_file}_{run_type}.strs", 'w') as f:
json.dump({'x': parser_x, 'y': parser_y, 'block': blocks, 'ht': ht, 'mirrored': mirrored,
'pole_1_blocks': pole_blocks, 'poles': len(self.geom.coil.coils[1].poles)},
f)
def saveBoundaryRepresentationFile(self, run_type):
self.occ.synchronize()
gmsh.write(f'{self.model_file}_{run_type}.brep')
gmsh.clear()
def loadBoundaryRepresentationFile(self, run_type):
gmsh.option.setString('Geometry.OCCTargetUnit', 'M') # set units to meters
gmsh.open(os.path.join(f'{self.model_file}_{run_type}.brep'))
def saveAuxiliaryFile(self, run_type):
Util.write_data_to_yaml(f'{self.model_file}_{run_type}.aux', self.md.model_dump())
@staticmethod
def findMidLayerPoint(bc_current, bc_next, center, mean_rad):
mid_layer = [(bc_current.x + bc_next.x) / 2, (bc_current.y + bc_next.y) / 2]
mid_rad = Func.points_distance(mid_layer, [center.x, center.y])
dist_from_mid = mean_rad - mid_rad
angle = Func.arc_angle_between_point_and_abscissa(mid_layer, [center.x, center.y])
mid_layer[0] += dist_from_mid * np.cos(angle)
mid_layer[1] += dist_from_mid * np.sin(angle)
return mid_layer
@staticmethod
def getMidLayerEndpoints(el_current, el_next, center, mid_layer_arc_pnt=None, coil_type='cos-theta', cable_type='Rutherford', is_for_mid_pole=False):
thin_shell_endpoints = {'higher': list, 'lower': list}
which_block = {'higher': str, 'lower': str}
angles = {'higher': float, 'lower': float}
# Check if the element crosses the x axis
angles_to_correct = []
correction_angle = 0
l_curr = Func.arc_angle_between_point_and_abscissa([el_current.iL.x, el_current.iL.y], center)
h_curr = Func.arc_angle_between_point_and_abscissa([el_current.iH.x, el_current.iH.y], center)
l_next = Func.arc_angle_between_point_and_abscissa([el_next.iL.x, el_next.iL.y], center)
h_next = Func.arc_angle_between_point_and_abscissa([el_next.iH.x, el_next.iH.y], center)
if abs(l_curr - h_curr) > np.pi:
angles_to_correct.append('current')
correction_angle = max(1.05 * (2 * np.pi - l_curr), correction_angle)
if abs(l_next - h_next) > np.pi:
angles_to_correct.append('next')
correction_angle = max(1.05 * (2 * np.pi - l_next), correction_angle)
for side in thin_shell_endpoints.keys():
if mid_layer_arc_pnt:
if side == 'higher':
mid_lyr_curr, mid_lyr_next = [el_current.oH, el_current.iH], [el_next.oH, el_next.iH]
else:
mid_lyr_curr, mid_lyr_next = [el_current.oL, el_current.iL], [el_next.oL, el_next.iL]
if cable_type in ['Mono', 'Ribbon']:
pnts_curr = Func.intersection_between_circle_and_line(
Func.line_through_two_points([mid_lyr_curr[0].x, mid_lyr_curr[0].y], [mid_lyr_curr[1].x, mid_lyr_curr[1].y]),
[center, mid_layer_arc_pnt])
pnt_curr = pnts_curr[0] if Func.points_distance(pnts_curr[0], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) <\
Func.points_distance(pnts_curr[1], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) else pnts_curr[1]
pnts_next = Func.intersection_between_circle_and_line(
Func.line_through_two_points([mid_lyr_next[0].x, mid_lyr_next[0].y], [mid_lyr_next[1].x, mid_lyr_next[1].y]),
[center, mid_layer_arc_pnt])
pnt_next = pnts_next[0] if Func.points_distance(pnts_next[0], [mid_lyr_next[0].x, mid_lyr_next[0].y]) <\
Func.points_distance(pnts_next[1], [mid_lyr_next[0].x, mid_lyr_next[0].y]) else pnts_next[1]
elif cable_type == 'Rutherford':
pnt_curr = Func.intersection_between_circle_and_line(
Func.line_through_two_points([mid_lyr_curr[0].x, mid_lyr_curr[0].y], [mid_lyr_curr[1].x, mid_lyr_curr[1].y]),
[center, mid_layer_arc_pnt], get_only_closest=True)[0]
pnt_next = Func.intersection_between_circle_and_line(
Func.line_through_two_points([mid_lyr_next[0].x, mid_lyr_next[0].y], [mid_lyr_next[1].x, mid_lyr_next[1].y]),
[center, mid_layer_arc_pnt], get_only_closest=True)[0]
else:
if cable_type == 'Rutherford':
if coil_type == 'common-block-coil':
mid_layer_x = (el_current.oH.x + el_next.iH.x) / 2
if side == 'higher':
pnt_curr, pnt_next = [mid_layer_x, el_current.iH.y], [mid_layer_x, el_next.iH.y]
else:
pnt_curr, pnt_next = [mid_layer_x, el_current.iL.y], [mid_layer_x, el_next.iL.y]
else:
mid_layer_y = (el_current.iH.y + el_next.iH.y) / 2 if is_for_mid_pole else (el_current.oH.y + el_next.iH.y) / 2
if side == 'higher':
pnt_curr, pnt_next = [el_current.iH.x, mid_layer_y], [el_next.iL.x if is_for_mid_pole else el_next.iH.x, mid_layer_y]
else:
pnt_curr, pnt_next = [el_current.iL.x, mid_layer_y], [el_next.iH.x if is_for_mid_pole else el_next.iL.x, mid_layer_y]
elif cable_type in ['Mono', 'Ribbon']:
pnt_curr = [(el_current.oH.x + el_next.iH.x) / 2, (el_current.oH.y + el_next.iH.y) / 2] if side == 'higher'\
else [(el_current.oL.x + el_next.iL.x) / 2, (el_current.oL.y + el_next.iL.y) / 2]
pnt_next = pnt_curr
angle_curr = Func.arc_angle_between_point_and_abscissa(pnt_curr, center)
angle_next = Func.arc_angle_between_point_and_abscissa(pnt_next, center)
if 'current' in angles_to_correct:
angle_curr = angle_curr + correction_angle - (2 * np.pi if side == 'lower' else 0)
elif 'next' in angles_to_correct:
if angle_curr < np.pi / 2: angle_curr += correction_angle
elif angle_curr > np.pi * 3 / 2: angle_curr = angle_curr + correction_angle - 2 * np.pi
if 'next' in angles_to_correct:
angle_next = angle_next + correction_angle - (2 * np.pi if side == 'lower' else 0)
elif 'current' in angles_to_correct:
if angle_next < np.pi / 2: angle_next += correction_angle
elif angle_next > np.pi * 3 / 2: angle_next = angle_next + correction_angle - 2 * np.pi
if abs(angle_curr - angle_next) < 1e-6:
thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'
elif angle_curr * (-1 if side == 'lower' else 1) < angle_next * (-1 if side == 'lower' else 1):
thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'
else:
thin_shell_endpoints[side], angles[side], which_block[side] = pnt_next, angle_next, 'next'
if angles['higher'] < angles['lower']: return None
else: return thin_shell_endpoints, which_block
def create_geom_dict(self, geometry_setting):
return {v: k in geometry_setting.areas for k, v in self.nc.items()}
def constructIronGeometry(self, symmetry, geometry_setting, run_type):
"""
Generates points, hyper lines, and curve loops for the iron yoke
"""
iron = self.geom.iron #roxie
if symmetry == 'xy':
self.md.geometries.iron.quadrants = {1: dM.Region()}
list_bnds = ['x_p', 'y_p']
elif symmetry == 'x':
self.md.geometries.iron.quadrants = {1: dM.Region(), 2: dM.Region()}
list_bnds = ['x_p', 'x_n']
elif symmetry == 'y':
self.md.geometries.iron.quadrants = {1: dM.Region(), 4: dM.Region()}
list_bnds = ['y_p', 'y_n']
else:
for k in self.nc.keys(): getattr(self.md.geometries, k).quadrants = {1: dM.Region(), 2: dM.Region(), 4: dM.Region(), 3: dM.Region()}
list_bnds = []
lc = 1e-2
geom_dict = self.create_geom_dict(geometry_setting)
for point_name, point in iron.key_points.items():
identifier = next((k for k in self.inv_nc.keys() if re.match(f'^{k}', point_name[2:])), None)
if not geom_dict.get(identifier, False): continue
quadrants = getattr(self.md.geometries, self.inv_nc[identifier]).quadrants #re.sub(r'\d+', '', point_name[2:])
if symmetry in ['x', 'xy']:
if point.y == 0.:
self.symmetric_bnds['x_p']['pnts'].append([point_name, point.x])
if symmetry in ['y', 'xy']:
if point.x == 0.:
self.symmetric_bnds['y_p']['pnts'].append([point_name, point.y])
quadrants[1].points[point_name] = self.occ.addPoint(point.x, point.y, 0, lc)
if symmetry in ['x', 'none']:
if point.x == 0.:
quadrants[2].points[point_name] = quadrants[1].points[point_name]
else:
quadrants[2].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]
self.occ.mirror([(0, quadrants[2].points[point_name])], 1, 0, 0, 0)
if point.y == 0. and symmetry == 'x':
self.symmetric_bnds['x_n']['pnts'].append([point_name, point.x])
if symmetry in ['y', 'none']:
if point.y == 0.:
quadrants[4].points[point_name] = quadrants[1].points[point_name]
else:
quadrants[4].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]
self.occ.mirror([(0, quadrants[4].points[point_name])], 0, 1, 0, 0)
if point.x == 0. and symmetry == 'y':
self.symmetric_bnds['y_n']['pnts'].append([point_name, point.y])
if symmetry == 'none':
if point.y == 0.:
quadrants[3].points[point_name] = quadrants[2].points[point_name]
elif point.x == 0.:
quadrants[3].points[point_name] = quadrants[4].points[point_name]
else:
quadrants[3].points[point_name] = self.occ.copy([(0, quadrants[2].points[point_name])])[0][1]
self.occ.mirror([(0, quadrants[3].points[point_name])], 0, 1, 0, 0)
mirror_x = [1, -1, -1, 1]
mirror_y = [1, 1, -1, -1]
symmetric_bnds_order = {'x': [], 'y': []}
sym_lines_tags = {'x_p': [], 'y_p': [], 'x_n': [], 'y_n': []}
for line_name, line in iron.hyper_lines.items():
identifier = next((k for k in self.inv_nc.keys() if re.match(f'^{k}', line_name[2:])), None)
if not geom_dict.get(identifier, False): continue
quadrants = getattr(self.md.geometries, self.inv_nc[identifier]).quadrants #re.sub(r'\d+', '', line_name[2:])
pt1 = iron.key_points[line.kp1]
pt2 = iron.key_points[line.kp2]
if line.type == 'line':
for quadrant, qq in quadrants.items():
if quadrant == 1:
qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
if pt1.y == 0. and pt2.y == 0. and 'x_p' in list_bnds:
self.symmetric_bnds['x_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)
sym_lines_tags['x_p'].append(qq.lines[line_name])
symmetric_bnds_order['x'].append(min(pt1.x, pt2.x))
elif pt1.x == 0. and pt2.x == 0. and 'y_p' in list_bnds:
self.symmetric_bnds['y_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)
sym_lines_tags['y_p'].append(qq.lines[line_name])
symmetric_bnds_order['y'].append(min(pt1.y, pt2.y))
elif quadrant == 2:
if pt1.x == 0. and pt2.x == 0.:
qq.lines[line_name] = quadrants[1].lines[line_name]
else:
qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
if pt1.y == 0. and pt2.y == 0. and 'x_n' in list_bnds:
self.symmetric_bnds['x_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)
sym_lines_tags['x_n'].append(qq.lines[line_name])
elif quadrant == 4:
if pt1.y == 0. and pt2.y == 0.:
qq.lines[line_name] = quadrants[1].lines[line_name]
else:
qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
if pt1.x == 0. and pt2.x == 0. and 'y_n' in list_bnds:
self.symmetric_bnds['y_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)
sym_lines_tags['y_n'].append(qq.lines[line_name])
else: # 3
if pt1.y == 0. and pt2.y == 0.:
qq.lines[line_name] = quadrants[2].lines[line_name]
elif pt1.x == 0. and pt2.x == 0.:
qq.lines[line_name] = quadrants[4].lines[line_name]
else:
qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
elif line.type == 'arc':
center = Func.arc_center_from_3_points([pt1.x, pt1.y],
[iron.key_points[line.kp3].x, iron.key_points[line.kp3].y],
[pt2.x, pt2.y])
new_point_name = 'kp' + line_name + '_center'
arc_coordinates1 = (pt1.x, pt1.y)
arc_coordinates2 = (pt2.x, pt2.y)
arc_coordinates3 = (iron.key_points[line.kp3].x, iron.key_points[line.kp3].y)
# This code addresses a meshing error in MQXA and MB_2COILS that occurs when an arc is defined on any of
# the axes. The issue arises because the function Func.arc_center_from_3_points does not return exactly
# zero but a value with a magnitude of approximately 10^-17 when the two points are placed on the axes.
# Consequently, when using the method self.occ.addCircleArc(), which only takes in three points without
# specifying a direction, a problem arises. The addCircleArc() function always creates the arc with the
# smallest angle. However, since center point can be slightly above or below the axis, the arc can
# inadvertently be drawn in the wrong quadrant, leading to an incorrect result.
# -----------------------
# Check that arcs with points on the x-axis are drawn in the first quadrant
if arc_coordinates3[1] > 0 and arc_coordinates2[1] == 0 and arc_coordinates1[1] == 0 and center[1] > 0:
quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], -center[1], 0)
# Check that arcs with points on the y-axis are drawn in the first quadrant
elif arc_coordinates3[0] > 0 and arc_coordinates2[0] == 0 and arc_coordinates1[0] == 0 and center[0] > 0:
quadrants[1].points[new_point_name] = self.occ.addPoint(-center[0], center[1], 0)
else:
quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], center[1], 0)
# -----------------------
# gmsh.model.setEntityName(0, gm.iron.quadrants[1].points[new_point_name], 'iron_' + new_point_name)
if symmetry in ['x', 'none']:
if center[0] == 0.:
quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]
else:
quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)
if symmetry in ['y', 'none']:
if center[1] == 0.:
quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]
else:
quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)
if symmetry == 'none':
if center[1] == 0.:
quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]
else:
quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]
self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)
for quadrant, qq in quadrants.items():
qq.lines[line_name] = self.occ.addCircleArc(
qq.points[line.kp1], qq.points[line.kp3], qq.points[line.kp2], center=False)
elif line.type == 'circle':
center = [(pt1.x + pt2.x) / 2, (pt1.y + pt2.y) / 2]
radius = (np.sqrt(np.square(pt1.x - center[0]) + np.square(pt1.y - center[1])) +
np.sqrt(np.square(pt2.x - center[0]) + np.square(pt2.y - center[1]))) / 2
for quadrant, qq in quadrants.items():
qq.lines[line_name] = self.occ.addCircle(
mirror_x[quadrant - 1] * center[0], mirror_y[quadrant - 1] * center[1], 0, radius)
qq.points['kp' + line_name] = len(qq.points) + 1
elif line.type == 'ellipticArc':
a, b = line.arg1, line.arg2
x1, y1 = pt1.x, pt1.y
x2, y2 = pt2.x, pt2.y
x3 = np.power(x1, 2.0)
y3 = np.power(y1, 2.0)
x4 = np.power(x2, 2.0)
y4 = np.power(y2, 2.0)
a2 = np.power(a, 2.0)
b2 = np.power(b, 2.0)
expression = -4.0 * a2 * b2 + a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 - 2.0 * b2 * x1 * x2 + b2 * x4
xc = x1 / 2.0 + x2 / 2.0 - a * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 -
2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (y1 - y2) / (2.0 * b)
yc = y1 / 2.0 + y2 / 2.0 + b * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3
- 2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (x1 - x2) / (2.0 * a)
center = self.occ.addPoint(xc, yc, 0, lc)
axis_point_a = self.occ.addPoint(xc + a, yc, 0, lc)
axis_point_b = self.occ.addPoint(xc, yc + b, 0, lc)
new_point_name = 'kp' + line_name + '_center'
new_axis_a_point_name = 'kp' + line_name + '_a'
new_axis_b_point_name = 'kp' + line_name + '_b'
quadrants[1].points[new_point_name] = center
quadrants[1].points[new_axis_a_point_name] = axis_point_a
quadrants[1].points[new_axis_b_point_name] = axis_point_b
if symmetry in ['x', 'none']:
if xc == 0.: # Least amount of possible points.
quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]
quadrants[2].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]
quadrants[2].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]
else:
quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
quadrants[2].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]
quadrants[2].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]
self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)
self.occ.mirror([(0, quadrants[2].points[new_axis_a_point_name])], 1, 0, 0, 0)
self.occ.mirror([(0, quadrants[2].points[new_axis_b_point_name])], 1, 0, 0, 0)
if symmetry in ['y', 'none']:
if yc == 0.:
quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]
quadrants[4].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]
quadrants[4].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]
else:
quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)
quadrants[4].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]
self.occ.mirror([(0, quadrants[4].points[new_axis_a_point_name])], 0, 1, 0, 0)
quadrants[4].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]
self.occ.mirror([(0, quadrants[4].points[new_axis_b_point_name])], 0, 1, 0, 0)
if symmetry == 'none':
if yc == 0.:
quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]
quadrants[3].points[new_axis_a_point_name] = quadrants[2].points[new_axis_a_point_name]
quadrants[3].points[new_axis_b_point_name] = quadrants[2].points[new_axis_b_point_name]
else:
quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]
self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)
quadrants[3].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_a_point_name])])[0][1]
self.occ.mirror([(0, quadrants[3].points[new_axis_a_point_name])], 0, 1, 0, 0)
quadrants[3].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_b_point_name])])[0][1]
self.occ.mirror([(0, quadrants[3].points[new_axis_b_point_name])], 0, 1, 0, 0)
for quadrant, qq in quadrants.items():
qq.lines[line_name] = self.occ.addEllipseArc(
qq.points[line.kp1], qq.points[new_point_name], qq.points[new_axis_a_point_name if a > b else new_axis_b_point_name],
qq.points[line.kp2])
else:
raise ValueError('Hyper line {} not supported'.format(line.type))
if symmetry != 'none':
quadrants = self.md.geometries.iron_yoke.quadrants
indexes = {'x_p': 1, 'y_p': 1, 'x_n': 1, 'y_n': 1}
self.md.geometries.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)
for sym in list_bnds:
if sym in ['x_p', 'y_p']:
quadrant = 1
elif sym == 'x_n':
quadrant = 2
else: # 'y_n'
quadrant = 4
sym_lines_tags[sym] = [x for _, x in sorted(zip(symmetric_bnds_order[sym[0]], sym_lines_tags[sym]))]
self.symmetric_bnds[sym]['pnts'].append(['center', 0])
self.symmetric_bnds[sym]['pnts'].sort(key=lambda x: x[1])
self.md.geometries.symmetric_boundaries.lines[sym + '_center'] = self.occ.addLine(
self.md.geometries.air_inf.points['center'], quadrants[quadrant].points[self.symmetric_bnds[sym]['pnts'][1][0]])
sym_lines_tags[sym].insert(0, self.md.geometries.symmetric_boundaries.lines[sym + '_center'])
for i, pnt in enumerate(self.symmetric_bnds[sym]['pnts'][1:-1]):
pnt_next = self.symmetric_bnds[sym]['pnts'][i + 2][0]
if not any(pnt[0] in s and pnt_next in s for s in self.symmetric_bnds[sym]['line_pnts']):
self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]] =\
self.occ.addLine(quadrants[quadrant].points[pnt[0]], quadrants[quadrant].points[pnt_next])
sym_lines_tags[sym].insert(indexes[sym], self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]])
indexes[sym] += 1
if symmetry == 'xy':
self.symmetric_loop_lines['x'] = sym_lines_tags['x_p']
sym_lines_tags['y_p'].reverse()
self.symmetric_loop_lines['y'] = sym_lines_tags['y_p']
elif symmetry == 'x':
sym_lines_tags['x_n'].reverse()
self.symmetric_loop_lines['x'] = sym_lines_tags['x_n'] + sym_lines_tags['x_p']
elif symmetry == 'y':
sym_lines_tags['y_p'].reverse()
self.symmetric_loop_lines['y'] = sym_lines_tags['y_p'] + sym_lines_tags['y_n']
# add all areas of each quadrant. Useful for brep curves and meshing
for key in geometry_setting.areas: # only consider areas that are implemented
quadrants = getattr(self.md.geometries, key).quadrants
for quadrant, qq in quadrants.items():
for area_name, area in iron.hyper_areas.items(): ## all areas
def _add_loop():
# prevent additional curveloop generation when Enforcing the TSA mapping on the collar
if (run_type == 'TH'
and self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping
and (area_name.startswith('arc') and not area_name.startswith('arch') or area_name.startswith('arp'))
): # need to disable the pole area too as it is linked to the same curve
qq.areas[area_name] = dM.Area() ## initialise area without loop
else:
qq.areas[area_name] = dM.Area(
loop=self.occ.addCurveLoop([qq.lines[line] for line in area.lines]))
if iron.hyper_areas[area_name].material not in getattr(self.md.domains.groups_entities, key) and \
iron.hyper_areas[area_name].material != 'BH_air': ## add the material to the keys
# for the collar region, it is possible to overwrite the material -> intercept it here
if key == 'collar' and (self.data.magnet.solve.collar.material != iron.hyper_areas[area_name].material) and self.data.magnet.solve.collar.material is not None:
logger.warning("Overwriting the collar material for area {} to {} ".format(area_name, self.data.magnet.solve.collar.material))
iron.hyper_areas[area_name].material = self.data.magnet.solve.collar.material
getattr(self.md.domains.groups_entities, key)[iron.hyper_areas[area_name].material] = []
identifier = next((k for k in geom_dict.keys() if re.match(f'^{k}', area_name[2:])),
None) # match key from geom_dict to the area name (see naming convention)
if key == self.inv_nc.get(identifier, None): # re.sub(r'\d+', '', area_name[2:]),
_add_loop() # adds arch to collar, because c is in the naming convention of the collar
elif area_name.startswith('arh') and key == 'iron_yoke': # if not previous but it is a hole, assume iron
_add_loop()
# define inner collar lines
def define_inner_collar():
"""
Defines the inner collar line used for the thermal TSA + for the A projection
"""
self.occ.synchronize()
# only works if the inner collar line is an arc -> just disable 'arc' and calc for all lines
# alternative method. Find all "arc" lines and then select the closest to the center
for quad, object in self.md.geometries.collar.quadrants.items():
arc_line_tags = [object.lines[name] for name in object.lines.keys() if
self.geom.iron.hyper_lines[name].type == 'arc']
closest_dist = 1000.
for tag in arc_line_tags:
x, y, _ = gmsh.model.getValue(1, tag, [0.5]) # pick one point on the arc
dist = np.sqrt(x ** 2 + y ** 2)
if dist < closest_dist:
closest_dist = dist
closest_line = tag ## assumes it is only one line per quadrant
self.md.geometries.collar.inner_boundary_tags[quad] = [closest_line]
def define_collar_cooling():
"""
Defines the cooling holes in the collar
"""
self.occ.synchronize()
line_names = [item for key in self.geom.iron.hyper_areas.keys() if 'ch' in key for item in self.geom.iron.hyper_areas[key].lines]
# line names are the same in each quadrant. Tags are unique
for quad, qq in self.md.geometries.collar.quadrants.items():
self.md.geometries.collar.cooling_tags.extend([qq.lines[line] for line in line_names])
# these tags are only used to be skipped for enfrocing_TSA_mapping
# we need the inner collar lines if we want to do TSA, so no need to define it
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA_new: define_inner_collar()
# we only need to specify the air holes if we want cooling OR if we enforce TSA nodes on the collar
if run_type == 'TH' and (self.data.magnet.solve.thermal.collar_cooling.enabled or self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping): define_collar_cooling()
def constructWedgeGeometry(self, use_TSA):
"""
Generates points, hyper lines, and curve loops for the wedges
"""
def _addMidLayerThinShellPoints(wedge_current):
def __addThinShellPoints(side_case, mid_layer_ts):
if side_case == 'outer':
mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +
Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2
else:
mean_rad_current = (Func.points_distance([wedge_current.iH.x, wedge_current.iH.y], wedge_center) +
Func.points_distance([wedge_current.iL.x, wedge_current.iL.y], wedge_center)) / 2
are_endpoints = {}
for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():
blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]
blk_next = wnd.blocks[blk_nr_next]
ht_list_next = (list(blk_next.half_turns.keys()) if blk_nr_next == list(wnd.blocks.keys())[0] else list(
reversed(blk_next.half_turns.keys())))
hh = blk_next.half_turns[ht_list_next[-1]].corners.bare
ll = blk_next.half_turns[ht_list_next[0]].corners.bare
bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
if side_case == 'outer':
block_list = self.md.geometries.coil.anticlockwise_order.coils[wedge.order_l.coil].layers[wedge.order_l.layer + 1]
blk_index = [blk.block for blk in block_list].index(blk_nr_next)
if blk_index + 1 == len(block_list): blk_index = -1
for blk in block_list[blk_index + 1:] + block_list[:blk_index + 1]:
if blk.winding == block_list[blk_index].winding:
ht_index = -1
break
elif blk.pole != block_list[blk_index].pole:
ht_index = 0
break
hh = blk_next.half_turns[ht_list_next[ht_index]].corners.bare
ll = blk_next.half_turns[ht_list_next[0 if ht_index == -1 else -1]].corners.bare
mean_rad_next = (Func.points_distance([hh.iH.x, hh.iH.y], wedge_center) +
Func.points_distance([ll.iL.x, ll.iL.y], wedge_center)) / 2
else:
mean_rad_next = (Func.points_distance([bc_next.oH.x, bc_next.oH.y], wedge_center) +
Func.points_distance([bc_next.oL.x, bc_next.oL.y], wedge_center)) / 2
mean_rad = (mean_rad_current + mean_rad_next) / 2
mid_layer = self.findMidLayerPoint(wedge_current.oH, bc_next.iH, wedge.corrected_center.outer, mean_rad)\
if side_case == 'outer' else self.findMidLayerPoint(wedge_current.iH, bc_next.oH, wedge.corrected_center.inner, mean_rad)
are_endpoints[wnd_nr] = self.getMidLayerEndpoints(wedge_current, bc_next, wedge_center, mid_layer_arc_pnt=mid_layer)
for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():
blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]
blk_next = wnd.blocks[blk_nr_next]
is_first_blk_next = blk_nr_next == list(wnd.blocks.keys())[0]
ht_list_next = (list(blk_next.half_turns.keys()) if is_first_blk_next else list(
reversed(blk_next.half_turns.keys())))
if are_endpoints[wnd_nr]: # this is empty if the wedge and the block are not radially adjacent
endpoints = are_endpoints[wnd_nr][0]
which_entity = are_endpoints[wnd_nr][1]
mid_layer_name = 'w' + str(wedge_nr) + '_' + str(blk_nr_next)
mid_layer_ts[mid_layer_name] = dM.Region()
ts_wdg = mid_layer_ts[mid_layer_name]
beg = ('w' + str(wedge_nr) if which_entity['lower'] == 'current' else str(ht_list_next[0])) + 'l'
ts_wdg.points[beg] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
ht_lower_angles = {}
for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):
for pnt1, pnt2, side in zip([[ht.corners.bare.iL.x, ht.corners.bare.iL.y], [ht.corners.bare.iH.x, ht.corners.bare.iH.y]],
[[ht.corners.bare.oL.x, ht.corners.bare.oL.y], [ht.corners.bare.oH.x, ht.corners.bare.oH.y]],
['l', 'h']):
line_pars_current = Func.line_through_two_points(pnt1, pnt2)
intersect_prev = Func.intersection_between_arc_and_line(
line_pars_current, [wedge_center, endpoints['higher'], endpoints['lower']])
if intersect_prev:
ts_wdg.points[str(ht_nr) + side] = self.occ.addPoint(intersect_prev[0][0], intersect_prev[0][1], 0)
elif side == 'l':
intrsc = Func.intersection_between_circle_and_line(line_pars_current, [wedge_center, endpoints['lower']], get_only_closest=True)[0]
ht_lower_angles[ht_nr] = Func.arc_angle_between_point_and_abscissa([intrsc[0], intrsc[1]], wedge_center)
end = ('w' + str(wedge_nr) if which_entity['higher'] == 'current' else str(ht_list_next[-1])) + 'h'
if all('w' in pnt_name for pnt_name in list(ts_wdg.points.keys())): # only one thin-shell 'within' the facing half-turn
wdg_angle_il = Func.arc_angle_between_point_and_abscissa([endpoints['lower'][0], endpoints['lower'][1]], wedge_center)
for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):
if ht_lower_angles[ht_nr] > wdg_angle_il: break
prev_nr = str(ht_nr)
end = prev_nr + 'h'
ts_wdg.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
# Create auxiliary thin shells for outliers
# if both corners belong to thin shells, continue
used_wdg_corners = [False, False]
for ep in are_endpoints.values():
if ep is not None:
if ep[1]['higher'] == 'current': used_wdg_corners[1] = True
if ep[1]['lower'] == 'current': used_wdg_corners[0] = True
if side_case == 'inner':
for ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.keys():
if ts_name[ts_name.index('_') + 1:] == 'w' + str(wedge_nr):
for ep_key, ep in are_endpoints_wdg[int(ts_name[1:ts_name.index('_')])].items():
if ep is not None:
if ep[1]['higher'] == 'next': used_wdg_corners[1] = True
if ep[1]['lower'] == 'next': used_wdg_corners[0] = True
else:
if wedge_nr in are_endpoints_wdg:
for ep in are_endpoints_wdg[wedge_nr].values():
if ep is not None:
if ep[1]['higher'] == 'current': used_wdg_corners[1] = True
if ep[1]['lower'] == 'current': used_wdg_corners[0] = True
if not used_wdg_corners[1]:
for wdg_nr, wdg in self.geom.wedges.items():
if blk_nr_next == wdg.order_l.block: used_wdg_corners[1] = True
if not used_wdg_corners[0]:
for wdg_nr, wdg in self.geom.wedges.items():
if blk_nr_next == wdg.order_h.block: used_wdg_corners[0] = True
if not all(used_wdg_corners):
def ___create_aux_mid_layer_point(ss, points):
mid_layer_ts_aux[mid_layer_name] = dM.Region()
circle_pnt = [endpoints[ss][0], endpoints[ss][1]]
inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),
[[wedge.corrected_center.outer.x, wedge.corrected_center.outer.y], circle_pnt], get_only_closest=True)[0]
mid_layer_ts_aux[mid_layer_name].points[str(wedge_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)
mid_layer_ts_aux[mid_layer_name].points['center'] = self.occ.addPoint(wedge_data[wedge_nr][1].x, wedge_data[wedge_nr][1].y, 0)
mid_layer_ts_aux[mid_layer_name].lines['w' + str(wedge_nr)] = 0
if which_entity['higher'] == 'current' and which_entity['lower'] != 'current':
___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],
[wedge_current.oL.x, wedge_current.oL.y]])
elif which_entity['higher'] != 'current' and which_entity['lower'] == 'current':
___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],
[wedge_current.oH.x, wedge_current.oH.y]])
else: # whole block 'within' the facing wedge
for wdg_nr, wdg in self.geom.wedges.items():
if blk_nr_next == wdg.order_h.block:
___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],
[wedge_current.oH.x, wedge_current.oH.y]])
break
elif blk_nr_next == wdg.order_l.block:
___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],
[wedge_current.oL.x, wedge_current.oL.y]])
break
pole = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole]
blk_list_current = list(pole.layers[wedge.order_l.layer].windings[wedge.order_l.winding].blocks.keys())
if wedge.order_l.layer < len(pole.layers):
__addThinShellPoints('outer', self.md.geometries.thin_shells.mid_layers_wdg_to_ht)
if wedge.order_l.layer > 1:
__addThinShellPoints('inner', self.md.geometries.thin_shells.mid_layers_ht_to_wdg)
wedges = self.md.geometries.wedges
mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
wedge_data = {}
wdgs_corners = {}
for wedge_nr, wedge in self.geom.wedges.items():
wdgs_corners[wedge_nr] = {}
corners = wdgs_corners[wedge_nr]
if wedge.order_l.coil not in wedges.coils:
wedges.coils[wedge.order_l.coil] = dM.WedgeLayer()
if wedge.order_l.layer not in wedges.coils[wedge.order_l.coil].layers:
wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer] = dM.WedgeRegion()
wedge_layer = wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer]
wedge_layer.wedges[wedge_nr] = dM.Region()
wedge_reg = wedge_layer.wedges[wedge_nr]
wedge_layer.block_prev[wedge_nr] = wedge.order_l.block
wedge_layer.block_next[wedge_nr] = wedge.order_h.block
wnd = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[
wedge.order_l.layer].windings[wedge.order_l.winding]
wnd_next = self.geom.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[
wedge.order_h.layer].windings[wedge.order_h.winding]
block = wnd.blocks[wedge.order_l.block]
block_next = wnd_next.blocks[wedge.order_h.block]
corners['last_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[
wedge.order_l.layer].windings[wedge.order_l.winding].blocks[wedge.order_l.block].half_turns.areas.keys())[-1])
corners['first_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[
wedge.order_h.layer].windings[wedge.order_h.winding].blocks[wedge.order_h.block].half_turns.areas.keys())[0])
ht_current = block.half_turns[corners['last_ht']].corners.bare
ht_next = block_next.half_turns[corners['first_ht']].corners.bare
d_current = self.data.conductors[wnd.conductor_name].cable.th_insulation_along_width * 2
d_next = self.data.conductors[wnd_next.conductor_name].cable.th_insulation_along_width * 2
for pnt_close, pnt_far, wdg_corner, d in zip([ht_current.iH, ht_current.oH, ht_next.iL, ht_next.oL],
[ht_current.iL, ht_current.oL, ht_next.iH, ht_next.oH],
['il', 'ol', 'ih', 'oh'], [d_current, d_current, d_next, d_next]):
if abs(pnt_far.x - pnt_close.x) > 0.:
m = (pnt_far.y - pnt_close.y) / (pnt_far.x - pnt_close.x)
b = pnt_close.y - m * pnt_close.x
root = np.sqrt(- pnt_close.x ** 2 * m ** 2 - 2 * pnt_close.x * b * m + 2 * pnt_close.x * pnt_close.y * m
- b ** 2 + 2 * b * pnt_close.y - pnt_close.y ** 2 + d ** 2 * m ** 2 + d ** 2)
pnt1_x = (pnt_close.x - b * m + pnt_close.y * m + root) / (m ** 2 + 1)
pnt1_y = m * pnt1_x + b
pnt2_x = (pnt_close.x - b * m + pnt_close.y * m - root) / (m ** 2 + 1)
pnt2_y = m * pnt2_x + b
corners[wdg_corner] = Coord(x=pnt1_x, y=pnt1_y) if Func.points_distance([pnt1_x, pnt1_y], [pnt_far.x, pnt_far.y]) >\
Func.points_distance([pnt_close.x, pnt_close.y], [pnt_far.x, pnt_far.y]) else Coord(x=pnt2_x, y=pnt2_y)
else:
bore_cnt_x = self.geom.coil.coils[wedge.order_l.coil].bore_center.x
pnt1_y, pnt2_y = pnt_close.y + d, pnt_close.y - d
corners[wdg_corner] = Coord(x=pnt_close.x,
y=pnt1_y if (wdg_corner[-1] == 'l' and pnt_close.x > bore_cnt_x) or
(wdg_corner[-1] == 'h' and pnt_close.x < bore_cnt_x) else pnt2_y)
wedge_reg.points[wdg_corner] = self.occ.addPoint(corners[wdg_corner].x, corners[wdg_corner].y, 0)
inner = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],
[corners['ih'].x, corners['ih'].y], [corners['il'].x, corners['il'].y])
outer = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],
[corners['oh'].x, corners['oh'].y], [corners['ol'].x, corners['ol'].y])
wedge_data[wedge_nr] = [Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']), wedge.corrected_center.outer]
wedge_reg.points['inner_center'] = self.occ.addPoint(inner[0], inner[1], 0)
wedge_reg.points['outer_center'] = self.occ.addPoint(outer[0], outer[1], 0)
wedge_reg.lines['h'] = self.occ.addLine(wedge_reg.points['ih'], wedge_reg.points['oh'])
wedge_reg.lines['l'] = self.occ.addLine(wedge_reg.points['il'], wedge_reg.points['ol'])
wedge_reg.lines['i'] = self.occ.addCircleArc(wedge_reg.points['ih'], wedge_reg.points['inner_center'], wedge_reg.points['il'])
wedge_reg.lines['o'] = self.occ.addCircleArc(wedge_reg.points['oh'], wedge_reg.points['outer_center'], wedge_reg.points['ol'])
"""
logger.warning("Using straight wedge geometry") # required for the projection
wedge_reg.lines['i'] = self.occ.addLine(wedge_reg.points['ih'], wedge_reg.points['il'])
wedge_reg.lines['o'] = self.occ.addLine(wedge_reg.points['oh'], wedge_reg.points['ol'])
"""
wedge_reg.areas[str(wedge_nr)] = dM.Area(loop=self.occ.addCurveLoop(
[wedge_reg.lines['i'], wedge_reg.lines['l'], wedge_reg.lines['o'], wedge_reg.lines['h']]))
if use_TSA:
# Wedge thin shells
mid_layer_ts = self.md.geometries.thin_shells.mid_layers_wdg_to_wdg
are_endpoints_wdg = {}
for coil_nr, coil in self.md.geometries.wedges.coils.items():
layer_list = list(coil.layers.keys())
for layer_nr, layer in coil.layers.items():
if layer_list.index(layer_nr) + 1 < len(layer_list):
for wedge_nr, wedge in layer.wedges.items():
are_endpoints_wdg[wedge_nr] = {}
are_endpoints = are_endpoints_wdg[wedge_nr]
wedge_current = wedge_data[wedge_nr][0]
wedge_center = [wedge_data[wedge_nr][1].x, wedge_data[wedge_nr][1].y]
mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +
Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2
for wdg_next_nr, wdg_next in coil.layers[layer_nr + 1].wedges.items():
if self.geom.wedges[wedge_nr].order_l.pole == self.geom.wedges[wdg_next_nr].order_l.pole:
wedge_next = wedge_data[wdg_next_nr][0]
mean_rad_next = (Func.points_distance([wedge_next.iH.x, wedge_next.iH.y], wedge_center) +
Func.points_distance([wedge_next.iL.x, wedge_next.iL.y], wedge_center)) / 2
mean_rad = (mean_rad_current + mean_rad_next) / 2
mid_layer = self.findMidLayerPoint(wedge_current.oH, wedge_next.iH, wedge_data[wedge_nr][1], mean_rad)
are_endpoints[wdg_next_nr] = self.getMidLayerEndpoints(wedge_current, wedge_next, wedge_center, mid_layer_arc_pnt=mid_layer)
if are_endpoints[wdg_next_nr]: # this is empty if the wedges are not radially adjacent
endpoints = are_endpoints[wdg_next_nr][0]
mid_layer_name = 'w' + str(wedge_nr) + '_w' + str(wdg_next_nr)
mid_layer_ts[mid_layer_name] = dM.Region()
ts = mid_layer_ts[mid_layer_name]
ts.points['center'] = self.occ.addPoint(wedge_center[0], wedge_center[1], 0)
ts.points['beg'] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
end = 'w' + str(wedge_nr if are_endpoints[wdg_next_nr][1] == 'current' else wdg_next_nr)
ts.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
# Half-turn thin shells
for wedge_nr, wedge in self.geom.wedges.items():
corners = wdgs_corners[wedge_nr]
# Mid layer lines
wedge_center = [self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y]
_addMidLayerThinShellPoints(Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']))
# Mid wedge-turn lines
mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn
for adj_blk, ht, inner, outer in zip([wedge.order_l, wedge.order_h], [corners['last_ht'], corners['first_ht']],
[corners['il'], corners['ih']], [corners['ol'], corners['oh']]):
mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)] = dM.Region()
ts = mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)]
ht_corners = self.geom.coil.coils[adj_blk.coil].poles[adj_blk.pole].layers[
adj_blk.layer].windings[adj_blk.winding].blocks[adj_blk.block].half_turns[ht].corners.bare
ht_corners_i = ht_corners.iH if ht == corners['last_ht'] else ht_corners.iL
ht_corners_o = ht_corners.oH if ht == corners['last_ht'] else ht_corners.oL
mid_inner = [(inner.x + ht_corners_i.x) / 2, (inner.y + ht_corners_i.y) / 2]
mid_outer = [(outer.x + ht_corners_o.x) / 2, (outer.y + ht_corners_o.y) / 2]
line_name = 'w' + str(wedge_nr) + '_' + str(ht)
ts.points[line_name + '_i'] = self.occ.addPoint(mid_inner[0], mid_inner[1], 0)
ts.points[line_name + '_o'] = self.occ.addPoint(mid_outer[0], mid_outer[1], 0)
def constructCoilGeometry(self, run_type):
"""
Generates points, hyper lines, and curve loops for the coil half-turns
"""
symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'
# Sub domains angles: first key means 'from 0 to x'; second key means 'from x to 2*pi'
if symmetry == 'xy':
angle_range = {'to': np.pi / 2, 'from': 2 * np.pi}
elif symmetry == 'x':
angle_range = {'to': np.pi, 'from': 2 * np.pi}
elif symmetry == 'y':
angle_range = {'to': np.pi / 2, 'from': 3 / 2 * np.pi}
elif symmetry == 'none':
angle_range = {'to': 2 * np.pi, 'from': 0}
else:
raise Exception('Symmetry plane not supported.')
def _addMidLayerThinShellPoints(pnt_params, ss, name, case):
endpnts, cnt = ts_endpoints[name]
if len(pnt_params) == 3: # line parameters (cos-theta Rutherford)
intersect[name] = Func.intersection_between_arc_and_line(pnt_params, [cnt, endpnts['higher'], endpnts['lower']])
if intersect[name]:
intersect[name] = intersect[name][0]
pnt_angle = Func.arc_angle_between_point_and_abscissa(intersect[name], cnt)
elif len(pnt_params) == 4: # points coordinates (cos-theta Mono)
wnd_next = list(pole.layers[layer_nr + (1 if case == 'current' else -1)].windings.keys())[
list(pole.layers[layer_nr].windings.keys()).index(winding_nr)]
blk_next = pole.layers[layer_nr + (1 if case == 'current' else -1)].windings[wnd_next].blocks[
int(ts_name[ts_name.index('_') + 1:] if case == 'current' else ts_name[:ts_name.index('_')])]
ht_next = blk_next.half_turns[list(blk_next.half_turns.keys() if is_first_blk else reversed(blk_next.half_turns.keys()))[ht_list.index(halfTurn_nr)]].corners.bare
coord_next = (ht_next.iL if ss == 'l' else ht_next.iH) if case == 'current' else (ht_next.oL if ss == 'l' else ht_next.oH)
pnt = [(pnt_params[2 if case == 'current' else 0] + coord_next.x) / 2, (pnt_params[3 if case == 'current' else 1] + coord_next.y) / 2]
pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)
pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)
pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)
intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None
else: # point coordinates (block-coil)
pnt = [endpnts['higher'][0], pnt_params[1]] if coil.type == 'common-block-coil' else [pnt_params[0], endpnts['higher'][1]]
if abs(endpnts['higher'][1]) > 1e-6:
pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)
pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)
pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)
else:
pnt_angle = abs(pnt_params[0])
pnt_angle_h = abs(endpnts['higher'][0])
pnt_angle_l = abs(endpnts['lower'][0])
intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None
if intersect[name]:
mid_layer_ts[name].mid_layers.points[str(halfTurn_nr) + ss] = \
self.occ.addPoint(intersect[name][0], intersect[name][1], 0)
mid_layer_ts[name].point_angles[str(halfTurn_nr) + ss] = Func.sig_dig(pnt_angle)
if len(pnt_params) == 2 and not intersect[name] and (abs(pnt_angle - pnt_angle_h) < 1e-6 or abs(pnt_angle - pnt_angle_l) < 1e-6):
intersect[name] = pnt
return intersect
def _addMidLayerThinShellGroup(cl, for_mid_pole=False, mid_coil=False):
is_first_blk_next = block_nr_next == list(winding_next.blocks.keys())[0]
if 'solenoid' in cl.type:
ht_list_next = list(reversed(block_next.half_turns.keys()) if layer_nr % 2 == 0 else list(block_next.half_turns.keys()))
elif cl.type == 'reversed-block-coil':
ht_list_next = (list(block_next.half_turns.keys()) if not is_first_blk_next else list(reversed(block_next.half_turns.keys())))
else:
ht_list_next = (list(block_next.half_turns.keys()) if is_first_blk_next else list(reversed(block_next.half_turns.keys())))
hh = block_next.half_turns[ht_list_next[-1]].corners.bare
ll = block_next.half_turns[ht_list_next[0]].corners.bare
bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
if 'block-coil' in cl.type or (cable_type_curr in ['Mono', 'Ribbon'] and not mid_coil):
center = [cl.bore_center.x, cl.bore_center.y]
are_endpoints = self.getMidLayerEndpoints(bc_current, bc_next, center, coil_type=cl.type, cable_type=cable_type_curr, is_for_mid_pole=for_mid_pole)
else:
mean_rad_next = (Func.points_distance([bc_next.iH.x, bc_next.iH.y], [cl.bore_center.x, cl.bore_center.y]) +
Func.points_distance([bc_next.iL.x, bc_next.iL.y], [cl.bore_center.x, cl.bore_center.y])) / 2
mean_rad = (mean_rad_current + mean_rad_next) / 2
mid_layer_h = self.findMidLayerPoint(bc_current.oH, bc_next.iH, cl.bore_center, mean_rad)
mid_layer_l = self.findMidLayerPoint(bc_current.oL, bc_next.iL, cl.bore_center, mean_rad)
mid_ht_next_i = int(len(ht_list_next) / 2) if len(ht_list_next) % 2 == 0 else round(len(ht_list_next) / 2)
mid_ht_next = block_next.half_turns[ht_list_next[mid_ht_next_i - 1]].corners.insulated
mid_layer_m = self.findMidLayerPoint(mid_ht_current.oH, mid_ht_next.iH, cl.bore_center, mean_rad)
center = Func.arc_center_from_3_points(mid_layer_h, mid_layer_m, mid_layer_l)
are_endpoints = self.getMidLayerEndpoints(bc_current, bc_next, center, mid_layer_arc_pnt=mid_layer_h, cable_type=cable_type_curr)
if are_endpoints: # this is empty if the blocks are not radially adjacent
endpoints = are_endpoints[0]
which_block = are_endpoints[1]
mid_layer_name = blk_nr + '_' + str(block_nr_next)
if for_mid_pole:
block_coil_mid_pole_next_blks_list[block_nr_next].append(mid_layer_name)
block_coil_ts_endpoints[mid_layer_name] = [endpoints, center]
else:
if block_nr_next not in list(next_blks_list.keys()):
next_blks_list[block_nr_next] = []
next_blks_list[block_nr_next].append(mid_layer_name)
ts_endpoints[mid_layer_name] = [endpoints, center]
mid_layer_ts[mid_layer_name] = dM.MidLayer()
mid_layer_ts[mid_layer_name].half_turn_lists[blk_nr] = ht_list
mid_layer_ts[mid_layer_name].half_turn_lists[str(block_nr_next)] = ht_list_next
beg = (str(ht_list[0]) if which_block['lower'] == 'current' else str(ht_list_next[0])) + 'l'
mid_layer_ts[mid_layer_name].mid_layers.points[beg] = \
self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
end = (str(ht_list[-1]) if which_block['higher'] == 'current' else str(ht_list_next[-1])) + 'h'
mid_layer_ts[mid_layer_name].mid_layers.points[end] = \
self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
if not for_mid_pole or (for_mid_pole and abs(endpoints['higher'][1]) > 1e-6):
mid_layer_ts[mid_layer_name].point_angles[beg] =\
Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['lower'], center))
mid_layer_ts[mid_layer_name].point_angles[end] =\
Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['higher'], center))
else:
mid_layer_ts[mid_layer_name].point_angles[beg] = abs(endpoints['lower'][0])
mid_layer_ts[mid_layer_name].point_angles[end] = abs(endpoints['higher'][0])
# Create anticlockwise order of blocks
present_blocks = []
block_corner_angles = {}
concentric_coils = self.md.geometries.coil.concentric_coils
acw_order = self.md.geometries.coil.anticlockwise_order.coils
self.md.geometries.coil.physical_order = self.geom.coil.physical_order
for coil_nr, coil in self.geom.coil.coils.items():
# if coil_nr not in block_corner_angles:
block_corner_angles[coil_nr] = {}
if (coil.bore_center.x, coil.bore_center.y) not in concentric_coils:
concentric_coils[(coil.bore_center.x, coil.bore_center.y)] = []
concentric_coils[(coil.bore_center.x, coil.bore_center.y)].append(coil_nr)
for pole_nr, pole in coil.poles.items():
for layer_nr, layer in pole.layers.items():
if layer_nr not in block_corner_angles[coil_nr]:
block_corner_angles[coil_nr][layer_nr] = {}
blk_angles = block_corner_angles[coil_nr][layer_nr]
for winding_nr, winding in layer.windings.items():
for block_nr, block in winding.blocks.items():
blk_angles[block_nr] = {'angle': Func.sig_dig(Func.arc_angle_between_point_and_abscissa(
[block.block_corners.iL.x, block.block_corners.iL.y],
[coil.bore_center.x, coil.bore_center.y])), 'keys': [pole_nr, winding_nr]}
higher_angle = Func.sig_dig(Func.arc_angle_between_point_and_abscissa(
[block.block_corners.iH.x, block.block_corners.iH.y],
[coil.bore_center.x, coil.bore_center.y]))
if ((blk_angles[block_nr]['angle'] <= angle_range['to'] and higher_angle <= angle_range['to']) or
(angle_range['from'] <= blk_angles[block_nr]['angle'] and angle_range['from'] <= higher_angle)):
present_blocks.append(block_nr)
for coil_nr, coil in block_corner_angles.items():
acw_order[coil_nr] = dM.LayerOrder()
for layer_nr, layer in coil.items():
acw_order[coil_nr].layers[layer_nr] = []
ordered_blocks = [[block_nr, block['angle'], block['keys']] for block_nr, block in layer.items()]
ordered_blocks.sort(key=lambda x: x[1])
for blk in ordered_blocks:
if blk[0] in present_blocks:
acw_order[coil_nr].layers[layer_nr].append(dM.AnticlockwiseOrder(pole=blk[2][0], winding=blk[2][1], block=blk[0]))
# Check if there are concentric coils
for bore_center, coils in concentric_coils.items():
if len(coils) > 1:
radii = []
for coil_nr in coils:
lyr = self.geom.coil.coils[coil_nr].poles[1].layers[1]
blk = list(lyr.windings.keys())[0]
radii.append([coil_nr, Func.points_distance(bore_center, [lyr.windings[blk].blocks[blk].block_corners.iL.x, lyr.windings[blk].blocks[blk].block_corners.iL.y])])
radii.sort(key=lambda x: x[1])
concentric_coils[bore_center] = [rad[0] for rad in radii]
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht
# Collect block couples for block-coil mid-pole thin shells
block_coil_mid_pole_next_blks_list = {}
block_coil_ts_endpoints = {}
for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
if self.geom.coil.coils[coil_nr].type in ['block-coil', 'reversed-block-coil']:
self.block_coil_mid_pole_blks[coil_nr] = []
first_lyr = list(coil.layers.keys())[0]
layer = coil.layers[first_lyr]
for nr, block_order in enumerate(layer):
blk_next_index = nr + 1 if nr + 1 < len(layer) else 0
if layer[blk_next_index].pole != block_order.pole:
self.block_coil_mid_pole_blks[coil_nr].append([block_order, layer[blk_next_index]])
block_coil_mid_pole_next_blks_list[layer[blk_next_index].block] = []
# Mid pole lines for block-coils
for coil_nr, coil in self.block_coil_mid_pole_blks.items():
coil_geom = self.geom.coil.coils[coil_nr]
for mid_pole in coil:
winding = self.geom.coil.coils[coil_nr].poles[mid_pole[0].pole].layers[1].windings[mid_pole[0].winding]
cable_type_curr = self.data.conductors[winding.conductor_name].cable.type
block_nr = mid_pole[0].block
blk_nr = str(block_nr)
block = winding.blocks[block_nr]
is_first_blk = block_nr == list(winding.blocks.keys())[0]
if coil_geom.type == 'reversed-block-coil':
ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))
else:
ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))
hh = block.half_turns[ht_list[-1]].corners.bare
ll = block.half_turns[ht_list[0]].corners.bare
bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
winding_next = self.geom.coil.coils[coil_nr].poles[mid_pole[1].pole].layers[1].windings[mid_pole[1].winding]
block_nr_next = mid_pole[1].block
block_next = winding_next.blocks[block_nr_next]
_addMidLayerThinShellGroup(coil_geom, for_mid_pole=True)
mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
self.md.geometries.coil.physical_order = self.geom.coil.physical_order
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
next_blks_list = block_coil_mid_pole_next_blks_list.copy()
ts_endpoints = block_coil_ts_endpoints.copy()
for coil_nr, coil in self.geom.coil.coils.items():
self.md.geometries.coil.coils[coil_nr] = dM.Pole()
coils = self.md.geometries.coil.coils[coil_nr]
coils.type = coil.type
coils.bore_center = coil.bore_center
for pole_nr, pole in coil.poles.items():
coils.poles[pole_nr] = dM.Layer()
poles = coils.poles[pole_nr]
for layer_nr, layer in pole.layers.items():
poles.layers[layer_nr] = dM.Winding()
layers = poles.layers[layer_nr]
for winding_nr, winding in layer.windings.items():
cable_type_curr = self.data.conductors[winding.conductor_name].cable.type
layers.windings[winding_nr] = dM.Block(conductor_name=winding.conductor_name, conductors_number=winding.conductors_number)
windings = layers.windings[winding_nr]
blk_list_current = list(winding.blocks.keys())
for block_nr, block in winding.blocks.items():
if block_nr in present_blocks:
blk_nr = str(block_nr)
windings.blocks[block_nr] = dM.BlockData(current_sign=block.current_sign)
hts = windings.blocks[block_nr].half_turns
is_first_blk = block_nr == list(winding.blocks.keys())[0]
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
if 'solenoid' in coil.type:
ht_list = (list(reversed(block.half_turns.keys()) if (layer_nr - 1) % 2 == 0 else list(block.half_turns.keys())))
elif coil.type == 'reversed-block-coil':
ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))
else:
ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))
hh = block.half_turns[ht_list[-1]].corners.bare
ll = block.half_turns[ht_list[0]].corners.bare
bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
# Mid layer lines
mean_rad_current = (Func.points_distance([bc_current.oH.x, bc_current.oH.y], [coil.bore_center.x, coil.bore_center.y]) +
Func.points_distance([bc_current.oL.x, bc_current.oL.y], [coil.bore_center.x, coil.bore_center.y])) / 2
mid_ht_current_i = int(len(ht_list) / 2) if len(ht_list) % 2 == 0 else round(len(ht_list) / 2)
mid_ht_current = block.half_turns[ht_list[mid_ht_current_i - 1]].corners.insulated
concentric_coil = concentric_coils[(coil.bore_center.x, coil.bore_center.y)]
if layer_nr < len(pole.layers):
for winding_nr_next, winding_next in pole.layers[layer_nr + 1].windings.items():
if cable_type_curr == 'Rutherford' or\
(cable_type_curr in ['Mono', 'Ribbon'] and
list(pole.layers[layer_nr + 1].windings.keys()).index(winding_nr_next) == list(layer.windings.keys()).index(winding_nr)):
blk_list_next = list(winding_next.blocks.keys())
block_nr_next = blk_list_next[blk_list_current.index(block_nr)]
block_next = winding_next.blocks[block_nr_next]
_addMidLayerThinShellGroup(coil)
elif concentric_coil.index(coil_nr) + 1 < len(concentric_coil):
coil_nr_next = concentric_coil[concentric_coil.index(coil_nr) + 1]
for pole_nr_next, pole_next in self.geom.coil.coils[coil_nr_next].poles.items():
for layer_nr_next, layer_next in pole_next.layers.items():
if layer_nr_next == 1:
for winding_nr_next, winding_next in layer_next.windings.items():
for block_nr_next, block_next in winding_next.blocks.items():
_addMidLayerThinShellGroup(coil, mid_coil=True)
else:
blk_ins = windings.blocks[block_nr].insulation
blk_ins.areas[blk_nr] = dM.Area()
if 'solenoid' in coil.type:
ht_items = (list(reversed(block.half_turns.items()) if layer_nr - 1 % 2 == 0 else list(block.half_turns.items())))
elif coil.type == 'reversed-block-coil':
ht_items = (block.half_turns.items() if not is_first_blk else reversed(block.half_turns.items()))
else:
ht_items = (block.half_turns.items() if is_first_blk else reversed(block.half_turns.items()))
for halfTurn_nr, halfTurn in ht_items:
ht_nr = str(halfTurn_nr)
ht = halfTurn.corners.insulated
hts.areas[ht_nr] = dM.Area()
ht_b = halfTurn.corners.bare
hts.points[ht_nr + 'ih'] = self.occ.addPoint(ht_b.iH.x, ht_b.iH.y, 0)
hts.points[ht_nr + 'il'] = self.occ.addPoint(ht_b.iL.x, ht_b.iL.y, 0)
hts.points[ht_nr + 'oh'] = self.occ.addPoint(ht_b.oH.x, ht_b.oH.y, 0)
hts.points[ht_nr + 'ol'] = self.occ.addPoint(ht_b.oL.x, ht_b.oL.y, 0)
hts.lines[ht_nr + 'i'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'il'])
hts.lines[ht_nr + 'o'] = self.occ.addLine(hts.points[ht_nr + 'oh'], hts.points[ht_nr + 'ol'])
hts.lines[ht_nr + 'l'] = self.occ.addLine(hts.points[ht_nr + 'il'], hts.points[ht_nr + 'ol'])
hts.lines[ht_nr + 'h'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'oh'])
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
intersection = {}
# Create mid layer points and compute their angle to the x-axis
for mid_lyr_type in ['current', 'previous']:
for pnt1, pnt2, side in zip(
[[ht_b.iH.x, ht_b.iH.y], [ht_b.iL.x, ht_b.iL.y]],
[[ht_b.oH.x, ht_b.oH.y], [ht_b.oL.x, ht_b.oL.y]], ['h', 'l']):
if (cable_type_curr in ['Mono', 'Ribbon'] and coil.type == 'cos-theta' and
(layer_nr < len(pole.layers) and mid_lyr_type == 'current' or layer_nr > 1 and mid_lyr_type == 'previous')):
pnts_input = pnt1 + pnt2
elif coil.type == 'cos-theta' and (cable_type_curr == 'Rutherford' or cable_type_curr in ['Mono', 'Ribbon'] and\
(layer_nr == len(pole.layers) and mid_lyr_type == 'current' or layer_nr == 1 and mid_lyr_type == 'previous')):
pnts_input = Func.line_through_two_points(pnt1, pnt2)
elif 'block-coil' in coil.type:
pnts_input = pnt1
intersect = {}
if mid_lyr_type == 'current':
# Current mid-layer
for ts_name in ts_endpoints.keys():
if blk_nr == ts_name[:ts_name.index('_')]:
_addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)
elif mid_lyr_type == 'previous':
# Previous mid-layer
if block_nr in next_blks_list:
for ts_name in next_blks_list[block_nr]:
_addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)
for key, value in intersect.items():
if key in intersection:
intersection[key][side] = value
else:
intersection[key] = {side: value}
# Search for half turns that face thin shells only partially
def __create_aux_mid_layer_point(ss, points):
mid_layer_ts_aux[key] = dM.Region()
if 'block-coil' in coil.type:
inter_pnt = [points[0], ts_endpoints[key][0][ss][1]]
else:
inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),
[ts_endpoints[key][1], ts_endpoints[key][0][ss]], get_only_closest=True)[0]
mid_layer_ts_aux[key].points[str(halfTurn_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)
mid_layer_ts_aux[key].lines[blk_nr] = 0
for key, value in intersection.items():
first_blk, second_blk = key.split('_')
if 'block-coil' in coil.type: #any(int(second_blk) == blk_order.block for blk_order in acw_order[coil_nr].layers[layer_nr]): # block-coil mid-pole case
if value['h'] and not value['l']:
__create_aux_mid_layer_point('lower', [ht_b.iL.x, ht_b.iL.y])
elif value['l'] and not value['h']:
__create_aux_mid_layer_point('higher', [ht_b.iH.x, ht_b.iH.y])
else:
relevant_blk = int(first_blk) if second_blk == blk_nr else int(second_blk)
if layer_nr == len(pole.layers) and blk_nr == first_blk:
lyr_blks = acw_order[coil_nr + 1].layers[1]
elif layer_nr == 1 and blk_nr == second_blk:
lyr_blks = acw_order[coil_nr - 1].layers[len(acw_order[coil_nr - 1].layers)]
else:
lyr_blks = acw_order[coil_nr].layers[layer_nr + (1 if first_blk == blk_nr else -1)]
for nr, block_order in enumerate(lyr_blks):
if block_order.block == relevant_blk:
block_order_curr = block_order
block_order_prev = lyr_blks[-1] if nr == 0 else lyr_blks[nr - 1]
block_order_next = lyr_blks[0] if nr + 1 == len(lyr_blks) else lyr_blks[nr + 1]
break
if value['h'] and not value['l'] and block_order_curr.winding == block_order_prev.winding:
__create_aux_mid_layer_point('lower', [[ht_b.iL.x, ht_b.iL.y], [ht_b.oL.x, ht_b.oL.y]])
elif value['l'] and not value['h'] and block_order_curr.winding == block_order_next.winding:
__create_aux_mid_layer_point('higher', [[ht_b.iH.x, ht_b.iH.y], [ht_b.oH.x, ht_b.oH.y]])
else:
blk_ins.points[ht_nr + 'ih'] = self.occ.addPoint(ht.iH.x, ht.iH.y, 0)
blk_ins.points[ht_nr + 'il'] = self.occ.addPoint(ht.iL.x, ht.iL.y, 0)
blk_ins.points[ht_nr + 'oh'] = self.occ.addPoint(ht.oH.x, ht.oH.y, 0)
blk_ins.points[ht_nr + 'ol'] = self.occ.addPoint(ht.oL.x, ht.oL.y, 0)
hts.areas[ht_nr].loop = self.occ.addCurveLoop(
[hts.lines[ht_nr + 'i'], # inner
hts.lines[ht_nr + 'l'], # lower
hts.lines[ht_nr + 'o'], # outer
hts.lines[ht_nr + 'h']]) # higher
# Build wire order of the insulation lines of the current block
if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
ht_list = list(hts.areas.keys())
ht_list.extend(list(reversed(ht_list))[1:])
self.blk_ins_lines[block_nr] = ['l']
for nr, ht_nr in enumerate(ht_list):
if nr + 1 == winding.conductors_number: # end of first round
self.blk_ins_lines[block_nr].extend([ht_nr + 'i', 'h', ht_nr + 'o'])
else:
if nr + 1 < winding.conductors_number: # within first round
self.blk_ins_lines[block_nr].extend([ht_nr + 'i', ht_nr + 'i' + ht_list[nr + 1]])
else: # within second round
self.blk_ins_lines[block_nr].extend([ht_nr + 'o' + ht_list[nr - 1], ht_nr + 'o'])
def constructInsulationGeometry(self):
"""
Generates points, hyper lines, and curve loops for the coil insulations
"""
def _createMidPoleLines(case, cnt=0):
if 'block-coil' in geom_coil.type:
if case == 'inner':
group.lines['mid_pole_' + case[0]] = self.occ.addLine(ins_pnt[first_ht_curr + case[0] + 'l'], ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
ordered_lines[group_nr].append(['mid_pole_' + case[0], (len(coil.layers) * 2) * 1e3 + 5e2, group.lines['mid_pole_' + case[0]]])
else:
group.lines['mid_pole_' + case[0]] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt_opposite[first_ht_prev + 'il'])
ordered_lines[group_nr].append(['mid_pole_' + case[0], 0, group.lines['mid_pole_' + case[0]]])
else:
ht_curr = geom_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[
block_order.block].half_turns[int(first_ht_curr)].corners.insulated
ht_prev = geom_coil.poles[block_order_prev.pole].layers[layer_nr].windings[block_order_prev.winding].blocks[
block_order_prev.block].half_turns[int(last_ht_prev)].corners.insulated
pnt_curr = [ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y]
pnt_prev = [ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y]
if Func.points_distance(pnt_curr, pnt_prev) > 1e-6:
correct_center = Func.corrected_arc_center([self.md.geometries.coil.coils[coil_nr].bore_center.x, self.md.geometries.coil.coils[coil_nr].bore_center.y],
[ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y],
[ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y])
ln_name = 'mid_pole_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]
group.lines[ln_name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],
self.occ.addPoint(correct_center[0], correct_center[1], 0),
ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
# self.occ.addLine(ins_pnt[first_ht_curr + case[0] + 'l'], ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
cnt += 1 if case == 'inner' else -1
ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])
return cnt
def _createMidWindingLines(case, cnt):
name = 'mid_wind_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]
# Create corrected center
blk1 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order.block)][0]].layers[
blks_info[str(block_order.block)][1]].windings[blks_info[str(block_order.block)][2]].blocks[int(str(block_order.block))]
blk2 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order_prev.block)][0]].layers[
blks_info[str(block_order_prev.block)][1]].windings[blks_info[str(block_order_prev.block)][2]].blocks[int(block_order_prev.block)]
pnt1 = blk1.half_turns[int(first_ht_curr)].corners.insulated.iL if case == 'inner' else blk1.half_turns[int(first_ht_curr)].corners.insulated.oL
pnt2 = blk2.half_turns[int(last_ht_prev)].corners.insulated.iH if case == 'inner' else blk2.half_turns[int(last_ht_prev)].corners.insulated.oH
outer_center = Func.corrected_arc_center([self.md.geometries.coil.coils[coil_nr].bore_center.x,
self.md.geometries.coil.coils[coil_nr].bore_center.y],
[pnt1.x, pnt1.y], [pnt2.x, pnt2.y])
group.lines[name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],
self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
cnt += 1 if case == 'inner' else -1
ordered_lines[group_nr].append([name, cnt, group.lines[name]])
return cnt
def _createInnerOuterLines(case, cnt):
# Create half turn lines
idxs = [1, round(len(self.blk_ins_lines[block_order.block]) / 2), 1] if case == 'inner'\
else [len(self.blk_ins_lines[block_order.block]) - 1, round(len(self.blk_ins_lines[block_order.block]) / 2), -1]
lns = self.blk_ins_lines[block_order.block][idxs[0]:idxs[1]:idxs[2]]
for ln_nr, ln_name in enumerate(lns):
skip_cnt = False
if ln_name[-1].isdigit():
try:
group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + 'h'],
ins_pnt[ln_name[ln_name.index(case[0]) + 1:] + case[0] + 'l'])
except:
skip_cnt = True
next_line = lns[ln_nr + 1]
pos = 'first' if next_line[:-1] == ln_name[:ln_name.index(case[0])] else 'second'
lns[ln_nr + 1] = next_line + (ln_name[ln_name.index(case[0]) + 1:] + 'l' if pos == 'first' else ln_name[:ln_name.index(case[0])] + 'h')
elif ln_name[-1] in ['i', 'o']:
group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name + 'l'], ins_pnt[ln_name + 'h'])
else:
group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + ln_name[-1]],
ins_pnt[ln_name[ln_name.index(case[0]) + 1:-1] + case[0] + ln_name[-1]])
if not skip_cnt:
cnt += 1 if case == 'inner' else -1
ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])
return cnt
def _computePointAngle(case):
points_angles = pa_next if case == 'outer' else pa_prev
current_ht_h = [current_ht.oH.x, current_ht.oH.y] if case == 'outer' else [current_ht.iH.x, current_ht.iH.y]
if ht_nr == 0:
current_ht_l = [current_ht.oL.x, current_ht.oL.y] if case == 'outer' else [current_ht.iL.x, current_ht.iL.y]
if 'block-coil' in geom_coil.type: current_ht_l[1] = 1 if current_ht_l[1] > 0 else -1
points_angles[str(block_order.block) + '_' + ht_name + 'l'] = Func.arc_angle_between_point_and_abscissa(current_ht_l, center)
if ht_nr == len(ht_list) - 1:
name = ht_name + 'h'
coord = current_ht_h
else: # for mid half turns, get the outer corner
next_ht_ins = geom_hts[int(ht_list[ht_nr + 1])].corners.insulated
next_ht = [next_ht_ins.oL.x, next_ht_ins.oL.y] if case == 'outer' else [next_ht_ins.iL.x, next_ht_ins.iL.y]
condition = (Func.points_distance(current_ht_h, center) > Func.points_distance(next_ht, center))\
if case == 'outer' else (Func.points_distance(current_ht_h, center) < Func.points_distance(next_ht, center))
if condition:
name = ht_name + 'h'
coord = current_ht_h
else:
name = ht_list[ht_nr + 1] + 'l'
coord = next_ht
if 'block-coil' in geom_coil.type: coord[1] = 1 if coord[1] > 0 else -1
points_angles[str(block_order.block) + '_' + name] = Func.arc_angle_between_point_and_abscissa(coord, center)
ins = self.md.geometries.insulation
for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
aux_coil = self.md.geometries.coil.coils[coil_nr]
geom_coil = self.geom.coil.coils[coil_nr]
groups = len(geom_coil.poles)
count = {}
ordered_lines = {}
points_angle = {}
blks_info = {}
ending_line = {}
center = [geom_coil.bore_center.x, geom_coil.bore_center.y]
if coil_nr not in ins.coils:
ins.coils[coil_nr] = dM.InsulationGroup()
ins_groups = ins.coils[coil_nr].group
for layer_nr, layer in coil.layers.items():
group_nr = 1
wnd_nr = len(aux_coil.poles[1].layers[layer_nr].windings)
ordered_layer = layer[wnd_nr:] + layer[:wnd_nr] if layer[0].pole != layer[-1].pole else layer
for nr, block_order in enumerate(ordered_layer):
blks_info[str(block_order.block)] = [block_order.pole, layer_nr, block_order.winding]
# Get previous block in anticlockwise order
block_order_prev = ordered_layer[-1] if nr == 0 else ordered_layer[nr - 1]
# Update insulation group
if block_order.winding == block_order_prev.winding:
group_nr = group_nr + 1 if group_nr < groups else 1
# Initialize dicts
if group_nr not in ins_groups:
ins_groups[group_nr] = dM.InsulationRegion()
points_angle[group_nr] = {}
ordered_lines[group_nr] = []
count[group_nr] = [0, (len(coil.layers) + 1) * 1e3]
group = ins_groups[group_nr].ins
ins_groups[group_nr].blocks.append([block_order.pole, layer_nr, block_order.winding, block_order.block])
# Find the wedge
if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
for wdg, blk in self.md.geometries.wedges.coils[coil_nr].layers[layer_nr].block_prev.items():
if blk == block_order_prev.block:
ins_groups[group_nr].wedges.append([layer_nr, wdg])
break
if layer_nr < len(coil.layers):
mid_layer_next = str(layer_nr) + '_' + str(layer_nr + 1)
if mid_layer_next not in points_angle[group_nr]:
points_angle[group_nr][mid_layer_next] = {}
pa_next = points_angle[group_nr][mid_layer_next]
if layer_nr > 1:
mid_layer_prev = str(layer_nr - 1) + '_' + str(layer_nr)
pa_prev = points_angle[group_nr][mid_layer_prev]
# Get point tags of insulation
ins_pnt = aux_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[
block_order.block].insulation.points
# Get relevant info for line names
first_ht_curr = self.blk_ins_lines[block_order.block][1][:-1]
last_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[
layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[-1]
ins_pnt_opposite = aux_coil.poles[block_order_prev.pole].layers[
layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].insulation.points
if 'cos-theta' == geom_coil.type:
# Create lower and higher angle lines
if block_order.winding == block_order_prev.winding:
group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])
ending_line[group_nr - 1 if group_nr > 1 else groups] =\
[ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]
# Create inner lines of insulation group
if layer_nr == 1:
if block_order.pole != block_order_prev.pole:
count[group_nr][0] = _createMidPoleLines('inner', count[group_nr][0])
if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
count[group_nr][0] = _createMidWindingLines('inner', count[group_nr][0])
count[group_nr][0] = _createInnerOuterLines('inner', count[group_nr][0])
# Create outer lines of insulation group
if layer_nr == len(coil.layers):
if block_order.pole != block_order_prev.pole:
count[group_nr][1] = _createMidPoleLines('outer', count[group_nr][1])
if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
count[group_nr][1] = _createMidWindingLines('outer', count[group_nr][1])
count[group_nr][1] = _createInnerOuterLines('outer', count[group_nr][1])
elif 'block-coil' in geom_coil.type:
last_ht_curr = self.blk_ins_lines[block_order.block][self.blk_ins_lines[block_order.block].index('h') - 1][:-1]
first_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[layer_nr].windings[
block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[0]
# Create lower and higher angle lines
if block_order.winding == block_order_prev.winding:
group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 4 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])
ending_line[group_nr - 1 if group_nr > 1 else groups] =\
[ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]
group.lines[str(layer_nr) + 'bh'] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt[last_ht_curr + 'oh'])
ordered_lines[group_nr].append([str(layer_nr) + 'bh', (len(coil.layers) * 2 + layer_nr) * 1e3, group.lines[str(layer_nr) + 'bh']])
# Create inner lines of insulation group
if block_order.pole != block_order_prev.pole:
if layer_nr == 1:
_createMidPoleLines('inner')
_createMidPoleLines('outer')
group.lines[str(layer_nr) + 'bl'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
ordered_lines[group_nr].append([str(layer_nr) + 'bl', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'bl']])
# Create outer lines of insulation group
if layer_nr == len(coil.layers):
count[group_nr][1] = _createInnerOuterLines(
'outer', (len(coil.layers) * 4 - layer_nr + 1) * 1e3 if block_order.winding == block_order_prev.winding else (len(coil.layers) + 1) * 1e3)
# Store info about the angle of each point in between layers
ht_list = list(aux_coil.poles[block_order.pole].layers[
layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns.areas.keys())
geom_hts = geom_coil.poles[block_order.pole].layers[
layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns
for ht_nr, ht_name in enumerate(ht_list): # half turns in anticlockwise order
current_ht = geom_hts[int(ht_name)].corners.insulated
if layer_nr < len(coil.layers): # if it's not the last layer, fetch all outer corners angles
_computePointAngle('outer')
if layer_nr > 1: # if it's not the first layer, fetch all inner corners angles
_computePointAngle('inner')
# Create closing lines
for grp_nr, grp in ending_line.items():
ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h'] = self.occ.addLine(grp[0], grp[1])
ordered_lines[grp_nr].append([str(layer_nr) + 'h', layer_nr * 1e3, ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h']])
# Create lines connecting different layers and generate closed loops
for group_nr, group in points_angle.items():
ins_group = ins_groups[group_nr].ins
for mid_l_name, mid_l in group.items():
first_layer = mid_l_name[:mid_l_name.index('_')]
# Correct angles if the group crosses the abscissa
max_angle = max(mid_l.values())
max_diff = max_angle - min(mid_l.values())
if max_diff > np.pi:
for pnt_name, angle in mid_l.items():
if angle < max_diff / 2:
mid_l[pnt_name] = angle + max_angle
# Order points according to angle
ordered_pnts = [[pnt_name, angle] for pnt_name, angle in mid_l.items()]
ordered_pnts.sort(key=lambda x: x[1])
ordered_names = [x[0] for x in ordered_pnts]
for case in ['beg', 'end']:
past_blocks = []
sides = ['l', 'o', 'h', 'l'] if case == 'beg' else ['h', 'i', 'l', 'h']
# count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2
for i in range(2 if 'block-coil' in geom_coil.type else 1):
count = int(first_layer) * 1e3 + 5e2 if i == 0 else (len(coil.layers) * 2 + int(first_layer)) * 1e3 + 5e2
if case == 'beg':
pnt_position = 0 if i == 0 else int(len(ordered_names) / 2)
else:
pnt_position = -1 if i == 0 else int(len(ordered_names) / 2 - 1)
first_block = ordered_names[pnt_position][:ordered_names[pnt_position].index('_')] # ordered_pnts[pnt_position][0][:ordered_pnts[pnt_position][0].index('_')] #
ordered_search_names = ordered_names[pnt_position::1 if case == 'beg' else -1]
for nr, pnt in enumerate(ordered_search_names[1:], 1): # enumerate(ordered_names if case == 'beg' else reversed(ordered_names)): #
current_blk = pnt[:pnt.index('_')]
ins_pnt = aux_coil.poles[blks_info[current_blk][0]].layers[blks_info[current_blk][1]].windings[
blks_info[current_blk][2]].blocks[int(current_blk)].insulation.points
prev_pnt = ordered_search_names[nr - 1] # ordered_pnts[nr - 1 if case == 'beg' else - nr][0] #
prev_blk = prev_pnt[:prev_pnt.index('_')]
start_pnt_name = prev_pnt[prev_pnt.index('_') + 1:-1] + ('o' if str(blks_info[prev_blk][1]) == first_layer else 'i')
ins_pnt_prev = aux_coil.poles[blks_info[prev_blk][0]].layers[blks_info[prev_blk][1]].windings[
blks_info[prev_blk][2]].blocks[int(prev_blk)].insulation.points
# Create lines when you find the first edge belonging to a block of the opposite layer
if blks_info[current_blk][1] != blks_info[first_block][1]:
pnt_tag_name = pnt[pnt.index('_') + 1:-1] + ('o' if str(blks_info[current_blk][1]) == first_layer else 'i') + ('l' if pnt[-1] == 'l' else 'h')
pnt_tag_name_opposite = start_pnt_name + ('l' if prev_pnt[-1] == 'l' else 'h')
opp_blk_ins_lines = self.blk_ins_lines[int(prev_blk)]
indexes = [opp_blk_ins_lines.index(start_pnt_name) + (1 if prev_pnt[-1] == sides[0] else 0),
len(opp_blk_ins_lines) if case == 'beg' else opp_blk_ins_lines.index('h'), 1] if start_pnt_name[-1] == sides[1]\
else [opp_blk_ins_lines.index(start_pnt_name) - (1 if prev_pnt[-1] == sides[0] else 0),
0 if case == 'beg' else opp_blk_ins_lines.index('h'), -1]
if case == 'beg':
if i == 0:
count = (len(coil.layers) * (4 if 'block-coil' in geom_coil.type else 2) - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])
else:
count = (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])
else:
count += 1 + abs(indexes[0] - indexes[1])
# Create all remaining lines of the current layer block
for line_name in opp_blk_ins_lines[indexes[0]:indexes[1]:indexes[2]]:
if 'block-coil' in geom_coil.type:
if not line_name[-1].isdigit():
ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])
count += 1 if (case == 'beg' and i == 1) or (case == 'end' and i == 0) else -1
ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
else:
if line_name[-1].isdigit():
ins_group.lines[line_name] = self.occ.addLine(
ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'],
ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'])
else:
ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])
count += 1 if case == 'beg' else -1 # if start_pnt_name[-1] == sides[1] else 1
ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
# Create mid layer line
if 'block-coil' in geom_coil.type:
count_rest = -abs(indexes[0] - indexes[1]) if (case == 'beg' and i == 1) or (case == 'end' and i == 0) else 1 + abs(indexes[0] - indexes[1])
else:
count_rest = -abs(indexes[0] - indexes[1]) if case == 'beg' else 1 + abs(indexes[0] - indexes[1])
line_name = 'mid_layer_' + mid_l_name + ('b' if i == 1 else '') + (
'_l' if case == 'beg' else '_h')
gmsh.model.occ.synchronize()
"""
The line that connects two layers may overlap with a new (pole) region. This can be avoided by
modifying the line to be an L shape by adding an intermediate point.
We add the point along the upper half-turn radial direction towards the center.
"""
p = ins_pnt[pnt_tag_name] # point to be extended
linetag = gmsh.model.getAdjacencies(0, ins_pnt[pnt_tag_name])[0][0] # line to be extended
# find the distance to move the point towards the center
coord1 = gmsh.model.getValue(0, p, [])
coord2 = gmsh.model.getValue(0, ins_pnt_prev[pnt_tag_name_opposite], [])
# this is the
r1 = np.sqrt(coord1[0] ** 2 + coord1[1] ** 2)
r2 = np.sqrt(coord2[0] ** 2 + coord2[1] ** 2)
distance = (r1 - r2) / 2
p1, p2 = [b[1] for b in gmsh.model.getBoundary([(1, linetag)], oriented=True)]
dir_vector = gmsh.model.getValue(0, p1, [])-gmsh.model.getValue(0, p2, [])
unit_vector = dir_vector / np.linalg.norm(dir_vector)
coord = gmsh.model.getValue(0, p, [])
X = self.occ.addPoint(
coord[0] + unit_vector[0] * distance,
coord[1] + unit_vector[1] * distance,
0)
gmsh.model.occ.synchronize()
ins_group.lines[line_name + "_A"] = self.occ.addLine(ins_pnt[pnt_tag_name], X)
ins_group.lines[line_name + "_B"] = self.occ.addLine(X, ins_pnt_prev[pnt_tag_name_opposite])
ordered_lines[group_nr].append(
[line_name+"_A", count + count_rest -1 if case == 'beg' else count + count_rest, ins_group.lines[line_name+"_A"]])
ordered_lines[group_nr].append(
[line_name+"_B", count + count_rest if case == 'beg' else count + count_rest -1, ins_group.lines[line_name+"_B"]])
""" # original code
line_name = 'mid_layer_' + mid_l_name + ('b' if i == 1 else '') + ('_l' if case == 'beg' else '_h')
ins_group.lines[line_name] = self.occ.addLine(ins_pnt[pnt_tag_name], ins_pnt_prev[pnt_tag_name_opposite])
ordered_lines[group_nr].append([line_name, count + count_rest, ins_group.lines[line_name]])
"""
break
# Create all edges of the first block sticking out completely todo: might have to be extended to multiple blocks
if current_blk != first_block and current_blk not in past_blocks:
def __createWedgeInsulation(cnt):
# Create the line connecting the blocks (where a wedge is)
line_name = self.blk_ins_lines[int(current_blk)][
(-1 if start_pnt_name[-1] == 'o' else 1) if case == 'beg'
else (round(len(self.blk_ins_lines[int(current_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1))]
line_name_prev = self.blk_ins_lines[int(prev_blk)][
(round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1)) if case == 'beg'
else (-1 if start_pnt_name[-1] == 'o' else 1)]
# Create corrected center
blk1 = geom_coil.poles[blks_info[prev_blk][0]].layers[
blks_info[prev_blk][1]].windings[blks_info[prev_blk][2]].blocks[int(prev_blk)]
blk2 = geom_coil.poles[blks_info[current_blk][0]].layers[
blks_info[current_blk][1]].windings[blks_info[current_blk][2]].blocks[int(current_blk)]
pnt1 = blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oH if case == 'beg'\
else blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oL
pnt2 = blk2.half_turns[int(line_name[:-1])].corners.insulated.oL if case == 'beg'\
else blk2.half_turns[int(line_name[:-1])].corners.insulated.oH
outer_center = Func.corrected_arc_center([aux_coil.bore_center.x, aux_coil.bore_center.y],
[pnt2.x, pnt2.y] if case == 'beg' else [pnt1.x, pnt1.y],
[pnt1.x, pnt1.y] if case == 'beg' else [pnt2.x, pnt2.y])
ins_group.lines[line_name_prev + line_name] =\
self.occ.addCircleArc(ins_pnt_prev[line_name_prev + sides[2]],
self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt[line_name + sides[3]])
ordered_lines[group_nr].append([line_name_prev + line_name, cnt, ins_group.lines[line_name_prev + line_name]])
count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2
past_blocks.append(current_blk)
indexes = [round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + 1,
len(self.blk_ins_lines[int(prev_blk)])] if str(blks_info[prev_blk][1]) == first_layer\
else [1, round(len(self.blk_ins_lines[int(prev_blk)]) / 2)]
if case == 'beg':
count += 1
__createWedgeInsulation(count)
lines = self.blk_ins_lines[int(prev_blk)][indexes[0]:indexes[1]]
side = 'o' if str(blks_info[prev_blk][1]) == first_layer else 'i'
for line_nr, line_name in enumerate(lines):
skip_count = False
if line_name[-1].isdigit():
try:
ins_group.lines[line_name] =\
self.occ.addLine(ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'],
ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'])
except: # points are too close to each other
skip_count = True
next_line = lines[line_nr + 1]
pnt1, pnt2 = line_name.split(side)
pos = 'first' if next_line[:-1] == pnt1 else 'second'
lines[line_nr + 1] = next_line + (pnt2 + 'l' if pos == 'first' else pnt1 + 'h')
elif line_name[-1] in ['i', 'o']:
ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'h'], ins_pnt_prev[line_name + 'l'])
else:
ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name[:line_name.index(side)] + side + line_name[-1]],
ins_pnt_prev[line_name[line_name.index(side) + 1:-1] + side + line_name[-1]])
if not skip_count:
count += 1 # if start_pnt_name[-1] == sides[1] else -1
ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
if case == 'end':
count += 1
__createWedgeInsulation(count)
# Generate closed loops
ordered_lines[group_nr].sort(key=lambda x: x[1])
area_name = str((coil_nr - 1) * len(ins_groups) + group_nr)
ins_group.areas[area_name] = dM.Area()
if len(points_angle) == 1:
ins_group.areas['inner_loop'] = dM.Area(loop=self.occ.addCurveLoop([ins_group.lines[line] for line in [x[0] for x in ordered_lines[group_nr]]
if 'i' in line and line[0].isdigit() or '_i' in line]))
ins_group.areas[area_name].loop = self.occ.addCurveLoop([ins_group.lines[line] for line in [x[0] for x in ordered_lines[group_nr]]
if 'o' in line and line[0].isdigit() or '_o' in line])
else:
ins_group.areas[area_name].loop = self.occ.addCurveLoop([ins_group.lines[line] for line in [x[0] for x in ordered_lines[group_nr]]])
def constructThinShells_poles(self):
ts_layer = self.md.geometries.thin_shells.pole_layers
ts_av_ins_thick = self.md.geometries.thin_shells.ins_thickness.poles
def _construct_thin_shell_corners_to_line(pnt1, pnt2, pole_line, name):
# use gmsh to calculate distance to a line
coord_a = gmsh.model.getClosestPoint(1, pole_line, coord=(pnt1[0], pnt1[1], 0))[0]
coord_b = gmsh.model.getClosestPoint(1, pole_line, coord=(pnt2[0], pnt2[1], 0))[0]
# draw new point at half the distance between iH and coord_a
new_i = self.occ.addPoint((pnt1[0] + coord_a[0]) / 2, (pnt1[1] + coord_a[1]) / 2, 0)
new_o = self.occ.addPoint((pnt2[0] + coord_b[0]) / 2, (pnt2[1] + coord_b[1]) / 2, 0)
ts_layer[name] = dM.Region()
self.occ.synchronize()
ts_layer[name].lines['1'] = self.occ.addLine(new_i, new_o)
self.occ.synchronize()
cond_name = next(iter(self.data.conductors.keys()))
other_material = 0.5*(self.data.conductors[cond_name].cable.th_insulation_along_height+self.data.conductors[cond_name].cable.th_insulation_along_width)#todo, better select which one
# distance -> Average between coord_a and iH and coord_b and oH AND remove the G10 thickness
ts_av_ins_thick[name] = float(0.5 * (np.sqrt((pnt1[0] - coord_a[0]) ** 2 + (pnt1[1] - coord_a[1]) ** 2) +
np.sqrt((pnt2[0] - coord_b[0]) ** 2 + (pnt2[1] - coord_b[1]) ** 2))) - other_material
return
def _find_line_closest_to_points(pnt1, pnt2, line_list_tags):
"""
Should work for any (pole) geometry. Given the half turn corner points pnt1 = [x1, y1], pnt2 = [x2, y2],
and a list of pole line tags, we need to select which one is on the opposite side to construct the thin shell line
thus, we search the closest line(s) to each corner point and then take the intersection of those two sets
"""
closest_lines = {0: [], 1: []}
for i, pnt in enumerate([pnt1, pnt2]):
min_dist = float('inf')
for line_tag in line_list_tags:
tag1, tag2 = gmsh.model.getAdjacencies(1, line_tag)[1]
start_pnt = gmsh.model.getValue(0, tag1, [])
end_pnt = gmsh.model.getValue(0, tag2, [])
v = end_pnt - start_pnt
w = pnt - start_pnt
c1 = np.dot(w, v)
c2 = np.dot(v, v)
# avoid extending the line
if c1 <= 0:
dist = np.linalg.norm(pnt - start_pnt) # Closest to p1
elif c2 <= c1:
dist = np.linalg.norm(pnt - end_pnt) # Closest to p2
else:
b = c1 / c2
Pb = start_pnt + b * v
dist = np.linalg.norm(pnt - Pb)
if np.isclose(min_dist, dist):
closest_lines[i].append(line_tag) # add to list
elif dist < min_dist:
min_dist = dist
closest_lines[i] = [line_tag] # overwrite
# return the intersection of closest lines
return list(set(closest_lines[0]).intersection(set(closest_lines[1])))[0] # should be only one line
def _split_lines_azimuthal_radial(line_tag_list):
alines = []
rlines = []
for tag in line_tag_list:
pointTags = gmsh.model.getAdjacencies(1, tag)[1] # get tag
p1 = gmsh.model.getValue(0, pointTags[0], [])
p2 = gmsh.model.getValue(0, pointTags[1], [])
dr = np.sqrt(p2[0] ** 2 + p2[1] ** 2) - np.sqrt(p1[0] ** 2 + p1[1] ** 2)
dt = (np.arctan2(p2[1], p2[0]) - np.arctan2(p1[1], p1[0])) * np.sqrt(p2[0] ** 2 + p2[1] ** 2) ## convert to length
if np.abs(dt)>np.abs(dr):
alines.append(tag)
else:
rlines.append(tag)
return alines, rlines
def _wedge_to_pole_lines():
"""
These lines do not exist, those that already exist contain the G10 insulation so we need more.
"""
for _, region in self.md.geometries.thin_shells.mid_layers_aux.items(): #those are the g10 layers, we need to duplicate this line to account for the kapton
# obtain the position of the lines and draw more on top of it
for name, line_tag in region.lines.items():
name = f"p{name}_r" # all radial lines
ts_layer[name] = dM.Region()
point_tag = gmsh.model.getBoundary([(1, line_tag)], oriented=False)
p1, p2 = point_tag[0][1], point_tag[1][1]
x1, y1, _ = gmsh.model.getValue(0, p1, [])
x2, y2, _ = gmsh.model.getValue(0, p2, [])
# create new line
p1_new = gmsh.model.occ.addPoint(x1, y1, 0.0)
p2_new = gmsh.model.occ.addPoint(x2, y2, 0.0)
self.occ.synchronize()
ts_layer[name].lines['1'] = self.occ.addLine(p1_new, p2_new)
self.occ.synchronize()
# thickness = distance between the wedge and the pole - G10 thickness
# how to approximate this thickness? -> use the same thickness as for the ht 108
"""
cond_name = next(iter(self.data.conductors.keys()))
other_material = 0.5 * (
self.data.conductors[cond_name].cable.th_insulation_along_height + self.data.conductors[
cond_name].cable.th_insulation_along_width)
# distance -> Average between coord_a and iH and coord_b and oH AND remove the G10 thickness ?
ts_av_ins_thick[name] = ???
"""
ts_av_ins_thick[name] = ts_av_ins_thick['p108_a'] #@emma hardcoded thickness, use the same as for this ht
max_layer = max(self.geom.coil.coils[1].poles[1].layers.keys())
# first, the lines connecting hts to the pole
for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items(): # coilnr is only 1 here
for layer_nr, layer in coil.layers.items(): # we need to add radial TS lines for both layers
for nr, blk_order in enumerate(layer):
block = self.geom.coil.coils[coil_nr].poles[blk_order.pole].layers[
layer_nr].windings[blk_order.winding].blocks[blk_order.block]
ht_list = list(self.md.geometries.coil.coils[coil_nr].poles[blk_order.pole].layers[
layer_nr].windings[blk_order.winding].blocks[
blk_order.block].half_turns.areas.keys())
blk_index_next = nr + 1 if nr + 1 < len(layer) else 0
block_order_next = layer[blk_index_next]
block_next = self.geom.coil.coils[coil_nr].poles[block_order_next.pole].layers[
layer_nr].windings[block_order_next.winding].blocks[block_order_next.block]
ht_list_next = list(self.md.geometries.coil.coils[coil_nr].poles[block_order_next.pole].layers[
layer_nr].windings[block_order_next.winding].blocks[
block_order_next.block].half_turns.areas.keys())
ht_last = int(ht_list[-1])
ht_next_first = int(ht_list_next[0])
# if winding nr is the same and pole is the same -> pole connection
if blk_order.pole == block_order_next.pole and blk_order.winding == block_order_next.winding:
# Create radial lines for pole connection
pole_lines_a, pole_lines_r = _split_lines_azimuthal_radial(
[v for qq in self.md.geometries.poles.quadrants.keys() for v in
self.md.geometries.poles.quadrants[qq].lines.values()])
# todo maybe this can be speed up if we select the correct quadrant
if layer_nr != max_layer:
# add azimuthal thin shells at the ht side 'o' (normals radial)
# Now we consider all the HT's from one block, this might not hold true in general
for ht in map(int, ht_list):
oH = [block.half_turns[ht].corners.bare.oH.x,
block.half_turns[ht].corners.bare.oH.y, 0.0]
oL = [block.half_turns[ht].corners.bare.oL.x,
block.half_turns[ht].corners.bare.oL.y, 0.0]
_construct_thin_shell_corners_to_line(
oH, oL,
_find_line_closest_to_points(oH, oL, pole_lines_a),
name=f"p{ht}_r")
for ht in map(int, ht_list_next):
oH = [block_next.half_turns[ht].corners.bare.oH.x,
block_next.half_turns[ht].corners.bare.oH.y, 0.0]
oL = [block_next.half_turns[ht].corners.bare.oL.x,
block_next.half_turns[ht].corners.bare.oL.y, 0.0]
_construct_thin_shell_corners_to_line(
oH, oL,
_find_line_closest_to_points(oH, oL, pole_lines_a),
name=f"p{ht}_r")
iH = [block.half_turns[ht_last].corners.bare.iH.x,
block.half_turns[ht_last].corners.bare.iH.y, 0.0]
oH = [block.half_turns[ht_last].corners.bare.oH.x,
block.half_turns[ht_last].corners.bare.oH.y, 0.0]
_construct_thin_shell_corners_to_line(
iH, oH,
_find_line_closest_to_points(iH, oH, pole_lines_r),
name=f"p{ht_last}_a")
iL = [block_next.half_turns[ht_next_first].corners.bare.iL.x,
block_next.half_turns[ht_next_first].corners.bare.iL.y, 0.0]
oL = [block_next.half_turns[ht_next_first].corners.bare.oL.x,
block_next.half_turns[ht_next_first].corners.bare.oL.y, 0.0]
_construct_thin_shell_corners_to_line(
iL, oL,
_find_line_closest_to_points(iL, oL, pole_lines_r),
name=f"p{ht_next_first}_a")
# second, the lines connecting wedge to the pole
_wedge_to_pole_lines()
def constructAdditionalThinShells(self):
def _create_lines_ht_alignment(ohx, ohy, olx, oly, ihx, ihy, ilx, ily):
def __line_circle_intersection(p1, p2):
x1, y1 = p1
x2, y2 = p2
dx = x2 - x1
dy = y2 - y1
A = dx**2 + dy**2
B = 2 * (dx*x1 + dy*y1)
C = x1**2 + y1**2 - R**2 # R collar is known from outer scope
disc = B**2 - 4*A*C
if disc < 0:
logger.warning(" No intersection between line and circle.")
return [] # no intersection
else:
sqrt_disc = np.sqrt(disc)
t1 = (-B + sqrt_disc) / (2*A)
t2 = (-B - sqrt_disc) / (2*A)
return [(x1 + t1*dx, y1 + t1*dy),
(x1 + t2*dx, y1 + t2*dy)]
mid_bare_o = np.mean([[ohx, ohy], [olx, oly]], axis=0)
mid_bare_i = np.mean([[ihx, ihy], [ilx, ily]], axis=0)
offsets = np.array([[ohx, ohy], [olx, oly]]) - mid_bare_o
quad = 1 if mid_bare_o[0] >= 0 and mid_bare_o[1] >= 0 else (
2 if mid_bare_o[0] <= 0 <= mid_bare_o[1] else (3 if mid_bare_o[0] <= 0 and mid_bare_o[1] <= 0 else 4))
inters = __line_circle_intersection(mid_bare_o, mid_bare_i)
quadrants = {
1: lambda x, y: x >= 0 and y >= 0,
2: lambda x, y: x <= 0 <= y,
3: lambda x, y: x <= 0 and y <= 0,
4: lambda x, y: x >= 0 >= y,
}
check = quadrants[quad]
for x, y in inters:
if check(x, y): # select the correct intersection based on the quadrant
xi, yi = x, y
break
# add a point halfway between (x,y) and (xi, yi)
dr = float(np.sqrt((xi - mid_bare_o[0]) ** 2 + (yi - mid_bare_o[1]) ** 2))
new_point_coords = [mid_bare_o[0] + 0.5 * (xi - mid_bare_o[0]), mid_bare_o[1] + 0.5 * (yi - mid_bare_o[1])]
for offset in offsets:
point_tag_final = self.occ.addPoint(new_point_coords[0] + offset[0], new_point_coords[1] + offset[1], 0)
self.occ.synchronize()
return dr, point_tag_final
def _embed_points_to_collar_curve():
def __add_point_in_curve(points, curve_tag):
return self.occ.fragment([(0, k) for k in points], [(1, curve_tag)], removeObject=True)[0]
### First, cut the collar line
self.occ.synchronize()
new_tags = {1: [], 2: [], 3: [], 4: []} # tuples
new_point_tags = {1: [], 2: [], 3: [], 4: []} # only tags
collar_size = self.data.magnet.mesh.thermal.collar.SizeMin
for ts_name, ts in self.md.geometries.thin_shells.collar_layers.items():
for name, line in ts.lines.items():
coords = gmsh.model.getValue(1, line, [i[0] for i in gmsh.model.getParametrizationBounds(1, line)])
t1 = np.arctan2(coords[0], coords[1])
t2 = np.arctan2(coords[3], coords[4])
quad = t1 // (np.pi / 2) + 1 if t1 > 0 else 4 + t1 // (np.pi / 2) + 1
curve_tag = self.md.geometries.collar.inner_boundary_tags[quad][0]
start, end = min(t1 % (np.pi / 2), t2 % (np.pi / 2)), max(t1 % (np.pi / 2),
t2 % (np.pi / 2)) # ensure start < end
tmp_coords = gmsh.model.getValue(1, line, [i[0] for i in gmsh.model.getParametrizationBounds(1, line)])
target_size = collar_size
elements = max(1, round(Func.points_distance(tmp_coords[:2], tmp_coords[3:-1]) / target_size))+1
if elements%2 == 1: elements += 1
para_coords = np.linspace(start, end, elements, endpoint=True)
for u in para_coords:
if quad == 1 or quad == 3:
u = np.pi / 2 - u ## magic
x, y, z = gmsh.model.getValue(1, curve_tag, [u]) # Evaluate point on curve
new_point_tags[quad].append(self.occ.addPoint(x, y, z)) # Add point coordinates to the list
for q in new_point_tags.keys():
new_tags[q].extend(__add_point_in_curve(new_point_tags[q], self.md.geometries.collar.inner_boundary_tags[q][0]))
# so for occ the old curve no longer exists, but in the objects the tag is still there
REMOVED_TAGS = [self.md.geometries.collar.inner_boundary_tags[q][0] for q in new_tags.keys()]
# update boundary tags
for quad, taglist in new_tags.items():
curvelist = [tag[1] for tag in taglist if tag[0]==1 ] # select the curves only
self.md.geometries.collar.inner_boundary_tags[quad] = curvelist # replace
for k, new_tag in enumerate(self.md.geometries.collar.inner_boundary_tags[quad]):
self.md.geometries.collar.quadrants[quad].lines[str(k)] = new_tag ## adding new lines
# REDEFINE THE CURVELOOP
loop_list = []
tmpdict = copy.deepcopy(self.md.geometries.collar.quadrants[quad].lines)
for tag_name, tag in tmpdict.items():
if tag in self.md.geometries.collar.cooling_tags: # skip holes
continue
if tag in REMOVED_TAGS:
# remove from dictionary
del self.md.geometries.collar.quadrants[quad].lines[tag_name]
continue
loop_list.append(tag)
self.occ.synchronize()
for area in self.md.geometries.collar.quadrants[quad].areas:
if area.startswith('arc') and not area.startswith('arch'): #collar region, not the holes
self.md.geometries.collar.quadrants[quad].areas[area] = dM.Area(
loop=self.occ.addCurveLoop(loop_list))
# Cut the pole lines
"""
idea, loop over the thin shell lines, then draw line from the HT corner to the line ends, then intersect it with the pole curve
just find intersection and select the shortest one ?
# not implemented, maybe not necessary
"""
for quad in [1, 2, 3, 4]:
# additionally add the pole area back
for area in self.md.geometries.poles.quadrants[quad].areas:
if area.startswith('arp'):
self.md.geometries.poles.quadrants[quad].areas[area] = dM.Area(
loop=self.occ.addCurveLoop(list(self.md.geometries.poles.quadrants[quad].lines.values())))
self.occ.synchronize()
# point still needs to be added
ts_layer = self.md.geometries.thin_shells.collar_layers
ts_av_ins_thick = self.md.geometries.thin_shells.ins_thickness.collar
collar_tag_dict = self.md.geometries.collar.inner_boundary_tags
enforce_TSA_mapping_collar = self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping # if True, cut the collar curve into segments, otherwise use the whole curve
center = self.occ.addPoint(0, 0, 0)
collar_x, collar_y, _ = gmsh.model.getValue(1, collar_tag_dict[1][0], [0])
R = np.sqrt(collar_x ** 2 + collar_y ** 2) # radius of collar curve
alignment = ['radial', 'ht'][1] # pick ht alignment
# COLLAR
for pid, pole in self.geom.coil.coils[1].poles.items():
layer_num = max(pole.layers.keys()) # outside layer
for wid, winding in pole.layers[layer_num].windings.items():
for block_idx in winding.blocks.keys():
block = winding.blocks[block_idx]
ht_nr_area = list(self.md.geometries.coil.coils[1].poles[pid].layers[layer_num].windings[wid].blocks[block_idx].half_turns.areas.keys())
i = 0
for ht_nr, ht in block.half_turns.items(): #ht_idx is not the same as number
ht_old = ht_nr_area[i]
i+=1
dr = 0.
if alignment == 'radial':
for (x, y), (x1, y1) in zip([[ht.corners.bare.oH.x, ht.corners.bare.oH.y], [ht.corners.bare.oL.x, ht.corners.bare.oL.y]],
[[ht.corners.bare.iH.x, ht.corners.bare.iH.y], [ht.corners.bare.iL.x, ht.corners.bare.iL.y]]): # uses the bare coordinates of the HT
t = np.arctan2(y, x)
quad = t // (np.pi / 2) + 1 if t > 0 else 4 + t // (np.pi / 2) + 1
collar_idx = collar_tag_dict[quad][0] # get the (debug: ONE) collar curve tag for the quadrant
dr_prev = dr
dummy = self.occ.addPoint(x, y, 0)
dr = self.occ.get_distance(0, dummy, 1, collar_idx)[0]
self.occ.synchronize()
self.occ.remove([(0, dummy)])
point_tag_final = self.occ.addPoint(x + 0.5 * dr * np.cos(t), y + 0.5 * dr * np.sin(t), 0)
elif alignment == 'ht': # distance to the next half turn (more relevant for coil to coil distances)
dr, point_tag_final = _create_lines_ht_alignment(ohx = ht.corners.bare.oH.x, ohy=ht.corners.bare.oH.y,
olx = ht.corners.bare.oL.x, oly=ht.corners.bare.oL.y,
ihx = ht.corners.bare.iH.x, ihy=ht.corners.bare.iH.y,
ilx = ht.corners.bare.iL.x, ily=ht.corners.bare.iL.y)
# save the line
name = f'{ht_nr}_x'
ts_layer[name] = dM.Region()
self.occ.synchronize()
ts_layer[name].lines['1'] = self.occ.addLine(point_tag_final-1, point_tag_final)
cond_name = next(iter(self.data.conductors.keys()))
if alignment == 'radial':
ts_av_ins_thick[name] = 0.5*(dr+dr_prev) - self.data.conductors[cond_name].cable.th_insulation_along_width # approx thickness, average - insulation
elif alignment == 'ht':
ts_av_ins_thick[name] = dr - self.data.conductors[cond_name].cable.th_insulation_along_width
self.occ.synchronize()
# WEDGES
if self.data.magnet.geometry.thermal.with_wedges:
for wedge_idx, wedge in self.md.geometries.wedges.coils[1].layers[max(self.md.geometries.wedges.coils[1].layers.keys())].wedges.items():
dr=0.
if alignment == 'radial':
raise NotImplementedError("Wedge radial alignment not implemented")
elif alignment == 'ht': # distance to the next half turn (more relevant for coil to coil distances)
ohx, ohy, _ = gmsh.model.getValue(0, wedge.points['oh'], [])
olx, oly, _ = gmsh.model.getValue(0, wedge.points['ol'], [])
ihx, ihy, _ = gmsh.model.getValue(0, wedge.points['ih'], [])
ilx, ily, _ = gmsh.model.getValue(0, wedge.points['il'], [])
dr, point_tag_final = _create_lines_ht_alignment(ohx=ohx, ohy=ohy, olx=olx, oly=oly,
ihx=ihx, ihy=ihy, ilx=ilx, ily=ily)
# find the smallest thickness
dr_b = 0.0
for x, y in [[ohx, ohy],
[olx, oly]]:
dummy = self.occ.addPoint(x, y, 0)
dr_a = dr_b
t = np.arctan2(y, x)
quad = t // (np.pi / 2) + 1 if t > 0 else 4 + t // (np.pi / 2) + 1
collar_idx = collar_tag_dict[quad][0]
dr_b = self.occ.get_distance(0, dummy, 1, collar_idx)[0]
self.occ.remove([(0, dummy)])
# save line
name = f'w{wedge_idx}_x'
ts_layer[name] = dM.Region()
ts_layer[name].lines['1'] = self.occ.addLine(point_tag_final - 1, point_tag_final) # make index line trivial (no distinction)
cond_name = next(iter(self.data.conductors.keys()))
#### ts_av_ins_thick[name] = dr - self.data.conductors[cond_name].cable.th_insulation_along_width
# This overshoots the thickness. The wedges are curved and the center between the corners (straight line) is further away from the collar than the curved wedge boundary.
# Maybe one could use gmsh to obtain the distance between the outer curve and the collar, but this should be close to taking the minimum distance of the cornerpoints as both curves are
# circle segments of (approximately) two concentric circles around the centre
logger = logging.getLogger('FiQuS')
logger.warning("Using alternative wedge insulation thickness approximation ")
ts_av_ins_thick[name] = min(dr_a,dr_b)- self.data.conductors[cond_name].cable.th_insulation_along_width # approx thickness, average - insulation
# POLES
if 'poles' in self.data.magnet.geometry.thermal.areas:
# generate additional thin shell lines
self.constructThinShells_poles()
#gmsh.fltk.run() constructed correctly
if enforce_TSA_mapping_collar:
self.occ.synchronize()
_embed_points_to_collar_curve() ## both coils and wedges
self.occ.remove([(0, center)]) # remove center point
self.occ.synchronize()
def constructThinShells(self, with_wedges):
# default
ins_th = self.md.geometries.thin_shells.ins_thickness
mid_pole_ts = self.md.geometries.thin_shells.mid_poles
mid_winding_ts = self.md.geometries.thin_shells.mid_windings
mid_turn_ts = self.md.geometries.thin_shells.mid_turn_blocks
# not default
mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht
mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
# Create mid-pole and mid-turn thin shells
for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
for layer_nr, layer in coil.layers.items():
for nr, blk_order in enumerate(layer):
block = self.geom.coil.coils[coil_nr].poles[blk_order.pole].layers[
layer_nr].windings[blk_order.winding].blocks[blk_order.block]
ht_list = list(self.md.geometries.coil.coils[coil_nr].poles[blk_order.pole].layers[
layer_nr].windings[blk_order.winding].blocks[blk_order.block].half_turns.areas.keys())
# Create mid-pole and mid-winding thin shells
blk_index_next = nr + 1 if nr + 1 < len(layer) else 0
block_order_next = layer[blk_index_next]
block_next = self.geom.coil.coils[coil_nr].poles[block_order_next.pole].layers[
layer_nr].windings[block_order_next.winding].blocks[block_order_next.block]
ht_list_next = list(self.md.geometries.coil.coils[coil_nr].poles[block_order_next.pole].layers[
layer_nr].windings[block_order_next.winding].blocks[block_order_next.block].half_turns.areas.keys())
ht_last = int(ht_list[-1])
ht_next_first = int(ht_list_next[0])
iH = [block.half_turns[ht_last].corners.bare.iH.x, block.half_turns[ht_last].corners.bare.iH.y]
iL = [block_next.half_turns[ht_next_first].corners.bare.iL.x, block_next.half_turns[ht_next_first].corners.bare.iL.y]
oH = [block.half_turns[ht_last].corners.bare.oH.x, block.half_turns[ht_last].corners.bare.oH.y]
oL = [block_next.half_turns[ht_next_first].corners.bare.oL.x, block_next.half_turns[ht_next_first].corners.bare.oL.y]
ts_name = str(blk_order.block) + '_' + str(block_order_next.block)
for ts, th, condition in zip([mid_pole_ts, mid_winding_ts], [ins_th.mid_pole, ins_th.mid_winding],
# ['_ly' + str(layer_nr), '_wd' + str(blk_order.winding) + '_wd' + str(block_order_next.winding)],
[self.geom.coil.coils[coil_nr].type == 'cos-theta' and block_order_next.pole != blk_order.pole,
(not with_wedges or not self.geom.wedges) and self.geom.coil.coils[coil_nr].type in
['cos-theta', 'common-block-coil'] and block_order_next.pole == blk_order.pole and block_order_next.winding != blk_order.winding]):
if condition:
ts[ts_name] = dM.Region()
ts[ts_name].points['i'] = self.occ.addPoint((iH[0] + iL[0]) / 2, (iH[1] + iL[1]) / 2, 0)
ts[ts_name].points['o'] = self.occ.addPoint((oH[0] + oL[0]) / 2, (oH[1] + oL[1]) / 2, 0)
ts[ts_name].lines[str(ht_last) + '_' + str(ht_next_first)] =\
self.occ.addLine(ts[ts_name].points['i'], ts[ts_name].points['o'])
# Get insulation thickness
th[ts_name] = Func.sig_dig((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)
# if 'cl' + str(coil_nr) + th_name not in th:
# th['cl' + str(coil_nr) + th_name] = float((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)
# Create mid-turn thin shells
mid_turn_ts[str(blk_order.block)] = dM.Region()
ts = mid_turn_ts[str(blk_order.block)]
for nr_ht, ht in enumerate(ht_list[:-1]):
line_name = ht + '_' + ht_list[nr_ht + 1]
current_ht = block.half_turns[int(ht)].corners.bare
next_ht = block.half_turns[int(ht_list[nr_ht + 1])].corners.bare
mid_inner = [(current_ht.iH.x + next_ht.iL.x) / 2, (current_ht.iH.y + next_ht.iL.y) / 2]
mid_outer = [(current_ht.oH.x + next_ht.oL.x) / 2, (current_ht.oH.y + next_ht.oL.y) / 2]
mid_length = Func.points_distance(mid_inner, mid_outer)
mid_line = Func.line_through_two_points(mid_inner, mid_outer)
points = {'inner': list, 'outer': list}
for case, current_h, current_l, next_h, next_l, mid_point in zip(
['inner', 'outer'], [current_ht.iH, current_ht.oH], [current_ht.iL, current_ht.oL],
[next_ht.iH, next_ht.oH], [next_ht.iL, next_ht.oL], [mid_outer, mid_inner]):
current_line = Func.line_through_two_points([current_h.x, current_h.y], [current_l.x, current_l.y])
next_line = Func.line_through_two_points([next_h.x, next_h.y], [next_l.x, next_l.y])
current_intersect = Func.intersection_between_two_lines(mid_line, current_line)
next_intersect = Func.intersection_between_two_lines(mid_line, next_line)
points[case] = current_intersect if Func.points_distance(
current_intersect, mid_point) < mid_length else next_intersect
ts.points[line_name + '_i'] = self.occ.addPoint(points['inner'][0], points['inner'][1], 0)
ts.points[line_name + '_o'] = self.occ.addPoint(points['outer'][0], points['outer'][1], 0)
ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])
# Create mid-layer thin shells
block_coil_mid_pole_list = [str(blks[0].block) + '_' + str(blks[1].block) for coil_nr, coil in self.block_coil_mid_pole_blks.items() for blks in coil]
for ts_name, ts in mid_layer_ts.items():
# Order mid-layer thin shell points according to their angle with respect to the x-axis to generate lines
blk1, blk2 = ts_name.split('_')
max_angle = max(ts.point_angles.values())
max_diff = max_angle - min(ts.point_angles.values())
if max_diff > np.pi:
for pnt_name, angle in ts.point_angles.items():
if angle < max_diff / 2:
ts.point_angles[pnt_name] = angle + max_angle
ordered_pnts = [[pnt_name, ts.point_angles[pnt_name], pnt] for pnt_name, pnt in ts.mid_layers.points.items()]
ordered_pnts.sort(key=lambda x: x[1])
for nr, pnt in enumerate(ordered_pnts[:-1]):
pnt_current = pnt[0]
pnt_next = ordered_pnts[nr + 1][0]
if ((pnt_current[-1] == 'l' and pnt_next[-1] == 'h' and ts_name not in block_coil_mid_pole_list) or # cos-theta
(ts_name in block_coil_mid_pole_list and
((pnt_current[-1] == pnt_next[-1] == 'h' and block_coil_mid_pole_list.index(ts_name) == 0) or # assumes a dipole block-coil
(pnt_current[-1] == pnt_next[-1] == 'l' and block_coil_mid_pole_list.index(ts_name) == 1) or # assumes a dipole block-coil
(pnt_current[:-1] == pnt_next[:-1])))):
if pnt_current[:-1] == pnt_next[:-1]:
relevant_blk = blk2 if int(pnt_current[:-1]) in ts.half_turn_lists[blk1] else blk1
if nr > 0:
iter_nr = nr - 1
while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr -= 1
line_name = ordered_pnts[iter_nr][0][:-1] + '_' + pnt_current[:-1]
else:
if len(ordered_pnts) == 2: # todo: get right ht from relevant_blk for 1-ht blocks
line_name = pnt_current[:-1] + '_' + str(ts.half_turn_lists[relevant_blk][0])
else:
iter_nr = nr + 1
while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr += 1
line_name = pnt_current[:-1] + '_' + ordered_pnts[iter_nr][0][:-1]
else:
line_name = pnt_current[:-1] + '_' + pnt_next[:-1]
# TODO look into why this exception handling is needed for SMC magnet with ESC coils - the following meshing stage does not work
try:
tag = self.occ.addLine(pnt[2], ordered_pnts[nr + 1][2])
ts.mid_layers.lines[line_name] = tag
except Exception as e:
ts.mid_layers.lines[line_name] = tag # this will be the last tag, i.e. from previously created line
x1, y1, z1 = gmsh.model.occ.getBoundingBox(1, pnt[2])[:3]
x2, y2, z2 = gmsh.model.occ.getBoundingBox(1, ordered_pnts[nr + 1][2])[:3]
distance = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)
logger.info(f"{e} {line_name} between point tags {pnt[2], ordered_pnts[nr + 1][2]} and distance between points {distance}")
if ts_name in mid_layer_ts_aux:
aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[0]
other_pnt = ordered_pnts[0 if aux_pnt[-1] == 'l' else -1]
other_pnt_coord = gmsh.model.getValue(0, other_pnt[2], [])[:2] # needs to be a new point
mid_layer_ts_aux[ts_name].points[other_pnt[0]] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)
line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]
try:
mid_layer_ts_aux[ts_name].lines[line_name] = \
self.occ.addLine(mid_layer_ts_aux[ts_name].points[aux_pnt], mid_layer_ts_aux[ts_name].points[other_pnt[0]])
except:
mid_layer_ts_aux[ts_name].lines.pop(line_name)
# Create wedge-to-block and block-to-wedge lines
for wdg_ts in [self.md.geometries.thin_shells.mid_layers_wdg_to_ht, self.md.geometries.thin_shells.mid_layers_ht_to_wdg]:
for ts_name, ts in wdg_ts.items():
pnt_list = list(ts.points.keys())
for nr, pnt in enumerate(pnt_list[:-1]):
if pnt[-1] == 'l' and pnt_list[nr + 1][-1] == 'h':
ts.lines[pnt[:-1] + '_' + pnt_list[nr + 1][:-1]] = self.occ.addLine(ts.points[pnt], ts.points[pnt_list[nr + 1]])
if ts_name in mid_layer_ts_aux:
aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[
1 if list(mid_layer_ts_aux[ts_name].points.keys()).index('center') == 0 else 0]
other_pnt = pnt_list[0 if aux_pnt[-1] == 'l' else -1]
other_pnt_coord = gmsh.model.getValue(0, ts.points[other_pnt], [])[:2] # needs to be a new point
mid_layer_ts_aux[ts_name].points[other_pnt] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)
line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]
mid_layer_ts_aux[ts_name].lines[line_name] = self.occ.addCircleArc(
mid_layer_ts_aux[ts_name].points[aux_pnt], mid_layer_ts_aux[ts_name].points['center'], mid_layer_ts_aux[ts_name].points[other_pnt])
# Create wedge-to-wedge lines
for ts_nr, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():
ts.lines[ts_nr] = self.occ.addCircleArc(ts.points['beg'], ts.points['center'], ts.points[list(ts.points.keys())[-1]])
# Create mid wedge-turn lines
mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn
for ts_nr, ts in mid_turn_ts.items():
line_name = list(ts.points.keys())[0][:-2]
ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])
# Get insulation thickness
for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
geom_coil = self.geom.coil.coils[coil_nr]
# Get block-coil mid-pole thickness
if coil_nr in self.block_coil_mid_pole_blks:
for blk_orders in self.block_coil_mid_pole_blks[coil_nr]:
block_y = geom_coil.poles[blk_orders[0].pole].layers[1].windings[blk_orders[0].winding].blocks[blk_orders[0].block].block_corners.iH.y
block_next_y = geom_coil.poles[blk_orders[1].pole].layers[1].windings[blk_orders[1].winding].blocks[blk_orders[1].block].block_corners.iH.y
ins_th.mid_layer[str(blk_orders[0].block) + '_' + str(blk_orders[1].block)] = Func.sig_dig(abs(block_y - block_next_y))
# Get mid-layer thickness by intersecting the line passing through i-o of the ht of one side with the line passing through l-h of the ht of the opposite side
for layer_nr, layer in coil.layers.items():
for blk_order in layer:
for ts_name, ts in mid_layer_ts.items():
blk1, blk2 = ts_name.split('_')
if blk1 == str(blk_order.block) and ts_name not in block_coil_mid_pole_list:
block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]
if layer_nr < len(coil.layers):
for blk_order_next in coil.layers[layer_nr + 1]:
if blk_order_next.block == int(blk2):
block_next = geom_coil.poles[blk_order_next.pole].layers[layer_nr + 1].windings[blk_order_next.winding].blocks[int(blk2)]
break
else:
for blk_order_next in self.md.geometries.coil.anticlockwise_order.coils[coil_nr + 1].layers[1]:
if blk_order_next.block == int(blk2):
block_next = self.geom.coil.coils[coil_nr + 1].poles[blk_order_next.pole].layers[1].windings[blk_order_next.winding].blocks[int(blk2)]
break
distances = []
lines = list(ts.mid_layers.lines.keys())
for line_name in [lines[0], lines[-1]]:
ht_1, ht_2 = int(line_name[:line_name.index('_')]), int(line_name[line_name.index('_') + 1:])
ht_char = {'low_p': ht_1, 'high_p': ht_2,
'current': ht_1 if ht_1 in ts.half_turn_lists[blk1] else ht_2,
'next': ht_2 if ht_1 in ts.half_turn_lists[blk1] else ht_1}
hts = {'current': block.half_turns[ht_char['current']].corners.bare,
'next': block_next.half_turns[ht_char['next']].corners.bare}
hts_p = {'low_p': [hts['current'].oL, hts['current'].iL] if ht_char['low_p'] == ht_char['current'] else [hts['next'].iL, hts['next'].oL],
'high_p': [hts['current'].oH, hts['current'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['next'].iH, hts['next'].oH],
'low_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['low_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH],
'high_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH]}
low_line = Func.line_through_two_points([hts_p['low_p'][0].x, hts_p['low_p'][0].y],
[hts_p['low_p'][1].x, hts_p['low_p'][1].y])
high_line = Func.line_through_two_points([hts_p['high_p'][0].x, hts_p['high_p'][0].y],
[hts_p['high_p'][1].x, hts_p['high_p'][1].y])
distances.extend([Func.points_distance([hts_p['low_p'][0].x, hts_p['low_p'][0].y], Func.intersection_between_two_lines(
low_line, Func.line_through_two_points([hts_p['low_p_opp'][0].x, hts_p['low_p_opp'][0].y], [hts_p['low_p_opp'][1].x, hts_p['low_p_opp'][1].y]))),
Func.points_distance([hts_p['high_p'][0].x, hts_p['high_p'][0].y], Func.intersection_between_two_lines(
high_line, Func.line_through_two_points([hts_p['high_p_opp'][0].x, hts_p['high_p_opp'][0].y], [hts_p['high_p_opp'][1].x, hts_p['high_p_opp'][1].y])))])
ins_th.mid_layer[ts_name] = Func.sig_dig(min(distances))
for ts_type, wdg_ts in enumerate([self.md.geometries.thin_shells.mid_layers_wdg_to_ht, self.md.geometries.thin_shells.mid_layers_ht_to_wdg]):
for ts_name, ts in wdg_ts.items():
wdg, blk = ts_name.split('_')
if blk == str(blk_order.block):
block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]
wedge = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + (1 if ts_type == 1 else -1)].wedges[int(wdg[1:])]
pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]
pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]
pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]
pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]
low_line = Func.line_through_two_points(pnt_il, pnt_ol)
high_line = Func.line_through_two_points(pnt_ih, pnt_oh)
el1_l, el2_l = list(ts.lines.keys())[0].split('_')
ht_l = block.half_turns[int(el1_l) if el2_l == wdg else int(el2_l)].corners.bare
el1_h, el2_h = list(ts.lines.keys())[-1].split('_')
ht_h = block.half_turns[int(el1_h) if el2_h == wdg else int(el2_h)].corners.bare
opp_line_l = Func.line_through_two_points([ht_l.iL.x, ht_l.iL.y], [ht_l.iH.x, ht_l.iH.y]) if ts_type == 0\
else Func.line_through_two_points([ht_l.oL.x, ht_l.oL.y], [ht_l.oH.x, ht_l.oH.y])
opp_line_h = Func.line_through_two_points([ht_h.iL.x, ht_h.iL.y], [ht_h.iH.x, ht_h.iH.y]) if ts_type == 0 \
else Func.line_through_two_points([ht_h.oL.x, ht_h.oL.y], [ht_h.oH.x, ht_h.oH.y])
ins_th.mid_layer[ts_name] = Func.sig_dig(
(Func.points_distance(pnt_ol if ts_type == 0 else pnt_il, Func.intersection_between_two_lines(low_line, opp_line_l)) +
Func.points_distance(pnt_oh if ts_type == 0 else pnt_ih, Func.intersection_between_two_lines(high_line, opp_line_h))) / 2)
for coil_nr, coil in self.md.geometries.wedges.coils.items():
# Get mid-layer thickness by intersecting the line passing through i-o of the wdg of one side with the line passing through l-h of the wdg of the opposite side
for layer_nr, layer in coil.layers.items():
for wedge_nr, wedge in layer.wedges.items():
for ts_name, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():
wdg1, wdg2 = ts_name[1:ts_name.index('_')], ts_name[ts_name.index('_') + 2:]
if wdg1 == str(wedge_nr):
wedge_next = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + 1].wedges[int(wdg2)]
# pnt_il_next = gmsh.model.getValue(0, wedge_next.points['il'], [])[:2]
# pnt_ih_next = gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2]
pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]
pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]
pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]
pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]
low_line = Func.line_through_two_points(pnt_il, pnt_ol)
high_line = Func.line_through_two_points(pnt_ih, pnt_oh)
opp_line = Func.line_through_two_points(gmsh.model.getValue(0, wedge_next.points['il'], [])[:2],
gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2])
ins_th.mid_layer[ts_name] = Func.sig_dig(
(Func.points_distance(pnt_ol, Func.intersection_between_two_lines(low_line, opp_line)) +
Func.points_distance(pnt_oh, Func.intersection_between_two_lines(high_line, opp_line))) / 2)
def buildDomains(self, run_type, symmetry):
"""
Generates plane surfaces from the curve loops
"""
iron = self.geom.iron
gm = self.md.geometries
geometry_setting = self.data.magnet.geometry.electromagnetics if run_type == 'EM' \
else self.data.magnet.geometry.thermal
with_wedges = geometry_setting.with_wedges
inv_nc = {v: k for k, v in self.nc.items()} #invert naming convention
for a in geometry_setting.areas: # a in ['iron_yoke', 'collar', ...]:
for quadrant, qq in getattr(gm, a).quadrants.items():
for area_name, area in qq.areas.items():
identifier = next((k for k in self.inv_nc.keys() if (k in area_name[2:])), None)#re.sub(r'\d+', '',area_name[2:])
if a == inv_nc.get(identifier, None): # ensure it is part of the iron yoke or collar (iron, collar)
build = True
loops = [area.loop]
for hole_key, hole in iron.hyper_holes.items():
if area_name == hole.areas[1]:
loops.append(qq.areas[hole.areas[0]].loop)
elif area_name == hole.areas[0]: #skip holes
area.surface = self.occ.addPlaneSurface(loops) # also build the holes. An existing curveloop without area is very annoying
build = False
if build:
area.surface = self.occ.addPlaneSurface(loops)
getattr(self.md.domains.groups_entities, a)[iron.hyper_areas[area_name].material].append(area.surface) ## save the material
# Build coil domains
for coil_nr, coil in gm.coil.coils.items():
for pole_nr, pole in coil.poles.items():
for layer_nr, layer in pole.layers.items():
for winding_nr, winding in layer.windings.items():
for block_key, block in winding.blocks.items():
for area_name, area in block.half_turns.areas.items():
area.surface = self.occ.addPlaneSurface([area.loop])
# Build wedges domains
if with_wedges:
for coil_nr, coil in gm.wedges.coils.items():
for layer_nr, layer in coil.layers.items():
for wedge_nr, wedge in layer.wedges.items():
wedge.areas[str(wedge_nr)].surface = self.occ.addPlaneSurface([wedge.areas[str(wedge_nr)].loop])
# Build insulation domains
if run_type == 'TH' and not geometry_setting.use_TSA:
for coil_nr, coil in gm.insulation.coils.items():
for group_nr, group in coil.group.items():
holes = []
for blk in group.blocks:
holes.extend([ht.loop for ht_nr, ht in gm.coil.coils[
coil_nr].poles[blk[0]].layers[blk[1]].windings[blk[2]].blocks[blk[3]].half_turns.areas.items()])
for wdg in group.wedges:
holes.extend([wedge.loop for wedge_nr, wedge in gm.wedges.coils[
coil_nr].layers[wdg[0]].wedges[wdg[1]].areas.items()])
if len(group.ins.areas) == 1:
for area_name, area in group.ins.areas.items():
area.surface = self.occ.addPlaneSurface([area.loop] + holes)
else:
for area_name, area in group.ins.areas.items():
if area_name.isdigit():
area.surface = self.occ.addPlaneSurface([area.loop] + holes + [group.ins.areas['inner_loop'].loop])
# Create and build air far field
if run_type == 'EM':
if 'iron_yoke' in geometry_setting.areas:
for i in iron.key_points:
gm.iron_yoke.max_radius = max(gm.iron_yoke.max_radius, max(iron.key_points[i].x, iron.key_points[i].y)) # this also contains other regions, e.g. collar but this has no effect
greatest_radius = gm.iron_yoke.max_radius
else: # no iron yoke data available
for coil_nr, coil in self.geom.coil.coils.items():
for pole_nr, pole in coil.poles.items():
first_winding = list(pole.layers[len(pole.layers)].windings.keys())[0]
first_block = list(pole.layers[len(pole.layers)].windings[first_winding].blocks)[0]
gm.coil.max_radius = max(abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.x),
abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.y),
gm.coil.max_radius)
greatest_radius = gm.coil.max_radius
radius_in = greatest_radius * (2.5 if 'iron_yoke' in geometry_setting.areas else 6)
radius_out = greatest_radius * (3.2 if 'iron_yoke' in geometry_setting.areas else 8)
air_inf_center_x, air_inf_center_y = 0, 0
for coil_nr, coil in self.md.geometries.coil.coils.items():
air_inf_center_x += coil.bore_center.x
air_inf_center_y += coil.bore_center.y
gm.air.points['bore_center' + str(coil_nr)] = self.occ.addPoint(coil.bore_center.x, coil.bore_center.y, 0.)
air_inf_center = [air_inf_center_x / len(self.md.geometries.coil.coils), air_inf_center_y / len(self.md.geometries.coil.coils)]
if symmetry == 'none':
gm.air_inf.lines['inner'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_in)
gm.air_inf.lines['outer'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_out)
gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['inner']]))
gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['outer']]))
gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop, gm.air_inf.areas['inner'].loop])
else:
pnt1 = [1, 0] if symmetry in ['xy', 'x'] else [0, -1]
pnt2 = [0, 1] if symmetry in ['xy', 'y'] else [-1, 0]
gm.air.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_in, pnt1[1] * radius_in, 0)
gm.air.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_in, pnt2[1] * radius_in, 0)
gm.air_inf.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_out, pnt1[1] * radius_out, 0)
gm.air_inf.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_out, pnt2[1] * radius_out, 0)
gm.air.lines['ln1'] = self.occ.addLine(gm.air.points['pnt1'], gm.air_inf.points['pnt1'])
gm.air.lines['ln2'] = self.occ.addLine(gm.air.points['pnt2'], gm.air_inf.points['pnt2'])
if not self.data.magnet.geometry.electromagnetics.with_iron_yoke:
gm.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)
gm.air_inf.lines['inner'] = self.occ.addCircleArc(gm.air.points['pnt2'], gm.air_inf.points['center'], gm.air.points['pnt1'])
gm.air_inf.lines['outer'] = self.occ.addCircleArc(gm.air_inf.points['pnt2'], gm.air_inf.points['center'], gm.air_inf.points['pnt1'])
if symmetry in ['xy', 'x']:
gm.air.lines['x_p'] = self.occ.addLine(self.md.geometries.air_inf.points['center'] if 'solenoid' in self.geom.coil.coils[1].type else
gm.iron.quadrants[1].points[self.symmetric_bnds['x_p']['pnts'][-1][0]], gm.air.points['pnt1'])
self.symmetric_loop_lines['x'].append(gm.air.lines['x_p'])
else: # y
gm.air.lines['y_n'] = self.occ.addLine(gm.iron.quadrants[4].points[self.symmetric_bnds['y_n']['pnts'][-1][0]], gm.air.points['pnt1'])
self.symmetric_loop_lines['y'].append(gm.air.lines['y_n'])
if symmetry in ['xy', 'y']:
gm.air.lines['y_p'] = self.occ.addLine(gm.iron.quadrants[1].points[self.symmetric_bnds['y_p']['pnts'][-1][0]], gm.air.points['pnt2'])
self.symmetric_loop_lines['y'].insert(0, gm.air.lines['y_p'])
else: # x
gm.air.lines['x_n'] = self.occ.addLine(self.md.geometries.air_inf.points['center'] if 'solenoid' in self.geom.coil.coils[1].type else
gm.iron.quadrants[2].points[self.symmetric_bnds['x_n']['pnts'][-1][0]], gm.air.points['pnt2'])
self.symmetric_loop_lines['x'].insert(0, gm.air.lines['x_n'])
inner_lines = self.symmetric_loop_lines['x'] + [gm.air_inf.lines['inner']] + self.symmetric_loop_lines['y']\
if symmetry == 'xy' else self.symmetric_loop_lines[symmetry] + [gm.air_inf.lines['inner']]
gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop(inner_lines))
gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop(
[gm.air.lines['ln1'], gm.air_inf.lines['outer'], gm.air.lines['ln2'], gm.air_inf.lines['inner']]))
gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop])
# self.md.domains.groups_entities.air_inf = [gm.air_inf.areas['outer'].surface]
gm.air_inf.areas['inner'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['inner'].loop])
self.occ.synchronize()
#self.gu.launch_interactive_GUI()
def fragment(self):
"""
Fragment and group air domains
"""
# Collect surfaces to be subtracted by background air
holes = []
# iron yoke and collar
group_keys = self.nc.keys()
for key in group_keys:
group = getattr(self.md.domains.groups_entities, key)
for _, surfaces in group.items():
holes.extend([(2, s) for s in surfaces])
# Coils
for coil_nr, coil in self.md.geometries.coil.coils.items():
for pole_nr, pole in coil.poles.items():
for layer_nr, layer in pole.layers.items():
for winding_nr, winding in layer.windings.items():
for block_key, block in winding.blocks.items():
for area_name, area in block.half_turns.areas.items():
holes.append((2, area.surface))
# Wedges
for coil_nr, coil in self.md.geometries.wedges.coils.items():
for layer_nr, layer in coil.layers.items():
for wedge_nr, wedge in layer.wedges.items():
for area_name, area in wedge.areas.items():
holes.append((2, area.surface))
# Insulation
# if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
# for coil_nr, coil in self.md.geometries.insulation.coils.items():
# for group_nr, group in coil.group.items():
# for area_name, area in group.ins.areas.items():
# holes.append((2, area.surface))
# Fragment
fragmented = self.occ.fragment([(2, self.md.geometries.air_inf.areas['inner'].surface)], holes)[1]
self.occ.synchronize()
self.md.domains.groups_entities.air = []
existing_domains = [e[0][1] for e in fragmented[1:]]
for e in fragmented[0]:
if e[1] not in existing_domains:
self.md.domains.groups_entities.air.append(e[1])
def updateTags(self, run_type, symmetry):
# Update half turn line tags
for coil_nr, coil in self.md.geometries.coil.coils.items():
for pole_nr, pole in coil.poles.items():
for layer_nr, layer in pole.layers.items():
for winding_nr, winding in layer.windings.items():
for block_key, block in winding.blocks.items():
hts = block.half_turns
# Get half turn ID numbers
area_list = list(hts.areas.keys())
for nr, ht_nr in enumerate(area_list):
first_tag = int(min(gmsh.model.getAdjacencies(2, hts.areas[ht_nr].surface)[1]))
hts.lines[ht_nr + 'i'] = first_tag
hts.lines[ht_nr + 'l'] = first_tag + 1
hts.lines[ht_nr + 'o'] = first_tag + 2
hts.lines[ht_nr + 'h'] = first_tag + 3
# Update collar tags
if run_type == "TH" and 'collar' in self.data.magnet.geometry.thermal.areas:
for quad, old_tags in self.md.geometries.collar.quadrants.items():
self.md.geometries.collar.inner_boundary_tags[quad] = [] # reset the inner boundary tags
new_tags = []
for name, area in self.md.geometries.collar.quadrants[quad].areas.items(): # arcol contains the boundaries of the holes too
if not re.match(r"^ar.h", name):
# the issue is that you don't know which line is which, e.g. which is the inner collar line
new_tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])
for k, name in enumerate(self.md.geometries.collar.quadrants[quad].lines.keys()):
self.md.geometries.collar.quadrants[quad].lines[name] = new_tags[k]
# Update inner collar tags
collar_lines = [self.md.geometries.collar.quadrants[quad].lines[name] for name in self.md.geometries.collar.quadrants[quad].lines.keys()]
closest_dist = 1000.
closest_lines = []
max_dist = 0.0
max_lines = []
# We assume that the middle of the curve of the collar is the closest one to the centre (0,0)
for tag in collar_lines:
##x, y, _ = gmsh.model.getValue(1, tag, [0.0]) # pick one point on the line
curve_type = gmsh.model.getType(1, tag)
if curve_type == 'Line': # find the middle of the line
tag1, tag2 = gmsh.model.getAdjacencies(1, tag)[1]
x1, y1, z1 = gmsh.model.getValue(0, tag1, [])
x2, y2, z2 = gmsh.model.getValue(0, tag2, [])
x = 0.5 * (x1 + x2)
y = 0.5 * (y1 + y2)
# take the average of the end points
dist = np.sqrt(x ** 2 + y ** 2)
elif curve_type == "Circle": # use any point on the circle (same distance because concentric with origin)
x, y, z = gmsh.model.getValue(1, tag, [0.5])
dist = np.sqrt(x ** 2 + y ** 2)
if dist < closest_dist+1e-10:
if dist < closest_dist-1e-10: # clear if new min is found
closest_dist = dist
closest_lines = []
closest_lines.append(tag)
if dist > max_dist-1e-10:
if dist > max_dist+1e-10: # clear if new max is found
max_dist = dist
max_lines = []
max_lines.append(tag)
self.md.geometries.collar.inner_boundary_tags[quad] = closest_lines
self.md.geometries.collar.outer_boundary_tags[quad] = max_lines
# outer collar tags, does not work because geom.iron has the old tags and this is not accurate anymore
"""
outer = [old_tags.lines[name] for name in old_tags.lines.keys() if
self.geom.iron.hyper_lines[name].type == 'line']
logger.info("outer tags", outer)
self.md.geometries.collar.outer_boundary_tags[quad] = outer
"""
# concerning the TSL, it seems impossible to get the tags of the lines, as they are not adjacent to a surface nor (yet) grouped in a physical group
# only way to ensure equal tags is by creating them in the same way as gmsh orders the lines. (e.g. all at the start or all at the end would work)
# we cannot swap order... soo lets hope that just a shift in numbers is sufficient
if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA_new and self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping:
shift = len(self.md.geometries.collar.inner_boundary_tags[1])*4 -4
if 'poles' in self.data.magnet.geometry.thermal.areas:
shift -= (1*4) # we shifted too much
## update TSA collar lines
for _, ts in self.md.geometries.thin_shells.collar_layers.items():
ts.lines['1'] += shift
for _, ts in self.md.geometries.thin_shells.pole_layers.items():
ts.lines['1'] += shift
atts = ['mid_layers_ht_to_ht', 'mid_layers_wdg_to_ht', 'mid_layers_ht_to_wdg', 'mid_layers_wdg_to_wdg', 'mid_poles', 'mid_windings', 'mid_turn_blocks' , 'mid_wedge_turn']
for at in atts:
for _, ts_region in getattr(self.md.geometries.thin_shells, at).items():
try: ts_region = ts_region.mid_layers
except AttributeError: pass
for key in ts_region.lines.keys():
ts_region.lines[key] += shift
# Update coil cooling tags
### no consistent way to get the tags of the cooling lines, so we assume that they are ordered in the same way as gmsh orders the lines
if self.data.magnet.solve.thermal.collar_cooling.enabled:
tags =[]
## if we want all cooling holes
if self.data.magnet.solve.thermal.collar_cooling.which == 'all':
for quad, region in self.md.geometries.collar.quadrants.items():
for name, area in region.areas.items():
if re.match(r"^ar.h", name):
tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])
else:
nr_applied_cooling = self.data.magnet.solve.thermal.collar_cooling.which
nr = 1
for _, quad_data in self.md.geometries.collar.quadrants.items():
for name, area in quad_data.areas.items():
if re.match(r"^ar.h", name):
if nr in nr_applied_cooling:
tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])
nr += 1
self.md.geometries.collar.cooling_tags = tags
# Update insulation line tags
if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
pass # todo
# Update wedge line tags
for coil_nr, coil in self.md.geometries.wedges.coils.items():
for layer_nr, layer in coil.layers.items():
for wedge_nr, wedge in layer.wedges.items():
lines_tags = list(gmsh.model.getAdjacencies(2, wedge.areas[str(wedge_nr)].surface)[1])
lines_tags.sort(key=lambda x: x)
wedge.lines['i'] = int(lines_tags[0])
wedge.lines['l'] = int(lines_tags[1])
wedge.lines['o'] = int(lines_tags[2])
wedge.lines['h'] = int(lines_tags[3])
if run_type == 'EM':
def _get_bnd_lines():
return [pair[1] for pair in self.occ.getEntitiesInBoundingBox(corner_min[0], corner_min[1], corner_min[2],
corner_max[0], corner_max[1], corner_max[2], dim=1)]
tol = 1e-6
# Update tags of air and air_inf arcs and their points
lines_tags = gmsh.model.getAdjacencies(2, self.md.geometries.air_inf.areas['outer'].surface)[1]
self.md.geometries.air_inf.lines['outer'] = int(lines_tags[0 if symmetry == 'none' else 1])
self.md.geometries.air_inf.lines['inner'] = int(lines_tags[1 if symmetry == 'none' else 3])
if symmetry == 'none': # todo: check if this holds for symmetric models too
for coil_nr, coil in self.md.geometries.coil.coils.items():
self.md.geometries.air.points['bore_center' + str(coil_nr)] += 2
else:
pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['outer'])[1])
indexes = [0, 1, 0] if 'x' in symmetry else [1, 0, 1]
pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] >\
gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]
self.md.geometries.air_inf.points['pnt1'] = int(pnt_tags[pnts[0]])
self.md.geometries.air_inf.points['pnt2'] = int(pnt_tags[pnts[1]])
pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['inner'])[1])
pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] > \
gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]
self.md.geometries.air.points['pnt1'] = int(pnt_tags[pnts[0]])
self.md.geometries.air.points['pnt2'] = int(pnt_tags[pnts[1]])
for coil_nr, coil in self.md.geometries.coil.coils.items():
self.md.geometries.air.points['bore_center' + str(coil_nr)] =(
self.occ.getEntitiesInBoundingBox(-tol + coil.bore_center.x, -tol + coil.bore_center.y, -tol,
tol + coil.bore_center.x, tol + coil.bore_center.y, tol, dim=0))[0][1]
# Group symmetry boundary lines per type
if symmetry == 'xy':
corner_min = [-tol, -tol, -tol]
corner_max = [gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0] + tol, tol, tol]
self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()
corner_max = [tol, gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1] + tol, tol]
self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()
elif symmetry == 'x':
x_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0]
corner_min = [- x_coord - tol, -tol, -tol]
corner_max = [x_coord + tol, tol, tol]
self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()
elif symmetry == 'y':
y_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1]
corner_min = [-tol, - y_coord - tol, -tol]
corner_max = [tol, y_coord + tol, tol]
self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()
def move_keypoints(self, keypoints, displacement, keypoint_names=None):
if keypoint_names is None:
keypoint_names = []
for name, hole in self.geom.iron.hyper_areas.items():
if not 'ch' in name: # ch -> collar hole
continue
line_names = hole.lines
keypoint_names.append(list(set([getattr(self.geom.iron.hyper_lines[line], kp_name)
for line in line_names for kp_name in ['kp1', 'kp2', 'kp3']])))
if type(displacement) == list:
list_displacement = displacement
elif str(displacement) == "0":
list_displacement = [[0.0, 0.0], [0.0, 0.0]]
elif str(displacement) == "1":
list_displacement = [[0.004, -0.015], [0.03, -0.025]]
elif str(displacement) == "2":
list_displacement = [[0.004, -0.015], [0.0, 0.0]]
elif str(displacement) == "3":
list_displacement = [[-0.035, 0.045], [-0.004, -0.0015]]
else:
raise ValueError("displacement_type not recognized")
for i, hole in enumerate(keypoint_names):
for name in hole:
if name is None:
continue
keypoints[name].x += list_displacement[i][0]
keypoints[name].y += list_displacement[i][1]
return keypoints