Coverage for fiqus/geom_generators/GeometryPancake3DUtils.py: 12%

246 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2025-03-14 02:31 +0100

1import math 

2from enum import Enum 

3from typing import List, Tuple, Dict 

4 

5import os 

6import numpy as np 

7import gmsh 

8import copy 

9 

10import matplotlib.pyplot as plt 

11from mpl_toolkits.mplot3d.art3d import Poly3DCollection 

12 

13from fiqus.parsers.ParserCOND import ParserCOND 

14 

15class direction(Enum): 

16 """ 

17 A class to specify direction easily. 

18 """ 

19 

20 ccw = 0 

21 cw = 1 

22 

23class coordinate(Enum): 

24 """ 

25 A class to specify coordinate types easily. 

26 """ 

27 

28 rectangular = 0 

29 cylindrical = 1 

30 spherical = 2 

31 

32class point: 

33 """ 

34 This is a class for creating points in GMSH. It supports rectangular and cylindrical 

35 coordinates. Moreover, vector operations are supported. 

36 

37 :param r0: x, r, or r (default: 0.0) 

38 :type r0: float, optional 

39 :param r1: y, theta, or theta (default: 0.0) 

40 :type r1: float, optional 

41 :param r2: z, z, or phi (default: 0.0) 

42 :type r2: float, optional 

43 :param type: coordinate type (default: coordinate.rectangular) 

44 :type type: coordinate, optional 

45 """ 

46 

47 def __init__(self, r0=0.0, r1=0.0, r2=0.0, type=coordinate.rectangular) -> None: 

48 

49 self.type = type # Store 'type' as an instance attribute 

50 

51 if type is coordinate.rectangular: 

52 self.x = r0 

53 self.y = r1 

54 self.z = r2 

55 

56 self.r = math.sqrt(self.x**2 + self.y**2) 

57 self.theta = math.atan2(self.y, self.x) 

58 elif type is coordinate.cylindrical: 

59 self.r = r0 

60 self.theta = r1 

61 self.x = self.r * math.cos(self.theta) 

62 self.y = self.r * math.sin(self.theta) 

63 self.z = r2 

64 elif type is coordinate.spherical: 

65 raise ValueError("Spherical coordinates are not supported yet!") 

66 else: 

67 raise ValueError("Improper coordinate type value!") 

68 

69 self.tag = gmsh.model.occ.addPoint(self.x, self.y, self.z) 

70 

71 def __repr__(self): 

72 """ 

73 Returns the string representation of the point. 

74 

75 :return: string representation of the point 

76 :rtype: str 

77 """ 

78 return "point(%r, %r, %r, %r)" % (self.x, self.y, self.z, self.type) 

79 

80 def __abs__(self): 

81 """ 

82 Returns the magnitude of the point vector. 

83 

84 :return: the magnitude of the point vector 

85 :rtype: float 

86 """ 

87 return math.hypot(self.x, self.y, self.z) 

88 

89 def __add__(self, other): 

90 """ 

91 Returns the summation of two point vectors. 

92 

93 :param other: point vector to be added 

94 :type other: point 

95 :return: the summation of two point vectors 

96 :rtype: point 

97 """ 

98 x = self.x + other.x 

99 y = self.y + other.y 

100 z = self.z + other.z 

101 return point(x, y, z, coordinate.rectangular) 

102 

103 def __mul__(self, scalar): 

104 """ 

105 Returns the product of a point vector and a scalar. 

106 

107 :param scalar: a scalar value 

108 :type scalar: float 

109 :return: point 

110 :rtype: point 

111 """ 

112 return point( 

113 self.x * scalar, 

114 self.y * scalar, 

115 self.z * scalar, 

116 coordinate.rectangular, 

117 ) 

118 

119 

120 

121class spiralCurve: 

