Coverage for fiqus/geom_generators/GeometryMultipole.py: 80%

1import copy

2import logging

3import os

4import gmsh

5import numpy as np

6import pandas as pd

7import json

9from fiqus.utils.Utils import GmshUtils

10from fiqus.utils.Utils import FilesAndFolders as Util

11from fiqus.utils.Utils import GeometricFunctions as Func

12from fiqus.data import DataFiQuS as dF

13from fiqus.data import DataMultipole as dM

14from fiqus.data.DataRoxieParser import FiQuSGeometry

15from fiqus.data.DataRoxieParser import Corner

16from fiqus.data.DataRoxieParser import Coord

17import re

19logger = logging.getLogger('FiQuS')

20class Geometry:

21 def __init__(self, data: dF.FDM() = None, geom: FiQuSGeometry() = None,

22 geom_folder: str = None, verbose: bool = False):

23 """

24 Class to generate geometry

25 :param data: FiQuS data model

26 :param geom: ROXIE geometry data

27 :param verbose: If True more information is printed in python console.

28 """

29 self.data: dF.FDM() = data

30 self.geom: FiQuSGeometry() = geom.Roxie_Data

32 # move cooling holes to a desired position

33 if self.data.magnet.solve.thermal.collar_cooling.move_cooling_holes:

34 self.geom.iron.key_points = self.move_keypoints(self.geom.iron.key_points, self.data.magnet.solve.thermal.collar_cooling.move_cooling_holes)

36 self.geom_folder = geom_folder

37 self.verbose: bool = verbose

39 self.md = dM.MultipoleData()

41 self.gu = GmshUtils(self.geom_folder, self.verbose)

42 self.gu.initialize(verbosity_Gmsh=self.data.run.verbosity_Gmsh)

43 self.occ = gmsh.model.occ

45 self.model_file = os.path.join(self.geom_folder, self.data.general.magnet_name)

47 self.blk_ins_lines = {} # for meshed insulation

48 self.ins_wire_lines = {} # for meshed insulation

49 self.block_coil_mid_pole_blks = {}

51 self.nc = {'collar': 'c', 'iron_yoke': 'i', 'poles': 'p'}

52 self.inv_nc = {v: k for k, v in self.nc.items()} #invert naming convention

54 if self.data.magnet.geometry.electromagnetics.symmetry != 'none':

55 self.symmetric_loop_lines = {'x': [], 'y': []}

56 self.symmetric_bnds = {'x_p': {'pnts': [], 'line_pnts': []}, 'y_p': {'pnts': [], 'line_pnts': []},

57 'x_n': {'pnts': [], 'line_pnts': []}, 'y_n': {'pnts': [], 'line_pnts': []}}

59 def clear(self):

60 self.md = dM.MultipoleData()

61 self.block_coil_mid_pole_blks = {}

62 gmsh.clear()

64 def ending_step(self, gui: bool = False):

65 if gui:

66 self.gu.launch_interactive_GUI()

67 else:

68 if gmsh.isInitialized():

69 gmsh.clear()

70 gmsh.finalize()

72 def saveHalfTurnCornerPositions(self):

73 self.occ.synchronize()

74 iH, iL, oH, oL, iHr, iLr, oHr, oLr = [], [], [], [], [], [], [], []

75 for po in self.geom.coil.physical_order:

76 block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[

77 po.winding].blocks[po.block]

78 for halfTurn_nr, halfTurn in block.half_turns.items():

79 ht = halfTurn.corners.insulated

80 ht_b = halfTurn.corners.bare

81 iHr.append([ht_b.iH.x, ht_b.iH.y])

82 iLr.append([ht_b.iL.x, ht_b.iL.y])

83 oHr.append([ht_b.oH.x, ht_b.oH.y])

84 oLr.append([ht_b.oL.x, ht_b.oL.y])

85 iH.append([ht.iH.x, ht.iH.y])

86 iL.append([ht.iL.x, ht.iL.y])

87 oH.append([ht.oH.x, ht.oH.y])

88 oL.append([ht.oL.x, ht.oL.y])

89 with open(f"{self.model_file}.crns", 'w') as f:

90 json.dump({'iH': iH, 'iL': iL, 'oH': oH, 'oL': oL,

91 'iHr': iHr, 'iLr': iLr, 'oHr': oHr, 'oLr': oLr}, f)

93 def saveStrandPositions(self, run_type):

94 symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'

95 ht_nr = 0

96 std_nr = 0

97 parser_x, parser_y, blocks, ht, std, pole_blocks = [], [], [], [], [], []

98 for po in self.geom.coil.physical_order:

99 block = self.geom.coil.coils[po.coil].poles[po.pole].layers[po.layer].windings[

100 po.winding].blocks[po.block]

101 if po.pole == 1: pole_blocks.append(po.block)

102 for halfTurn_nr, halfTurn in block.half_turns.items():

103 ht_nr += 1

104 for strand_group_nr, strand_group in halfTurn.strand_groups.items():

105 for strand_nr, strand in strand_group.strand_positions.items():

106 std_nr += 1

107 blocks.append(po.block)

108 ht.append(ht_nr)

109 std.append(std_nr)

110 parser_x.append(strand.x)

111 parser_y.append(strand.y)

112 mirrored = {}

113 condition = {2: [1, -1], 3: [1, 1], 4: [-1, 1]}

114 if symmetry == 'xy': mirroring = {2: [-1, 1], 3: [-1, -1], 4: [1, -1]}

115 elif symmetry == 'x': mirroring = {3: [1, -1], 4: [1, -1]}

116 elif symmetry == 'y': mirroring = {2: [-1, 1], 3: [-1, 1]}

117 else: mirroring = {}

118 if mirroring:

119 df = pd.DataFrame({'parser_x': parser_x, 'parser_y': parser_y}, index=std)

120 for qdr, mrr in mirroring.items():

121 subdf = df[(condition[qdr][0] * df['parser_x'] < 0) & (condition[qdr][1] * df['parser_y'] < 0)]

122 for strand, x, y in zip(subdf.index, subdf['parser_x'], subdf['parser_y']):

123 mirrored[strand] = df[(df['parser_x'] == mrr[0] * x) & (df['parser_y'] == mrr[1] * y)].index.item()

124 with open(f"{self.model_file}_{run_type}.strs", 'w') as f:

125 json.dump({'x': parser_x, 'y': parser_y, 'block': blocks, 'ht': ht, 'mirrored': mirrored,

126 'pole_1_blocks': pole_blocks, 'poles': len(self.geom.coil.coils[1].poles)},

127 f)

128

129 def saveBoundaryRepresentationFile(self, run_type):

130 self.occ.synchronize()

131 gmsh.write(f'{self.model_file}_{run_type}.brep')

132 gmsh.clear()

133

134 def loadBoundaryRepresentationFile(self, run_type):

135 gmsh.option.setString('Geometry.OCCTargetUnit', 'M') # set units to meters

136 gmsh.open(os.path.join(f'{self.model_file}_{run_type}.brep'))

137

138 def saveAuxiliaryFile(self, run_type):

139 Util.write_data_to_yaml(f'{self.model_file}_{run_type}.aux', self.md.model_dump())

140

141 @staticmethod

142 def findMidLayerPoint(bc_current, bc_next, center, mean_rad):

143 mid_layer = [(bc_current.x + bc_next.x) / 2, (bc_current.y + bc_next.y) / 2]

144 mid_rad = Func.points_distance(mid_layer, [center.x, center.y])

145 dist_from_mid = mean_rad - mid_rad

146 angle = Func.arc_angle_between_point_and_abscissa(mid_layer, [center.x, center.y])

147 mid_layer[0] += dist_from_mid * np.cos(angle)

148 mid_layer[1] += dist_from_mid * np.sin(angle)

149 return mid_layer

150

151 @staticmethod

152 def getMidLayerEndpoints(el_current, el_next, center, mid_layer_arc_pnt=None, coil_type='cos-theta', cable_type='Rutherford', is_for_mid_pole=False):

153 thin_shell_endpoints = {'higher': list, 'lower': list}

154 which_block = {'higher': str, 'lower': str}

155 angles = {'higher': float, 'lower': float}

156 # Check if the element crosses the x axis

157 angles_to_correct = []

158 correction_angle = 0

159 l_curr = Func.arc_angle_between_point_and_abscissa([el_current.iL.x, el_current.iL.y], center)

160 h_curr = Func.arc_angle_between_point_and_abscissa([el_current.iH.x, el_current.iH.y], center)

161 l_next = Func.arc_angle_between_point_and_abscissa([el_next.iL.x, el_next.iL.y], center)

162 h_next = Func.arc_angle_between_point_and_abscissa([el_next.iH.x, el_next.iH.y], center)

163 if abs(l_curr - h_curr) > np.pi:

164 angles_to_correct.append('current')

165 correction_angle = max(1.05 * (2 * np.pi - l_curr), correction_angle)

166 if abs(l_next - h_next) > np.pi:

167 angles_to_correct.append('next')

168 correction_angle = max(1.05 * (2 * np.pi - l_next), correction_angle)

169 for side in thin_shell_endpoints.keys():

170 if mid_layer_arc_pnt:

171 if side == 'higher':

172 mid_lyr_curr, mid_lyr_next = [el_current.oH, el_current.iH], [el_next.oH, el_next.iH]

173 else:

174 mid_lyr_curr, mid_lyr_next = [el_current.oL, el_current.iL], [el_next.oL, el_next.iL]

175 if cable_type in ['Mono', 'Ribbon']:

176 pnts_curr = Func.intersection_between_circle_and_line(

177 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]),

178 [center, mid_layer_arc_pnt])

179 pnt_curr = pnts_curr[0] if Func.points_distance(pnts_curr[0], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) <\

180 Func.points_distance(pnts_curr[1], [mid_lyr_curr[0].x, mid_lyr_curr[0].y]) else pnts_curr[1]

181 pnts_next = Func.intersection_between_circle_and_line(

182 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]),

183 [center, mid_layer_arc_pnt])

184 pnt_next = pnts_next[0] if Func.points_distance(pnts_next[0], [mid_lyr_next[0].x, mid_lyr_next[0].y]) <\

185 Func.points_distance(pnts_next[1], [mid_lyr_next[0].x, mid_lyr_next[0].y]) else pnts_next[1]

186 elif cable_type == 'Rutherford':

187 pnt_curr = Func.intersection_between_circle_and_line(

188 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]),

189 [center, mid_layer_arc_pnt], get_only_closest=True)[0]

190 pnt_next = Func.intersection_between_circle_and_line(

191 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]),

192 [center, mid_layer_arc_pnt], get_only_closest=True)[0]

193 else:

194 if cable_type == 'Rutherford':

195 if coil_type == 'common-block-coil':

196 mid_layer_x = (el_current.oH.x + el_next.iH.x) / 2

197 if side == 'higher':

198 pnt_curr, pnt_next = [mid_layer_x, el_current.iH.y], [mid_layer_x, el_next.iH.y]

199 else:

200 pnt_curr, pnt_next = [mid_layer_x, el_current.iL.y], [mid_layer_x, el_next.iL.y]

201 else:

202 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

203 if side == 'higher':

204 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]

205 else:

206 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]

207 elif cable_type in ['Mono', 'Ribbon']:

208 pnt_curr = [(el_current.oH.x + el_next.iH.x) / 2, (el_current.oH.y + el_next.iH.y) / 2] if side == 'higher'\

209 else [(el_current.oL.x + el_next.iL.x) / 2, (el_current.oL.y + el_next.iL.y) / 2]

210 pnt_next = pnt_curr

211 angle_curr = Func.arc_angle_between_point_and_abscissa(pnt_curr, center)

212 angle_next = Func.arc_angle_between_point_and_abscissa(pnt_next, center)

213 if 'current' in angles_to_correct:

214 angle_curr = angle_curr + correction_angle - (2 * np.pi if side == 'lower' else 0)

215 elif 'next' in angles_to_correct:

216 if angle_curr < np.pi / 2: angle_curr += correction_angle

217 elif angle_curr > np.pi * 3 / 2: angle_curr = angle_curr + correction_angle - 2 * np.pi

218 if 'next' in angles_to_correct:

219 angle_next = angle_next + correction_angle - (2 * np.pi if side == 'lower' else 0)

220 elif 'current' in angles_to_correct:

221 if angle_next < np.pi / 2: angle_next += correction_angle

222 elif angle_next > np.pi * 3 / 2: angle_next = angle_next + correction_angle - 2 * np.pi

223 if abs(angle_curr - angle_next) < 1e-6:

224 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'

225 elif angle_curr * (-1 if side == 'lower' else 1) < angle_next * (-1 if side == 'lower' else 1):

226 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_curr, angle_curr, 'current'

227 else:

228 thin_shell_endpoints[side], angles[side], which_block[side] = pnt_next, angle_next, 'next'

229 if angles['higher'] < angles['lower']: return None

230 else: return thin_shell_endpoints, which_block

231

232 def create_geom_dict(self, geometry_setting):

233 return {v: k in geometry_setting.areas for k, v in self.nc.items()}

234

235 def constructIronGeometry(self, symmetry, geometry_setting, run_type):

236 """

237 Generates points, hyper lines, and curve loops for the iron yoke

238 """

239 iron = self.geom.iron #roxie

240

241 if symmetry == 'xy':

242 self.md.geometries.iron.quadrants = {1: dM.Region()}

243 list_bnds = ['x_p', 'y_p']

244 elif symmetry == 'x':

245 self.md.geometries.iron.quadrants = {1: dM.Region(), 2: dM.Region()}

246 list_bnds = ['x_p', 'x_n']

247 elif symmetry == 'y':

248 self.md.geometries.iron.quadrants = {1: dM.Region(), 4: dM.Region()}

249 list_bnds = ['y_p', 'y_n']

250 else:

251 for k in self.nc.keys(): getattr(self.md.geometries, k).quadrants = {1: dM.Region(), 2: dM.Region(), 4: dM.Region(), 3: dM.Region()}

252 list_bnds = []

253

254 lc = 1e-2

255 geom_dict = self.create_geom_dict(geometry_setting)

256

257 for point_name, point in iron.key_points.items():

258 identifier = next((k for k in self.inv_nc.keys() if re.match(f'^{k}', point_name[2:])), None)

259 if not geom_dict.get(identifier, False): continue

260 quadrants = getattr(self.md.geometries, self.inv_nc[identifier]).quadrants #re.sub(r'\d+', '', point_name[2:])

261 if symmetry in ['x', 'xy']:

262 if point.y == 0.:

263 self.symmetric_bnds['x_p']['pnts'].append([point_name, point.x])

264 if symmetry in ['y', 'xy']:

265 if point.x == 0.:

266 self.symmetric_bnds['y_p']['pnts'].append([point_name, point.y])

267 quadrants[1].points[point_name] = self.occ.addPoint(point.x, point.y, 0, lc)

268 if symmetry in ['x', 'none']:

269 if point.x == 0.:

270 quadrants[2].points[point_name] = quadrants[1].points[point_name]

271 else:

272 quadrants[2].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]

273 self.occ.mirror([(0, quadrants[2].points[point_name])], 1, 0, 0, 0)

274 if point.y == 0. and symmetry == 'x':

275 self.symmetric_bnds['x_n']['pnts'].append([point_name, point.x])

276 if symmetry in ['y', 'none']:

277 if point.y == 0.:

278 quadrants[4].points[point_name] = quadrants[1].points[point_name]

279 else:

280 quadrants[4].points[point_name] = self.occ.copy([(0, quadrants[1].points[point_name])])[0][1]

281 self.occ.mirror([(0, quadrants[4].points[point_name])], 0, 1, 0, 0)

282 if point.x == 0. and symmetry == 'y':

283 self.symmetric_bnds['y_n']['pnts'].append([point_name, point.y])

284 if symmetry == 'none':

285 if point.y == 0.:

286 quadrants[3].points[point_name] = quadrants[2].points[point_name]

287 elif point.x == 0.:

288 quadrants[3].points[point_name] = quadrants[4].points[point_name]

289 else:

290 quadrants[3].points[point_name] = self.occ.copy([(0, quadrants[2].points[point_name])])[0][1]

291 self.occ.mirror([(0, quadrants[3].points[point_name])], 0, 1, 0, 0)

292

293 mirror_x = [1, -1, -1, 1]

294 mirror_y = [1, 1, -1, -1]

295 symmetric_bnds_order = {'x': [], 'y': []}

296 sym_lines_tags = {'x_p': [], 'y_p': [], 'x_n': [], 'y_n': []}

297 for line_name, line in iron.hyper_lines.items():

298 identifier = next((k for k in self.inv_nc.keys() if re.match(f'^{k}', line_name[2:])), None)

299 if not geom_dict.get(identifier, False): continue

300 quadrants = getattr(self.md.geometries, self.inv_nc[identifier]).quadrants #re.sub(r'\d+', '', line_name[2:])

301 pt1 = iron.key_points[line.kp1]

302 pt2 = iron.key_points[line.kp2]

303 if line.type == 'line':

304 for quadrant, qq in quadrants.items():

305 if quadrant == 1:

306 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])

307 if pt1.y == 0. and pt2.y == 0. and 'x_p' in list_bnds:

308 self.symmetric_bnds['x_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)

309 sym_lines_tags['x_p'].append(qq.lines[line_name])

310 symmetric_bnds_order['x'].append(min(pt1.x, pt2.x))

311 elif pt1.x == 0. and pt2.x == 0. and 'y_p' in list_bnds:

312 self.symmetric_bnds['y_p']['line_pnts'].append(line.kp1 + '_' + line.kp2)

313 sym_lines_tags['y_p'].append(qq.lines[line_name])

314 symmetric_bnds_order['y'].append(min(pt1.y, pt2.y))

315 elif quadrant == 2:

316 if pt1.x == 0. and pt2.x == 0.:

317 qq.lines[line_name] = quadrants[1].lines[line_name]

318 else:

319 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])

320 if pt1.y == 0. and pt2.y == 0. and 'x_n' in list_bnds:

321 self.symmetric_bnds['x_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)

322 sym_lines_tags['x_n'].append(qq.lines[line_name])

323 elif quadrant == 4:

324 if pt1.y == 0. and pt2.y == 0.:

325 qq.lines[line_name] = quadrants[1].lines[line_name]

326 else:

327 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])

328 if pt1.x == 0. and pt2.x == 0. and 'y_n' in list_bnds:

329 self.symmetric_bnds['y_n']['line_pnts'].append(line.kp1 + '_' + line.kp2)

330 sym_lines_tags['y_n'].append(qq.lines[line_name])

331 else: # 3

332 if pt1.y == 0. and pt2.y == 0.:

333 qq.lines[line_name] = quadrants[2].lines[line_name]

334 elif pt1.x == 0. and pt2.x == 0.:

335 qq.lines[line_name] = quadrants[4].lines[line_name]

336 else:

337 qq.lines[line_name] = self.occ.addLine(qq.points[line.kp1], qq.points[line.kp2])

338

339 elif line.type == 'arc':

340 center = Func.arc_center_from_3_points([pt1.x, pt1.y],

341 [iron.key_points[line.kp3].x, iron.key_points[line.kp3].y],

342 [pt2.x, pt2.y])

343 new_point_name = 'kp' + line_name + '_center'

344 arc_coordinates1 = (pt1.x, pt1.y)

345 arc_coordinates2 = (pt2.x, pt2.y)

346 arc_coordinates3 = (iron.key_points[line.kp3].x, iron.key_points[line.kp3].y)

347

348 # This code addresses a meshing error in MQXA and MB_2COILS that occurs when an arc is defined on any of

349 # the axes. The issue arises because the function Func.arc_center_from_3_points does not return exactly

350 # zero but a value with a magnitude of approximately 10^-17 when the two points are placed on the axes.

351 # Consequently, when using the method self.occ.addCircleArc(), which only takes in three points without

352 # specifying a direction, a problem arises. The addCircleArc() function always creates the arc with the

353 # smallest angle. However, since center point can be slightly above or below the axis, the arc can

354 # inadvertently be drawn in the wrong quadrant, leading to an incorrect result.

355 # -----------------------

356 # Check that arcs with points on the x-axis are drawn in the first quadrant

357 if arc_coordinates3[1] > 0 and arc_coordinates2[1] == 0 and arc_coordinates1[1] == 0 and center[1] > 0:

358 quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], -center[1], 0)

359 # Check that arcs with points on the y-axis are drawn in the first quadrant

360 elif arc_coordinates3[0] > 0 and arc_coordinates2[0] == 0 and arc_coordinates1[0] == 0 and center[0] > 0:

361 quadrants[1].points[new_point_name] = self.occ.addPoint(-center[0], center[1], 0)

362 else:

363 quadrants[1].points[new_point_name] = self.occ.addPoint(center[0], center[1], 0)

364 # -----------------------

365 # gmsh.model.setEntityName(0, gm.iron.quadrants[1].points[new_point_name], 'iron_' + new_point_name)

366 if symmetry in ['x', 'none']:

367 if center[0] == 0.:

368 quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]

369 else:

370 quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]

371 self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)

372 if symmetry in ['y', 'none']:

373 if center[1] == 0.:

374 quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]

375 else:

376 quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]

377 self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)

378 if symmetry == 'none':

379 if center[1] == 0.:

380 quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]

381 else:

382 quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]

383 self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)

384

385 for quadrant, qq in quadrants.items():

386 qq.lines[line_name] = self.occ.addCircleArc(

387 qq.points[line.kp1], qq.points[line.kp3], qq.points[line.kp2], center=False)

