Coverage for fiqus/post_processors/PostProcessCCT.py: 92%

290 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2024-12-27 02:39 +0100

1import os 

2import math 

3import csv 

4 

5import gmsh 

6import json 

7import numpy as np 

8import pandas as pd 

9from pathlib import Path 

10 

11from fiqus.geom_generators.GeometryCCT import Winding, FQPL 

12from fiqus.data.DataWindingsCCT import WindingsInformation 

13from fiqus.parsers.ParserDAT import ParserDAT 

14from fiqus.utils.Utils import GmshUtils, FilesAndFolders 

15 

16 

17class Post_Process: 

18 def __init__(self, fdm, verbose=True): 

19 """ 

20 Class to cct models postprocessing 

21 :param fdm: FiQuS data model 

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

23 """ 

24 self.cctdm = fdm.magnet 

25 self.model_folder = os.path.join(os.getcwd()) 

26 self.magnet_name = fdm.general.magnet_name 

27 

28 self.verbose = verbose 

29 self.gu = GmshUtils(self.model_folder, self.verbose) 

30 self.gu.initialize() 

31 self.masks_fqpls = {} 

32 self.pos_names = [] 

33 for variable, volume, file_ext in zip(self.cctdm.postproc.variables, self.cctdm.postproc.volumes, self.cctdm.postproc.file_exts): 

34 self.pos_names.append(f'{variable}_{volume}.{file_ext}') 

35 self.pos_name = self.pos_names[0] 

36 self.field_map_3D = os.path.join(self.model_folder, 'field_map_3D.csv') 

37 self.model_file = self.field_map_3D 

38 self.geom_folder = Path(self.model_folder).parent.parent 

39 # csv output definition for fields 

40 self.csv_ds = {'ds [m]': None} # length of each hexahedron 

41 self.csv_sAve = {'sAvePositions [m]': None} # cumulative positions along the electrical order of the center of each hexahedron 

42 self.csv_3Dsurf = { # basically this is the mapping of coordinates of corners to output. This dictionary defines which corner of hexahedra gets 

43 # used and which coordinate for each name like 'x3Dsurf 1'. Look into Hexahedron.py base class for drawing of corner numbering 

44 'x3Dsurf 1 [m]': {5: 'x'}, 

45 'x3Dsurf 2 [m]': {6: 'x'}, 

46 'y3Dsurf 1 [m]': {5: 'y'}, 

47 'y3Dsurf 2 [m]': {6: 'y'}, 

48 'z3Dsurf 1 [m]': {5: 'z'}, 

49 'z3Dsurf 2 [m]': {6: 'z'}, 

50 'x3Dsurf_2 1 [m]': {7: 'x'}, 

51 'x3Dsurf_2 2 [m]': {6: 'x'}, 

52 'y3Dsurf_2 1 [m]': {7: 'y'}, 

53 'y3Dsurf_2 2 [m]': {6: 'y'}, 

54 'z3Dsurf_2 1 [m]': {7: 'z'}, 

55 'z3Dsurf_2 2 [m]': {6: 'z'}} 

56 self.csv_nodeToTurns = {'nodeToHalfTurn [-]': None, # hexahedron (node) mapping to "Electrical" turns (electrical order) 

57 'nodeToPhysicalTurn [-]': None} # hexahedron (node) mapping to "Physical" turns (this is used for defining thermal connections between turns) 

58 self.csv_TransportCurrent = {'TransportCurrent [A]': None} # Transport current in the turn (not channel) for which the magnetic field was calculated. 

59 # The above does not need to be changed for running models at different current. 

60 self.csv_Bs = { 

61 # lists with magnetic field flux density components along cartesian coordinates (x, y, z) and l, h, w i.e. length, height and width of the channel (or wires if aligned with channel). 

62 # x, y, z is used for display, l, h, w are used in LEDET for Ic scaling (well only h and w for perpendicular field). 

63 'Bx [T]': 'Vx', 

64 'By [T]': 'Vy', 

65 'Bz [T]': 'Vz', 

66 'Bl [T]': 'Vl', 

67 'Bh [T]': 'Vh', 

68 'Bw [T]': 'Vw', 

69 } 