122 """ 

123 A class to create a spiral curves parallel to XY plane in GMSH. The curve is defined 

124 by a spline and it is divided into sub-curves. Sub-curves are used because it makes 

125 the geometry creation process easier. 

126 

127 :param innerRadius: inner radius 

128 :type innerRadius: float 

129 :param gap: gap after each turn 

130 :type gap: float 

131 :param turns: number of turns 

132 :type turns: float 

133 :param z: z coordinate 

134 :type z: float 

135 :param initialTheta: initial theta angle in radians 

136 :type initialTheta: float 

137 :param direction: direction of the spiral (default: direction.ccw) 

138 :type direction: direction, optional 

139 :param cutPlaneNormal: normal vector of the plane that will cut the spiral curve 

140 (default: None) 

141 :type cutPlaneNormal: tuple[float, float, float], optional 

142 """ 

143 

144 # If the number of points used per turn is n, then the number of sections per turn 

145 # is n-1. They set the resolution of the spiral curve. It sets the limit of the 

146 # precision of the float number of turns that can be used to create the spiral 

147 # curve. The value below might be modified in Geometry.__init__ method. 

148 sectionsPerTurn = 16 

149 

150 # There will be curvesPerTurn curve entities per turn. It will be effectively the 

151 # number of volumes per turn in the end. The value below might be modified in 

152 # Geometry.__init__ method. 

153 curvesPerTurn = 2 

154 

155 def __init__( 

156 self, 

157 innerRadius, 

158 gap, 

159 turns, 

160 z, 

161 initialTheta, 

162 transitionNotchAngle, 

163 geo, 

164 direction=direction.ccw, # TODO code this to understand cw direction 

165 cutPlaneNormal=Tuple[float, float, float] 

166 ) -> None: 

167 spt = self.sectionsPerTurn # just to make the code shorter 

168 self.turnRes = 1 / spt # turn resolution 

169 cpt = self.curvesPerTurn # just to make the code shorter 

170 self.turns = turns 

171 self.geo = geo 

172 # ============================================================================= 

173 # GENERATING POINTS STARTS ==================================================== 

174 # ============================================================================= 

175 print('The theta is') 

176 print(initialTheta) 

177 

178 print('This is the status of the cutPlane') 

179 print(cutPlaneNormal) 

180 # Calculate the coordinates of the points that define the spiral curve: 

181 if direction is direction.ccw: 

182 # If the spiral is counter-clockwise, the initial theta angle decreases, 

183 # and r increases as the theta angle decreases. 

184 multiplier = 1 

185 elif direction is direction.cw: 

186 # If the spiral is clockwise, the initial theta angle increases, and r 

187 # increases as the theta angle increases. 

188 multiplier = -1 

189 

190 NofPointsPerTurn = int(spt + 1) 

191 thetaArrays = [] 

192 for turn in range(1, int(self.turns) + 1): 

193 thetaArrays.append( 

194 np.linspace( 

195 initialTheta + (turn - 1) * 2 * math.pi * multiplier, 

196 initialTheta + (turn) * 2 * math.pi * multiplier, 

197 NofPointsPerTurn, 

198 ) 

199 ) 

200 

201 thetaArrays.append( 

202 np.linspace( 

203 initialTheta + (turn) * 2 * math.pi * multiplier, 

204 initialTheta + (self.turns) * 2 * math.pi * multiplier, 

205 round(spt * (self.turns - turn) + 1), 

206 ) 

207 ) 

208 

209 if cutPlaneNormal is not None: 

210 # If the cutPlaneNormal is specified, the spiral curve will be cut by a 

211 # plane that is normal to the cutPlaneNormal vector and passes through the 

212 # origin. 

213 

214 alpha = math.atan2(cutPlaneNormal[1], cutPlaneNormal[0]) - math.pi / 2 

215 alpha2 = alpha + math.pi 

216 

217 listOfBreakPoints = [] 

218 for turn in range(1, int(self.turns) + 2): 

219 breakPoint1 = alpha + (turn - 1) * 2 * math.pi * multiplier 

220 breakPoint2 = alpha2 + (turn - 1) * 2 * math.pi * multiplier 

221 if ( 

222 breakPoint1 > initialTheta 

223 and breakPoint1 < initialTheta + 2 * math.pi * self.turns 

224 ): 