388

389 elif line.type == 'circle':

390 center = [(pt1.x + pt2.x) / 2, (pt1.y + pt2.y) / 2]

391 radius = (np.sqrt(np.square(pt1.x - center[0]) + np.square(pt1.y - center[1])) +

392 np.sqrt(np.square(pt2.x - center[0]) + np.square(pt2.y - center[1]))) / 2

393

394 for quadrant, qq in quadrants.items():

395 qq.lines[line_name] = self.occ.addCircle(

396 mirror_x[quadrant - 1] * center[0], mirror_y[quadrant - 1] * center[1], 0, radius)

397 qq.points['kp' + line_name] = len(qq.points) + 1

398

399 elif line.type == 'ellipticArc':

400 a, b = line.arg1, line.arg2

401 x1, y1 = pt1.x, pt1.y

402 x2, y2 = pt2.x, pt2.y

403 x3 = np.power(x1, 2.0)

404 y3 = np.power(y1, 2.0)

405 x4 = np.power(x2, 2.0)

406 y4 = np.power(y2, 2.0)

407 a2 = np.power(a, 2.0)

408 b2 = np.power(b, 2.0)

409 expression = -4.0 * a2 * b2 + a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 - 2.0 * b2 * x1 * x2 + b2 * x4

410 xc = x1 / 2.0 + x2 / 2.0 - a * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3 -

411 2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (y1 - y2) / (2.0 * b)

412 yc = y1 / 2.0 + y2 / 2.0 + b * np.power(- expression / (a2 * y3 - 2.0 * a2 * y1 * y2 + a2 * y4 + b2 * x3

413 - 2.0 * b2 * x1 * x2 + b2 * x4), 0.5) * (x1 - x2) / (2.0 * a)

414

415 center = self.occ.addPoint(xc, yc, 0, lc)

416 axis_point_a = self.occ.addPoint(xc + a, yc, 0, lc)

417 axis_point_b = self.occ.addPoint(xc, yc + b, 0, lc)

418

419 new_point_name = 'kp' + line_name + '_center'

420 new_axis_a_point_name = 'kp' + line_name + '_a'

421 new_axis_b_point_name = 'kp' + line_name + '_b'

422

423 quadrants[1].points[new_point_name] = center

424 quadrants[1].points[new_axis_a_point_name] = axis_point_a

425 quadrants[1].points[new_axis_b_point_name] = axis_point_b

426

427 if symmetry in ['x', 'none']:

428 if xc == 0.: # Least amount of possible points.

429 quadrants[2].points[new_point_name] = quadrants[1].points[new_point_name]

430 quadrants[2].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]

431 quadrants[2].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]

432 else:

433 quadrants[2].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]

434 quadrants[2].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]

435 quadrants[2].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]

436 self.occ.mirror([(0, quadrants[2].points[new_point_name])], 1, 0, 0, 0)

437 self.occ.mirror([(0, quadrants[2].points[new_axis_a_point_name])], 1, 0, 0, 0)

438 self.occ.mirror([(0, quadrants[2].points[new_axis_b_point_name])], 1, 0, 0, 0)

439 if symmetry in ['y', 'none']:

440 if yc == 0.:

441 quadrants[4].points[new_point_name] = quadrants[1].points[new_point_name]

442 quadrants[4].points[new_axis_a_point_name] = quadrants[1].points[new_axis_a_point_name]

443 quadrants[4].points[new_axis_b_point_name] = quadrants[1].points[new_axis_b_point_name]

444 else:

445 quadrants[4].points[new_point_name] = self.occ.copy([(0, quadrants[1].points[new_point_name])])[0][1]

446 self.occ.mirror([(0, quadrants[4].points[new_point_name])], 0, 1, 0, 0)

447 quadrants[4].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_a_point_name])])[0][1]

448 self.occ.mirror([(0, quadrants[4].points[new_axis_a_point_name])], 0, 1, 0, 0)

449 quadrants[4].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[1].points[new_axis_b_point_name])])[0][1]

450 self.occ.mirror([(0, quadrants[4].points[new_axis_b_point_name])], 0, 1, 0, 0)

451 if symmetry == 'none':

452 if yc == 0.:

453 quadrants[3].points[new_point_name] = quadrants[2].points[new_point_name]

454 quadrants[3].points[new_axis_a_point_name] = quadrants[2].points[new_axis_a_point_name]

455 quadrants[3].points[new_axis_b_point_name] = quadrants[2].points[new_axis_b_point_name]

456 else:

457 quadrants[3].points[new_point_name] = self.occ.copy([(0, quadrants[2].points[new_point_name])])[0][1]

458 self.occ.mirror([(0, quadrants[3].points[new_point_name])], 0, 1, 0, 0)

459 quadrants[3].points[new_axis_a_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_a_point_name])])[0][1]

460 self.occ.mirror([(0, quadrants[3].points[new_axis_a_point_name])], 0, 1, 0, 0)

461 quadrants[3].points[new_axis_b_point_name] = self.occ.copy([(0, quadrants[2].points[new_axis_b_point_name])])[0][1]

462 self.occ.mirror([(0, quadrants[3].points[new_axis_b_point_name])], 0, 1, 0, 0)

463

464 for quadrant, qq in quadrants.items():

465 qq.lines[line_name] = self.occ.addEllipseArc(

466 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],

467 qq.points[line.kp2])

468

469 else:

470 raise ValueError('Hyper line {} not supported'.format(line.type))

471

472 if symmetry != 'none':

473 quadrants = self.md.geometries.iron_yoke.quadrants

474 indexes = {'x_p': 1, 'y_p': 1, 'x_n': 1, 'y_n': 1}

475 self.md.geometries.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)

476 for sym in list_bnds:

477 if sym in ['x_p', 'y_p']:

478 quadrant = 1

479 elif sym == 'x_n':

480 quadrant = 2

481 else: # 'y_n'

482 quadrant = 4

483 sym_lines_tags[sym] = [x for _, x in sorted(zip(symmetric_bnds_order[sym[0]], sym_lines_tags[sym]))]

484

485 self.symmetric_bnds[sym]['pnts'].append(['center', 0])

486 self.symmetric_bnds[sym]['pnts'].sort(key=lambda x: x[1])

487 self.md.geometries.symmetric_boundaries.lines[sym + '_center'] = self.occ.addLine(

488 self.md.geometries.air_inf.points['center'], quadrants[quadrant].points[self.symmetric_bnds[sym]['pnts'][1][0]])

489 sym_lines_tags[sym].insert(0, self.md.geometries.symmetric_boundaries.lines[sym + '_center'])

490 for i, pnt in enumerate(self.symmetric_bnds[sym]['pnts'][1:-1]):

491 pnt_next = self.symmetric_bnds[sym]['pnts'][i + 2][0]

492 if not any(pnt[0] in s and pnt_next in s for s in self.symmetric_bnds[sym]['line_pnts']):

493 self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]] =\

494 self.occ.addLine(quadrants[quadrant].points[pnt[0]], quadrants[quadrant].points[pnt_next])

495 sym_lines_tags[sym].insert(indexes[sym], self.md.geometries.symmetric_boundaries.lines[sym + '_' + pnt[0]])

496 indexes[sym] += 1

497 if symmetry == 'xy':

498 self.symmetric_loop_lines['x'] = sym_lines_tags['x_p']

499 sym_lines_tags['y_p'].reverse()

500 self.symmetric_loop_lines['y'] = sym_lines_tags['y_p']

501 elif symmetry == 'x':

502 sym_lines_tags['x_n'].reverse()

503 self.symmetric_loop_lines['x'] = sym_lines_tags['x_n'] + sym_lines_tags['x_p']

504 elif symmetry == 'y':

505 sym_lines_tags['y_p'].reverse()

506 self.symmetric_loop_lines['y'] = sym_lines_tags['y_p'] + sym_lines_tags['y_n']

507

508 # add all areas of each quadrant. Useful for brep curves and meshing

509 for key in geometry_setting.areas: # only consider areas that are implemented

510 quadrants = getattr(self.md.geometries, key).quadrants

511 for quadrant, qq in quadrants.items():

512 for area_name, area in iron.hyper_areas.items(): ## all areas

513 def _add_loop():

514 # prevent additional curveloop generation when Enforcing the TSA mapping on the collar

515 if (run_type == 'TH'

516 and self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping

517 and (area_name.startswith('arc') and not area_name.startswith('arch') or area_name.startswith('arp'))

518 ): # need to disable the pole area too as it is linked to the same curve

519 qq.areas[area_name] = dM.Area() ## initialise area without loop

520 else:

521 qq.areas[area_name] = dM.Area(

522 loop=self.occ.addCurveLoop([qq.lines[line] for line in area.lines]))

523

524 if iron.hyper_areas[area_name].material not in getattr(self.md.domains.groups_entities, key) and \

525 iron.hyper_areas[area_name].material != 'BH_air': ## add the material to the keys

526 # for the collar region, it is possible to overwrite the material -> intercept it here

527 if key == 'collar' and (self.data.magnet.solve.collar.material != iron.hyper_areas[area_name].material) and self.data.magnet.solve.collar.material is not None:

528 logger.warning("Overwriting the collar material for area {} to {} ".format(area_name, self.data.magnet.solve.collar.material))

529 iron.hyper_areas[area_name].material = self.data.magnet.solve.collar.material

530 getattr(self.md.domains.groups_entities, key)[iron.hyper_areas[area_name].material] = []

531

532 identifier = next((k for k in geom_dict.keys() if re.match(f'^{k}', area_name[2:])),

533 None) # match key from geom_dict to the area name (see naming convention)

534

535 if key == self.inv_nc.get(identifier, None): # re.sub(r'\d+', '', area_name[2:]),

536 _add_loop() # adds arch to collar, because c is in the naming convention of the collar

537 elif area_name.startswith('arh') and key == 'iron_yoke': # if not previous but it is a hole, assume iron

538 _add_loop()

539

540 # define inner collar lines

541 def define_inner_collar():

542 """

543 Defines the inner collar line used for the thermal TSA + for the A projection

544 """

545 self.occ.synchronize()

546 # only works if the inner collar line is an arc -> just disable 'arc' and calc for all lines

547 # alternative method. Find all "arc" lines and then select the closest to the center

548 for quad, object in self.md.geometries.collar.quadrants.items():

549 arc_line_tags = [object.lines[name] for name in object.lines.keys() if

550 self.geom.iron.hyper_lines[name].type == 'arc']

551 closest_dist = 1000.

552 for tag in arc_line_tags:

553 x, y, _ = gmsh.model.getValue(1, tag, [0.5]) # pick one point on the arc

554 dist = np.sqrt(x ** 2 + y ** 2)

555 if dist < closest_dist:

556 closest_dist = dist

557 closest_line = tag ## assumes it is only one line per quadrant

558 self.md.geometries.collar.inner_boundary_tags[quad] = [closest_line]

559 def define_collar_cooling():

560 """

561 Defines the cooling holes in the collar

562 """

563 self.occ.synchronize()

564 line_names = [item for key in self.geom.iron.hyper_areas.keys() if 'ch' in key for item in self.geom.iron.hyper_areas[key].lines]

565 # line names are the same in each quadrant. Tags are unique

566 for quad, qq in self.md.geometries.collar.quadrants.items():

567 self.md.geometries.collar.cooling_tags.extend([qq.lines[line] for line in line_names])

568 # these tags are only used to be skipped for enfrocing_TSA_mapping

569

570 # we need the inner collar lines if we want to do TSA, so no need to define it

571 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA_new: define_inner_collar()

572 # we only need to specify the air holes if we want cooling OR if we enforce TSA nodes on the collar

573 if run_type == 'TH' and (self.data.magnet.solve.thermal.collar_cooling.enabled or self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping): define_collar_cooling()

574

575

576

577 def constructWedgeGeometry(self, use_TSA):

578 """

579 Generates points, hyper lines, and curve loops for the wedges

580 """

581 def _addMidLayerThinShellPoints(wedge_current):

582 def __addThinShellPoints(side_case, mid_layer_ts):

583 if side_case == 'outer':

584 mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +

585 Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2

586 else:

587 mean_rad_current = (Func.points_distance([wedge_current.iH.x, wedge_current.iH.y], wedge_center) +

588 Func.points_distance([wedge_current.iL.x, wedge_current.iL.y], wedge_center)) / 2

589 are_endpoints = {}

590 for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():

591 blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]

592 blk_next = wnd.blocks[blk_nr_next]

593 ht_list_next = (list(blk_next.half_turns.keys()) if blk_nr_next == list(wnd.blocks.keys())[0] else list(

594 reversed(blk_next.half_turns.keys())))

595 hh = blk_next.half_turns[ht_list_next[-1]].corners.bare

596 ll = blk_next.half_turns[ht_list_next[0]].corners.bare

597 bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)

598 if side_case == 'outer':

599 block_list = self.md.geometries.coil.anticlockwise_order.coils[wedge.order_l.coil].layers[wedge.order_l.layer + 1]

600 blk_index = [blk.block for blk in block_list].index(blk_nr_next)

601 if blk_index + 1 == len(block_list): blk_index = -1

602 for blk in block_list[blk_index + 1:] + block_list[:blk_index + 1]:

603 if blk.winding == block_list[blk_index].winding:

604 ht_index = -1

605 break

606 elif blk.pole != block_list[blk_index].pole:

607 ht_index = 0

608 break

609 hh = blk_next.half_turns[ht_list_next[ht_index]].corners.bare

610 ll = blk_next.half_turns[ht_list_next[0 if ht_index == -1 else -1]].corners.bare

611 mean_rad_next = (Func.points_distance([hh.iH.x, hh.iH.y], wedge_center) +

612 Func.points_distance([ll.iL.x, ll.iL.y], wedge_center)) / 2

613 else:

614 mean_rad_next = (Func.points_distance([bc_next.oH.x, bc_next.oH.y], wedge_center) +

615 Func.points_distance([bc_next.oL.x, bc_next.oL.y], wedge_center)) / 2

616 mean_rad = (mean_rad_current + mean_rad_next) / 2

617 mid_layer = self.findMidLayerPoint(wedge_current.oH, bc_next.iH, wedge.corrected_center.outer, mean_rad)\

618 if side_case == 'outer' else self.findMidLayerPoint(wedge_current.iH, bc_next.oH, wedge.corrected_center.inner, mean_rad)

619 are_endpoints[wnd_nr] = self.getMidLayerEndpoints(wedge_current, bc_next, wedge_center, mid_layer_arc_pnt=mid_layer)

620 for wnd_nr, wnd in pole.layers[wedge.order_l.layer + (1 if side_case == 'outer' else -1)].windings.items():

621 blk_nr_next = list(wnd.blocks.keys())[blk_list_current.index(wedge.order_l.block)]

622 blk_next = wnd.blocks[blk_nr_next]

623 is_first_blk_next = blk_nr_next == list(wnd.blocks.keys())[0]

624 ht_list_next = (list(blk_next.half_turns.keys()) if is_first_blk_next else list(

625 reversed(blk_next.half_turns.keys())))

626 if are_endpoints[wnd_nr]: # this is empty if the wedge and the block are not radially adjacent

627 endpoints = are_endpoints[wnd_nr][0]

628 which_entity = are_endpoints[wnd_nr][1]

629 mid_layer_name = 'w' + str(wedge_nr) + '_' + str(blk_nr_next)

630 mid_layer_ts[mid_layer_name] = dM.Region()

631 ts_wdg = mid_layer_ts[mid_layer_name]

632 beg = ('w' + str(wedge_nr) if which_entity['lower'] == 'current' else str(ht_list_next[0])) + 'l'

633 ts_wdg.points[beg] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)

634 ht_lower_angles = {}

635 for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):

636 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]],

637 [[ht.corners.bare.oL.x, ht.corners.bare.oL.y], [ht.corners.bare.oH.x, ht.corners.bare.oH.y]],

638 ['l', 'h']):

639 line_pars_current = Func.line_through_two_points(pnt1, pnt2)

640 intersect_prev = Func.intersection_between_arc_and_line(

641 line_pars_current, [wedge_center, endpoints['higher'], endpoints['lower']])

642 if intersect_prev:

643 ts_wdg.points[str(ht_nr) + side] = self.occ.addPoint(intersect_prev[0][0], intersect_prev[0][1], 0)

644 elif side == 'l':

645 intrsc = Func.intersection_between_circle_and_line(line_pars_current, [wedge_center, endpoints['lower']], get_only_closest=True)[0]

646 ht_lower_angles[ht_nr] = Func.arc_angle_between_point_and_abscissa([intrsc[0], intrsc[1]], wedge_center)

647 end = ('w' + str(wedge_nr) if which_entity['higher'] == 'current' else str(ht_list_next[-1])) + 'h'

648 if all('w' in pnt_name for pnt_name in list(ts_wdg.points.keys())): # only one thin-shell 'within' the facing half-turn

649 wdg_angle_il = Func.arc_angle_between_point_and_abscissa([endpoints['lower'][0], endpoints['lower'][1]], wedge_center)

650 for ht_nr, ht in (blk_next.half_turns.items() if is_first_blk_next else reversed(blk_next.half_turns.items())):

651 if ht_lower_angles[ht_nr] > wdg_angle_il: break

652 prev_nr = str(ht_nr)

653 end = prev_nr + 'h'

654 ts_wdg.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)

655

656 # Create auxiliary thin shells for outliers

657 # if both corners belong to thin shells, continue

658 used_wdg_corners = [False, False]

659 for ep in are_endpoints.values():

660 if ep is not None:

661 if ep[1]['higher'] == 'current': used_wdg_corners[1] = True

662 if ep[1]['lower'] == 'current': used_wdg_corners[0] = True

663 if side_case == 'inner':

664 for ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.keys():

665 if ts_name[ts_name.index('_') + 1:] == 'w' + str(wedge_nr):

666 for ep_key, ep in are_endpoints_wdg[int(ts_name[1:ts_name.index('_')])].items():

667 if ep is not None:

668 if ep[1]['higher'] == 'next': used_wdg_corners[1] = True

669 if ep[1]['lower'] == 'next': used_wdg_corners[0] = True

670 else:

671 if wedge_nr in are_endpoints_wdg:

672 for ep in are_endpoints_wdg[wedge_nr].values():

673 if ep is not None:

674 if ep[1]['higher'] == 'current': used_wdg_corners[1] = True

675 if ep[1]['lower'] == 'current': used_wdg_corners[0] = True

676 if not used_wdg_corners[1]:

677 for wdg_nr, wdg in self.geom.wedges.items():

678 if blk_nr_next == wdg.order_l.block: used_wdg_corners[1] = True

679 if not used_wdg_corners[0]:

680 for wdg_nr, wdg in self.geom.wedges.items():

681 if blk_nr_next == wdg.order_h.block: used_wdg_corners[0] = True

682 if not all(used_wdg_corners):

683 def ___create_aux_mid_layer_point(ss, points):

684 mid_layer_ts_aux[mid_layer_name] = dM.Region()

685 circle_pnt = [endpoints[ss][0], endpoints[ss][1]]

686 inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),

687 [[wedge.corrected_center.outer.x, wedge.corrected_center.outer.y], circle_pnt], get_only_closest=True)[0]

688 mid_layer_ts_aux[mid_layer_name].points[str(wedge_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)

689 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)

690 mid_layer_ts_aux[mid_layer_name].lines['w' + str(wedge_nr)] = 0

691 if which_entity['higher'] == 'current' and which_entity['lower'] != 'current':

692 ___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],

693 [wedge_current.oL.x, wedge_current.oL.y]])

694 elif which_entity['higher'] != 'current' and which_entity['lower'] == 'current':

695 ___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],

696 [wedge_current.oH.x, wedge_current.oH.y]])

697 else: # whole block 'within' the facing wedge

698 for wdg_nr, wdg in self.geom.wedges.items():

699 if blk_nr_next == wdg.order_h.block:

700 ___create_aux_mid_layer_point('higher', [[wedge_current.iH.x, wedge_current.iH.y],

701 [wedge_current.oH.x, wedge_current.oH.y]])

702 break

703 elif blk_nr_next == wdg.order_l.block:

704 ___create_aux_mid_layer_point('lower', [[wedge_current.iL.x, wedge_current.iL.y],

705 [wedge_current.oL.x, wedge_current.oL.y]])

706 break

707

708 pole = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole]

709 blk_list_current = list(pole.layers[wedge.order_l.layer].windings[wedge.order_l.winding].blocks.keys())

710 if wedge.order_l.layer < len(pole.layers):

711 __addThinShellPoints('outer', self.md.geometries.thin_shells.mid_layers_wdg_to_ht)

712 if wedge.order_l.layer > 1:

713 __addThinShellPoints('inner', self.md.geometries.thin_shells.mid_layers_ht_to_wdg)

714

715 wedges = self.md.geometries.wedges

716 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux

717 wedge_data = {}

718

719 wdgs_corners = {}

720 for wedge_nr, wedge in self.geom.wedges.items():

721 wdgs_corners[wedge_nr] = {}

722 corners = wdgs_corners[wedge_nr]

723 if wedge.order_l.coil not in wedges.coils:

724 wedges.coils[wedge.order_l.coil] = dM.WedgeLayer()

725 if wedge.order_l.layer not in wedges.coils[wedge.order_l.coil].layers:

726 wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer] = dM.WedgeRegion()

727 wedge_layer = wedges.coils[wedge.order_l.coil].layers[wedge.order_l.layer]

728 wedge_layer.wedges[wedge_nr] = dM.Region()