70 self.headers_dict = {**self.csv_ds, **self.csv_sAve, **self.csv_3Dsurf, **self.csv_nodeToTurns, **self.csv_TransportCurrent, **self.csv_Bs} 

71 self.data_for_csv = {} 

72 self.precision_for_csv = {} 

73 # -------- make keys in data for csv ---- 

74 for key in list(self.headers_dict.keys()): 

75 self.precision_for_csv[key] = '%.6f' # data will be written to csv output file with 6 decimal places, i.e. down to um and uT level. 

76 for key in list(self.csv_ds.keys()): 

77 self.data_for_csv[key] = [] 

78 for key in list(self.csv_sAve.keys()): 

79 self.data_for_csv[key] = [] 

80 for key in list(self.csv_3Dsurf.keys()): 

81 self.data_for_csv[key] = [] 

82 for key in list(self.csv_nodeToTurns.keys()): 

83 self.data_for_csv[key] = [] 

84 for key in list(self.csv_TransportCurrent.keys()): 

85 self.data_for_csv[key] = [] 

86 for key in list(self.csv_Bs.keys()): 

87 self.data_for_csv[key] = [] 

88 self.distance_key = list(self.csv_ds.keys())[0] 

89 self.sAve_key = list(self.csv_sAve.keys())[0] 

90 

91 @staticmethod 

92 def _get_fields(coord_center, normals): 

93 """ 

94 Helper funciton to probe magnetic field solution form gmsh view and calculate magnetic field along the height, width and length of the channel section (hex element) 

95 :param coord_center: dictionary with 'x', 'y' and 'z' coordinates stored as lists with items for each hexahedra. Magnetic field is probed at these locations. 

96 :param normals: dictionary with normals along the height, width and length stored per coordinate direction 'n_x', 'n_y' and 'n_z' and lists with items for each hexahedra. 

97 :return: tuple with list of Bxs, Bys, Bzs, Bhs, Bws, Bls 

98 """ 

99 Bxs = [] 

100 Bys = [] 

101 Bzs = [] 

102 Bhs = [] 

103 Bws = [] 

104 Bls = [] 

105 view_tag = gmsh.view.getTags()[0] 

106 for v, _ in enumerate(coord_center['x']): 

107 field = gmsh.view.probe(view_tag, coord_center['x'][v], coord_center['y'][v], coord_center['z'][v])[0] 

108 Bx = field[0] 

109 By = field[1] 

110 Bz = field[2] 

111 Bmod = math.sqrt(Bx ** 2 + By ** 2 + Bz ** 2) 

112 Bh = abs(Bx * normals['normals_h']['n_x'][v] + By * normals['normals_h']['n_y'][v] + Bz * normals['normals_h']['n_z'][v]) 

113 Bw = abs(Bx * normals['normals_w']['n_x'][v] + By * normals['normals_w']['n_y'][v] + Bz * normals['normals_w']['n_z'][v]) 

114 Bl = math.sqrt(abs(Bmod ** 2 - Bh ** 2 - Bw ** 2)) # sometimes it is very small and negative, so force it to be positive before taking the root. 

115 for arr, val in zip([Bxs, Bys, Bzs, Bhs, Bws, Bls], [Bx, By, Bz, Bh, Bw, Bl]): 

116 arr.append(val) 

117 return Bxs, Bys, Bzs, Bhs, Bws, Bls 

118 

119 @staticmethod 

120 def _plot_fields_in_views(p_name, global_channel_pos, coord_center, Bxs, Bys, Bzs): 