225 listOfBreakPoints.append(breakPoint1) 

226 if ( 

227 breakPoint2 > initialTheta 

228 and breakPoint2 < initialTheta + 2 * math.pi * self.turns 

229 ): 

230 listOfBreakPoints.append(breakPoint2) 

231 

232 thetaArrays.append(np.array(listOfBreakPoints)) 

233 

234 theta = np.concatenate(thetaArrays) 

235 theta = np.round(theta, 10) 

236 theta = np.unique(theta) 

237 theta = np.sort(theta) 

238 theta = theta[::multiplier] 

239 

240 r = innerRadius + (theta - initialTheta) / (2 * math.pi) * (gap) * multiplier 

241 z = np.ones(theta.shape) * z 

242 

243 # Create the points and store their tags: 

244 points = [] # point objects 

245 pointTags = [] # point tags 

246 breakPointObjectsDueToCutPlane = [] # only used if cutPlaneNormal is not None 

247 breakPointTagsDueToCutPlane = [] # only used if cutPlaneNormal is not None 

248 pointObjectsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None 

249 pointTagsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None 

250 breakPointObjectsDueToTransition = [] 

251 breakPointTagsDueToTransition = [] 

252 coordinateList = [] 

253 

254 for j in range(len(theta)): 

255 pointObject = point(r[j], theta[j], z[j], coordinate.cylindrical) 

256 [x_c, y_c, z_c] = [r[j], theta[j], z[j]] 

257 coordinateList.append([x_c, y_c, z_c]) 

258 points.append(pointObject) 

259 pointTags.append(pointObject.tag) 

260 if cutPlaneNormal is not None: 

261 if theta[j] in listOfBreakPoints: 

262 breakPointObjectsDueToCutPlane.append(pointObject) 

263 breakPointTagsDueToCutPlane.append(pointObject.tag) 

264 else: 

265 pointObjectsWithoutBreakPoints.append(pointObject) 

266 pointTagsWithoutBreakPoints.append(pointObject.tag) 

267 

268 # identify if the point is a break point due to the layer transition: 

269 angle1 = initialTheta + (2 * math.pi - transitionNotchAngle) * multiplier 

270 angle2 = ( 

271 initialTheta 

272 + ((self.turns % 1) * 2 * math.pi + transitionNotchAngle) * multiplier 

273 ) 

274 if math.isclose( 

275 math.fmod(theta[j], 2 * math.pi), angle1, abs_tol=1e-6 

276 ) or math.isclose(math.fmod(theta[j], 2 * math.pi), angle2, abs_tol=1e-6): 

277 breakPointObjectsDueToTransition.append(pointObject) 

278 breakPointTagsDueToTransition.append(pointObject.tag) 

279 

280 # Logic to break the points up into relevant geom coordinates 

281 # Brick points structure (for X-Y plane only for now): 

282 # [[[x1, y1, z1], [x2, y2, z2], [x3, y3, z3], [x4, y4, z4]], ...] 

283 # Theoretically, very easy to extend to 8 points knowing the height of the 

284 

285 

286 # Defining the coordinate lists to which the points are to be added 

287 

288 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number 

289 winding_1 = [] 

290 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number 

291 winding_2 = [] 

292 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number 

293 winding_3 = [] 

294 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number 

295 winding_4 = [] 

296 

297 for i in range(len(theta)-1): # range is reduced as no brick can be created starting at the last point 

298 # print(i) 

299 # Assuming theta is a numpy array and you're looking for the index of a value close to pi 

300 value_to_find = theta[i]+2*np.pi 

301 tolerance = 1e-10 # Define a small tolerance 

302 # Find indices where the condition is true 

303 indices = np.where(np.abs(theta - value_to_find) < tolerance)[0] 

304 z = self.geo.wi.h/2 

305 z_p = z 

306 z_n = -z 

307 

308 if len(indices) > 0: 

309 windingUpIndex = indices[0] # Take the first index if there are multiple matches 

310 

311 try: 

312 x_1 = r[i] * np.cos(theta[i]) 

