Coverage for fiqus/mesh_generators/MeshMultipole.py: 86%
978 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-11-03 01:32 +0000
1import os
2import gmsh
3import copy
4import numpy as np
5import json
6import logging
8from fiqus.utils.Utils import GmshUtils
9from fiqus.utils.Utils import FilesAndFolders as Util
10from fiqus.utils.Utils import GeometricFunctions as Func
11from fiqus.data import DataFiQuS as dF
12from fiqus.data import DataMultipole as dM
13from fiqus.data import RegionsModelFiQuS as rM
15logger = logging.getLogger(__name__)
17class Mesh:
18 def __init__(self, data: dF.FDM() = None, mesh_folder: str = None,
19 verbose: bool = False):
20 """
21 Class to generate mesh
22 :param data: FiQuS data model
23 :param verbose: If True more information is printed in python console.
24 """
25 self.data: dF.FDM() = data
26 self.mesh_folder = mesh_folder
27 self.verbose: bool = verbose
29 self.md = dM.MultipoleData()
30 self.rc = dM.MultipoleRegionCoordinate()
31 self.rm = rM.RegionsModel()
32 self.strands = None
34 self.gu = GmshUtils(self.mesh_folder, self.verbose)
35 self.gu.initialize(verbosity_Gmsh=self.data.run.verbosity_Gmsh)
36 self.occ = gmsh.model.occ
37 self.mesh = gmsh.model.mesh
39 self.brep_iron_curves = {1: set(), 2: set(), 3: set(), 4: set()}
40 self.mesh_parameters = dict.fromkeys(['SJ', 'SICN', 'SIGE', 'Gamma', 'nodes'])
41 self.geom_files = os.path.join(os.path.dirname(self.mesh_folder), self.data.general.magnet_name)
42 self.model_file = os.path.join(self.mesh_folder, self.data.general.magnet_name)
44 # Insulation sequence involving cable insulation only (turn-to-turn, outlying conductor edge)
45 self.ins_type_cond = {}
46 # Insulation sequence involving quench heaters (outlying or mid-layer/pole)
47 qh_keys = {key: {} for key in range(1, self.data.quench_protection.quench_heaters.N_strips + 1)}
48 self.ins_type_qh = {'internal_double': {}, 'internal': copy.deepcopy(qh_keys), 'external': copy.deepcopy(qh_keys)}
49 # Insulation sequence between blocks (layer-to-layer, pole-to-pole)
50 self.ins_type = {'mid_pole': {}, 'mid_winding': {}, 'mid_layer': {}, 'aux': {}}
51 self.qh_data, self.wedge_cond = {}, {}
53 self.colors = {'wedges': [86, 180, 233], # sky blue
54 'insul': [119, 136, 153], # light slate grey
55 'half_turns_pos': [213, 94, 0], # vermilion
56 'half_turns_neg': [255, 136, 42], # light vermilion
57 'air': [240, 228, 66], # yellow
58 'air_inf': [220, 208, 46], # dark yellow
59 # yoke
60 'BHiron1': [0, 114, 178], # blue
61 'BHiron2': [0, 158, 115], # bluish green
62 'BHiron4': [86, 180, 233], # sky blue
63 # key
64 'BHiron3': [220, 208, 46], # dark yellow
65 # [230, 159, 0], # orange
66 # collar
67 'BHiron5': [204, 121, 167], # hopbush
68 'BHiron6': [0, 114, 178], # blue
69 'BHiron7': [204, 121, 167]} # reddish purple
71 def clear(self):
72 self.md = dM.MultipoleData()
73 self.rc = dM.MultipoleRegionCoordinate()
74 self.rm = rM.RegionsModel()
75 gmsh.clear()
77 def ending_step(self, gui: bool = False):
78 if gui:
79 self.gu.launch_interactive_GUI()
80 else:
81 gmsh.clear()
82 gmsh.finalize()
84 def loadStrandPositions(self, run_type):
85 self.strands = json.load(open(f"{self.geom_files}_{run_type}.strs"))
87 def loadAuxiliaryFile(self, run_type):
88 self.md = Util.read_data_from_yaml(f"{self.geom_files}_{run_type}.aux", dM.MultipoleData)
90 def updateAuxiliaryFile(self, run_type):
91 Util.write_data_to_yaml(f'{self.model_file}_{run_type}.aux', self.md.model_dump())
93 def saveMeshFile(self, run_type):
94 gmsh.write(f'{self.model_file}_{run_type}.msh')
96 def saveRegionFile(self, run_type):
97 Util.write_data_to_yaml(f'{self.model_file}_{run_type}.reg', self.rm.model_dump())
99 def saveRegionCoordinateFile(self, run_type):
100 Util.write_data_to_yaml(f'{self.model_file}_{run_type}.reco', self.rc.model_dump())
102 def getIronCurvesTags(self):
103 for quadrant, qq in self.md.geometries.iron.quadrants.items():
104 for _, area in qq.areas.items():
105 if area.surface:
106 self.brep_iron_curves[quadrant] |= set(gmsh.model.getAdjacencies(2, area.surface)[1])
108 def defineMesh(self, geometry, mesh, run_type):
109 thresholds = []
110 self.occ.synchronize()
112 if mesh.conductors.field.enabled:
113 distance_conductors = self.mesh.field.add("Distance")
114 self.mesh.field.setNumbers(distance_conductors, "CurvesList",
115 [line for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items() for layer_nr, layer in coil.layers.items()
116 for _, block_order in enumerate(layer) for _, line in self.md.geometries.coil.coils[coil_nr].poles[block_order.pole].layers[
117 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns.lines.items()]
118 )
119 self.mesh.field.setNumber(distance_conductors, "Sampling", 100)
121 threshold_conductors = self.mesh.field.add("Threshold")
122 self.mesh.field.setNumber(threshold_conductors, "InField", distance_conductors)
123 self.mesh.field.setNumber(threshold_conductors, "SizeMin", mesh.conductors.field.SizeMin)
124 self.mesh.field.setNumber(threshold_conductors, "SizeMax", mesh.conductors.field.SizeMax)
125 self.mesh.field.setNumber(threshold_conductors, "DistMin", mesh.conductors.field.DistMin)
126 self.mesh.field.setNumber(threshold_conductors, "DistMax", mesh.conductors.field.DistMax)
127 self.mesh.field.setNumber(threshold_conductors, "StopAtDistMax", 1)
128 thresholds.append(threshold_conductors)
130 if mesh.wedges.field.enabled:
131 distance_wedges = self.mesh.field.add("Distance")
132 self.mesh.field.setNumbers(distance_wedges, "CurvesList",
133 [line for _, coil in self.md.geometries.wedges.coils.items() for _, layer in coil.layers.items() for _, wdg in layer.wedges.items() for _, line in wdg.lines.items()]
134 )
135 self.mesh.field.setNumber(distance_wedges, "Sampling", 100)
137 # raise Exception(f"cannot set threshold for wedges field: {[line for _, coil in self.md.geometries.wedges.coils.items() for _, layer in coil.layers.items() for _, wdg in layer.wedges.items() for _, line in wdg.lines.items()]}")
139 threshold_wedges = self.mesh.field.add("Threshold")
140 self.mesh.field.setNumber(threshold_wedges, "InField", distance_wedges)
141 self.mesh.field.setNumber(threshold_wedges, "SizeMin", mesh.wedges.field.SizeMin)
142 self.mesh.field.setNumber(threshold_wedges, "SizeMax", mesh.wedges.field.SizeMax)
143 self.mesh.field.setNumber(threshold_wedges, "DistMin", mesh.wedges.field.DistMin)
144 self.mesh.field.setNumber(threshold_wedges, "DistMax", mesh.wedges.field.DistMax)
145 self.mesh.field.setNumber(threshold_wedges, "StopAtDistMax", 1)
146 thresholds.append(threshold_wedges)
148 if geometry.with_iron_yoke:
149 distance_iron = self.mesh.field.add("Distance")
150 self.mesh.field.setNumbers(distance_iron, "CurvesList", [line for _, qq in self.brep_iron_curves.items() for line in qq])
151 self.mesh.field.setNumber(distance_iron, "Sampling", 100)
153 if mesh.iron_field.enabled:
155 threshold_iron = self.mesh.field.add("Threshold")
156 self.mesh.field.setNumber(threshold_iron, "InField", distance_iron)
157 self.mesh.field.setNumber(threshold_iron, "SizeMin", mesh.iron_field.SizeMin)
158 self.mesh.field.setNumber(threshold_iron, "SizeMax", mesh.iron_field.SizeMax)
159 self.mesh.field.setNumber(threshold_iron, "DistMin", mesh.iron_field.DistMin)
160 self.mesh.field.setNumber(threshold_iron, "DistMax", mesh.iron_field.DistMax)
161 # the iron field is typically rather coarse so that 'overwriting' other fields is not a problem
162 # its distMax however would need to be large for StopAtDistMax to be effective
163 #self.mesh.field.setNumber(threshold_iron, "StopAtDistMax", 1)
164 thresholds.append(threshold_iron)
166 if run_type == 'EM' and mesh.bore_field.enabled:
168 distance_bore = self.mesh.field.add("Distance")
169 self.mesh.field.setNumbers(distance_bore, "PointsList", [pnt for pnt_name, pnt in self.md.geometries.air.points.items() if 'bore' in pnt_name])
170 self.mesh.field.setNumber(distance_bore, "Sampling", 100)
172 threshold_bore = self.mesh.field.add("Threshold")
173 self.mesh.field.setNumber(threshold_bore, "InField", distance_bore)
174 self.mesh.field.setNumber(threshold_bore, "SizeMin", mesh.bore_field.SizeMin)
175 self.mesh.field.setNumber(threshold_bore, "SizeMax", mesh.bore_field.SizeMax)
176 self.mesh.field.setNumber(threshold_bore, "DistMin", mesh.bore_field.DistMin)
177 self.mesh.field.setNumber(threshold_bore, "DistMax", mesh.bore_field.DistMax)
178 self.mesh.field.setNumber(threshold_bore, "StopAtDistMax", 1)
179 thresholds.append(threshold_bore)
181 insulation_mesh_fields = []
182 if run_type == 'TH':
183 for coil_nr, coil in self.md.geometries.insulation.coils.items():
184 for _, group in coil.group.items():
185 # Areas
186 for area_name, area in group.ins.areas.items():
187 if area_name.isdigit():
188 insulation_mesh_fields.append(self.mesh.field.add("Constant"))
189 insulation_mesh_field = insulation_mesh_fields[-1]
190 self.mesh.field.setNumbers(insulation_mesh_field, "SurfacesList", [area.surface])
191 self.mesh.field.setNumber(insulation_mesh_field, "VIn", mesh.insulation.global_size)
193 background = self.mesh.field.add("Min")
194 self.mesh.field.setNumbers(background, "FieldsList", thresholds + insulation_mesh_fields)
195 self.mesh.field.setAsBackgroundMesh(background)
197 # Apply transfinite curves and potentially surfaces to conductors and wedges
198 if mesh.conductors.transfinite.enabled_for in ['curves', 'curves_and_surfaces']:
199 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
200 for layer_nr, layer in coil.layers.items():
201 for _, block_order in enumerate(layer):
202 winding = self.md.geometries.coil.coils[coil_nr].poles[block_order.pole].layers[
203 layer_nr].windings[block_order.winding]
204 cable = self.data.conductors[winding.conductor_name].cable
205 for line_key, line in winding.blocks[block_order.block].half_turns.lines.items():
206 if mesh.model_dump().get('isothermal_conductors', False):
207 elements = 1
208 elif any([i in line_key for i in ['i', 'o']]):
209 elements = max(1, round(cable.bare_cable_height_mean / mesh.conductors.transfinite.curve_target_size_height))
210 elif any([i in line_key for i in ['l', 'h']]):
211 elements = max(1, round(cable.bare_cable_width / mesh.conductors.transfinite.curve_target_size_width))
213 self.mesh.setTransfiniteCurve(line, elements)
214 if mesh.conductors.transfinite.enabled_for=='curves_and_surfaces' or mesh.model_dump().get('isothermal_conductors', False):
215 for _, area in winding.blocks[block_order.block].half_turns.areas.items():
216 self.mesh.setTransfiniteSurface(area.surface)
217 self.mesh.setRecombine(2, area.surface)
219 if 'insulation' in mesh.model_dump() and 'TSA' in mesh.model_dump()["insulation"]:
220 # Apply transfinite curves to thin shell lines
221 if geometry.use_TSA:
222 gts = self.md.geometries.thin_shells
223 conductor_target_sizes = {'width': mesh.conductors.transfinite.curve_target_size_width, 'height': mesh.conductors.transfinite.curve_target_size_height}
224 wedge_target_sizes = {'width': mesh.wedges.transfinite.curve_target_size_width, 'height': mesh.wedges.transfinite.curve_target_size_height}
225 for ts_group, side in zip([gts.mid_layers_ht_to_ht, gts.mid_layers_ht_to_wdg, gts.mid_layers_wdg_to_ht, gts.mid_layers_wdg_to_wdg,
226 gts.mid_poles, gts.mid_windings, gts.mid_turn_blocks, gts.mid_wedge_turn, gts.mid_layers_aux],
227 ['height', 'height', 'height', 'height', 'width', 'width', 'width', 'width', 'height']):
228 for ts_name, ts in ts_group.items():
229 for _, line in ts.lines.items() if isinstance(ts, dM.Region) else ts.mid_layers.lines.items():
230 if mesh.isothermal_conductors or mesh.isothermal_wedges:
231 elements = 1
232 else:
233 coords = gmsh.model.getValue(1, line, [i[0] for i in gmsh.model.getParametrizationBounds(1, line)])
234 target_size = wedge_target_sizes[side] if ts_name.count('w') == 2 else conductor_target_sizes[side]
235 elements = max(1, round(Func.points_distance(coords[:2], coords[3:-1]) / target_size))
237 # it's a wedge
238 if ts_name.count('w') == 2 and mesh.wedges.transfinite.enabled_for in ['curves', 'curves_and_surfaces']:
239 self.mesh.setTransfiniteCurve(line, elements)
240 elif ts_name.count('w') != 2 and mesh.conductors.transfinite.enabled_for in ['curves', 'curves_and_surfaces']:
241 self.mesh.setTransfiniteCurve(line, elements)
242 # COMMENTED since this overwrites also the cable transfinite meshes and in general,
243 # restricting the insulation boundaries to be transfinite seems very restrictive due to their complex geometry
244 # Apply transfinite curves to insulation boundaries
245 # else:
246 # for coil_nr, coil in self.md.geometries.insulation.coils.items():
247 # for group_nr, group in coil.group.items():
248 # cable_height = self.data.conductors[self.md.geometries.coil.coils[coil_nr].poles[
249 # group.blocks[0][0]].layers[group.blocks[0][1]].windings[group.blocks[0][2]].conductor_name].cable.bare_cable_height_mean
250 # for line in [bnd[1] for bnd in gmsh.model.getBoundary(
251 # [(2, list(group.ins.areas.values())[0].surface)], # +
252 # # [(2, ht.surface) for blk_order in group.blocks for ht in
253 # # self.md.geometries.coil.coils[coil_nr].poles[blk_order[0]].layers[blk_order[1]].windings[blk_order[2]].blocks[blk_order[3]].half_turns.areas.values()] +
254 # # [(2, self.md.geometries.wedges.coils[coil_nr].layers[wdg_order[0]].wedges[wdg_order[1]].areas[str(wdg_order[1])].surface) for wdg_order in group.wedges],
255 # combined=True, oriented=False)]:
256 # pnts = gmsh.model.getAdjacencies(1, line)[1]
257 # length = Func.points_distance(gmsh.model.getValue(0, pnts[0], []), gmsh.model.getValue(0, pnts[1], []))
258 # self.mesh.setTransfiniteCurve(line,max(1, round(length / (mesh.conductor_target_sizes.width if length > 2 * cable_height else mesh.conductor_target_sizes.height))))
260 # Apply transfinite curves to wedges
261 if mesh.wedges.transfinite.enabled_for in ['curves', 'curves_and_surfaces']:
262 if geometry.with_wedges:
263 for coil_nr, coil in self.md.geometries.wedges.coils.items():
264 for layer_nr, layer in coil.layers.items():
265 for _, wedge in layer.wedges.items():
266 pnts = gmsh.model.getAdjacencies(1, wedge.lines['i'])[1]
267 inner_height = Func.points_distance(gmsh.model.getValue(0, pnts[0], []), gmsh.model.getValue(0, pnts[1], []))
268 pnts = gmsh.model.getAdjacencies(1, wedge.lines['l'])[1]
269 width = Func.points_distance(gmsh.model.getValue(0, pnts[0], []), gmsh.model.getValue(0, pnts[1], []))
270 pnts = gmsh.model.getAdjacencies(1, wedge.lines['o'])[1]
271 outer_height = Func.points_distance(gmsh.model.getValue(0, pnts[0], []), gmsh.model.getValue(0, pnts[1], []))
272 for line_key, line in wedge.lines.items():
273 if mesh.model_dump().get('isothermal_wedges', False):
274 elements = 1
275 elif 'i' in line_key:
276 elements = max(1, round(inner_height / mesh.wedges.transfinite.curve_target_size_height))
277 elif 'o' in line_key:
278 elements = max(1, round((inner_height if mesh.wedges.transfinite.enabled_for in ['curves', 'curves_and_surfaces']
279 else outer_height) / mesh.wedges.transfinite.curve_target_size_height))
280 elif any([i in line_key for i in ['l', 'h']]):
281 elements = max(1, round(width / mesh.wedges.transfinite.curve_target_size_width))
282 if mesh.wedges.transfinite.enabled_for in ['curves', 'curves_and_surfaces']:
283 self.mesh.setTransfiniteCurve(line, elements)
284 if mesh.wedges.transfinite.enabled_for=='curves_and_surfaces' or mesh.model_dump().get('isothermal_wedges', False):
285 self.mesh.setTransfiniteSurface(list(wedge.areas.values())[0].surface)
286 self.mesh.setRecombine(2, list(wedge.areas.values())[0].surface)
288 def createPhysicalGroups(self, geometry):
289 """
290 Creates physical groups by grouping the mirrored entities according to the Roxie domains
291 """
292 offset: int = 1 if 'symmetry' in geometry.model_dump() else int(1e6)
293 pg_tag = offset
295 # Create physical groups of iron yoke regions and block insulation
296 pg = self.md.domains.physical_groups
297 for group_name, surfaces in self.md.domains.groups_entities.iron.items():
298 pg.iron.surfaces[group_name] = gmsh.model.addPhysicalGroup(2, surfaces, pg_tag)
299 gmsh.model.setPhysicalName(2, pg.iron.surfaces[group_name], group_name)
300 color = self.colors[group_name]
301 gmsh.model.setColor([(2, i) for i in surfaces], color[0], color[1], color[2])
302 pg_tag += 1
304 # Create the physical group of air infinite
305 if 'symmetry' in geometry.model_dump():
306 gmsh.model.setPhysicalName(0, gmsh.model.addPhysicalGroup(
307 0, [pnt for pnt_name, pnt in self.md.geometries.air.points.items() if 'bore_field' in pnt_name], pg_tag), 'bore_centers')
308 pg_tag += 1
309 pg.air_inf_bnd = gmsh.model.addPhysicalGroup(1, [self.md.geometries.air_inf.lines['outer']], pg_tag)
310 gmsh.model.setPhysicalName(1, pg.air_inf_bnd, 'air_inf_bnd')
311 pg_tag += 1
312 pg.air_inf = gmsh.model.addPhysicalGroup(2, [self.md.geometries.air_inf.areas['outer'].surface], pg_tag)
313 gmsh.model.setPhysicalName(2, pg.air_inf, 'air_inf')
314 gmsh.model.setColor([(2, self.md.geometries.air_inf.areas['outer'].surface)], self.colors['air_inf'][0], self.colors['air_inf'][1], self.colors['air_inf'][2])
315 pg_tag += 1
316 pg.air = gmsh.model.addPhysicalGroup(2, self.md.domains.groups_entities.air, pg_tag)
317 gmsh.model.setPhysicalName(2, pg.air, 'air')
318 gmsh.model.setColor([(2, i) for i in self.md.domains.groups_entities.air], self.colors['air'][0], self.colors['air'][1], self.colors['air'][2])
319 pg_tag += 1
321 # Create physical groups of half turns
322 lyr_list_group = []
323 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
324 lyr_list_group.extend(['cl' + str(coil_nr) + 'ly' + str(lyr) for lyr in list(coil.layers.keys())])
325 for layer_nr, layer in coil.layers.items():
326 ht_list_group = []
327 for nr, block_order in enumerate(layer):
328 wnd = self.md.geometries.coil.coils[coil_nr].poles[block_order.pole].layers[
329 layer_nr].windings[block_order.winding]
330 block = wnd.blocks[block_order.block]
331 block_nr = block_order.block
332 pg.blocks[block_nr] = dM.PoweredBlock()
333 blk_pg = pg.blocks[block_nr]
334 blk_pg.current_sign = block.current_sign
335 blk_pg.conductor = wnd.conductor_name
336 color = self.colors['half_turns_pos'] if block.current_sign > 0 else self.colors['half_turns_neg']
337 ht_list = list(block.half_turns.areas.keys())
338 ht_list_group.extend(ht_list)
339 if nr + 1 < len(layer):
340 if layer[nr + 1].pole == block_order.pole and layer[nr + 1].winding != block_order.winding:
341 ht_list_group.append('w' + str(nr))
342 # Create 2D physical groups of half turns
343 for ht_key, ht in block.half_turns.areas.items():
344 blk_pg.half_turns[int(ht_key)] = dM.PoweredGroup()
345 ht_pg = blk_pg.half_turns[int(ht_key)]
346 # Create physical group and assign name and color
347 ht_pg.tag = gmsh.model.addPhysicalGroup(2, [ht.surface], pg_tag)
348 gmsh.model.setPhysicalName(2, ht_pg.tag, ht_key)
349 gmsh.model.setColor([(2, ht.surface)], color[0], color[1], color[2])
350 pg_tag += 1
351 # Assign thin-shell group
352 # the check for reversed block coil is not tested well
353 if geometry.model_dump().get('correct_block_coil_tsa_checkered_scheme', False) and self.md.geometries.coil.coils[coil_nr].type == 'reversed-block-coil':
354 azimuthal = 'a1' if list(wnd.blocks.keys()).index(block_nr) % 2 == 0 else 'a2'
355 else:
356 azimuthal = 'a1' if lyr_list_group.index('cl' + str(coil_nr) + 'ly' + str(layer_nr)) % 2 == 0 else 'a2'
357 radial = 'r1' if ht_list_group.index(ht_key) % 2 == 0 else 'r2'
358 ht_pg.group = radial + '_' + azimuthal
359 if geometry.model_dump().get('use_TSA', False):
360 # Create 1D physical groups of thin shells
361 for ht_line_key, ht_line in block.half_turns.lines.items():
362 ht_nr = ht_line_key[:-1]
363 # Create half turn line groups
364 line_pg = gmsh.model.addPhysicalGroup(1, [ht_line], pg_tag)
365 gmsh.model.setPhysicalName(1, line_pg, ht_line_key)
366 color = [0, 0, 0] if blk_pg.half_turns[int(ht_nr)].group[:2] == 'r1' else [150, 150, 150]
367 gmsh.model.setColor([(1, ht_line)], color[0], color[1], color[2])
368 pg_tag += 1
369 # Store thin shell tags
370 blk_pg.half_turns[int(ht_nr)].lines[ht_line_key[-1]] = line_pg
372 # Create points region for projection
373 if 'use_TSA' in geometry.model_dump():
374 self.md.domains.physical_groups.half_turns_points = gmsh.model.addPhysicalGroup(
375 0, [gmsh.model.getAdjacencies(1, gmsh.model.getAdjacencies(2, ht.surface)[1][0])[1][0]
376 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items()
377 for layer_nr, layer in coil.layers.items() for block_order in layer
378 for ht_key, ht in self.md.geometries.coil.coils[coil_nr].poles[block_order.pole].layers[
379 layer_nr].windings[block_order.winding].blocks[block_order.block].half_turns.areas.items()], int(4e6))
380 gmsh.model.setPhysicalName(0, self.md.domains.physical_groups.half_turns_points, 'points')
382 # Create physical groups of insulations
383 if not geometry.model_dump().get('use_TSA', True):
384 for coil_nr, coil in self.md.geometries.insulation.coils.items():
385 for group_nr, group in coil.group.items():
386 # Areas
387 for area_name, area in group.ins.areas.items():
388 if area_name.isdigit():
389 pg.insulations.surfaces[area_name] = gmsh.model.addPhysicalGroup(2, [area.surface], pg_tag)
390 gmsh.model.setPhysicalName(2, pg.insulations.surfaces[area_name], 'insul_' + area_name)
391 gmsh.model.setColor([(2, area.surface)], self.colors['insul'][0], self.colors['insul'][1], self.colors['insul'][2])
392 pg_tag += 1
394 # Boundaries
395 area_name = list(group.ins.areas.keys())[0] # todo: test for Mono
396 pg.insulations.curves['ext' + area_name] = gmsh.model.addPhysicalGroup(
397 1, [bnd[1] for bnd in gmsh.model.getBoundary([(2, list(group.ins.areas.values())[0].surface)], combined=False, oriented=False)[:len(group.ins.lines)]], pg_tag)
398 gmsh.model.setPhysicalName(1, pg.insulations.curves['ext' + area_name], 'insul_ext' + area_name)
399 pg_tag += 1
400 # todo: NOT COMPLETED: would work if the tags were updated in the Geometry script after saving and loading brep
401 # side_lines = {'i': [], 'o': [], 'l': [], 'h': []}
402 # for line_name, line in group.ins.lines.items():
403 # side_lines[line_name[-1] if line_name[-1].isalpha() else sorted(line_name)[-1]].append(line)
404 # for side, side_line in side_lines.items():
405 # pg.insulations.curves[str(group_nr) + side] = gmsh.model.addPhysicalGroup(1, side_line, pg_tag)
406 # gmsh.model.setPhysicalName(1, pg.insulations.curves[str(group_nr) + side], str(group_nr) + side)
407 # pg_tag += 1
409 # Create physical groups of wedges
410 for coil_nr, coil in self.md.geometries.wedges.coils.items():
411 for layer_nr, layer in coil.layers.items():
412 for wedge_nr, wedge in layer.wedges.items():
413 pg.wedges[wedge_nr] = dM.WedgeGroup()
414 wedge_pg = pg.wedges[wedge_nr]
415 wedge_pg.tag = gmsh.model.addPhysicalGroup(2, [wedge.areas[str(wedge_nr)].surface], pg_tag)
416 gmsh.model.setPhysicalName(2, wedge_pg.tag, 'w' + str(wedge_nr))
417 gmsh.model.setColor([(2, wedge.areas[str(wedge_nr)].surface)],
418 self.colors['wedges'][0], self.colors['wedges'][1], self.colors['wedges'][2])
419 pg_tag += 1
420 # Assign thin-shell group
421 prev_block_hts = pg.blocks[layer.block_prev[wedge_nr]].half_turns
422 if len(list(prev_block_hts.keys())) > 1:
423 wedge_pg.group = prev_block_hts[list(prev_block_hts.keys())[-2]].group
424 else:
425 prev_group = prev_block_hts[list(prev_block_hts.keys())[0]].group
426 wedge_pg.group = ('r1' if prev_group[1] == '2' else 'r2') + prev_group[prev_group.index('_'):]
427 if geometry.model_dump().get('use_TSA', False):
428 # Create 1D physical groups of thin shells
429 for line_key, line in wedge.lines.items():
430 wedge_pg.lines[line_key] = gmsh.model.addPhysicalGroup(1, [line], pg_tag)
431 gmsh.model.setPhysicalName(1, wedge_pg.lines[line_key], 'w' + str(wedge_nr) + line_key)
432 color = [0, 0, 0] if wedge_pg.group[:2] == 'r1' else [150, 150, 150]
433 gmsh.model.setColor([(1, line)], color[0], color[1], color[2])
434 pg_tag += 1
436 # Create physical groups of thin shells
437 if geometry.model_dump().get('use_TSA', False):
438 gts = self.md.geometries.thin_shells
439 # Create physical groups of block mid-layer lines
440 block_coil_flag = False
441 for ts_name, ts in gts.mid_layers_ht_to_ht.items():
442 blk_i, blk_o = ts_name[:ts_name.index('_')], ts_name[ts_name.index('_') + 1:]
443 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
444 if (any(int(blk_i) == blk_order.block for blk_order in self.md.geometries.coil.anticlockwise_order.coils[coil_nr].layers[1]) and
445 any(int(blk_o) == blk_order.block for blk_order in self.md.geometries.coil.anticlockwise_order.coils[coil_nr].layers[1])): # block-coil mid-pole case
446 block_coil_flag = True
447 for line_name, line in ts.mid_layers.lines.items():
448 line_pg = gmsh.model.addPhysicalGroup(1, [line], pg_tag)
449 gmsh.model.setPhysicalName(1, line_pg, line_name)
450 pg_tag += 1
451 ht1, ht2 = int(line_name[:line_name.index('_')]), int(line_name[line_name.index('_') + 1:])
452 if ht1 in ts.half_turn_lists[blk_i]:
453 if block_coil_flag:
454 pg.blocks[int(blk_i)].half_turns[ht1].mid_layer_lines.inner[line_name] = line_pg
455 else:
456 pg.blocks[int(blk_i)].half_turns[ht1].mid_layer_lines.outer[line_name] = line_pg
457 pg.blocks[int(blk_o)].half_turns[ht2].mid_layer_lines.inner[line_name] = line_pg
458 else:
459 pg.blocks[int(blk_o)].half_turns[ht1].mid_layer_lines.inner[line_name] = line_pg
460 if block_coil_flag:
461 pg.blocks[int(blk_i)].half_turns[ht2].mid_layer_lines.inner[line_name] = line_pg
462 else:
463 pg.blocks[int(blk_i)].half_turns[ht2].mid_layer_lines.outer[line_name] = line_pg
465 # Create physical groups of wedge-to-block mid-layer lines
466 for ts_name, ts in gts.mid_layers_wdg_to_ht.items():
467 wdg, blk = ts_name.split('_')
468 for line_name, line in ts.lines.items():
469 line_pg = gmsh.model.addPhysicalGroup(1, [line], pg_tag)
470 gmsh.model.setPhysicalName(1, line_pg, line_name)
471 pg_tag += 1
472 pg.wedges[int(wdg[1:])].mid_layer_lines.outer[line_name] = line_pg
473 ht = line_name[:line_name.index('_')] if line_name[line_name.index('_') + 1:] == wdg else line_name[line_name.index('_') + 1:]
474 pg.blocks[int(blk)].half_turns[int(ht)].mid_layer_lines.inner[line_name] = line_pg
476 # Create physical groups of block-to-wedge mid-layer lines
477 for ts_name, ts in gts.mid_layers_ht_to_wdg.items():
478 wdg, blk = ts_name.split('_')
479 for line_name, line in ts.lines.items():
480 line_pg = gmsh.model.addPhysicalGroup(1, [line], pg_tag)
481 gmsh.model.setPhysicalName(1, line_pg, line_name)
482 pg_tag += 1
483 pg.wedges[int(wdg[1:])].mid_layer_lines.inner[line_name] = line_pg
484 ht = line_name[:line_name.index('_')] if line_name[line_name.index('_') + 1:] == wdg else line_name[line_name.index('_') + 1:]
485 pg.blocks[int(blk)].half_turns[int(ht)].mid_layer_lines.outer[line_name] = line_pg
487 # Create physical groups of wedge-to-wedge mid-layer lines
488 for ts_name, ts in gts.mid_layers_wdg_to_wdg.items():
489 wdg_i, wdg_o = ts_name[1:ts_name.index('_')], ts_name[ts_name.index('_') + 2:]
490 line_pg = gmsh.model.addPhysicalGroup(1, [ts.lines[list(ts.lines.keys())[0]]], pg_tag)
491 gmsh.model.setPhysicalName(1, line_pg, ts_name)
492 pg_tag += 1
493 pg.wedges[int(wdg_i)].mid_layer_lines.outer[ts_name] = line_pg
494 pg.wedges[int(wdg_o)].mid_layer_lines.inner[ts_name] = line_pg
496 # Create non-mid-layer thin shells
497 for ts_group_name, ts_group in zip(['mid_pole', 'mid_winding', 'mid_turn', 'mid_turn', 'aux'],
498 [gts.mid_poles, gts.mid_windings, gts.mid_turn_blocks, gts.mid_wedge_turn, gts.mid_layers_aux]):
499 for ts_name, ts in ts_group.items():
500 if '_' in ts_name:
501 el_name1, el_name2 = ts_name.split('_')
502 el_name1, el_name2 = el_name1.strip('w'), el_name2.strip('w')
503 else: # mid_turn_blocks
504 el_name1, el_name2 = ts_name, ts_name
505 for line_name, line in ts.lines.items():
506 line_pg = gmsh.model.addPhysicalGroup(1, [line], pg_tag)
507 gmsh.model.setPhysicalName(1, line_pg, line_name)
508 pg_tag += 1
509 if ts_group_name == 'aux':
510 ht1 = int(list(ts.points.keys())[0][:-1])
511 else:
512 ht1, ht2 = line_name.split('_')
513 pg.blocks[int(el_name2)].half_turns[int(ht2)].__dict__[ts_group_name + '_lines'][line_name] = line_pg
514 if ts_name[0] == 'w':
515 pg.wedges[int(el_name1)].__dict__[ts_group_name + '_lines'][line_name] = line_pg
516 else:
517 pg.blocks[int(el_name1)].half_turns[int(ht1)].__dict__[ts_group_name + '_lines'][line_name] = line_pg
519 # Create physical groups of symmetric boundaries
520 if geometry.model_dump().get('symmetry', 'none') != 'none':
521 line_tags_normal_free, line_tags_tangent_free = [], []
522 if geometry.symmetry == 'xy':
523 if len(self.md.geometries.coil.coils[1].poles) == 2:
524 line_tags_normal_free = self.md.domains.groups_entities.symmetric_boundaries.y
525 line_tags_tangent_free = self.md.domains.groups_entities.symmetric_boundaries.x
526 elif len(self.md.geometries.coil.coils[1].poles) == 4: # todo: think about other multi-pole types
527 line_tags_tangent_free = self.md.domains.groups_entities.symmetric_boundaries.x +\
528 self.md.domains.groups_entities.symmetric_boundaries.y
529 elif geometry.symmetry == 'x':
530 if 'solenoid' in self.md.geometries.coil.coils[1].type:
531 line_tags_normal_free = self.md.domains.groups_entities.symmetric_boundaries.x
532 else:
533 line_tags_tangent_free = self.md.domains.groups_entities.symmetric_boundaries.x
534 elif geometry.symmetry == 'y':
535 if len(self.md.geometries.coil.coils[1].poles) == 2:
536 line_tags_normal_free = self.md.domains.groups_entities.symmetric_boundaries.y
537 elif len(self.md.geometries.coil.coils[1].poles) == 4:
538 line_tags_tangent_free = self.md.domains.groups_entities.symmetric_boundaries.y
539 if line_tags_normal_free:
540 pg.symmetric_boundaries.normal_free = gmsh.model.addPhysicalGroup(1, line_tags_normal_free, pg_tag)
541 gmsh.model.setPhysicalName(1, pg.symmetric_boundaries.normal_free, 'normal_free_bnd')
542 pg_tag += 1
543 if line_tags_tangent_free:
544 pg.symmetric_boundaries.tangential_free = gmsh.model.addPhysicalGroup(1, line_tags_tangent_free, pg_tag)
545 gmsh.model.setPhysicalName(1, pg.symmetric_boundaries.tangential_free, 'tangent_free_bnd')
546 pg_tag += 1
548 def rearrangeThinShellsData(self):
549 pg = self.md.domains.physical_groups
550 ins_th = self.md.geometries.thin_shells.ins_thickness
551 qh = self.data.quench_protection.quench_heaters
553 def _store_qh_data(position, thin_shell, ts_tag):
554 qh_ins = self.ins_type_qh[position][qh.ids[ht_index]]
555 if thin_shell not in qh_ins: qh_ins[thin_shell] = []
556 qh_ins[thin_shell].append(ts_tag)
557 if thin_shell not in self.qh_data[qh.ids[ht_index]]:
558 self.qh_data[qh.ids[ht_index]][thin_shell] = {'conductor': blk_pg.conductor, 'ht_side': qh_side}
560 def _store_ts_tags(pg_el, geom_ts_name='', geom_ts_name2=None, ts_grp='', lines='', lines_side=None):
561 geom_ts = self.md.geometries.thin_shells.model_dump()[geom_ts_name]
562 for ln_name, ln_tag in (pg_el.model_dump()[lines][lines_side] if lines_side else pg_el.model_dump()[lines]).items():
563 for ts_name in ts_groups[ts_grp]:
564 if ts_name in geom_ts:
565 if ln_name in (geom_ts[ts_name][geom_ts_name2]['lines'] if geom_ts_name2 else geom_ts[ts_name]['lines']):
566 if ts_name not in self.ins_type[ts_grp]: self.ins_type[ts_grp][ts_name] = []
567 self.ins_type[ts_grp][ts_name].append(ln_tag)
568 break
570 # Collect thin shell tags
571 # Half turns
572 for blk_pg_nr, blk_pg in pg.blocks.items():
573 ts_groups = {'mid_wedge_turn': [ts_name for ts_name in self.md.geometries.thin_shells.mid_wedge_turn if blk_pg_nr == int(ts_name.split('_')[1])],
574 'aux': [ts_name for ts_name in ins_th.mid_layer if str(blk_pg_nr) in ts_name.split('_')],
575 'mid_winding': [ts_name for ts_name in ins_th.mid_winding if blk_pg_nr == int(ts_name.split('_')[0])],
576 'mid_pole': [ts_name for ts_name in ins_th.mid_pole if blk_pg_nr == int(ts_name.split('_')[0])],
577 'mid_layer': [ts_name for ts_name in ins_th.mid_layer if ts_name[0] != 'w' and blk_pg_nr == int(ts_name.split('_')[0])
578 or ts_name[0] == 'w' and ts_name.split('_')[1][0] != 'w' and blk_pg_nr == int(ts_name.split('_')[1])],
579 'mid_layer_qh_i': [ts_name for ts_name in ins_th.mid_layer if ts_name.split('_')[1][0] != 'w' and blk_pg_nr == int(ts_name.split('_')[1])]}
580 ht_list = list(blk_pg.half_turns.keys())
581 self.ins_type_cond[str(blk_pg_nr)] = {'inner': [], 'outer': [], 'higher': [], 'lower': [], 'mid_turn': []}
582 for ht_nr, ht in blk_pg.half_turns.items():
583 if ht_nr in qh.turns and qh.N_strips > 0: # check if a quench heater strip touches the current half-turn
584 ht_index = qh.turns.index(ht_nr)
585 qh_side = qh.turns_sides[ht_index]
586 if qh.ids[ht_index] not in self.qh_data: self.qh_data[qh.ids[ht_index]] = {}
587 else: qh_side = ''
589 # find conductor type of ht adjacent to wedge
590 if ht_list.index(ht_nr) == len(ht_list) - 1 and ht.mid_turn_lines:
591 for line_name, line_tag in ht.mid_turn_lines.items():
592 for ts_name in ts_groups['mid_wedge_turn']:
593 if line_name in self.md.geometries.thin_shells.mid_wedge_turn[ts_name].lines:
594 self.wedge_cond[int(ts_name[1:ts_name.index('_')])] = blk_pg.conductor
595 break
597 self.ins_type_cond[str(blk_pg_nr)]['mid_turn'].extend(list(ht.mid_turn_lines.values())) # mid-turn insulation
599 if ht.aux_lines: # outer mid-layer insulation
600 _store_ts_tags(ht, geom_ts_name='mid_layers_aux', ts_grp='aux', lines='aux_lines')
602 if ht.mid_layer_lines.inner and qh_side == 'i':
603 for line_name, line_tag in ht.mid_layer_lines.inner.items():
604 for ts_name in ts_groups['mid_layer_qh_i']:
605 if ts_name in self.md.geometries.thin_shells.mid_layers_ht_to_ht:
606 ts_lines = self.md.geometries.thin_shells.mid_layers_ht_to_ht[ts_name].mid_layers.lines
607 elif ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_ht:
608 ts_lines = self.md.geometries.thin_shells.mid_layers_wdg_to_ht[ts_name].lines
609 else: ts_lines = []
610 if line_name in ts_lines:
611 _store_qh_data('internal', ts_name, line_tag)
612 break
613 elif ht.mid_layer_lines.inner: # block-coil (!) inner mid-layer insulation
614 _store_ts_tags(ht, geom_ts_name='mid_layers_ht_to_ht', geom_ts_name2='mid_layers', ts_grp='mid_layer', lines='mid_layer_lines', lines_side='inner')
615 elif not ht.mid_layer_lines.inner and qh_side == 'i': # quench heater inner insulation
616 _store_qh_data('external', blk_pg_nr, ht.lines['i'])
617 elif not ht.mid_layer_lines.inner: # inner insulation
618 self.ins_type_cond[str(blk_pg_nr)]['inner'].append(ht.lines['i'])
620 if ht.mid_layer_lines.outer and qh_side == 'o': # mid-layer quench heater insulation
621 for line_name, line_tag in ht.mid_layer_lines.outer.items():
622 for ts_name in ts_groups['mid_layer']:
623 if ts_name in self.md.geometries.thin_shells.mid_layers_ht_to_ht:
624 ts_lines = self.md.geometries.thin_shells.mid_layers_ht_to_ht[ts_name].mid_layers.lines
625 elif ts_name in self.md.geometries.thin_shells.mid_layers_ht_to_wdg:
626 ts_lines = self.md.geometries.thin_shells.mid_layers_ht_to_wdg[ts_name].lines
627 else: ts_lines = []
628 if line_name in ts_lines:
629 _store_qh_data('internal', ts_name, line_tag)
630 break
631 elif ht.mid_layer_lines.outer: # mid-layer insulation
632 for line_name, line_tag in ht.mid_layer_lines.outer.items():
633 for ts_name in ts_groups['mid_layer']:
634 if ts_name in self.md.geometries.thin_shells.mid_layers_ht_to_ht:
635 ts_lines = self.md.geometries.thin_shells.mid_layers_ht_to_ht[ts_name].mid_layers.lines
636 elif ts_name in self.md.geometries.thin_shells.mid_layers_ht_to_wdg:
637 ts_lines = self.md.geometries.thin_shells.mid_layers_ht_to_wdg[ts_name].lines
638 else: ts_lines = []
639 if line_name in ts_lines:
640 if ts_name not in self.ins_type['mid_layer']: self.ins_type['mid_layer'][ts_name] = []
641 self.ins_type['mid_layer'][ts_name].append(line_tag)
642 break
643 elif not ht.mid_layer_lines.outer and qh_side == 'o': # quench heater outer insulation
644 _store_qh_data('external', blk_pg_nr, ht.lines['o'])
645 else: # outer insulation
646 self.ins_type_cond[str(blk_pg_nr)]['outer'].append(ht.lines['o'])
648 # mid-pole insulation
649 _store_ts_tags(ht, geom_ts_name='mid_poles', ts_grp='mid_pole', lines='mid_pole_lines')
651 # mid-winding insulation
652 _store_ts_tags(ht, geom_ts_name='mid_windings', ts_grp='mid_winding', lines='mid_winding_lines')
654 if ht_list.index(ht_nr) == 0 and len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1: # lower angle external side insulation
655 if qh_side == 'l': _store_qh_data('external', blk_pg_nr, ht.lines['l'])
656 else: self.ins_type_cond[str(blk_pg_nr)]['lower'].append(ht.lines['l'])
657 if ht_list.index(ht_nr) == len(ht_list) - 1 and len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1: # higher angle external side insulation
658 if qh_side == 'h': _store_qh_data('external', blk_pg_nr, ht.lines['h'])
659 else: self.ins_type_cond[str(blk_pg_nr)]['higher'].append(ht.lines['h'])
661 # Wedges
662 for wdg_pg_nr, wdg_pg in pg.wedges.items():
663 ts_groups = {'aux': [ts_name for ts_name in ins_th.mid_layer if 'w' + str(wdg_pg_nr) in ts_name.split('_')],
664 'mid_layer': [ts_name for ts_name in ins_th.mid_layer if 'w' + str(wdg_pg_nr) == ts_name.split('_')[0]]}
665 self.ins_type_cond['w' + str(wdg_pg_nr)] = {'inner': [], 'outer': [], 'higher': [], 'lower': [], 'mid_turn': []}
666 if wdg_pg.aux_lines: _store_ts_tags(wdg_pg, geom_ts_name='mid_layers_aux', ts_grp='aux', lines='aux_lines')
667 if not wdg_pg.mid_layer_lines.inner: self.ins_type_cond['w' + str(wdg_pg_nr)]['inner'].append(wdg_pg.lines['i'])
668 if wdg_pg.mid_layer_lines.outer:
669 for line_name, line_tag in wdg_pg.mid_layer_lines.outer.items():
670 for ts_name in ts_groups['mid_layer']:
671 if ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_ht:
672 ts_lines = self.md.geometries.thin_shells.mid_layers_wdg_to_ht[ts_name].lines
673 elif ts_name in self.md.geometries.thin_shells.mid_layers_wdg_to_wdg:
674 ts_lines = self.md.geometries.thin_shells.mid_layers_wdg_to_wdg[ts_name].lines
675 else: ts_lines = []
676 if line_name in ts_lines:
677 if ts_name not in self.ins_type['mid_layer']: self.ins_type['mid_layer'][ts_name] = []
678 self.ins_type['mid_layer'][ts_name].append(line_tag)
679 break
680 else: self.ins_type_cond['w' + str(wdg_pg_nr)]['outer'].append(wdg_pg.lines['o'])
682 # Collect common thin shells for double qh mid-layers
683 for qh_nr, ts_groups in self.ins_type_qh['internal'].items():
684 for qh_nr2, ts_groups2 in self.ins_type_qh['internal'].items():
685 if qh_nr != qh_nr2:
686 common_ts_groups = list(set(ts_groups) & set(ts_groups2))
687 for ts in common_ts_groups:
688 tags, tags2 = ts_groups[ts], ts_groups2[ts]
689 common_tags = list(set(tags) & set(tags2))
690 for tag in common_tags:
691 tags.remove(tag), tags2.remove(tag)
692 if self.qh_data[qh_nr2][ts]['ht_side'] == 'i': qh_name = str(qh_nr) + '_' + str(qh_nr2)
693 else: qh_name = str(qh_nr2) + '_' + str(qh_nr)
694 if qh_name not in self.ins_type_qh['internal_double']: self.ins_type_qh['internal_double'][qh_name] = {}
695 qh_ins_id = self.ins_type_qh['internal_double'][qh_name]
696 if ts not in qh_ins_id: qh_ins_id[ts] = []
697 qh_ins_id[ts].append(tag)
699 def assignRegionsTags(self, geometry, mesh):
700 def _get_input_insulation_data(i_name, i_type=None):
701 ow_idx = next((index for index, couple in enumerate(self.data.magnet.solve.thermal.insulation_TSA.block_to_block.blocks_connection_overwrite)
702 if all(element in couple for element in i_name.split('_'))), None)
703 if i_type == 'mid_winding': mid_mat, mid_th = [self.data.magnet.solve.wedges.material], []
704 elif ow_idx is not None:
705 mid_mat = self.data.magnet.solve.thermal.insulation_TSA.block_to_block.materials_overwrite[ow_idx]
706 mid_th = self.data.magnet.solve.thermal.insulation_TSA.block_to_block.thicknesses_overwrite[ow_idx]
707 else: mid_mat, mid_th = [self.data.magnet.solve.thermal.insulation_TSA.block_to_block.material], []
708 return mid_mat, mid_th
710 def _compute_insulation_thicknesses(tot_th, known_ins_th):
711 if not mid_thicknesses:
712 mid_lyrs = [Func.sig_dig(tot_th - known_ins_th) / len(mid_materials)] * len(mid_materials)
713 elif None in mid_thicknesses:
714 input_ths = sum([th for th in mid_thicknesses if th is not None])
715 mid_lyrs = [th if th is not None else Func.sig_dig(tot_th - known_ins_th - input_ths) / mid_thicknesses.count(None) for th in mid_thicknesses]
716 else: mid_lyrs = mid_thicknesses
717 zeros = [nbr for nbr, th in enumerate(mid_lyrs) if th < 1e-8]
718 if tot_th - known_ins_th - sum(mid_lyrs) < -1e-8:
719 raise ValueError("Layer-to-layer insulation exceeds the space between blocks: check 'solve'->'insulation_TSA'->'block_to_block'->'thicknesses_overwrite'")
720 else: return mid_lyrs, zeros
722 pg = self.md.domains.physical_groups
723 qh = self.data.quench_protection.quench_heaters
725 # Air and air far field
726 if 'bore_field' in mesh.model_dump():
727 self.rm.air_far_field.vol.radius_out = float(abs(max(gmsh.model.getValue(0, gmsh.model.getAdjacencies(
728 1, self.md.geometries.air_inf.lines['outer'])[1][0], []), key=abs)))
729 self.rm.air_far_field.vol.radius_in = float(abs(max(gmsh.model.getValue(0, gmsh.model.getAdjacencies(
730 1, self.md.geometries.air_inf.lines['inner'])[1][0], []), key=abs)))
731 self.rm.air.vol.name = "Air"
732 self.rm.air.vol.number = pg.air
733 self.rm.air_far_field.vol.names = ["AirInf"]
734 self.rm.air_far_field.vol.numbers = [pg.air_inf]
735 self.rm.air_far_field.surf.name = "Surface_Inf"
736 self.rm.air_far_field.surf.number = pg.air_inf_bnd
737 if geometry.model_dump().get('symmetry', 'none') != 'none':
738 self.rm.boundaries.symmetry.normal_free.name = 'normal_free_bnd'
739 self.rm.boundaries.symmetry.normal_free.number = pg.symmetric_boundaries.normal_free
740 self.rm.boundaries.symmetry.tangential_free.name = 'tangent_free_bnd'
741 self.rm.boundaries.symmetry.tangential_free.number = pg.symmetric_boundaries.tangential_free
743 if 'use_TSA' in geometry.model_dump():
744 self.rm.projection_points.name = 'projection_points'
745 self.rm.projection_points.number = self.md.domains.physical_groups.half_turns_points
747 # Initialize lists
748 for group in ['r1_a1', 'r2_a1', 'r1_a2', 'r2_a2']:
749 self.rm.powered[group] = rM.Powered()
750 self.rm.powered[group].vol.names = []
751 self.rm.powered[group].vol.numbers = []
752 for cond_name in self.data.conductors.keys(): self.rm.powered[group].conductors[cond_name] = []
753 if geometry.with_wedges:
754 self.rm.induced[group] = rM.Induced()
755 self.rm.induced[group].vol.names = []
756 self.rm.induced[group].vol.numbers = []
757 if 'bore_field' in mesh.model_dump():
758 initial_current = self.data.power_supply.I_initial
759 self.rm.powered[group].vol.currents = []
760 if geometry.model_dump().get('use_TSA', False):
761 self.rm.powered[group].surf_in.names = []
762 self.rm.powered[group].surf_in.numbers = []
763 self.rm.powered[group].surf_out.names = []
764 self.rm.powered[group].surf_out.numbers = []
765 if geometry.with_wedges:
766 self.rm.induced[group].surf_in.names = []
767 self.rm.induced[group].surf_in.numbers = []
768 self.rm.induced[group].surf_out.names = []
769 self.rm.induced[group].surf_out.numbers = []
770 if geometry.with_iron_yoke:
771 self.rm.iron_yoke.vol.names = []
772 self.rm.iron_yoke.vol.numbers = []
773 if geometry.model_dump().get('use_TSA', False):
774 unique_thin_shells = []
775 self.rm.thin_shells.second_group_is_next['azimuthally'] = []
776 self.rm.thin_shells.second_group_is_next['radially'] = []
777 self.rm.thin_shells.normals_directed['azimuthally'] = []
778 self.rm.thin_shells.normals_directed['radially'] = []
779 else:
780 self.rm.insulator.vol.names = []
781 self.rm.insulator.vol.numbers = []
782 self.rm.insulator.surf.names = []
783 self.rm.insulator.surf.numbers = []
785 if geometry.model_dump().get('use_TSA', False):
786 # Categorize insulation types
787 min_h = mesh.insulation.global_size
788 # min_h = 1
789 # for conductor in self.data.conductors.keys():
790 # min_h = min([self.data.conductors[conductor].cable.th_insulation_along_height,
791 # self.data.conductors[conductor].cable.th_insulation_along_width, min_h])
792 min_h_QH = mesh.insulation.TSA.global_size_QH if mesh.insulation.TSA.global_size_QH else min_h
794 # Conductor insulation layers
795 for el, ins in self.ins_type_cond.items():
796 cond = self.data.conductors[self.wedge_cond[int(el[1:])]].cable if 'w' in el else self.data.conductors[pg.blocks[int(el)].conductor].cable
797 for ins_side, tags in ins.items():
798 if tags:
799 side_ins_type = [cond.material_insulation]
800 if ins_side in ['inner', 'outer']: side_ins = [cond.th_insulation_along_width]
801 elif ins_side in ['higher', 'lower']: side_ins = [cond.th_insulation_along_height]
802 else: # mid_turn
803 side_ins = [cond.th_insulation_along_height, cond.th_insulation_along_height]
804 side_ins_type.append(cond.material_insulation)
805 if ins_side[0] in 'iohl' and el + ins_side[0] in self.data.magnet.solve.thermal.insulation_TSA.exterior.blocks:
806 add_mat_idx = self.data.magnet.solve.thermal.insulation_TSA.exterior.blocks.index(el + ins_side[0])
807 side_ins.extend(self.data.magnet.solve.thermal.insulation_TSA.exterior.thicknesses_append[add_mat_idx])
808 side_ins_type.extend(self.data.magnet.solve.thermal.insulation_TSA.exterior.materials_append[add_mat_idx])
809 if ins_side[0] in 'il': side_ins.reverse(), side_ins_type.reverse()
810 self.rm.thin_shells.insulation_types.layers_number.append(0)
811 self.rm.thin_shells.insulation_types.thin_shells.append(list(set(tags)))
812 self.rm.thin_shells.insulation_types.thicknesses.append([])
813 self.rm.thin_shells.insulation_types.layers_material.append([])
814 for nr, ins_lyr in enumerate(side_ins):
815 tsa_layers = max(mesh.insulation.TSA.minimum_discretizations, round(ins_lyr / min_h))
816 self.rm.thin_shells.insulation_types.layers_number[-1] += tsa_layers
817 self.rm.thin_shells.insulation_types.thicknesses[-1].extend([ins_lyr / tsa_layers] * tsa_layers)
818 self.rm.thin_shells.insulation_types.layers_material[-1].extend([side_ins_type[nr]] * tsa_layers)
820 # Mid-pole, mid-winding, and mid-layer insulation layers
821 ins_th_dict = self.md.geometries.thin_shells.ins_thickness.model_dump()
822 for ins_type, ins in self.ins_type.items():
823 for ins_name, tags in ins.items():
824 # Get conductors insulation
825 if ins_name.count('w') == 2:
826 el1, el2 = int(ins_name[1:ins_name.index('_')]), int(ins_name[ins_name.index('_') + 2:])
827 cond1 = self.data.conductors[self.wedge_cond[el1]].cable
828 cond2 = self.data.conductors[self.wedge_cond[el2]].cable
829 elif ins_name.count('w') == 1:
830 el1, el2 = int(ins_name[1:ins_name.index('_')]), int(ins_name[ins_name.index('_') + 1:])
831 cond1 = self.data.conductors[self.wedge_cond[el1]].cable
832 cond2 = self.data.conductors[pg.blocks[el2].conductor].cable
833 if ins_name in self.md.geometries.thin_shells.mid_layers_ht_to_wdg: cond1, cond2 = cond2, cond1
834 else:
835 el1, el2 = int(ins_name[:ins_name.index('_')]), int(ins_name[ins_name.index('_') + 1:])
836 cond1 = self.data.conductors[pg.blocks[el1].conductor].cable
837 cond2 = self.data.conductors[pg.blocks[el2].conductor].cable
838 if ins_type in ['mid_layer', 'aux']:
839 cond_ins1, cond_ins2 = cond1.th_insulation_along_width, cond2.th_insulation_along_width
840 else:
841 cond_ins1, cond_ins2 = cond1.th_insulation_along_height, cond2.th_insulation_along_height
842 # Get insulation layer thickness
843 mid_materials, mid_thicknesses = _get_input_insulation_data(ins_name, i_type=ins_type)
844 ins_thickness = ins_th_dict['mid_layer'][ins_name] / 2 if ins_type == 'aux' else ins_th_dict[ins_type][ins_name]
845 mid_lyr_th, null_idx = _compute_insulation_thicknesses(ins_thickness, cond_ins1 + cond_ins2)
846 for idx in null_idx: mid_lyr_th.pop(idx), mid_materials.pop(idx)
847 side_ins = [cond_ins1] + mid_lyr_th + [cond_ins2]
848 # Get insulation materials
849 side_ins_type = [cond1.material_insulation] + mid_materials + [cond2.material_insulation]
850 # Initialize sublists
851 self.rm.thin_shells.insulation_types.layers_number.append(0)
852 self.rm.thin_shells.insulation_types.thin_shells.append(list(set(tags)))
853 self.rm.thin_shells.insulation_types.thicknesses.append([])
854 self.rm.thin_shells.insulation_types.layers_material.append([])
855 for nr, ins_lyr in enumerate(side_ins):
856 tsa_layers = max(mesh.insulation.TSA.minimum_discretizations, round(ins_lyr / min_h))
857 self.rm.thin_shells.insulation_types.layers_number[-1] += tsa_layers
858 self.rm.thin_shells.insulation_types.thicknesses[-1].extend([ins_lyr / tsa_layers] * tsa_layers)
859 self.rm.thin_shells.insulation_types.layers_material[-1].extend([side_ins_type[nr]] * tsa_layers)
861 # Quench heater insulation layers
862 for ins_type, ins in self.ins_type_qh.items():
863 for qh_nr, ts_groups in ins.items():
864 for ts_name, tags in ts_groups.items():
865 if ins_type != 'external': mid_materials, mid_thicknesses = _get_input_insulation_data(ts_name)
866 if tags:
867 if ins_type == 'external':
868 data = self.qh_data[qh_nr]
869 if str(ts_name) + data[ts_name]['ht_side'] in self.data.magnet.solve.thermal.insulation_TSA.exterior.blocks:
870 add_mat_idx = self.data.magnet.solve.thermal.insulation_TSA.exterior.blocks.index(str(ts_name) + data[ts_name]['ht_side'])
871 additional_ths = self.data.magnet.solve.thermal.insulation_TSA.exterior.thicknesses_append[add_mat_idx]
872 additional_mats = self.data.magnet.solve.thermal.insulation_TSA.exterior.materials_append[add_mat_idx]
873 else: additional_ths, additional_mats = [], []
874 cond = self.data.conductors[data[ts_name]['conductor']].cable
875 ht_ins_th = cond.th_insulation_along_width if data[ts_name]['ht_side'] in 'io' else cond.th_insulation_along_height
876 side_ins = [ht_ins_th] + [h_ins for h_ins in qh.h_ins[qh_nr - 1]] + [qh.h[qh_nr - 1]] + [h_ground_ins for h_ground_ins in qh.h_ground_ins[qh_nr - 1]] + additional_ths
877 side_ins_type = [cond.material_insulation] + [type_ins for type_ins in qh.type_ins[qh_nr - 1]] + ['SS'] +\
878 [type_ground_ins for type_ground_ins in qh.type_ground_ins[qh_nr - 1]] + additional_mats
879 if data[ts_name]['ht_side'] in 'il': side_ins.reverse(), side_ins_type.reverse()
880 qh_list = [qh_nr]
881 elif ins_type == 'internal':
882 data = self.qh_data[qh_nr]
883 cond = self.data.conductors[data[ts_name]['conductor']].cable
884 cond2 = self.data.conductors[self.wedge_cond[int(ts_name.split('_')[0][1:])] if 'w' in ts_name
885 else pg.blocks[int(ts_name.split('_')[1]) if data[ts_name]['ht_side'] == 'o' else int(ts_name.split('_')[0])].conductor].cable
886 side_ins_qh = [h_ins for h_ins in qh.h_ins[qh_nr - 1]] + [qh.h[qh_nr - 1]] + [h_ground_ins for h_ground_ins in qh.h_ground_ins[qh_nr - 1]]
887 mid_lyr_th, null_idx = _compute_insulation_thicknesses(
888 ins_th_dict['mid_layer'][ts_name], sum([cond.th_insulation_along_width, cond2.th_insulation_along_width] + side_ins_qh))
889 for idx in null_idx: mid_lyr_th.pop(idx), mid_materials.pop(idx)
890 side_ins = [cond.th_insulation_along_width] + side_ins_qh + mid_lyr_th + [cond2.th_insulation_along_width]
891 side_ins_type = [cond.material_insulation] + [type_ins for type_ins in qh.type_ins[qh_nr - 1]] + ['SS'] +\
892 [type_ground_ins for type_ground_ins in qh.type_ground_ins[qh_nr - 1]] + mid_materials + [cond2.material_insulation]
893 if data[ts_name]['ht_side'] == 'i': side_ins.reverse(), side_ins_type.reverse()
894 qh_list = [qh_nr]
895 elif ins_type == 'internal_double':
896 qh_nr1, qh_nr2 = int(qh_nr.split('_')[0]), int(qh_nr.split('_')[1])
897 data, data2 = self.qh_data[qh_nr1], self.qh_data[qh_nr2]
898 cond = self.data.conductors[data[ts_name]['conductor']].cable
899 cond2 = self.data.conductors[data2[ts_name]['conductor']].cable
900 side_ins_qh = [h_ins for h_ins in qh.h_ins[qh_nr1 - 1]] + [qh.h[qh_nr1 - 1]] + [h_ground_ins for h_ground_ins in qh.h_ground_ins[qh_nr1 - 1]]
901 side_ins_qh2 = [h_ground_ins for h_ground_ins in qh.h_ground_ins[qh_nr2 - 1][::-1]] + [qh.h[qh_nr2 - 1]] + [h_ins for h_ins in qh.h_ins[qh_nr2 - 1][::-1]]
902 mid_lyr_th, null_idx = _compute_insulation_thicknesses(
903 ins_th_dict['mid_layer'][ts_name], sum([cond.th_insulation_along_width, cond2.th_insulation_along_width] + side_ins_qh + side_ins_qh2))
904 for idx in null_idx: mid_lyr_th.pop(idx), mid_materials.pop(idx)
905 side_ins = [cond.th_insulation_along_width] + side_ins_qh + mid_lyr_th + side_ins_qh2 + [cond2.th_insulation_along_width]
906 side_ins_type = [cond.material_insulation] + [type_ins for type_ins in qh.type_ins[qh_nr1 - 1]] + ['SS'] +\
907 [type_ground_ins for type_ground_ins in qh.type_ground_ins[qh_nr1 - 1]] + mid_materials +\
908 [type_ground_ins for type_ground_ins in qh.type_ground_ins[qh_nr2 - 1][::-1]] + ['SS'] +\
909 [type_ins for type_ins in qh.type_ins[qh_nr2 - 1][::-1]] + [cond2.material_insulation]
910 qh_list = [qh_nr1, qh_nr2]
911 qh_labels = [1 if m == 'SS' else None for m in side_ins_type]
912 ss_indexes = [index for index, value in enumerate(qh_labels) if value == 1]
913 for nr, idx in enumerate(ss_indexes): qh_labels[idx] = qh_list[nr]
914 self.rm.thin_shells.quench_heaters.layers_number.append(0)
915 self.rm.thin_shells.quench_heaters.thin_shells.append(list(set(tags)))
916 self.rm.thin_shells.quench_heaters.thicknesses.append([])
917 self.rm.thin_shells.quench_heaters.layers_material.append([])
918 self.rm.thin_shells.quench_heaters.label.append([])
919 for nr, ins_lyr in enumerate(side_ins):
920 tsa_layers = max(mesh.insulation.TSA.minimum_discretizations_QH, round(ins_lyr / min_h_QH))
921 self.rm.thin_shells.quench_heaters.layers_number[-1] += tsa_layers
922 self.rm.thin_shells.quench_heaters.thicknesses[-1].extend([ins_lyr / tsa_layers] * tsa_layers)
923 self.rm.thin_shells.quench_heaters.layers_material[-1].extend([side_ins_type[nr]] * tsa_layers)
924 self.rm.thin_shells.quench_heaters.label[-1].extend([qh_labels[nr]] * tsa_layers)
926 # Powered
927 for blk_nr, blk in pg.blocks.items():
928 ht_list = list(blk.half_turns.keys())
929 for ht_nr, ht in blk.half_turns.items():
930 ht_name = f"ht{ht_nr}_{'EM' if 'bore_field' in mesh.model_dump() else 'TH'}"
931 self.rm.powered[ht.group].conductors[blk.conductor].append(ht_name)
932 self.rm.powered[ht.group].vol.names.append(ht_name)
933 self.rm.powered[ht.group].vol.numbers.append(ht.tag)
934 if 'bore_field' in mesh.model_dump():
935 self.rm.powered[ht.group].vol.currents.append(initial_current * (1 if blk.current_sign > 0 else -1))
936 if geometry.model_dump().get('use_TSA', False):
937 for line in ['l', 'i', 'o', 'h']:
938 # Bare edges
939 self.rm.powered[ht.group].surf_in.names.append(ht_name + line)
940 self.rm.powered[ht.group].surf_in.numbers.append(ht.lines[line])
941 if line in 'io': self.rm.thin_shells.normals_directed['radially'].append(ht.lines[line])
942 else: self.rm.thin_shells.normals_directed['azimuthally'].append(ht.lines[line])
943 # Auxiliary thin shells
944 for line_name, line_tag in ht.aux_lines.items():
945 self.rm.powered[ht.group].surf_in.names.append(line_name)
946 self.rm.powered[ht.group].surf_in.numbers.append(line_tag)
947 self.rm.thin_shells.normals_directed['radially'].append(line_tag)
948 # Thin shells
949 for line_name, line_tag in dict(ht.mid_layer_lines.inner, **ht.mid_layer_lines.outer).items():
950 self.rm.powered[ht.group].surf_out.names.append(line_name)
951 self.rm.powered[ht.group].surf_out.numbers.append(line_tag)
952 self.rm.thin_shells.normals_directed['radially'].append(line_tag)
953 for line_name, line_tag in dict(ht.mid_pole_lines, **ht.mid_winding_lines, **ht.mid_turn_lines).items():
954 self.rm.powered[ht.group].surf_out.names.append(line_name)
955 self.rm.powered[ht.group].surf_out.numbers.append(line_tag)
956 self.rm.thin_shells.normals_directed['azimuthally'].append(line_tag)
958 # Which thin shells or exterior conductor edges precede a second group (r2 or a2)
959 if ht.group[1] == '2':
960 # mid-turn thin shells precede r2
961 for line_name, line_tag in ht.mid_turn_lines.items():
962 if (ht_list.index(ht_nr) != 0 and int(line_name[line_name.index('_') + 1:]) == ht_nr) or\
963 (ht_list.index(ht_nr) == 0 and 'w' in line_name):
964 self.rm.thin_shells.second_group_is_next['azimuthally'].append(line_tag)
965 # mid-pole thin shells precede r2
966 if ht_list.index(ht_nr) == 0 and ht.mid_pole_lines:
967 self.rm.thin_shells.second_group_is_next['azimuthally'].append(list(ht.mid_pole_lines.values())[0])
968 # conductor edges precede r2
969 elif ht_list.index(ht_nr) == 0 and len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1:
970 self.rm.thin_shells.second_group_is_next['azimuthally'].append(ht.lines['l'])
971 # mid-winding thin shells precede r2
972 for line_name, line_tag in ht.mid_winding_lines.items():
973 if int(line_name[line_name.index('_') + 1:]) == ht_nr:
974 self.rm.thin_shells.second_group_is_next['azimuthally'].append(line_tag)
975 elif ht_list.index(ht_nr) == len(ht_list) - 1:
976 # mid-turn thin shells precede r2
977 for line_name, line_tag in ht.mid_turn_lines.items():
978 if 'w' in line_name:
979 self.rm.thin_shells.second_group_is_next['azimuthally'].append(line_tag)
980 # conductor edges precede r2
981 if len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1:
982 self.rm.thin_shells.second_group_is_next['azimuthally'].append(ht.lines['h'])
983 if ht.group[4] == '2':
984 # mid-layer thin shells precede a2
985 for line_name, line_tag in ht.mid_layer_lines.inner.items():
986 self.rm.thin_shells.second_group_is_next['radially'].append(line_tag)
987 elif not ht.mid_layer_lines.outer:
988 # conductor edges precede a2
989 self.rm.thin_shells.second_group_is_next['radially'].append(ht.lines['o'])
991 if geometry.model_dump().get('use_TSA', False):
992 for group in ['r1_a1', 'r2_a1', 'r1_a2', 'r2_a2']:
993 unique_thin_shells.extend(self.rm.powered[group].surf_out.numbers)
995 # Wedges
996 if geometry.with_wedges:
997 for wdg_nr, wdg in pg.wedges.items():
998 wdg_name = f"w{wdg_nr}_{'EM' if 'bore_field' in mesh.model_dump() else 'TH'}"
999 self.rm.induced[wdg.group].vol.names.append(wdg_name)
1000 self.rm.induced[wdg.group].vol.numbers.append(wdg.tag)
1001 if geometry.model_dump().get('use_TSA', False):
1002 # Bare edges
1003 for line in ['l', 'i', 'o', 'h']:
1004 self.rm.induced[wdg.group].surf_in.names.append(wdg_name + line)
1005 self.rm.induced[wdg.group].surf_in.numbers.append(wdg.lines[line])
1006 if line in 'io':
1007 self.rm.thin_shells.normals_directed['radially'].append(wdg.lines[line])
1008 else:
1009 self.rm.thin_shells.normals_directed['azimuthally'].append(wdg.lines[line])
1010 # Auxiliary thin shells
1011 for line_name, line_tag in wdg.aux_lines.items():
1012 self.rm.induced[wdg.group].surf_in.names.append(line_name)
1013 self.rm.induced[wdg.group].surf_in.numbers.append(line_tag)
1014 self.rm.thin_shells.normals_directed['radially'].append(line_tag)
1015 # Thin shells
1016 for line_name, line_tag in dict(wdg.mid_layer_lines.inner, **wdg.mid_layer_lines.outer).items():
1017 self.rm.induced[wdg.group].surf_out.names.append(line_name)
1018 self.rm.induced[wdg.group].surf_out.numbers.append(line_tag)
1019 self.rm.thin_shells.normals_directed['radially'].append(line_tag)
1020 for line_name, line_tag in wdg.mid_turn_lines.items():
1021 self.rm.induced[wdg.group].surf_out.names.append(line_name)
1022 self.rm.induced[wdg.group].surf_out.numbers.append(line_tag)
1023 self.rm.thin_shells.normals_directed['azimuthally'].append(line_tag)
1024 # Which thin shells or exterior conductor edges precede a second group (r2 or a2)
1025 if wdg.group[4] == '2':
1026 for line_name, line_tag in wdg.mid_layer_lines.inner.items():
1027 if line_name.count('w') == 2:
1028 self.rm.thin_shells.second_group_is_next['radially'].append(line_tag)
1029 elif not wdg.mid_layer_lines.outer:
1030 self.rm.thin_shells.second_group_is_next['radially'].append(wdg.lines['o'])
1031 if geometry.model_dump().get('use_TSA', False):
1032 for group in ['r1_a1', 'r2_a1', 'r1_a2', 'r2_a2']:
1033 unique_thin_shells.extend(self.rm.induced[group].surf_out.numbers)
1035 # Unique mid layers
1036 if geometry.model_dump().get('use_TSA', False):
1037 self.rm.thin_shells.mid_turns_layers_poles = list(set(unique_thin_shells))
1039 # Insulation
1040 for group_name, surface in pg.insulations.surfaces.items():
1041 self.rm.insulator.vol.names.append('ins' + group_name)
1042 self.rm.insulator.vol.numbers.append(surface)
1043 if 'insulation' in mesh.model_dump() and 'TSA' in mesh.model_dump()["insulation"]:
1044 for group_name, curve in pg.insulations.curves.items():
1045 self.rm.insulator.surf.names.append('ins' + group_name)
1046 self.rm.insulator.surf.numbers.append(curve)
1048 # Iron
1049 for group_name, surface in pg.iron.surfaces.items():
1050 self.rm.iron_yoke.vol.names.append(group_name)
1051 self.rm.iron_yoke.vol.numbers.append(surface)
1053 # Boundary conditions
1054 if 'insulation' in mesh.model_dump() and 'TSA' in mesh.model_dump()["insulation"]:
1055 # Initialize lists
1056 for bc_data, bc_rm in zip(self.data.magnet.solve.thermal.overwrite_boundary_conditions, self.rm.boundaries.thermal): # b.c. type
1057 bc_rm[1].bc.names = []
1058 bc_rm[1].bc.numbers = []
1059 if bc_data[0] == 'cooling':
1060 bc_rm[1].bc.values = []
1061 for group in ['1_r1_a1', '2_r1_a1', '1_r2_a1', '2_r2_a1', '1_r1_a2', '2_r1_a2', '1_r2_a2', '2_r2_a2']:
1062 bc_rm[1].groups[group] = []
1063 else:
1064 bc_rm[1].bc.value = []
1065 for group in ['r1_a1', 'r2_a1', 'r1_a2', 'r2_a2']:
1066 bc_rm[1].groups[group] = []
1068 # Apply general cooling and adiabatic
1069 if self.data.magnet.solve.thermal.He_cooling.enabled:
1070 cooling_side = {'i': any(coil_side in self.data.magnet.solve.thermal.He_cooling.sides for coil_side in ['inner', 'external']),
1071 'o': any(coil_side in self.data.magnet.solve.thermal.He_cooling.sides for coil_side in ['outer', 'external']),
1072 'hl': self.data.magnet.solve.thermal.He_cooling.sides == 'external'}
1073 else:
1074 cooling_side = {'i': False, 'o': False, 'hl': False}
1076 def __assign_bnd_tag(el, name, side, bc_type, tag=None):
1077 line_tag = tag if tag else el.lines[side]
1078 bnd_list_names[bc_type].append(name + side)
1079 bnd_list_numbers[bc_type].append(line_tag)
1080 if side in 'io': new_group = el.group[:3] + 'a1' if el.group[4] == '2' else el.group[:3] + 'a2'
1081 else: new_group = 'r1' + el.group[2:] if el.group[1] == '2' else 'r2' + el.group[2:]
1082 bc_rm[bc_type].groups[new_group].append(line_tag)
1083 for group_name, group in self.rm.thin_shells.groups.items():
1084 if line_tag in group:
1085 bc_rm[bc_type].groups[el.group[0] + '_' + new_group].append(line_tag)
1086 break
1088 bc_rm = {'Robin': self.rm.boundaries.thermal.cooling, 'Neumann': self.rm.boundaries.thermal.heat_flux}
1089 bnd_list_names = {'Robin': [], 'Neumann': []}
1090 bnd_list_numbers = {'Robin': [], 'Neumann': []}
1091 if geometry.model_dump().get('use_TSA', False):
1092 # Half turn boundaries
1093 for coil_nr, coil in self.md.geometries.coil.anticlockwise_order.coils.items():
1094 for lyr_nr, orders in coil.layers.items():
1095 for order in orders:
1096 ht_list = list(pg.blocks[order.block].half_turns.keys())
1097 for ht_nr, ht in pg.blocks[order.block].half_turns.items():
1098 if not ht.mid_layer_lines.inner:
1099 __assign_bnd_tag(ht, 'ht' + str(ht_nr), 'i', 'Robin' if cooling_side['i'] else 'Neumann')
1100 if not ht.mid_layer_lines.outer:
1101 __assign_bnd_tag(ht, 'ht' + str(ht_nr), 'o', 'Robin' if cooling_side['o'] else 'Neumann')
1102 if ht_list.index(ht_nr) == 0 and len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1:
1103 __assign_bnd_tag(ht, 'ht' + str(ht_nr), 'l', 'Robin' if cooling_side['hl'] else 'Neumann')
1104 if ht_list.index(ht_nr) == len(ht_list) - 1 and len(ht.mid_turn_lines) + len(ht.mid_winding_lines) + len(ht.mid_pole_lines) == 1:
1105 __assign_bnd_tag(ht, 'ht' + str(ht_nr), 'h', 'Robin' if cooling_side['hl'] else 'Neumann')
1106 if ht.aux_lines:
1107 __assign_bnd_tag(ht, 'ht' + str(ht_nr), 'o', 'Robin' if cooling_side['hl'] else 'Neumann', list(ht.aux_lines.values())[0])
1108 # Wedge boundaries
1109 for wdg_nr, wdg in pg.wedges.items():
1110 if not wdg.mid_layer_lines.inner:
1111 __assign_bnd_tag(wdg, 'wd' + str(wdg_nr), 'i', 'Robin' if cooling_side['i'] else 'Neumann')
1112 if not wdg.mid_layer_lines.outer:
1113 __assign_bnd_tag(wdg, 'wd' + str(wdg_nr), 'o', 'Robin' if cooling_side['o'] else 'Neumann')
1114 if wdg.aux_lines:
1115 __assign_bnd_tag(wdg, 'wd' + str(wdg_nr), 'o', 'Robin' if cooling_side['hl'] else 'Neumann', list(wdg.aux_lines.values())[0])
1116 else: # insulation case
1117 for curves_group, tag in pg.insulations.curves.items():
1118 if self.data.magnet.solve.thermal.He_cooling.enabled:
1119 if self.data.magnet.solve.thermal.He_cooling.sides == 'external': bc_type = 'Robin'
1120 else:
1121 raise ValueError(f"Cooling side '{self.data.magnet.solve.thermal.He_cooling.sides}' is not supported for meshed insulation models.")
1122 # bc_type = 'Robin' if (self.data.magnet.solve.thermal.He_cooling.sides == 'external' or
1123 # ('inner' in self.data.magnet.solve.thermal.He_cooling.sides and curves_group[-1] == 'i') or
1124 # ('outer' in self.data.magnet.solve.thermal.He_cooling.sides and curves_group[-1] == 'o')) else 'Neumann'
1125 else:
1126 bc_type = 'Neumann'
1127 bnd_list_names[bc_type].append('ins' + curves_group)
1128 bnd_list_numbers[bc_type].append(tag)
1129 if self.data.magnet.solve.thermal.He_cooling.enabled:
1130 bc_rm['Robin'].bc.names.append(bnd_list_names['Robin'])
1131 bc_rm['Robin'].bc.numbers.append(bnd_list_numbers['Robin'])
1132 bc_rm['Robin'].bc.values.append([self.data.magnet.solve.thermal.He_cooling.heat_transfer_coefficient,
1133 self.data.magnet.solve.thermal.init_temperature])
1134 if bnd_list_names['Neumann']:
1135 bc_rm['Neumann'].bc.names.append(bnd_list_names['Neumann'])
1136 bc_rm['Neumann'].bc.numbers.append(bnd_list_numbers['Neumann'])
1137 bc_rm['Neumann'].bc.value.append(0.)
1139 # Apply specific boundary conditions
1140 for bc_data, bc_rm in zip(self.data.magnet.solve.thermal.overwrite_boundary_conditions, self.rm.boundaries.thermal): # b.c. type
1141 # bc_data is a tuple like: ('temperature', {'const_T1': boundaries, value)})
1142 # bc_rm is a tuple like: ('temperature', DirichletCondition(names, numbers, value))
1144 for _, bc in bc_data[1].items(): # all boundary conditions of one b.c. type (e.g., Dirichlet with different temperatures)
1145 bnd_list_names = []
1146 bnd_list_numbers = []
1147 if geometry.model_dump().get('use_TSA', False):
1148 for bnd in bc.boundaries: # all boundaries of one boundary condition
1149 if bnd[0] == 'w':
1150 if not geometry.with_wedges:
1151 raise Exception('Wedge regions are disabled.')
1152 # Fetch the physical group of the wedge
1153 pg_el = pg.wedges[int(bnd[1:-1])]
1154 name = 'wd' + bnd[1:]
1155 else:
1156 # Fetch the physical group of the half turn
1157 ht_index = self.strands['ht'].index(int(bnd[:-1]))
1158 pg_el = pg.blocks[self.strands['block'][ht_index]].half_turns[int(bnd[:-1])]
1159 name = 'ht' + bnd
1160 line_pg_tag = pg_el.lines[bnd[-1]]
1161 bnd_list_names.append(name)
1162 bnd_list_numbers.append(line_pg_tag)
1163 # Find the half turn group this boundary is assigned to and take the complementary
1164 if bnd[-1] in 'io':
1165 new_group = pg_el.group[:3] + 'a1' if pg_el.group[4] == '2' else pg_el.group[:3] + 'a2'
1166 else: # ['l', 'h'] todo: if applied to an inner line (i.e., not domain boundaries), extra code needed because the line would belong to two groups
1167 new_group = 'r1' + pg_el.group[2:] if pg_el.group[1] == '2' else 'r2' + pg_el.group[2:]
1168 # Overwrite general cooling and adiabatic condition
1169 if self.data.magnet.solve.thermal.He_cooling.enabled:
1170 if name in self.rm.boundaries.thermal.cooling.bc.names[0]:
1171 bnd_idx = self.rm.boundaries.thermal.cooling.bc.names[0].index(name)
1172 self.rm.boundaries.thermal.cooling.bc.names[0].pop(bnd_idx)
1173 bnd_idx = self.rm.boundaries.thermal.cooling.bc.numbers[0].pop(bnd_idx)
1174 self.rm.boundaries.thermal.cooling.groups[new_group].pop(
1175 self.rm.boundaries.thermal.cooling.groups[new_group].index(bnd_idx))
1176 if self.data.magnet.solve.thermal.He_cooling.sides != 'external':
1177 if name in self.rm.boundaries.thermal.heat_flux.bc.names[0]:
1178 bnd_idx = self.rm.boundaries.thermal.heat_flux.bc.names[0].index(name)
1179 self.rm.boundaries.thermal.heat_flux.bc.names[0].pop(bnd_idx)
1180 bnd_idx = self.rm.boundaries.thermal.heat_flux.bc.numbers[0].pop(bnd_idx)
1181 self.rm.boundaries.thermal.heat_flux.groups[new_group].pop(
1182 self.rm.boundaries.thermal.heat_flux.groups[new_group].index(bnd_idx))
1183 # Assign the tag
1184 bc_rm[1].groups[new_group].append(line_pg_tag)
1185 # Extra grouping for Robin virtual shells
1186 if bc_data[0] == 'cooling':
1187 for group_name, group in self.rm.thin_shells.groups.items():
1188 if line_pg_tag in group:
1189 bc_rm[1].groups[group_name[0] + '_' + new_group].append(line_pg_tag)
1190 break
1191 else: # the b.c. are assigned to insulation boundaries instead
1192 pass # todo: not supported yet
1193 bc_rm[1].bc.names.append(bnd_list_names)
1194 bc_rm[1].bc.numbers.append(bnd_list_numbers)
1195 if bc_data[0] == 'cooling':
1196 bc_rm[1].bc.values.append([bc.heat_transfer_coefficient, self.data.magnet.solve.thermal.init_temperature])
1197 elif bc_data[0] == 'temperature':
1198 bc_rm[1].bc.value.append(bc.const_temperature)
1199 elif bc_data[0] == 'heat_flux':
1200 bc_rm[1].bc.value.append(bc.const_heat_flux)
1202 def setMeshOptions(self, run_type):
1203 """
1204 Meshes the generated domain
1205 """
1206 mesh = self.data.magnet.mesh
1207 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
1208 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
1209 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
1210 gmsh.option.setNumber("Mesh.Algorithm", 6)
1211 gmsh.option.setNumber("Mesh.Optimize", 1)
1212 gmsh.option.setNumber("Mesh.ElementOrder", 1)
1214 def generateMesh(self):
1215 gmsh.option.setNumber("General.Terminal", self.verbose)
1216 self.mesh.generate(2)
1217 # self.mesh.removeDuplicateNodes()
1218 # self.occ.synchronize()
1219 # self.gu.launch_interactive_GUI()
1221 def checkMeshQuality(self):
1222 tags = self.mesh.getElements(2)[1][0]
1224 self.mesh_parameters['SJ'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSJ'))
1225 self.mesh_parameters['SICN'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSICN'))
1226 self.mesh_parameters['SIGE'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='minSIGE'))
1227 self.mesh_parameters['Gamma'] = min(self.mesh.getElementQualities(elementTags=tags, qualityName='gamma'))
1228 self.mesh_parameters['nodes'] = len(self.mesh.getNodes()[0])
1230 # gmsh.plugin.setNumber("AnalyseMeshQuality", "JacobianDeterminant", 1)
1231 # gmsh.plugin.setNumber("AnalyseMeshQuality", "CreateView", 100)
1232 # test = gmsh.plugin.run("AnalyseMeshQuality")
1233 # test2 = gmsh.view.getModelData(test, test)
1235 # gmsh.logger.getLastError()
1236 # gmsh.logger.get()
1238 def saveClosestNeighboursList(self):
1240 def _closest_node_on_reference(origin_points, reference_points):
1241 """
1242 Compute list of lists with each list representing one original node (3 entries for 3 coordinates) and the closest
1243 reference mesh node (3 entries for 3 coordinates).
1244 :param origin_points: list of origin point as list of 3 floats per origin node
1245 :param reference_points: list of coordinates of reference mesh points as list of 3 floats per reference mesh node
1246 :return: returns list of lists with the closest reference mesh node per origin node
1247 """
1248 closest_node = []
1249 for x in range(0, len(origin_points), 3):
1250 origin_point = origin_points[x:x + 3]
1251 # Compute distance list between origin point and reference point list
1252 dist_lst = [Func.points_distance(origin_point, reference_points[y:y + 3]) for y in range(0, len(reference_points), 3)]
1253 min_idx = 3 * np.argmin(dist_lst)
1254 closest_node.append(origin_point + reference_points[min_idx:min_idx + 3])
1255 return closest_node
1257 def _get_closest_nodes(side):
1258 origin_list = self.mesh.getNodesForPhysicalGroup(1, mid_layer)[1].tolist()
1259 reference_list = self.mesh.getNodesForPhysicalGroup(1, el.lines[side])[1].tolist()
1260 self.rc.neighbouring_nodes.groups[group].extend(
1261 [node for node_list in _closest_node_on_reference(origin_list, reference_list) for node in node_list])
1263 logger.info(f"Info : {self.data.general.magnet_name} - F i n d i n g C l o s e s t N e i g h b o u r s . . .")
1264 logger.info(f"Info : Finding closest reference nodes ...")
1266 self.rc.neighbouring_nodes.groups['1_1'] = []
1267 self.rc.neighbouring_nodes.groups['2_1'] = []
1268 self.rc.neighbouring_nodes.groups['1_2'] = []
1269 self.rc.neighbouring_nodes.groups['2_2'] = []
1271 for blk_nr, blk in self.md.domains.physical_groups.blocks.items():
1272 for ht_nr, el in blk.half_turns.items():
1273 ht_list = list(blk.half_turns.keys())
1274 group = el.group[1] + '_' + el.group[-1]
1275 for line_name, mid_layer in el.mid_layer_lines.inner.items(): _get_closest_nodes('i')
1276 for line_name, mid_layer in el.mid_layer_lines.outer.items(): _get_closest_nodes('o')
1277 for line_name, mid_layer in el.aux_lines.items(): _get_closest_nodes('o')
1278 for line_name, mid_layer in el.mid_pole_lines.items(): _get_closest_nodes('l' if ht_list.index(ht_nr) == 0 else 'h')
1279 for line_name, mid_layer in el.mid_winding_lines.items(): _get_closest_nodes('l' if ht_list.index(ht_nr) == 0 else 'h')
1280 for line_name, mid_layer in el.mid_turn_lines.items():
1281 high = ht_list.index(ht_nr) == len(ht_list) - 1 if 'w' in line_name\
1282 else int(line_name[:line_name.index('_')]) == ht_nr
1283 _get_closest_nodes('h' if high else 'l')
1284 for wdg_nr, el in self.md.domains.physical_groups.wedges.items():
1285 group = el.group[1] + '_' + el.group[-1]
1286 for line_name, mid_layer in el.mid_layer_lines.inner.items(): _get_closest_nodes('i')
1287 for line_name, mid_layer in el.mid_layer_lines.outer.items(): _get_closest_nodes('o')
1288 for line_name, mid_layer in el.aux_lines.items(): _get_closest_nodes('o')
1289 for line_name, mid_layer in el.mid_turn_lines.items():
1290 _get_closest_nodes('l' if line_name == list(el.mid_turn_lines.keys())[0] else 'h')
1292 logger.info(f"Info : {self.data.general.magnet_name} - E n d F i n d i n g C l o s e s t N e i g h b o u r s")
1294 def selectMeshNodes(self, elements: str):
1296 logger.info(f"Info : {self.data.general.magnet_name} - S e l e c t i n g M e s h N o d e s . . .")
1297 logger.info(f"Info : Selecting a mesh node per isothermal {elements[:-1]} ...")
1299 if elements == 'conductors':
1300 bare_mesh = {'1_1': self.rm.powered['r1_a1'].vol.numbers, '2_1': self.rm.powered['r2_a1'].vol.numbers,
1301 '1_2': self.rm.powered['r1_a2'].vol.numbers, '2_2': self.rm.powered['r2_a2'].vol.numbers}
1302 groups = self.rc.isothermal_nodes.conductors
1304 # dir + robin + mid_layers
1305 # potentially all thin shell lines if easier
1306 line_tags = self.rm.boundaries.thermal.temperature.groups['r1_a1'] + \
1307 self.rm.boundaries.thermal.temperature.groups['r1_a2'] + \
1308 self.rm.boundaries.thermal.temperature.groups['r2_a1'] + \
1309 self.rm.boundaries.thermal.temperature.groups['r2_a2'] + \
1310 self.rm.boundaries.thermal.cooling.groups['r1_a1'] + \
1311 self.rm.boundaries.thermal.cooling.groups['r1_a2'] + \
1312 self.rm.boundaries.thermal.cooling.groups['r2_a1'] + \
1313 self.rm.boundaries.thermal.cooling.groups['r2_a2'] + \
1314 self.rm.thin_shells.mid_turns_layers_poles
1315 for tag in line_tags:
1316 coords = list(self.mesh.getNodesForPhysicalGroup(dim=1, tag=tag)[1])[:3]
1317 self.rc.isothermal_nodes.thin_shells[tag] = [float(coords[0]), float(coords[1]), float(coords[2])]
1319 elif elements == 'wedges':
1320 bare_mesh = {'1_1': self.rm.induced['r1_a1'].vol.numbers, '2_1': self.rm.induced['r2_a1'].vol.numbers,
1321 '1_2': self.rm.induced['r1_a2'].vol.numbers, '2_2': self.rm.induced['r2_a2'].vol.numbers}
1322 groups = self.rc.isothermal_nodes.wedges
1323 else:
1324 bare_mesh = {}
1325 groups = {}
1327 for group, tags_list in bare_mesh.items():
1328 groups[group] = {}
1329 for tag in tags_list:
1330 coords = list(self.mesh.getNodesForPhysicalGroup(dim=2, tag=tag)[1])[:3]
1331 groups[group][tag] = [float(coords[0]), float(coords[1]), float(coords[2])]
1333 for tag in self.rm.boundaries.thermal.cooling.groups['r' + group[0] + '_' + 'a' + group[-1]]:
1334 coords = list(self.mesh.getNodesForPhysicalGroup(dim=1, tag=tag)[1])[:3]
1335 groups[group][tag] = [float(coords[0]), float(coords[1]), float(coords[2])]
1337 logger.info(f"Info : {self.data.general.magnet_name} - E n d S e l e c t i n g M e s h N o d e s")
1338 # import time
1339 # time.sleep(100)
1340 # TODO: add to pro file automatically
1342 # self.occ.synchronize()
1343 # self.gu.launch_interactive_GUI()