729 wedge_reg = wedge_layer.wedges[wedge_nr]

730 wedge_layer.block_prev[wedge_nr] = wedge.order_l.block

731 wedge_layer.block_next[wedge_nr] = wedge.order_h.block

732 wnd = self.geom.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[

733 wedge.order_l.layer].windings[wedge.order_l.winding]

734 wnd_next = self.geom.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[

735 wedge.order_h.layer].windings[wedge.order_h.winding]

736 block = wnd.blocks[wedge.order_l.block]

737 block_next = wnd_next.blocks[wedge.order_h.block]

738 corners['last_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_l.coil].poles[wedge.order_l.pole].layers[

739 wedge.order_l.layer].windings[wedge.order_l.winding].blocks[wedge.order_l.block].half_turns.areas.keys())[-1])

740 corners['first_ht'] = int(list(self.md.geometries.coil.coils[wedge.order_h.coil].poles[wedge.order_h.pole].layers[

741 wedge.order_h.layer].windings[wedge.order_h.winding].blocks[wedge.order_h.block].half_turns.areas.keys())[0])

742 ht_current = block.half_turns[corners['last_ht']].corners.bare

743 ht_next = block_next.half_turns[corners['first_ht']].corners.bare

744 d_current = self.data.conductors[wnd.conductor_name].cable.th_insulation_along_width * 2

745 d_next = self.data.conductors[wnd_next.conductor_name].cable.th_insulation_along_width * 2

746 for pnt_close, pnt_far, wdg_corner, d in zip([ht_current.iH, ht_current.oH, ht_next.iL, ht_next.oL],

747 [ht_current.iL, ht_current.oL, ht_next.iH, ht_next.oH],

748 ['il', 'ol', 'ih', 'oh'], [d_current, d_current, d_next, d_next]):

749 if abs(pnt_far.x - pnt_close.x) > 0.:

750 m = (pnt_far.y - pnt_close.y) / (pnt_far.x - pnt_close.x)

751 b = pnt_close.y - m * pnt_close.x

752 root = np.sqrt(- pnt_close.x ** 2 * m ** 2 - 2 * pnt_close.x * b * m + 2 * pnt_close.x * pnt_close.y * m

753 - b ** 2 + 2 * b * pnt_close.y - pnt_close.y ** 2 + d ** 2 * m ** 2 + d ** 2)

754 pnt1_x = (pnt_close.x - b * m + pnt_close.y * m + root) / (m ** 2 + 1)

755 pnt1_y = m * pnt1_x + b

756 pnt2_x = (pnt_close.x - b * m + pnt_close.y * m - root) / (m ** 2 + 1)

757 pnt2_y = m * pnt2_x + b

758 corners[wdg_corner] = Coord(x=pnt1_x, y=pnt1_y) if Func.points_distance([pnt1_x, pnt1_y], [pnt_far.x, pnt_far.y]) >\

759 Func.points_distance([pnt_close.x, pnt_close.y], [pnt_far.x, pnt_far.y]) else Coord(x=pnt2_x, y=pnt2_y)

760 else:

761 bore_cnt_x = self.geom.coil.coils[wedge.order_l.coil].bore_center.x

762 pnt1_y, pnt2_y = pnt_close.y + d, pnt_close.y - d

763 corners[wdg_corner] = Coord(x=pnt_close.x,

764 y=pnt1_y if (wdg_corner[-1] == 'l' and pnt_close.x > bore_cnt_x) or

765 (wdg_corner[-1] == 'h' and pnt_close.x < bore_cnt_x) else pnt2_y)

766 wedge_reg.points[wdg_corner] = self.occ.addPoint(corners[wdg_corner].x, corners[wdg_corner].y, 0)

767 inner = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,

768 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],

769 [corners['ih'].x, corners['ih'].y], [corners['il'].x, corners['il'].y])

770 outer = Func.corrected_arc_center([self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,

771 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y],

772 [corners['oh'].x, corners['oh'].y], [corners['ol'].x, corners['ol'].y])

773 wedge_data[wedge_nr] = [Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']), wedge.corrected_center.outer]

774 wedge_reg.points['inner_center'] = self.occ.addPoint(inner[0], inner[1], 0)

775 wedge_reg.points['outer_center'] = self.occ.addPoint(outer[0], outer[1], 0)

776 wedge_reg.lines['h'] = self.occ.addLine(wedge_reg.points['ih'], wedge_reg.points['oh'])

777 wedge_reg.lines['l'] = self.occ.addLine(wedge_reg.points['il'], wedge_reg.points['ol'])

778 wedge_reg.lines['i'] = self.occ.addCircleArc(wedge_reg.points['ih'], wedge_reg.points['inner_center'], wedge_reg.points['il'])

779 wedge_reg.lines['o'] = self.occ.addCircleArc(wedge_reg.points['oh'], wedge_reg.points['outer_center'], wedge_reg.points['ol'])

780 """

781 logger.warning("Using straight wedge geometry") # required for the projection

782 wedge_reg.lines['i'] = self.occ.addLine(wedge_reg.points['ih'], wedge_reg.points['il'])

783 wedge_reg.lines['o'] = self.occ.addLine(wedge_reg.points['oh'], wedge_reg.points['ol'])

784 """

785 wedge_reg.areas[str(wedge_nr)] = dM.Area(loop=self.occ.addCurveLoop(

786 [wedge_reg.lines['i'], wedge_reg.lines['l'], wedge_reg.lines['o'], wedge_reg.lines['h']]))

787

788 if use_TSA:

789 # Wedge thin shells

790 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_wdg_to_wdg

791 are_endpoints_wdg = {}

792 for coil_nr, coil in self.md.geometries.wedges.coils.items():

793 layer_list = list(coil.layers.keys())

794 for layer_nr, layer in coil.layers.items():

795 if layer_list.index(layer_nr) + 1 < len(layer_list):

796 for wedge_nr, wedge in layer.wedges.items():

797 are_endpoints_wdg[wedge_nr] = {}

798 are_endpoints = are_endpoints_wdg[wedge_nr]

799 wedge_current = wedge_data[wedge_nr][0]

800 wedge_center = [wedge_data[wedge_nr][1].x, wedge_data[wedge_nr][1].y]

801 mean_rad_current = (Func.points_distance([wedge_current.oH.x, wedge_current.oH.y], wedge_center) +

802 Func.points_distance([wedge_current.oL.x, wedge_current.oL.y], wedge_center)) / 2

803 for wdg_next_nr, wdg_next in coil.layers[layer_nr + 1].wedges.items():

804 if self.geom.wedges[wedge_nr].order_l.pole == self.geom.wedges[wdg_next_nr].order_l.pole:

805 wedge_next = wedge_data[wdg_next_nr][0]

806 mean_rad_next = (Func.points_distance([wedge_next.iH.x, wedge_next.iH.y], wedge_center) +

807 Func.points_distance([wedge_next.iL.x, wedge_next.iL.y], wedge_center)) / 2

808 mean_rad = (mean_rad_current + mean_rad_next) / 2

809 mid_layer = self.findMidLayerPoint(wedge_current.oH, wedge_next.iH, wedge_data[wedge_nr][1], mean_rad)

810 are_endpoints[wdg_next_nr] = self.getMidLayerEndpoints(wedge_current, wedge_next, wedge_center, mid_layer_arc_pnt=mid_layer)

811 if are_endpoints[wdg_next_nr]: # this is empty if the wedges are not radially adjacent

812 endpoints = are_endpoints[wdg_next_nr][0]

813 mid_layer_name = 'w' + str(wedge_nr) + '_w' + str(wdg_next_nr)

814 mid_layer_ts[mid_layer_name] = dM.Region()

815 ts = mid_layer_ts[mid_layer_name]

816 ts.points['center'] = self.occ.addPoint(wedge_center[0], wedge_center[1], 0)

817 ts.points['beg'] = self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)

818 end = 'w' + str(wedge_nr if are_endpoints[wdg_next_nr][1] == 'current' else wdg_next_nr)

819 ts.points[end] = self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)

820

821 # Half-turn thin shells

822 for wedge_nr, wedge in self.geom.wedges.items():

823 corners = wdgs_corners[wedge_nr]

824 # Mid layer lines

825 wedge_center = [self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.x,

826 self.md.geometries.coil.coils[wedge.order_l.coil].bore_center.y]

827 _addMidLayerThinShellPoints(Corner(iH=corners['ih'], oH=corners['oh'], iL=corners['il'], oL=corners['ol']))

828 # Mid wedge-turn lines

829 mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn

830 for adj_blk, ht, inner, outer in zip([wedge.order_l, wedge.order_h], [corners['last_ht'], corners['first_ht']],

831 [corners['il'], corners['ih']], [corners['ol'], corners['oh']]):

832 mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)] = dM.Region()

833 ts = mid_turn_ts['w' + str(wedge_nr) + '_' + str(adj_blk.block)]

834 ht_corners = self.geom.coil.coils[adj_blk.coil].poles[adj_blk.pole].layers[

835 adj_blk.layer].windings[adj_blk.winding].blocks[adj_blk.block].half_turns[ht].corners.bare

836 ht_corners_i = ht_corners.iH if ht == corners['last_ht'] else ht_corners.iL

837 ht_corners_o = ht_corners.oH if ht == corners['last_ht'] else ht_corners.oL

838 mid_inner = [(inner.x + ht_corners_i.x) / 2, (inner.y + ht_corners_i.y) / 2]

839 mid_outer = [(outer.x + ht_corners_o.x) / 2, (outer.y + ht_corners_o.y) / 2]

840 line_name = 'w' + str(wedge_nr) + '_' + str(ht)

841 ts.points[line_name + '_i'] = self.occ.addPoint(mid_inner[0], mid_inner[1], 0)

842 ts.points[line_name + '_o'] = self.occ.addPoint(mid_outer[0], mid_outer[1], 0)

843

844 def constructCoilGeometry(self, run_type):

845 """

846 Generates points, hyper lines, and curve loops for the coil half-turns

847 """

848 symmetry = self.data.magnet.geometry.electromagnetics.symmetry if run_type == 'EM' else 'none'

849 # Sub domains angles: first key means 'from 0 to x'; second key means 'from x to 2*pi'

850 if symmetry == 'xy':

851 angle_range = {'to': np.pi / 2, 'from': 2 * np.pi}

852 elif symmetry == 'x':

853 angle_range = {'to': np.pi, 'from': 2 * np.pi}

854 elif symmetry == 'y':

855 angle_range = {'to': np.pi / 2, 'from': 3 / 2 * np.pi}

856 elif symmetry == 'none':

857 angle_range = {'to': 2 * np.pi, 'from': 0}

858 else:

859 raise Exception('Symmetry plane not supported.')

860

861 def _addMidLayerThinShellPoints(pnt_params, ss, name, case):

862 endpnts, cnt = ts_endpoints[name]

863 if len(pnt_params) == 3: # line parameters (cos-theta Rutherford)

864 intersect[name] = Func.intersection_between_arc_and_line(pnt_params, [cnt, endpnts['higher'], endpnts['lower']])

865 if intersect[name]:

866 intersect[name] = intersect[name][0]

867 pnt_angle = Func.arc_angle_between_point_and_abscissa(intersect[name], cnt)

868 elif len(pnt_params) == 4: # points coordinates (cos-theta Mono)

869 wnd_next = list(pole.layers[layer_nr + (1 if case == 'current' else -1)].windings.keys())[

870 list(pole.layers[layer_nr].windings.keys()).index(winding_nr)]

871 blk_next = pole.layers[layer_nr + (1 if case == 'current' else -1)].windings[wnd_next].blocks[

872 int(ts_name[ts_name.index('_') + 1:] if case == 'current' else ts_name[:ts_name.index('_')])]

873 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

874 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)

875 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]

876 pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)

877 pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)

878 pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)

879 intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None

880 else: # point coordinates (block-coil)

881 pnt = [endpnts['higher'][0], pnt_params[1]] if coil.type == 'common-block-coil' else [pnt_params[0], endpnts['higher'][1]]

882 if abs(endpnts['higher'][1]) > 1e-6:

883 pnt_angle = Func.arc_angle_between_point_and_abscissa(pnt, cnt)

884 pnt_angle_h = Func.arc_angle_between_point_and_abscissa(endpnts['higher'], cnt)

885 pnt_angle_l = Func.arc_angle_between_point_and_abscissa(endpnts['lower'], cnt)

886 else:

887 pnt_angle = abs(pnt_params[0])

888 pnt_angle_h = abs(endpnts['higher'][0])

889 pnt_angle_l = abs(endpnts['lower'][0])

890 intersect[name] = pnt if pnt_angle_h > pnt_angle > pnt_angle_l else None

891 if intersect[name]:

892 mid_layer_ts[name].mid_layers.points[str(halfTurn_nr) + ss] = \

893 self.occ.addPoint(intersect[name][0], intersect[name][1], 0)

894 mid_layer_ts[name].point_angles[str(halfTurn_nr) + ss] = Func.sig_dig(pnt_angle)

895 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):

896 intersect[name] = pnt

897 return intersect

898

899 def _addMidLayerThinShellGroup(cl, for_mid_pole=False, mid_coil=False):

900 is_first_blk_next = block_nr_next == list(winding_next.blocks.keys())[0]

901 if 'solenoid' in cl.type:

902 ht_list_next = list(reversed(block_next.half_turns.keys()) if layer_nr % 2 == 0 else list(block_next.half_turns.keys()))

903 elif cl.type == 'reversed-block-coil':

904 ht_list_next = (list(block_next.half_turns.keys()) if not is_first_blk_next else list(reversed(block_next.half_turns.keys())))

905 else:

906 ht_list_next = (list(block_next.half_turns.keys()) if is_first_blk_next else list(reversed(block_next.half_turns.keys())))

907 hh = block_next.half_turns[ht_list_next[-1]].corners.bare

908 ll = block_next.half_turns[ht_list_next[0]].corners.bare

909 bc_next = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)

910 if 'block-coil' in cl.type or (cable_type_curr in ['Mono', 'Ribbon'] and not mid_coil):

911 center = [cl.bore_center.x, cl.bore_center.y]

912 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)

913 else:

914 mean_rad_next = (Func.points_distance([bc_next.iH.x, bc_next.iH.y], [cl.bore_center.x, cl.bore_center.y]) +

915 Func.points_distance([bc_next.iL.x, bc_next.iL.y], [cl.bore_center.x, cl.bore_center.y])) / 2

916 mean_rad = (mean_rad_current + mean_rad_next) / 2

917 mid_layer_h = self.findMidLayerPoint(bc_current.oH, bc_next.iH, cl.bore_center, mean_rad)

918 mid_layer_l = self.findMidLayerPoint(bc_current.oL, bc_next.iL, cl.bore_center, mean_rad)

919 mid_ht_next_i = int(len(ht_list_next) / 2) if len(ht_list_next) % 2 == 0 else round(len(ht_list_next) / 2)

920 mid_ht_next = block_next.half_turns[ht_list_next[mid_ht_next_i - 1]].corners.insulated

921 mid_layer_m = self.findMidLayerPoint(mid_ht_current.oH, mid_ht_next.iH, cl.bore_center, mean_rad)

922 center = Func.arc_center_from_3_points(mid_layer_h, mid_layer_m, mid_layer_l)

923 are_endpoints = self.getMidLayerEndpoints(bc_current, bc_next, center, mid_layer_arc_pnt=mid_layer_h, cable_type=cable_type_curr)

924 if are_endpoints: # this is empty if the blocks are not radially adjacent

925 endpoints = are_endpoints[0]

926 which_block = are_endpoints[1]

927 mid_layer_name = blk_nr + '_' + str(block_nr_next)

928 if for_mid_pole:

929 block_coil_mid_pole_next_blks_list[block_nr_next].append(mid_layer_name)

930 block_coil_ts_endpoints[mid_layer_name] = [endpoints, center]

931 else:

932 if block_nr_next not in list(next_blks_list.keys()):

933 next_blks_list[block_nr_next] = []

934 next_blks_list[block_nr_next].append(mid_layer_name)

935 ts_endpoints[mid_layer_name] = [endpoints, center]

936 mid_layer_ts[mid_layer_name] = dM.MidLayer()

937 mid_layer_ts[mid_layer_name].half_turn_lists[blk_nr] = ht_list

938 mid_layer_ts[mid_layer_name].half_turn_lists[str(block_nr_next)] = ht_list_next

939 beg = (str(ht_list[0]) if which_block['lower'] == 'current' else str(ht_list_next[0])) + 'l'

940 mid_layer_ts[mid_layer_name].mid_layers.points[beg] = \

941 self.occ.addPoint(endpoints['lower'][0], endpoints['lower'][1], 0)

942 end = (str(ht_list[-1]) if which_block['higher'] == 'current' else str(ht_list_next[-1])) + 'h'

943 mid_layer_ts[mid_layer_name].mid_layers.points[end] = \

944 self.occ.addPoint(endpoints['higher'][0], endpoints['higher'][1], 0)

945 if not for_mid_pole or (for_mid_pole and abs(endpoints['higher'][1]) > 1e-6):

946 mid_layer_ts[mid_layer_name].point_angles[beg] =\

947 Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['lower'], center))

948 mid_layer_ts[mid_layer_name].point_angles[end] =\

949 Func.sig_dig(Func.arc_angle_between_point_and_abscissa(endpoints['higher'], center))

950 else:

951 mid_layer_ts[mid_layer_name].point_angles[beg] = abs(endpoints['lower'][0])

952 mid_layer_ts[mid_layer_name].point_angles[end] = abs(endpoints['higher'][0])

953

954 # Create anticlockwise order of blocks

955 present_blocks = []

956 block_corner_angles = {}

957 concentric_coils = self.md.geometries.coil.concentric_coils

958 acw_order = self.md.geometries.coil.anticlockwise_order.coils

959 self.md.geometries.coil.physical_order = self.geom.coil.physical_order

960 for coil_nr, coil in self.geom.coil.coils.items():

961 # if coil_nr not in block_corner_angles:

962 block_corner_angles[coil_nr] = {}

963 if (coil.bore_center.x, coil.bore_center.y) not in concentric_coils:

964 concentric_coils[(coil.bore_center.x, coil.bore_center.y)] = []

965 concentric_coils[(coil.bore_center.x, coil.bore_center.y)].append(coil_nr)

966 for pole_nr, pole in coil.poles.items():

967 for layer_nr, layer in pole.layers.items():

968 if layer_nr not in block_corner_angles[coil_nr]:

969 block_corner_angles[coil_nr][layer_nr] = {}

970 blk_angles = block_corner_angles[coil_nr][layer_nr]

971 for winding_nr, winding in layer.windings.items():

972 for block_nr, block in winding.blocks.items():

973 blk_angles[block_nr] = {'angle': Func.sig_dig(Func.arc_angle_between_point_and_abscissa(

974 [block.block_corners.iL.x, block.block_corners.iL.y],

975 [coil.bore_center.x, coil.bore_center.y])), 'keys': [pole_nr, winding_nr]}

976 higher_angle = Func.sig_dig(Func.arc_angle_between_point_and_abscissa(

977 [block.block_corners.iH.x, block.block_corners.iH.y],

978 [coil.bore_center.x, coil.bore_center.y]))

979 if ((blk_angles[block_nr]['angle'] <= angle_range['to'] and higher_angle <= angle_range['to']) or

980 (angle_range['from'] <= blk_angles[block_nr]['angle'] and angle_range['from'] <= higher_angle)):

981 present_blocks.append(block_nr)

982 for coil_nr, coil in block_corner_angles.items():

983 acw_order[coil_nr] = dM.LayerOrder()

984 for layer_nr, layer in coil.items():

985 acw_order[coil_nr].layers[layer_nr] = []

986 ordered_blocks = [[block_nr, block['angle'], block['keys']] for block_nr, block in layer.items()]

987 ordered_blocks.sort(key=lambda x: x[1])

988 for blk in ordered_blocks:

989 if blk[0] in present_blocks:

990 acw_order[coil_nr].layers[layer_nr].append(dM.AnticlockwiseOrder(pole=blk[2][0], winding=blk[2][1], block=blk[0]))

991

992 # Check if there are concentric coils

993 for bore_center, coils in concentric_coils.items():

994 if len(coils) > 1:

995 radii = []

996 for coil_nr in coils:

997 lyr = self.geom.coil.coils[coil_nr].poles[1].layers[1]

998 blk = list(lyr.windings.keys())[0]

999 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])])

1000 radii.sort(key=lambda x: x[1])

1001 concentric_coils[bore_center] = [rad[0] for rad in radii]

1002

1003 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:

1004 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht

1005 # Collect block couples for block-coil mid-pole thin shells

1006 block_coil_mid_pole_next_blks_list = {}

1007 block_coil_ts_endpoints = {}

1008 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():

1009 if self.geom.coil.coils[coil_nr].type in ['block-coil', 'reversed-block-coil']:

1010 self.block_coil_mid_pole_blks[coil_nr] = []

1011 first_lyr = list(coil.layers.keys())[0]

1012 layer = coil.layers[first_lyr]

1013 for nr, block_order in enumerate(layer):

1014 blk_next_index = nr + 1 if nr + 1 < len(layer) else 0

1015 if layer[blk_next_index].pole != block_order.pole:

1016 self.block_coil_mid_pole_blks[coil_nr].append([block_order, layer[blk_next_index]])

1017 block_coil_mid_pole_next_blks_list[layer[blk_next_index].block] = []

1018 # Mid pole lines for block-coils

1019 for coil_nr, coil in self.block_coil_mid_pole_blks.items():