313 y_1 = r[i] * np.sin(theta[i]) 

314 x_2 = r[i + 1] * np.cos(theta[i + 1]) 

315 y_2 = r[i + 1] * np.sin(theta[i + 1]) 

316 x_3 = r[windingUpIndex] * np.cos(theta[windingUpIndex]) 

317 y_3 = r[windingUpIndex] * np.sin(theta[windingUpIndex]) 

318 x_4 = r[windingUpIndex + 1] * np.cos(theta[windingUpIndex + 1]) 

319 y_4 = r[windingUpIndex + 1] * np.sin(theta[windingUpIndex + 1]) 

320 

321 addPoints = [[x_1, y_1, z_n], [x_2, y_2, z_n], [x_4, y_4, z_n], [x_3, y_3, z_n], [x_1, y_1, z_p], [x_2, y_2, z_p], [x_4, y_4, z_p], [x_3, y_3, z_p]] 

322 

323 k = int(theta[i] / (2 * np.pi)) 

324 

325 angle_in_current_turn = round(theta[i] % (2 * np.pi), 6) 

326 angle_second_brick_point = round(theta[i+1] % (2 * np.pi), 6) 

327 

328 if (((round(angle_in_current_turn, 6) <= round(np.pi, 6)) or (round(angle_in_current_turn, 6) == round(2*np.pi, 6))) and (round(angle_second_brick_point, 6) <= round(np.pi, 6))): 

329 k = int(round(theta[i] / (2 * np.pi))) 

330 if k % 2 == 0: 

331 winding_1.append(addPoints) 

332 else: 

333 winding_3.append(addPoints) 

334 

335 if (round(angle_in_current_turn, 6) >= round(np.pi, 6) and (round(angle_second_brick_point, 6) >= round(np.pi, 6) or round(angle_second_brick_point, 6) == 0.0)): 

336 

337 if k % 2 == 0: 

338 winding_2.append(addPoints) 

339 else: 

340 winding_4.append(addPoints) 

341 

342 except IndexError: 

343 print('All of the winding conductor points have been found') 

344 

345 x_coords = [] 

346 y_coords = [] 

347 z_coords = [] 

348 

349 # Writing the conductor file 

350 windingPointList = [winding_1, winding_2, winding_3, winding_4] 

351 

352 if self.geo.conductorWrite: 

353 dict_cond_sample = {0: {'SHAPE': 'BR8', 'XCENTRE': '0.0', 'YCENTRE': '0.0', 'ZCENTRE': '0.0', 'PHI1': '0.0', 'THETA1': '0.0', 'PSI1': '0.0', 'XCEN2': '0.0', 'YCEN2': '0.0', 'ZCEN2': '0.0', 'THETA2': '0.0', 'PHI2': '0.0', 'PSI2': '0.0', 'XP1': '-0.879570', 'YP1': '-0.002940', 'ZP1': '-1.131209', 'XP2': '-0.879570', 'YP2': '0.002940', 'ZP2': '-1.131209', 'XP3': '-0.881381', 'YP3': '0.002940', 'ZP3': '-1.114205', 'XP4': '-0.881381', 'YP4': '-0.002940', 'ZP4': '-1.114205', 'XP5': '-0.861227', 'YP5': '-0.002972', 'ZP5': '-1.129183', 'XP6': '-0.861208', 'YP6': '0.002908', 'ZP6': '-1.129182', 'XP7': '-0.863294', 'YP7': '0.002912', 'ZP7': '-1.112210', 'XP8': '-0.863313', 'YP8': '-0.002968', 'ZP8': '-1.112211', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}} 

354 

355 # Use a dictionary to manage dict_cond_1, dict_cond_2, etc. 

356 dict_conds = {} 

357 

358 k = 1 # Start from 1 for dict_cond_1, dict_cond_2, ... 

359 for winding in windingPointList: 

360 n = 0 

361 dict_conds[k] = {} 

362 for brick in winding: 

363 dict_conds[k][n] = copy.deepcopy(dict_cond_sample[0]) 

