Coverage for fiqus/geom_generators/GeometryMultipole.py: 79%
1505 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-03-14 02:31 +0100
1import os
2import gmsh
3import numpy as np
4import pandas as pd
5import json
7from fiqus.utils.Utils import GmshUtils
8from fiqus.utils.Utils import FilesAndFolders as Util
9from fiqus.utils.Utils import GeometricFunctions as Func
10from fiqus.data import DataFiQuS as dF
11from fiqus.data import DataMultipole as dM
12from fiqus.data.DataRoxieParser import Corner
13from fiqus.data.DataRoxieParser import Coord
16class Geometry:
17 def __init__(self, data: dF.FDM() = None, geom: dF.FiQuSGeometry() = None,
18 geom_folder: str = None, verbose: bool = False):
19 """
20 Class to generate geometry
21 :param data: FiQuS data model
22 :param geom: ROXIE geometry data
23 :param verbose: If True more information is printed in python console.
24 """
25 self.data: dF.FDM() = data
26 self.geom: dF.RoxieData() = geom.Roxie_Data
27 self.geom_folder = geom_folder
28 self.verbose: bool = verbose
30 self.md = dM.MultipoleData()
32 self.gu = GmshUtils(self.geom_folder, self.verbose)
33 self.gu.initialize(verbosity_Gmsh=self.data.run.verbosity_Gmsh)
34 self.occ = gmsh.model.occ
36 self.model_file = os.path.join(self.geom_folder, self.data.general.magnet_name)
38 self.blk_ins_lines = {} # for meshed insulation
39 self.ins_wire_lines = {} # for meshed insulation
40 self.block_coil_mid_pole_blks = {}
42 if self.data.magnet.geometry.electromagnetics.symmetry != 'none':
43 self.symmetric_loop_lines = {'x': [], 'y': []}
44 self.symmetric_bnds = {'x_p': {'pnts': [], 'line_pnts': []}, 'y_p': {'pnts': [], 'line_pnts': []},
45 'x_n': {'pnts': [], 'line_pnts': []}, 'y_n': {'pnts': [], 'line_pnts': []}}
47 def clear(self):
48 self.md = dM.MultipoleData()
49 self.block_coil_mid_pole_blks = {}
50 gmsh.clear()
52 def ending_step(self, gui: bool = False):
53 if gui:
54 self.gu.launch_interactive_GUI()
55 else:
56 gmsh.clear()
57 gmsh.finalize()
59 def saveHalfTurnCornerPositions(self):
60 self.occ.synchronize()
61 iH, iL, oH, oL, iHr, iLr, oHr, oLr = [], [], [], [], [], [], [], []
62 for po in self.geom.coil.physical_order:
63 block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[
64 po.winding].blocks[po.block]
65 for halfTurn_nr, halfTurn in block.half_turns.items():
66 ht = halfTurn.corners.insulated
67 ht_b = halfTurn.corners.bare
68 iHr.append([ht_b.iH.x, ht_b.iH.y])
69 iLr.append([ht_b.iL.x, ht_b.iL.y])
70 oHr.append([ht_b.oH.x, ht_b.oH.y])
71 oLr.append([ht_b.oL.x, ht_b.oL.y])
72 iH.append([ht.iH.x, ht.iH.y])
73 iL.append([ht.iL.x, ht.iL.y])
74 oH.append([ht.oH.x, ht.oH.y])
75 oL.append([ht.oL.x, ht.oL.y])
76 json.dump({'iH': iH, 'iL': iL, 'oH': oH, 'oL': oL,
77 'iHr': iHr, 'iLr': iLr, 'oHr': oHr, 'oLr': oLr}, open(f"{self.model_file}.crns", 'w'))
79 def saveStrandPositions(self, run_type):
80 symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'
81 ht_nr = 0
82 std_nr = 0
83 parser_x, parser_y, blocks, ht, std, pole_blocks = [], [], [], [], [], []
84 for po in self.geom.coil.physical_order:
85 block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[
86 po.winding].blocks[po.block]
87 if po.pole == 1: pole_blocks.append(po.block)
88 for halfTurn_nr, halfTurn in block.half_turns.items():
89 ht_nr += 1
90 for strand_group_nr, strand_group in halfTurn.strand_groups.items():
91 for strand_nr, strand in strand_group.strand_positions.items():
92 std_nr += 1
93 blocks.append(po.block)
94 ht.append(ht_nr)
95 std.append(std_nr)
96 parser_x.append(strand.x)
97 parser_y.append(strand.y)
98 mirrored = {}
99 condition = {2: [1, -1], 3: [1, 1], 4: [-1, 1]}
100 if symmetry == 'xy': mirroring = {2: [-1, 1], 3: [-1, -1], 4: [1, -1]}
101 elif symmetry == 'x': mirroring = {3: [1, -1], 4: [1, -1]}
102 elif symmetry == 'y': mirroring = {2: [-1, 1], 3: [-1, 1]}
103 else: mirroring = {}
104 if mirroring:
105 df = pd.DataFrame({'parser_x': parser_x, 'parser_y': parser_y}, index=std)
106 for qdr, mrr in mirroring.items():
107 subdf = df[(condition[qdr][0] * df['parser_x'] < 0) & (condition[qdr][1] * df['parser_y'] < 0)]
108 for strand, x, y in zip(subdf.index, subdf['parser_x'], subdf['parser_y']):
109 mirrored[strand] = df[(df['parser_x'] == mrr[0] * x) & (df['parser_y'] == mrr[1] * y)].index.item()
110 json.dump({'x': parser_x, 'y': parser_y, 'block': blocks, 'ht': ht, 'mirrored': mirrored,
111 'pole_1_blocks': pole_blocks, 'poles': len(self.geom.coil.coils[1].poles)},
112 open(f"{self.model_file}_{run_type}.strs", 'w'))
114 def saveBoundaryRepresentationFile(self, run_type):
115 self.occ.synchronize()
116 gmsh.write(f'{self.model_file}_{run_type}.brep')
117 gmsh.clear()
119 def loadBoundaryRepresentationFile(self, run_type):
120 gmsh.option.setString('Geometry.OCCTargetUnit', 'M') # set units to meters
121 gmsh.open(os.path.join(f'{self.model_file}_{run_type}.brep'))
123 def saveAuxiliaryFile(self, run_type):
124 Util.write_data_to_yaml(f'{self.model_file}_{run_type}.aux', self.md.dict())
126 @staticmethod
127 def findMidLayerPoint(bc_current, bc_next, center, mean_rad):
128 mid_layer = [(bc_current.x + bc_next.x) / 2, (bc_current.y + bc_next.y) / 2]
129 mid_rad = Func.points_distance(mid_layer, [center.x, center.y])
130 dist_from_mid = mean_rad - mid_rad
131 angle = Func.arc_angle_between_point_and_abscissa(mid_layer, [center.x, center.y])
132 mid_layer[0] += dist_from_mid * np.cos(angle)
133 mid_layer[1] += dist_from_mid * np.sin(angle)
134 return mid_layer
136 @staticmethod
137 def getMidLayerEndpoints(el_current, el_next, center, mid_layer_arc_pnt=None, coil_type='cos-theta', cable_type='Rutherford', is_for_mid_pole=False):
138 thin_shell_endpoints = {'higher': list, 'lower': list}
139 which_block = {'higher': str, 'lower': str}
140 angles = {'higher': float, 'lower': float}
141 # Check if the element crosses the x axis
142 angles_to_correct = []
143 correction_angle = 0
144 l_curr = Func.arc_angle_between_point_and_abscissa([el_current.iL.x, el_current.iL.y], center)
145 h_curr = Func.arc_angle_between_point_and_abscissa([el_current.iH.x, el_current.iH.y], center)
146 l_next = Func.arc_angle_between_point_and_abscissa([el_next.iL.x, el_next.iL.y], center)
147 h_next = Func.arc_angle_between_point_and_abscissa([el_next.iH.x, el_next.iH.y], center)
148 if abs(l_curr - h_curr) > np.pi:
149 angles_to_correct.append('current')
150 correction_angle = max(1.05 * (2 * np.pi - l_curr), correction_angle)
151 if abs(l_next - h_next) > np.pi:
152 angles_to_correct.append('next')
153 correction_angle = max(1.05 * (2 * np.pi - l_next), correction_angle)
154 for side in thin_shell_endpoints.keys():
155 if mid_layer_arc_pnt:
156 if side == 'higher':
157 mid_lyr_curr, mid_lyr_next = [el_current.oH, el_current.iH], [el_next.oH, el_next.iH]
158 else:
159 mid_lyr_curr, mid_lyr_next = [el_current.oL, el_current.iL], [el_next.oL, el_next.iL]
160 if cable_type in ['Mono', 'Ribbon']:
161 pnts_curr = Func.intersection_between_circle_and_line(
162 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]),
163 [center, mid_layer_arc_pnt])
164 pnt_curr = pnts_curr[0] if Func.points_distance(pnts_curr[0], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) <\
165 Func.points_distance(pnts_curr[1], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) else pnts_curr[1]
166 pnts_next = Func.intersection_between_circle_and_line(
167 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]),
168 [center, mid_layer_arc_pnt])
169 pnt_next = pnts_next[0] if Func.points_distance(pnts_next[0], [mid_lyr_next[0].x, mid_lyr_next[0].y]) <\
170 Func.points_distance(pnts_next[1], [mid_lyr_next[0].x, mid_lyr_next[0].y]) else pnts_next[1]
171 elif cable_type == 'Rutherford':
172 pnt_curr = Func.intersection_between_circle_and_line(
173 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]),
174 [center, mid_layer_arc_pnt], get_only_closest=True)[0]
175 pnt_next = Func.intersection_between_circle_and_line(
176 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]),
177 [center, mid_layer_arc_pnt], get_only_closest=True)[0]
178 else:
179 if cable_type == 'Rutherford':
180 if coil_type == 'common-block-coil':
181 mid_layer_x = (el_current.oH.x + el_next.iH.x) / 2
182 if side == 'higher':
183 pnt_curr, pnt_next = [mid_layer_x, el_current.iH.y], [mid_layer_x, el_next.iH.y]
184 else:
185 pnt_curr, pnt_next = [mid_layer_x, el_current.iL.y], [mid_layer_x, el_next.iL.y]
186 else:
187 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
188 if side == 'higher':
189 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]
190 else:
191 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]
192 elif cable_type in ['Mono', 'Ribbon']:
193 pnt_curr = [(el_current.oH.x + el_next.iH.x) / 2, (el_current.oH.y + el_next.iH.y) / 2] if side == 'higher'\
194 else [(el_current.oL.x + el_next.iL.x) / 2, (el_current.oL.y + el_next.iL.y) / 2]
195 pnt_next = pnt_curr
196 angle_curr = Func.arc_angle_between_point_and_abscissa(pnt_curr, center)
197 angle_next = Func.arc_angle_between_point_and_abscissa(pnt_next, center)
198 if 'current' in angles_to_correct:
199 angle_curr = angle_curr + correction_angle - (2 * np.pi if side == 'lower' else 0)
200 elif 'next' in angles_to_correct:
201 if angle_curr < np.pi / 2: angle_curr += correction_angle
202 elif angle_curr > np.pi * 3 / 2: angle_curr = angle_curr + correction_angle - 2 * np.pi
203 if 'next' in angles_to_correct:
204 angle_next = angle_next + correction_angle - (2 * np.pi if side == 'lower' else 0)
205 elif 'current' in angles_to_correct:
206 if angle_next < np.pi / 2: angle_next += correction_angle
207 elif angle_next > np.pi * 3 / 2: angle_next = angle_next + correction_angle - 2 * np.pi
208 if abs(angle_curr - angle_next) < 1e-6: # todo: check if needed
209 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'
210 elif angle_curr * (-1 if side == 'lower' else 1) < angle_next * (-1 if side == 'lower' else 1):
211 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'
212 else:
213 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_next, angle_next, 'next'
214 if angles['higher'] < angles['lower']: return None
215 else: return thin_shell_endpoints, which_block
217 def constructIronGeometry(self, symmetry):
218 """
219 Generates points, hyper lines, and curve loops for the iron yoke
220 """
221 iron = self.geom.iron
222 if symmetry == 'xy':
223 self.md.geometries.iron.quadrants = {1: dM.Region()}
224 list_bnds = ['x_p', 'y_p']
225 elif symmetry == 'x':
226 self.md.geometries.iron.quadrants = {1: dM.Region(), 2: dM.Region()}
227 list_bnds = ['x_p', 'x_n']
228 elif symmetry == 'y':
229 self.md.geometries.iron.quadrants = {1: dM.Region(), 4: dM.Region()}
230 list_bnds = ['y_p', 'y_n']
231 else:
232 self.md.geometries.iron.quadrants = {1: dM.Region(), 2: dM.Region(), 4: dM.Region(), 3: dM.Region()}
233 list_bnds = []
234 quadrants = self.md.geometries.iron.quadrants
236 lc = 1e-2
237 for point_name, point in iron.key_points.items():
238 if symmetry in ['x', 'xy']:
239 if point.y == 0.:
240 self.symmetric_bnds['x_p']['pnts'].append([point_name, point.x])
241 if symmetry in ['y', 'xy']:
242 if point.x == 0.:
243 self.symmetric_bnds['y_p']['pnts'].append([point_name, point.y])
244 quadrants[1].points[point_name] = self.occ.addPoint(point.x, point.y, 0, lc)
245 if symmetry in ['x', 'none']:
246 if point.x == 0.:
247 quadrants[2].points[point_name] = quadrants[1].points[point_name]
248 else:
249 quadrants[2].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]
250 self.occ.mirror([(0, quadrants[2].points[point_name])], 1, 0, 0, 0)
251 if point.y == 0. and symmetry == 'x':
252 self.symmetric_bnds['x_n']['pnts'].append([point_name, point.x])
253 if symmetry in ['y', 'none']:
254 if point.y == 0.:
255 quadrants[4].points[point_name] = quadrants[1].points[point_name]
256 else:
257 quadrants[4].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]
258 self.occ.mirror([(0, quadrants[4].points[point_name])], 0, 1, 0, 0)
259 if point.x == 0. and symmetry == 'y':
260 self.symmetric_bnds['y_n']['pnts'].append([point_name, point.y])
261 if symmetry == 'none':
262 if point.y == 0.:
263 quadrants[3].points[point_name] = quadrants[2].points[point_name]
264 elif point.x == 0.:
265 quadrants[3].points[point_name] = quadrants[4].points[point_name]
266 else:
267 quadrants[3].points[point_name] = self.occ.copy([(0, quadrants[2].points[point_name])])[0][1]
268 self.occ.mirror([(0, quadrants[3].points[point_name])], 0, 1, 0, 0)
270 mirror_x = [1, -1, -1, 1]
271 mirror_y = [1, 1, -1, -1]
272 symmetric_bnds_order = {'x': [], 'y': []}
273 sym_lines_tags = {'x_p': [], 'y_p': [], 'x_n': [], 'y_n': []}
274 for line_name, line in iron.hyper_lines.items():
275 pt1 = iron.key_points[line.kp1]
276 pt2 = iron.key_points[line.kp2]
277 if line.type == 'line':
278 for quadrant, qq in quadrants.items():
279 if quadrant == 1:
280 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
281 if pt1.y == 0. and pt2.y == 0. and 'x_p' in list_bnds:
282 self.symmetric_bnds['x_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)
283 sym_lines_tags['x_p'].append(qq.lines[line_name])
284 symmetric_bnds_order['x'].append(min(pt1.x, pt2.x))
285 elif pt1.x == 0. and pt2.x == 0. and 'y_p' in list_bnds:
286 self.symmetric_bnds['y_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)
287 sym_lines_tags['y_p'].append(qq.lines[line_name])
288 symmetric_bnds_order['y'].append(min(pt1.y, pt2.y))
289 elif quadrant == 2:
290 if pt1.x == 0. and pt2.x == 0.:
291 qq.lines[line_name] = quadrants[1].lines[line_name]
292 else:
293 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
294 if pt1.y == 0. and pt2.y == 0. and 'x_n' in list_bnds:
295 self.symmetric_bnds['x_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)
296 sym_lines_tags['x_n'].append(qq.lines[line_name])
297 elif quadrant == 4:
298 if pt1.y == 0. and pt2.y == 0.:
299 qq.lines[line_name] = quadrants[1].lines[line_name]
300 else:
301 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
302 if pt1.x == 0. and pt2.x == 0. and 'y_n' in list_bnds:
303 self.symmetric_bnds['y_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)
304 sym_lines_tags['y_n'].append(qq.lines[line_name])
305 else: # 3
306 if pt1.y == 0. and pt2.y == 0.:
307 qq.lines[line_name] = quadrants[2].lines[line_name]
308 elif pt1.x == 0. and pt2.x == 0.:
309 qq.lines[line_name] = quadrants[4].lines[line_name]
310 else:
311 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])
313 elif line.type == 'arc':
314 center = Func.arc_center_from_3_points([pt1.x, pt1.y],
315 [iron.key_points[line.kp3].x, iron.key_points[line.kp3].y],
316 [pt2.x, pt2.y])
317 new_point_name = 'kp' + line_name + '_center'
318 arc_coordinates1 = (pt1.x, pt1.y)
319 arc_coordinates2 = (pt2.x, pt2.y)
320 arc_coordinates3 = (iron.key_points[line.kp3].x, iron.key_points[line.kp3].y)
322 # This code addresses a meshing error in MQXA and MB_2COILS that occurs when an arc is defined on any of
323 # the axes. The issue arises because the function Func.arc_center_from_3_points does not return exactly
324 # zero but a value with a magnitude of approximately 10^-17 when the two points are placed on the axes.
325 # Consequently, when using the method self.occ.addCircleArc(), which only takes in three points without
326 # specifying a direction, a problem arises. The addCircleArc() function always creates the arc with the
327 # smallest angle. However, since center point can be slightly above or below the axis, the arc can
328 # inadvertently be drawn in the wrong quadrant, leading to an incorrect result.
329 # -----------------------
330 # Check that arcs with points on the x-axis are drawn in the first quadrant
331 if arc_coordinates3[1] > 0 and arc_coordinates2[1] == 0 and arc_coordinates1[1] == 0 and center[1] > 0:
332 quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], -center[1], 0)
333 # Check that arcs with points on the y-axis are drawn in the first quadrant
334 elif arc_coordinates3[0] > 0 and arc_coordinates2[0] == 0 and arc_coordinates1[0] == 0 and center[0] > 0:
335 quadrants[1].points[new_point_name] = self.occ.addPoint(-center[0], center[1], 0)
336 else:
337 quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], center[1], 0)
338 # -----------------------
339 # gmsh.model.setEntityName(0, gm.iron.quadrants[1].points[new_point_name], 'iron_' + new_point_name)
340 if symmetry in ['x', 'none']:
341 if center[0] == 0.:
342 quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]
343 else:
344 quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
345 self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)
346 if symmetry in ['y', 'none']:
347 if center[1] == 0.:
348 quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]
349 else:
350 quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
351 self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)
352 if symmetry == 'none':
353 if center[1] == 0.:
354 quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]
355 else:
356 quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]
357 self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)
359 for quadrant, qq in quadrants.items():
360 qq.lines[line_name] = self.occ.addCircleArc(
361 qq.points[line.kp1], qq.points[new_point_name], qq.points[line.kp2])
363 elif line.type == 'circle':
364 center = [(pt1.x + pt2.x) / 2, (pt1.y + pt2.y) / 2]
365 radius = (np.sqrt(np.square(pt1.x - center[0]) + np.square(pt1.y - center[1])) +
366 np.sqrt(np.square(pt2.x - center[0]) + np.square(pt2.y - center[1]))) / 2
368 for quadrant, qq in quadrants.items():
369 qq.lines[line_name] = self.occ.addCircle(
370 mirror_x[quadrant - 1] * center[0], mirror_y[quadrant - 1] * center[1], 0, radius)
371 qq.points['kp' + line_name] = len(qq.points) + 1
373 elif line.type == 'ellipticArc':
374 a, b = line.arg1, line.arg2
375 x1, y1 = pt1.x, pt1.y
376 x2, y2 = pt2.x, pt2.y
377 x3 = np.power(x1, 2.0)
378 y3 = np.power(y1, 2.0)
379 x4 = np.power(x2, 2.0)
380 y4 = np.power(y2, 2.0)
381 a2 = np.power(a, 2.0)
382 b2 = np.power(b, 2.0)
383 expression = -4.0 * a2 * b2 + a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 - 2.0 * b2 * x1 * x2 + b2 * x4
384 xc = x1 / 2.0 + x2 / 2.0 - a * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 -
385 2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (y1 - y2) / (2.0 * b)
386 yc = y1 / 2.0 + y2 / 2.0 + b * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3
387 - 2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (x1 - x2) / (2.0 * a)
389 center = self.occ.addPoint(xc, yc, 0, lc)
390 axis_point_a = self.occ.addPoint(xc + a, yc, 0, lc)
391 axis_point_b = self.occ.addPoint(xc, yc + b, 0, lc)
393 new_point_name = 'kp' + line_name + '_center'
394 new_axis_a_point_name = 'kp' + line_name + '_a'
395 new_axis_b_point_name = 'kp' + line_name + '_b'
397 quadrants[1].points[new_point_name] = center
398 quadrants[1].points[new_axis_a_point_name] = axis_point_a
399 quadrants[1].points[new_axis_b_point_name] = axis_point_b
401 if symmetry in ['x', 'none']:
402 if xc == 0.: # Least amount of possible points.
403 quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]
404 quadrants[2].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]
405 quadrants[2].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]
406 else:
407 quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
408 quadrants[2].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]
409 quadrants[2].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]
410 self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)
411 self.occ.mirror([(0, quadrants[2].points[new_axis_a_point_name])], 1, 0, 0, 0)
412 self.occ.mirror([(0, quadrants[2].points[new_axis_b_point_name])], 1, 0, 0, 0)
413 if symmetry in ['y', 'none']:
414 if yc == 0.:
415 quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]
416 quadrants[4].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]
417 quadrants[4].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]
418 else:
419 quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]
420 self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)
421 quadrants[4].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]
422 self.occ.mirror([(0, quadrants[4].points[new_axis_a_point_name])], 0, 1, 0, 0)
423 quadrants[4].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]
424 self.occ.mirror([(0, quadrants[4].points[new_axis_b_point_name])], 0, 1, 0, 0)
425 if symmetry == 'none':
426 if yc == 0.:
427 quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]
428 quadrants[3].points[new_axis_a_point_name] = quadrants[2].points[new_axis_a_point_name]
429 quadrants[3].points[new_axis_b_point_name] = quadrants[2].points[new_axis_b_point_name]
430 else:
431 quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]
432 self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)
433 quadrants[3].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_a_point_name])])[0][1]
434 self.occ.mirror([(0, quadrants[3].points[new_axis_a_point_name])], 0, 1, 0, 0)
435 quadrants[3].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_b_point_name])])[0][1]
436 self.occ.mirror([(0, quadrants[3].points[new_axis_b_point_name])], 0, 1, 0, 0)
438 for quadrant, qq in quadrants.items():
439 qq.lines[line_name] = self.occ.addEllipseArc(
440 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],
441 qq.points[line.kp2])
443 else:
444 raise ValueError('Hyper line {} not supported'.format(line.type))
446 if symmetry != 'none':
447 indexes = {'x_p': 1, 'y_p': 1, 'x_n': 1, 'y_n': 1}
448 self.md.geometries.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)
449 for sym in list_bnds:
450 if sym in ['x_p', 'y_p']:
451 quadrant = 1
452 elif sym == 'x_n':
453 quadrant = 2
454 else: # 'y_n'
455 quadrant = 4
456 sym_lines_tags[sym] = [x for _, x in sorted(zip(symmetric_bnds_order[sym[0]], sym_lines_tags[sym]))]
458 self.symmetric_bnds[sym]['pnts'].append(['center', 0])
459 self.symmetric_bnds[sym]['pnts'].sort(key=lambda x: x[1])
460 self.md.geometries.symmetric_boundaries.lines[sym + '_center'] = self.occ.addLine(
461 self.md.geometries.air_inf.points['center'], quadrants[quadrant].points[self.symmetric_bnds[sym]['pnts'][1][0]])
462 sym_lines_tags[sym].insert(0, self.md.geometries.symmetric_boundaries.lines[sym + '_center'])
463 for i, pnt in enumerate(self.symmetric_bnds[sym]['pnts'][1:-1]):
464 pnt_next = self.symmetric_bnds[sym]['pnts'][i + 2][0]
465 if not any(pnt[0] in s and pnt_next in s for s in self.symmetric_bnds[sym]['line_pnts']):
466 self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]] =\
467 self.occ.addLine(quadrants[quadrant].points[pnt[0]], quadrants[quadrant].points[pnt_next])
468 sym_lines_tags[sym].insert(indexes[sym], self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]])
469 indexes[sym] += 1
470 if symmetry == 'xy':
471 self.symmetric_loop_lines['x'] = sym_lines_tags['x_p']
472 sym_lines_tags['y_p'].reverse()
473 self.symmetric_loop_lines['y'] = sym_lines_tags['y_p']
474 elif symmetry == 'x':
475 sym_lines_tags['x_n'].reverse()
476 self.symmetric_loop_lines['x'] = sym_lines_tags['x_n'] + sym_lines_tags['x_p']
477 elif symmetry == 'y':
478 sym_lines_tags['y_p'].reverse()
479 self.symmetric_loop_lines['y'] = sym_lines_tags['y_p'] + sym_lines_tags['y_n']
481 for quadrant, qq in quadrants.items():
482 for area_name, area in iron.hyper_areas.items():
483 qq.areas[area_name] = dM.Area(loop=self.occ.addCurveLoop([qq.lines[line] for line in area.lines]))
484 if (iron.hyper_areas[area_name].material not in self.md.domains.groups_entities.iron and
485 iron.hyper_areas[area_name].material != 'BH_air'):
486 self.md.domains.groups_entities.iron[iron.hyper_areas[area_name].material] = []
488 def constructWedgeGeometry(self, use_TSA):
489 """
490 Generates points, hyper lines, and curve loops for the wedges
491 """
492 def _addMidLayerThinShellPoints(wedge_current):
493 def __addThinShellPoints(side_case, mid_layer_ts):
494 if side_case == 'outer':
495 mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +
496 Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2
497 else:
498 mean_rad_current = (Func.points_distance([wedge_current.iH.x, wedge_current.iH.y], wedge_center) +
499 Func.points_distance([wedge_current.iL.x, wedge_current.iL.y], wedge_center)) / 2
500 are_endpoints = {}
501 for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():
502 blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]
503 blk_next = wnd.blocks[blk_nr_next]
504 ht_list_next = (list(blk_next.half_turns.keys()) if blk_nr_next == list(wnd.blocks.keys())[0] else list(
505 reversed(blk_next.half_turns.keys())))
506 hh = blk_next.half_turns[ht_list_next[-1]].corners.bare
507 ll = blk_next.half_turns[ht_list_next[0]].corners.bare
508 bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
509 if side_case == 'outer':
510 block_list = self.md.geometries.coil.anticlockwise_order.coils[wedge.order_l.coil].layers[wedge.order_l.layer + 1]
511 blk_index = [blk.block for blk in block_list].index(blk_nr_next)
512 if blk_index + 1 == len(block_list): blk_index = -1
513 for blk in block_list[blk_index + 1:] + block_list[:blk_index + 1]:
514 if blk.winding == block_list[blk_index].winding:
515 ht_index = -1
516 break
517 elif blk.pole != block_list[blk_index].pole:
518 ht_index = 0
519 break
520 hh = blk_next.half_turns[ht_list_next[ht_index]].corners.bare
521 ll = blk_next.half_turns[ht_list_next[0 if ht_index == -1 else -1]].corners.bare
522 mean_rad_next = (Func.points_distance([hh.iH.x, hh.iH.y], wedge_center) +
523 Func.points_distance([ll.iL.x, ll.iL.y], wedge_center)) / 2
524 else:
525 mean_rad_next = (Func.points_distance([bc_next.oH.x, bc_next.oH.y], wedge_center) +
526 Func.points_distance([bc_next.oL.x, bc_next.oL.y], wedge_center)) / 2
527 mean_rad = (mean_rad_current + mean_rad_next) / 2
528 mid_layer = self.findMidLayerPoint(wedge_current.oH, bc_next.iH, wedge.corrected_center.outer, mean_rad)\
529 if side_case == 'outer' else self.findMidLayerPoint(wedge_current.iH, bc_next.oH, wedge.corrected_center.inner, mean_rad)
530 are_endpoints[wnd_nr] = self.getMidLayerEndpoints(wedge_current, bc_next, wedge_center, mid_layer_arc_pnt=mid_layer)
531 for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():
532 blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]
533 blk_next = wnd.blocks[blk_nr_next]
534 is_first_blk_next = blk_nr_next == list(wnd.blocks.keys())[0]
535 ht_list_next = (list(blk_next.half_turns.keys()) if is_first_blk_next else list(
536 reversed(blk_next.half_turns.keys())))
537 if are_endpoints[wnd_nr]: # this is empty if the wedge and the block are not radially adjacent
538 endpoints = are_endpoints[wnd_nr][0]
539 which_entity = are_endpoints[wnd_nr][1]
540 mid_layer_name = 'w' + str(wedge_nr) + '_' + str(blk_nr_next)
541 mid_layer_ts[mid_layer_name] = dM.Region()
542 ts_wdg = mid_layer_ts[mid_layer_name]
543 beg = ('w' + str(wedge_nr) if which_entity['lower'] == 'current' else str(ht_list_next[0])) + 'l'
544 ts_wdg.points[beg] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
545 ht_lower_angles = {}
546 for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):
547 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]],
548 [[ht.corners.bare.oL.x, ht.corners.bare.oL.y], [ht.corners.bare.oH.x, ht.corners.bare.oH.y]],
549 ['l', 'h']):
550 line_pars_current = Func.line_through_two_points(pnt1, pnt2)
551 intersect_prev = Func.intersection_between_arc_and_line(
552 line_pars_current, [wedge_center, endpoints['higher'], endpoints['lower']])
553 if intersect_prev:
554 ts_wdg.points[str(ht_nr) + side] = self.occ.addPoint(intersect_prev[0][0], intersect_prev[0][1], 0)
555 elif side == 'l':
556 intrsc = Func.intersection_between_circle_and_line(line_pars_current, [wedge_center, endpoints['lower']], get_only_closest=True)[0]
557 ht_lower_angles[ht_nr] = Func.arc_angle_between_point_and_abscissa([intrsc[0], intrsc[1]], wedge_center)
558 end = ('w' + str(wedge_nr) if which_entity['higher'] == 'current' else str(ht_list_next[-1])) + 'h'
559 if all('w' in pnt_name for pnt_name in list(ts_wdg.points.keys())): # only one thin-shell 'within' the facing half-turn
560 wdg_angle_il = Func.arc_angle_between_point_and_abscissa([endpoints['lower'][0], endpoints['lower'][1]], wedge_center)
561 for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):
562 if ht_lower_angles[ht_nr] > wdg_angle_il: break
563 prev_nr = str(ht_nr)
564 end = prev_nr + 'h'
565 ts_wdg.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
567 # Create auxiliary thin shells for outliers
568 # if both corners belong to thin shells, continue
569 used_wdg_corners = [False, False]
570 for ep in are_endpoints.values():
571 if ep is not None:
572 if ep[1]['higher'] == 'current': used_wdg_corners[1] = True
573 if ep[1]['lower'] == 'current': used_wdg_corners[0] = True
574 if side_case == 'inner':
575 for ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.keys():
576 if ts_name[ts_name.index('_') + 1:] == 'w' + str(wedge_nr):
577 for ep_key, ep in are_endpoints_wdg[int(ts_name[1:ts_name.index('_')])].items():
578 if ep is not None:
579 if ep[1]['higher'] == 'next': used_wdg_corners[1] = True
580 if ep[1]['lower'] == 'next': used_wdg_corners[0] = True
581 else:
582 if wedge_nr in are_endpoints_wdg:
583 for ep in are_endpoints_wdg[wedge_nr].values():
584 if ep is not None:
585 if ep[1]['higher'] == 'current': used_wdg_corners[1] = True
586 if ep[1]['lower'] == 'current': used_wdg_corners[0] = True
587 if not used_wdg_corners[1]:
588 for wdg_nr, wdg in self.geom.wedges.items():
589 if blk_nr_next == wdg.order_l.block: used_wdg_corners[1] = True
590 if not used_wdg_corners[0]:
591 for wdg_nr, wdg in self.geom.wedges.items():
592 if blk_nr_next == wdg.order_h.block: used_wdg_corners[0] = True
593 if not all(used_wdg_corners):
594 def ___create_aux_mid_layer_point(ss, points):
595 mid_layer_ts_aux[mid_layer_name] = dM.Region()
596 circle_pnt = [endpoints[ss][0], endpoints[ss][1]]
597 inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),
598 [[wedge.corrected_center.outer.x, wedge.corrected_center.outer.y], circle_pnt], get_only_closest=True)[0]
599 mid_layer_ts_aux[mid_layer_name].points[str(wedge_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)
600 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)
601 mid_layer_ts_aux[mid_layer_name].lines['w' + str(wedge_nr)] = 0
602 if which_entity['higher'] == 'current' and which_entity['lower'] != 'current':
603 ___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],
604 [wedge_current.oL.x, wedge_current.oL.y]])
605 elif which_entity['higher'] != 'current' and which_entity['lower'] == 'current':
606 ___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],
607 [wedge_current.oH.x, wedge_current.oH.y]])
608 else: # whole block 'within' the facing wedge
609 for wdg_nr, wdg in self.geom.wedges.items():
610 if blk_nr_next == wdg.order_h.block:
611 ___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],
612 [wedge_current.oH.x, wedge_current.oH.y]])
613 break
614 elif blk_nr_next == wdg.order_l.block:
615 ___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],
616 [wedge_current.oL.x, wedge_current.oL.y]])
617 break
619 pole = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole]
620 blk_list_current = list(pole.layers[wedge.order_l.layer].windings[wedge.order_l.winding].blocks.keys())
621 if wedge.order_l.layer < len(pole.layers):
622 __addThinShellPoints('outer', self.md.geometries.thin_shells.mid_layers_wdg_to_ht)
623 if wedge.order_l.layer > 1:
624 __addThinShellPoints('inner', self.md.geometries.thin_shells.mid_layers_ht_to_wdg)
626 wedges = self.md.geometries.wedges
627 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
628 wedge_data = {}
630 wdgs_corners = {}
631 for wedge_nr, wedge in self.geom.wedges.items():
632 wdgs_corners[wedge_nr] = {}
633 corners = wdgs_corners[wedge_nr]
634 if wedge.order_l.coil not in wedges.coils:
635 wedges.coils[wedge.order_l.coil] = dM.WedgeLayer()
636 if wedge.order_l.layer not in wedges.coils[wedge.order_l.coil].layers:
637 wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer] = dM.WedgeRegion()
638 wedge_layer = wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer]
639 wedge_layer.wedges[wedge_nr] = dM.Region()
640 wedge_reg = wedge_layer.wedges[wedge_nr]
641 wedge_layer.block_prev[wedge_nr] = wedge.order_l.block
642 wedge_layer.block_next[wedge_nr] = wedge.order_h.block
643 wnd = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[
644 wedge.order_l.layer].windings[wedge.order_l.winding]
645 wnd_next = self.geom.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[
646 wedge.order_h.layer].windings[wedge.order_h.winding]
647 block = wnd.blocks[wedge.order_l.block]
648 block_next = wnd_next.blocks[wedge.order_h.block]
649 corners['last_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[
650 wedge.order_l.layer].windings[wedge.order_l.winding].blocks[wedge.order_l.block].half_turns.areas.keys())[-1])
651 corners['first_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[
652 wedge.order_h.layer].windings[wedge.order_h.winding].blocks[wedge.order_h.block].half_turns.areas.keys())[0])
653 ht_current = block.half_turns[corners['last_ht']].corners.bare
654 ht_next = block_next.half_turns[corners['first_ht']].corners.bare
655 d_current = self.data.conductors[wnd.conductor_name].cable.th_insulation_along_width * 2
656 d_next = self.data.conductors[wnd_next.conductor_name].cable.th_insulation_along_width * 2
657 for pnt_close, pnt_far, wdg_corner, d in zip([ht_current.iH, ht_current.oH, ht_next.iL, ht_next.oL],
658 [ht_current.iL, ht_current.oL, ht_next.iH, ht_next.oH],
659 ['il', 'ol', 'ih', 'oh'], [d_current, d_current, d_next, d_next]):
660 if abs(pnt_far.x - pnt_close.x) > 0.:
661 m = (pnt_far.y - pnt_close.y) / (pnt_far.x - pnt_close.x)
662 b = pnt_close.y - m * pnt_close.x
663 root = np.sqrt(- pnt_close.x ** 2 * m ** 2 - 2 * pnt_close.x * b * m + 2 * pnt_close.x * pnt_close.y * m
664 - b ** 2 + 2 * b * pnt_close.y - pnt_close.y ** 2 + d ** 2 * m ** 2 + d ** 2)
665 pnt1_x = (pnt_close.x - b * m + pnt_close.y * m + root) / (m ** 2 + 1)
666 pnt1_y = m * pnt1_x + b
667 pnt2_x = (pnt_close.x - b * m + pnt_close.y * m - root) / (m ** 2 + 1)
668 pnt2_y = m * pnt2_x + b
669 corners[wdg_corner] = Coord(x=pnt1_x, y=pnt1_y) if Func.points_distance([pnt1_x, pnt1_y], [pnt_far.x, pnt_far.y]) >\
670 Func.points_distance([pnt_close.x, pnt_close.y], [pnt_far.x, pnt_far.y]) else Coord(x=pnt2_x, y=pnt2_y)
671 else:
672 bore_cnt_x = self.geom.coil.coils[wedge.order_l.coil].bore_center.x
673 pnt1_y, pnt2_y = pnt_close.y + d, pnt_close.y - d
674 corners[wdg_corner] = Coord(x=pnt_close.x,
675 y=pnt1_y if (wdg_corner[-1] == 'l' and pnt_close.x > bore_cnt_x) or
676 (wdg_corner[-1] == 'h' and pnt_close.x < bore_cnt_x) else pnt2_y)
677 wedge_reg.points[wdg_corner] = self.occ.addPoint(corners[wdg_corner].x, corners[wdg_corner].y, 0)
678 inner = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
679 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],
680 [corners['ih'].x, corners['ih'].y], [corners['il'].x, corners['il'].y])
681 outer = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
682 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],
683 [corners['oh'].x, corners['oh'].y], [corners['ol'].x, corners['ol'].y])
684 wedge_data[wedge_nr] = [Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']), wedge.corrected_center.outer]
685 wedge_reg.points['inner_center'] = self.occ.addPoint(inner[0], inner[1], 0)
686 wedge_reg.points['outer_center'] = self.occ.addPoint(outer[0], outer[1], 0)
687 wedge_reg.lines['h'] = self.occ.addLine(wedge_reg.points['ih'], wedge_reg.points['oh'])
688 wedge_reg.lines['l'] = self.occ.addLine(wedge_reg.points['il'], wedge_reg.points['ol'])
689 wedge_reg.lines['i'] = self.occ.addCircleArc(wedge_reg.points['ih'], wedge_reg.points['inner_center'], wedge_reg.points['il'])
690 wedge_reg.lines['o'] = self.occ.addCircleArc(wedge_reg.points['oh'], wedge_reg.points['outer_center'], wedge_reg.points['ol'])
691 wedge_reg.areas[str(wedge_nr)] = dM.Area(loop=self.occ.addCurveLoop(
692 [wedge_reg.lines['i'], wedge_reg.lines['l'], wedge_reg.lines['o'], wedge_reg.lines['h']]))
694 if use_TSA:
695 # Wedge thin shells
696 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_wdg_to_wdg
697 are_endpoints_wdg = {}
698 for coil_nr, coil in self.md.geometries.wedges.coils.items():
699 layer_list = list(coil.layers.keys())
700 for layer_nr, layer in coil.layers.items():
701 if layer_list.index(layer_nr) + 1 < len(layer_list):
702 for wedge_nr, wedge in layer.wedges.items():
703 are_endpoints_wdg[wedge_nr] = {}
704 are_endpoints = are_endpoints_wdg[wedge_nr]
705 wedge_current = wedge_data[wedge_nr][0]
706 wedge_center = [wedge_data[wedge_nr][1].x, wedge_data[wedge_nr][1].y]
707 mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +
708 Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2
709 for wdg_next_nr, wdg_next in coil.layers[layer_nr + 1].wedges.items():
710 if self.geom.wedges[wedge_nr].order_l.pole == self.geom.wedges[wdg_next_nr].order_l.pole:
711 wedge_next = wedge_data[wdg_next_nr][0]
712 mean_rad_next = (Func.points_distance([wedge_next.iH.x, wedge_next.iH.y], wedge_center) +
713 Func.points_distance([wedge_next.iL.x, wedge_next.iL.y], wedge_center)) / 2
714 mean_rad = (mean_rad_current + mean_rad_next) / 2
715 mid_layer = self.findMidLayerPoint(wedge_current.oH, wedge_next.iH, wedge_data[wedge_nr][1], mean_rad)
716 are_endpoints[wdg_next_nr] = self.getMidLayerEndpoints(wedge_current, wedge_next, wedge_center, mid_layer_arc_pnt=mid_layer)
717 if are_endpoints[wdg_next_nr]: # this is empty if the wedges are not radially adjacent
718 endpoints = are_endpoints[wdg_next_nr][0]
719 mid_layer_name = 'w' + str(wedge_nr) + '_w' + str(wdg_next_nr)
720 mid_layer_ts[mid_layer_name] = dM.Region()
721 ts = mid_layer_ts[mid_layer_name]
722 ts.points['center'] = self.occ.addPoint(wedge_center[0], wedge_center[1], 0)
723 ts.points['beg'] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
724 end = 'w' + str(wedge_nr if are_endpoints[wdg_next_nr][1] == 'current' else wdg_next_nr)
725 ts.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
727 # Half-turn thin shells
728 for wedge_nr, wedge in self.geom.wedges.items():
729 corners = wdgs_corners[wedge_nr]
730 # Mid layer lines
731 wedge_center = [self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,
732 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y]
733 _addMidLayerThinShellPoints(Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']))
734 # Mid wedge-turn lines
735 mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn
736 for adj_blk, ht, inner, outer in zip([wedge.order_l, wedge.order_h], [corners['last_ht'], corners['first_ht']],
737 [corners['il'], corners['ih']], [corners['ol'], corners['oh']]):
738 mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)] = dM.Region()
739 ts = mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)]
740 ht_corners = self.geom.coil.coils[adj_blk.coil].poles[adj_blk.pole].layers[
741 adj_blk.layer].windings[adj_blk.winding].blocks[adj_blk.block].half_turns[ht].corners.bare
742 ht_corners_i = ht_corners.iH if ht == corners['last_ht'] else ht_corners.iL
743 ht_corners_o = ht_corners.oH if ht == corners['last_ht'] else ht_corners.oL
744 mid_inner = [(inner.x + ht_corners_i.x) / 2, (inner.y + ht_corners_i.y) / 2]
745 mid_outer = [(outer.x + ht_corners_o.x) / 2, (outer.y + ht_corners_o.y) / 2]
746 line_name = 'w' + str(wedge_nr) + '_' + str(ht)
747 ts.points[line_name + '_i'] = self.occ.addPoint(mid_inner[0], mid_inner[1], 0)
748 ts.points[line_name + '_o'] = self.occ.addPoint(mid_outer[0], mid_outer[1], 0)
750 def constructCoilGeometry(self, run_type):
751 """
752 Generates points, hyper lines, and curve loops for the coil half-turns
753 """
754 symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'
755 # Sub domains angles: first key means 'from 0 to x'; second key means 'from x to 2*pi'
756 if symmetry == 'xy':
757 angle_range = {'to': np.pi / 2, 'from': 2 * np.pi}
758 elif symmetry == 'x':
759 angle_range = {'to': np.pi, 'from': 2 * np.pi}
760 elif symmetry == 'y':
761 angle_range = {'to': np.pi / 2, 'from': 3 / 2 * np.pi}
762 elif symmetry == 'none':
763 angle_range = {'to': 2 * np.pi, 'from': 0}
764 else:
765 raise Exception('Symmetry plane not supported.')
767 def _addMidLayerThinShellPoints(pnt_params, ss, name, case):
768 endpnts, cnt = ts_endpoints[name]
769 if len(pnt_params) == 3: # line parameters (cos-theta Rutherford)
770 intersect[name] = Func.intersection_between_arc_and_line(pnt_params, [cnt, endpnts['higher'], endpnts['lower']])
771 if intersect[name]:
772 intersect[name] = intersect[name][0]
773 pnt_angle = Func.arc_angle_between_point_and_abscissa(intersect[name], cnt)
774 elif len(pnt_params) == 4: # points coordinates (cos-theta Mono)
775 wnd_next = list(pole.layers[layer_nr + (1 if case == 'current' else -1)].windings.keys())[
776 list(pole.layers[layer_nr].windings.keys()).index(winding_nr)]
777 blk_next = pole.layers[layer_nr + (1 if case == 'current' else -1)].windings[wnd_next].blocks[
778 int(ts_name[ts_name.index('_') + 1:] if case == 'current' else ts_name[:ts_name.index('_')])]
779 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
780 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)
781 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]
782 pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)
783 pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)
784 pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)
785 intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None
786 else: # point coordinates (block-coil)
787 pnt = [endpnts['higher'][0], pnt_params[1]] if coil.type == 'common-block-coil' else [pnt_params[0], endpnts['higher'][1]]
788 if abs(endpnts['higher'][1]) > 1e-6:
789 pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)
790 pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)
791 pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)
792 else:
793 pnt_angle = abs(pnt_params[0])
794 pnt_angle_h = abs(endpnts['higher'][0])
795 pnt_angle_l = abs(endpnts['lower'][0])
796 intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None
797 if intersect[name]:
798 mid_layer_ts[name].mid_layers.points[str(halfTurn_nr) + ss] = \
799 self.occ.addPoint(intersect[name][0], intersect[name][1], 0)
800 mid_layer_ts[name].point_angles[str(halfTurn_nr) + ss] = Func.sig_dig(pnt_angle)
801 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):
802 intersect[name] = pnt
803 return intersect
805 def _addMidLayerThinShellGroup(cl, for_mid_pole=False, mid_coil=False):
806 is_first_blk_next = block_nr_next == list(winding_next.blocks.keys())[0]
807 if 'solenoid' in cl.type:
808 ht_list_next = list(reversed(block_next.half_turns.keys()) if layer_nr % 2 == 0 else list(block_next.half_turns.keys()))
809 elif cl.type == 'reversed-block-coil':
810 ht_list_next = (list(block_next.half_turns.keys()) if not is_first_blk_next else list(reversed(block_next.half_turns.keys())))
811 else:
812 ht_list_next = (list(block_next.half_turns.keys()) if is_first_blk_next else list(reversed(block_next.half_turns.keys())))
813 hh = block_next.half_turns[ht_list_next[-1]].corners.bare
814 ll = block_next.half_turns[ht_list_next[0]].corners.bare
815 bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
816 if 'block-coil' in cl.type or (cable_type_curr in ['Mono', 'Ribbon'] and not mid_coil):
817 center = [cl.bore_center.x, cl.bore_center.y]
818 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)
819 else:
820 mean_rad_next = (Func.points_distance([bc_next.iH.x, bc_next.iH.y], [cl.bore_center.x, cl.bore_center.y]) +
821 Func.points_distance([bc_next.iL.x, bc_next.iL.y], [cl.bore_center.x, cl.bore_center.y])) / 2
822 mean_rad = (mean_rad_current + mean_rad_next) / 2
823 mid_layer_h = self.findMidLayerPoint(bc_current.oH, bc_next.iH, cl.bore_center, mean_rad)
824 mid_layer_l = self.findMidLayerPoint(bc_current.oL, bc_next.iL, cl.bore_center, mean_rad)
825 mid_ht_next_i = int(len(ht_list_next) / 2) if len(ht_list_next) % 2 == 0 else round(len(ht_list_next) / 2)
826 mid_ht_next = block_next.half_turns[ht_list_next[mid_ht_next_i - 1]].corners.insulated
827 mid_layer_m = self.findMidLayerPoint(mid_ht_current.oH, mid_ht_next.iH, cl.bore_center, mean_rad)
828 center = Func.arc_center_from_3_points(mid_layer_h, mid_layer_m, mid_layer_l)
829 are_endpoints = self.getMidLayerEndpoints(bc_current, bc_next, center, mid_layer_arc_pnt=mid_layer_h, cable_type=cable_type_curr)
830 if are_endpoints: # this is empty if the blocks are not radially adjacent
831 endpoints = are_endpoints[0]
832 which_block = are_endpoints[1]
833 mid_layer_name = blk_nr + '_' + str(block_nr_next)
834 if for_mid_pole:
835 block_coil_mid_pole_next_blks_list[block_nr_next].append(mid_layer_name)
836 block_coil_ts_endpoints[mid_layer_name] = [endpoints, center]
837 else:
838 if block_nr_next not in list(next_blks_list.keys()):
839 next_blks_list[block_nr_next] = []
840 next_blks_list[block_nr_next].append(mid_layer_name)
841 ts_endpoints[mid_layer_name] = [endpoints, center]
842 mid_layer_ts[mid_layer_name] = dM.MidLayer()
843 mid_layer_ts[mid_layer_name].half_turn_lists[blk_nr] = ht_list
844 mid_layer_ts[mid_layer_name].half_turn_lists[str(block_nr_next)] = ht_list_next
845 beg = (str(ht_list[0]) if which_block['lower'] == 'current' else str(ht_list_next[0])) + 'l'
846 mid_layer_ts[mid_layer_name].mid_layers.points[beg] = \
847 self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)
848 end = (str(ht_list[-1]) if which_block['higher'] == 'current' else str(ht_list_next[-1])) + 'h'
849 mid_layer_ts[mid_layer_name].mid_layers.points[end] = \
850 self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)
851 if not for_mid_pole or (for_mid_pole and abs(endpoints['higher'][1]) > 1e-6):
852 mid_layer_ts[mid_layer_name].point_angles[beg] =\
853 Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['lower'], center))
854 mid_layer_ts[mid_layer_name].point_angles[end] =\
855 Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['higher'], center))
856 else:
857 mid_layer_ts[mid_layer_name].point_angles[beg] = abs(endpoints['lower'][0])
858 mid_layer_ts[mid_layer_name].point_angles[end] = abs(endpoints['higher'][0])
860 # Create anticlockwise order of blocks
861 present_blocks = []
862 block_corner_angles = {}
863 concentric_coils = self.md.geometries.coil.concentric_coils
864 acw_order = self.md.geometries.coil.anticlockwise_order.coils
865 self.md.geometries.coil.physical_order = self.geom.coil.physical_order
866 for coil_nr, coil in self.geom.coil.coils.items():
867 # if coil_nr not in block_corner_angles:
868 block_corner_angles[coil_nr] = {}
869 if (coil.bore_center.x, coil.bore_center.y) not in concentric_coils:
870 concentric_coils[(coil.bore_center.x, coil.bore_center.y)] = []
871 concentric_coils[(coil.bore_center.x, coil.bore_center.y)].append(coil_nr)
872 for pole_nr, pole in coil.poles.items():
873 for layer_nr, layer in pole.layers.items():
874 if layer_nr not in block_corner_angles[coil_nr]:
875 block_corner_angles[coil_nr][layer_nr] = {}
876 blk_angles = block_corner_angles[coil_nr][layer_nr]
877 for winding_nr, winding in layer.windings.items():
878 for block_nr, block in winding.blocks.items():
879 blk_angles[block_nr] = {'angle': Func.sig_dig(Func.arc_angle_between_point_and_abscissa(
880 [block.block_corners.iL.x, block.block_corners.iL.y],
881 [coil.bore_center.x, coil.bore_center.y])), 'keys': [pole_nr, winding_nr]}
882 higher_angle = Func.sig_dig(Func.arc_angle_between_point_and_abscissa(
883 [block.block_corners.iH.x, block.block_corners.iH.y],
884 [coil.bore_center.x, coil.bore_center.y]))
885 if ((blk_angles[block_nr]['angle'] <= angle_range['to'] and higher_angle <= angle_range['to']) or
886 (angle_range['from'] <= blk_angles[block_nr]['angle'] and angle_range['from'] <= higher_angle)):
887 present_blocks.append(block_nr)
888 for coil_nr, coil in block_corner_angles.items():
889 acw_order[coil_nr] = dM.LayerOrder()
890 for layer_nr, layer in coil.items():
891 acw_order[coil_nr].layers[layer_nr] = []
892 ordered_blocks = [[block_nr, block['angle'], block['keys']] for block_nr, block in layer.items()]
893 ordered_blocks.sort(key=lambda x: x[1])
894 for blk in ordered_blocks:
895 if blk[0] in present_blocks:
896 acw_order[coil_nr].layers[layer_nr].append(dM.AnticlockwiseOrder(pole=blk[2][0], winding=blk[2][1], block=blk[0]))
898 # Check if there are concentric coils
899 for bore_center, coils in concentric_coils.items():
900 if len(coils) > 1:
901 radii = []
902 for coil_nr in coils:
903 lyr = self.geom.coil.coils[coil_nr].poles[1].layers[1]
904 blk = list(lyr.windings.keys())[0]
905 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])])
906 radii.sort(key=lambda x: x[1])
907 concentric_coils[bore_center] = [rad[0] for rad in radii]
909 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
910 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht
911 # Collect block couples for block-coil mid-pole thin shells
912 block_coil_mid_pole_next_blks_list = {}
913 block_coil_ts_endpoints = {}
914 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
915 if self.geom.coil.coils[coil_nr].type in ['block-coil', 'reversed-block-coil']:
916 self.block_coil_mid_pole_blks[coil_nr] = []
917 first_lyr = list(coil.layers.keys())[0]
918 layer = coil.layers[first_lyr]
919 for nr, block_order in enumerate(layer):
920 blk_next_index = nr + 1 if nr + 1 < len(layer) else 0
921 if layer[blk_next_index].pole != block_order.pole:
922 self.block_coil_mid_pole_blks[coil_nr].append([block_order, layer[blk_next_index]])
923 block_coil_mid_pole_next_blks_list[layer[blk_next_index].block] = []
924 # Mid pole lines for block-coils
925 for coil_nr, coil in self.block_coil_mid_pole_blks.items():
926 coil_geom = self.geom.coil.coils[coil_nr]
927 for mid_pole in coil:
928 winding = self.geom.coil.coils[coil_nr].poles[mid_pole[0].pole].layers[1].windings[mid_pole[0].winding]
929 cable_type_curr = self.data.conductors[winding.conductor_name].cable.type
930 block_nr = mid_pole[0].block
931 blk_nr = str(block_nr)
932 block = winding.blocks[block_nr]
933 is_first_blk = block_nr == list(winding.blocks.keys())[0]
934 if coil_geom.type == 'reversed-block-coil':
935 ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))
936 else:
937 ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))
938 hh = block.half_turns[ht_list[-1]].corners.bare
939 ll = block.half_turns[ht_list[0]].corners.bare
940 bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
941 winding_next = self.geom.coil.coils[coil_nr].poles[mid_pole[1].pole].layers[1].windings[mid_pole[1].winding]
942 block_nr_next = mid_pole[1].block
943 block_next = winding_next.blocks[block_nr_next]
944 _addMidLayerThinShellGroup(coil_geom, for_mid_pole=True)
946 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
947 self.md.geometries.coil.physical_order = self.geom.coil.physical_order
948 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
949 next_blks_list = block_coil_mid_pole_next_blks_list.copy()
950 ts_endpoints = block_coil_ts_endpoints.copy()
951 for coil_nr, coil in self.geom.coil.coils.items():
952 self.md.geometries.coil.coils[coil_nr] = dM.Pole()
953 coils = self.md.geometries.coil.coils[coil_nr]
954 coils.type = coil.type
955 coils.bore_center = coil.bore_center
956 for pole_nr, pole in coil.poles.items():
957 coils.poles[pole_nr] = dM.Layer()
958 poles = coils.poles[pole_nr]
959 for layer_nr, layer in pole.layers.items():
960 poles.layers[layer_nr] = dM.Winding()
961 layers = poles.layers[layer_nr]
962 for winding_nr, winding in layer.windings.items():
963 cable_type_curr = self.data.conductors[winding.conductor_name].cable.type
964 layers.windings[winding_nr] = dM.Block(conductor_name=winding.conductor_name, conductors_number=winding.conductors_number)
965 windings = layers.windings[winding_nr]
966 blk_list_current = list(winding.blocks.keys())
967 for block_nr, block in winding.blocks.items():
968 if block_nr in present_blocks:
969 blk_nr = str(block_nr)
970 windings.blocks[block_nr] = dM.BlockData(current_sign=block.current_sign)
971 hts = windings.blocks[block_nr].half_turns
972 is_first_blk = block_nr == list(winding.blocks.keys())[0]
973 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
974 if 'solenoid' in coil.type:
975 ht_list = (list(reversed(block.half_turns.keys()) if (layer_nr - 1) % 2 == 0 else list(block.half_turns.keys())))
976 elif coil.type == 'reversed-block-coil':
977 ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))
978 else:
979 ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))
980 hh = block.half_turns[ht_list[-1]].corners.bare
981 ll = block.half_turns[ht_list[0]].corners.bare
982 bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)
983 # Mid layer lines
984 mean_rad_current = (Func.points_distance([bc_current.oH.x, bc_current.oH.y], [coil.bore_center.x, coil.bore_center.y]) +
985 Func.points_distance([bc_current.oL.x, bc_current.oL.y], [coil.bore_center.x, coil.bore_center.y])) / 2
986 mid_ht_current_i = int(len(ht_list) / 2) if len(ht_list) % 2 == 0 else round(len(ht_list) / 2)
987 mid_ht_current = block.half_turns[ht_list[mid_ht_current_i - 1]].corners.insulated
988 concentric_coil = concentric_coils[(coil.bore_center.x, coil.bore_center.y)]
989 if layer_nr < len(pole.layers):
990 for winding_nr_next, winding_next in pole.layers[layer_nr + 1].windings.items():
991 if cable_type_curr == 'Rutherford' or\
992 (cable_type_curr in ['Mono', 'Ribbon'] and
993 list(pole.layers[layer_nr + 1].windings.keys()).index(winding_nr_next) == list(layer.windings.keys()).index(winding_nr)):
994 blk_list_next = list(winding_next.blocks.keys())
995 block_nr_next = blk_list_next[blk_list_current.index(block_nr)]
996 block_next = winding_next.blocks[block_nr_next]
997 _addMidLayerThinShellGroup(coil)
998 elif concentric_coil.index(coil_nr) + 1 < len(concentric_coil):
999 coil_nr_next = concentric_coil[concentric_coil.index(coil_nr) + 1]
1000 for pole_nr_next, pole_next in self.geom.coil.coils[coil_nr_next].poles.items():
1001 for layer_nr_next, layer_next in pole_next.layers.items():
1002 if layer_nr_next == 1:
1003 for winding_nr_next, winding_next in layer_next.windings.items():
1004 for block_nr_next, block_next in winding_next.blocks.items():
1005 _addMidLayerThinShellGroup(coil, mid_coil=True)
1006 else:
1007 blk_ins = windings.blocks[block_nr].insulation
1008 blk_ins.areas[blk_nr] = dM.Area()
1010 if 'solenoid' in coil.type:
1011 ht_items = (list(reversed(block.half_turns.items()) if layer_nr - 1 % 2 == 0 else list(block.half_turns.items())))
1012 elif coil.type == 'reversed-block-coil':
1013 ht_items = (block.half_turns.items() if not is_first_blk else reversed(block.half_turns.items()))
1014 else:
1015 ht_items = (block.half_turns.items() if is_first_blk else reversed(block.half_turns.items()))
1016 for halfTurn_nr, halfTurn in ht_items:
1017 ht_nr = str(halfTurn_nr)
1018 ht = halfTurn.corners.insulated
1019 hts.areas[ht_nr] = dM.Area()
1020 ht_b = halfTurn.corners.bare
1022 hts.points[ht_nr + 'ih'] = self.occ.addPoint(ht_b.iH.x, ht_b.iH.y, 0)
1023 hts.points[ht_nr + 'il'] = self.occ.addPoint(ht_b.iL.x, ht_b.iL.y, 0)
1024 hts.points[ht_nr + 'oh'] = self.occ.addPoint(ht_b.oH.x, ht_b.oH.y, 0)
1025 hts.points[ht_nr + 'ol'] = self.occ.addPoint(ht_b.oL.x, ht_b.oL.y, 0)
1027 hts.lines[ht_nr + 'i'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'il'])
1028 hts.lines[ht_nr + 'o'] = self.occ.addLine(hts.points[ht_nr + 'oh'], hts.points[ht_nr + 'ol'])
1029 hts.lines[ht_nr + 'l'] = self.occ.addLine(hts.points[ht_nr + 'il'], hts.points[ht_nr + 'ol'])
1030 hts.lines[ht_nr + 'h'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'oh'])
1032 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:
1033 intersection = {}
1034 # Create mid layer points and compute their angle to the x-axis
1035 for mid_lyr_type in ['current', 'previous']:
1036 for pnt1, pnt2, side in zip(
1037 [[ht_b.iH.x, ht_b.iH.y], [ht_b.iL.x, ht_b.iL.y]],
1038 [[ht_b.oH.x, ht_b.oH.y], [ht_b.oL.x, ht_b.oL.y]], ['h', 'l']):
1039 if (cable_type_curr in ['Mono', 'Ribbon'] and coil.type == 'cos-theta' and
1040 (layer_nr < len(pole.layers) and mid_lyr_type == 'current' or layer_nr > 1 and mid_lyr_type == 'previous')):
1041 pnts_input = pnt1 + pnt2
1042 elif coil.type == 'cos-theta' and (cable_type_curr == 'Rutherford' or cable_type_curr in ['Mono', 'Ribbon'] and\
1043 (layer_nr == len(pole.layers) and mid_lyr_type == 'current' or layer_nr == 1 and mid_lyr_type == 'previous')):
1044 pnts_input = Func.line_through_two_points(pnt1, pnt2)
1045 elif 'block-coil' in coil.type:
1046 pnts_input = pnt1
1047 intersect = {}
1048 if mid_lyr_type == 'current':
1049 # Current mid-layer
1050 for ts_name in ts_endpoints.keys():
1051 if blk_nr == ts_name[:ts_name.index('_')]:
1052 _addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)
1053 elif mid_lyr_type == 'previous':
1054 # Previous mid-layer
1055 if block_nr in next_blks_list:
1056 for ts_name in next_blks_list[block_nr]:
1057 _addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)
1058 for key, value in intersect.items():
1059 if key in intersection:
1060 intersection[key][side] = value
1061 else:
1062 intersection[key] = {side: value}
1064 # Search for half turns that face thin shells only partially
1065 def __create_aux_mid_layer_point(ss, points):
1066 mid_layer_ts_aux[key] = dM.Region()
1067 if 'block-coil' in coil.type:
1068 inter_pnt = [points[0], ts_endpoints[key][0][ss][1]]
1069 else:
1070 inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),
1071 [ts_endpoints[key][1], ts_endpoints[key][0][ss]], get_only_closest=True)[0]
1072 mid_layer_ts_aux[key].points[str(halfTurn_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)
1073 mid_layer_ts_aux[key].lines[blk_nr] = 0
1074 for key, value in intersection.items():
1075 first_blk, second_blk = key.split('_')
1076 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
1077 if value['h'] and not value['l']:
1078 __create_aux_mid_layer_point('lower', [ht_b.iL.x, ht_b.iL.y])
1079 elif value['l'] and not value['h']:
1080 __create_aux_mid_layer_point('higher', [ht_b.iH.x, ht_b.iH.y])
1081 else:
1082 relevant_blk = int(first_blk) if second_blk == blk_nr else int(second_blk)
1083 if layer_nr == len(pole.layers) and blk_nr == first_blk:
1084 lyr_blks = acw_order[coil_nr + 1].layers[1]
1085 elif layer_nr == 1 and blk_nr == second_blk:
1086 lyr_blks = acw_order[coil_nr - 1].layers[len(acw_order[coil_nr - 1].layers)]
1087 else:
1088 lyr_blks = acw_order[coil_nr].layers[layer_nr + (1 if first_blk == blk_nr else -1)]
1089 for nr, block_order in enumerate(lyr_blks):
1090 if block_order.block == relevant_blk:
1091 block_order_curr = block_order
1092 block_order_prev = lyr_blks[-1] if nr == 0 else lyr_blks[nr - 1]
1093 block_order_next = lyr_blks[0] if nr + 1 == len(lyr_blks) else lyr_blks[nr + 1]
1094 break
1095 if value['h'] and not value['l'] and block_order_curr.winding == block_order_prev.winding:
1096 __create_aux_mid_layer_point('lower', [[ht_b.iL.x, ht_b.iL.y], [ht_b.oL.x, ht_b.oL.y]])
1097 elif value['l'] and not value['h'] and block_order_curr.winding == block_order_next.winding:
1098 __create_aux_mid_layer_point('higher', [[ht_b.iH.x, ht_b.iH.y], [ht_b.oH.x, ht_b.oH.y]])
1099 else:
1100 blk_ins.points[ht_nr + 'ih'] = self.occ.addPoint(ht.iH.x, ht.iH.y, 0)
1101 blk_ins.points[ht_nr + 'il'] = self.occ.addPoint(ht.iL.x, ht.iL.y, 0)
1102 blk_ins.points[ht_nr + 'oh'] = self.occ.addPoint(ht.oH.x, ht.oH.y, 0)
1103 blk_ins.points[ht_nr + 'ol'] = self.occ.addPoint(ht.oL.x, ht.oL.y, 0)
1105 hts.areas[ht_nr].loop = self.occ.addCurveLoop(
1106 [hts.lines[ht_nr + 'i'], # inner
1107 hts.lines[ht_nr + 'l'], # lower
1108 hts.lines[ht_nr + 'o'], # outer
1109 hts.lines[ht_nr + 'h']]) # higher
1111 # Build wire order of the insulation lines of the current block
1112 if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
1113 ht_list = list(hts.areas.keys())
1114 ht_list.extend(list(reversed(ht_list))[1:])
1115 self.blk_ins_lines[block_nr] = ['l']
1116 for nr, ht_nr in enumerate(ht_list):
1117 if nr + 1 == winding.conductors_number: # end of first round
1118 self.blk_ins_lines[block_nr].extend([ht_nr + 'i', 'h', ht_nr + 'o'])
1119 else:
1120 if nr + 1 < winding.conductors_number: # within first round
1121 self.blk_ins_lines[block_nr].extend([ht_nr + 'i', ht_nr + 'i' + ht_list[nr + 1]])
1122 else: # within second round
1123 self.blk_ins_lines[block_nr].extend([ht_nr + 'o' + ht_list[nr - 1], ht_nr + 'o'])
1125 def constructInsulationGeometry(self):
1126 """
1127 Generates points, hyper lines, and curve loops for the coil insulations
1128 """
1129 def _createMidPoleLines(case, cnt=0):
1130 if 'block-coil' in geom_coil.type:
1131 if case == 'inner':
1132 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'])
1133 ordered_lines[group_nr].append(['mid_pole_' + case[0], (len(coil.layers) * 2) * 1e3 + 5e2, group.lines['mid_pole_' + case[0]]])
1134 else:
1135 group.lines['mid_pole_' + case[0]] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt_opposite[first_ht_prev + 'il'])
1136 ordered_lines[group_nr].append(['mid_pole_' + case[0], 0, group.lines['mid_pole_' + case[0]]])
1137 else:
1138 ht_curr = geom_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[
1139 block_order.block].half_turns[int(first_ht_curr)].corners.insulated
1140 ht_prev = geom_coil.poles[block_order_prev.pole].layers[layer_nr].windings[block_order_prev.winding].blocks[
1141 block_order_prev.block].half_turns[int(last_ht_prev)].corners.insulated
1142 pnt_curr = [ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y]
1143 pnt_prev = [ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y]
1144 if Func.points_distance(pnt_curr, pnt_prev) > 1e-6:
1145 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],
1146 [ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y],
1147 [ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y])
1148 ln_name = 'mid_pole_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]
1149 group.lines[ln_name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],
1150 self.occ.addPoint(correct_center[0], correct_center[1], 0),
1151 ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
1152 # self.occ.addLine(ins_pnt[first_ht_curr + case[0] + 'l'], ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
1153 cnt += 1 if case == 'inner' else -1
1154 ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])
1155 return cnt
1157 def _createMidWindingLines(case, cnt):
1158 name = 'mid_wind_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]
1159 # Create corrected center
1160 blk1 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order.block)][0]].layers[
1161 blks_info[str(block_order.block)][1]].windings[blks_info[str(block_order.block)][2]].blocks[int(str(block_order.block))]
1162 blk2 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order_prev.block)][0]].layers[
1163 blks_info[str(block_order_prev.block)][1]].windings[blks_info[str(block_order_prev.block)][2]].blocks[int(block_order_prev.block)]
1164 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
1165 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
1166 outer_center = Func.corrected_arc_center([self.md.geometries.coil.coils[coil_nr].bore_center.x,
1167 self.md.geometries.coil.coils[coil_nr].bore_center.y],
1168 [pnt1.x, pnt1.y], [pnt2.x, pnt2.y])
1169 group.lines[name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],
1170 self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt_opposite[last_ht_prev + case[0] + 'h'])
1171 cnt += 1 if case == 'inner' else -1
1172 ordered_lines[group_nr].append([name, cnt, group.lines[name]])
1173 return cnt
1175 def _createInnerOuterLines(case, cnt):
1176 # Create half turn lines
1177 idxs = [1, round(len(self.blk_ins_lines[block_order.block]) / 2), 1] if case == 'inner'\
1178 else [len(self.blk_ins_lines[block_order.block]) - 1, round(len(self.blk_ins_lines[block_order.block]) / 2), -1]
1179 lns = self.blk_ins_lines[block_order.block][idxs[0]:idxs[1]:idxs[2]]
1180 for ln_nr, ln_name in enumerate(lns):
1181 skip_cnt = False
1182 if ln_name[-1].isdigit():
1183 try:
1184 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + 'h'],
1185 ins_pnt[ln_name[ln_name.index(case[0]) + 1:] + case[0] + 'l'])
1186 except:
1187 skip_cnt = True
1188 next_line = lns[ln_nr + 1]
1189 pos = 'first' if next_line[:-1] == ln_name[:ln_name.index(case[0])] else 'second'
1190 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')
1191 elif ln_name[-1] in ['i', 'o']:
1192 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name + 'l'], ins_pnt[ln_name + 'h'])
1193 else:
1194 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + ln_name[-1]],
1195 ins_pnt[ln_name[ln_name.index(case[0]) + 1:-1] + case[0] + ln_name[-1]])
1196 if not skip_cnt:
1197 cnt += 1 if case == 'inner' else -1
1198 ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])
1199 return cnt
1201 def _computePointAngle(case):
1202 points_angles = pa_next if case == 'outer' else pa_prev
1203 current_ht_h = [current_ht.oH.x, current_ht.oH.y] if case == 'outer' else [current_ht.iH.x, current_ht.iH.y]
1204 if ht_nr == 0:
1205 current_ht_l = [current_ht.oL.x, current_ht.oL.y] if case == 'outer' else [current_ht.iL.x, current_ht.iL.y]
1206 if 'block-coil' in geom_coil.type: current_ht_l[1] = 1 if current_ht_l[1] > 0 else -1
1207 points_angles[str(block_order.block) + '_' + ht_name + 'l'] = Func.arc_angle_between_point_and_abscissa(current_ht_l, center)
1208 if ht_nr == len(ht_list) - 1:
1209 name = ht_name + 'h'
1210 coord = current_ht_h
1211 else: # for mid half turns, get the outer corner
1212 next_ht_ins = geom_hts[int(ht_list[ht_nr + 1])].corners.insulated
1213 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]
1214 condition = (Func.points_distance(current_ht_h, center) > Func.points_distance(next_ht, center))\
1215 if case == 'outer' else (Func.points_distance(current_ht_h, center) < Func.points_distance(next_ht, center))
1216 if condition:
1217 name = ht_name + 'h'
1218 coord = current_ht_h
1219 else:
1220 name = ht_list[ht_nr + 1] + 'l'
1221 coord = next_ht
1222 if 'block-coil' in geom_coil.type: coord[1] = 1 if coord[1] > 0 else -1
1223 points_angles[str(block_order.block) + '_' + name] = Func.arc_angle_between_point_and_abscissa(coord, center)
1225 ins = self.md.geometries.insulation
1226 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
1227 aux_coil = self.md.geometries.coil.coils[coil_nr]
1228 geom_coil = self.geom.coil.coils[coil_nr]
1229 groups = len(geom_coil.poles)
1230 count = {}
1231 ordered_lines = {}
1232 points_angle = {}
1233 blks_info = {}
1234 ending_line = {}
1235 center = [geom_coil.bore_center.x, geom_coil.bore_center.y]
1236 if coil_nr not in ins.coils:
1237 ins.coils[coil_nr] = dM.InsulationGroup()
1238 ins_groups = ins.coils[coil_nr].group
1239 for layer_nr, layer in coil.layers.items():
1240 group_nr = 1
1241 wnd_nr = len(aux_coil.poles[1].layers[layer_nr].windings)
1242 ordered_layer = layer[wnd_nr:] + layer[:wnd_nr] if layer[0].pole != layer[-1].pole else layer
1243 for nr, block_order in enumerate(ordered_layer):
1244 blks_info[str(block_order.block)] = [block_order.pole, layer_nr, block_order.winding]
1245 # Get previous block in anticlockwise order
1246 block_order_prev = ordered_layer[-1] if nr == 0 else ordered_layer[nr - 1]
1247 # Update insulation group
1248 if block_order.winding == block_order_prev.winding:
1249 group_nr = group_nr + 1 if group_nr < groups else 1
1250 # Initialize dicts
1251 if group_nr not in ins_groups:
1252 ins_groups[group_nr] = dM.InsulationRegion()
1253 points_angle[group_nr] = {}
1254 ordered_lines[group_nr] = []
1255 count[group_nr] = [0, (len(coil.layers) + 1) * 1e3]
1256 group = ins_groups[group_nr].ins
1257 ins_groups[group_nr].blocks.append([block_order.pole, layer_nr, block_order.winding, block_order.block])
1258 # Find the wedge
1259 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
1260 for wdg, blk in self.md.geometries.wedges.coils[coil_nr].layers[layer_nr].block_prev.items():
1261 if blk == block_order_prev.block:
1262 ins_groups[group_nr].wedges.append([layer_nr, wdg])
1263 break
1264 if layer_nr < len(coil.layers):
1265 mid_layer_next = str(layer_nr) + '_' + str(layer_nr + 1)
1266 if mid_layer_next not in points_angle[group_nr]:
1267 points_angle[group_nr][mid_layer_next] = {}
1268 pa_next = points_angle[group_nr][mid_layer_next]
1269 if layer_nr > 1:
1270 mid_layer_prev = str(layer_nr - 1) + '_' + str(layer_nr)
1271 pa_prev = points_angle[group_nr][mid_layer_prev]
1272 # Get point tags of insulation
1273 ins_pnt = aux_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[
1274 block_order.block].insulation.points
1275 # Get relevant info for line names
1276 first_ht_curr = self.blk_ins_lines[block_order.block][1][:-1]
1277 last_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[
1278 layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[-1]
1279 ins_pnt_opposite = aux_coil.poles[block_order_prev.pole].layers[
1280 layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].insulation.points
1281 if 'cos-theta' == geom_coil.type:
1282 # Create lower and higher angle lines
1283 if block_order.winding == block_order_prev.winding:
1284 group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
1285 ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])
1286 ending_line[group_nr - 1 if group_nr > 1 else groups] =\
1287 [ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]
1288 # Create inner lines of insulation group
1289 if layer_nr == 1:
1290 if block_order.pole != block_order_prev.pole:
1291 count[group_nr][0] = _createMidPoleLines('inner', count[group_nr][0])
1292 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
1293 count[group_nr][0] = _createMidWindingLines('inner', count[group_nr][0])
1294 count[group_nr][0] = _createInnerOuterLines('inner', count[group_nr][0])
1295 # Create outer lines of insulation group
1296 if layer_nr == len(coil.layers):
1297 if block_order.pole != block_order_prev.pole:
1298 count[group_nr][1] = _createMidPoleLines('outer', count[group_nr][1])
1299 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:
1300 count[group_nr][1] = _createMidWindingLines('outer', count[group_nr][1])
1301 count[group_nr][1] = _createInnerOuterLines('outer', count[group_nr][1])
1302 elif 'block-coil' in geom_coil.type:
1303 last_ht_curr = self.blk_ins_lines[block_order.block][self.blk_ins_lines[block_order.block].index('h') - 1][:-1]
1304 first_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[layer_nr].windings[
1305 block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[0]
1306 # Create lower and higher angle lines
1307 if block_order.winding == block_order_prev.winding:
1308 group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
1309 ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 4 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])
1310 ending_line[group_nr - 1 if group_nr > 1 else groups] =\
1311 [ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]
1312 group.lines[str(layer_nr) + 'bh'] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt[last_ht_curr + 'oh'])
1313 ordered_lines[group_nr].append([str(layer_nr) + 'bh', (len(coil.layers) * 2 + layer_nr) * 1e3, group.lines[str(layer_nr) + 'bh']])
1314 # Create inner lines of insulation group
1315 if block_order.pole != block_order_prev.pole:
1316 if layer_nr == 1:
1317 _createMidPoleLines('inner')
1318 _createMidPoleLines('outer')
1319 group.lines[str(layer_nr) + 'bl'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])
1320 ordered_lines[group_nr].append([str(layer_nr) + 'bl', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'bl']])
1321 # Create outer lines of insulation group
1322 if layer_nr == len(coil.layers):
1323 count[group_nr][1] = _createInnerOuterLines(
1324 'outer', (len(coil.layers) * 4 - layer_nr + 1) * 1e3 if block_order.winding == block_order_prev.winding else (len(coil.layers) + 1) * 1e3)
1325 # Store info about the angle of each point in between layers
1326 ht_list = list(aux_coil.poles[block_order.pole].layers[
1327 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns.areas.keys())
1328 geom_hts = geom_coil.poles[block_order.pole].layers[
1329 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns
1330 for ht_nr, ht_name in enumerate(ht_list): # half turns in anticlockwise order
1331 current_ht = geom_hts[int(ht_name)].corners.insulated
1332 if layer_nr < len(coil.layers): # if it's not the last layer, fetch all outer corners angles
1333 _computePointAngle('outer')
1334 if layer_nr > 1: # if it's not the first layer, fetch all inner corners angles
1335 _computePointAngle('inner')
1336 # Create closing lines
1337 for grp_nr, grp in ending_line.items():
1338 ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h'] = self.occ.addLine(grp[0], grp[1])
1339 ordered_lines[grp_nr].append([str(layer_nr) + 'h', layer_nr * 1e3, ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h']])
1340 # Create lines connecting different layers and generate closed loops
1341 for group_nr, group in points_angle.items():
1342 ins_group = ins_groups[group_nr].ins
1343 for mid_l_name, mid_l in group.items():
1344 first_layer = mid_l_name[:mid_l_name.index('_')]
1345 # Correct angles if the group crosses the abscissa
1346 max_angle = max(mid_l.values())
1347 max_diff = max_angle - min(mid_l.values())
1348 if max_diff > np.pi:
1349 for pnt_name, angle in mid_l.items():
1350 if angle < max_diff / 2:
1351 mid_l[pnt_name] = angle + max_angle
1352 # Order points according to angle
1353 ordered_pnts = [[pnt_name, angle] for pnt_name, angle in mid_l.items()]
1354 ordered_pnts.sort(key=lambda x: x[1])
1355 ordered_names = [x[0] for x in ordered_pnts]
1356 for case in ['beg', 'end']:
1357 past_blocks = []
1358 sides = ['l', 'o', 'h', 'l'] if case == 'beg' else ['h', 'i', 'l', 'h']
1359 # count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2
1360 for i in range(2 if 'block-coil' in geom_coil.type else 1):
1361 count = int(first_layer) * 1e3 + 5e2 if i == 0 else (len(coil.layers) * 2 + int(first_layer)) * 1e3 + 5e2
1362 if case == 'beg':
1363 pnt_position = 0 if i == 0 else int(len(ordered_names) / 2)
1364 else:
1365 pnt_position = -1 if i == 0 else int(len(ordered_names) / 2 - 1)
1366 first_block = ordered_names[pnt_position][:ordered_names[pnt_position].index('_')] # ordered_pnts[pnt_position][0][:ordered_pnts[pnt_position][0].index('_')] #
1367 ordered_search_names = ordered_names[pnt_position::1 if case == 'beg' else -1]
1368 for nr, pnt in enumerate(ordered_search_names[1:], 1): # enumerate(ordered_names if case == 'beg' else reversed(ordered_names)): #
1369 current_blk = pnt[:pnt.index('_')]
1370 ins_pnt = aux_coil.poles[blks_info[current_blk][0]].layers[blks_info[current_blk][1]].windings[
1371 blks_info[current_blk][2]].blocks[int(current_blk)].insulation.points
1372 prev_pnt = ordered_search_names[nr - 1] # ordered_pnts[nr - 1 if case == 'beg' else - nr][0] #
1373 prev_blk = prev_pnt[:prev_pnt.index('_')]
1374 start_pnt_name = prev_pnt[prev_pnt.index('_') + 1:-1] + ('o' if str(blks_info[prev_blk][1]) == first_layer else 'i')
1375 ins_pnt_prev = aux_coil.poles[blks_info[prev_blk][0]].layers[blks_info[prev_blk][1]].windings[
1376 blks_info[prev_blk][2]].blocks[int(prev_blk)].insulation.points
1377 # Create lines when you find the first edge belonging to a block of the opposite layer
1378 if blks_info[current_blk][1] != blks_info[first_block][1]:
1379 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')
1380 pnt_tag_name_opposite = start_pnt_name + ('l' if prev_pnt[-1] == 'l' else 'h')
1381 opp_blk_ins_lines = self.blk_ins_lines[int(prev_blk)]
1382 indexes = [opp_blk_ins_lines.index(start_pnt_name) + (1 if prev_pnt[-1] == sides[0] else 0),
1383 len(opp_blk_ins_lines) if case == 'beg' else opp_blk_ins_lines.index('h'), 1] if start_pnt_name[-1] == sides[1]\
1384 else [opp_blk_ins_lines.index(start_pnt_name) - (1 if prev_pnt[-1] == sides[0] else 0),
1385 0 if case == 'beg' else opp_blk_ins_lines.index('h'), -1]
1386 if case == 'beg':
1387 if i == 0:
1388 count = (len(coil.layers) * (4 if 'block-coil' in geom_coil.type else 2) - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])
1389 else:
1390 count = (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])
1391 else:
1392 count += 1 + abs(indexes[0] - indexes[1])
1393 # Create all remaining lines of the current layer block
1394 for line_name in opp_blk_ins_lines[indexes[0]:indexes[1]:indexes[2]]:
1395 if 'block-coil' in geom_coil.type:
1396 if not line_name[-1].isdigit():
1397 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])
1398 count += 1 if (case == 'beg' and i == 1) or (case == 'end' and i == 0) else -1
1399 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
1400 else:
1401 if line_name[-1].isdigit():
1402 ins_group.lines[line_name] = self.occ.addLine(
1403 ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'],
1404 ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'])
1405 else:
1406 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])
1407 count += 1 if case == 'beg' else -1 # if start_pnt_name[-1] == sides[1] else 1
1408 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
1409 # Create mid layer line
1410 if 'block-coil' in geom_coil.type:
1411 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])
1412 else:
1413 count_rest = -abs(indexes[0] - indexes[1]) if case == 'beg' else 1 + abs(indexes[0] - indexes[1])
1414 line_name = 'mid_layer_' + mid_l_name + ('b' if i == 1 else '') + ('_l' if case == 'beg' else '_h')
1415 ins_group.lines[line_name] = self.occ.addLine(ins_pnt[pnt_tag_name], ins_pnt_prev[pnt_tag_name_opposite])
1416 ordered_lines[group_nr].append([line_name, count + count_rest, ins_group.lines[line_name]])
1417 break
1418 # Create all edges of the first block sticking out completely todo: might have to be extended to multiple blocks
1419 if current_blk != first_block and current_blk not in past_blocks:
1420 def __createWedgeInsulation(cnt):
1421 # Create the line connecting the blocks (where a wedge is)
1422 line_name = self.blk_ins_lines[int(current_blk)][
1423 (-1 if start_pnt_name[-1] == 'o' else 1) if case == 'beg'
1424 else (round(len(self.blk_ins_lines[int(current_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1))]
1425 line_name_prev = self.blk_ins_lines[int(prev_blk)][
1426 (round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1)) if case == 'beg'
1427 else (-1 if start_pnt_name[-1] == 'o' else 1)]
1428 # Create corrected center
1429 blk1 = geom_coil.poles[blks_info[prev_blk][0]].layers[
1430 blks_info[prev_blk][1]].windings[blks_info[prev_blk][2]].blocks[int(prev_blk)]
1431 blk2 = geom_coil.poles[blks_info[current_blk][0]].layers[
1432 blks_info[current_blk][1]].windings[blks_info[current_blk][2]].blocks[int(current_blk)]
1433 pnt1 = blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oH if case == 'beg'\
1434 else blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oL
1435 pnt2 = blk2.half_turns[int(line_name[:-1])].corners.insulated.oL if case == 'beg'\
1436 else blk2.half_turns[int(line_name[:-1])].corners.insulated.oH
1437 outer_center = Func.corrected_arc_center([aux_coil.bore_center.x, aux_coil.bore_center.y],
1438 [pnt2.x, pnt2.y] if case == 'beg' else [pnt1.x, pnt1.y],
1439 [pnt1.x, pnt1.y] if case == 'beg' else [pnt2.x, pnt2.y])
1440 ins_group.lines[line_name_prev + line_name] =\
1441 self.occ.addCircleArc(ins_pnt_prev[line_name_prev + sides[2]],
1442 self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt[line_name + sides[3]])
1443 ordered_lines[group_nr].append([line_name_prev + line_name, cnt, ins_group.lines[line_name_prev + line_name]])
1445 count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2
1446 past_blocks.append(current_blk)
1447 indexes = [round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + 1,
1448 len(self.blk_ins_lines[int(prev_blk)])] if str(blks_info[prev_blk][1]) == first_layer\
1449 else [1, round(len(self.blk_ins_lines[int(prev_blk)]) / 2)]
1450 if case == 'beg':
1451 count += 1
1452 __createWedgeInsulation(count)
1453 lines = self.blk_ins_lines[int(prev_blk)][indexes[0]:indexes[1]]
1454 side = 'o' if str(blks_info[prev_blk][1]) == first_layer else 'i'
1455 for line_nr, line_name in enumerate(lines):
1456 skip_count = False
1457 if line_name[-1].isdigit():
1458 try:
1459 ins_group.lines[line_name] =\
1460 self.occ.addLine(ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'],
1461 ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'])
1462 except: # points are too close to each other
1463 skip_count = True
1464 next_line = lines[line_nr + 1]
1465 pnt1, pnt2 = line_name.split(side)
1466 pos = 'first' if next_line[:-1] == pnt1 else 'second'
1467 lines[line_nr + 1] = next_line + (pnt2 + 'l' if pos == 'first' else pnt1 + 'h')
1468 elif line_name[-1] in ['i', 'o']:
1469 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'h'], ins_pnt_prev[line_name + 'l'])
1470 else:
1471 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name[:line_name.index(side)] + side + line_name[-1]],
1472 ins_pnt_prev[line_name[line_name.index(side) + 1:-1] + side + line_name[-1]])
1473 if not skip_count:
1474 count += 1 # if start_pnt_name[-1] == sides[1] else -1
1475 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])
1476 if case == 'end':
1477 count += 1
1478 __createWedgeInsulation(count)
1480 # Generate closed loops
1481 ordered_lines[group_nr].sort(key=lambda x: x[1])
1482 area_name = str((coil_nr - 1) * len(ins_groups) + group_nr)
1483 ins_group.areas[area_name] = dM.Area()
1484 if len(points_angle) == 1:
1485 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]]
1486 if 'i' in line and line[0].isdigit() or '_i' in line]))
1487 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]]
1488 if 'o' in line and line[0].isdigit() or '_o' in line])
1489 else:
1490 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]]])
1492 def constructThinShells(self, with_wedges):
1493 ins_th = self.md.geometries.thin_shells.ins_thickness
1494 mid_pole_ts = self.md.geometries.thin_shells.mid_poles
1495 mid_winding_ts = self.md.geometries.thin_shells.mid_windings
1496 mid_turn_ts = self.md.geometries.thin_shells.mid_turn_blocks
1497 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht
1498 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux
1500 # Create mid-pole and mid-turn thin shells
1501 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
1502 for layer_nr, layer in coil.layers.items():
1503 for nr, blk_order in enumerate(layer):
1504 block = self.geom.coil.coils[coil_nr].poles[blk_order.pole].layers[
1505 layer_nr].windings[blk_order.winding].blocks[blk_order.block]
1506 ht_list = list(self.md.geometries.coil.coils[coil_nr].poles[blk_order.pole].layers[
1507 layer_nr].windings[blk_order.winding].blocks[blk_order.block].half_turns.areas.keys())
1508 # Create mid-pole and mid-winding thin shells
1509 blk_index_next = nr + 1 if nr + 1 < len(layer) else 0
1510 block_order_next = layer[blk_index_next]
1511 block_next = self.geom.coil.coils[coil_nr].poles[block_order_next.pole].layers[
1512 layer_nr].windings[block_order_next.winding].blocks[block_order_next.block]
1513 ht_list_next = list(self.md.geometries.coil.coils[coil_nr].poles[block_order_next.pole].layers[
1514 layer_nr].windings[block_order_next.winding].blocks[block_order_next.block].half_turns.areas.keys())
1515 ht_last = int(ht_list[-1])
1516 ht_next_first = int(ht_list_next[0])
1517 iH = [block.half_turns[ht_last].corners.bare.iH.x, block.half_turns[ht_last].corners.bare.iH.y]
1518 iL = [block_next.half_turns[ht_next_first].corners.bare.iL.x, block_next.half_turns[ht_next_first].corners.bare.iL.y]
1519 oH = [block.half_turns[ht_last].corners.bare.oH.x, block.half_turns[ht_last].corners.bare.oH.y]
1520 oL = [block_next.half_turns[ht_next_first].corners.bare.oL.x, block_next.half_turns[ht_next_first].corners.bare.oL.y]
1521 ts_name = str(blk_order.block) + '_' + str(block_order_next.block)
1522 for ts, th, condition in zip([mid_pole_ts, mid_winding_ts], [ins_th.mid_pole, ins_th.mid_winding],
1523 # ['_ly' + str(layer_nr), '_wd' + str(blk_order.winding) + '_wd' + str(block_order_next.winding)],
1524 [self.geom.coil.coils[coil_nr].type == 'cos-theta' and block_order_next.pole != blk_order.pole,
1525 (not with_wedges or not self.geom.wedges) and self.geom.coil.coils[coil_nr].type in
1526 ['cos-theta', 'common-block-coil'] and block_order_next.pole == blk_order.pole and block_order_next.winding != blk_order.winding]):
1527 if condition:
1528 ts[ts_name] = dM.Region()
1529 ts[ts_name].points['i'] = self.occ.addPoint((iH[0] + iL[0]) / 2, (iH[1] + iL[1]) / 2, 0)
1530 ts[ts_name].points['o'] = self.occ.addPoint((oH[0] + oL[0]) / 2, (oH[1] + oL[1]) / 2, 0)
1531 ts[ts_name].lines[str(ht_last) + '_' + str(ht_next_first)] =\
1532 self.occ.addLine(ts[ts_name].points['i'], ts[ts_name].points['o'])
1533 # Get insulation thickness
1534 th[ts_name] = Func.sig_dig((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)
1535 # if 'cl' + str(coil_nr) + th_name not in th:
1536 # th['cl' + str(coil_nr) + th_name] = float((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)
1537 # Create mid-turn thin shells
1538 mid_turn_ts[str(blk_order.block)] = dM.Region()
1539 ts = mid_turn_ts[str(blk_order.block)]
1540 for nr_ht, ht in enumerate(ht_list[:-1]):
1541 line_name = ht + '_' + ht_list[nr_ht + 1]
1542 current_ht = block.half_turns[int(ht)].corners.bare
1543 next_ht = block.half_turns[int(ht_list[nr_ht + 1])].corners.bare
1544 mid_inner = [(current_ht.iH.x + next_ht.iL.x) / 2, (current_ht.iH.y + next_ht.iL.y) / 2]
1545 mid_outer = [(current_ht.oH.x + next_ht.oL.x) / 2, (current_ht.oH.y + next_ht.oL.y) / 2]
1546 mid_length = Func.points_distance(mid_inner, mid_outer)
1547 mid_line = Func.line_through_two_points(mid_inner, mid_outer)
1548 points = {'inner': list, 'outer': list}
1549 for case, current_h, current_l, next_h, next_l, mid_point in zip(
1550 ['inner', 'outer'], [current_ht.iH, current_ht.oH], [current_ht.iL, current_ht.oL],
1551 [next_ht.iH, next_ht.oH], [next_ht.iL, next_ht.oL], [mid_outer, mid_inner]):
1552 current_line = Func.line_through_two_points([current_h.x, current_h.y], [current_l.x, current_l.y])
1553 next_line = Func.line_through_two_points([next_h.x, next_h.y], [next_l.x, next_l.y])
1554 current_intersect = Func.intersection_between_two_lines(mid_line, current_line)
1555 next_intersect = Func.intersection_between_two_lines(mid_line, next_line)
1556 points[case] = current_intersect if Func.points_distance(
1557 current_intersect, mid_point) < mid_length else next_intersect
1558 ts.points[line_name + '_i'] = self.occ.addPoint(points['inner'][0], points['inner'][1], 0)
1559 ts.points[line_name + '_o'] = self.occ.addPoint(points['outer'][0], points['outer'][1], 0)
1560 ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])
1562 # Create mid-layer thin shells
1563 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]
1564 for ts_name, ts in mid_layer_ts.items():
1565 # Order mid-layer thin shell points according to their angle with respect to the x-axis to generate lines
1566 blk1, blk2 = ts_name.split('_')
1567 max_angle = max(ts.point_angles.values())
1568 max_diff = max_angle - min(ts.point_angles.values())
1569 if max_diff > np.pi:
1570 for pnt_name, angle in ts.point_angles.items():
1571 if angle < max_diff / 2:
1572 ts.point_angles[pnt_name] = angle + max_angle
1573 ordered_pnts = [[pnt_name, ts.point_angles[pnt_name], pnt] for pnt_name, pnt in ts.mid_layers.points.items()]
1574 ordered_pnts.sort(key=lambda x: x[1])
1575 for nr, pnt in enumerate(ordered_pnts[:-1]):
1576 pnt_current = pnt[0]
1577 pnt_next = ordered_pnts[nr + 1][0]
1578 if ((pnt_current[-1] == 'l' and pnt_next[-1] == 'h' and ts_name not in block_coil_mid_pole_list) or # cos-theta
1579 (ts_name in block_coil_mid_pole_list and
1580 ((pnt_current[-1] == pnt_next[-1] == 'h' and block_coil_mid_pole_list.index(ts_name) == 0) or # assumes a dipole block-coil
1581 (pnt_current[-1] == pnt_next[-1] == 'l' and block_coil_mid_pole_list.index(ts_name) == 1) or # assumes a dipole block-coil
1582 (pnt_current[:-1] == pnt_next[:-1])))):
1583 if pnt_current[:-1] == pnt_next[:-1]:
1584 relevant_blk = blk2 if int(pnt_current[:-1]) in ts.half_turn_lists[blk1] else blk1
1585 if nr > 0:
1586 iter_nr = nr - 1
1587 while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr -= 1
1588 line_name = ordered_pnts[iter_nr][0][:-1] + '_' + pnt_current[:-1]
1589 else:
1590 if len(ordered_pnts) == 2: # todo: get right ht from relevant_blk for 1-ht blocks
1591 line_name = pnt_current[:-1] + '_' + str(ts.half_turn_lists[relevant_blk][0])
1592 else:
1593 iter_nr = nr + 1
1594 while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr += 1
1595 line_name = pnt_current[:-1] + '_' + ordered_pnts[iter_nr][0][:-1]
1596 else:
1597 line_name = pnt_current[:-1] + '_' + pnt_next[:-1]
1598 ts.mid_layers.lines[line_name] = self.occ.addLine(pnt[2], ordered_pnts[nr + 1][2])
1599 if ts_name in mid_layer_ts_aux:
1600 aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[0]
1601 other_pnt = ordered_pnts[0 if aux_pnt[-1] == 'l' else -1]
1602 other_pnt_coord = gmsh.model.getValue(0, other_pnt[2], [])[:2] # needs to be a new point
1603 mid_layer_ts_aux[ts_name].points[other_pnt[0]] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)
1604 line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]
1605 try:
1606 mid_layer_ts_aux[ts_name].lines[line_name] = \
1607 self.occ.addLine(mid_layer_ts_aux[ts_name].points[aux_pnt], mid_layer_ts_aux[ts_name].points[other_pnt[0]])
1608 except:
1609 mid_layer_ts_aux[ts_name].lines.pop(line_name)
1611 # Create wedge-to-block and block-to-wedge lines
1612 for wdg_ts in [self.md.geometries.thin_shells.mid_layers_wdg_to_ht, self.md.geometries.thin_shells.mid_layers_ht_to_wdg]:
1613 for ts_name, ts in wdg_ts.items():
1614 pnt_list = list(ts.points.keys())
1615 for nr, pnt in enumerate(pnt_list[:-1]):
1616 if pnt[-1] == 'l' and pnt_list[nr + 1][-1] == 'h':
1617 ts.lines[pnt[:-1] + '_' + pnt_list[nr + 1][:-1]] = self.occ.addLine(ts.points[pnt], ts.points[pnt_list[nr + 1]])
1618 if ts_name in mid_layer_ts_aux:
1619 aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[
1620 1 if list(mid_layer_ts_aux[ts_name].points.keys()).index('center') == 0 else 0]
1621 other_pnt = pnt_list[0 if aux_pnt[-1] == 'l' else -1]
1622 other_pnt_coord = gmsh.model.getValue(0, ts.points[other_pnt], [])[:2] # needs to be a new point
1623 mid_layer_ts_aux[ts_name].points[other_pnt] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)
1624 line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]
1625 mid_layer_ts_aux[ts_name].lines[line_name] = self.occ.addCircleArc(
1626 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])
1628 # Create wedge-to-wedge lines
1629 for ts_nr, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():
1630 ts.lines[ts_nr] = self.occ.addCircleArc(ts.points['beg'], ts.points['center'], ts.points[list(ts.points.keys())[-1]])
1632 # Create mid wedge-turn lines
1633 mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn
1634 for ts_nr, ts in mid_turn_ts.items():
1635 line_name = list(ts.points.keys())[0][:-2]
1636 ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])
1638 # Get insulation thickness
1639 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
1640 geom_coil = self.geom.coil.coils[coil_nr]
1641 # Get block-coil mid-pole thickness
1642 if coil_nr in self.block_coil_mid_pole_blks:
1643 for blk_orders in self.block_coil_mid_pole_blks[coil_nr]:
1644 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
1645 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
1646 ins_th.mid_layer[str(blk_orders[0].block) + '_' + str(blk_orders[1].block)] = Func.sig_dig(abs(block_y - block_next_y))
1647 # 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
1648 for layer_nr, layer in coil.layers.items():
1649 for blk_order in layer:
1650 for ts_name, ts in mid_layer_ts.items():
1651 blk1, blk2 = ts_name.split('_')
1652 if blk1 == str(blk_order.block) and ts_name not in block_coil_mid_pole_list:
1653 block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]
1654 if layer_nr < len(coil.layers):
1655 for blk_order_next in coil.layers[layer_nr + 1]:
1656 if blk_order_next.block == int(blk2):
1657 block_next = geom_coil.poles[blk_order_next.pole].layers[layer_nr + 1].windings[blk_order_next.winding].blocks[int(blk2)]
1658 break
1659 else:
1660 for blk_order_next in self.md.geometries.coil.anticlockwise_order.coils[coil_nr + 1].layers[1]:
1661 if blk_order_next.block == int(blk2):
1662 block_next = self.geom.coil.coils[coil_nr + 1].poles[blk_order_next.pole].layers[1].windings[blk_order_next.winding].blocks[int(blk2)]
1663 break
1664 distances = []
1665 lines = list(ts.mid_layers.lines.keys())
1666 for line_name in [lines[0], lines[-1]]:
1667 ht_1, ht_2 = int(line_name[:line_name.index('_')]), int(line_name[line_name.index('_') + 1:])
1668 ht_char = {'low_p': ht_1, 'high_p': ht_2,
1669 'current': ht_1 if ht_1 in ts.half_turn_lists[blk1] else ht_2,
1670 'next': ht_2 if ht_1 in ts.half_turn_lists[blk1] else ht_1}
1671 hts = {'current': block.half_turns[ht_char['current']].corners.bare,
1672 'next': block_next.half_turns[ht_char['next']].corners.bare}
1673 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],
1674 'high_p': [hts['current'].oH, hts['current'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['next'].iH, hts['next'].oH],
1675 'low_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['low_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH],
1676 'high_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH]}
1677 low_line = Func.line_through_two_points([hts_p['low_p'][0].x, hts_p['low_p'][0].y],
1678 [hts_p['low_p'][1].x, hts_p['low_p'][1].y])
1679 high_line = Func.line_through_two_points([hts_p['high_p'][0].x, hts_p['high_p'][0].y],
1680 [hts_p['high_p'][1].x, hts_p['high_p'][1].y])
1681 distances.extend([Func.points_distance([hts_p['low_p'][0].x, hts_p['low_p'][0].y], Func.intersection_between_two_lines(
1682 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]))),
1683 Func.points_distance([hts_p['high_p'][0].x, hts_p['high_p'][0].y], Func.intersection_between_two_lines(
1684 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])))])
1685 ins_th.mid_layer[ts_name] = Func.sig_dig(min(distances))
1686 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]):
1687 for ts_name, ts in wdg_ts.items():
1688 wdg, blk = ts_name.split('_')
1689 if blk == str(blk_order.block):
1690 block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]
1691 wedge = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + (1 if ts_type == 1 else -1)].wedges[int(wdg[1:])]
1692 pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]
1693 pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]
1694 pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]
1695 pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]
1696 low_line = Func.line_through_two_points(pnt_il, pnt_ol)
1697 high_line = Func.line_through_two_points(pnt_ih, pnt_oh)
1698 el1_l, el2_l = list(ts.lines.keys())[0].split('_')
1699 ht_l = block.half_turns[int(el1_l) if el2_l == wdg else int(el2_l)].corners.bare
1700 el1_h, el2_h = list(ts.lines.keys())[-1].split('_')
1701 ht_h = block.half_turns[int(el1_h) if el2_h == wdg else int(el2_h)].corners.bare
1702 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\
1703 else Func.line_through_two_points([ht_l.oL.x, ht_l.oL.y], [ht_l.oH.x, ht_l.oH.y])
1704 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 \
1705 else Func.line_through_two_points([ht_h.oL.x, ht_h.oL.y], [ht_h.oH.x, ht_h.oH.y])
1706 ins_th.mid_layer[ts_name] = Func.sig_dig(
1707 (Func.points_distance(pnt_ol if ts_type == 0 else pnt_il, Func.intersection_between_two_lines(low_line, opp_line_l)) +
1708 Func.points_distance(pnt_oh if ts_type == 0 else pnt_ih, Func.intersection_between_two_lines(high_line, opp_line_h))) / 2)
1710 for coil_nr, coil in self.md.geometries.wedges.coils.items():
1711 # 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
1712 for layer_nr, layer in coil.layers.items():
1713 for wedge_nr, wedge in layer.wedges.items():
1714 for ts_name, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():
1715 wdg1, wdg2 = ts_name[1:ts_name.index('_')], ts_name[ts_name.index('_') + 2:]
1716 if wdg1 == str(wedge_nr):
1717 wedge_next = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + 1].wedges[int(wdg2)]
1718 # pnt_il_next = gmsh.model.getValue(0, wedge_next.points['il'], [])[:2]
1719 # pnt_ih_next = gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2]
1720 pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]
1721 pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]
1722 pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]
1723 pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]
1724 low_line = Func.line_through_two_points(pnt_il, pnt_ol)
1725 high_line = Func.line_through_two_points(pnt_ih, pnt_oh)
1726 opp_line = Func.line_through_two_points(gmsh.model.getValue(0, wedge_next.points['il'], [])[:2],
1727 gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2])
1728 ins_th.mid_layer[ts_name] = Func.sig_dig(
1729 (Func.points_distance(pnt_ol, Func.intersection_between_two_lines(low_line, opp_line)) +
1730 Func.points_distance(pnt_oh, Func.intersection_between_two_lines(high_line, opp_line))) / 2)
1732 def buildDomains(self, run_type, symmetry):
1733 """
1734 Generates plane surfaces from the curve loops
1735 """
1736 iron = self.geom.iron
1737 gm = self.md.geometries
1738 with_iron_yoke = self.data.magnet.geometry.electromagnetics.with_iron_yoke if run_type == 'EM'\
1739 else self.data.magnet.geometry.thermal.with_iron_yoke
1740 with_wedges = self.data.magnet.geometry.electromagnetics.with_wedges if run_type == 'EM' \
1741 else self.data.magnet.geometry.thermal.with_wedges
1743 # Build iron yoke domains
1744 if with_iron_yoke:
1745 for quadrant, qq in gm.iron.quadrants.items():
1746 for area_name, area in qq.areas.items():
1747 build = True
1748 loops = [area.loop]
1749 for hole_key, hole in iron.hyper_holes.items():
1750 if area_name == hole.areas[1]:
1751 loops.append(qq.areas[hole.areas[0]].loop)
1752 elif area_name == hole.areas[0]: # or iron.hyper_areas[area_name].material == 'BH_air':
1753 build = False
1754 if build:
1755 area.surface = self.occ.addPlaneSurface(loops)
1756 # Group areas per material type
1757 self.md.domains.groups_entities.iron[iron.hyper_areas[area_name].material].append(area.surface)
1759 # Build coil domains
1760 for coil_nr, coil in gm.coil.coils.items():
1761 for pole_nr, pole in coil.poles.items():
1762 for layer_nr, layer in pole.layers.items():
1763 for winding_nr, winding in layer.windings.items():
1764 for block_key, block in winding.blocks.items():
1765 for area_name, area in block.half_turns.areas.items():
1766 area.surface = self.occ.addPlaneSurface([area.loop])
1768 # Build wedges domains
1769 if with_wedges:
1770 for coil_nr, coil in gm.wedges.coils.items():
1771 for layer_nr, layer in coil.layers.items():
1772 for wedge_nr, wedge in layer.wedges.items():
1773 wedge.areas[str(wedge_nr)].surface = self.occ.addPlaneSurface([wedge.areas[str(wedge_nr)].loop])
1775 # Build insulation domains
1776 if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
1777 for coil_nr, coil in gm.insulation.coils.items():
1778 for group_nr, group in coil.group.items():
1779 holes = []
1780 for blk in group.blocks:
1781 holes.extend([ht.loop for ht_nr, ht in gm.coil.coils[
1782 coil_nr].poles[blk[0]].layers[blk[1]].windings[blk[2]].blocks[blk[3]].half_turns.areas.items()])
1783 for wdg in group.wedges:
1784 holes.extend([wedge.loop for wedge_nr, wedge in gm.wedges.coils[
1785 coil_nr].layers[wdg[0]].wedges[wdg[1]].areas.items()])
1786 if len(group.ins.areas) == 1:
1787 for area_name, area in group.ins.areas.items():
1788 area.surface = self.occ.addPlaneSurface([area.loop] + holes)
1789 else:
1790 for area_name, area in group.ins.areas.items():
1791 if area_name.isdigit():
1792 area.surface = self.occ.addPlaneSurface([area.loop] + holes + [group.ins.areas['inner_loop'].loop])
1794 # Create and build air far field
1795 if run_type == 'EM':
1796 if self.data.magnet.geometry.electromagnetics.with_iron_yoke:
1797 for i in iron.key_points:
1798 gm.iron.max_radius = max(gm.iron.max_radius, max(iron.key_points[i].x, iron.key_points[i].y))
1799 greatest_radius = gm.iron.max_radius
1800 else: # no iron yoke data available
1801 for coil_nr, coil in self.geom.coil.coils.items():
1802 for pole_nr, pole in coil.poles.items():
1803 first_winding = list(pole.layers[len(pole.layers)].windings.keys())[0]
1804 first_block = list(pole.layers[len(pole.layers)].windings[first_winding].blocks)[0]
1805 gm.coil.max_radius = max(abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.x),
1806 abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.y),
1807 gm.coil.max_radius)
1808 greatest_radius = gm.coil.max_radius
1809 radius_in = greatest_radius * (2.5 if self.data.magnet.geometry.electromagnetics.with_iron_yoke else 6)
1810 radius_out = greatest_radius * (3.2 if self.data.magnet.geometry.electromagnetics.with_iron_yoke else 8)
1811 air_inf_center_x, air_inf_center_y = 0, 0
1812 for coil_nr, coil in self.md.geometries.coil.coils.items():
1813 air_inf_center_x += coil.bore_center.x
1814 air_inf_center_y += coil.bore_center.y
1815 gm.air.points['bore_center' + str(coil_nr)] = self.occ.addPoint(coil.bore_center.x, coil.bore_center.y, 0.)
1816 air_inf_center = [air_inf_center_x / len(self.md.geometries.coil.coils), air_inf_center_y / len(self.md.geometries.coil.coils)]
1817 if symmetry == 'none':
1818 gm.air_inf.lines['inner'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_in)
1819 gm.air_inf.lines['outer'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_out)
1820 gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['inner']]))
1821 gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['outer']]))
1822 gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop, gm.air_inf.areas['inner'].loop])
1823 else:
1824 pnt1 = [1, 0] if symmetry in ['xy', 'x'] else [0, -1]
1825 pnt2 = [0, 1] if symmetry in ['xy', 'y'] else [-1, 0]
1826 gm.air.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_in, pnt1[1] * radius_in, 0)
1827 gm.air.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_in, pnt2[1] * radius_in, 0)
1828 gm.air_inf.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_out, pnt1[1] * radius_out, 0)
1829 gm.air_inf.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_out, pnt2[1] * radius_out, 0)
1830 gm.air.lines['ln1'] = self.occ.addLine(gm.air.points['pnt1'], gm.air_inf.points['pnt1'])
1831 gm.air.lines['ln2'] = self.occ.addLine(gm.air.points['pnt2'], gm.air_inf.points['pnt2'])
1832 if not self.data.magnet.geometry.electromagnetics.with_iron_yoke:
1833 gm.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)
1834 gm.air_inf.lines['inner'] = self.occ.addCircleArc(gm.air.points['pnt2'], gm.air_inf.points['center'], gm.air.points['pnt1'])
1835 gm.air_inf.lines['outer'] = self.occ.addCircleArc(gm.air_inf.points['pnt2'], gm.air_inf.points['center'], gm.air_inf.points['pnt1'])
1837 if symmetry in ['xy', 'x']:
1838 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
1839 gm.iron.quadrants[1].points[self.symmetric_bnds['x_p']['pnts'][-1][0]], gm.air.points['pnt1'])
1840 self.symmetric_loop_lines['x'].append(gm.air.lines['x_p'])
1841 else: # y
1842 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'])
1843 self.symmetric_loop_lines['y'].append(gm.air.lines['y_n'])
1844 if symmetry in ['xy', 'y']:
1845 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'])
1846 self.symmetric_loop_lines['y'].insert(0, gm.air.lines['y_p'])
1847 else: # x
1848 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
1849 gm.iron.quadrants[2].points[self.symmetric_bnds['x_n']['pnts'][-1][0]], gm.air.points['pnt2'])
1850 self.symmetric_loop_lines['x'].insert(0, gm.air.lines['x_n'])
1852 inner_lines = self.symmetric_loop_lines['x'] + [gm.air_inf.lines['inner']] + self.symmetric_loop_lines['y']\
1853 if symmetry == 'xy' else self.symmetric_loop_lines[symmetry] + [gm.air_inf.lines['inner']]
1854 gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop(inner_lines))
1855 gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop(
1856 [gm.air.lines['ln1'], gm.air_inf.lines['outer'], gm.air.lines['ln2'], gm.air_inf.lines['inner']]))
1857 gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop])
1858 # self.md.domains.groups_entities.air_inf = [gm.air_inf.areas['outer'].surface]
1859 gm.air_inf.areas['inner'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['inner'].loop])
1861 # self.occ.synchronize()
1862 # self.gu.launch_interactive_GUI()
1864 def fragment(self):
1865 """
1866 Fragment and group air domains
1867 """
1868 # Collect surfaces to be subtracted by background air
1869 holes = []
1871 # Iron
1872 for group_name, surfaces in self.md.domains.groups_entities.iron.items():
1873 holes.extend([(2, s) for s in surfaces])
1874 # Coils
1875 for coil_nr, coil in self.md.geometries.coil.coils.items():
1876 for pole_nr, pole in coil.poles.items():
1877 for layer_nr, layer in pole.layers.items():
1878 for winding_nr, winding in layer.windings.items():
1879 for block_key, block in winding.blocks.items():
1880 for area_name, area in block.half_turns.areas.items():
1881 holes.append((2, area.surface))
1882 # Wedges
1883 for coil_nr, coil in self.md.geometries.wedges.coils.items():
1884 for layer_nr, layer in coil.layers.items():
1885 for wedge_nr, wedge in layer.wedges.items():
1886 for area_name, area in wedge.areas.items():
1887 holes.append((2, area.surface))
1888 # Insulation
1889 # if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
1890 # for coil_nr, coil in self.md.geometries.insulation.coils.items():
1891 # for group_nr, group in coil.group.items():
1892 # for area_name, area in group.ins.areas.items():
1893 # holes.append((2, area.surface))
1895 # Fragment
1896 fragmented = self.occ.fragment([(2, self.md.geometries.air_inf.areas['inner'].surface)], holes)[1]
1897 self.occ.synchronize()
1899 self.md.domains.groups_entities.air = []
1900 existing_domains = [e[0][1] for e in fragmented[1:]]
1901 for e in fragmented[0]:
1902 if e[1] not in existing_domains:
1903 self.md.domains.groups_entities.air.append(e[1])
1905 def updateTags(self, run_type, symmetry):
1907 # Update half turn line tags
1908 for coil_nr, coil in self.md.geometries.coil.coils.items():
1909 for pole_nr, pole in coil.poles.items():
1910 for layer_nr, layer in pole.layers.items():
1911 for winding_nr, winding in layer.windings.items():
1912 for block_key, block in winding.blocks.items():
1913 hts = block.half_turns
1914 # Get half turn ID numbers
1915 area_list = list(hts.areas.keys())
1916 for nr, ht_nr in enumerate(area_list):
1917 first_tag = int(min(gmsh.model.getAdjacencies(2, hts.areas[ht_nr].surface)[1]))
1918 hts.lines[ht_nr + 'i'] = first_tag
1919 hts.lines[ht_nr + 'l'] = first_tag + 1
1920 hts.lines[ht_nr + 'o'] = first_tag + 2
1921 hts.lines[ht_nr + 'h'] = first_tag + 3
1923 # Update insulation line tags
1924 if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:
1925 pass # todo
1927 # Update wedge line tags
1928 for coil_nr, coil in self.md.geometries.wedges.coils.items():
1929 for layer_nr, layer in coil.layers.items():
1930 for wedge_nr, wedge in layer.wedges.items():
1931 lines_tags = list(gmsh.model.getAdjacencies(2, wedge.areas[str(wedge_nr)].surface)[1])
1932 lines_tags.sort(key=lambda x: x)
1933 wedge.lines['i'] = int(lines_tags[0])
1934 wedge.lines['l'] = int(lines_tags[1])
1935 wedge.lines['o'] = int(lines_tags[2])
1936 wedge.lines['h'] = int(lines_tags[3])
1938 if run_type == 'EM':
1939 def _get_bnd_lines():
1940 return [pair[1] for pair in self.occ.getEntitiesInBoundingBox(corner_min[0], corner_min[1], corner_min[2],
1941 corner_max[0], corner_max[1], corner_max[2], dim=1)]
1943 tol = 1e-6
1944 # Update tags of air and air_inf arcs and their points
1945 lines_tags = gmsh.model.getAdjacencies(2, self.md.geometries.air_inf.areas['outer'].surface)[1]
1946 self.md.geometries.air_inf.lines['outer'] = int(lines_tags[0 if symmetry == 'none' else 1])
1947 self.md.geometries.air_inf.lines['inner'] = int(lines_tags[1 if symmetry == 'none' else 3])
1948 if symmetry == 'none': # todo: check if this holds for symmetric models too
1949 for coil_nr, coil in self.md.geometries.coil.coils.items():
1950 self.md.geometries.air.points['bore_center' + str(coil_nr)] += 2
1951 else:
1952 pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['outer'])[1])
1953 indexes = [0, 1, 0] if 'x' in symmetry else [1, 0, 1]
1954 pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] >\
1955 gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]
1956 self.md.geometries.air_inf.points['pnt1'] = int(pnt_tags[pnts[0]])
1957 self.md.geometries.air_inf.points['pnt2'] = int(pnt_tags[pnts[1]])
1958 pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['inner'])[1])
1959 pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] > \
1960 gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]
1961 self.md.geometries.air.points['pnt1'] = int(pnt_tags[pnts[0]])
1962 self.md.geometries.air.points['pnt2'] = int(pnt_tags[pnts[1]])
1963 for coil_nr, coil in self.md.geometries.coil.coils.items():
1964 self.md.geometries.air.points['bore_center' + str(coil_nr)] =(
1965 self.occ.getEntitiesInBoundingBox(-tol + coil.bore_center.x, -tol + coil.bore_center.y, -tol,
1966 tol + coil.bore_center.x, tol + coil.bore_center.y, tol, dim=0))[0][1]
1968 # Group symmetry boundary lines per type
1969 if symmetry == 'xy':
1970 corner_min = [-tol, -tol, -tol]
1971 corner_max = [gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0] + tol, tol, tol]
1972 self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()
1973 corner_max = [tol, gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1] + tol, tol]
1974 self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()
1975 elif symmetry == 'x':
1976 x_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0]
1977 corner_min = [- x_coord - tol, -tol, -tol]
1978 corner_max = [x_coord + tol, tol, tol]
1979 self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()
1980 elif symmetry == 'y':
1981 y_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1]
1982 corner_min = [-tol, - y_coord - tol, -tol]
1983 corner_max = [tol, y_coord + tol, tol]
1984 self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()
1986 # self.occ.synchronize()
1987 # self.gu.launch_interactive_GUI()