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
5import os
6import numpy as np
7import gmsh
8import copy
10import matplotlib.pyplot as plt
11from mpl_toolkits.mplot3d.art3d import Poly3DCollection
13from fiqus.parsers.ParserCOND import ParserCOND
15class direction(Enum):
16 """
17 A class to specify direction easily.
18 """
20 ccw = 0
21 cw = 1
23class coordinate(Enum):
24 """
25 A class to specify coordinate types easily.
26 """
28 rectangular = 0
29 cylindrical = 1
30 spherical = 2
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.
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 """
47 def __init__(self, r0=0.0, r1=0.0, r2=0.0, type=coordinate.rectangular) -> None:
49 self.type = type # Store 'type' as an instance attribute
51 if type is coordinate.rectangular:
52 self.x = r0
53 self.y = r1
54 self.z = r2
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!")
69 self.tag = gmsh.model.occ.addPoint(self.x, self.y, self.z)
71 def __repr__(self):
72 """
73 Returns the string representation of the point.
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)
80 def __abs__(self):
81 """
82 Returns the magnitude of the point vector.
84 :return: the magnitude of the point vector
85 :rtype: float
86 """
87 return math.hypot(self.x, self.y, self.z)
89 def __add__(self, other):
90 """
91 Returns the summation of two point vectors.
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)
103 def __mul__(self, scalar):
104 """
105 Returns the product of a point vector and a scalar.
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 )
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.
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 """
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
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
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)
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
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 )
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 )
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.
214 alpha = math.atan2(cutPlaneNormal[1], cutPlaneNormal[0]) - math.pi / 2
215 alpha2 = alpha + math.pi
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)
232 thetaArrays.append(np.array(listOfBreakPoints))
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]
240 r = innerRadius + (theta - initialTheta) / (2 * math.pi) * (gap) * multiplier
241 z = np.ones(theta.shape) * z
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 = []
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)
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)
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
286 # Defining the coordinate lists to which the points are to be added
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 = []
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
308 if len(indices) > 0:
309 windingUpIndex = indices[0] # Take the first index if there are multiple matches
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])
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]]
323 k = int(theta[i] / (2 * np.pi))
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)
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)
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)):
337 if k % 2 == 0:
338 winding_2.append(addPoints)
339 else:
340 winding_4.append(addPoints)
342 except IndexError:
343 print('All of the winding conductor points have been found')
345 x_coords = []
346 y_coords = []
347 z_coords = []
349 # Writing the conductor file
350 windingPointList = [winding_1, winding_2, winding_3, winding_4]
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'}}
355 # Use a dictionary to manage dict_cond_1, dict_cond_2, etc.
356 dict_conds = {}
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
371 target_dir = os.path.join(os.getcwd(), 'tests', '_outputs', 'parsers')
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)
385 # =============================================================================
386 # GENERATING POINTS ENDS ======================================================
387 # =============================================================================
389 # =============================================================================
390 # GENERATING SPLINES STARTS ===================================================
391 # =============================================================================
393 # Create the spline with the points:
394 spline = gmsh.model.occ.addSpline(pointTags)
396 # Split the spline into sub-curves:
397 sectionsPerCurve = int(spt / cpt)
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
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
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
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]
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]
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
460 if direction is direction.ccw:
461 self.splineTurns.append((startTurn, endTurn))
462 else:
463 self.splineTurns.append((-startTurn, -endTurn))
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])
469 # =============================================================================
470 # GENERATING SPLINES ENDS =====================================================
471 # =============================================================================
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]
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.
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]
501 def getPointTag(self, turn, endPoint=False):
502 """
503 Returns the point object at a specific turn.
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!")
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]
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.
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]