121 """ 

122 Add gmsh list data view with name 'B cartesian {p_name}, {global_channel_pos}' at coordinates stored in coord_center and magnetic field values from list Bxs, Bys, Bzs 

123 :param p_name: string with powered region name 

124 :param global_channel_pos: integer with wire position in the channel 

125 :param coord_center: dictionary with 'x', 'y' and 'z' coordinates stored as lists with items for each hexahedra. Magnetic field is plotted at these locations. 

126 :param Bxs: list of magnetic field along the x axis to use in the view 

127 :param Bys: list of magnetic field along the y axis to use in the view 

128 :param Bzs: list of magnetic field along the z axis to use in the view 

129 :return: none, adds view to currently initialized gmsh model and synchronizes 

130 """ 

131 data_cartesian = [] 

132 for v, _ in enumerate(coord_center['x']): 

133 data_cartesian.append(coord_center['x'][v]) 

134 data_cartesian.append(coord_center['y'][v]) 

135 data_cartesian.append(coord_center['z'][v]) 

136 data_cartesian.append(Bxs[v]) 

137 data_cartesian.append(Bys[v]) 

138 data_cartesian.append(Bzs[v]) 

139 gmsh.view.addListData(gmsh.view.add(f'B cartesian {p_name}, {global_channel_pos}'), "VP", len(data_cartesian) // 6, data_cartesian) 

140 gmsh.model.occ.synchronize() 

141 

142 @staticmethod 

143 def _plot_turns_in_view(data_turns, coord_center, turns_for_ch_pos): 

144 """ 

145 Add gmsh list data view with name 'Turn numbering' at coordinates stored in coord_center and values from turns_for_ch_pos 

146 :param data_turns: This is list to which the data gets appended so at the end one view with all the turns labeled is created with looping through turns in the channels. 

147 :param coord_center: dictionary with 'x', 'y' and 'z' coordinates stored as lists with items for each hexahedra. Magnetic field is plotted at these locations. 

148 :param turns_for_ch_pos: Lists with turn numbers to plot in the view. 

149 :return: none, adds view to currently initialized gmsh model and synchronizes 

150 """ 

151 for v, _ in enumerate(coord_center['x']): 

152 data_turns.append(coord_center['x'][v]) 

153 data_turns.append(coord_center['y'][v]) 

154 data_turns.append(coord_center['z'][v]) 

155 data_turns.append(turns_for_ch_pos[v]) 

156 gmsh.view.addListData(gmsh.view.add('Turn numbering'), "SP", len(data_turns) // 4, data_turns) 

157 gmsh.model.occ.synchronize() 

158 

159 @staticmethod 

160 def _check_normals_match_corners(coord_center, normals_dict, p_name): 

161 """ 

162 This is error checking function. Coord_centers are calculated in this class 'on the 'fly' form cctdm file. Normals are generated form geom_generators before meshing, then saved to file and now loaded here. 

163 There is a chance that thses could get 'out of sync'. A basic check of the length (i.e. number of hexahedra involved is performed) of these lists is performed. 

164 This will only detect if number of turns or discretization of each powered geom_generators has changed. 

165 :param coord_center: list with coordinate centers. The content is not used only list length is queried. 

166 :param normals_dict: list with normals read form json file. The content is not used only list length is queried. 

167 :param p_name: string with powered volume name. Only used for error message display to say which region is not matching. 

168 :return: None, prints error to the python console, if any. 

169 """ 

170 for cord in ['x', 'y', 'z']: 

171 if coord_center[cord].size != len(normals_dict['normals_h'][cord]) or coord_center[cord].size != len(normals_dict['normals_w'][cord]): 

172 raise ValueError(f'Number of volumes in normals does not match the magnet definition for {p_name} winding!') 

173 

174 @staticmethod 

175 def _calc_coord_center_and_ds(corners_lists): 

176 """ 

177 Calculates a geometrical centre of each hexahedron base on its corners coordinates 

178 :param corners_lists: lists of lists corresponding to a list of hexahedrons with each list containing 8 corner coordinates for each of the 8 points in the corners of hexahedron 

179 :return: coord_center: dictionary with 'x', 'y' and 'z' coordinates stored as lists with items for each hexahedra. and ds 

180 """ 

181 coord_center = {} 

182 for cord in ['x', 'y', 'z']: 

183 coord_center[cord] = np.mean(corners_lists[cord], axis=0) # coordinates center 

184 surf_close = {} 

185 surf_far = {} 

186 for cord in ['x', 'y', 'z']: 

187 surf_close[cord] = np.mean(corners_lists[cord][0:4], axis=0) # coordinates surface close 

188 surf_far[cord] = np.mean(corners_lists[cord][4:8], axis=0) # coordinates surface far 

189 ds = np.sqrt(np.square(surf_close['x'] - surf_far['x']) + np.square(surf_close['y'] - surf_far['y']) + np.square(surf_close['z'] - surf_far['z'])) 

190 return coord_center, ds 

191 

192 def _load_to_fields_csv_ds(self, ds, h3Dsurf, nodeToHalfTurn, nodeToPhysicalTurn, current, Bxs, Bys, Bzs, Bls, Bhs, Bws, s_ave): 

193 """ 

194 Method to load data to csv dictionary to be dumped to a csv file. 

195 :param ds: length of each hexahedron 

196 :param h3Dsurf: dict with coordinates to export 

197 :param nodeToHalfTurn: hexahedron (node) mapping to "Electrical" turns (electrical order) 

198 :param nodeToPhysicalTurn: hexahedron (node) mapping to "Physical" turns (this is used for defining thermal connections between turns) 

199 :param current: Transport current in the turn (not channel) for which the magnetic field was calculated. It does not need to be changed for running models at different current. 

200 :param Bxs: list of magnetic flux density component along x direction 

201 :param Bys: list of magnetic flux density component along y direction 

202 :param Bzs: list of magnetic flux density component along z direction 

203 :param Bls: list of magnetic flux density component along l (length) of wire direction 

204 :param Bhs: list of magnetic flux density component along h (height) of wire direction 

205 :param Bws: list of magnetic flux density component along w (wight) of wire direction 

206 :param s_ave: cumulative positions along the electrical order of the center of each hexahedron 

207 :return: nothing, appends data to this class attribute data_for_csv 

208 """ 

209 for arr, key in zip([ds, s_ave], list(self.csv_ds.keys()) + list(self.csv_sAve.keys())): 

210 self.data_for_csv[key].extend(arr) 

211 for key, dict_what in self.csv_3Dsurf.items(): 

212 for p_n, p_c in dict_what.items(): 

213 self.data_for_csv[key].extend(h3Dsurf[p_c][p_n]) 

214 for arr, key in zip([Bxs, Bys, Bzs, Bls, Bhs, Bws, current, nodeToHalfTurn, nodeToPhysicalTurn], list(self.csv_Bs.keys()) + list(self.csv_TransportCurrent.keys()) + list(self.csv_nodeToTurns.keys())): 

215 self.data_for_csv[key].extend(arr) 

216 

217 def _save_fields_csv(self): 

218 """ 

219 Helper method to save csv file from data_for_csv class attribute. 

220 :return: nothing, saves csv to file with name magnet_name(from yaml input file)_total_n_turns.csv 

221 """ 

222 # csv_file_path = os.path.join(self.model_folder, f'{self.magnet_name[:self.magnet_name.index("_")]}_{int(total_n_turns)}.csv') 

223 

224 df = pd.DataFrame(self.data_for_csv) 

225 for column in df: 

226 df[column] = df[column].map(lambda x: self.precision_for_csv[column] % x) 

227 # put NaN at the end of some columns to match what LEDET needs. 

228 df.loc[df.index[-1], self.distance_key] = 'NaN' 

229 df.loc[df.index[-1], self.sAve_key] = 'NaN' 

230 # df.loc[df.index[-1], list(self.csv_nodeToTurns.keys())[0]] = 'NaN' 

231 for h in list(self.csv_nodeToTurns.keys()): 

232 df.loc[df.index[-1], h] = 'NaN' 

233 df.loc[df.index[-1], list(self.csv_TransportCurrent.keys())[0]] = 'NaN' 

234 for h in list(self.csv_Bs.keys()): 

235 df.loc[df.index[-1], h] = 'NaN' 

236 df.to_csv(self.field_map_3D, index=False, sep=',') # , float_format='%.7f') 

237 

238 @staticmethod 

239 def _trim_array(array, mask): 

240 return list(np.array(array)[mask]) 

241 

242 @staticmethod 

243 def _all_equal(iterator): 

244 iterator = iter(iterator) 

245 try: 

246 first = next(iterator) 

247 except StopIteration: 

248 return True 

249 return all(first == x for x in iterator) 

250 

251 def postporcess_inductance(self): 

252 """ 

253 Function to postprocess .dat file with inductance saved by GetDP to selfMutualInductanceMatrix.csv required by LEDET for CCT magnet solve 

254 :return: Nothing, writes selfMutualInductanceMatrix.csv file on disc. 

255 :rtype: None 

256 """ 

257 inductance_file_from_getdp = os.path.join(self.model_folder, 'Inductance.dat') 

258 channels_inductance = ParserDAT(inductance_file_from_getdp).pqv # postprocessed quantity value (pqv) for Inductance.dat contains magnet self-inductance (channel, not wire turns) 

259 

260 total_n_turns_channel = 0 # including fqpls 

261 for n_turns in self.cctdm.geometry.windings.n_turnss: 

262 total_n_turns_channel = total_n_turns_channel + n_turns 

263 if self._all_equal(self.cctdm.postproc.windings_wwns) and self._all_equal(self.cctdm.postproc.windings_whns): 

264 number_of_turns_in_channel = self.cctdm.postproc.windings_wwns[0] * self.cctdm.postproc.windings_whns[0] 

265 total_n_turns_channel = int(total_n_turns_channel) 

266 total_n_turns_channel_and_fqpls = total_n_turns_channel 

267 

268 #windings_inductance = channels_inductance * total_n_turns_windings**2 

269 windings_inductance = channels_inductance * number_of_turns_in_channel**3 # scaling with square for number of turns in the channel as the model has current in the channel but not number of turns 

270 

271 

272 geometry_folder = Path(self.model_folder).parent.parent 

273 winding_information_file = os.path.join(geometry_folder, f"{self.magnet_name}.wi") 

274 cctwi = FilesAndFolders.read_data_from_yaml(winding_information_file, WindingsInformation) 

275 windings_avg_length = cctwi.windings_avg_length 

276 windings_inductance_per_m = windings_inductance / windings_avg_length 

277 

278 print(f"Channels self-inductance: {channels_inductance} H") 

279 print(f"Number of turns in channel: {number_of_turns_in_channel} H") 

280 print(f"Total number of channel turns: {total_n_turns_channel} turns") 

281 print(f"Total number of inductance blocks: {total_n_turns_channel}") 

282 print(f"Windings self-inductance: {windings_inductance} H") 

283 print(f"Windings self-inductance per meter: {windings_inductance_per_m} H/m") 

284 print(f"Magnetic length: {windings_avg_length} m") 

285 

286 fqpls_dummy_inductance_block = 1e-9 # H/m (H per meter of magnet!) 

287 for _ in self.cctdm.geometry.fqpls.names: 

288 total_n_turns_channel_and_fqpls = total_n_turns_channel_and_fqpls + 1 

289 windings_inductance_per_m = windings_inductance_per_m - fqpls_dummy_inductance_block # subtracting fqpls as they will be added later to the matrix 

290 

291 win_ind_for_matrix = windings_inductance_per_m / total_n_turns_channel**2 

292 M_matrix = np.ones(shape=(total_n_turns_channel_and_fqpls, total_n_turns_channel_and_fqpls)) * win_ind_for_matrix 

293 M_matrix[:, total_n_turns_channel:total_n_turns_channel_and_fqpls] = 0.0 # mutual inductance of windings to fqpl set to zero for the columns 

294 M_matrix[total_n_turns_channel:total_n_turns_channel_and_fqpls, :] = 0.0 # mutual inductance of windings to fqpl set to zero for the rows 

295 for ij in range(total_n_turns_channel, total_n_turns_channel_and_fqpls, 1): 

296 M_matrix[ij, ij] = fqpls_dummy_inductance_block # self-inductance of fqpls inductance block 

297 # below code is the same as in BuilderLEDET in steam_sdk 

298 csv_write_path = os.path.join(self.model_folder, 'selfMutualInductanceMatrix.csv') 

299 print(f'Saving square matrix of size {total_n_turns_channel_and_fqpls}x{total_n_turns_channel_and_fqpls} into: {csv_write_path}') 

300 with open(csv_write_path, 'w', newline='') as file: 

301 reader = csv.writer(file) 

302 reader.writerow(["# Extended self mutual inductance matrix [H/m]"]) 

303 for i in range(M_matrix.shape[0]): 

304 reader.writerow(M_matrix[i]) 

305 

306 def postprocess_fields(self, gui=False): 

307 """ 

308 Methods to calculated 'virtual' turn positions in the channels and output information for about them for LEDET as a csv file. 

309 :param gui: if True, graphical user interface (gui) of gmsh is shown with views created 

310 :return: nothing, writes csv to file 

311 """ 

312 winding_order = self.cctdm.postproc.winding_order 

313 total_n_turns = 0 

314 for n_turns, wwn, whn in zip(self.cctdm.geometry.windings.n_turnss, self.cctdm.postproc.windings_wwns, self.cctdm.postproc.windings_whns): 

315 total_n_turns = total_n_turns + n_turns * wwn * whn 

316 

317 gmsh.open(self.pos_name) 

318 global_channel_pos = 1 

319 # data_turns = [] # used in self._plot_turns_in_view() - do not delete 

320 csv_data_dict = {} 

321 s_max = 0 

322 for w_i, w_name in enumerate(self.cctdm.geometry.windings.names): # + self.cctwi.f_names: 

323 normals_dict = json.load(open(os.path.join(self.geom_folder, f'{w_name}.normals'))) 

324 winding_obj = Winding(self.cctdm, w_i, post_proc=True) 

325 ww2 = self.cctdm.geometry.windings.wwws[w_i] / 2 # half of channel width 

326 wh2 = self.cctdm.geometry.windings.wwhs[w_i] / 2 # half of channel height 

327 wwns = self.cctdm.postproc.windings_wwns[w_i] # number of wires in width direction 

328 whns = self.cctdm.postproc.windings_whns[w_i] # number of wires in height direction 

329 wsw = self.cctdm.geometry.windings.wwws[w_i] / wwns # wire size in width 

330 wsh = self.cctdm.geometry.windings.wwhs[w_i] / whns # wire size in height 

331 for i_w in range(wwns): 

332 for i_h in range(whns): 

333 corner_i = 0 

334 corners_lists = {'x': [], 'y': [], 'z': []} 

335 for far in [False, True]: 

336 for ii_w, ii_h in zip([0, 0, 1, 1], [0, 1, 1, 0]): 

337 wt = (-ww2 + (ii_w + i_w) * wsw) / ww2 

338 ht = (-wh2 + (ii_h + i_h) * wsh) / wh2 

339 corner_i += 1 

340 xs, ys, zs, turns_from_rotation = winding_obj.calc_turns_corner_coords(wt, ht, far=far) 

341 zs = [z - winding_obj.z_corr for z in zs] 

342 corners_lists['x'].append(xs) 

343 corners_lists['y'].append(ys) 

344 corners_lists['z'].append(zs) 

345 nodeToPhysicalTurn = [(global_channel_pos - 1) * self.cctdm.geometry.windings.n_turnss[w_i] + t for t in turns_from_rotation] 

346 nodeToHalfTurn = [winding_order.index(global_channel_pos) * self.cctdm.geometry.windings.n_turnss[w_i] + t for t in turns_from_rotation] 

347 coord_center, ds = self._calc_coord_center_and_ds(corners_lists) 

348 self._check_normals_match_corners(coord_center, normals_dict, w_name) 

349 # if global_channel_pos == 1: 

350 # s_ave_elem_len = np.copy(ds) 

351 # s_ave_elem_len[0] = s_ave_elem_len[0] / 2 

352 # s_ave = np.cumsum(s_ave_elem_len) 

353 # else: 

354 # s_ave = np.cumsum(ds) + s_max 

355 # s_max = np.max(s_ave) 

356 current = [abs(self.cctdm.solve.windings.currents[w_i]) / (wwns * whns)] * len(nodeToHalfTurn) 

357 Bxs, Bys, Bzs, Bhs, Bws, Bls = self._get_fields(coord_center, normals_dict) 

358 if gui: 

359 self._plot_fields_in_views(w_name, global_channel_pos, coord_center, Bxs, Bys, Bzs) 

360 # self._plot_turns_in_view(data_turns, coord_center, nodeToHalfTurn) 

361 csv_data_dict[winding_order.index(global_channel_pos) + 1] = ds, corners_lists, nodeToHalfTurn, nodeToPhysicalTurn, current, Bxs, Bys, Bzs, Bls, Bhs, Bws 

362 global_channel_pos += 1 

363 for f_i, f_name in enumerate(self.cctdm.geometry.fqpls.names): 

364 # ----- calculate centers of coordinates for hexes of fqpl 

365 fqpl_obj = FQPL(self.cctdm, f_i) 

366 corners_lists = {'x': [], 'y': [], 'z': []} 

367 for corner_i in range(8): 

368 for cord in ['x', 'y', 'z']: 

369 corners_lists[cord].append([v_dict[corner_i][cord] for v_dict in fqpl_obj.hexes.values()]) 

370 coord_center, ds = self._calc_coord_center_and_ds(corners_lists) 

371 normals_dict = json.load(open(os.path.join(self.geom_folder, f'{f_name}.normals'))) 

372 self._check_normals_match_corners(coord_center, normals_dict, f_name) 

373 Bxs, Bys, Bzs, Bhs, Bws, Bls = self._get_fields(coord_center, normals_dict) 

374 self.masks_fqpls[f_i] = [] # dictionary to keep masks for trimming one end of fqpls that is long to extend to the end of air region, but not needed that long in simulations 

375 for z_close, z_far in zip(corners_lists['z'][0], corners_lists['z'][-1]): 

376 if z_close > -self.cctdm.geometry.fqpls.z_ends[f_i] * self.cctdm.postproc.fqpl_export_trim_tol[f_i] and z_far > -self.cctdm.geometry.fqpls.z_ends[f_i] * self.cctdm.postproc.fqpl_export_trim_tol[f_i]: 

377 self.masks_fqpls[f_i].append(True) 

378 else: 

379 self.masks_fqpls[f_i].append(False) 

380 self.masks_fqpls[f_i] = np.array(self.masks_fqpls[f_i], dtype=bool) 

381 for cord in ['x', 'y', 'z']: 

382 for corner_i, corner in enumerate(corners_lists[cord]): 

383 corners_lists[cord][corner_i] = list(np.array(corner)[self.masks_fqpls[f_i]]) 

384 lists_to_trim = [Bxs, Bys, Bzs, Bhs, Bws, Bls, ds] 

385 for idx, array in enumerate(lists_to_trim): 

386 lists_to_trim[idx] = self._trim_array(array, self.masks_fqpls[f_i]) 

387 [Bxs, Bys, Bzs, Bhs, Bws, Bls, ds] = lists_to_trim 

388 for cord in ['x', 'y', 'z']: 

389 coord_center[cord] = self._trim_array(coord_center[cord], self.masks_fqpls[f_i]) 

390 # s_ave = np.cumsum(ds) + s_max 

391 # s_max = np.max(s_ave) 

392 go_dict = {'from': 0, 'to': len(Bxs)//2} # fqpl 'go' part 

393 ret_dict = {'from': len(Bxs)//2, 'to': len(Bxs)} 

394 for go_ret in [go_dict, ret_dict]: 

395 total_n_turns += 1 

396 ds_half = ds[go_ret['from']:go_ret['to']] 

397 corners_lists_half = {} 

398 for cord in ['x', 'y', 'z']: 

399 corners_lists_half[cord] = [] 

400 for corner_i, corner in enumerate(corners_lists[cord]): 

401 corners_lists_half[cord].append(corner[go_ret['from']:go_ret['to']]) 

402 nodeToHalfTurn = [total_n_turns] * len(Bxs[go_ret['from']:go_ret['to']]) # total expected number of turns + one + fqpl number 

403 nodeToPhysicalTurn = [total_n_turns] * len(Bxs[go_ret['from']:go_ret['to']]) 

404 current = [self.cctdm.solve.fqpls.currents[f_i]] * len(Bxs[go_ret['from']:go_ret['to']]) 

405 Bxs_half = Bxs[go_ret['from']:go_ret['to']] 

406 Bys_half = Bys[go_ret['from']:go_ret['to']] 

407 Bzs_half = Bzs[go_ret['from']:go_ret['to']] 

408 Bls_half = Bls[go_ret['from']:go_ret['to']] 

409 Bhs_half = Bhs[go_ret['from']:go_ret['to']] 

410 Bws_half = Bws[go_ret['from']:go_ret['to']] 

411 #s_ave_half = s_ave[go_ret['from']:go_ret['to']] 

412 csv_data_dict[total_n_turns] = ds_half, corners_lists_half, nodeToHalfTurn, nodeToPhysicalTurn, current, Bxs_half, Bys_half, Bzs_half, Bls_half, Bhs_half, Bws_half 

413 if gui: 

414 self._plot_fields_in_views(f_name, global_channel_pos, coord_center, Bxs, Bys, Bzs) 

415 global_channel_pos += 1 

416 

417 ds_global = [] 

418 csv_data_sorded = dict(sorted(csv_data_dict.items())).values()# sorted list of tuples with the csv data 

419 for tuple_with_values in csv_data_sorded: 

420 ds_global.append(tuple_with_values[0]) # 0 index in tuple is for ds, changing to list to easily extend. 

421 

422 s_ave_global = [] 

423 s_max = -ds_global[0][0]/2 # make a s_max to a negative half of the length of the first hexahedra. This is to set the first s_ave max to a half of the first ds element 

424 for ds_i, ds_np in enumerate(ds_global): 

425 s_ave_np = np.cumsum(ds_np) + s_max 

426 s_ave_global.append(s_ave_np) 

427 s_max = np.max(s_ave_np) 

428 

429 for sorted_v_tuple, s_ave_np in zip(csv_data_sorded, s_ave_global): 

430 self._load_to_fields_csv_ds(*sorted_v_tuple, s_ave_np) # this is to avoid rearranging s_ave by electrical order 

431 

432 #el_order_half_turns = [int(t) for t in list(pd.unique(self.data_for_csv['nodeToPhysicalTurn [-]']))] 

433 #json.dump({'el_order_half_turns': el_order_half_turns}, open(os.path.join(self.model_folder, "el_order_half_turns.json"), 'w')) 

434 self._save_fields_csv() 

435 if gui: 

436 self.gu.launch_interactive_GUI() 

437 else: 

438 gmsh.clear() 

439 gmsh.finalize() 

440