1020 coil_geom = self.geom.coil.coils[coil_nr]

1021 for mid_pole in coil:

1022 winding = self.geom.coil.coils[coil_nr].poles[mid_pole[0].pole].layers[1].windings[mid_pole[0].winding]

1023 cable_type_curr = self.data.conductors[winding.conductor_name].cable.type

1024 block_nr = mid_pole[0].block

1025 blk_nr = str(block_nr)

1026 block = winding.blocks[block_nr]

1027 is_first_blk = block_nr == list(winding.blocks.keys())[0]

1028 if coil_geom.type == 'reversed-block-coil':

1029 ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))

1030 else:

1031 ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))

1032 hh = block.half_turns[ht_list[-1]].corners.bare

1033 ll = block.half_turns[ht_list[0]].corners.bare

1034 bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)

1035 winding_next = self.geom.coil.coils[coil_nr].poles[mid_pole[1].pole].layers[1].windings[mid_pole[1].winding]

1036 block_nr_next = mid_pole[1].block

1037 block_next = winding_next.blocks[block_nr_next]

1038 _addMidLayerThinShellGroup(coil_geom, for_mid_pole=True)

1039

1040 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux

1041 self.md.geometries.coil.physical_order = self.geom.coil.physical_order

1042 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:

1043 next_blks_list = block_coil_mid_pole_next_blks_list.copy()

1044 ts_endpoints = block_coil_ts_endpoints.copy()

1045 for coil_nr, coil in self.geom.coil.coils.items():

1046 self.md.geometries.coil.coils[coil_nr] = dM.Pole()

1047 coils = self.md.geometries.coil.coils[coil_nr]

1048 coils.type = coil.type

1049 coils.bore_center = coil.bore_center

1050 for pole_nr, pole in coil.poles.items():

1051 coils.poles[pole_nr] = dM.Layer()

1052 poles = coils.poles[pole_nr]

1053 for layer_nr, layer in pole.layers.items():

1054 poles.layers[layer_nr] = dM.Winding()

1055 layers = poles.layers[layer_nr]

1056 for winding_nr, winding in layer.windings.items():

1057 cable_type_curr = self.data.conductors[winding.conductor_name].cable.type

1058 layers.windings[winding_nr] = dM.Block(conductor_name=winding.conductor_name, conductors_number=winding.conductors_number)

1059 windings = layers.windings[winding_nr]

1060 blk_list_current = list(winding.blocks.keys())

1061 for block_nr, block in winding.blocks.items():

1062 if block_nr in present_blocks:

1063 blk_nr = str(block_nr)

1064 windings.blocks[block_nr] = dM.BlockData(current_sign=block.current_sign)

1065 hts = windings.blocks[block_nr].half_turns

1066 is_first_blk = block_nr == list(winding.blocks.keys())[0]

1067 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:

1068 if 'solenoid' in coil.type:

1069 ht_list = (list(reversed(block.half_turns.keys()) if (layer_nr - 1) % 2 == 0 else list(block.half_turns.keys())))

1070 elif coil.type == 'reversed-block-coil':

1071 ht_list = (list(block.half_turns.keys()) if not is_first_blk else list(reversed(block.half_turns.keys())))

1072 else:

1073 ht_list = (list(block.half_turns.keys()) if is_first_blk else list(reversed(block.half_turns.keys())))

1074 hh = block.half_turns[ht_list[-1]].corners.bare

1075 ll = block.half_turns[ht_list[0]].corners.bare

1076 bc_current = Corner(oH=hh.oH, iH=hh.iH, oL=ll.oL, iL=ll.iL)

1077 # Mid layer lines

1078 mean_rad_current = (Func.points_distance([bc_current.oH.x, bc_current.oH.y], [coil.bore_center.x, coil.bore_center.y]) +

1079 Func.points_distance([bc_current.oL.x, bc_current.oL.y], [coil.bore_center.x, coil.bore_center.y])) / 2

1080 mid_ht_current_i = int(len(ht_list) / 2) if len(ht_list) % 2 == 0 else round(len(ht_list) / 2)

1081 mid_ht_current = block.half_turns[ht_list[mid_ht_current_i - 1]].corners.insulated

1082 concentric_coil = concentric_coils[(coil.bore_center.x, coil.bore_center.y)]

1083 if layer_nr < len(pole.layers):

1084 for winding_nr_next, winding_next in pole.layers[layer_nr + 1].windings.items():

1085 if cable_type_curr == 'Rutherford' or\

1086 (cable_type_curr in ['Mono', 'Ribbon'] and

1087 list(pole.layers[layer_nr + 1].windings.keys()).index(winding_nr_next) == list(layer.windings.keys()).index(winding_nr)):

1088 blk_list_next = list(winding_next.blocks.keys())

1089 block_nr_next = blk_list_next[blk_list_current.index(block_nr)]

1090 block_next = winding_next.blocks[block_nr_next]

1091 _addMidLayerThinShellGroup(coil)

1092 elif concentric_coil.index(coil_nr) + 1 < len(concentric_coil):

1093 coil_nr_next = concentric_coil[concentric_coil.index(coil_nr) + 1]

1094 for pole_nr_next, pole_next in self.geom.coil.coils[coil_nr_next].poles.items():

1095 for layer_nr_next, layer_next in pole_next.layers.items():

1096 if layer_nr_next == 1:

1097 for winding_nr_next, winding_next in layer_next.windings.items():

1098 for block_nr_next, block_next in winding_next.blocks.items():

1099 _addMidLayerThinShellGroup(coil, mid_coil=True)

1100 else:

1101 blk_ins = windings.blocks[block_nr].insulation

1102 blk_ins.areas[blk_nr] = dM.Area()

1103

1104 if 'solenoid' in coil.type:

1105 ht_items = (list(reversed(block.half_turns.items()) if layer_nr - 1 % 2 == 0 else list(block.half_turns.items())))

1106 elif coil.type == 'reversed-block-coil':

1107 ht_items = (block.half_turns.items() if not is_first_blk else reversed(block.half_turns.items()))

1108 else:

1109 ht_items = (block.half_turns.items() if is_first_blk else reversed(block.half_turns.items()))

1110 for halfTurn_nr, halfTurn in ht_items:

1111 ht_nr = str(halfTurn_nr)

1112 ht = halfTurn.corners.insulated

1113 hts.areas[ht_nr] = dM.Area()

1114 ht_b = halfTurn.corners.bare

1115

1116 hts.points[ht_nr + 'ih'] = self.occ.addPoint(ht_b.iH.x, ht_b.iH.y, 0)

1117 hts.points[ht_nr + 'il'] = self.occ.addPoint(ht_b.iL.x, ht_b.iL.y, 0)

1118 hts.points[ht_nr + 'oh'] = self.occ.addPoint(ht_b.oH.x, ht_b.oH.y, 0)

1119 hts.points[ht_nr + 'ol'] = self.occ.addPoint(ht_b.oL.x, ht_b.oL.y, 0)

1120

1121 hts.lines[ht_nr + 'i'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'il'])

1122 hts.lines[ht_nr + 'o'] = self.occ.addLine(hts.points[ht_nr + 'oh'], hts.points[ht_nr + 'ol'])

1123 hts.lines[ht_nr + 'l'] = self.occ.addLine(hts.points[ht_nr + 'il'], hts.points[ht_nr + 'ol'])

1124 hts.lines[ht_nr + 'h'] = self.occ.addLine(hts.points[ht_nr + 'ih'], hts.points[ht_nr + 'oh'])

1125

1126 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA:

1127 intersection = {}

1128 # Create mid layer points and compute their angle to the x-axis

1129 for mid_lyr_type in ['current', 'previous']:

1130 for pnt1, pnt2, side in zip(

1131 [[ht_b.iH.x, ht_b.iH.y], [ht_b.iL.x, ht_b.iL.y]],

1132 [[ht_b.oH.x, ht_b.oH.y], [ht_b.oL.x, ht_b.oL.y]], ['h', 'l']):

1133 if (cable_type_curr in ['Mono', 'Ribbon'] and coil.type == 'cos-theta' and

1134 (layer_nr < len(pole.layers) and mid_lyr_type == 'current' or layer_nr > 1 and mid_lyr_type == 'previous')):

1135 pnts_input = pnt1 + pnt2

1136 elif coil.type == 'cos-theta' and (cable_type_curr == 'Rutherford' or cable_type_curr in ['Mono', 'Ribbon'] and\

1137 (layer_nr == len(pole.layers) and mid_lyr_type == 'current' or layer_nr == 1 and mid_lyr_type == 'previous')):

1138 pnts_input = Func.line_through_two_points(pnt1, pnt2)

1139 elif 'block-coil' in coil.type:

1140 pnts_input = pnt1

1141 intersect = {}

1142 if mid_lyr_type == 'current':

1143 # Current mid-layer

1144 for ts_name in ts_endpoints.keys():

1145 if blk_nr == ts_name[:ts_name.index('_')]:

1146 _addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)

1147 elif mid_lyr_type == 'previous':

1148 # Previous mid-layer

1149 if block_nr in next_blks_list:

1150 for ts_name in next_blks_list[block_nr]:

1151 _addMidLayerThinShellPoints(pnts_input, side, ts_name, mid_lyr_type)

1152 for key, value in intersect.items():

1153 if key in intersection:

1154 intersection[key][side] = value

1155 else:

1156 intersection[key] = {side: value}

1157

1158 # Search for half turns that face thin shells only partially

1159 def __create_aux_mid_layer_point(ss, points):

1160 mid_layer_ts_aux[key] = dM.Region()

1161 if 'block-coil' in coil.type:

1162 inter_pnt = [points[0], ts_endpoints[key][0][ss][1]]

1163 else:

1164 inter_pnt = Func.intersection_between_circle_and_line(Func.line_through_two_points(points[0], points[1]),

1165 [ts_endpoints[key][1], ts_endpoints[key][0][ss]], get_only_closest=True)[0]

1166 mid_layer_ts_aux[key].points[str(halfTurn_nr) + ss[0]] = self.occ.addPoint(inter_pnt[0], inter_pnt[1], 0)

1167 mid_layer_ts_aux[key].lines[blk_nr] = 0

1168 for key, value in intersection.items():

1169 first_blk, second_blk = key.split('_')

1170 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

1171 if value['h'] and not value['l']:

1172 __create_aux_mid_layer_point('lower', [ht_b.iL.x, ht_b.iL.y])

1173 elif value['l'] and not value['h']:

1174 __create_aux_mid_layer_point('higher', [ht_b.iH.x, ht_b.iH.y])

1175 else:

1176 relevant_blk = int(first_blk) if second_blk == blk_nr else int(second_blk)

1177 if layer_nr == len(pole.layers) and blk_nr == first_blk:

1178 lyr_blks = acw_order[coil_nr + 1].layers[1]

1179 elif layer_nr == 1 and blk_nr == second_blk:

1180 lyr_blks = acw_order[coil_nr - 1].layers[len(acw_order[coil_nr - 1].layers)]

1181 else:

1182 lyr_blks = acw_order[coil_nr].layers[layer_nr + (1 if first_blk == blk_nr else -1)]

1183 for nr, block_order in enumerate(lyr_blks):

1184 if block_order.block == relevant_blk:

1185 block_order_curr = block_order

1186 block_order_prev = lyr_blks[-1] if nr == 0 else lyr_blks[nr - 1]

1187 block_order_next = lyr_blks[0] if nr + 1 == len(lyr_blks) else lyr_blks[nr + 1]

1188 break

1189 if value['h'] and not value['l'] and block_order_curr.winding == block_order_prev.winding:

1190 __create_aux_mid_layer_point('lower', [[ht_b.iL.x, ht_b.iL.y], [ht_b.oL.x, ht_b.oL.y]])

1191 elif value['l'] and not value['h'] and block_order_curr.winding == block_order_next.winding:

1192 __create_aux_mid_layer_point('higher', [[ht_b.iH.x, ht_b.iH.y], [ht_b.oH.x, ht_b.oH.y]])

1193 else:

1194 blk_ins.points[ht_nr + 'ih'] = self.occ.addPoint(ht.iH.x, ht.iH.y, 0)

1195 blk_ins.points[ht_nr + 'il'] = self.occ.addPoint(ht.iL.x, ht.iL.y, 0)

1196 blk_ins.points[ht_nr + 'oh'] = self.occ.addPoint(ht.oH.x, ht.oH.y, 0)

1197 blk_ins.points[ht_nr + 'ol'] = self.occ.addPoint(ht.oL.x, ht.oL.y, 0)

1198

1199 hts.areas[ht_nr].loop = self.occ.addCurveLoop(

1200 [hts.lines[ht_nr + 'i'], # inner

1201 hts.lines[ht_nr + 'l'], # lower

1202 hts.lines[ht_nr + 'o'], # outer

1203 hts.lines[ht_nr + 'h']]) # higher

1204

1205 # Build wire order of the insulation lines of the current block

1206 if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:

1207 ht_list = list(hts.areas.keys())

1208 ht_list.extend(list(reversed(ht_list))[1:])

1209 self.blk_ins_lines[block_nr] = ['l']

1210 for nr, ht_nr in enumerate(ht_list):

1211 if nr + 1 == winding.conductors_number: # end of first round

1212 self.blk_ins_lines[block_nr].extend([ht_nr + 'i', 'h', ht_nr + 'o'])

1213 else:

1214 if nr + 1 < winding.conductors_number: # within first round

1215 self.blk_ins_lines[block_nr].extend([ht_nr + 'i', ht_nr + 'i' + ht_list[nr + 1]])

1216 else: # within second round

1217 self.blk_ins_lines[block_nr].extend([ht_nr + 'o' + ht_list[nr - 1], ht_nr + 'o'])

1218

1219 def constructInsulationGeometry(self):

1220 """

1221 Generates points, hyper lines, and curve loops for the coil insulations

1222 """

1223 def _createMidPoleLines(case, cnt=0):

1224 if 'block-coil' in geom_coil.type:

1225 if case == 'inner':

1226 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'])

1227 ordered_lines[group_nr].append(['mid_pole_' + case[0], (len(coil.layers) * 2) * 1e3 + 5e2, group.lines['mid_pole_' + case[0]]])

1228 else:

1229 group.lines['mid_pole_' + case[0]] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt_opposite[first_ht_prev + 'il'])

1230 ordered_lines[group_nr].append(['mid_pole_' + case[0], 0, group.lines['mid_pole_' + case[0]]])

1231 else:

1232 ht_curr = geom_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[

1233 block_order.block].half_turns[int(first_ht_curr)].corners.insulated

1234 ht_prev = geom_coil.poles[block_order_prev.pole].layers[layer_nr].windings[block_order_prev.winding].blocks[

1235 block_order_prev.block].half_turns[int(last_ht_prev)].corners.insulated

1236 pnt_curr = [ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y]

1237 pnt_prev = [ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y]

1238 if Func.points_distance(pnt_curr, pnt_prev) > 1e-6:

1239 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],

1240 [ht_curr.iL.x, ht_curr.iL.y] if case == 'inner' else [ht_curr.oL.x, ht_curr.oL.y],

1241 [ht_prev.iH.x, ht_prev.iH.y] if case == 'inner' else [ht_prev.oH.x, ht_prev.oH.y])

1242 ln_name = 'mid_pole_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]

1243 group.lines[ln_name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],

1244 self.occ.addPoint(correct_center[0], correct_center[1], 0),

1245 ins_pnt_opposite[last_ht_prev + case[0] + 'h'])

1246 # self.occ.addLine(ins_pnt[first_ht_curr + case[0] + 'l'], ins_pnt_opposite[last_ht_prev + case[0] + 'h'])

1247 cnt += 1 if case == 'inner' else -1

1248 ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])

1249 return cnt

1250

1251 def _createMidWindingLines(case, cnt):

1252 name = 'mid_wind_' + str(block_order_prev.block) + '_' + str(block_order.block) + '_' + case[0]

1253 # Create corrected center

1254 blk1 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order.block)][0]].layers[

1255 blks_info[str(block_order.block)][1]].windings[blks_info[str(block_order.block)][2]].blocks[int(str(block_order.block))]

1256 blk2 = self.geom.coil.coils[coil_nr].poles[blks_info[str(block_order_prev.block)][0]].layers[

1257 blks_info[str(block_order_prev.block)][1]].windings[blks_info[str(block_order_prev.block)][2]].blocks[int(block_order_prev.block)]

1258 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

1259 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

1260 outer_center = Func.corrected_arc_center([self.md.geometries.coil.coils[coil_nr].bore_center.x,

1261 self.md.geometries.coil.coils[coil_nr].bore_center.y],

1262 [pnt1.x, pnt1.y], [pnt2.x, pnt2.y])

1263 group.lines[name] = self.occ.addCircleArc(ins_pnt[first_ht_curr + case[0] + 'l'],

1264 self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt_opposite[last_ht_prev + case[0] + 'h'])

1265 cnt += 1 if case == 'inner' else -1

1266 ordered_lines[group_nr].append([name, cnt, group.lines[name]])

1267 return cnt

1268

1269 def _createInnerOuterLines(case, cnt):

1270 # Create half turn lines

1271 idxs = [1, round(len(self.blk_ins_lines[block_order.block]) / 2), 1] if case == 'inner'\

1272 else [len(self.blk_ins_lines[block_order.block]) - 1, round(len(self.blk_ins_lines[block_order.block]) / 2), -1]

1273 lns = self.blk_ins_lines[block_order.block][idxs[0]:idxs[1]:idxs[2]]

1274 for ln_nr, ln_name in enumerate(lns):

1275 skip_cnt = False

1276 if ln_name[-1].isdigit():

1277 try:

1278 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + 'h'],

1279 ins_pnt[ln_name[ln_name.index(case[0]) + 1:] + case[0] + 'l'])

1280 except:

1281 skip_cnt = True

1282 next_line = lns[ln_nr + 1]

1283 pos = 'first' if next_line[:-1] == ln_name[:ln_name.index(case[0])] else 'second'

1284 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')

1285 elif ln_name[-1] in ['i', 'o']:

1286 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name + 'l'], ins_pnt[ln_name + 'h'])

1287 else:

1288 group.lines[ln_name] = self.occ.addLine(ins_pnt[ln_name[:ln_name.index(case[0])] + case[0] + ln_name[-1]],

1289 ins_pnt[ln_name[ln_name.index(case[0]) + 1:-1] + case[0] + ln_name[-1]])

1290 if not skip_cnt:

1291 cnt += 1 if case == 'inner' else -1

1292 ordered_lines[group_nr].append([ln_name, cnt, group.lines[ln_name]])

1293 return cnt

1294

1295 def _computePointAngle(case):

1296 points_angles = pa_next if case == 'outer' else pa_prev

1297 current_ht_h = [current_ht.oH.x, current_ht.oH.y] if case == 'outer' else [current_ht.iH.x, current_ht.iH.y]

1298 if ht_nr == 0:

1299 current_ht_l = [current_ht.oL.x, current_ht.oL.y] if case == 'outer' else [current_ht.iL.x, current_ht.iL.y]

1300 if 'block-coil' in geom_coil.type: current_ht_l[1] = 1 if current_ht_l[1] > 0 else -1

1301 points_angles[str(block_order.block) + '_' + ht_name + 'l'] = Func.arc_angle_between_point_and_abscissa(current_ht_l, center)

1302 if ht_nr == len(ht_list) - 1:

1303 name = ht_name + 'h'

1304 coord = current_ht_h

1305 else: # for mid half turns, get the outer corner

1306 next_ht_ins = geom_hts[int(ht_list[ht_nr + 1])].corners.insulated

1307 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]

1308 condition = (Func.points_distance(current_ht_h, center) > Func.points_distance(next_ht, center))\

1309 if case == 'outer' else (Func.points_distance(current_ht_h, center) < Func.points_distance(next_ht, center))

1310 if condition:

1311 name = ht_name + 'h'

1312 coord = current_ht_h

1313 else:

1314 name = ht_list[ht_nr + 1] + 'l'

1315 coord = next_ht

1316 if 'block-coil' in geom_coil.type: coord[1] = 1 if coord[1] > 0 else -1

1317 points_angles[str(block_order.block) + '_' + name] = Func.arc_angle_between_point_and_abscissa(coord, center)

1318

1319 ins = self.md.geometries.insulation

1320 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():

1321 aux_coil = self.md.geometries.coil.coils[coil_nr]

1322 geom_coil = self.geom.coil.coils[coil_nr]

1323 groups = len(geom_coil.poles)

1324 count = {}

1325 ordered_lines = {}

1326 points_angle = {}

1327 blks_info = {}

1328 ending_line = {}

1329 center = [geom_coil.bore_center.x, geom_coil.bore_center.y]

1330 if coil_nr not in ins.coils:

1331 ins.coils[coil_nr] = dM.InsulationGroup()

1332 ins_groups = ins.coils[coil_nr].group

1333 for layer_nr, layer in coil.layers.items():

1334 group_nr = 1

1335 wnd_nr = len(aux_coil.poles[1].layers[layer_nr].windings)

1336 ordered_layer = layer[wnd_nr:] + layer[:wnd_nr] if layer[0].pole != layer[-1].pole else layer

1337 for nr, block_order in enumerate(ordered_layer):

1338 blks_info[str(block_order.block)] = [block_order.pole, layer_nr, block_order.winding]

1339 # Get previous block in anticlockwise order

1340 block_order_prev = ordered_layer[-1] if nr == 0 else ordered_layer[nr - 1]

1341 # Update insulation group

1342 if block_order.winding == block_order_prev.winding:

1343 group_nr = group_nr + 1 if group_nr < groups else 1

1344 # Initialize dicts

1345 if group_nr not in ins_groups:

1346 ins_groups[group_nr] = dM.InsulationRegion()

1347 points_angle[group_nr] = {}

1348 ordered_lines[group_nr] = []

1349 count[group_nr] = [0, (len(coil.layers) + 1) * 1e3]

1350 group = ins_groups[group_nr].ins

1351 ins_groups[group_nr].blocks.append([block_order.pole, layer_nr, block_order.winding, block_order.block])

1352 # Find the wedge

1353 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:

1354 for wdg, blk in self.md.geometries.wedges.coils[coil_nr].layers[layer_nr].block_prev.items():

1355 if blk == block_order_prev.block:

1356 ins_groups[group_nr].wedges.append([layer_nr, wdg])

1357 break

1358 if layer_nr < len(coil.layers):

1359 mid_layer_next = str(layer_nr) + '_' + str(layer_nr + 1)

1360 if mid_layer_next not in points_angle[group_nr]:

1361 points_angle[group_nr][mid_layer_next] = {}

1362 pa_next = points_angle[group_nr][mid_layer_next]

1363 if layer_nr > 1:

1364 mid_layer_prev = str(layer_nr - 1) + '_' + str(layer_nr)

1365 pa_prev = points_angle[group_nr][mid_layer_prev]

1366 # Get point tags of insulation

1367 ins_pnt = aux_coil.poles[block_order.pole].layers[layer_nr].windings[block_order.winding].blocks[

1368 block_order.block].insulation.points

1369 # Get relevant info for line names

1370 first_ht_curr = self.blk_ins_lines[block_order.block][1][:-1]

1371 last_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[

1372 layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[-1]

1373 ins_pnt_opposite = aux_coil.poles[block_order_prev.pole].layers[

1374 layer_nr].windings[block_order_prev.winding].blocks[block_order_prev.block].insulation.points

1375 if 'cos-theta' == geom_coil.type:

1376 # Create lower and higher angle lines

1377 if block_order.winding == block_order_prev.winding:

1378 group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])

1379 ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])

1380 ending_line[group_nr - 1 if group_nr > 1 else groups] =\

1381 [ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]

1382 # Create inner lines of insulation group

1383 if layer_nr == 1:

1384 if block_order.pole != block_order_prev.pole:

1385 count[group_nr][0] = _createMidPoleLines('inner', count[group_nr][0])

1386 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:

1387 count[group_nr][0] = _createMidWindingLines('inner', count[group_nr][0])

1388 count[group_nr][0] = _createInnerOuterLines('inner', count[group_nr][0])

1389 # Create outer lines of insulation group

1390 if layer_nr == len(coil.layers):

1391 if block_order.pole != block_order_prev.pole:

1392 count[group_nr][1] = _createMidPoleLines('outer', count[group_nr][1])

1393 if block_order.pole == block_order_prev.pole and block_order.winding != block_order_prev.winding:

1394 count[group_nr][1] = _createMidWindingLines('outer', count[group_nr][1])

1395 count[group_nr][1] = _createInnerOuterLines('outer', count[group_nr][1])

1396 elif 'block-coil' in geom_coil.type:

1397 last_ht_curr = self.blk_ins_lines[block_order.block][self.blk_ins_lines[block_order.block].index('h') - 1][:-1]

1398 first_ht_prev = list(aux_coil.poles[block_order_prev.pole].layers[layer_nr].windings[

1399 block_order_prev.winding].blocks[block_order_prev.block].half_turns.areas.keys())[0]

1400 # Create lower and higher angle lines

1401 if block_order.winding == block_order_prev.winding:

1402 group.lines[str(layer_nr) + 'l'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])

1403 ordered_lines[group_nr].append([str(layer_nr) + 'l', (len(coil.layers) * 4 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'l']])

1404 ending_line[group_nr - 1 if group_nr > 1 else groups] =\

1405 [ins_pnt_opposite[last_ht_prev + 'ih'], ins_pnt_opposite[last_ht_prev + 'oh']]

1406 group.lines[str(layer_nr) + 'bh'] = self.occ.addLine(ins_pnt[last_ht_curr + 'ih'], ins_pnt[last_ht_curr + 'oh'])

1407 ordered_lines[group_nr].append([str(layer_nr) + 'bh', (len(coil.layers) * 2 + layer_nr) * 1e3, group.lines[str(layer_nr) + 'bh']])

1408 # Create inner lines of insulation group

1409 if block_order.pole != block_order_prev.pole:

1410 if layer_nr == 1:

1411 _createMidPoleLines('inner')

1412 _createMidPoleLines('outer')

1413 group.lines[str(layer_nr) + 'bl'] = self.occ.addLine(ins_pnt[first_ht_curr + 'il'], ins_pnt[first_ht_curr + 'ol'])

1414 ordered_lines[group_nr].append([str(layer_nr) + 'bl', (len(coil.layers) * 2 - layer_nr + 1) * 1e3, group.lines[str(layer_nr) + 'bl']])

1415 # Create outer lines of insulation group

1416 if layer_nr == len(coil.layers):

1417 count[group_nr][1] = _createInnerOuterLines(

1418 'outer', (len(coil.layers) * 4 - layer_nr + 1) * 1e3 if block_order.winding == block_order_prev.winding else (len(coil.layers) + 1) * 1e3)

1419 # Store info about the angle of each point in between layers

1420 ht_list = list(aux_coil.poles[block_order.pole].layers[

1421 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns.areas.keys())

1422 geom_hts = geom_coil.poles[block_order.pole].layers[

1423 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns

1424 for ht_nr, ht_name in enumerate(ht_list): # half turns in anticlockwise order

1425 current_ht = geom_hts[int(ht_name)].corners.insulated

1426 if layer_nr < len(coil.layers): # if it's not the last layer, fetch all outer corners angles

1427 _computePointAngle('outer')

1428 if layer_nr > 1: # if it's not the first layer, fetch all inner corners angles

1429 _computePointAngle('inner')

1430 # Create closing lines

1431 for grp_nr, grp in ending_line.items():

1432 ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h'] = self.occ.addLine(grp[0], grp[1])

1433 ordered_lines[grp_nr].append([str(layer_nr) + 'h', layer_nr * 1e3, ins_groups[grp_nr].ins.lines[str(layer_nr) + 'h']])

1434 # Create lines connecting different layers and generate closed loops

1435 for group_nr, group in points_angle.items():

1436 ins_group = ins_groups[group_nr].ins

1437 for mid_l_name, mid_l in group.items():

1438 first_layer = mid_l_name[:mid_l_name.index('_')]

1439 # Correct angles if the group crosses the abscissa

1440 max_angle = max(mid_l.values())

1441 max_diff = max_angle - min(mid_l.values())

1442 if max_diff > np.pi:

1443 for pnt_name, angle in mid_l.items():

1444 if angle < max_diff / 2:

1445 mid_l[pnt_name] = angle + max_angle

1446 # Order points according to angle

1447 ordered_pnts = [[pnt_name, angle] for pnt_name, angle in mid_l.items()]

1448 ordered_pnts.sort(key=lambda x: x[1])

1449 ordered_names = [x[0] for x in ordered_pnts]

1450 for case in ['beg', 'end']:

1451 past_blocks = []

1452 sides = ['l', 'o', 'h', 'l'] if case == 'beg' else ['h', 'i', 'l', 'h']

1453 # count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2

1454 for i in range(2 if 'block-coil' in geom_coil.type else 1):

1455 count = int(first_layer) * 1e3 + 5e2 if i == 0 else (len(coil.layers) * 2 + int(first_layer)) * 1e3 + 5e2

1456 if case == 'beg':

1457 pnt_position = 0 if i == 0 else int(len(ordered_names) / 2)

1458 else:

1459 pnt_position = -1 if i == 0 else int(len(ordered_names) / 2 - 1)

1460 first_block = ordered_names[pnt_position][:ordered_names[pnt_position].index('_')] # ordered_pnts[pnt_position][0][:ordered_pnts[pnt_position][0].index('_')] #

1461 ordered_search_names = ordered_names[pnt_position::1 if case == 'beg' else -1]

1462 for nr, pnt in enumerate(ordered_search_names[1:], 1): # enumerate(ordered_names if case == 'beg' else reversed(ordered_names)): #

1463 current_blk = pnt[:pnt.index('_')]

1464 ins_pnt = aux_coil.poles[blks_info[current_blk][0]].layers[blks_info[current_blk][1]].windings[

1465 blks_info[current_blk][2]].blocks[int(current_blk)].insulation.points

1466 prev_pnt = ordered_search_names[nr - 1] # ordered_pnts[nr - 1 if case == 'beg' else - nr][0] #

1467 prev_blk = prev_pnt[:prev_pnt.index('_')]

1468 start_pnt_name = prev_pnt[prev_pnt.index('_') + 1:-1] + ('o' if str(blks_info[prev_blk][1]) == first_layer else 'i')

1469 ins_pnt_prev = aux_coil.poles[blks_info[prev_blk][0]].layers[blks_info[prev_blk][1]].windings[

1470 blks_info[prev_blk][2]].blocks[int(prev_blk)].insulation.points

1471 # Create lines when you find the first edge belonging to a block of the opposite layer

1472 if blks_info[current_blk][1] != blks_info[first_block][1]:

1473 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')

1474 pnt_tag_name_opposite = start_pnt_name + ('l' if prev_pnt[-1] == 'l' else 'h')

1475 opp_blk_ins_lines = self.blk_ins_lines[int(prev_blk)]

1476 indexes = [opp_blk_ins_lines.index(start_pnt_name) + (1 if prev_pnt[-1] == sides[0] else 0),

1477 len(opp_blk_ins_lines) if case == 'beg' else opp_blk_ins_lines.index('h'), 1] if start_pnt_name[-1] == sides[1]\

1478 else [opp_blk_ins_lines.index(start_pnt_name) - (1 if prev_pnt[-1] == sides[0] else 0),

1479 0 if case == 'beg' else opp_blk_ins_lines.index('h'), -1]

1480 if case == 'beg':

1481 if i == 0:

1482 count = (len(coil.layers) * (4 if 'block-coil' in geom_coil.type else 2) - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])

1483 else:

1484 count = (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2 - abs(indexes[0] - indexes[1])

1485 else:

1486 count += 1 + abs(indexes[0] - indexes[1])

1487 # Create all remaining lines of the current layer block

1488 for line_name in opp_blk_ins_lines[indexes[0]:indexes[1]:indexes[2]]:

1489 if 'block-coil' in geom_coil.type:

1490 if not line_name[-1].isdigit():

1491 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])

1492 count += 1 if (case == 'beg' and i == 1) or (case == 'end' and i == 0) else -1

1493 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])

1494 else:

1495 if line_name[-1].isdigit():

1496 ins_group.lines[line_name] = self.occ.addLine(

1497 ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'],

1498 ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'])

1499 else:

1500 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'l'], ins_pnt_prev[line_name + 'h'])

1501 count += 1 if case == 'beg' else -1 # if start_pnt_name[-1] == sides[1] else 1

1502 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])

1503 # Create mid layer line

1504 if 'block-coil' in geom_coil.type:

1505 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])

1506 else:

1507 count_rest = -abs(indexes[0] - indexes[1]) if case == 'beg' else 1 + abs(indexes[0] - indexes[1])

1508

1509 line_name = 'mid_layer_' + mid_l_name + ('b' if i == 1 else '') + (

1510 '_l' if case == 'beg' else '_h')

1511

1512 gmsh.model.occ.synchronize()

1513 """

1514 The line that connects two layers may overlap with a new (pole) region. This can be avoided by

1515 modifying the line to be an L shape by adding an intermediate point.

1516 We add the point along the upper half-turn radial direction towards the center.

1517 """

1518 p = ins_pnt[pnt_tag_name] # point to be extended

1519 linetag = gmsh.model.getAdjacencies(0, ins_pnt[pnt_tag_name])[0][0] # line to be extended

1520 # find the distance to move the point towards the center

1521 coord1 = gmsh.model.getValue(0, p, [])

1522 coord2 = gmsh.model.getValue(0, ins_pnt_prev[pnt_tag_name_opposite], [])

1523 # this is the

1524 r1 = np.sqrt(coord1[0] ** 2 + coord1[1] ** 2)

1525 r2 = np.sqrt(coord2[0] ** 2 + coord2[1] ** 2)

1526 distance = (r1 - r2) / 2

1527

1528 p1, p2 = [b[1] for b in gmsh.model.getBoundary([(1, linetag)], oriented=True)]

1529 dir_vector = gmsh.model.getValue(0, p1, [])-gmsh.model.getValue(0, p2, [])

1530 unit_vector = dir_vector / np.linalg.norm(dir_vector)

1531 coord = gmsh.model.getValue(0, p, [])

1532 X = self.occ.addPoint(

1533 coord[0] + unit_vector[0] * distance,

1534 coord[1] + unit_vector[1] * distance,

1535 0)

1536 gmsh.model.occ.synchronize()

1537

1538 ins_group.lines[line_name + "_A"] = self.occ.addLine(ins_pnt[pnt_tag_name], X)

1539 ins_group.lines[line_name + "_B"] = self.occ.addLine(X, ins_pnt_prev[pnt_tag_name_opposite])

1540

1541 ordered_lines[group_nr].append(

1542 [line_name+"_A", count + count_rest -1 if case == 'beg' else count + count_rest, ins_group.lines[line_name+"_A"]])

1543 ordered_lines[group_nr].append(

1544 [line_name+"_B", count + count_rest if case == 'beg' else count + count_rest -1, ins_group.lines[line_name+"_B"]])

1545

1546 """ # original code

1547 line_name = 'mid_layer_' + mid_l_name + ('b' if i == 1 else '') + ('_l' if case == 'beg' else '_h')

1548 ins_group.lines[line_name] = self.occ.addLine(ins_pnt[pnt_tag_name], ins_pnt_prev[pnt_tag_name_opposite])

1549 ordered_lines[group_nr].append([line_name, count + count_rest, ins_group.lines[line_name]])

1550 """

1551 break

1552 # Create all edges of the first block sticking out completely todo: might have to be extended to multiple blocks

1553 if current_blk != first_block and current_blk not in past_blocks:

1554 def __createWedgeInsulation(cnt):

1555 # Create the line connecting the blocks (where a wedge is)

1556 line_name = self.blk_ins_lines[int(current_blk)][

1557 (-1 if start_pnt_name[-1] == 'o' else 1) if case == 'beg'

1558 else (round(len(self.blk_ins_lines[int(current_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1))]

1559 line_name_prev = self.blk_ins_lines[int(prev_blk)][

1560 (round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + (1 if start_pnt_name[-1] == 'o' else -1)) if case == 'beg'

1561 else (-1 if start_pnt_name[-1] == 'o' else 1)]

1562 # Create corrected center

1563 blk1 = geom_coil.poles[blks_info[prev_blk][0]].layers[

1564 blks_info[prev_blk][1]].windings[blks_info[prev_blk][2]].blocks[int(prev_blk)]

1565 blk2 = geom_coil.poles[blks_info[current_blk][0]].layers[

1566 blks_info[current_blk][1]].windings[blks_info[current_blk][2]].blocks[int(current_blk)]

1567 pnt1 = blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oH if case == 'beg'\

1568 else blk1.half_turns[int(line_name_prev[:-1])].corners.insulated.oL

1569 pnt2 = blk2.half_turns[int(line_name[:-1])].corners.insulated.oL if case == 'beg'\

1570 else blk2.half_turns[int(line_name[:-1])].corners.insulated.oH

1571 outer_center = Func.corrected_arc_center([aux_coil.bore_center.x, aux_coil.bore_center.y],

1572 [pnt2.x, pnt2.y] if case == 'beg' else [pnt1.x, pnt1.y],

1573 [pnt1.x, pnt1.y] if case == 'beg' else [pnt2.x, pnt2.y])

1574 ins_group.lines[line_name_prev + line_name] =\

1575 self.occ.addCircleArc(ins_pnt_prev[line_name_prev + sides[2]],

1576 self.occ.addPoint(outer_center[0], outer_center[1], 0), ins_pnt[line_name + sides[3]])

1577 ordered_lines[group_nr].append([line_name_prev + line_name, cnt, ins_group.lines[line_name_prev + line_name]])

1578

1579 count = int(first_layer) * 1e3 + 5e2 if case == 'end' else (len(coil.layers) * 2 - int(first_layer)) * 1e3 + 5e2

1580 past_blocks.append(current_blk)

1581 indexes = [round(len(self.blk_ins_lines[int(prev_blk)]) / 2) + 1,

1582 len(self.blk_ins_lines[int(prev_blk)])] if str(blks_info[prev_blk][1]) == first_layer\

1583 else [1, round(len(self.blk_ins_lines[int(prev_blk)]) / 2)]

1584 if case == 'beg':

1585 count += 1

1586 __createWedgeInsulation(count)

1587 lines = self.blk_ins_lines[int(prev_blk)][indexes[0]:indexes[1]]

1588 side = 'o' if str(blks_info[prev_blk][1]) == first_layer else 'i'

1589 for line_nr, line_name in enumerate(lines):

1590 skip_count = False

1591 if line_name[-1].isdigit():

1592 try:

1593 ins_group.lines[line_name] =\

1594 self.occ.addLine(ins_pnt_prev[line_name[line_name.index(start_pnt_name[-1]) + 1:] + start_pnt_name[-1] + 'l'],

1595 ins_pnt_prev[line_name[:line_name.index(start_pnt_name[-1])] + start_pnt_name[-1] + 'h'])

1596 except: # points are too close to each other

1597 skip_count = True

1598 next_line = lines[line_nr + 1]

1599 pnt1, pnt2 = line_name.split(side)

1600 pos = 'first' if next_line[:-1] == pnt1 else 'second'

1601 lines[line_nr + 1] = next_line + (pnt2 + 'l' if pos == 'first' else pnt1 + 'h')

1602 elif line_name[-1] in ['i', 'o']:

1603 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name + 'h'], ins_pnt_prev[line_name + 'l'])

1604 else:

1605 ins_group.lines[line_name] = self.occ.addLine(ins_pnt_prev[line_name[:line_name.index(side)] + side + line_name[-1]],

1606 ins_pnt_prev[line_name[line_name.index(side) + 1:-1] + side + line_name[-1]])

1607 if not skip_count:

1608 count += 1 # if start_pnt_name[-1] == sides[1] else -1

1609 ordered_lines[group_nr].append([line_name, count, ins_group.lines[line_name]])

1610 if case == 'end':

1611 count += 1

1612 __createWedgeInsulation(count)

1613

1614 # Generate closed loops

1615 ordered_lines[group_nr].sort(key=lambda x: x[1])

1616 area_name = str((coil_nr - 1) * len(ins_groups) + group_nr)

1617 ins_group.areas[area_name] = dM.Area()

1618 if len(points_angle) == 1:

1619 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]]

1620 if 'i' in line and line[0].isdigit() or '_i' in line]))

1621 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]]

1622 if 'o' in line and line[0].isdigit() or '_o' in line])

1623 else:

1624 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]]])

1625

1626 def constructThinShells_poles(self):

1627 ts_layer = self.md.geometries.thin_shells.pole_layers

1628 ts_av_ins_thick = self.md.geometries.thin_shells.ins_thickness.poles

1629

1630 def _construct_thin_shell_corners_to_line(pnt1, pnt2, pole_line, name):

1631 # use gmsh to calculate distance to a line

1632 coord_a = gmsh.model.getClosestPoint(1, pole_line, coord=(pnt1[0], pnt1[1], 0))[0]

1633 coord_b = gmsh.model.getClosestPoint(1, pole_line, coord=(pnt2[0], pnt2[1], 0))[0]

1634 # draw new point at half the distance between iH and coord_a

1635 new_i = self.occ.addPoint((pnt1[0] + coord_a[0]) / 2, (pnt1[1] + coord_a[1]) / 2, 0)

1636 new_o = self.occ.addPoint((pnt2[0] + coord_b[0]) / 2, (pnt2[1] + coord_b[1]) / 2, 0)

1637

1638 ts_layer[name] = dM.Region()

1639

1640 self.occ.synchronize()

1641 ts_layer[name].lines['1'] = self.occ.addLine(new_i, new_o)

1642 self.occ.synchronize()

1643

1644 cond_name = next(iter(self.data.conductors.keys()))

1645 other_material = 0.5*(self.data.conductors[cond_name].cable.th_insulation_along_height+self.data.conductors[cond_name].cable.th_insulation_along_width)#todo, better select which one

1646 # distance -> Average between coord_a and iH and coord_b and oH AND remove the G10 thickness

1647 ts_av_ins_thick[name] = float(0.5 * (np.sqrt((pnt1[0] - coord_a[0]) ** 2 + (pnt1[1] - coord_a[1]) ** 2) +

1648 np.sqrt((pnt2[0] - coord_b[0]) ** 2 + (pnt2[1] - coord_b[1]) ** 2))) - other_material

1649 return

1650

1651 def _find_line_closest_to_points(pnt1, pnt2, line_list_tags):

1652 """

1653 Should work for any (pole) geometry. Given the half turn corner points pnt1 = [x1, y1], pnt2 = [x2, y2],

1654 and a list of pole line tags, we need to select which one is on the opposite side to construct the thin shell line

1655 thus, we search the closest line(s) to each corner point and then take the intersection of those two sets

1656 """

1657 closest_lines = {0: [], 1: []}

1658 for i, pnt in enumerate([pnt1, pnt2]):