364 for pointIndex in range(8): 

365 dict_conds[k][n][f'XP{pointIndex + 1}'] = str(brick[pointIndex][0]) 

366 dict_conds[k][n][f'YP{pointIndex + 1}'] = str(brick[pointIndex][1]) 

367 dict_conds[k][n][f'ZP{pointIndex + 1}'] = str(brick[pointIndex][2]) 

368 n += 1 

369 k += 1 

370 

371 target_dir = os.path.join(os.getcwd(), 'tests', '_outputs', 'parsers') 

372 

373 # Ensure the target directory exists 

374 os.makedirs(target_dir, exist_ok=True) 

375 for i in range (1, 5): 

376 # Define the output file path 

377 out_file = os.path.join(target_dir, f'winding_{i}.cond') 

378 input_dict = dict_conds[i] 

379 list_of_shapes = ['BR8'] 

380 for shape in list_of_shapes: 

381 pc = ParserCOND() 

382 print(f'the input dictionary is {input_dict}') 

383 pc.write_cond(input_dict, out_file) 

384 

385 # ============================================================================= 

386 # GENERATING POINTS ENDS ====================================================== 

387 # ============================================================================= 

388 

389 # ============================================================================= 

390 # GENERATING SPLINES STARTS =================================================== 

391 # ============================================================================= 

392 

393 # Create the spline with the points: 

394 spline = gmsh.model.occ.addSpline(pointTags) 

395 

396 # Split the spline into sub-curves: 

397 sectionsPerCurve = int(spt / cpt) 

398 

399 # Create a list of point tags that will split the spline: 

400 # Basically, they are the points to divide the spirals into sectionsPerCurve 

401 # turns. However, some other points are also included to support the float 

402 # number of turns. It is best to visually look at the divisions with 

403 # gmsh.fltk.run() to understand why the split points are chosen the way they are 

404 # selected. 

405 if cutPlaneNormal is None: 

406 pointObjectsWithoutBreakPoints = points 

407 pointTagsWithoutBreakPoints = pointTags 

408 

409 splitPointTags = list( 

410 set(pointTagsWithoutBreakPoints[:-1:sectionsPerCurve]) 

411 | set(pointTagsWithoutBreakPoints[-spt - 1 :: -spt]) 

412 | set(breakPointTagsDueToCutPlane) 

413 | set(breakPointTagsDueToTransition) 

414 ) 

415 splitPointTags = sorted(splitPointTags) 

416 # Remove the first element of the list (starting point): 

417 _, *splitPointTags = splitPointTags 

418 

419 # Also create a list of corresponding point objects: 

420 splitPoints = list( 

421 set(pointObjectsWithoutBreakPoints[:-1:sectionsPerCurve]) 

422 | set(pointObjectsWithoutBreakPoints[-spt - 1 :: -spt]) 

423 | set(breakPointObjectsDueToCutPlane) 

424 | set(breakPointObjectsDueToTransition) 

425 ) 

426 splitPoints = sorted(splitPoints, key=lambda x: x.tag) 

427 # Remove the first element of the list (starting point): 

428 _, *splitPoints = splitPoints 

429 

430 # Split the spline: 

431 dims = [0] * len(splitPointTags) 

432 _, splines = gmsh.model.occ.fragment( 

433 [(1, spline)], 

434 list(zip(dims, splitPointTags)), 

435 removeObject=True, 

436 removeTool=True, 

437 ) 

438 splines = splines[0] 

439 self.splineTags = [j for _, j in splines] 

440 

441 # Note the turn number of each spline. This will be used in getSplineTag and 

442 # getSplineTags methods. 

443 self.splineTurns = [] 

444 for i in range(len(self.splineTags)): 

445 if i == 0: 

446 startPoint = points[0] 

447 endPoint = splitPoints[0] 

448 elif i == len(self.splineTags) - 1: 

449 startPoint = splitPoints[-1] 

450 endPoint = points[-1] 

451 else: 

452 startPoint = splitPoints[i - 1] 

453 endPoint = splitPoints[i] 

454 

