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