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

1import os 

2import gmsh 

3from numpy import square as square 

4from numpy import sqrt as sqrt 

5import json 

6 

7from fiqus.utils.Utils import GmshUtils 

8from fiqus.utils.Utils import FilesAndFolders as Util 

9from fiqus.utils.Utils import GeometricFunctions as Func 

10from fiqus.data import DataFiQuS as dF 

11from fiqus.data import DataMultipole as dM 

12 

13 

14class Geometry: 

15 def __init__(self, data: dF.FDM() = None, geom: dF.FiQuSGeometry() = None, sett: dF.FiQuSSettings() = None, 

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

17 """ 

18 Class to generate geometry 

19 :param data: FiQuS data model 

20 :param geom: ROXIE geometry data 

21 :param sett: settings data model 

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

23 """ 

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

25 self.geom: dF.FiQuSGeometry() = geom 

26 self.sett: dF.FiQuSSettings() = sett 

27 self.geom_folder = geom_folder 

28 self.verbose: bool = verbose 

29 

30 self.md = dM.MultipoleData() 

31 

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

33 self.gu.initialize() 

34 self.occ = gmsh.model.occ 

35 

36 self.model_file = f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.brep" 

37 self.iL, self.iR, self.oL, self.oR, self.iLr, self.iRr, self.oLr, self.oRr = [], [], [], [], [], [], [], [] 

38 

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

40 if gui: 

41 self.gu.launch_interactive_GUI() 

42 else: 

43 gmsh.clear() 

44 gmsh.finalize() 

45 

46 def saveHalfTurnCornerPositions(self): 

47 json.dump({'iL': self.iL, 'iR': self.iR, 'oL': self.oL, 'oR': self.oR, 

48 'iLr': self.iLr, 'iRr': self.iRr, 'oLr': self.oLr, 'oRr': self.oRr}, 

49 open(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.crns", 'w')) 

50 

51 def saveStrandPositions(self): 

52 ht_nr = 0 

53 parser_x, parser_y, blocks, ht = [], [], [], [] 

54 for eo in self.geom.Roxie_Data.coil.electrical_order: 

55 block = self.geom.Roxie_Data.coil.coils[eo.coil].poles[eo.pole].layers[eo.layer].windings[ 

56 eo.winding].blocks[eo.block] 

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

58 ht_nr += 1 

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

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

61 blocks.append(eo.block) 

62 ht.append(ht_nr) 

63 parser_x.append(strand.x) 

64 parser_y.append(strand.y) 

65 json.dump({'x': parser_x, 'y': parser_y, 'block': blocks, 'ht': ht}, 

66 open(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.strs", 'w')) 

67 

68 def saveBoundaryRepresentationFile(self): 

69 self.occ.synchronize() 

70 gmsh.write(self.model_file) 

71 gmsh.clear() 

72 

73 def loadBoundaryRepresentationFile(self): 

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

75 gmsh.open(self.model_file) 

76 

77 def saveAuxiliaryFile(self): 

78 Util.write_data_to_yaml(f"{os.path.join(self.geom_folder, self.data.general.magnet_name)}.aux", self.md.dict()) 

79 

80 def constructIronGeometry(self): 

81 """ 

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

83 """ 

84 iron = self.geom.Roxie_Data.iron 

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

86 quadrants = self.md.geometries.iron.quadrants 

87 

88 lc = 1e-2 

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

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

91 

92 if point.x == 0.: 

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

94 else: 

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

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

97 if point.y == 0.: 

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

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

100 else: 

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

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

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

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

105 

106 mirror_x = [1, -1, -1, 1] 

107 mirror_y = [1, 1, -1, -1] 

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

109 if line.type == 'line': 

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

111 if quadrant == 1: 

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

113 elif quadrant == 2: 

114 if iron.key_points[line.kp1].x == 0. and iron.key_points[line.kp2].x == 0.: 

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

116 else: 

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

118 elif quadrant == 3: 

119 if iron.key_points[line.kp1].y == 0. and iron.key_points[line.kp2].y == 0.: 

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

121 else: 

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

123 else: 

124 if iron.key_points[line.kp1].y == 0. and iron.key_points[line.kp2].y == 0.: 

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

126 else: 

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

128 

129 elif line.type == 'arc': 

130 center = Func.arc_center_from_3_points( 

131 iron.key_points[line.kp1], iron.key_points[line.kp3], iron.key_points[line.kp2]) 

132 new_point_name = 'kp' + line_name + '_center' 

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

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

135 if center[0] == 0.: 

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

137 else: 

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

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

140 if center[1] == 0.: 

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

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

143 else: 

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

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

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

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

148 

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

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

151 qq.points[line.kp1], qq.points[new_point_name], qq.points[line.kp2]) 

152 

153 elif line.type == 'circle': 

154 pt1 = iron.key_points[line.kp1] 

155 pt2 = iron.key_points[line.kp2] 

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

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

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

159 

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

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

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

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

164 

165 else: 

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

167 

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

169 for area_name, area in iron.hyper_areas.items(): 

170 qq.areas[area_name] = dM.Area(loop=self.occ.addCurveLoop([qq.lines[line] for line in area.lines])) 

171 if (iron.hyper_areas[area_name].material not in self.md.domains.groups_surfaces and 

172 iron.hyper_areas[area_name].material != 'BH_air'): 

173 self.md.domains.groups_surfaces[iron.hyper_areas[area_name].material] = [] 

174 

175 def constructWedgeGeometry(self): 

176 """ 

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

178 """ 

179 wedges = self.geom.Roxie_Data.wedges 

180 for i in wedges: 

181 kp0_inner = self.occ.addPoint(wedges[i].corrected_center.inner.x, wedges[i].corrected_center.inner.y, 0) 

182 kp0_outer = self.occ.addPoint(wedges[i].corrected_center.outer.x, wedges[i].corrected_center.outer.y, 0) 

183 arg = [self.occ.addPoint(wedges[i].corners.iL.x, wedges[i].corners.iL.y, 0), 

184 self.occ.addPoint(wedges[i].corners.iR.x, wedges[i].corners.iR.y, 0), 

185 self.occ.addPoint(wedges[i].corners.oL.x, wedges[i].corners.oL.y, 0), 

186 self.occ.addPoint(wedges[i].corners.oR.x, wedges[i].corners.oR.y, 0)] 

187 

188 left = self.occ.addLine(arg[0], arg[2]) 

189 right = self.occ.addLine(arg[1], arg[3]) 

190 inner = self.occ.addCircleArc(arg[0], kp0_inner, arg[1]) 

191 outer = self.occ.addCircleArc(arg[2], kp0_outer, arg[3]) 

192 self.md.geometries.wedges.areas[i] = dM.Area(loop=self.occ.addCurveLoop([inner, right, outer, left])) 

193 

194 def constructCoilGeometry(self): 

195 """ 

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

197 """ 

198 frac = 1 # self.frac if self.frac else 1 

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

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

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

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

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

204 poles = coils.poles[pole_nr] 

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

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

207 layers = poles.layers[layer_nr] 

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

209 layers.windings[winding_nr] = dM.Block(conductor_name=winding.conductor_name, 

210 conductors_number=winding.conductors_number) 

211 windings = layers.windings[winding_nr] 

212 ll = self.sett.Model_Data_GS.conductors[winding.conductor_name].cable.bare_cable_height_mean 

213 for block_key, block in winding.blocks.items(): 

214 windings.blocks[block_key] = dM.BlockData(current_sign=block.current_sign) 

215 hts = windings.blocks[block_key].half_turns 

216 

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

218 ht = str(halfTurn_nr) 

219 hts.areas[ht] = dM.Area() 

220 hf_current = halfTurn.corners.insulated 

221 if halfTurn_nr == 1: 

222 bc = block.block_corners 

223 if (hf_current.iL.y < 0 and block_key == max(winding.blocks.keys())) or \ 

224 (hf_current.iL.y > 0 and block_key == min(winding.blocks.keys())): 

225 hts.points[ht + 'i'] = self.occ.addPoint(bc.iR.x, bc.iR.y, 0, ll / frac) 

226 hts.points[ht + 'o'] = self.occ.addPoint(bc.oR.x, bc.oR.y, 0, ll / frac) 

227 else: 

228 hts.points[ht + 'i'] = self.occ.addPoint(bc.iL.x, bc.iL.y, 0, ll / frac) 

229 hts.points[ht + 'o'] = self.occ.addPoint(bc.oL.x, bc.oL.y, 0, ll / frac) 

230 hts.lines[ht + 'r'] = self.occ.addLine(hts.points[ht + 'i'], hts.points[ht + 'o']) 

231 

232 if halfTurn_nr == winding.conductors_number: 

233 bc = block.block_corners 

234 

235 if (hf_current.iL.y < 0 and block_key == max(winding.blocks.keys())) or \ 

236 (hf_current.iL.y > 0 and block_key == min(winding.blocks.keys())): 

237 hts.points['end_i'] = self.occ.addPoint(bc.iL.x, bc.iL.y, 0, ll/frac) 

238 hts.points['end_o'] = self.occ.addPoint(bc.oL.x, bc.oL.y, 0, ll/frac) 

239 else: 

240 hts.points['end_i'] = self.occ.addPoint(bc.iR.x, bc.iR.y, 0, ll/frac) 

241 hts.points['end_o'] = self.occ.addPoint(bc.oR.x, bc.oR.y, 0, ll/frac) 

242 

243 hts.lines[ht + 'i'] = self.occ.addLine(hts.points['end_i'], hts.points[ht + 'i']) 

244 hts.lines[ht + 'o'] = self.occ.addLine(hts.points['end_o'], hts.points[ht + 'o']) 

245 hts.lines['end'] = self.occ.addLine(hts.points['end_i'], hts.points['end_o']) 

246 

247 # For plotting only 

248 if pole_nr == 1: 

249 self.occ.synchronize() 

250 self.iL.append(list(gmsh.model.getValue(0, hts.points['end_i'], [])[:-1])) 

251 self.oL.append(list(gmsh.model.getValue(0, hts.points['end_o'], [])[:-1])) 

252 # if block_key == 1 or block_key == 3 or block_key == 5 or block_key == 6: 

253 # radius = np.sqrt(np.square(iL[0] - coil.bore_center.x) + 

254 # np.square(iL[1] - coil.bore_center.y)) 

255 # self.ax.add_patch( 

256 # patches.Circle((coil.bore_center.x, coil.bore_center.y), radius=radius, 

257 # color='b', fill=False)) 

258 # radius = np.sqrt(np.square(oL[0] - coil.bore_center.x) + 

259 # np.square(oL[1] - coil.bore_center.y)) 

260 # self.ax.add_patch( 

261 # patches.Circle((coil.bore_center.x, coil.bore_center.y), radius=radius, 

262 # color='c', fill=False)) 

263 ################ 

264 

265 left = hts.lines['end'] 

266 else: 

267 hf_next = block.half_turns[halfTurn_nr + 1].corners.insulated 

268 next_ht = str(halfTurn_nr + 1) 

269 mid_point_i = [(hf_next.iR.x + hf_current.iL.x) / 2, 

270 (hf_next.iR.y + hf_current.iL.y) / 2] 

271 mid_point_o = [(hf_next.oR.x + hf_current.oL.x) / 2, 

272 (hf_next.oR.y + hf_current.oL.y) / 2] 

273 hts.points[next_ht + 'i'] = self.occ.addPoint(mid_point_i[0], mid_point_i[1], 0, 

274 ll / frac) 

275 hts.points[next_ht + 'o'] = self.occ.addPoint(mid_point_o[0], mid_point_o[1], 0, 

276 ll / frac) 

277 

278 hts.lines[ht + 'i'] = \ 

279 self.occ.addLine(hts.points[next_ht + 'i'], hts.points[ht + 'i']) 

280 hts.lines[ht + 'o'] = \ 

281 self.occ.addLine(hts.points[next_ht + 'o'], hts.points[ht + 'o']) 

282 hts.lines[next_ht + 'r'] = \ 

283 self.occ.addLine(hts.points[next_ht + 'i'], hts.points[next_ht + 'o']) 

284 left = hts.lines[next_ht + 'r'] 

285 

286 # For plotting only 

287 if pole_nr == 1: 

288 self.occ.synchronize() 

289 self.iL.append(list(gmsh.model.getValue(0, hts.points[next_ht + 'i'], [])[:-1])) 

290 self.oL.append(list(gmsh.model.getValue(0, hts.points[next_ht + 'o'], [])[:-1])) 

291 if pole_nr == 1: 

292 self.iR.append(list(gmsh.model.getValue(0, hts.points[ht + 'i'], [])[:-1])) 

293 self.oR.append(list(gmsh.model.getValue(0, hts.points[ht + 'o'], [])[:-1])) 

294 hti = halfTurn.corners.insulated 

295 self.iLr.append([hti.iL.x, hti.iL.y]) 

296 self.iRr.append([hti.iR.x, hti.iR.y]) 

297 self.oLr.append([hti.oL.x, hti.oL.y]) 

298 self.oRr.append([hti.oR.x, hti.oR.y]) 

299 ################ 

300 

301 hts.areas[ht].loop = self.occ.addCurveLoop([hts.lines[ht + 'i'], # inner 

302 hts.lines[ht + 'r'], # right 

303 hts.lines[ht + 'o'], # outer 

304 left]) # left 

305 

306 self.saveHalfTurnCornerPositions() 

307 

308 def buildDomains(self): 

309 """ 

310 Generates plane surfaces from the curve loops 

311 """ 

312 iron = self.geom.Roxie_Data.iron 

313 gm = self.md.geometries 

314 

315 for i in iron.key_points: 

316 gm.iron.max_radius = max(gm.iron.max_radius, max(iron.key_points[i].x, iron.key_points[i].y)) 

317 

318 # Create and build air far field 

319 radius_in = gm.iron.max_radius * (2.5 if self.data.magnet.geometry.with_iron_yoke else 6) 

320 gm.air_inf.lines['inner'] = self.occ.addCircle(0., 0., 0., radius_in) 

321 gm.air_inf.areas['inner'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['inner']])) 

322 gm.air_inf.areas['inner'].surface = self.occ.addPlaneSurface([self.md.geometries.air_inf.areas['inner'].loop]) 

323 radius_out = gm.iron.max_radius * (3.2 if self.data.magnet.geometry.with_iron_yoke else 8) 

324 gm.air_inf.lines['outer'] = self.occ.addCircle(0., 0., 0., radius_out) 

325 gm.air_inf.areas['outer'] = dM.Area(loop=self.occ.addCurveLoop([gm.air_inf.lines['outer']])) 

326 gm.air_inf.areas['outer'].surface = self.occ.addPlaneSurface([gm.air_inf.areas['outer'].loop, 

327 gm.air_inf.areas['inner'].loop]) 

328 self.md.domains.groups_surfaces['air_inf'] = [gm.air_inf.areas['outer'].surface] 

329 

330 # Build iron yoke domains 

331 if self.data.magnet.geometry.with_iron_yoke: 

332 for quadrant, qq in gm.iron.quadrants.items(): 

333 for area_name, area in qq.areas.items(): 

334 build = True 

335 loops = [area.loop] 

336 for hole_key, hole in iron.hyper_holes.items(): 

337 if area_name == hole.areas[1]: 

338 loops.append(qq.areas[hole.areas[0]].loop) 

339 elif area_name == hole.areas[0]: # or iron.hyper_areas[area_name].material == 'BH_air': 

340 build = False 

341 if build: 

342 area.surface = self.occ.addPlaneSurface(loops) 

343 self.md.domains.groups_surfaces[iron.hyper_areas[area_name].material].append(area.surface) 

344 

345 # Build coil domains 

346 for coil_nr, coil in gm.coil.coils.items(): 

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

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

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

350 # cable = self.set.Model_Data_GS.conductors[winding.conductor_name].cable 

351 # ht_area = cable.bare_cable_height_mean * cable.bare_cable_width 

352 for block_key, block in winding.blocks.items(): 

353 current_sign = '_neg' if block.current_sign < 0 else '_pos' 

354 for area_name, area in block.half_turns.areas.items(): 

355 ht_name = 'block' + str(block_key) + '_ht' + str(area_name) + current_sign 

356 # self.conductor_areas[ht_name] = ht_area 

357 area.surface = self.occ.addPlaneSurface([area.loop]) 

358 self.md.domains.groups_surfaces[ht_name] = [area.surface] 

359 

360 # Build wedges domains 

361 for wedge_nr, wedge in gm.wedges.areas.items(): 

362 wedge.surface = self.occ.addPlaneSurface([wedge.loop]) 

363 self.md.domains.groups_surfaces['wedge_' + str(wedge_nr)] = [wedge.surface] 

364 

365 def updateTags(self): 

366 self.md.geometries.coil.electrical_order = self.geom.Roxie_Data.coil.electrical_order 

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

368 self.md.geometries.coil.coils[coil_nr].bore_center = coil.bore_center 

369 

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

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

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

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

374 for block_key, block in winding.blocks.items(): 

375 hts = block.half_turns 

376 for ht_nr, ht in hts.areas.items(): 

377 lines_tags = gmsh.model.getAdjacencies(2, ht.surface)[1] 

378 first_tag = int(min(lines_tags)) 

379 hts.lines[ht_nr + ('i' if ht_nr == '1' else 'r')] = first_tag 

380 hts.lines[ht_nr + ('r' if ht_nr == '1' else 'i')] = first_tag + 1 

381 hts.lines[ht_nr + 'o'] = first_tag + 2 

382 if int(ht_nr) == len(hts.areas): 

383 hts.lines['end'] = first_tag + 3 

384 

385 lines_tags = gmsh.model.getAdjacencies(2, self.md.geometries.air_inf.areas['outer'].surface)[1] 

386 self.md.geometries.air_inf.lines['outer'] = int(lines_tags[0]) 

387 self.md.geometries.air_inf.lines['inner'] = int(lines_tags[1])