1659 min_dist = float('inf')

1660 for line_tag in line_list_tags:

1661 tag1, tag2 = gmsh.model.getAdjacencies(1, line_tag)[1]

1662 start_pnt = gmsh.model.getValue(0, tag1, [])

1663 end_pnt = gmsh.model.getValue(0, tag2, [])

1664 v = end_pnt - start_pnt

1665 w = pnt - start_pnt

1666 c1 = np.dot(w, v)

1667 c2 = np.dot(v, v)

1668 # avoid extending the line

1669 if c1 <= 0:

1670 dist = np.linalg.norm(pnt - start_pnt) # Closest to p1

1671 elif c2 <= c1:

1672 dist = np.linalg.norm(pnt - end_pnt) # Closest to p2

1673 else:

1674 b = c1 / c2

1675 Pb = start_pnt + b * v

1676 dist = np.linalg.norm(pnt - Pb)

1677

1678 if np.isclose(min_dist, dist):

1679 closest_lines[i].append(line_tag) # add to list

1680 elif dist < min_dist:

1681 min_dist = dist

1682 closest_lines[i] = [line_tag] # overwrite

1683

1684 # return the intersection of closest lines

1685 return list(set(closest_lines[0]).intersection(set(closest_lines[1])))[0] # should be only one line

1686

1687 def _split_lines_azimuthal_radial(line_tag_list):

1688 alines = []

1689 rlines = []

1690 for tag in line_tag_list:

1691 pointTags = gmsh.model.getAdjacencies(1, tag)[1] # get tag

1692 p1 = gmsh.model.getValue(0, pointTags[0], [])

1693 p2 = gmsh.model.getValue(0, pointTags[1], [])

1694

1695 dr = np.sqrt(p2[0] ** 2 + p2[1] ** 2) - np.sqrt(p1[0] ** 2 + p1[1] ** 2)

1696 dt = (np.arctan2(p2[1], p2[0]) - np.arctan2(p1[1], p1[0])) * np.sqrt(p2[0] ** 2 + p2[1] ** 2) ## convert to length

1697 if np.abs(dt)>np.abs(dr):

1698 alines.append(tag)

1699 else:

1700 rlines.append(tag)

1701 return alines, rlines

1702

1703 def _wedge_to_pole_lines():

1704 """

1705 These lines do not exist, those that already exist contain the G10 insulation so we need more.

1706 """

1707 for _, region in self.md.geometries.thin_shells.mid_layers_aux.items(): #those are the g10 layers, we need to duplicate this line to account for the kapton

1708 # obtain the position of the lines and draw more on top of it

1709 for name, line_tag in region.lines.items():

1710 name = f"p{name}_r" # all radial lines

1711 ts_layer[name] = dM.Region()

1712

1713 point_tag = gmsh.model.getBoundary([(1, line_tag)], oriented=False)

1714 p1, p2 = point_tag[0][1], point_tag[1][1]

1715 x1, y1, _ = gmsh.model.getValue(0, p1, [])

1716 x2, y2, _ = gmsh.model.getValue(0, p2, [])

1717 # create new line

1718 p1_new = gmsh.model.occ.addPoint(x1, y1, 0.0)

1719 p2_new = gmsh.model.occ.addPoint(x2, y2, 0.0)

1720 self.occ.synchronize()

1721 ts_layer[name].lines['1'] = self.occ.addLine(p1_new, p2_new)

1722 self.occ.synchronize()

1723

1724 # thickness = distance between the wedge and the pole - G10 thickness

1725 # how to approximate this thickness? -> use the same thickness as for the ht 108

1726 """

1727 cond_name = next(iter(self.data.conductors.keys()))

1728 other_material = 0.5 * (

1729 self.data.conductors[cond_name].cable.th_insulation_along_height + self.data.conductors[

1730 cond_name].cable.th_insulation_along_width)

1731 # distance -> Average between coord_a and iH and coord_b and oH AND remove the G10 thickness ?

1732 ts_av_ins_thick[name] = ???

1733 """

1734 ts_av_ins_thick[name] = ts_av_ins_thick['p108_a'] #@emma hardcoded thickness, use the same as for this ht

1735

1736 max_layer = max(self.geom.coil.coils[1].poles[1].layers.keys())

1737 # first, the lines connecting hts to the pole

1738 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items(): # coilnr is only 1 here

1739 for layer_nr, layer in coil.layers.items(): # we need to add radial TS lines for both layers

1740 for nr, blk_order in enumerate(layer):

1741 block = self.geom.coil.coils[coil_nr].poles[blk_order.pole].layers[

1742 layer_nr].windings[blk_order.winding].blocks[blk_order.block]

1743 ht_list = list(self.md.geometries.coil.coils[coil_nr].poles[blk_order.pole].layers[

1744 layer_nr].windings[blk_order.winding].blocks[

1745 blk_order.block].half_turns.areas.keys())

1746

1747 blk_index_next = nr + 1 if nr + 1 < len(layer) else 0

1748 block_order_next = layer[blk_index_next]

1749 block_next = self.geom.coil.coils[coil_nr].poles[block_order_next.pole].layers[

1750 layer_nr].windings[block_order_next.winding].blocks[block_order_next.block]

1751 ht_list_next = list(self.md.geometries.coil.coils[coil_nr].poles[block_order_next.pole].layers[

1752 layer_nr].windings[block_order_next.winding].blocks[

1753 block_order_next.block].half_turns.areas.keys())

1754 ht_last = int(ht_list[-1])

1755 ht_next_first = int(ht_list_next[0])

1756 # if winding nr is the same and pole is the same -> pole connection

1757

1758 if blk_order.pole == block_order_next.pole and blk_order.winding == block_order_next.winding:

1759 # Create radial lines for pole connection

1760 pole_lines_a, pole_lines_r = _split_lines_azimuthal_radial(

1761 [v for qq in self.md.geometries.poles.quadrants.keys() for v in

1762 self.md.geometries.poles.quadrants[qq].lines.values()])

1763 # todo maybe this can be speed up if we select the correct quadrant

1764 if layer_nr != max_layer:

1765 # add azimuthal thin shells at the ht side 'o' (normals radial)

1766 # Now we consider all the HT's from one block, this might not hold true in general

1767

1768 for ht in map(int, ht_list):

1769 oH = [block.half_turns[ht].corners.bare.oH.x,

1770 block.half_turns[ht].corners.bare.oH.y, 0.0]

1771 oL = [block.half_turns[ht].corners.bare.oL.x,

1772 block.half_turns[ht].corners.bare.oL.y, 0.0]

1773 _construct_thin_shell_corners_to_line(

1774 oH, oL,

1775 _find_line_closest_to_points(oH, oL, pole_lines_a),

1776 name=f"p{ht}_r")

1777

1778 for ht in map(int, ht_list_next):

1779 oH = [block_next.half_turns[ht].corners.bare.oH.x,

1780 block_next.half_turns[ht].corners.bare.oH.y, 0.0]

1781 oL = [block_next.half_turns[ht].corners.bare.oL.x,

1782 block_next.half_turns[ht].corners.bare.oL.y, 0.0]

1783 _construct_thin_shell_corners_to_line(

1784 oH, oL,

1785 _find_line_closest_to_points(oH, oL, pole_lines_a),

1786 name=f"p{ht}_r")

1787

1788 iH = [block.half_turns[ht_last].corners.bare.iH.x,

1789 block.half_turns[ht_last].corners.bare.iH.y, 0.0]

1790 oH = [block.half_turns[ht_last].corners.bare.oH.x,

1791 block.half_turns[ht_last].corners.bare.oH.y, 0.0]

1792 _construct_thin_shell_corners_to_line(

1793 iH, oH,

1794 _find_line_closest_to_points(iH, oH, pole_lines_r),

1795 name=f"p{ht_last}_a")

1796

1797 iL = [block_next.half_turns[ht_next_first].corners.bare.iL.x,

1798 block_next.half_turns[ht_next_first].corners.bare.iL.y, 0.0]

1799 oL = [block_next.half_turns[ht_next_first].corners.bare.oL.x,

1800 block_next.half_turns[ht_next_first].corners.bare.oL.y, 0.0]

1801 _construct_thin_shell_corners_to_line(

1802 iL, oL,

1803 _find_line_closest_to_points(iL, oL, pole_lines_r),

1804 name=f"p{ht_next_first}_a")

1805 # second, the lines connecting wedge to the pole

1806 _wedge_to_pole_lines()

1807 def constructAdditionalThinShells(self):

1808 def _create_lines_ht_alignment(ohx, ohy, olx, oly, ihx, ihy, ilx, ily):

1809 def __line_circle_intersection(p1, p2):

1810 x1, y1 = p1

1811 x2, y2 = p2

1812 dx = x2 - x1

1813 dy = y2 - y1

1814

1815 A = dx**2 + dy**2

1816 B = 2 * (dx*x1 + dy*y1)

1817 C = x1**2 + y1**2 - R**2 # R collar is known from outer scope

1818

1819 disc = B**2 - 4*A*C

1820 if disc < 0:

1821 logger.warning(" No intersection between line and circle.")

1822 return [] # no intersection

1823 else:

1824 sqrt_disc = np.sqrt(disc)

1825 t1 = (-B + sqrt_disc) / (2*A)

1826 t2 = (-B - sqrt_disc) / (2*A)

1827 return [(x1 + t1*dx, y1 + t1*dy),

1828 (x1 + t2*dx, y1 + t2*dy)]

1829

1830 mid_bare_o = np.mean([[ohx, ohy], [olx, oly]], axis=0)

1831 mid_bare_i = np.mean([[ihx, ihy], [ilx, ily]], axis=0)

1832 offsets = np.array([[ohx, ohy], [olx, oly]]) - mid_bare_o

1833 quad = 1 if mid_bare_o[0] >= 0 and mid_bare_o[1] >= 0 else (

1834 2 if mid_bare_o[0] <= 0 <= mid_bare_o[1] else (3 if mid_bare_o[0] <= 0 and mid_bare_o[1] <= 0 else 4))

1835 inters = __line_circle_intersection(mid_bare_o, mid_bare_i)

1836 quadrants = {

1837 1: lambda x, y: x >= 0 and y >= 0,

1838 2: lambda x, y: x <= 0 <= y,

1839 3: lambda x, y: x <= 0 and y <= 0,

1840 4: lambda x, y: x >= 0 >= y,

1841 }

1842 check = quadrants[quad]

1843 for x, y in inters:

1844 if check(x, y): # select the correct intersection based on the quadrant

1845 xi, yi = x, y

1846 break

1847

1848 # add a point halfway between (x,y) and (xi, yi)

1849 dr = float(np.sqrt((xi - mid_bare_o[0]) ** 2 + (yi - mid_bare_o[1]) ** 2))

1850 new_point_coords = [mid_bare_o[0] + 0.5 * (xi - mid_bare_o[0]), mid_bare_o[1] + 0.5 * (yi - mid_bare_o[1])]

1851 for offset in offsets:

1852 point_tag_final = self.occ.addPoint(new_point_coords[0] + offset[0], new_point_coords[1] + offset[1], 0)

1853 self.occ.synchronize()

1854 return dr, point_tag_final

1855 def _embed_points_to_collar_curve():

1856 def __add_point_in_curve(points, curve_tag):

1857 return self.occ.fragment([(0, k) for k in points], [(1, curve_tag)], removeObject=True)[0]

1858

1859 ### First, cut the collar line

1860 self.occ.synchronize()

1861 new_tags = {1: [], 2: [], 3: [], 4: []} # tuples

1862 new_point_tags = {1: [], 2: [], 3: [], 4: []} # only tags

1863 collar_size = self.data.magnet.mesh.thermal.collar.SizeMin

1864 for ts_name, ts in self.md.geometries.thin_shells.collar_layers.items():

1865 for name, line in ts.lines.items():

1866 coords = gmsh.model.getValue(1, line, [i[0] for i in gmsh.model.getParametrizationBounds(1, line)])

1867 t1 = np.arctan2(coords[0], coords[1])

1868 t2 = np.arctan2(coords[3], coords[4])

1869 quad = t1 // (np.pi / 2) + 1 if t1 > 0 else 4 + t1 // (np.pi / 2) + 1

1870 curve_tag = self.md.geometries.collar.inner_boundary_tags[quad][0]

1871 start, end = min(t1 % (np.pi / 2), t2 % (np.pi / 2)), max(t1 % (np.pi / 2),

1872 t2 % (np.pi / 2)) # ensure start < end

1873 tmp_coords = gmsh.model.getValue(1, line, [i[0] for i in gmsh.model.getParametrizationBounds(1, line)])

1874 target_size = collar_size

1875 elements = max(1, round(Func.points_distance(tmp_coords[:2], tmp_coords[3:-1]) / target_size))+1

1876 if elements%2 == 1: elements += 1

1877 para_coords = np.linspace(start, end, elements, endpoint=True)

1878

1879 for u in para_coords:

1880 if quad == 1 or quad == 3:

1881 u = np.pi / 2 - u ## magic

1882 x, y, z = gmsh.model.getValue(1, curve_tag, [u]) # Evaluate point on curve

1883 new_point_tags[quad].append(self.occ.addPoint(x, y, z)) # Add point coordinates to the list

1884

1885 for q in new_point_tags.keys():

1886 new_tags[q].extend(__add_point_in_curve(new_point_tags[q], self.md.geometries.collar.inner_boundary_tags[q][0]))

1887 # so for occ the old curve no longer exists, but in the objects the tag is still there

1888

1889 REMOVED_TAGS = [self.md.geometries.collar.inner_boundary_tags[q][0] for q in new_tags.keys()]

1890

1891 # update boundary tags

1892 for quad, taglist in new_tags.items():

1893 curvelist = [tag[1] for tag in taglist if tag[0]==1 ] # select the curves only

1894 self.md.geometries.collar.inner_boundary_tags[quad] = curvelist # replace

1895 for k, new_tag in enumerate(self.md.geometries.collar.inner_boundary_tags[quad]):

1896 self.md.geometries.collar.quadrants[quad].lines[str(k)] = new_tag ## adding new lines

1897

1898 # REDEFINE THE CURVELOOP

1899 loop_list = []

1900 tmpdict = copy.deepcopy(self.md.geometries.collar.quadrants[quad].lines)

1901 for tag_name, tag in tmpdict.items():

1902 if tag in self.md.geometries.collar.cooling_tags: # skip holes

1903 continue

1904 if tag in REMOVED_TAGS:

1905 # remove from dictionary

1906 del self.md.geometries.collar.quadrants[quad].lines[tag_name]

1907 continue

1908 loop_list.append(tag)

1909 self.occ.synchronize()

1910

1911 for area in self.md.geometries.collar.quadrants[quad].areas:

1912 if area.startswith('arc') and not area.startswith('arch'): #collar region, not the holes

1913 self.md.geometries.collar.quadrants[quad].areas[area] = dM.Area(

1914 loop=self.occ.addCurveLoop(loop_list))

1915

1916 # Cut the pole lines

1917 """

1918 idea, loop over the thin shell lines, then draw line from the HT corner to the line ends, then intersect it with the pole curve

1919 just find intersection and select the shortest one ?

1920

1921 # not implemented, maybe not necessary

1922 """

1923 for quad in [1, 2, 3, 4]:

1924 # additionally add the pole area back

1925 for area in self.md.geometries.poles.quadrants[quad].areas:

1926 if area.startswith('arp'):

1927 self.md.geometries.poles.quadrants[quad].areas[area] = dM.Area(

1928 loop=self.occ.addCurveLoop(list(self.md.geometries.poles.quadrants[quad].lines.values())))

1929

1930 self.occ.synchronize()

1931

1932 # point still needs to be added

1933 ts_layer = self.md.geometries.thin_shells.collar_layers

1934 ts_av_ins_thick = self.md.geometries.thin_shells.ins_thickness.collar

1935 collar_tag_dict = self.md.geometries.collar.inner_boundary_tags

1936 enforce_TSA_mapping_collar = self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping # if True, cut the collar curve into segments, otherwise use the whole curve

1937

1938 center = self.occ.addPoint(0, 0, 0)

1939 collar_x, collar_y, _ = gmsh.model.getValue(1, collar_tag_dict[1][0], [0])

1940 R = np.sqrt(collar_x ** 2 + collar_y ** 2) # radius of collar curve

1941

1942 alignment = ['radial', 'ht'][1] # pick ht alignment

1943

1944 # COLLAR

1945 for pid, pole in self.geom.coil.coils[1].poles.items():

1946 layer_num = max(pole.layers.keys()) # outside layer

1947 for wid, winding in pole.layers[layer_num].windings.items():

1948 for block_idx in winding.blocks.keys():

1949 block = winding.blocks[block_idx]

1950 ht_nr_area = list(self.md.geometries.coil.coils[1].poles[pid].layers[layer_num].windings[wid].blocks[block_idx].half_turns.areas.keys())

1951 i = 0

1952 for ht_nr, ht in block.half_turns.items(): #ht_idx is not the same as number

1953 ht_old = ht_nr_area[i]

1954 i+=1

1955 dr = 0.

1956 if alignment == 'radial':

1957 for (x, y), (x1, y1) in zip([[ht.corners.bare.oH.x, ht.corners.bare.oH.y], [ht.corners.bare.oL.x, ht.corners.bare.oL.y]],

1958 [[ht.corners.bare.iH.x, ht.corners.bare.iH.y], [ht.corners.bare.iL.x, ht.corners.bare.iL.y]]): # uses the bare coordinates of the HT

1959 t = np.arctan2(y, x)

1960 quad = t // (np.pi / 2) + 1 if t > 0 else 4 + t // (np.pi / 2) + 1

1961 collar_idx = collar_tag_dict[quad][0] # get the (debug: ONE) collar curve tag for the quadrant

1962

1963 dr_prev = dr

1964

1965 dummy = self.occ.addPoint(x, y, 0)

1966 dr = self.occ.get_distance(0, dummy, 1, collar_idx)[0]

1967 self.occ.synchronize()

1968 self.occ.remove([(0, dummy)])

1969 point_tag_final = self.occ.addPoint(x + 0.5 * dr * np.cos(t), y + 0.5 * dr * np.sin(t), 0)

1970

1971 elif alignment == 'ht': # distance to the next half turn (more relevant for coil to coil distances)

1972 dr, point_tag_final = _create_lines_ht_alignment(ohx = ht.corners.bare.oH.x, ohy=ht.corners.bare.oH.y,

1973 olx = ht.corners.bare.oL.x, oly=ht.corners.bare.oL.y,

1974 ihx = ht.corners.bare.iH.x, ihy=ht.corners.bare.iH.y,

1975 ilx = ht.corners.bare.iL.x, ily=ht.corners.bare.iL.y)

1976 # save the line

1977 name = f'{ht_nr}_x'

1978 ts_layer[name] = dM.Region()

1979 self.occ.synchronize()

1980 ts_layer[name].lines['1'] = self.occ.addLine(point_tag_final-1, point_tag_final)

1981 cond_name = next(iter(self.data.conductors.keys()))

1982 if alignment == 'radial':

1983 ts_av_ins_thick[name] = 0.5*(dr+dr_prev) - self.data.conductors[cond_name].cable.th_insulation_along_width # approx thickness, average - insulation

1984 elif alignment == 'ht':

1985 ts_av_ins_thick[name] = dr - self.data.conductors[cond_name].cable.th_insulation_along_width

1986 self.occ.synchronize()

1987

1988 # WEDGES

1989 if self.data.magnet.geometry.thermal.with_wedges:

1990 for wedge_idx, wedge in self.md.geometries.wedges.coils[1].layers[max(self.md.geometries.wedges.coils[1].layers.keys())].wedges.items():

1991 dr=0.

1992 if alignment == 'radial':

1993 raise NotImplementedError("Wedge radial alignment not implemented")

1994

1995 elif alignment == 'ht': # distance to the next half turn (more relevant for coil to coil distances)

1996 ohx, ohy, _ = gmsh.model.getValue(0, wedge.points['oh'], [])

1997 olx, oly, _ = gmsh.model.getValue(0, wedge.points['ol'], [])

1998 ihx, ihy, _ = gmsh.model.getValue(0, wedge.points['ih'], [])

1999 ilx, ily, _ = gmsh.model.getValue(0, wedge.points['il'], [])

2000 dr, point_tag_final = _create_lines_ht_alignment(ohx=ohx, ohy=ohy, olx=olx, oly=oly,

2001 ihx=ihx, ihy=ihy, ilx=ilx, ily=ily)

2002

2003 # find the smallest thickness

2004 dr_b = 0.0

2005 for x, y in [[ohx, ohy],

2006 [olx, oly]]:

2007 dummy = self.occ.addPoint(x, y, 0)

2008 dr_a = dr_b

2009 t = np.arctan2(y, x)

2010 quad = t // (np.pi / 2) + 1 if t > 0 else 4 + t // (np.pi / 2) + 1

2011 collar_idx = collar_tag_dict[quad][0]

2012 dr_b = self.occ.get_distance(0, dummy, 1, collar_idx)[0]

2013 self.occ.remove([(0, dummy)])

2014

2015 # save line

2016 name = f'w{wedge_idx}_x'

2017 ts_layer[name] = dM.Region()

2018 ts_layer[name].lines['1'] = self.occ.addLine(point_tag_final - 1, point_tag_final) # make index line trivial (no distinction)

2019

2020 cond_name = next(iter(self.data.conductors.keys()))

2021 #### ts_av_ins_thick[name] = dr - self.data.conductors[cond_name].cable.th_insulation_along_width