455 startTurn = (startPoint.theta - initialTheta) / (2 * math.pi) 

456 startTurn = round(startTurn / self.turnRes) * self.turnRes 

457 endTurn = (endPoint.theta - initialTheta) / (2 * math.pi) 

458 endTurn = round(endTurn / self.turnRes) * self.turnRes 

459 

460 if direction is direction.ccw: 

461 self.splineTurns.append((startTurn, endTurn)) 

462 else: 

463 self.splineTurns.append((-startTurn, -endTurn)) 

464 

465 # Check if splineTurn tuples starts with the small turn number: 

466 for i in range(len(self.splineTurns)): 

467 self.splineTurns[i] = sorted(self.splineTurns[i]) 

468 

469 # ============================================================================= 

470 # GENERATING SPLINES ENDS ===================================================== 

471 # ============================================================================= 

472 

473 # Find start and end points of the spiral curve: 

474 gmsh.model.occ.synchronize() # synchronize the model to make getBoundary work 

475 self.startPointTag = gmsh.model.getBoundary([(1, self.getSplineTag(0))])[1][1] 

476 self.endPointTag = gmsh.model.getBoundary( 

477 [(1, self.getSplineTag(self.turns, endPoint=True))] 

478 )[1][1] 

479 

480 def getSplineTag(self, turn, endPoint=False): 

481 """ 

482 Returns the spline tag at a specific turn. It returns the spline tag of the 

483 section that is on the turn except its end point. 

484 

485 :param turn: turn number (it can be a float) 

486 :type turn: float 

487 :param endPoint: if True, return the spline tag of the section that is on the 

488 turn including its end point but not its start point (default: False) 

489 :type endPoint: bool, optional 

490 :return: spline tag 

491 """ 

492 if endPoint: 

493 for i, r in enumerate(self.splineTurns): 

494 if r[0] + (self.turnRes / 2) < turn <= r[1] + (self.turnRes / 2): 

495 return self.splineTags[i] 

496 else: 

497 for i, r in enumerate(self.splineTurns): 

498 if r[0] - (self.turnRes / 2) <= turn < r[1] - (self.turnRes / 2): 

499 return self.splineTags[i] 

500 

501 def getPointTag(self, turn, endPoint=False): 

502 """ 

503 Returns the point object at a specific turn. 

504 

505 :param turn: turn number (it can be a float) 

506 :type turn: float 

507 :return: point object 

508 :rtype: point 

509 """ 

510 if turn < 0 or turn > self.turns: 

511 raise ValueError("Turn number is out of range!") 

512 

513 if turn == 0: 

514 return self.startPointTag 

515 elif turn == self.turns: 

516 return self.endPointTag 

517 else: 

518 curveTag = self.getSplineTag(turn, endPoint=endPoint) 

519 if endPoint: 

520 points = gmsh.model.getBoundary([(1, curveTag)]) 

521 return points[1][1] 

522 else: 

523 points = gmsh.model.getBoundary([(1, curveTag)]) 

524 return points[0][1] 

525 

526 def getSplineTags(self, turnStart=None, turnEnd=None): 

527 """ 

528 Get the spline tags from a specific turn to another specific turn. If turnStart 

529 and turnEnd are not specified, it returns all the spline tags. 

530 

531 :param turnStart: start turn number (it can be a float) (default: None) 

532 :type turnStart: float, optional 

533 :param turnEnd: end turn number (it can be a float) (default: None) 

534 :type turnEnd: float, optional 

535 :return: spline tags 

536 :rtype: list[int] 

537 """ 

538 if turnStart is None and turnEnd is None: 

539 return self.splineTags 

540 elif turnStart is None or turnEnd is None: 

541 raise ValueError( 

542 "turnStart and turnEnd must be both specified or both not specified." 

543 " You specified only one of them." 

544 ) 

545 else: 

546 start = self.splineTags.index(self.getSplineTag(turnStart, False)) 

547 end = self.splineTags.index(self.getSplineTag(turnEnd, True)) + 1 

548 return self.splineTags[start:end] 

549