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