2022 # This overshoots the thickness. The wedges are curved and the center between the corners (straight line) is further away from the collar than the curved wedge boundary.

2023 # Maybe one could use gmsh to obtain the distance between the outer curve and the collar, but this should be close to taking the minimum distance of the cornerpoints as both curves are

2024 # circle segments of (approximately) two concentric circles around the centre

2025

2026 logger = logging.getLogger('FiQuS')

2027 logger.warning("Using alternative wedge insulation thickness approximation ")

2028 ts_av_ins_thick[name] = min(dr_a,dr_b)- self.data.conductors[cond_name].cable.th_insulation_along_width # approx thickness, average - insulation

2029

2030 # POLES

2031 if 'poles' in self.data.magnet.geometry.thermal.areas:

2032 # generate additional thin shell lines

2033 self.constructThinShells_poles()

2034

2035 #gmsh.fltk.run() constructed correctly

2036

2037 if enforce_TSA_mapping_collar:

2038 self.occ.synchronize()

2039 _embed_points_to_collar_curve() ## both coils and wedges

2040

2041 self.occ.remove([(0, center)]) # remove center point

2042 self.occ.synchronize()

2043 def constructThinShells(self, with_wedges):

2044 # default

2045 ins_th = self.md.geometries.thin_shells.ins_thickness

2046 mid_pole_ts = self.md.geometries.thin_shells.mid_poles

2047 mid_winding_ts = self.md.geometries.thin_shells.mid_windings

2048 mid_turn_ts = self.md.geometries.thin_shells.mid_turn_blocks

2049 # not default

2050 mid_layer_ts = self.md.geometries.thin_shells.mid_layers_ht_to_ht

2051 mid_layer_ts_aux = self.md.geometries.thin_shells.mid_layers_aux

2052

2053 # Create mid-pole and mid-turn thin shells

2054 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():

2055 for layer_nr, layer in coil.layers.items():

2056 for nr, blk_order in enumerate(layer):

2057 block = self.geom.coil.coils[coil_nr].poles[blk_order.pole].layers[

2058 layer_nr].windings[blk_order.winding].blocks[blk_order.block]

2059 ht_list = list(self.md.geometries.coil.coils[coil_nr].poles[blk_order.pole].layers[

2060 layer_nr].windings[blk_order.winding].blocks[blk_order.block].half_turns.areas.keys())

2061 # Create mid-pole and mid-winding thin shells

2062 blk_index_next = nr + 1 if nr + 1 < len(layer) else 0

2063 block_order_next = layer[blk_index_next]

2064 block_next = self.geom.coil.coils[coil_nr].poles[block_order_next.pole].layers[

2065 layer_nr].windings[block_order_next.winding].blocks[block_order_next.block]

2066 ht_list_next = list(self.md.geometries.coil.coils[coil_nr].poles[block_order_next.pole].layers[

2067 layer_nr].windings[block_order_next.winding].blocks[block_order_next.block].half_turns.areas.keys())

2068 ht_last = int(ht_list[-1])

2069 ht_next_first = int(ht_list_next[0])

2070 iH = [block.half_turns[ht_last].corners.bare.iH.x, block.half_turns[ht_last].corners.bare.iH.y]

2071 iL = [block_next.half_turns[ht_next_first].corners.bare.iL.x, block_next.half_turns[ht_next_first].corners.bare.iL.y]

2072 oH = [block.half_turns[ht_last].corners.bare.oH.x, block.half_turns[ht_last].corners.bare.oH.y]

2073 oL = [block_next.half_turns[ht_next_first].corners.bare.oL.x, block_next.half_turns[ht_next_first].corners.bare.oL.y]

2074 ts_name = str(blk_order.block) + '_' + str(block_order_next.block)

2075

2076 for ts, th, condition in zip([mid_pole_ts, mid_winding_ts], [ins_th.mid_pole, ins_th.mid_winding],

2077 # ['_ly' + str(layer_nr), '_wd' + str(blk_order.winding) + '_wd' + str(block_order_next.winding)],

2078 [self.geom.coil.coils[coil_nr].type == 'cos-theta' and block_order_next.pole != blk_order.pole,

2079 (not with_wedges or not self.geom.wedges) and self.geom.coil.coils[coil_nr].type in

2080 ['cos-theta', 'common-block-coil'] and block_order_next.pole == blk_order.pole and block_order_next.winding != blk_order.winding]):

2081 if condition:

2082 ts[ts_name] = dM.Region()

2083 ts[ts_name].points['i'] = self.occ.addPoint((iH[0] + iL[0]) / 2, (iH[1] + iL[1]) / 2, 0)

2084 ts[ts_name].points['o'] = self.occ.addPoint((oH[0] + oL[0]) / 2, (oH[1] + oL[1]) / 2, 0)

2085 ts[ts_name].lines[str(ht_last) + '_' + str(ht_next_first)] =\

2086 self.occ.addLine(ts[ts_name].points['i'], ts[ts_name].points['o'])

2087 # Get insulation thickness

2088 th[ts_name] = Func.sig_dig((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)

2089 # if 'cl' + str(coil_nr) + th_name not in th:

2090 # th['cl' + str(coil_nr) + th_name] = float((Func.points_distance(iH, iL) + Func.points_distance(oH, oL)) / 2)

2091 # Create mid-turn thin shells

2092 mid_turn_ts[str(blk_order.block)] = dM.Region()

2093 ts = mid_turn_ts[str(blk_order.block)]

2094 for nr_ht, ht in enumerate(ht_list[:-1]):

2095 line_name = ht + '_' + ht_list[nr_ht + 1]

2096 current_ht = block.half_turns[int(ht)].corners.bare

2097 next_ht = block.half_turns[int(ht_list[nr_ht + 1])].corners.bare

2098 mid_inner = [(current_ht.iH.x + next_ht.iL.x) / 2, (current_ht.iH.y + next_ht.iL.y) / 2]

2099 mid_outer = [(current_ht.oH.x + next_ht.oL.x) / 2, (current_ht.oH.y + next_ht.oL.y) / 2]

2100 mid_length = Func.points_distance(mid_inner, mid_outer)

2101 mid_line = Func.line_through_two_points(mid_inner, mid_outer)

2102 points = {'inner': list, 'outer': list}

2103 for case, current_h, current_l, next_h, next_l, mid_point in zip(

2104 ['inner', 'outer'], [current_ht.iH, current_ht.oH], [current_ht.iL, current_ht.oL],

2105 [next_ht.iH, next_ht.oH], [next_ht.iL, next_ht.oL], [mid_outer, mid_inner]):

2106 current_line = Func.line_through_two_points([current_h.x, current_h.y], [current_l.x, current_l.y])

2107 next_line = Func.line_through_two_points([next_h.x, next_h.y], [next_l.x, next_l.y])

2108 current_intersect = Func.intersection_between_two_lines(mid_line, current_line)

2109 next_intersect = Func.intersection_between_two_lines(mid_line, next_line)

2110 points[case] = current_intersect if Func.points_distance(

2111 current_intersect, mid_point) < mid_length else next_intersect

2112 ts.points[line_name + '_i'] = self.occ.addPoint(points['inner'][0], points['inner'][1], 0)

2113 ts.points[line_name + '_o'] = self.occ.addPoint(points['outer'][0], points['outer'][1], 0)

2114 ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])

2115

2116 # Create mid-layer thin shells

2117 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]

2118 for ts_name, ts in mid_layer_ts.items():

2119 # Order mid-layer thin shell points according to their angle with respect to the x-axis to generate lines

2120 blk1, blk2 = ts_name.split('_')

2121 max_angle = max(ts.point_angles.values())

2122 max_diff = max_angle - min(ts.point_angles.values())

2123 if max_diff > np.pi:

2124 for pnt_name, angle in ts.point_angles.items():

2125 if angle < max_diff / 2:

2126 ts.point_angles[pnt_name] = angle + max_angle

2127 ordered_pnts = [[pnt_name, ts.point_angles[pnt_name], pnt] for pnt_name, pnt in ts.mid_layers.points.items()]

2128 ordered_pnts.sort(key=lambda x: x[1])

2129 for nr, pnt in enumerate(ordered_pnts[:-1]):

2130 pnt_current = pnt[0]

2131 pnt_next = ordered_pnts[nr + 1][0]

2132 if ((pnt_current[-1] == 'l' and pnt_next[-1] == 'h' and ts_name not in block_coil_mid_pole_list) or # cos-theta

2133 (ts_name in block_coil_mid_pole_list and

2134 ((pnt_current[-1] == pnt_next[-1] == 'h' and block_coil_mid_pole_list.index(ts_name) == 0) or # assumes a dipole block-coil

2135 (pnt_current[-1] == pnt_next[-1] == 'l' and block_coil_mid_pole_list.index(ts_name) == 1) or # assumes a dipole block-coil

2136 (pnt_current[:-1] == pnt_next[:-1])))):

2137 if pnt_current[:-1] == pnt_next[:-1]:

2138 relevant_blk = blk2 if int(pnt_current[:-1]) in ts.half_turn_lists[blk1] else blk1

2139 if nr > 0:

2140 iter_nr = nr - 1

2141 while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr -= 1

2142 line_name = ordered_pnts[iter_nr][0][:-1] + '_' + pnt_current[:-1]

2143 else:

2144 if len(ordered_pnts) == 2: # todo: get right ht from relevant_blk for 1-ht blocks

2145 line_name = pnt_current[:-1] + '_' + str(ts.half_turn_lists[relevant_blk][0])

2146 else:

2147 iter_nr = nr + 1

2148 while int(ordered_pnts[iter_nr][0][:-1]) not in ts.half_turn_lists[relevant_blk]: iter_nr += 1

2149 line_name = pnt_current[:-1] + '_' + ordered_pnts[iter_nr][0][:-1]

2150 else:

2151 line_name = pnt_current[:-1] + '_' + pnt_next[:-1]

2152

2153 # TODO look into why this exception handling is needed for SMC magnet with ESC coils - the following meshing stage does not work

2154 try:

2155 tag = self.occ.addLine(pnt[2], ordered_pnts[nr + 1][2])

2156 ts.mid_layers.lines[line_name] = tag

2157 except Exception as e:

2158 ts.mid_layers.lines[line_name] = tag # this will be the last tag, i.e. from previously created line

2159 x1, y1, z1 = gmsh.model.occ.getBoundingBox(1, pnt[2])[:3]

2160 x2, y2, z2 = gmsh.model.occ.getBoundingBox(1, ordered_pnts[nr + 1][2])[:3]

2161 distance = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2 + (z2 - z1) ** 2)

2162 logger.info(f"{e} {line_name} between point tags {pnt[2], ordered_pnts[nr + 1][2]} and distance between points {distance}")

2163

2164 if ts_name in mid_layer_ts_aux:

2165 aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[0]

2166 other_pnt = ordered_pnts[0 if aux_pnt[-1] == 'l' else -1]

2167 other_pnt_coord = gmsh.model.getValue(0, other_pnt[2], [])[:2] # needs to be a new point

2168 mid_layer_ts_aux[ts_name].points[other_pnt[0]] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)

2169 line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]

2170 try:

2171 mid_layer_ts_aux[ts_name].lines[line_name] = \

2172 self.occ.addLine(mid_layer_ts_aux[ts_name].points[aux_pnt], mid_layer_ts_aux[ts_name].points[other_pnt[0]])

2173 except:

2174 mid_layer_ts_aux[ts_name].lines.pop(line_name)

2175

2176 # Create wedge-to-block and block-to-wedge lines

2177 for wdg_ts in [self.md.geometries.thin_shells.mid_layers_wdg_to_ht, self.md.geometries.thin_shells.mid_layers_ht_to_wdg]:

2178 for ts_name, ts in wdg_ts.items():

2179 pnt_list = list(ts.points.keys())

2180 for nr, pnt in enumerate(pnt_list[:-1]):

2181 if pnt[-1] == 'l' and pnt_list[nr + 1][-1] == 'h':

2182 ts.lines[pnt[:-1] + '_' + pnt_list[nr + 1][:-1]] = self.occ.addLine(ts.points[pnt], ts.points[pnt_list[nr + 1]])

2183 if ts_name in mid_layer_ts_aux:

2184 aux_pnt = list(mid_layer_ts_aux[ts_name].points.keys())[

2185 1 if list(mid_layer_ts_aux[ts_name].points.keys()).index('center') == 0 else 0]

2186 other_pnt = pnt_list[0 if aux_pnt[-1] == 'l' else -1]

2187 other_pnt_coord = gmsh.model.getValue(0, ts.points[other_pnt], [])[:2] # needs to be a new point

2188 mid_layer_ts_aux[ts_name].points[other_pnt] = self.occ.addPoint(other_pnt_coord[0], other_pnt_coord[1], 0)

2189 line_name = list(mid_layer_ts_aux[ts_name].lines.keys())[0]

2190 mid_layer_ts_aux[ts_name].lines[line_name] = self.occ.addCircleArc(

2191 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])

2192

2193 # Create wedge-to-wedge lines

2194 for ts_nr, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():

2195 ts.lines[ts_nr] = self.occ.addCircleArc(ts.points['beg'], ts.points['center'], ts.points[list(ts.points.keys())[-1]])

2196

2197 # Create mid wedge-turn lines

2198 mid_turn_ts = self.md.geometries.thin_shells.mid_wedge_turn

2199 for ts_nr, ts in mid_turn_ts.items():

2200 line_name = list(ts.points.keys())[0][:-2]

2201 ts.lines[line_name] = self.occ.addLine(ts.points[line_name + '_i'], ts.points[line_name + '_o'])

2202

2203 # Get insulation thickness

2204 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():

2205 geom_coil = self.geom.coil.coils[coil_nr]

2206 # Get block-coil mid-pole thickness

2207 if coil_nr in self.block_coil_mid_pole_blks:

2208 for blk_orders in self.block_coil_mid_pole_blks[coil_nr]:

2209 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

2210 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

2211 ins_th.mid_layer[str(blk_orders[0].block) + '_' + str(blk_orders[1].block)] = Func.sig_dig(abs(block_y - block_next_y))

2212 # 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

2213 for layer_nr, layer in coil.layers.items():

2214 for blk_order in layer:

2215 for ts_name, ts in mid_layer_ts.items():

2216 blk1, blk2 = ts_name.split('_')

2217 if blk1 == str(blk_order.block) and ts_name not in block_coil_mid_pole_list:

2218 block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]

2219 if layer_nr < len(coil.layers):

2220 for blk_order_next in coil.layers[layer_nr + 1]:

2221 if blk_order_next.block == int(blk2):

2222 block_next = geom_coil.poles[blk_order_next.pole].layers[layer_nr + 1].windings[blk_order_next.winding].blocks[int(blk2)]

2223 break

2224 else:

2225 for blk_order_next in self.md.geometries.coil.anticlockwise_order.coils[coil_nr + 1].layers[1]:

2226 if blk_order_next.block == int(blk2):

2227 block_next = self.geom.coil.coils[coil_nr + 1].poles[blk_order_next.pole].layers[1].windings[blk_order_next.winding].blocks[int(blk2)]

2228 break

2229 distances = []

2230 lines = list(ts.mid_layers.lines.keys())

2231 for line_name in [lines[0], lines[-1]]:

2232 ht_1, ht_2 = int(line_name[:line_name.index('_')]), int(line_name[line_name.index('_') + 1:])

2233 ht_char = {'low_p': ht_1, 'high_p': ht_2,

2234 'current': ht_1 if ht_1 in ts.half_turn_lists[blk1] else ht_2,

2235 'next': ht_2 if ht_1 in ts.half_turn_lists[blk1] else ht_1}

2236 hts = {'current': block.half_turns[ht_char['current']].corners.bare,

2237 'next': block_next.half_turns[ht_char['next']].corners.bare}

2238 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],

2239 'high_p': [hts['current'].oH, hts['current'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['next'].iH, hts['next'].oH],

2240 'low_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['low_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH],

2241 'high_p_opp': [hts['next'].iL, hts['next'].iH] if ht_char['high_p'] == ht_char['current'] else [hts['current'].oL, hts['current'].oH]}

2242 low_line = Func.line_through_two_points([hts_p['low_p'][0].x, hts_p['low_p'][0].y],

2243 [hts_p['low_p'][1].x, hts_p['low_p'][1].y])

2244 high_line = Func.line_through_two_points([hts_p['high_p'][0].x, hts_p['high_p'][0].y],

2245 [hts_p['high_p'][1].x, hts_p['high_p'][1].y])

2246 distances.extend([Func.points_distance([hts_p['low_p'][0].x, hts_p['low_p'][0].y], Func.intersection_between_two_lines(

2247 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]))),

2248 Func.points_distance([hts_p['high_p'][0].x, hts_p['high_p'][0].y], Func.intersection_between_two_lines(

2249 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])))])

2250 ins_th.mid_layer[ts_name] = Func.sig_dig(min(distances))

2251 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]):

2252 for ts_name, ts in wdg_ts.items():

2253 wdg, blk = ts_name.split('_')

2254 if blk == str(blk_order.block):

2255 block = geom_coil.poles[blk_order.pole].layers[layer_nr].windings[blk_order.winding].blocks[blk_order.block]

2256 wedge = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + (1 if ts_type == 1 else -1)].wedges[int(wdg[1:])]

2257 pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]

2258 pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]

2259 pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]

2260 pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]

2261 low_line = Func.line_through_two_points(pnt_il, pnt_ol)

2262 high_line = Func.line_through_two_points(pnt_ih, pnt_oh)

2263 el1_l, el2_l = list(ts.lines.keys())[0].split('_')

2264 ht_l = block.half_turns[int(el1_l) if el2_l == wdg else int(el2_l)].corners.bare

2265 el1_h, el2_h = list(ts.lines.keys())[-1].split('_')

2266 ht_h = block.half_turns[int(el1_h) if el2_h == wdg else int(el2_h)].corners.bare

2267 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\

2268 else Func.line_through_two_points([ht_l.oL.x, ht_l.oL.y], [ht_l.oH.x, ht_l.oH.y])

2269 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 \

2270 else Func.line_through_two_points([ht_h.oL.x, ht_h.oL.y], [ht_h.oH.x, ht_h.oH.y])

2271 ins_th.mid_layer[ts_name] = Func.sig_dig(

2272 (Func.points_distance(pnt_ol if ts_type == 0 else pnt_il, Func.intersection_between_two_lines(low_line, opp_line_l)) +

2273 Func.points_distance(pnt_oh if ts_type == 0 else pnt_ih, Func.intersection_between_two_lines(high_line, opp_line_h))) / 2)

2274

2275 for coil_nr, coil in self.md.geometries.wedges.coils.items():

2276 # 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

2277 for layer_nr, layer in coil.layers.items():

2278 for wedge_nr, wedge in layer.wedges.items():

2279 for ts_name, ts in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg.items():

2280 wdg1, wdg2 = ts_name[1:ts_name.index('_')], ts_name[ts_name.index('_') + 2:]

2281 if wdg1 == str(wedge_nr):

2282 wedge_next = self.md.geometries.wedges.coils[coil_nr].layers[layer_nr + 1].wedges[int(wdg2)]

2283 # pnt_il_next = gmsh.model.getValue(0, wedge_next.points['il'], [])[:2]

2284 # pnt_ih_next = gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2]

2285 pnt_il = gmsh.model.getValue(0, wedge.points['il'], [])[:2]

2286 pnt_ol = gmsh.model.getValue(0, wedge.points['ol'], [])[:2]

2287 pnt_ih = gmsh.model.getValue(0, wedge.points['ih'], [])[:2]

2288 pnt_oh = gmsh.model.getValue(0, wedge.points['oh'], [])[:2]

2289 low_line = Func.line_through_two_points(pnt_il, pnt_ol)

2290 high_line = Func.line_through_two_points(pnt_ih, pnt_oh)

2291 opp_line = Func.line_through_two_points(gmsh.model.getValue(0, wedge_next.points['il'], [])[:2],

2292 gmsh.model.getValue(0, wedge_next.points['ih'], [])[:2])

2293 ins_th.mid_layer[ts_name] = Func.sig_dig(

2294 (Func.points_distance(pnt_ol, Func.intersection_between_two_lines(low_line, opp_line)) +

2295 Func.points_distance(pnt_oh, Func.intersection_between_two_lines(high_line, opp_line))) / 2)

2296

2297 def buildDomains(self, run_type, symmetry):

2298 """

2299 Generates plane surfaces from the curve loops

2300 """

2301 iron = self.geom.iron

2302 gm = self.md.geometries

2303 geometry_setting = self.data.magnet.geometry.electromagnetics if run_type == 'EM' \

2304 else self.data.magnet.geometry.thermal

2305

2306 with_wedges = geometry_setting.with_wedges

2307

2308 inv_nc = {v: k for k, v in self.nc.items()} #invert naming convention

2309 for a in geometry_setting.areas: # a in ['iron_yoke', 'collar', ...]:

2310 for quadrant, qq in getattr(gm, a).quadrants.items():

2311 for area_name, area in qq.areas.items():

2312 identifier = next((k for k in self.inv_nc.keys() if (k in area_name[2:])), None)#re.sub(r'\d+', '',area_name[2:])

2313 if a == inv_nc.get(identifier, None): # ensure it is part of the iron yoke or collar (iron, collar)

2314 build = True

2315 loops = [area.loop]

2316 for hole_key, hole in iron.hyper_holes.items():

2317 if area_name == hole.areas[1]:

2318 loops.append(qq.areas[hole.areas[0]].loop)

2319 elif area_name == hole.areas[0]: #skip holes

2320 area.surface = self.occ.addPlaneSurface(loops) # also build the holes. An existing curveloop without area is very annoying

2321 build = False

2322 if build:

2323 area.surface = self.occ.addPlaneSurface(loops)

2324 getattr(self.md.domains.groups_entities, a)[iron.hyper_areas[area_name].material].append(area.surface) ## save the material

2325

2326 # Build coil domains

2327 for coil_nr, coil in gm.coil.coils.items():

2328 for pole_nr, pole in coil.poles.items():

2329 for layer_nr, layer in pole.layers.items():

2330 for winding_nr, winding in layer.windings.items():

2331 for block_key, block in winding.blocks.items():

2332 for area_name, area in block.half_turns.areas.items():

2333 area.surface = self.occ.addPlaneSurface([area.loop])

2334

2335 # Build wedges domains

2336 if with_wedges:

2337 for coil_nr, coil in gm.wedges.coils.items():

2338 for layer_nr, layer in coil.layers.items():

2339 for wedge_nr, wedge in layer.wedges.items():

2340 wedge.areas[str(wedge_nr)].surface = self.occ.addPlaneSurface([wedge.areas[str(wedge_nr)].loop])

2341

2342 # Build insulation domains

2343 if run_type == 'TH' and not geometry_setting.use_TSA:

2344 for coil_nr, coil in gm.insulation.coils.items():

2345 for group_nr, group in coil.group.items():

2346 holes = []

2347 for blk in group.blocks:

2348 holes.extend([ht.loop for ht_nr, ht in gm.coil.coils[

2349 coil_nr].poles[blk[0]].layers[blk[1]].windings[blk[2]].blocks[blk[3]].half_turns.areas.items()])

2350 for wdg in group.wedges:

2351 holes.extend([wedge.loop for wedge_nr, wedge in gm.wedges.coils[

2352 coil_nr].layers[wdg[0]].wedges[wdg[1]].areas.items()])

2353 if len(group.ins.areas) == 1:

2354 for area_name, area in group.ins.areas.items():

2355 area.surface = self.occ.addPlaneSurface([area.loop] + holes)

2356 else:

2357 for area_name, area in group.ins.areas.items():

2358 if area_name.isdigit():

2359 area.surface = self.occ.addPlaneSurface([area.loop] + holes + [group.ins.areas['inner_loop'].loop])

2360

2361 # Create and build air far field

2362 if run_type == 'EM':

2363 if 'iron_yoke' in geometry_setting.areas:

2364 for i in iron.key_points:

2365 gm.iron_yoke.max_radius = max(gm.iron_yoke.max_radius, max(iron.key_points[i].x, iron.key_points[i].y)) # this also contains other regions, e.g. collar but this has no effect

2366 greatest_radius = gm.iron_yoke.max_radius

2367 else: # no iron yoke data available

2368 for coil_nr, coil in self.geom.coil.coils.items():

2369 for pole_nr, pole in coil.poles.items():

2370 first_winding = list(pole.layers[len(pole.layers)].windings.keys())[0]

2371 first_block = list(pole.layers[len(pole.layers)].windings[first_winding].blocks)[0]

2372 gm.coil.max_radius = max(abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.x),

2373 abs(pole.layers[len(pole.layers)].windings[first_winding].blocks[first_block].block_corners.oL.y),

2374 gm.coil.max_radius)

2375 greatest_radius = gm.coil.max_radius

2376 radius_in = greatest_radius * (2.5 if 'iron_yoke' in geometry_setting.areas else 6)

2377 radius_out = greatest_radius * (3.2 if 'iron_yoke' in geometry_setting.areas else 8)

2378 air_inf_center_x, air_inf_center_y = 0, 0

2379 for coil_nr, coil in self.md.geometries.coil.coils.items():

2380 air_inf_center_x += coil.bore_center.x

2381 air_inf_center_y += coil.bore_center.y

2382 gm.air.points['bore_center' + str(coil_nr)] = self.occ.addPoint(coil.bore_center.x, coil.bore_center.y, 0.)

2383 air_inf_center = [air_inf_center_x / len(self.md.geometries.coil.coils), air_inf_center_y / len(self.md.geometries.coil.coils)]

2384 if symmetry == 'none':

2385 gm.air_inf.lines['inner'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_in)

2386 gm.air_inf.lines['outer'] = self.occ.addCircle(air_inf_center[0], air_inf_center[1], 0., radius_out)

2387 gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['inner']]))

2388 gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['outer']]))

2389 gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop, gm.air_inf.areas['inner'].loop])

2390 else:

2391 pnt1 = [1, 0] if symmetry in ['xy', 'x'] else [0, -1]

2392 pnt2 = [0, 1] if symmetry in ['xy', 'y'] else [-1, 0]

2393 gm.air.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_in, pnt1[1] * radius_in, 0)

2394 gm.air.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_in, pnt2[1] * radius_in, 0)

2395 gm.air_inf.points['pnt1'] = self.occ.addPoint(pnt1[0] * radius_out, pnt1[1] * radius_out, 0)

2396 gm.air_inf.points['pnt2'] = self.occ.addPoint(pnt2[0] * radius_out, pnt2[1] * radius_out, 0)

2397 gm.air.lines['ln1'] = self.occ.addLine(gm.air.points['pnt1'], gm.air_inf.points['pnt1'])

2398 gm.air.lines['ln2'] = self.occ.addLine(gm.air.points['pnt2'], gm.air_inf.points['pnt2'])

2399 if not self.data.magnet.geometry.electromagnetics.with_iron_yoke:

2400 gm.air_inf.points['center'] = self.occ.addPoint(0, 0, 0)

2401 gm.air_inf.lines['inner'] = self.occ.addCircleArc(gm.air.points['pnt2'], gm.air_inf.points['center'], gm.air.points['pnt1'])

2402 gm.air_inf.lines['outer'] = self.occ.addCircleArc(gm.air_inf.points['pnt2'], gm.air_inf.points['center'], gm.air_inf.points['pnt1'])

2403

2404 if symmetry in ['xy', 'x']:

2405 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

2406 gm.iron.quadrants[1].points[self.symmetric_bnds['x_p']['pnts'][-1][0]], gm.air.points['pnt1'])

2407 self.symmetric_loop_lines['x'].append(gm.air.lines['x_p'])

2408 else: # y

2409 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'])

2410 self.symmetric_loop_lines['y'].append(gm.air.lines['y_n'])

2411 if symmetry in ['xy', 'y']:

2412 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'])

2413 self.symmetric_loop_lines['y'].insert(0, gm.air.lines['y_p'])

2414 else: # x

2415 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

2416 gm.iron.quadrants[2].points[self.symmetric_bnds['x_n']['pnts'][-1][0]], gm.air.points['pnt2'])

2417 self.symmetric_loop_lines['x'].insert(0, gm.air.lines['x_n'])

2418

2419 inner_lines = self.symmetric_loop_lines['x'] + [gm.air_inf.lines['inner']] + self.symmetric_loop_lines['y']\

2420 if symmetry == 'xy' else self.symmetric_loop_lines[symmetry] + [gm.air_inf.lines['inner']]

2421 gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop(inner_lines))

2422 gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop(

2423 [gm.air.lines['ln1'], gm.air_inf.lines['outer'], gm.air.lines['ln2'], gm.air_inf.lines['inner']]))

2424 gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop])

2425 # self.md.domains.groups_entities.air_inf = [gm.air_inf.areas['outer'].surface]

2426 gm.air_inf.areas['inner'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['inner'].loop])

2427

2428 self.occ.synchronize()

2429 #self.gu.launch_interactive_GUI()

2430

2431 def fragment(self):

2432 """

2433 Fragment and group air domains

2434 """

2435 # Collect surfaces to be subtracted by background air

2436 holes = []

2437

2438 # iron yoke and collar

2439 group_keys = self.nc.keys()

2440

2441 for key in group_keys:

2442 group = getattr(self.md.domains.groups_entities, key)

2443 for _, surfaces in group.items():

2444 holes.extend([(2, s) for s in surfaces])

2445

2446 # Coils

2447 for coil_nr, coil in self.md.geometries.coil.coils.items():

2448 for pole_nr, pole in coil.poles.items():

2449 for layer_nr, layer in pole.layers.items():

2450 for winding_nr, winding in layer.windings.items():

2451 for block_key, block in winding.blocks.items():

2452 for area_name, area in block.half_turns.areas.items():

2453 holes.append((2, area.surface))

2454 # Wedges

2455 for coil_nr, coil in self.md.geometries.wedges.coils.items():

2456 for layer_nr, layer in coil.layers.items():

2457 for wedge_nr, wedge in layer.wedges.items():

2458 for area_name, area in wedge.areas.items():

2459 holes.append((2, area.surface))

2460 # Insulation

2461 # if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:

2462 # for coil_nr, coil in self.md.geometries.insulation.coils.items():

2463 # for group_nr, group in coil.group.items():

2464 # for area_name, area in group.ins.areas.items():

2465 # holes.append((2, area.surface))

2466

2467 # Fragment

2468 fragmented = self.occ.fragment([(2, self.md.geometries.air_inf.areas['inner'].surface)], holes)[1]

2469 self.occ.synchronize()

2470

2471 self.md.domains.groups_entities.air = []

2472 existing_domains = [e[0][1] for e in fragmented[1:]]

2473 for e in fragmented[0]:

2474 if e[1] not in existing_domains:

2475 self.md.domains.groups_entities.air.append(e[1])

2476

2477 def updateTags(self, run_type, symmetry):

2478 # Update half turn line tags

2479 for coil_nr, coil in self.md.geometries.coil.coils.items():

2480 for pole_nr, pole in coil.poles.items():

2481 for layer_nr, layer in pole.layers.items():

2482 for winding_nr, winding in layer.windings.items():

2483 for block_key, block in winding.blocks.items():

2484 hts = block.half_turns

2485 # Get half turn ID numbers

2486 area_list = list(hts.areas.keys())

2487 for nr, ht_nr in enumerate(area_list):

2488 first_tag = int(min(gmsh.model.getAdjacencies(2, hts.areas[ht_nr].surface)[1]))

2489 hts.lines[ht_nr + 'i'] = first_tag

2490 hts.lines[ht_nr + 'l'] = first_tag + 1

2491 hts.lines[ht_nr + 'o'] = first_tag + 2

2492 hts.lines[ht_nr + 'h'] = first_tag + 3

2493

2494 # Update collar tags

2495 if run_type == "TH" and 'collar' in self.data.magnet.geometry.thermal.areas:

2496 for quad, old_tags in self.md.geometries.collar.quadrants.items():

2497 self.md.geometries.collar.inner_boundary_tags[quad] = [] # reset the inner boundary tags

2498 new_tags = []

2499 for name, area in self.md.geometries.collar.quadrants[quad].areas.items(): # arcol contains the boundaries of the holes too

2500 if not re.match(r"^ar.h", name):

2501 # the issue is that you don't know which line is which, e.g. which is the inner collar line

2502 new_tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])

2503

2504 for k, name in enumerate(self.md.geometries.collar.quadrants[quad].lines.keys()):

2505 self.md.geometries.collar.quadrants[quad].lines[name] = new_tags[k]

2506

2507 # Update inner collar tags

2508 collar_lines = [self.md.geometries.collar.quadrants[quad].lines[name] for name in self.md.geometries.collar.quadrants[quad].lines.keys()]

2509 closest_dist = 1000.

2510 closest_lines = []

2511 max_dist = 0.0

2512 max_lines = []

2513 # We assume that the middle of the curve of the collar is the closest one to the centre (0,0)

2514 for tag in collar_lines:

2515 ##x, y, _ = gmsh.model.getValue(1, tag, [0.0]) # pick one point on the line

2516 curve_type = gmsh.model.getType(1, tag)

2517 if curve_type == 'Line': # find the middle of the line

2518 tag1, tag2 = gmsh.model.getAdjacencies(1, tag)[1]

2519 x1, y1, z1 = gmsh.model.getValue(0, tag1, [])

2520 x2, y2, z2 = gmsh.model.getValue(0, tag2, [])

2521 x = 0.5 * (x1 + x2)

2522 y = 0.5 * (y1 + y2)

2523 # take the average of the end points

2524 dist = np.sqrt(x ** 2 + y ** 2)

2525 elif curve_type == "Circle": # use any point on the circle (same distance because concentric with origin)

2526 x, y, z = gmsh.model.getValue(1, tag, [0.5])

2527 dist = np.sqrt(x ** 2 + y ** 2)

2528

2529 if dist < closest_dist+1e-10:

2530 if dist < closest_dist-1e-10: # clear if new min is found

2531 closest_dist = dist

2532 closest_lines = []

2533 closest_lines.append(tag)

2534 if dist > max_dist-1e-10:

2535 if dist > max_dist+1e-10: # clear if new max is found

2536 max_dist = dist

2537 max_lines = []

2538 max_lines.append(tag)

2539 self.md.geometries.collar.inner_boundary_tags[quad] = closest_lines

2540 self.md.geometries.collar.outer_boundary_tags[quad] = max_lines

2541

2542 # outer collar tags, does not work because geom.iron has the old tags and this is not accurate anymore

2543 """

2544 outer = [old_tags.lines[name] for name in old_tags.lines.keys() if

2545 self.geom.iron.hyper_lines[name].type == 'line']

2546 logger.info("outer tags", outer)

2547 self.md.geometries.collar.outer_boundary_tags[quad] = outer

2548 """

2549

2550 # concerning the TSL, it seems impossible to get the tags of the lines, as they are not adjacent to a surface nor (yet) grouped in a physical group

2551 # only way to ensure equal tags is by creating them in the same way as gmsh orders the lines. (e.g. all at the start or all at the end would work)

2552 # we cannot swap order... soo lets hope that just a shift in numbers is sufficient

2553 if run_type == 'TH' and self.data.magnet.geometry.thermal.use_TSA_new and self.data.magnet.mesh.thermal.collar.Enforce_TSA_mapping:

2554 shift = len(self.md.geometries.collar.inner_boundary_tags[1])*4 -4

2555

2556 if 'poles' in self.data.magnet.geometry.thermal.areas:

2557 shift -= (1*4) # we shifted too much

2558

2559 ## update TSA collar lines

2560 for _, ts in self.md.geometries.thin_shells.collar_layers.items():

2561 ts.lines['1'] += shift

2562 for _, ts in self.md.geometries.thin_shells.pole_layers.items():

2563 ts.lines['1'] += shift

2564

2565 atts = ['mid_layers_ht_to_ht', 'mid_layers_wdg_to_ht', 'mid_layers_ht_to_wdg', 'mid_layers_wdg_to_wdg', 'mid_poles', 'mid_windings', 'mid_turn_blocks' , 'mid_wedge_turn']

2566 for at in atts:

2567 for _, ts_region in getattr(self.md.geometries.thin_shells, at).items():

2568 try: ts_region = ts_region.mid_layers

2569 except AttributeError: pass

2570 for key in ts_region.lines.keys():

2571 ts_region.lines[key] += shift

2572

2573 # Update coil cooling tags

2574 ### no consistent way to get the tags of the cooling lines, so we assume that they are ordered in the same way as gmsh orders the lines

2575 if self.data.magnet.solve.thermal.collar_cooling.enabled:

2576 tags =[]

2577 ## if we want all cooling holes

2578 if self.data.magnet.solve.thermal.collar_cooling.which == 'all':

2579 for quad, region in self.md.geometries.collar.quadrants.items():

2580 for name, area in region.areas.items():

2581 if re.match(r"^ar.h", name):

2582 tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])

2583 else:

2584 nr_applied_cooling = self.data.magnet.solve.thermal.collar_cooling.which

2585 nr = 1

2586 for _, quad_data in self.md.geometries.collar.quadrants.items():

2587 for name, area in quad_data.areas.items():

2588 if re.match(r"^ar.h", name):

2589 if nr in nr_applied_cooling:

2590 tags.extend([int(k) for k in gmsh.model.getAdjacencies(2, area.surface)[1]])

2591 nr += 1

2592 self.md.geometries.collar.cooling_tags = tags

2593

2594

2595 # Update insulation line tags

2596 if run_type == 'TH' and not self.data.magnet.geometry.thermal.use_TSA:

2597 pass # todo

2598

2599 # Update wedge line tags

2600 for coil_nr, coil in self.md.geometries.wedges.coils.items():

2601 for layer_nr, layer in coil.layers.items():

2602 for wedge_nr, wedge in layer.wedges.items():

2603 lines_tags = list(gmsh.model.getAdjacencies(2, wedge.areas[str(wedge_nr)].surface)[1])

2604 lines_tags.sort(key=lambda x: x)

2605 wedge.lines['i'] = int(lines_tags[0])

2606 wedge.lines['l'] = int(lines_tags[1])

2607 wedge.lines['o'] = int(lines_tags[2])

2608 wedge.lines['h'] = int(lines_tags[3])

2609

2610 if run_type == 'EM':

2611 def _get_bnd_lines():

2612 return [pair[1] for pair in self.occ.getEntitiesInBoundingBox(corner_min[0], corner_min[1], corner_min[2],

2613 corner_max[0], corner_max[1], corner_max[2], dim=1)]

2614

2615 tol = 1e-6

2616 # Update tags of air and air_inf arcs and their points

2617 lines_tags = gmsh.model.getAdjacencies(2, self.md.geometries.air_inf.areas['outer'].surface)[1]

2618 self.md.geometries.air_inf.lines['outer'] = int(lines_tags[0 if symmetry == 'none' else 1])

2619 self.md.geometries.air_inf.lines['inner'] = int(lines_tags[1 if symmetry == 'none' else 3])

2620 if symmetry == 'none': # todo: check if this holds for symmetric models too

2621 for coil_nr, coil in self.md.geometries.coil.coils.items():

2622 self.md.geometries.air.points['bore_center' + str(coil_nr)] += 2

2623 else:

2624 pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['outer'])[1])

2625 indexes = [0, 1, 0] if 'x' in symmetry else [1, 0, 1]

2626 pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] >\

2627 gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]

2628 self.md.geometries.air_inf.points['pnt1'] = int(pnt_tags[pnts[0]])

2629 self.md.geometries.air_inf.points['pnt2'] = int(pnt_tags[pnts[1]])

2630 pnt_tags = list(gmsh.model.getAdjacencies(1, self.md.geometries.air_inf.lines['inner'])[1])

2631 pnts = [0, 1] if gmsh.model.getValue(0, pnt_tags[indexes[0]], [])[indexes[2]] > \

2632 gmsh.model.getValue(0, pnt_tags[indexes[1]], [])[indexes[2]] else [1, 0]

2633 self.md.geometries.air.points['pnt1'] = int(pnt_tags[pnts[0]])

2634 self.md.geometries.air.points['pnt2'] = int(pnt_tags[pnts[1]])

2635 for coil_nr, coil in self.md.geometries.coil.coils.items():

2636 self.md.geometries.air.points['bore_center' + str(coil_nr)] =(

2637 self.occ.getEntitiesInBoundingBox(-tol + coil.bore_center.x, -tol + coil.bore_center.y, -tol,

2638 tol + coil.bore_center.x, tol + coil.bore_center.y, tol, dim=0))[0][1]

2639

2640 # Group symmetry boundary lines per type

2641 if symmetry == 'xy':

2642 corner_min = [-tol, -tol, -tol]

2643 corner_max = [gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0] + tol, tol, tol]

2644 self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()

2645 corner_max = [tol, gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1] + tol, tol]

2646 self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()

2647 elif symmetry == 'x':

2648 x_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt1'], [])[0]

2649 corner_min = [- x_coord - tol, -tol, -tol]

2650 corner_max = [x_coord + tol, tol, tol]

2651 self.md.domains.groups_entities.symmetric_boundaries.x = _get_bnd_lines()

2652 elif symmetry == 'y':

2653 y_coord = gmsh.model.getValue(0, self.md.geometries.air_inf.points['pnt2'], [])[1]

2654 corner_min = [-tol, - y_coord - tol, -tol]

2655 corner_max = [tol, y_coord + tol, tol]

2656 self.md.domains.groups_entities.symmetric_boundaries.y = _get_bnd_lines()

2657

2658 def move_keypoints(self, keypoints, displacement, keypoint_names=None):

2659 if keypoint_names is None:

2660 keypoint_names = []

2661 for name, hole in self.geom.iron.hyper_areas.items():

2662 if not 'ch' in name: # ch -> collar hole

2663 continue

2664 line_names = hole.lines

2665 keypoint_names.append(list(set([getattr(self.geom.iron.hyper_lines[line], kp_name)

2666 for line in line_names for kp_name in ['kp1', 'kp2', 'kp3']])))

2667 if type(displacement) == list:

2668 list_displacement = displacement

2669 elif str(displacement) == "0":

2670 list_displacement = [[0.0, 0.0], [0.0, 0.0]]

2671 elif str(displacement) == "1":

2672 list_displacement = [[0.004, -0.015], [0.03, -0.025]]

2673 elif str(displacement) == "2":

2674 list_displacement = [[0.004, -0.015], [0.0, 0.0]]

2675 elif str(displacement) == "3":

2676 list_displacement = [[-0.035, 0.045], [-0.004, -0.0015]]

2677 else:

2678 raise ValueError("displacement_type not recognized")

2679 for i, hole in enumerate(keypoint_names):

2680 for name in hole:

2681 if name is None:

2682 continue

2683 keypoints[name].x += list_displacement[i][0]

2684 keypoints[name].y += list_displacement[i][1]

2685 return keypoints