Coverage for fiqus/geom_generators/GeometryPancake3D.py: 82%
1473 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-03-14 02:31 +0100
« prev ^ index » next coverage.py v7.4.4, created at 2025-03-14 02:31 +0100
1import math
2import logging
3from enum import Enum
4from inspect import currentframe, getframeinfo
5from typing import List, Tuple, Dict
6import operator
8import os
9import json
10import timeit
11import numpy as np
12import gmsh
14from fiqus.utils.Utils import GmshUtils
15from fiqus.utils.Utils import FilesAndFolders
16from fiqus.mains.MainPancake3D import Base
17from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry
19logger = logging.getLogger(__name__)
21# coordinateList = []
23def findSurfacesWithNormalsOnXYPlane(dimTags):
24 result = []
25 for dimTag in dimTags:
26 surfaceNormal = gmsh.model.getNormal(dimTag[1], [0.5, 0.5])
27 if abs(surfaceNormal[2]) < 1e-6:
28 result.append(dimTag)
30 return result
33def findOuterOnes(dimTags, findInnerOnes=False):
34 """
35 Finds the outermost surface/curve/point in a list of dimTags. The outermost means
36 the furthest from the origin.
37 """
38 dim = dimTags[0][0]
39 if dim == 2:
40 distances = []
41 for dimTag in dimTags:
42 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])
43 for curve in curves:
44 curve = list(curve)
45 pointTags = gmsh.model.getBoundary(
46 [(1, curveTag) for curveTag in curve],
47 oriented=False,
48 combined=False,
49 )
50 # Get the positions of the points:
51 points = []
52 for dimTag in pointTags:
53 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]
54 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]
55 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))
56 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))
58 distances.append(
59 max([point[0] ** 2 + point[1] ** 2 for point in points])
60 )
61 elif dim == 1:
62 distances = []
63 for dimTag in dimTags:
64 pointTags = gmsh.model.getBoundary(
65 [dimTag],
66 oriented=False,
67 combined=False,
68 )
69 # Get the positions of the points:
70 points = []
71 for dimTag in pointTags:
72 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]
73 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]
74 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))
75 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))
77 distances.append(max([point[0] ** 2 + point[1] ** 2 for point in points]))
79 if findInnerOnes:
80 goalDistance = min(distances)
81 else:
82 goalDistance = max(distances)
84 result = []
85 for distance, dimTag in zip(distances, dimTags):
86 # Return all the dimTags with the hoal distance:
87 if math.isclose(distance, goalDistance, abs_tol=1e-6):
88 result.append(dimTag)
90 return result
93class dimTags:
94 """
95 This is a class for (dim, tag) tuple lists. DimTags are heavily used in GMSH, and
96 dimTags class makes it easier to store and manipulate them.
98 Every dimTags instance with the save option True will be stored in the
99 dimTagsStorage class. dimTags instances with the save option False will not be
100 stored. dimTagsStorage class stores all the dimTags with the corresponding names so
101 that later the dimTagsStorage class can be used to create the volume information
102 file. The volume information file is required in the meshing stage to create the
103 appropriate regions.
105 If parentName is specified, during the volume information file generation, another
106 key that is equal to parentName is created, and all the dimTags are also added
107 there. For example, air consists of many parts such as gap, outer tube, inner
108 cylinder parts, and their dimTags object all have the parentName of self.geo.air.name
109 (user input of the air name) so that in the volume information file, they are all
110 combined under the air name as well.
112 :param name: name of the dimTags object (default: None)
113 :type name: str, optional
114 :param parentName: name of the parent dimTags object (default: None)
115 :type parentName: str, optional
116 :param save: True if the instance to be stored in dimTagsStorage class, False
117 otherwise (default: False)
118 :type save: bool, optional
119 """
121 point = 0
122 curve = 1
123 surface = 2
124 volume = 3
126 def storageUpdateRequired(func):
127 """
128 A decorator for dimTags class. It will update the dimTagsStorage class if the
129 save option is True and the decorated function is called. Use this decorator
130 for every function that changes the dimTags instance so that the storage gets
131 updated.
133 :param func: function to be decorated
134 :type func: function
135 :return: decorated function
136 :rtype: function
137 """
139 def wrapper(self, *args, **kwargs):
140 func(self, *args, **kwargs)
142 if self.save:
143 # Update the dimTagsStorage:
144 dimTagsStorage.updateDimTags(self)
146 return wrapper
148 def __init__(
149 self,
150 name: str = None,
151 parentName: str = None,
152 save: bool = False,
153 ):
154 self.name = name
156 dimTagsObjects = dimTagsStorage.getDimTagsObject(name)
157 if dimTagsObjects != []:
158 dimTagsObject = dimTagsObjects[0]
159 self.physicalTag = dimTagsObject.physicalTag
160 # To store points, curves, surfaces, and volumes separately:
161 self.dimTags = dimTagsObject.dimTags
162 self.dimTagsForPG = dimTagsObject.dimTagsForPG
163 self.allDimTags = dimTagsObject.allDimTags
164 self.parentName = dimTagsObject.parentName
165 self.save = dimTagsObject.save
166 else:
167 self.physicalTag = None
168 # To store points, curves, surfaces, and volumes separately:
169 self.dimTags = [[] for _ in range(4)]
170 self.dimTagsForPG = [[] for _ in range(4)] # dim tags for physical groups
171 self.allDimTags = []
172 self.parentName = parentName
173 self.save = save
174 if self.save:
175 dimTagsStorage.updateDimTags(self)
177 @storageUpdateRequired
178 def add(
179 self,
180 dimTagsList: List[Tuple[int, int]],
181 dimTagsListForPG: List[Tuple[int, int]] = None,
182 ):
183 """
184 Adds a list of (dim, tag) tuples to the dimTags object.
186 dimTagsListForPG is also accepted as an argument because sometimes, the stored
187 dimTags and the current dimTags in the geometry generation can be different. For
188 example, if volume 61 is deleted, the other volume tags (62, 63, ...) won't
189 shift back. However, after saving the geometry as a BREP file and rereading it,
190 the volume tags will be shifted back. In this case, the stored dimTags should be
191 shifted as well. But to create the correct physical region in the
192 geometry-creating process (which won't be saved in the BREP file, just used for
193 debugging purposes), another optional dimTagsListForPG argument is accepted.
195 :param dimTagsList: list of (dim, tag) tuples
196 :type dimTagsList: list[tuple[int, int]]
197 :param dimTagsListForPG: list of (dim, tag) tuples for physical groups
198 (default: None). If dimTagsListForPG is None, dimTagsList will be used for
199 physical groups as well.
200 :type dimTagsListForPG: list[tuple[int, int]], optional
202 """
203 if not isinstance(dimTagsList, list):
204 dimTagsList = [dimTagsList]
206 if not all(isinstance(element, tuple) for element in dimTagsList):
207 raise TypeError("Dim tags must be a list of tuples!")
209 # Sometimes, negative entities can be added for topology.
210 for i, v in enumerate(dimTagsList):
211 if v[1] < 0:
212 dimTagsList[i] = (v[0], -v[1])
214 # Add dim tags if they are not already added:
215 for v in dimTagsList:
216 if v not in self.allDimTags:
217 self.dimTags[v[0]].append(v)
218 if dimTagsListForPG is None:
219 self.dimTagsForPG[v[0]].append(v)
220 self.allDimTags.append(v)
222 if dimTagsListForPG is not None:
223 if not isinstance(dimTagsListForPG, list):
224 dimTagsListForPG = [dimTagsListForPG]
226 if not all(isinstance(element, tuple) for element in dimTagsListForPG):
227 raise TypeError("Dim tags must be a list of tuples!")
229 for i, v in enumerate(dimTagsListForPG):
230 if v[1] < 0:
231 dimTagsListForPG[i] = (v[0], -v[1])
232 for v in dimTagsListForPG:
233 if v not in self.dimTagsForPG:
234 self.dimTagsForPG[v[0]].append(v)
236 def addWithTags(self, dim: int, tags: List[int], tagsForPG: List[int] = None):
237 """
238 Adds a list of tags with a specific dimension to the dimTags object. The
239 explanation of the tagsForPG argument is given in the add() method.
241 :param dim: dimension of the tags
242 :type dim: int
243 :param tags: list of tags
244 :type tags: list[tuple[int, int]]
245 :param tagsForPG: list of tags for physical groups (default: None)
246 :type tagsForPG: list[tuple[int, int]], optional
248 """
249 if not isinstance(tags, list):
250 tags = [tags]
252 if not isinstance(dim, int):
253 raise TypeError("Dimension must be an integer!")
255 if not all(isinstance(element, int) for element in tags):
256 raise TypeError("Tags must be a list of integers!")
258 dims = [dim] * len(tags)
259 dimTagsList = list(zip(dims, tags))
261 if tagsForPG is not None:
262 if not isinstance(tagsForPG, list):
263 tagsForPG = [tagsForPG]
265 if not all(isinstance(element, int) for element in tagsForPG):
266 raise TypeError("Tags must be a list of integers!")
268 dimTagsListForPG = list(zip(dims, tagsForPG))
269 self.add(dimTagsList, dimTagsListForPG)
270 else:
271 self.add(dimTagsList)
273 def getDimTags(self, dim: int = None) -> List[Tuple[int, int]]:
274 """
275 Returns the stored list of (dim, tag) tuples with a specific dimension. If dim
276 is not specified, returns all the (dim, tag) tuples.
278 :param dim: dimension of the tags to be returned (default: None)
279 :type dim: int, optional
280 :return: list of (dim, tag) tuples
281 :rtype: list[tuple[int, int]]
282 """
283 if dim is None:
284 return self.dimTags[0] + self.dimTags[1] + self.dimTags[2] + self.dimTags[3]
285 else:
286 return self.dimTags[dim]
288 def getDimTagsForPG(self, dim: int = None) -> List[Tuple[int, int]]:
289 """
290 Returns the stored list of (dim, tag) tuples for physical groups with a specific
291 dimension. If dim is not specified, returns all the (dim, tag) tuples for
292 physical groups.
294 :param dim: dimension of the tags to be returned (default: None)
295 :type dim: int, optional
296 :return: list of (dim, tag) tuples for physical groups
297 :rtype: list[tuple[int, int]]
298 """
299 if dim is None:
300 return (
301 self.dimTagsForPG[0]
302 + self.dimTagsForPG[1]
303 + self.dimTagsForPG[2]
304 + self.dimTagsForPG[3]
305 )
306 else:
307 return self.dimTagsForPG[dim]
309 def getTags(self, dim: int, forPhysicalGroup=False) -> List[int]:
310 """
311 Returns the stored list of tags with a specific dimension.
313 :param dim: dimension of the tags to be returned
314 :type dim: int
315 :return: list of tags
316 :rtype: list[int]
317 """
318 if forPhysicalGroup:
319 return [v[1] for v in self.dimTagsForPG[dim]]
320 else:
321 return [v[1] for v in self.dimTags[dim]]
323 def getExtrusionTop(self, dim=3):
324 """
325 Returns the top surfaces, lines, or points of an extrusion operation if the
326 dimTags object contains the tags of an extrusion operation.
327 gmsh.model.occ.extrusion() function returns all the entities that are created
328 by the extrusion as a dim tags list. The first element is always the top
329 surface, the second is the volume. However, when more than one surface is
330 extruded, extracting the top surfaces is not trivial. This function returns
331 the top surfaces of an extrusion operation. It does that by finding the entities
332 right before the volume entities. If dim is 2, the function will return the top
333 curves of an extrusion operation. If dim is 1, the function will return the top
334 points of an extrusion.
336 :param dim: dimension of the entity that is being created (default: 3)
337 :type dim: int, optional
338 :return: list of (dim, tag) tuples of the top entities of an extrusion operation
339 :rtype: list[tuple[int, int]]
340 """
342 topSurfaces = []
343 for index, dimTag in enumerate(self.allDimTags):
344 if dimTag[0] == dim:
345 topSurfaces.append(self.allDimTags[index - 1])
347 return topSurfaces
349 def getExtrusionSide(self, dim=3):
350 """
351 Returns the side surfaces, lines, or points of an extrusion operation if the
352 dimTags object contains the tags of an extrusion operation.
353 gmsh.model.occ.extrusion() function returns all the entities that are created
354 by the extrusion as a dim tags list. The first element is always the top
355 surface, the second is the volume. The other elements are the side surfaces.
356 However, when more than one surface is extruded, extracting the side surfaces
357 is not trivial. This function returns the side surfaces of an extrusion
358 operation. It does that by finding returning all the entities except the top
359 surface and the volume. If dim is 2, the function will return the side curves of
360 an extrusion operation.
362 :param dim: dimension of the entity that is being created (default: 3)
363 :type dim: int, optional
364 :return: list of (dim, tag) tuples of the side entities of an extrusion operation
365 :rtype: list[tuple[int, int]]
366 """
367 sideSurfaces = []
368 sideSurfaceStartIndex = None
369 for index, dimTag in enumerate(self.allDimTags):
370 if dimTag[0] == dim:
371 if sideSurfaceStartIndex is not None:
372 sideSurfaces.append(
373 self.allDimTags[sideSurfaceStartIndex : index - 1]
374 )
375 sideSurfaceStartIndex = index + 1
376 else:
377 sideSurfaceStartIndex = index + 1
379 sideSurfaces.append(self.allDimTags[sideSurfaceStartIndex:])
381 return sideSurfaces
383 def __add__(self, other):
384 """
385 Adds two dimTags objects and returns a new dimTags object with the same save and
386 name attirbues of the first dimTags object.
388 It might cause bugs because of the recursive behavior of the
389 @storageUpdateRequired decorator. Use with caution. Currently only used by
390 dimTagsStorage.updateDimTags method.
392 :param other: dimTags object to be added
393 :type other: dimTags
394 :return: dimTags object with the sum of the two dimTags objects
395 :rtype: dimTags
396 """
397 result = dimTags()
398 result.name = self.name
399 result.parentName = self.parentName
400 result.physicalTag = self.physicalTag
401 result.dimTags = self.dimTags
402 result.dimTagsForPG = self.dimTagsForPG
403 result.allDimTags = self.allDimTags
405 result.add(other.allDimTags)
406 result.save = self.save
407 return result
409 def __repr__(self):
410 """
411 Returns the string representation of the dimTags object. If the dimTags object
412 is saved, it will return "SAVED: name". If the dimTags object is not saved, it
413 will return "NOT SAVED: name".
415 dimTags objects are used as dictionary keys throughout the code. This
416 representation makes debugging easier.
418 :return: string representation of the dimTags object
419 :rtype: str
420 """
421 if self.save:
422 return "SAVED: " + self.name
423 else:
424 return "NOT SAVED: " + self.name
427class dimTagsStorage:
428 """
429 This is a global class to store the dimTags of important entities in the model.
430 Every dimTags instance with self.save = True will be stored in this class. Later,
431 the storage will be used to generate the volume information (*.vi) file. *.vi file
432 will be used for generating the physical regions in the meshing part.
434 Users should not use this class directly. Instead, they should use the dimTags
435 class. If they assign save = True to the dimTags instance, it will be stored in this
436 class.
438 This class is a singleton class. It means that there will be only one instance of
439 this class in the whole module. This is done to be able to use the same storage
440 throughout this module. See the singleton design pattern for more information.
441 """
443 __instance = None
444 __dimTagsDict = {} # Dictionary with the names of the dimTags objects as keys and
445 # dimTags objects as values
447 def __new__(cls):
448 if cls.__instance is None:
449 cls.__instance = super().__new__(cls)
451 return cls.__instance
453 @classmethod
454 def updateDimTags(cls, dimTagsObject: dimTags):
455 """
456 Either adds or updates the dimTags object in the storage.
458 :param dimTags: dimTags object to be added or updated.
459 :type dimTags: dimTags
461 """
462 if dimTagsObject.name in cls.__dimTagsDict:
463 newDimTags = dimTagsObject + cls.__dimTagsDict[dimTagsObject.name]
464 cls.__dimTagsDict[dimTagsObject.name] = newDimTags
465 else:
466 cls.__dimTagsDict[dimTagsObject.name] = dimTagsObject
468 # @classmethod
469 # def updateDimTagsFromDict(cls, name: str, dimTagsList: List[Tuple[int, int]]):
470 # """
471 # Updates or adds dimTags from a list of (dim, tag) tuples.
472 #
473 # :param name: Name of the dimTags entry to update or add.
474 # :type name: str
475 # :param dimTagsList: List of (dim, tag) tuples to be associated with the name.
476 # :type dimTagsList: List[Tuple[int, int]]
477 # """
478 # # Check if the entry exists; if so, update it, otherwise add a new entry
479 # if name in cls.__dimTagsDict:
480 # existingDimTags = cls.__dimTagsDict[
481 # name].dimTags # Assuming dimTags object has a dimTags attribute storing the list of tuples
482 # updatedDimTags = existingDimTags + dimTagsList
483 # cls.__dimTagsDict[name].dimTags = updatedDimTags # Update the list of tuples
484 # else:
485 # # Create a new dimTags object (this step depends on how your dimTags class is structured)
486 # newDimTagsObject = dimTags(name=name, dimTags=dimTagsList) # Assuming such a constructor exists
487 # cls.__dimTagsDict[name] = newDimTagsObject
489 @classmethod
490 def getDimTagsObject(cls, names: List[str]):
491 """
492 Returns the dimTags object with the given names.
494 :param names: names of the dimTags objects.
495 :type names: list[str]
496 :return: dimTags objects with the given name.
497 :rtype: list[dimTags]
498 """
499 if not isinstance(names, list):
500 names = [names]
502 dimTagsObjects = []
503 for name in names:
504 if name in cls.__dimTagsDict:
505 dimTagsObjects.append(cls.__dimTagsDict[name])
507 return dimTagsObjects
509 @classmethod
510 def getDimTags(cls, names: List[str], dim: int = None) -> List[Tuple[int, int]]:
511 """
512 Returns the stored list of (dim, tag) tuples with dimension a specific dimenions
513 and names. If dim is not specified, all the stored (dim, tag) tuples under the
514 given names will be returned.
516 :param names: names of the dimTags object that will be returned
517 :type names: list[str]
518 :param dim: dimension of the (dim, tag) tuples to be returned (default: None).
519 If dim is None, all the stored (dim, tag) tuples under the name names will
520 be returned.
521 :type dim: int, optional
522 :return: list of (dim, tag) tuples
523 """
524 if not isinstance(names, list):
525 names = [names]
527 dimTagsResult = []
528 for name in names:
529 dimTagsResult.extend(cls.__dimTagsDict[name].getDimTags(dim))
531 return dimTagsResult
533 @classmethod
534 def getTags(cls, names: List[str], dim: int) -> List[int]:
535 """
536 Returns the stored list of tags with dimension a specific dimension and names.
538 :param names: names of the dimTags objects
539 :type names: list[str]
540 :param dim: dimension of the tags to be returned
541 :type dim: int
542 :return: list of tags
543 :rtype: list[int]
544 """
545 dimTags = cls.getDimTags(names, dim)
546 tags = [dimTag[1] for dimTag in dimTags]
548 return tags
550 @classmethod
551 def getDimTagsDict(
552 cls, forPhysicalGroups=False
553 ) -> Dict[str, List[Tuple[int, int]]]:
554 """
555 Returns a dictionary with the names of the dimTags objects as keys and the
556 stored list of (dim, tag) tuples as values. This method is used to generate the
557 .vi file. If forPhysicalGroups is True, the dimTags for physical groups will be
558 returned instead of the dimTags.
560 :param forPhysicalGroups: True if the dimTags for physical groups should be
561 returned, False otherwise (default: False)
562 :type forPhysicalGroups: bool, optional
563 :return: dictionary with the names of the dimTags objects as keys and the
564 stored list of (dim, tag) tuples as values
565 :rtype: dict[str, list[tuple[int, int]]]
566 """
567 dictionary = {}
568 for name, dimTagsObject in cls.__dimTagsDict.items():
569 if dimTagsObject.parentName is not None:
570 if dimTagsObject.parentName in dictionary:
571 if forPhysicalGroups:
572 dictionary[dimTagsObject.parentName].extend(
573 dimTagsObject.getDimTagsForPG()
574 )
575 else:
576 dictionary[dimTagsObject.parentName].extend(
577 dimTagsObject.getDimTags()
578 )
579 else:
580 if forPhysicalGroups:
581 dictionary[dimTagsObject.parentName] = (
582 dimTagsObject.getDimTagsForPG()
583 )
584 else:
585 dictionary[dimTagsObject.parentName] = (
586 dimTagsObject.getDimTags()
587 )
588 if forPhysicalGroups:
589 dictionary[name] = dimTagsObject.getDimTagsForPG()
590 else:
591 dictionary[name] = dimTagsObject.getDimTags()
593 return dictionary
595 @classmethod
596 def getAllStoredDimTags(cls) -> List[Tuple[int, int]]:
597 """
598 Returns a list of all the stored (dim, tag) tuples, regardless of the name of
599 the dimTags object (i.e. all the dimTags objects are merged into one list).
601 :return: list of all the stored (dim, tag) tuples.
602 :rtype: list[tuple[int, int]]
603 """
604 AllStoredDimTags = []
605 for name, dimTagsObject in cls.__dimTagsDict.items():
606 AllStoredDimTags.extend(dimTagsObject.getDimTags())
608 return AllStoredDimTags
610 @classmethod
611 def clear(cls):
612 """
613 Clears the dimTagsStorage class.
616 """
617 cls.__instance = None
618 cls.__dimTagsDict = (
619 {}
620 ) # Dictionary with the names of the dimTags objects as keys and
621 # dimTags objects as values
624class coordinate(Enum):
625 """
626 A class to specify coordinate types easily.
627 """
629 rectangular = 0
630 cylindrical = 1
631 spherical = 2
634class direction(Enum):
635 """
636 A class to specify direction easily.
637 """
639 ccw = 0
640 cw = 1
643class point:
644 """
645 This is a class for creating points in GMSH. It supports rectangular and cylindrical
646 coordinates. Moreover, vector operations are supported.
648 :param r0: x, r, or r (default: 0.0)
649 :type r0: float, optional
650 :param r1: y, theta, or theta (default: 0.0)
651 :type r1: float, optional
652 :param r2: z, z, or phi (default: 0.0)
653 :type r2: float, optional
654 :param type: coordinate type (default: coordinate.rectangular)
655 :type type: coordinate, optional
656 """
658 def __init__(self, r0=0.0, r1=0.0, r2=0.0, type=coordinate.rectangular) -> None:
660 self.type = type # Store 'type' as an instance attribute
662 if type is coordinate.rectangular:
663 self.x = r0
664 self.y = r1
665 self.z = r2
667 self.r = math.sqrt(self.x**2 + self.y**2)
668 self.theta = math.atan2(self.y, self.x)
669 elif type is coordinate.cylindrical:
670 self.r = r0
671 self.theta = r1
672 self.x = self.r * math.cos(self.theta)
673 self.y = self.r * math.sin(self.theta)
674 self.z = r2
675 elif type is coordinate.spherical:
676 raise ValueError("Spherical coordinates are not supported yet!")
677 else:
678 raise ValueError("Improper coordinate type value!")
680 self.tag = gmsh.model.occ.addPoint(self.x, self.y, self.z)
682 def __repr__(self):
683 """
684 Returns the string representation of the point.
686 :return: string representation of the point
687 :rtype: str
688 """
689 return "point(%r, %r, %r, %r)" % (self.x, self.y, self.z, self.type)
691 def __abs__(self):
692 """
693 Returns the magnitude of the point vector.
695 :return: the magnitude of the point vector
696 :rtype: float
697 """
698 return math.hypot(self.x, self.y, self.z)
700 def __add__(self, other):
701 """
702 Returns the summation of two point vectors.
704 :param other: point vector to be added
705 :type other: point
706 :return: the summation of two point vectors
707 :rtype: point
708 """
709 x = self.x + other.x
710 y = self.y + other.y
711 z = self.z + other.z
712 return point(x, y, z, coordinate.rectangular)
714 def __mul__(self, scalar):
715 """
716 Returns the product of a point vector and a scalar.
718 :param scalar: a scalar value
719 :type scalar: float
720 :return: point
721 :rtype: point
722 """
723 return point(
724 self.x * scalar,
725 self.y * scalar,
726 self.z * scalar,
727 coordinate.rectangular,
728 )
731class spiralCurve():
732 """
733 A class to create a spiral curves parallel to XY plane in GMSH. The curve is defined
734 by a spline and it is divided into sub-curves. Sub-curves are used because it makes
735 the geometry creation process easier.
737 :param innerRadius: inner radius
738 :type innerRadius: float
739 :param gap: gap after each turn
740 :type gap: float
741 :param turns: number of turns
742 :type turns: float
743 :param z: z coordinate
744 :type z: float
745 :param initialTheta: initial theta angle in radians
746 :type initialTheta: float
747 :param direction: direction of the spiral (default: direction.ccw)
748 :type direction: direction, optional
749 :param cutPlaneNormal: normal vector of the plane that will cut the spiral curve
750 (default: None)
751 :type cutPlaneNormal: tuple[float, float, float], optional
752 """
754 # If the number of points used per turn is n, then the number of sections per turn
755 # is n-1. They set the resolution of the spiral curve. It sets the limit of the
756 # precision of the float number of turns that can be used to create the spiral
757 # curve. The value below might be modified in Geometry.__init__ method.
758 sectionsPerTurn = 16
760 # There will be curvesPerTurn curve entities per turn. It will be effectively the
761 # number of volumes per turn in the end. The value below might be modified in
762 # Geometry.__init__ method.
763 curvesPerTurn = 2
765 def __init__(
766 self,
767 innerRadius,
768 gap,
769 turns,
770 z,
771 initialTheta,
772 transitionNotchAngle,
773 direction=direction.ccw,
774 cutPlaneNormal=Tuple[float, float, float],
775 ) -> None:
776 spt = self.sectionsPerTurn # just to make the code shorter
777 self.turnRes = 1 / spt # turn resolution
778 cpt = self.curvesPerTurn # just to make the code shorter
779 self.turns = turns
781 # =============================================================================
782 # GENERATING POINTS STARTS ====================================================
783 # =============================================================================
785 # Calculate the coordinates of the points that define the spiral curve:
786 if direction is direction.ccw:
787 # If the spiral is counter-clockwise, the initial theta angle decreases,
788 # and r increases as the theta angle decreases.
789 multiplier = 1
790 elif direction is direction.cw:
791 # If the spiral is clockwise, the initial theta angle increases, and r
792 # increases as the theta angle increases.
793 multiplier = -1
795 NofPointsPerTurn = int(spt + 1)
796 thetaArrays = []
797 for turn in range(1, int(self.turns) + 1):
798 thetaArrays.append(
799 np.linspace(
800 initialTheta + (turn - 1) * 2 * math.pi * multiplier,
801 initialTheta + (turn) * 2 * math.pi * multiplier,
802 NofPointsPerTurn,
803 )
804 )
806 thetaArrays.append(
807 np.linspace(
808 initialTheta + (turn) * 2 * math.pi * multiplier,
809 initialTheta + (self.turns) * 2 * math.pi * multiplier,
810 round(spt * (self.turns - turn) + 1),
811 )
812 )
814 if cutPlaneNormal is not None:
815 # If the cutPlaneNormal is specified, the spiral curve will be cut by a
816 # plane that is normal to the cutPlaneNormal vector and passes through the
817 # origin.
819 alpha = math.atan2(cutPlaneNormal[1], cutPlaneNormal[0]) - math.pi / 2
820 alpha2 = alpha + math.pi
822 listOfBreakPoints = []
823 for turn in range(1, int(self.turns) + 2):
824 breakPoint1 = alpha + (turn - 1) * 2 * math.pi * multiplier
825 breakPoint2 = alpha2 + (turn - 1) * 2 * math.pi * multiplier
826 if (
827 breakPoint1 > initialTheta
828 and breakPoint1 < initialTheta + 2 * math.pi * self.turns
829 ):
830 listOfBreakPoints.append(breakPoint1)
831 if (
832 breakPoint2 > initialTheta
833 and breakPoint2 < initialTheta + 2 * math.pi * self.turns
834 ):
835 listOfBreakPoints.append(breakPoint2)
837 thetaArrays.append(np.array(listOfBreakPoints))
839 theta = np.concatenate(thetaArrays)
840 theta = np.round(theta, 10)
841 theta = np.unique(theta)
842 theta = np.sort(theta)
843 theta = theta[::multiplier]
845 r = innerRadius + (theta - initialTheta) / (2 * math.pi) * (gap) * multiplier
846 z = np.ones(theta.shape) * z
848 # Create the points and store their tags:
849 points = [] # point objects
850 pointTags = [] # point tags
851 breakPointObjectsDueToCutPlane = [] # only used if cutPlaneNormal is not None
852 breakPointTagsDueToCutPlane = [] # only used if cutPlaneNormal is not None
853 pointObjectsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None
854 pointTagsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None
855 breakPointObjectsDueToTransition = []
856 breakPointTagsDueToTransition = []
857 coordinateList = []
859 for j in range(len(theta)):
860 pointObject = point(r[j], theta[j], z[j], coordinate.cylindrical)
861 [x_c, y_c, z_c] = [r[j], theta[j], z[j]]
862 #print([x_c, y_c, z_c])
863 coordinateList.append([x_c, y_c, z_c])
864 points.append(pointObject)
865 pointTags.append(pointObject.tag)
866 if cutPlaneNormal is not None:
867 if theta[j] in listOfBreakPoints:
868 breakPointObjectsDueToCutPlane.append(pointObject)
869 breakPointTagsDueToCutPlane.append(pointObject.tag)
870 else:
871 pointObjectsWithoutBreakPoints.append(pointObject)
872 pointTagsWithoutBreakPoints.append(pointObject.tag)
874 # identify if the point is a break point due to the layer transition:
875 angle1 = initialTheta + (2 * math.pi - transitionNotchAngle) * multiplier
876 angle2 = (
877 initialTheta
878 + ((self.turns % 1) * 2 * math.pi + transitionNotchAngle) * multiplier
879 )
880 if math.isclose(
881 math.fmod(theta[j], 2 * math.pi), angle1, abs_tol=1e-6
882 ) or math.isclose(math.fmod(theta[j], 2 * math.pi), angle2, abs_tol=1e-6):
883 breakPointObjectsDueToTransition.append(pointObject)
884 breakPointTagsDueToTransition.append(pointObject.tag)
886 # Plotter
887 # x_coords = [coord[0] for coord in coordinateList]
888 # y_coords = [coord[1] for coord in coordinateList]
889 # z_coords = [coord[2] for coord in coordinateList]
891 # print(f'number of divisions {self.Pancake3DMeshWinding.ane}')
892 # print(self.wi.ane)
893 # Creating the 3D plot
894 # fig = plt.figure()
895 # ax = fig.add_subplot(111, projection='3d')
897 # Plotting the coordinates
898 # ax.scatter(x_coords, y_coords, z_coords, c='r', marker='o')
899 #
900 # # Setting labels
901 # ax.set_xlabel('X Label')
902 # ax.set_ylabel('Y Label')
903 # ax.set_zlabel('Z Label')
904 #
905 # plt.show()
907 # Logic to break the points up into relevant geom coordinates
908 # Brick points structure (for X-Y plane only for now):
909 # [[[x1, y1, z1], [x2, y2, z2], [x3, y3, z3], [x4, y4, z4]], ...]
910 # Theoretically, very easy to extend to 8 points knowing the height of the
912 # Defining the coordinate lists to which the points are to be added
914 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number
915 winding_1 = []
916 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number
917 winding_2 = []
918 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number
919 winding_3 = []
920 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number
921 winding_4 = []
922 #print(theta[10])
923 # heightPancake = self.geo.winding.height
924 # print(heightPancake)
925 for i in range(len(theta)-1): # range is reduced as no brick can be created starting at the last point
926 # Assuming theta is a numpy array and you're looking for the index of a value close to pi
927 value_to_find = theta[i]+np.pi
928 tolerance = 1e-10 # Define a small tolerance
929 # Find indices where the condition is true
930 indices = np.where(np.abs(theta - value_to_find) < tolerance)[0]
931 if len(indices) > 0:
932 windingUpIndex = indices[0] # Take the first index if there are multiple matches
933 try:
934 x_1 = r[i] * np.cos(theta[i])
935 y_1 = r[i] * np.sin(theta[i])
936 z_g = z[i]
937 x_2 = r[i+1] * np.cos(theta[i+1])
938 y_2 = r[i+1] * np.sin(theta[i+1])
939 x_3 = r[windingUpIndex] * np.cos(theta[windingUpIndex])
940 y_3 = r[windingUpIndex] * np.sin(theta[windingUpIndex])
941 x_4 = r[windingUpIndex+1] * np.cos(theta[windingUpIndex+1])
942 y_4 = r[windingUpIndex+1] * np.sin(theta[windingUpIndex+1])
943 addPoints = [[x_1, y_1, z_g], [x_2, y_2, z_g], [x_3, y_3, z_g], [x_4, y_4, z_g]]
944 k = theta[i]//(2*np.pi)
945 if (theta[i] <= np.pi*(k+1)):
946 # print('winding 1 or 3')
947 if (k%2 == 0):
948 # print('winding 1')
949 winding_1.append(addPoints)
950 else:
951 # print('winding 3')
952 winding_3.append(addPoints)
954 if (theta[i] >= np.pi*(k+1)):
955 # print('winding 2 or 4')
956 if (k%2 == 0):
957 # print('winding 2')
958 winding_2.append(addPoints)
959 else:
960 # print('winding 4')
961 winding_4.append(addPoints)
962 except IndexError:
963 print('All of the winding conductor points have been found')
965 # print(winding_1)
966 # print(winding_2)
967 # print(winding_3)
968 # print(winding_4)
970 # Plotter
971 # x_coords = [coord[0] for coord in winding_1]
972 # y_coords = [coord[1] for coord in winding_1]
973 # z_coords = [coord[2] for coord in winding_1]
974 #
975 # # Creating the 3D plot
976 # fig = plt.figure()
977 # ax = fig.add_subplot(111, projection='3d')
978 #
979 # # Plotting the coordinates
980 # ax.scatter(x_coords, y_coords, z_coords, c='r', marker='o')
981 #
982 # # Setting labels
983 # ax.set_xlabel('X Label')
984 # ax.set_ylabel('Y Label')
985 # ax.set_zlabel('Z Label')
986 #
987 # plt.show()
989 # if True:
990 # indexPoint = 1
991 # rangeUpdated = 0
992 # dict_cond = {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'}, 1: {'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.861227', 'YP1': '-0.002972', 'ZP1': '-1.129183', 'XP2': '-0.861208', 'YP2': '0.002908', 'ZP2': '-1.129182', 'XP3': '-0.863294', 'YP3': '0.002912', 'ZP3': '-1.112210', 'XP4': '-0.863313', 'YP4': '-0.002968', 'ZP4': '-1.112211', 'XP5': '-0.842917', 'YP5': '-0.003066', 'ZP5': '-1.126858', 'XP6': '-0.842880', 'YP6': '0.002814', 'ZP6': '-1.126858', 'XP7': '-0.845242', 'YP7': '0.002830', 'ZP7': '-1.109922', 'XP8': '-0.845278', 'YP8': '-0.003050', 'ZP8': '-1.109922', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}, 2: {'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.842917', 'YP1': '-0.003066', 'ZP1': '-1.126858', 'XP2': '-0.842880', 'YP2': '0.002814', 'ZP2': '-1.126858', 'XP3': '-0.845242', 'YP3': '0.002830', 'ZP3': '-1.109922', 'XP4': '-0.845278', 'YP4': '-0.003050', 'ZP4': '-1.109922', 'XP5': '-0.824646', 'YP5': '-0.003216', 'ZP5': '-1.124235', 'XP6': '-0.824593', 'YP6': '0.002664', 'ZP6': '-1.124239', 'XP7': '-0.827229', 'YP7': '0.002698', 'ZP7': '-1.107343', 'XP8': '-0.827282', 'YP8': '-0.003181', 'ZP8': '-1.107339', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}}
993 # # print(dict_cond)
994 # for brick in dict_cond:
995 # for pointIndex in range (rangeUpdated, rangeUpdated+7):
996 # dict_cond[brick][f'XP{indexPoint}'] = str(coordinateList[pointIndex][0])
997 # dict_cond[brick][f'YP{indexPoint}'] = str(coordinateList[pointIndex][1])
998 # dict_cond[brick][f'ZP{indexPoint}'] = str(coordinateList[pointIndex][2])
999 # indexPoint+=1
1000 # indexPoint = 1
1001 # rangeUpdated = rangeUpdated + 8
1003 # writing COND.json file
1004 # Define the path for the JSON file, one directory up from the current script
1005 # json_file_path = os.path.join(os.path.dirname(os.getcwd()), "BR8.json")
1007 # Function to print the contents of a JSON file
1008 # def print_json_contents(path):
1009 # try:
1010 # with open(path, 'r') as file:
1011 # data = json.load(file)
1012 # print(json.dumps(data, indent=4))
1013 # except FileNotFoundError:
1014 # print("File not found.")
1015 # except json.JSONDecodeError:
1016 # print("File is empty or contains non-JSON conforming data.")
1017 #
1018 # # Print current contents
1019 # print("Current contents of BR8.json:")
1020 # print_json_contents(json_file_path)
1022 # Overwrite the JSON file
1023 # with open(json_file_path, 'w') as file:
1024 # json.dump(dict_cond, file, indent=4)
1025 #
1026 # print("\nContents of BR8.json after overwriting:")
1027 # print_json_contents(json_file_path)
1028 # writing the .cond file
1030 # p = tuh.Paths('tests/parsers', '')
1031 # FilesAndFolders.prep_folder(p.model_folder)
1033 # Specify the target directory relative to the current working directory
1034 target_dir = os.path.join(os.getcwd(), 'tests', '_outputs', 'parsers')
1036 # Ensure the target directory exists
1037 os.makedirs(target_dir, exist_ok=True)
1039 # Define the output file path
1040 # out_file = os.path.join(target_dir, 'BR8.cond')
1041 # list_of_shapes = ['BR8']
1042 # for shape in list_of_shapes:
1043 # pc = ParserCOND()
1044 # input_dict = dict_cond
1045 # pc.write_cond(input_dict, out_file)
1046 # print('path')
1047 # print(out_file)
1048 # print('hello world')
1049 # =============================================================================
1050 # GENERATING POINTS ENDS ======================================================
1051 # =============================================================================
1053 # =============================================================================
1054 # GENERATING SPLINES STARTS ===================================================
1055 # =============================================================================
1057 # Create the spline with the points:
1058 spline = gmsh.model.occ.addSpline(pointTags)
1060 # Split the spline into sub-curves:
1061 sectionsPerCurve = int(spt / cpt)
1063 # Create a list of point tags that will split the spline:
1064 # Basically, they are the points to divide the spirals into sectionsPerCurve
1065 # turns. However, some other points are also included to support the float
1066 # number of turns. It is best to visually look at the divisions with
1067 # gmsh.fltk.run() to understand why the split points are chosen the way they are
1068 # selected.
1070 if cutPlaneNormal is None:
1071 pointObjectsWithoutBreakPoints = points
1072 pointTagsWithoutBreakPoints = pointTags
1074 splitPointTags = list(
1075 set(pointTagsWithoutBreakPoints[:-1:sectionsPerCurve])
1076 | set(pointTagsWithoutBreakPoints[-spt - 1 :: -spt])
1077 | set(breakPointTagsDueToCutPlane)
1078 | set(breakPointTagsDueToTransition)
1079 )
1080 splitPointTags = sorted(splitPointTags)
1081 # Remove the first element of the list (starting point):
1082 _, *splitPointTags = splitPointTags
1084 # Also create a list of corresponding point objects:
1085 splitPoints = list(
1086 set(pointObjectsWithoutBreakPoints[:-1:sectionsPerCurve])
1087 | set(pointObjectsWithoutBreakPoints[-spt - 1 :: -spt])
1088 | set(breakPointObjectsDueToCutPlane)
1089 | set(breakPointObjectsDueToTransition)
1090 )
1091 splitPoints = sorted(splitPoints, key=lambda x: x.tag)
1092 # Remove the first element of the list (starting point):
1093 _, *splitPoints = splitPoints
1095 # Split the spline:
1096 dims = [0] * len(splitPointTags)
1097 _, splines = gmsh.model.occ.fragment(
1098 [(1, spline)],
1099 list(zip(dims, splitPointTags)),
1100 removeObject=True,
1101 removeTool=True,
1102 )
1103 splines = splines[0]
1104 self.splineTags = [j for _, j in splines]
1106 # Note the turn number of each spline. This will be used in getSplineTag and
1107 # getSplineTags methods.
1108 self.splineTurns = []
1109 for i in range(len(self.splineTags)):
1110 if i == 0:
1111 startPoint = points[0]
1112 endPoint = splitPoints[0]
1113 elif i == len(self.splineTags) - 1:
1114 startPoint = splitPoints[-1]
1115 endPoint = points[-1]
1116 else:
1117 startPoint = splitPoints[i - 1]
1118 endPoint = splitPoints[i]
1120 startTurn = (startPoint.theta - initialTheta) / (2 * math.pi)
1121 startTurn = round(startTurn / self.turnRes) * self.turnRes
1122 endTurn = (endPoint.theta - initialTheta) / (2 * math.pi)
1123 endTurn = round(endTurn / self.turnRes) * self.turnRes
1125 if direction is direction.ccw:
1126 self.splineTurns.append((startTurn, endTurn))
1127 else:
1128 self.splineTurns.append((-startTurn, -endTurn))
1130 # Check if splineTurn tuples starts with the small turn number:
1131 for i in range(len(self.splineTurns)):
1132 self.splineTurns[i] = sorted(self.splineTurns[i])
1134 # =============================================================================
1135 # GENERATING SPLINES ENDS =====================================================
1136 # =============================================================================
1138 # Find start and end points of the spiral curve:
1139 gmsh.model.occ.synchronize() # synchronize the model to make getBoundary work
1140 self.startPointTag = gmsh.model.getBoundary([(1, self.getSplineTag(0))])[1][1]
1141 self.endPointTag = gmsh.model.getBoundary(
1142 [(1, self.getSplineTag(self.turns, endPoint=True))]
1143 )[1][1]
1145 def getSplineTag(self, turn, endPoint=False):
1146 """
1147 Returns the spline tag at a specific turn. It returns the spline tag of the
1148 section that is on the turn except its end point.
1150 :param turn: turn number (it can be a float)
1151 :type turn: float
1152 :param endPoint: if True, return the spline tag of the section that is on the
1153 turn including its end point but not its start point (default: False)
1154 :type endPoint: bool, optional
1155 :return: spline tag
1156 """
1157 if endPoint:
1158 for i, r in enumerate(self.splineTurns):
1159 if r[0] + (self.turnRes / 2) < turn <= r[1] + (self.turnRes / 2):
1160 return self.splineTags[i]
1161 else:
1162 for i, r in enumerate(self.splineTurns):
1163 if r[0] - (self.turnRes / 2) <= turn < r[1] - (self.turnRes / 2):
1164 return self.splineTags[i]
1166 def getPointTag(self, turn, endPoint=False):
1167 """
1168 Returns the point object at a specific turn.
1170 :param turn: turn number (it can be a float)
1171 :type turn: float
1172 :return: point object
1173 :rtype: point
1174 """
1175 if turn < 0 or turn > self.turns:
1176 raise ValueError("Turn number is out of range!")
1178 if turn == 0:
1179 return self.startPointTag
1180 elif turn == self.turns:
1181 return self.endPointTag
1182 else:
1183 curveTag = self.getSplineTag(turn, endPoint=endPoint)
1184 if endPoint:
1185 points = gmsh.model.getBoundary([(1, curveTag)])
1186 return points[1][1]
1187 else:
1188 points = gmsh.model.getBoundary([(1, curveTag)])
1189 return points[0][1]
1191 def getSplineTags(self, turnStart=None, turnEnd=None):
1192 """
1193 Get the spline tags from a specific turn to another specific turn. If turnStart
1194 and turnEnd are not specified, it returns all the spline tags.
1196 :param turnStart: start turn number (it can be a float) (default: None)
1197 :type turnStart: float, optional
1198 :param turnEnd: end turn number (it can be a float) (default: None)
1199 :type turnEnd: float, optional
1200 :return: spline tags
1201 :rtype: list[int]
1202 """
1203 if turnStart is None and turnEnd is None:
1204 return self.splineTags
1205 elif turnStart is None or turnEnd is None:
1206 raise ValueError(
1207 "turnStart and turnEnd must be both specified or both not specified."
1208 " You specified only one of them."
1209 )
1210 else:
1211 start = self.splineTags.index(self.getSplineTag(turnStart, False))
1212 end = self.splineTags.index(self.getSplineTag(turnEnd, True)) + 1
1213 return self.splineTags[start:end]
1216class spiralSurface:
1217 """
1218 This is a class to create a spiral surface parallel to the XY plane in GMSH. If
1219 thinShellApproximation is set to False, it creates two spiral surfaces parallel to
1220 the XY plane, and their inner and outer curve loops in GMSH. One of the surfaces is
1221 the main surface specified, which is the winding surface, and the other is the
1222 contact layer surface (the gap between the winding surface). If thinShellApproximation
1223 is set to True, it creates only one spiral surface that touches each other
1224 (conformal).
1226 Note that surfaces are subdivided depending on the spiral curve divisions because
1227 this is required for the thin-shell approximation. Otherwise, when
1228 thinShellApproximation is set to True, it would be a disc rather than a spiral since
1229 it touches each other. However, this can be avoided by dividing the surfaces into
1230 small parts and making them conformal. Dividing the surfaces is not necessary when
1231 thinShellApproximation is set to false, but in order to use the same logic with TSA,
1232 it is divided anyway.
1234 :param innerRadius: inner radius
1235 :type innerRadius: float
1236 :param thickness: thickness
1237 :type thickness: float
1238 :param contactLayerThickness: contact layer thickness
1239 :type contactLayerThickness: float
1240 :param turns: number of turns
1241 :type turns: float
1242 :param z: z coordinate
1243 :type z: float
1244 :param initialTheta: initial theta angle in radians
1245 :type initialTheta: float
1246 :param spiralDirection: direction of the spiral (default: direction.ccw)
1247 :type spiralDirection: direction, optional
1248 :param thinShellApproximation: if True, the thin shell approximation is used
1249 (default: False)
1250 :type thinShellApproximation: bool, optional
1251 :param cutPlaneNormal: normal vector of the plane that will cut the spiral surface
1252 (default: None)
1253 :type cutPlaneNormal: tuple[float, float, float], optional
1254 """
1256 def __init__(
1257 self,
1258 innerRadius,
1259 thickness,
1260 contactLayerThickness,
1261 turns,
1262 z,
1263 initialTheta,
1264 transitionNotchAngle,
1265 spiralDirection=direction.ccw,
1266 thinShellApproximation=False,
1267 cutPlaneNormal=None,
1268 ) -> None:
1269 r_i = innerRadius
1270 t = thickness
1271 theta_i = initialTheta
1272 self.theta_i = theta_i
1274 self.surfaceTags = []
1275 self.contactLayerSurfaceTags = []
1277 self.direction = spiralDirection
1278 self.tsa = thinShellApproximation
1279 self.transitionNotchAngle = transitionNotchAngle
1280 # =============================================================================
1281 # GENERATING SPIRAL CURVES STARTS =============================================
1282 # =============================================================================
1283 if thinShellApproximation:
1284 # Create only one spiral curve because the thin shell approximation is used:
1285 # Winding thickness is increased slightly to ensure that the outer radius
1286 # would be the same without the thin shell approximation.
1287 turns = (
1288 turns + 1
1289 ) # for TSA, spiral has (n+1) turns but spiral surface has n
1290 spiral = spiralCurve(
1291 r_i,
1292 t + contactLayerThickness * (turns - 1) / (turns),
1293 turns,
1294 z,
1295 theta_i,
1296 transitionNotchAngle,
1297 spiralDirection,
1298 cutPlaneNormal=cutPlaneNormal,
1299 )
1301 # These are created to be able to use the same code with TSA and without TSA:
1302 innerSpiral = spiral
1303 outerSpiral = spiral
1304 else:
1305 # Create two spiral curves because the thin shell approximation is not used:
1306 innerSpiral = spiralCurve(
1307 r_i - contactLayerThickness,
1308 t + contactLayerThickness,
1309 turns + 1,
1310 z,
1311 theta_i,
1312 transitionNotchAngle,
1313 spiralDirection,
1314 cutPlaneNormal=cutPlaneNormal,
1315 )
1316 outerSpiral = spiralCurve(
1317 r_i,
1318 t + contactLayerThickness,
1319 turns + 1,
1320 z,
1321 theta_i,
1322 transitionNotchAngle,
1323 spiralDirection,
1324 cutPlaneNormal=cutPlaneNormal,
1325 )
1327 self.innerSpiral = innerSpiral
1328 self.outerSpiral = outerSpiral
1329 self.turns = turns
1330 # =============================================================================
1331 # GENERATING SPIRAL CURVES ENDS ===============================================
1332 # =============================================================================
1334 # =============================================================================
1335 # GENERATING SURFACES STARTS ==================================================
1336 # =============================================================================
1337 endLines = []
1338 endInsLines = []
1340 # This is used to check if all the contact layers are finished:
1341 allContactLayersAreFinished = False
1343 # Itterate over the spline tags:
1344 for i in range(len(innerSpiral.splineTags)):
1345 if thinShellApproximation:
1346 # The current spline will be the inner spline:
1347 innerSplineTag = spiral.splineTags[i]
1349 # Find the spline tag of the outer spline by finding the spline tag of
1350 # the spline that is exactly on the next turn:
1352 # Note the turn number of the current spline's start point:
1353 startTurn = spiral.splineTurns[i][0]
1355 if startTurn + 1 + 1e-4 > turns:
1356 # If the current spline is on the outer surface, break the loop,
1357 # because the whole surface is finished:
1358 break
1360 # Find the outer spline tag:
1361 isItBroken = True
1362 for j, turnTuple in enumerate(spiral.splineTurns):
1363 # Equality can not be checked with == because of the floating point
1364 # errors:
1365 if abs(turnTuple[0] - 1 - startTurn) < 1e-4:
1366 outerSplineTag = spiral.splineTags[j]
1367 isItBroken = False
1368 break
1370 if isItBroken:
1371 raise RuntimeError(
1372 "Something went wrong while creating the spiral surface. Outer"
1373 f" spline tag of {innerSplineTag} could not be found for TSA."
1374 )
1376 else:
1377 # Store the tags of the current splines:
1378 innerSplineTag = innerSpiral.splineTags[i]
1379 outerSplineTag = outerSpiral.splineTags[i]
1381 # The current outer spline will be the inner spline of the
1382 # contact layer:
1383 innerInsSplineTag = outerSpiral.splineTags[i]
1385 # Find the spline tag of the contact layer's outer spline by finding the
1386 # spline tag of the spline that is exactly on the next turn:
1388 # Note the turn number of the current spline's start point:
1389 startTurn = outerSpiral.splineTurns[i][0]
1391 if startTurn + 1 + 1e-4 > turns + 1:
1392 # If the current spline is on the outer surface, note that all the
1393 # contact layers are finished:
1394 allContactLayersAreFinished = True
1396 # Find the contact layer's outer spline tag:
1397 for j, turnTuple in enumerate(innerSpiral.splineTurns):
1398 if math.isclose(turnTuple[0], 1 + startTurn, abs_tol=1e-6):
1399 outerInsSplineTag = innerSpiral.splineTags[j]
1400 break
1402 # Create the lines to connect the two splines so that a surface can be
1403 # created:
1405 # Create start line:
1406 if i == 0:
1407 # If it is the first spline, start line should be created.
1409 # Create points:
1410 isStartPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]
1411 if thinShellApproximation:
1412 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]
1413 else:
1414 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]
1416 # Create the line:
1417 startLine = gmsh.model.occ.addLine(osStartPoint, isStartPoint)
1418 firstStartLine = startLine
1420 # Create lines for the contact layer if the thin shell approximation is not
1421 # used:
1422 if not thinShellApproximation and not allContactLayersAreFinished:
1423 isInsStartPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[
1424 1
1425 ][1]
1426 osInsStartPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[
1427 0
1428 ][1]
1430 # Create the line:
1431 startInsLine = gmsh.model.occ.addLine(
1432 osInsStartPoint, isInsStartPoint
1433 )
1434 firstInsStartLine = startInsLine
1436 else:
1437 # If it is not the first spline, the start line is the end line of the
1438 # previous surface. This guarantees that the surfaces are connected
1439 # (conformality).
1440 startLine = endLines[i - 1]
1442 # Do the same for the contact layer if the thin shell approximation is not
1443 # used:
1444 if not thinShellApproximation and not allContactLayersAreFinished:
1445 startInsLine = endInsLines[i - 1]
1447 # Create end line:
1449 # Create points:
1450 # The ifs are used because getBoundary is not consistent somehow.
1451 if i == 0:
1452 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[0][1]
1453 else:
1454 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]
1456 if (not i == 0) or thinShellApproximation:
1457 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]
1458 else:
1459 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]
1461 # Create the line:
1462 endLine = gmsh.model.occ.addLine(isEndPoint, osEndPoint)
1463 endLines.append(endLine)
1465 # Create lines for the contact layer if the thin shell approximation is not
1466 # used:
1467 if not thinShellApproximation and not allContactLayersAreFinished:
1468 if i == 0:
1469 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[0][
1470 1
1471 ]
1472 else:
1473 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[1][
1474 1
1475 ]
1477 osInsEndPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[1][1]
1479 # Create the line:
1480 endInsLine = gmsh.model.occ.addLine(isInsEndPoint, osInsEndPoint)
1481 endInsLines.append(endInsLine)
1483 # Create the surface:
1484 curveLoop = gmsh.model.occ.addCurveLoop(
1485 [startLine, innerSplineTag, endLine, outerSplineTag]
1486 )
1487 self.surfaceTags.append(gmsh.model.occ.addPlaneSurface([curveLoop]))
1489 # Create the surface for the contact layer if the thin shell approximation is
1490 # not used:
1491 if not thinShellApproximation and not allContactLayersAreFinished:
1492 curveLoop = gmsh.model.occ.addCurveLoop(
1493 [startInsLine, innerInsSplineTag, endInsLine, outerInsSplineTag]
1494 )
1495 self.contactLayerSurfaceTags.append(
1496 gmsh.model.occ.addPlaneSurface([curveLoop])
1497 )
1499 # =============================================================================
1500 # GENERATING SURFACES ENDS ====================================================
1501 # =============================================================================
1503 # =============================================================================
1504 # GENERATING CURVE LOOPS STARTS ===============================================
1505 # =============================================================================
1507 # Create the inner and outer curve loops (for both TSA == True and TSA == False):
1509 # VERY IMPORTANT NOTES ABOUT THE DIRECTION OF THE CURVE LOOPS
1510 # 1- GMSH doesn't like duplicates. Or the user doesn't like duplicates if they
1511 # want conformality. Actually, it's a positive thing about debugging because
1512 # you can always use `Geometry.remove_all_duplicates()` to see if there are
1513 # any duplicates. If there are, the problem is found. Solve it.
1514 # 2- The problem arises when one uses surface loops or curve loops. Because even
1515 # if you think there are no duplicates, GMSH/OCC might create some during
1516 # addCurveLoops and addSurfaceLoops operations. Even though
1517 # `geometry.remove_all_duplicates()` tells that there are duplicates, the
1518 # user doesn't suspect about addCurveLoop and addSurfaceLoop at first,
1519 # because it doesn't make sense.
1520 # 3- How you put the curves in the curve loops is very important! The same curve
1521 # loop with the same lines might cause problems if the user puts them in a
1522 # different order. For example, to create a plane surface with two curve
1523 # loops, the direction of the curve loops should be the same. That's why the
1524 # code has both innerCurveLoopTag and innerOppositeCurveLoopTag (the same
1525 # thing for the outer curve loop).
1527 # create the transition layer (notch):
1528 # Inner curve loop:
1529 notchStartPoint = innerSpiral.getPointTag(
1530 1 - transitionNotchAngle / (2 * math.pi)
1531 )
1532 notchLeftPoint = innerSpiral.getPointTag(0)
1533 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1534 notchRightLine = innerSpiral.getSplineTag(
1535 1 - transitionNotchAngle / (2 * math.pi)
1536 )
1538 if thinShellApproximation:
1539 innerStartCurves = [firstStartLine]
1540 else:
1541 innerStartCurves = [firstInsStartLine, firstStartLine]
1543 if thinShellApproximation:
1545 notchCurveLoop = gmsh.model.occ.addCurveLoop(
1546 [notchLeftLine, notchRightLine] + innerStartCurves
1547 )
1548 self.innerNotchSurfaceTags = [
1549 gmsh.model.occ.addPlaneSurface([notchCurveLoop])
1550 ]
1551 else:
1552 notchMiddlePoint = outerSpiral.getPointTag(0)
1553 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)
1555 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1556 [notchLeftLine, notchMiddleLine, firstStartLine]
1557 )
1558 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(
1559 [notchMiddleLine, notchRightLine, firstInsStartLine]
1560 )
1561 self.innerNotchSurfaceTags = [
1562 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1563 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),
1564 ]
1566 lines = innerSpiral.getSplineTags(
1567 0, 1 - transitionNotchAngle / (2 * math.pi)
1568 ) # The first turn of the spline
1569 innerCurves = lines + [notchLeftLine]
1570 self.innerNotchLeftLine = notchLeftLine
1572 self.innerStartCurves = innerStartCurves
1574 self.innerCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves)
1575 self.innerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves[::-1])
1577 # Outer curve loop:
1578 # The last turn of the spline:
1579 if thinShellApproximation:
1580 notchStartPoint = innerSpiral.getPointTag(
1581 self.turns + transitionNotchAngle / (2 * math.pi) - 1
1582 )
1583 notchLeftPoint = innerSpiral.getPointTag(self.turns)
1584 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1585 notchRightLine = innerSpiral.getSplineTag(
1586 self.turns - 1 + transitionNotchAngle / (2 * math.pi), endPoint=True
1587 )
1588 else:
1589 notchStartPoint = outerSpiral.getPointTag(
1590 self.turns + transitionNotchAngle / (2 * math.pi)
1591 )
1592 notchMiddlePoint = innerSpiral.getPointTag(self.turns + 1)
1593 notchLeftPoint = outerSpiral.getPointTag(self.turns + 1)
1594 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1595 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)
1596 notchRightLine = outerSpiral.getSplineTag(
1597 self.turns + transitionNotchAngle / (2 * math.pi), self.turns
1598 )
1599 if thinShellApproximation:
1600 lines = outerSpiral.getSplineTags(turns - 1, turns)
1601 else:
1602 lines = outerSpiral.getSplineTags(turns, turns + 1)
1604 if thinShellApproximation:
1605 outerEndCurves = [endLines[-1]]
1606 else:
1607 outerEndCurves = [endInsLines[-1], endLines[-1]]
1609 if thinShellApproximation:
1610 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1611 [notchLeftLine, notchRightLine, endLines[-1]]
1612 )
1613 self.outerNotchSurfaceTags = [
1614 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1615 ]
1616 else:
1617 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1618 [notchLeftLine, notchMiddleLine, endLines[-1]]
1619 )
1620 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(
1621 [notchMiddleLine, notchRightLine, endInsLines[-1]]
1622 )
1623 self.outerNotchSurfaceTags = [
1624 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1625 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),
1626 ]
1628 if thinShellApproximation:
1629 lines = innerSpiral.getSplineTags(
1630 self.turns - 1 + transitionNotchAngle / (2 * math.pi), self.turns
1631 ) # The first turn of the spline
1632 else:
1633 lines = outerSpiral.getSplineTags(
1634 self.turns + transitionNotchAngle / (2 * math.pi), self.turns + 1
1635 )
1636 outerCurves = lines + [notchLeftLine]
1637 self.outerNotchLeftLine = notchLeftLine
1639 self.outerEndCurves = outerEndCurves
1641 self.outerCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves)
1642 self.outerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves[::-1])
1643 # =============================================================================
1644 # GENERATING CURVE LOOPS ENDS =================================================
1645 # =============================================================================
1647 if not thinShellApproximation:
1648 surfaceTags = self.surfaceTags
1649 self.surfaceTags = self.contactLayerSurfaceTags
1650 self.contactLayerSurfaceTags = surfaceTags
1652 def getInnerRightPointTag(self):
1653 """
1654 Returns the point tag of the inner left point.
1656 :return: point tag
1657 :rtype: int
1658 """
1659 return self.innerSpiral.getPointTag(0)
1661 def getInnerUpperPointTag(self):
1662 """
1663 Returns the point tag of the inner right point.
1665 :return: point tag
1666 :rtype: int
1667 """
1668 if self.direction is direction.ccw:
1669 return self.innerSpiral.getPointTag(0.25)
1670 else:
1671 return self.innerSpiral.getPointTag(0.75)
1673 def getInnerLeftPointTag(self):
1674 """
1675 Returns the point tag of the inner upper point.
1677 :return: point tag
1678 :rtype: int
1679 """
1680 return self.innerSpiral.getPointTag(0.5)
1682 def getInnerLowerPointTag(self):
1683 """
1684 Returns the point tag of the inner lower point.
1686 :return: point tag
1687 :rtype: int
1688 """
1689 if self.direction is direction.ccw:
1690 return self.innerSpiral.getPointTag(0.75)
1691 else:
1692 return self.innerSpiral.getPointTag(0.25)
1694 def getOuterRightPointTag(self):
1695 """
1696 Returns the point tag of the outer left point.
1698 :return: point tag
1699 :rtype: int
1700 """
1701 if self.tsa:
1702 turns = self.turns
1703 else:
1704 turns = self.turns + 1
1705 return self.outerSpiral.getPointTag(turns, endPoint=False)
1707 def getOuterLowerPointTag(self):
1708 """
1709 Returns the point tag of the outer right point.
1711 :return: point tag
1712 :rtype: int
1713 """
1714 if self.tsa:
1715 turns = self.turns
1716 else:
1717 turns = self.turns + 1
1718 if self.direction is direction.ccw:
1719 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)
1720 else:
1721 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)
1723 def getOuterLeftPointTag(self):
1724 """
1725 Returns the point tag of the outer upper point.
1727 :return: point tag
1728 :rtype: int
1729 """
1730 if self.tsa:
1731 turns = self.turns
1732 else:
1733 turns = self.turns + 1
1734 return self.outerSpiral.getPointTag(turns - 0.5, endPoint=False)
1736 def getOuterUpperPointTag(self):
1737 """
1738 Returns the point tag of the outer lower point.
1740 :return: point tag
1741 :rtype: int
1742 """
1743 if self.tsa:
1744 turns = self.turns
1745 else:
1746 turns = self.turns + 1
1747 if self.direction is direction.ccw:
1748 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)
1749 else:
1750 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)
1752 def getInnerUpperRightCurves(self):
1753 """
1754 Returns the curve tags of the upper right curves.
1756 :return: curve tags
1757 :rtype: list[int]
1758 """
1759 if self.direction is direction.ccw:
1760 curves = self.innerSpiral.getSplineTags(0, 0.25)
1761 else:
1762 lines = self.innerSpiral.getSplineTags(
1763 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)
1764 ) # The first turn of the spline
1765 lines = lines + [self.innerNotchLeftLine]
1767 return lines
1769 return curves
1771 def getInnerUpperLeftCurves(self):
1772 """
1773 Returns the curve tags of the upper left curves.
1775 :return: curve tags
1776 :rtype: list[int]
1777 """
1778 if self.direction is direction.ccw:
1779 curves = self.innerSpiral.getSplineTags(0.25, 0.5)
1780 else:
1781 curves = self.innerSpiral.getSplineTags(0.5, 0.75)
1783 return curves
1785 def getInnerLowerLeftCurves(self):
1786 """
1787 Returns the curve tags of the lower left curves.
1789 :return: curve tags
1790 :rtype: list[int]
1791 """
1792 if self.direction is direction.ccw:
1793 curves = self.innerSpiral.getSplineTags(0.5, 0.75)
1794 else:
1795 curves = self.innerSpiral.getSplineTags(0.25, 0.5)
1797 return curves
1799 def getInnerLowerRightCurves(self):
1800 """
1801 Returns the curve tags of the lower right curves.
1803 :return: curve tags
1804 :rtype: list[int]
1805 """
1806 if self.direction is direction.ccw:
1807 lines = self.innerSpiral.getSplineTags(
1808 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)
1809 ) # The first turn of the spline
1810 lines = lines + [self.innerNotchLeftLine]
1812 return lines
1813 else:
1814 curves = self.innerSpiral.getSplineTags(0, 0.25)
1816 return curves
1818 def getOuterUpperRightCurves(self):
1819 """
1820 Returns the curve tags of the upper right curves.
1822 :return: curve tags
1823 :rtype: list[int]
1824 """
1825 if self.tsa:
1826 turns = self.turns
1827 else:
1828 turns = self.turns + 1
1830 if self.direction is direction.ccw:
1831 if self.tsa:
1832 lines = self.innerSpiral.getSplineTags(
1833 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),
1834 self.turns - 0.75,
1835 ) # The first turn of the spline
1836 else:
1837 lines = self.outerSpiral.getSplineTags(
1838 self.turns + self.transitionNotchAngle / (2 * math.pi),
1839 self.turns + 1 - 0.75,
1840 )
1841 lines = lines + [self.outerNotchLeftLine]
1843 return lines
1844 else:
1845 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)
1847 return curves
1849 def getOuterUpperLeftCurves(self):
1850 """
1851 Returns the curve tags of the lower right curves.
1853 :return: curve tags
1854 :rtype: list[int]
1855 """
1856 if self.tsa:
1857 turns = self.turns
1858 else:
1859 turns = self.turns + 1
1860 if self.direction is direction.ccw:
1861 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)
1862 else:
1863 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)
1865 return curves
1867 def getOuterLowerLeftCurves(self):
1868 """
1869 Returns the curve tags of the lower left curves.
1871 :return: curve tags
1872 :rtype: list[int]
1873 """
1874 if self.tsa:
1875 turns = self.turns
1876 else:
1877 turns = self.turns + 1
1878 if self.direction is direction.ccw:
1879 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)
1880 else:
1881 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)
1883 return curves
1885 def getOuterLowerRightCurves(self):
1886 """
1887 Returns the curve tags of the upper left curves.
1889 :return: curve tags
1890 :rtype: list[int]
1891 """
1892 if self.tsa:
1893 turns = self.turns
1894 else:
1895 turns = self.turns + 1
1896 if self.direction is direction.ccw:
1897 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)
1898 else:
1899 if self.tsa:
1900 lines = self.innerSpiral.getSplineTags(
1901 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),
1902 self.turns - 0.75,
1903 ) # The first turn of the spline
1904 else:
1905 lines = self.outerSpiral.getSplineTags(
1906 self.turns + self.transitionNotchAngle / (2 * math.pi),
1907 self.turns + 1 - 0.75,
1908 )
1909 lines = lines + [self.outerNotchLeftLine]
1911 return lines
1913 return curves
1915 def getInnerStartCurves(self):
1916 """
1917 Returns the curve tags of the start curves.
1919 :return: curve tags
1920 :rtype: list[int]
1921 """
1922 return self.innerStartCurves
1924 def getOuterEndCurves(self):
1925 """
1926 Returns the curve tags of the end curves.
1928 :return: curve tags
1929 :rtype: list[int]
1930 """
1931 return self.outerEndCurves
1934class circleWithFourCurves:
1935 def __init__(
1936 self,
1937 x,
1938 y,
1939 z,
1940 radius,
1941 upperRightTag=None,
1942 upperLeftTag=None,
1943 lowerLeftTag=None,
1944 lowerRightTag=None,
1945 leftPointTag=None,
1946 rightPointTag=None,
1947 upperPointTag=None,
1948 lowerPointTag=None,
1949 ):
1950 self.x = x
1951 self.y = y
1952 self.z = z
1953 self.r = radius
1954 if upperRightTag is None:
1955 dummyCircle = gmsh.model.occ.addCircle(self.x, self.y, self.z, self.r)
1956 self.leftPointTag = gmsh.model.occ.addPoint(self.x - self.r, self.y, self.z)
1957 self.rightPointTag = gmsh.model.occ.addPoint(
1958 self.x + self.r, self.y, self.z
1959 )
1960 self.upperPointTag = gmsh.model.occ.addPoint(
1961 self.x, self.y + self.r, self.z
1962 )
1963 self.lowerPointTag = gmsh.model.occ.addPoint(
1964 self.x, self.y - self.r, self.z
1965 )
1967 fragmentResults = gmsh.model.occ.fragment(
1968 [(1, dummyCircle)],
1969 [
1970 (0, self.leftPointTag),
1971 (0, self.rightPointTag),
1972 (0, self.upperPointTag),
1973 (0, self.lowerPointTag),
1974 ],
1975 )[0]
1976 linesDimTags = [dimTag for dimTag in fragmentResults if dimTag[0] == 1]
1978 self.upperRightTag = linesDimTags[0][1]
1979 self.upperLeftTag = linesDimTags[1][1]
1980 self.lowerLeftTag = linesDimTags[2][1]
1981 self.lowerRightTag = linesDimTags[3][1]
1982 else:
1983 self.upperRightTag = upperRightTag
1984 self.upperLeftTag = upperLeftTag
1985 self.lowerLeftTag = lowerLeftTag
1986 self.lowerRightTag = lowerRightTag
1988 self.leftPointTag = leftPointTag
1989 self.rightPointTag = rightPointTag
1990 self.upperPointTag = upperPointTag
1991 self.lowerPointTag = lowerPointTag
1994class outerAirSurface:
1995 def __init__(
1996 self,
1997 outerRadius,
1998 innerRadius,
1999 type="cylinder",
2000 divideIntoFourParts=False,
2001 divideTerminalPartIntoFourParts=False,
2002 ):
2003 self.surfaceTags = []
2005 self.divideIntoFourParts = divideIntoFourParts
2006 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts
2008 # for cylinder:
2009 self.shellTags = []
2011 # for cuboid:
2012 self.shellTagsPart1 = []
2013 self.shellTagsPart2 = []
2014 self.shellTagsPart3 = []
2015 self.shellTagsPart4 = []
2017 self.type = type
2018 self.outerRadius = outerRadius
2019 self.innerRadius = innerRadius
2021 def createFromScratch(self, z, shellTransformation=False, shellRadius=None):
2022 self.z = z
2024 if self.divideIntoFourParts:
2025 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)
2026 else:
2027 if self.divideTerminalPartIntoFourParts:
2028 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)
2029 innerCL = gmsh.model.occ.addCurveLoop(
2030 [
2031 self.innerCircle.upperRightTag,
2032 self.innerCircle.upperLeftTag,
2033 self.innerCircle.lowerLeftTag,
2034 self.innerCircle.lowerRightTag,
2035 ]
2036 )
2037 else:
2038 innerCL = gmsh.model.occ.addCircle(0, 0, z, self.innerRadius)
2039 innerCL = gmsh.model.occ.addCurveLoop([innerCL])
2041 if self.type == "cylinder" and self.divideIntoFourParts:
2042 outerCircle = circleWithFourCurves(0, 0, z, self.outerRadius)
2044 leftLineTag = gmsh.model.occ.addLine(
2045 outerCircle.leftPointTag, self.innerCircle.leftPointTag
2046 )
2047 rightLineTag = gmsh.model.occ.addLine(
2048 outerCircle.rightPointTag, self.innerCircle.rightPointTag
2049 )
2050 upperLineTag = gmsh.model.occ.addLine(
2051 outerCircle.upperPointTag, self.innerCircle.upperPointTag
2052 )
2053 lowerLineTag = gmsh.model.occ.addLine(
2054 outerCircle.lowerPointTag, self.innerCircle.lowerPointTag
2055 )
2057 # Create surfaces:
2058 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2059 [
2060 outerCircle.upperRightTag,
2061 rightLineTag,
2062 self.innerCircle.upperRightTag,
2063 -upperLineTag,
2064 ]
2065 )
2066 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2067 self.surfaceTags.append(self.upperRightTag)
2069 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2070 [
2071 outerCircle.upperLeftTag,
2072 leftLineTag,
2073 self.innerCircle.upperLeftTag,
2074 -upperLineTag,
2075 ]
2076 )
2077 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2078 self.surfaceTags.append(self.upperLeftTag)
2080 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2081 [
2082 outerCircle.lowerLeftTag,
2083 leftLineTag,
2084 self.innerCircle.lowerLeftTag,
2085 -lowerLineTag,
2086 ]
2087 )
2088 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2089 self.surfaceTags.append(self.lowerLeftTag)
2091 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2092 [
2093 outerCircle.lowerRightTag,
2094 rightLineTag,
2095 self.innerCircle.lowerRightTag,
2096 -lowerLineTag,
2097 ]
2098 )
2099 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2100 self.surfaceTags.append(self.lowerRightTag)
2102 if shellTransformation:
2103 shellOuterCircle = circleWithFourCurves(0, 0, z, shellRadius)
2104 shellLeftLineTag = gmsh.model.occ.addLine(
2105 shellOuterCircle.leftPointTag, outerCircle.leftPointTag
2106 )
2107 shellRightLineTag = gmsh.model.occ.addLine(
2108 shellOuterCircle.rightPointTag, outerCircle.rightPointTag
2109 )
2110 shellUpperLineTag = gmsh.model.occ.addLine(
2111 shellOuterCircle.upperPointTag, outerCircle.upperPointTag
2112 )
2113 shellLowerLineTag = gmsh.model.occ.addLine(
2114 shellOuterCircle.lowerPointTag, outerCircle.lowerPointTag
2115 )
2117 # Create surfaces:
2118 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2119 [
2120 shellOuterCircle.upperRightTag,
2121 shellRightLineTag,
2122 outerCircle.upperRightTag,
2123 -shellUpperLineTag,
2124 ]
2125 )
2126 self.upperRightTag = gmsh.model.occ.addPlaneSurface(
2127 [upperRightCurveLoop]
2128 )
2129 self.shellTags.append(self.upperRightTag)
2131 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2132 [
2133 shellOuterCircle.upperLeftTag,
2134 shellLeftLineTag,
2135 outerCircle.upperLeftTag,
2136 -shellUpperLineTag,
2137 ]
2138 )
2139 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2140 self.shellTags.append(self.upperLeftTag)
2142 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2143 [
2144 shellOuterCircle.lowerLeftTag,
2145 shellLeftLineTag,
2146 outerCircle.lowerLeftTag,
2147 -shellLowerLineTag,
2148 ]
2149 )
2150 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2151 self.shellTags.append(self.lowerLeftTag)
2153 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2154 [
2155 shellOuterCircle.lowerRightTag,
2156 shellRightLineTag,
2157 outerCircle.lowerRightTag,
2158 -shellLowerLineTag,
2159 ]
2160 )
2161 self.lowerRightTag = gmsh.model.occ.addPlaneSurface(
2162 [lowerRightCurveLoop]
2163 )
2164 self.shellTags.append(self.lowerRightTag)
2166 elif self.type == "cylinder" and not self.divideIntoFourParts:
2167 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.outerRadius)
2168 outerCL = gmsh.model.occ.addCurveLoop([outerCL])
2170 if shellTransformation:
2171 shellOuterCL = gmsh.model.occ.addCircle(0, 0, z, shellRadius)
2172 shellOuterCL = gmsh.model.occ.addCurveLoop([shellOuterCL])
2174 shellSurfaceTag = gmsh.model.occ.addPlaneSurface(
2175 [shellOuterCL, outerCL]
2176 )
2177 self.shellTags.append(shellSurfaceTag)
2179 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])
2180 self.surfaceTags.append(surfaceTag)
2182 elif self.type == "cuboid":
2183 # LL: lower left
2184 # LR: lower right
2185 # UR: upper right
2186 # UL: upper left
2187 airLLpointTag = point(-self.outerRadius, -self.outerRadius, z).tag
2188 airLRpointTag = point(self.outerRadius, -self.outerRadius, z).tag
2189 airURpointTag = point(self.outerRadius, self.outerRadius, z).tag
2190 airULpointTag = point(-self.outerRadius, self.outerRadius, z).tag
2192 # LH: lower horizontal
2193 # UH: upper horizontal
2194 # LV: left vertical
2195 # RV: right vertical
2196 airLHlineTag = gmsh.model.occ.addLine(airLLpointTag, airLRpointTag)
2197 airRVLineTag = gmsh.model.occ.addLine(airLRpointTag, airURpointTag)
2198 airUHLineTag = gmsh.model.occ.addLine(airURpointTag, airULpointTag)
2199 airLVLineTag = gmsh.model.occ.addLine(airULpointTag, airLLpointTag)
2201 outerCL = gmsh.model.occ.addCurveLoop(
2202 [airLHlineTag, airRVLineTag, airUHLineTag, airLVLineTag]
2203 )
2205 if self.divideIntoFourParts:
2206 innerCL = gmsh.model.occ.addCurveLoop(
2207 [
2208 self.innerCircle.upperRightTag,
2209 self.innerCircle.lowerRightTag,
2210 self.innerCircle.lowerLeftTag,
2211 self.innerCircle.upperLeftTag,
2212 ]
2213 )
2215 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])
2216 self.surfaceTags.append(surfaceTag)
2218 if shellTransformation:
2219 # LL: lower left
2220 # LR: lower right
2221 # UR: upper right
2222 # UL: upper left
2223 shellLLpointTag = point(
2224 -shellRadius,
2225 -shellRadius,
2226 z,
2227 ).tag
2228 shellLRpointTag = point(
2229 shellRadius,
2230 -shellRadius,
2231 z,
2232 ).tag
2233 shellURpointTag = point(
2234 shellRadius,
2235 shellRadius,
2236 z,
2237 ).tag
2238 shellULpointTag = point(
2239 -shellRadius,
2240 shellRadius,
2241 z,
2242 ).tag
2244 # LH: lower horizontal
2245 # UH: upper horizontal
2246 # LV: left vertical
2247 # RV: right vertical
2248 shellLHlineTag = gmsh.model.occ.addLine(
2249 shellLLpointTag, shellLRpointTag
2250 )
2251 shellRVLineTag = gmsh.model.occ.addLine(
2252 shellLRpointTag, shellURpointTag
2253 )
2254 shellUHLineTag = gmsh.model.occ.addLine(
2255 shellURpointTag, shellULpointTag
2256 )
2257 shellLVLineTag = gmsh.model.occ.addLine(
2258 shellULpointTag, shellLLpointTag
2259 )
2261 shellLowerLeftLineTag = gmsh.model.occ.addLine(
2262 shellLLpointTag, airLLpointTag
2263 )
2264 shellLowerRightLineTag = gmsh.model.occ.addLine(
2265 shellLRpointTag, airLRpointTag
2266 )
2267 shellUpperLeftLineTag = gmsh.model.occ.addLine(
2268 shellULpointTag, airULpointTag
2269 )
2270 shellUpperRightLineTag = gmsh.model.occ.addLine(
2271 shellURpointTag, airURpointTag
2272 )
2274 # Shell lower surface:
2275 shellLowerPSTag = gmsh.model.occ.addCurveLoop(
2276 [
2277 shellLowerLeftLineTag,
2278 airLHlineTag,
2279 shellLowerRightLineTag,
2280 shellLHlineTag,
2281 ]
2282 )
2283 shellLowerPSTag = gmsh.model.occ.addPlaneSurface([shellLowerPSTag])
2284 self.shellTagsPart1.append(shellLowerPSTag)
2286 # Shell right surface:
2287 shellRightPSTag = gmsh.model.occ.addCurveLoop(
2288 [
2289 shellLowerRightLineTag,
2290 airRVLineTag,
2291 shellUpperRightLineTag,
2292 shellRVLineTag,
2293 ]
2294 )
2295 shellRightPSTag = gmsh.model.occ.addPlaneSurface([shellRightPSTag])
2296 self.shellTagsPart2.append(shellRightPSTag)
2298 # Shell upper surface:
2299 shellUpperPSTag = gmsh.model.occ.addCurveLoop(
2300 [
2301 shellUpperLeftLineTag,
2302 airUHLineTag,
2303 shellUpperRightLineTag,
2304 shellUHLineTag,
2305 ]
2306 )
2307 shellUpperPSTag = gmsh.model.occ.addPlaneSurface([shellUpperPSTag])
2308 self.shellTagsPart3.append(shellUpperPSTag)
2310 # Shell left surface:
2311 shellLeftPSTag = gmsh.model.occ.addCurveLoop(
2312 [
2313 shellLowerLeftLineTag,
2314 airLVLineTag,
2315 shellUpperLeftLineTag,
2316 shellLVLineTag,
2317 ]
2318 )
2319 shellLeftPSTag = gmsh.model.occ.addPlaneSurface([shellLeftPSTag])
2320 self.shellTagsPart4.append(shellLeftPSTag)
2322 def setPrecreatedSurfaceTags(
2323 self,
2324 surfaceTags,
2325 cylinderShellTags=None,
2326 cuboidShellTags1=None,
2327 cuboidShellTags2=None,
2328 cuboidShellTags3=None,
2329 cuboidShellTags4=None,
2330 ):
2331 if not isinstance(surfaceTags, list):
2332 raise TypeError("surfaceTags must be a list.")
2334 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2]
2335 self.surfaceTags.extend(surfaceTags)
2337 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:
2338 # Create innerCircle object from the tags:
2339 curves = gmsh.model.getBoundary(
2340 [(2, tag) for tag in surfaceTags], oriented=False
2341 )
2342 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)
2343 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]
2344 innerCurveTags.sort()
2345 upperRightCurve = innerCurveTags[0]
2346 upperLeftCurve = innerCurveTags[1]
2347 lowerLeftCurve = innerCurveTags[2]
2348 lowerRightCurve = innerCurveTags[3]
2350 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2351 pointTags = [dimTag[1] for dimTag in points]
2352 pointTags.sort()
2353 upperPointTag = pointTags[0]
2354 leftPointTag = pointTags[1]
2356 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2357 pointTags = [dimTag[1] for dimTag in points]
2358 pointTags.sort()
2359 rightPointTag = pointTags[0]
2360 lowerPointTag = pointTags[1]
2362 self.innerCircle = circleWithFourCurves(
2363 0,
2364 0,
2365 self.z,
2366 self.outerRadius,
2367 upperRightTag=upperRightCurve,
2368 upperLeftTag=upperLeftCurve,
2369 lowerLeftTag=lowerLeftCurve,
2370 lowerRightTag=lowerRightCurve,
2371 leftPointTag=leftPointTag,
2372 rightPointTag=rightPointTag,
2373 upperPointTag=upperPointTag,
2374 lowerPointTag=lowerPointTag,
2375 )
2377 if cylinderShellTags is not None:
2378 self.shellTags.extend(cylinderShellTags)
2380 if cuboidShellTags1 is not None:
2381 self.shellTagsPart1.extend(cuboidShellTags1)
2382 self.shellTagsPart2.extend(cuboidShellTags2)
2383 self.shellTagsPart3.extend(cuboidShellTags3)
2384 self.shellTagsPart4.extend(cuboidShellTags4)
2386 def getInnerCL(self):
2387 # checked!
2388 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])
2389 # curves = list(curves)
2390 # curves = list(curves[1])
2392 # innerCL = gmsh.model.occ.addCurveLoop(curves)
2394 gmsh.model.occ.synchronize() # don't delete this line
2395 curves = gmsh.model.getBoundary(
2396 [(2, tag) for tag in self.surfaceTags], oriented=False
2397 )
2398 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)
2399 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]
2401 innerCL = gmsh.model.occ.addCurveLoop(innerCurveTags)
2402 return innerCL
2405class outerTerminalSurface:
2406 def __init__(
2407 self,
2408 outerRadius,
2409 tubeThickness,
2410 divideIntoFourParts=False,
2411 ):
2412 self.tubeSurfaceTags = []
2413 self.nontubeSurfaceTags = []
2415 self.divideIntoFourParts = divideIntoFourParts
2417 self.outerRadius = outerRadius
2418 self.tubeThickness = tubeThickness
2420 def createNontubePartWithMiddleCircleAndWinding(
2421 self, middleCircle: circleWithFourCurves, winding: spiralSurface
2422 ):
2423 leftLineTag = gmsh.model.occ.addLine(
2424 middleCircle.leftPointTag, winding.getOuterLeftPointTag()
2425 )
2426 rightLineTag = gmsh.model.occ.addLine(
2427 middleCircle.rightPointTag, winding.getOuterRightPointTag()
2428 )
2429 upperLineTag = gmsh.model.occ.addLine(
2430 middleCircle.upperPointTag, winding.getOuterUpperPointTag()
2431 )
2432 lowerLineTag = gmsh.model.occ.addLine(
2433 middleCircle.lowerPointTag, winding.getOuterLowerPointTag()
2434 )
2436 # Create surfaces for the nontube part:
2437 if winding.direction is direction.ccw:
2438 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2439 winding.getOuterUpperRightCurves()
2440 # + winding.getOuterEndCurves()
2441 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]
2442 )
2443 else:
2444 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2445 winding.getOuterUpperRightCurves()
2446 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]
2447 )
2448 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2449 self.nontubeSurfaceTags.append(self.upperRightTag)
2451 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2452 winding.getOuterUpperLeftCurves()
2453 + [leftLineTag, middleCircle.upperLeftTag, upperLineTag]
2454 )
2455 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2456 self.nontubeSurfaceTags.append(self.upperLeftTag)
2458 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2459 winding.getOuterLowerLeftCurves()
2460 + [leftLineTag, middleCircle.lowerLeftTag, lowerLineTag]
2461 )
2462 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2463 self.nontubeSurfaceTags.append(self.lowerLeftTag)
2465 if winding.direction is direction.ccw:
2466 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2467 winding.getOuterLowerRightCurves()
2468 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]
2469 )
2470 else:
2471 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2472 winding.getOuterLowerRightCurves()
2473 # + winding.getOuterEndCurves()
2474 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]
2475 )
2476 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2477 self.nontubeSurfaceTags.append(self.lowerRightTag)
2479 def createWithOuterAirAndWinding(
2480 self, outerAir: outerAirSurface, winding: spiralSurface, pancakeIndex
2481 ):
2482 # Tube part:
2483 z = outerAir.z
2485 if self.divideIntoFourParts:
2486 outerCircle = outerAir.innerCircle
2487 middleCircle = circleWithFourCurves(
2488 0, 0, z, self.outerRadius - self.tubeThickness
2489 )
2491 leftLineTag = gmsh.model.occ.addLine(
2492 outerCircle.leftPointTag, middleCircle.leftPointTag
2493 )
2494 rightLineTag = gmsh.model.occ.addLine(
2495 outerCircle.rightPointTag, middleCircle.rightPointTag
2496 )
2497 upperLineTag = gmsh.model.occ.addLine(
2498 outerCircle.upperPointTag, middleCircle.upperPointTag
2499 )
2500 lowerLineTag = gmsh.model.occ.addLine(
2501 outerCircle.lowerPointTag, middleCircle.lowerPointTag
2502 )
2504 # Create surfaces for the tube part:
2505 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2506 [
2507 outerCircle.upperRightTag,
2508 rightLineTag,
2509 middleCircle.upperRightTag,
2510 -upperLineTag,
2511 ]
2512 )
2513 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2514 self.tubeSurfaceTags.append(self.upperRightTag)
2516 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2517 [
2518 outerCircle.upperLeftTag,
2519 leftLineTag,
2520 middleCircle.upperLeftTag,
2521 -upperLineTag,
2522 ]
2523 )
2524 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2525 self.tubeSurfaceTags.append(self.upperLeftTag)
2527 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2528 [
2529 outerCircle.lowerLeftTag,
2530 leftLineTag,
2531 middleCircle.lowerLeftTag,
2532 -lowerLineTag,
2533 ]
2534 )
2535 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2536 self.tubeSurfaceTags.append(self.lowerLeftTag)
2538 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2539 [
2540 outerCircle.lowerRightTag,
2541 rightLineTag,
2542 middleCircle.lowerRightTag,
2543 -lowerLineTag,
2544 ]
2545 )
2546 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2547 self.tubeSurfaceTags.append(self.lowerRightTag)
2549 else:
2550 outerCL = outerAir.getInnerCL()
2552 middleCL = gmsh.model.occ.addCircle(
2553 0, 0, z, self.outerRadius - self.tubeThickness
2554 )
2555 middleCL = gmsh.model.occ.addCurveLoop([middleCL])
2557 tubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
2558 self.tubeSurfaceTags.append(tubeSurface)
2560 # Nontube part:
2561 if self.divideIntoFourParts:
2562 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2564 else:
2565 # Inner and outer curve loops. Sometimes the opposite curve loops are used
2566 # because the other one comes from the self.contactTerminalSurface. To create
2567 # a valid surface, the topology of the curve loops should be consistent. See the
2568 # note in the spiralSurface class.
2569 if pancakeIndex % 2 == 0:
2570 innerCL = winding.outerCurveLoopTag
2571 elif pancakeIndex % 2 == 1:
2572 innerCL = winding.outerOppositeCurveLoopTag
2574 # potential bug (curve order might be wrong)
2575 if self.divideIntoFourParts:
2576 middleCL = gmsh.model.occ.addCurveLoop(
2577 [
2578 middleCircle.upperRightTag,
2579 middleCircle.upperLeftTag,
2580 middleCircle.lowerLeftTag,
2581 middleCircle.lowerRightTag,
2582 ]
2583 )
2585 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2586 self.nontubeSurfaceTags.append(nontubeSurface)
2588 def createWithWindingAndTubeTags(
2589 self, winding: spiralSurface, tubeTags, pancakeIndex
2590 ):
2591 if not isinstance(tubeTags, list):
2592 raise TypeError("tubeTags must be a list.")
2594 self.tubeSurfaceTags.extend(tubeTags)
2596 middleCurves = gmsh.model.getBoundary(
2597 [(2, tag) for tag in tubeTags], oriented=False
2598 )
2599 middleCurveDimTags = findOuterOnes(middleCurves, findInnerOnes=True)
2600 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]
2602 if self.divideIntoFourParts:
2603 # Create middleCircle object from the tags:
2604 middleCurveTags.sort()
2605 upperRightCurve = middleCurveTags[0]
2606 upperLeftCurve = middleCurveTags[1]
2607 lowerLeftCurve = middleCurveTags[2]
2608 lowerRightCurve = middleCurveTags[3]
2610 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2611 pointTags = [dimTag[1] for dimTag in points]
2612 pointTags.sort()
2613 upperPointTag = pointTags[0]
2614 leftPointTag = pointTags[1]
2616 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2617 pointTags = [dimTag[1] for dimTag in points]
2618 pointTags.sort()
2619 rightPointTag = pointTags[0]
2620 lowerPointTag = pointTags[1]
2622 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]
2623 middleCircle = circleWithFourCurves(
2624 0,
2625 0,
2626 z,
2627 self.outerRadius - self.tubeThickness,
2628 upperRightTag=upperRightCurve,
2629 upperLeftTag=upperLeftCurve,
2630 lowerLeftTag=lowerLeftCurve,
2631 lowerRightTag=lowerRightCurve,
2632 leftPointTag=leftPointTag,
2633 rightPointTag=rightPointTag,
2634 upperPointTag=upperPointTag,
2635 lowerPointTag=lowerPointTag,
2636 )
2638 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2639 else:
2640 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)
2642 if pancakeIndex % 2 == 0:
2643 innerCL = winding.outerCurveLoopTag
2644 elif pancakeIndex % 2 == 1:
2645 innerCL = winding.outerOppositeCurveLoopTag
2647 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2648 self.nontubeSurfaceTags.append(nontubeSurface)
2651class innerAirSurface:
2652 def __init__(
2653 self, radius, divideIntoFourParts=False, divideTerminalPartIntoFourParts=False
2654 ):
2655 self.surfaceTags = []
2657 self.divideIntoFourParts = divideIntoFourParts
2658 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts
2660 self.radius = radius
2662 def createFromScratch(self, z):
2663 self.z = z
2664 if self.divideIntoFourParts:
2665 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)
2667 originTag = point(0, 0, z).tag
2669 leftLineTag = gmsh.model.occ.addLine(
2670 self.outerCircle.leftPointTag, originTag
2671 )
2672 rightLineTag = gmsh.model.occ.addLine(
2673 self.outerCircle.rightPointTag, originTag
2674 )
2675 upperLineTag = gmsh.model.occ.addLine(
2676 self.outerCircle.upperPointTag, originTag
2677 )
2678 lowerLineTag = gmsh.model.occ.addLine(
2679 self.outerCircle.lowerPointTag, originTag
2680 )
2682 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2683 [self.outerCircle.upperRightTag, rightLineTag, -upperLineTag]
2684 )
2685 upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2686 self.surfaceTags.append(upperRightTag)
2688 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2689 [self.outerCircle.upperLeftTag, leftLineTag, -upperLineTag]
2690 )
2691 upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2692 self.surfaceTags.append(upperLeftTag)
2694 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2695 [self.outerCircle.lowerLeftTag, leftLineTag, -lowerLineTag]
2696 )
2697 lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2698 self.surfaceTags.append(lowerLeftTag)
2700 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2701 [self.outerCircle.lowerRightTag, rightLineTag, -lowerLineTag]
2702 )
2703 lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2704 self.surfaceTags.append(lowerRightTag)
2706 else:
2707 if self.divideTerminalPartIntoFourParts:
2708 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)
2709 outerCL = gmsh.model.occ.addCurveLoop(
2710 [
2711 self.outerCircle.upperRightTag,
2712 self.outerCircle.upperLeftTag,
2713 self.outerCircle.lowerLeftTag,
2714 self.outerCircle.lowerRightTag,
2715 ]
2716 )
2717 else:
2718 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.radius)
2719 outerCL = gmsh.model.occ.addCurveLoop([outerCL])
2721 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL])
2722 self.surfaceTags.append(surfaceTag)
2724 def setPrecreatedSurfaceTags(self, surfaceTags):
2725 if not isinstance(surfaceTags, list):
2726 raise TypeError("surfaceTags must be a list.")
2728 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2] # potential bug
2729 self.surfaceTags = []
2730 self.surfaceTags.extend(surfaceTags)
2732 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:
2733 # Create outerCirle object from the tags:
2734 curves = gmsh.model.getBoundary(
2735 [(2, tag) for tag in surfaceTags], oriented=False
2736 )
2737 outerCurveDimTags = findOuterOnes(curves)
2738 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]
2739 outerCurveTags.sort()
2740 upperRightCurve = outerCurveTags[0]
2741 upperLeftCurve = outerCurveTags[1]
2742 lowerLeftCurve = outerCurveTags[2]
2743 lowerRightCurve = outerCurveTags[3]
2745 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2746 pointTags = [dimTag[1] for dimTag in points]
2747 pointTags.sort()
2748 upperPointTag = pointTags[0]
2749 leftPointTag = pointTags[1]
2751 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2752 pointTags = [dimTag[1] for dimTag in points]
2753 pointTags.sort()
2754 rightPointTag = pointTags[0]
2755 lowerPointTag = pointTags[1]
2757 self.outerCircle = circleWithFourCurves(
2758 0,
2759 0,
2760 self.z,
2761 self.radius,
2762 upperRightTag=upperRightCurve,
2763 upperLeftTag=upperLeftCurve,
2764 lowerLeftTag=lowerLeftCurve,
2765 lowerRightTag=lowerRightCurve,
2766 leftPointTag=leftPointTag,
2767 rightPointTag=rightPointTag,
2768 upperPointTag=upperPointTag,
2769 lowerPointTag=lowerPointTag,
2770 )
2772 def getOuterCL(self):
2773 # checked!
2774 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])
2775 # curves = list(curves)
2776 # curves = [int(curves[0])]
2778 # outerCL = gmsh.model.occ.addCurveLoop([curves[0]])
2780 # return outerCL
2782 curves = gmsh.model.getBoundary(
2783 [(2, tag) for tag in self.surfaceTags], oriented=False
2784 )
2785 outerCurveDimTags = findOuterOnes(curves)
2786 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]
2788 outerCL = gmsh.model.occ.addCurveLoop(outerCurveTags)
2789 return outerCL
2792class innerTerminalSurface:
2793 def __init__(self, innerRadius, tubeThickness, divideIntoFourParts=False):
2794 self.tubeSurfaceTags = []
2795 self.nontubeSurfaceTags = []
2797 self.divideIntoFourParts = divideIntoFourParts
2799 self.innerRadius = innerRadius
2800 self.tubeThickness = tubeThickness
2802 def createNontubePartWithMiddleCircleAndWinding(
2803 self, middleCircle: circleWithFourCurves, winding: spiralSurface
2804 ):
2805 leftLineTag = gmsh.model.occ.addLine(
2806 winding.getInnerLeftPointTag(), middleCircle.leftPointTag
2807 )
2808 rightLineTag = gmsh.model.occ.addLine(
2809 winding.getInnerRightPointTag(), middleCircle.rightPointTag
2810 )
2811 upperLineTag = gmsh.model.occ.addLine(
2812 winding.getInnerUpperPointTag(), middleCircle.upperPointTag
2813 )
2814 lowerLineTag = gmsh.model.occ.addLine(
2815 winding.getInnerLowerPointTag(), middleCircle.lowerPointTag
2816 )
2818 # Create surfaces for the nontube part:
2819 if winding.direction is direction.ccw:
2820 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2821 winding.getInnerUpperRightCurves()
2822 + [
2823 upperLineTag,
2824 middleCircle.upperRightTag,
2825 rightLineTag,
2826 ]
2827 )
2828 else:
2829 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2830 winding.getInnerUpperRightCurves()
2831 + [
2832 upperLineTag,
2833 middleCircle.upperRightTag,
2834 rightLineTag,
2835 ]
2836 # + winding.getInnerStartCurves()
2837 )
2838 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2839 self.nontubeSurfaceTags.append(self.upperRightTag)
2841 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2842 winding.getInnerUpperLeftCurves()
2843 + [
2844 upperLineTag,
2845 middleCircle.upperLeftTag,
2846 leftLineTag,
2847 ]
2848 )
2849 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2850 self.nontubeSurfaceTags.append(self.upperLeftTag)
2852 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2853 winding.getInnerLowerLeftCurves()
2854 + [
2855 lowerLineTag,
2856 middleCircle.lowerLeftTag,
2857 leftLineTag,
2858 ]
2859 )
2860 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2861 self.nontubeSurfaceTags.append(self.lowerLeftTag)
2863 if winding.direction is direction.ccw:
2864 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2865 winding.getInnerLowerRightCurves()
2866 + [
2867 lowerLineTag,
2868 middleCircle.lowerRightTag,
2869 rightLineTag,
2870 ]
2871 # + winding.getInnerStartCurves()
2872 )
2873 else:
2874 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2875 winding.getInnerLowerRightCurves()
2876 + [
2877 lowerLineTag,
2878 middleCircle.lowerRightTag,
2879 rightLineTag,
2880 ]
2881 )
2882 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2883 self.nontubeSurfaceTags.append(self.lowerRightTag)
2885 def createWithInnerAirAndWinding(
2886 self, innerAir: innerAirSurface, winding: spiralSurface, pancakeIndex
2887 ):
2888 z = innerAir.z
2890 # Tube part:
2891 if self.divideIntoFourParts:
2892 innerCircle = innerAir.outerCircle
2893 middleCircle = circleWithFourCurves(
2894 0, 0, z, self.innerRadius + self.tubeThickness
2895 )
2897 leftLineTag = gmsh.model.occ.addLine(
2898 middleCircle.leftPointTag, innerCircle.leftPointTag
2899 )
2900 rightLineTag = gmsh.model.occ.addLine(
2901 middleCircle.rightPointTag, innerCircle.rightPointTag
2902 )
2903 upperLineTag = gmsh.model.occ.addLine(
2904 middleCircle.upperPointTag, innerCircle.upperPointTag
2905 )
2906 lowerLineTag = gmsh.model.occ.addLine(
2907 middleCircle.lowerPointTag, innerCircle.lowerPointTag
2908 )
2910 # Create surfaces for the tube part:
2911 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2912 [
2913 middleCircle.upperRightTag,
2914 rightLineTag,
2915 innerCircle.upperRightTag,
2916 -upperLineTag,
2917 ]
2918 )
2919 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2920 self.tubeSurfaceTags.append(self.upperRightTag)
2922 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2923 [
2924 middleCircle.upperLeftTag,
2925 leftLineTag,
2926 innerCircle.upperLeftTag,
2927 -upperLineTag,
2928 ]
2929 )
2930 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2931 self.tubeSurfaceTags.append(self.upperLeftTag)
2933 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2934 [
2935 middleCircle.lowerLeftTag,
2936 leftLineTag,
2937 innerCircle.lowerLeftTag,
2938 -lowerLineTag,
2939 ]
2940 )
2941 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2942 self.tubeSurfaceTags.append(self.lowerLeftTag)
2944 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2945 [
2946 middleCircle.lowerRightTag,
2947 rightLineTag,
2948 innerCircle.lowerRightTag,
2949 -lowerLineTag,
2950 ]
2951 )
2952 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2953 self.tubeSurfaceTags.append(self.lowerRightTag)
2955 else:
2956 innerCL = innerAir.getOuterCL()
2958 middleCL = gmsh.model.occ.addCircle(
2959 0, 0, z, self.innerRadius + self.tubeThickness
2960 )
2961 middleCL = gmsh.model.occ.addCurveLoop([middleCL])
2963 tubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2964 self.tubeSurfaceTags.append(tubeSurface)
2966 # Nontube part:
2967 if self.divideIntoFourParts:
2968 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2970 else:
2971 # Inner and outer curve loops. Sometimes the opposite curve loops are used
2972 # because the other one comes from the self.contactTerminalSurface. To create
2973 # a valid surface, the topology of the curve loops should be consistent. See the
2974 # note in the spiralSurface class.
2975 if pancakeIndex == 0:
2976 outerCL = winding.innerCurveLoopTag
2978 elif pancakeIndex % 2 == 0:
2979 outerCL = winding.innerOppositeCurveLoopTag
2981 elif pancakeIndex % 2 == 1:
2982 outerCL = winding.innerCurveLoopTag
2984 # potential bug (curve order might be wrong)
2985 if self.divideIntoFourParts:
2986 middleCL = gmsh.model.occ.addCurveLoop(
2987 [
2988 middleCircle.upperRightTag,
2989 middleCircle.upperLeftTag,
2990 middleCircle.lowerLeftTag,
2991 middleCircle.lowerRightTag,
2992 ]
2993 )
2995 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
2996 self.nontubeSurfaceTags.append(nontubeSurface)
2998 def createWithWindingAndTubeTags(
2999 self, winding: spiralSurface, tubeTags, pancakeIndex
3000 ):
3001 if not isinstance(tubeTags, list):
3002 raise TypeError("tubeTags must be a list.")
3004 self.tubeSurfaceTags.extend(tubeTags)
3006 middleCurves = gmsh.model.getBoundary(
3007 [(2, tag) for tag in tubeTags], oriented=False
3008 )
3009 middleCurveDimTags = findOuterOnes(middleCurves)
3010 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]
3012 if self.divideIntoFourParts:
3013 # Create middleCircle object from the tags:
3014 middleCurveTags.sort()
3015 upperRightCurve = middleCurveTags[0]
3016 upperLeftCurve = middleCurveTags[1]
3017 lowerLeftCurve = middleCurveTags[2]
3018 lowerRightCurve = middleCurveTags[3]
3020 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
3021 pointTags = [dimTag[1] for dimTag in points]
3022 pointTags.sort()
3023 upperPointTag = pointTags[0]
3024 leftPointTag = pointTags[1]
3026 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
3027 pointTags = [dimTag[1] for dimTag in points]
3028 pointTags.sort()
3029 rightPointTag = pointTags[0]
3030 lowerPointTag = pointTags[1]
3032 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]
3033 middleCircle = circleWithFourCurves(
3034 0,
3035 0,
3036 z,
3037 self.innerRadius + self.tubeThickness,
3038 upperRightTag=upperRightCurve,
3039 upperLeftTag=upperLeftCurve,
3040 lowerLeftTag=lowerLeftCurve,
3041 lowerRightTag=lowerRightCurve,
3042 leftPointTag=leftPointTag,
3043 rightPointTag=rightPointTag,
3044 upperPointTag=upperPointTag,
3045 lowerPointTag=lowerPointTag,
3046 )
3048 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
3049 else:
3050 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)
3052 if pancakeIndex == 0:
3053 outerCL = winding.innerCurveLoopTag
3055 elif pancakeIndex % 2 == 0:
3056 outerCL = winding.innerOppositeCurveLoopTag
3058 elif pancakeIndex % 2 == 1:
3059 outerCL = winding.innerCurveLoopTag
3061 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
3062 self.nontubeSurfaceTags.append(nontubeSurface)
3065class pancakeCoilsWithAir:
3066 """
3067 A class to create Pancake3D coil. With this class, any number of pancake coils stack
3068 can be created in GMSH. Moreover, the class also creates some parts of the air
3069 volume as well. It creates the inner cylinder air volume and the outer tube air
3070 volume. However, the air between the pancake coils is not created. It is created in
3071 the gapAir class.
3073 self.fundamentalSurfaces are the surfaces at the bottom of each pancake coil. They
3074 are created one by one with self.generateFundamentalSurfaces() method. The first
3075 created self.fundamentalSurfaces are the first pancake coil's bottom surfaces. Those
3076 surfaces include the outer air shell surface, outer air tube surface, outer
3077 terminal's outer tube part, outer terminal's touching part, winding surfaces,
3078 contact layer surfaces, inner terminal's touching part, inner terminal's inner tube
3079 part, and inner air disc. Terminals are divided into two because they can only be
3080 connected with the perfect tubes since each pancake coil is rotated in a different
3081 direction.
3083 For the first pancake coil, self.fundamentalSurfaces are extruded downwards to
3084 connect the terminal to the end of the geometry at the bottom.
3086 Then self.fundamentalSurfaces are extruded upwards to create the first pancake coil
3087 with self.extrudeWindingPart method. The method returns the extrusion's top
3088 surfaces, which are saved in the topSurfaces variable.
3090 Then those topSurfaces variable is given to another method, self.extrudeGapPart, and
3091 they are further extruded upwards up to the bottom of the next pancake coil.
3092 However, only the air tube, air cylinder, and connection terminal (the perfect inner
3093 terminal tube or outer terminal tube) are extruded. Otherwise, conformality would be
3094 impossible. The gaps are filled with air in gapAir class with fragment operation
3095 later. Then the top surfaces are returned by the method and saved in
3096 self.contactSurfaces variable.
3098 Then using the self.contactSurfaces, self.generateFundamentalSurfaces method creates
3099 the new fundamental surfaces. All the surfaces from the self.contactSurfaces are
3100 used in the new self.fundamentalSurfaces variable to avoid surface duplication.
3102 The logic goes until the last topSurfaces are extruded upwards to connect the last
3103 terminal to the top of the geometry.
3105 Every pancake coil's rotation direction is different each time. Otherwise, their
3106 magnetic fields would neutralize each other.
3108 The first and second pancake coils are connected with the inner terminal. Then the
3109 second and the third pancake coils are connected with the outer terminal. And so on.
3111 :param geometryData: geometry information
3112 """
3114 def __init__(self, geometryData, meshData) -> None:
3115 logger.info("Generating pancake coils has been started.")
3116 start_time = timeit.default_timer()
3118 # Data:
3119 self.geo = geometryData
3120 self.mesh = meshData
3122 # ==============================================================================
3123 # CREATING VOLUME STORAGES STARTS ==============================================
3124 # ==============================================================================
3125 # Air shell (they will be empty if shellTransformation == False):
3126 # For cylinder type:
3127 self.airShellVolume = dimTags(name=self.geo.air.shellVolumeName, save=True)
3129 # For cuboid type:
3130 self.airShellVolumePart1 = dimTags(
3131 name=self.geo.air.shellVolumeName + "-Part1", save=True
3132 )
3133 self.airShellVolumePart2 = dimTags(
3134 name=self.geo.air.shellVolumeName + "-Part2", save=True
3135 )
3136 self.airShellVolumePart3 = dimTags(
3137 name=self.geo.air.shellVolumeName + "-Part3", save=True
3138 )
3139 self.airShellVolumePart4 = dimTags(
3140 name=self.geo.air.shellVolumeName + "-Part4", save=True
3141 )
3143 # Outer air tube volume (actually it is not a tube if the air type is cuboid):
3144 self.outerAirTubeVolume = dimTags(
3145 name=self.geo.air.name + "-OuterTube", save=True, parentName=self.geo.air.name
3146 )
3148 # Outer terminal's outer tube part:
3149 self.outerTerminalTubeVolume = dimTags(
3150 name=self.geo.terminals.outer.name + "-Tube", save=True, parentName=self.geo.terminals.outer.name
3151 )
3153 # Outer terminal's volume that touches the winding:
3154 self.outerTerminalTouchingVolume = dimTags(
3155 name=self.geo.terminals.outer.name + "-Touching",
3156 save=True,
3157 parentName=self.geo.terminals.outer.name,
3158 )
3160 # Inner terminal's volume that touches the winding:
3161 self.innerTerminalTouchingVolume = dimTags(
3162 name=self.geo.terminals.inner.name + "-Touching",
3163 save=True,
3164 parentName=self.geo.terminals.inner.name,
3165 )
3167 # Inner terminal's inner tube part:
3168 self.innerTerminalTubeVolume = dimTags(
3169 name=self.geo.terminals.inner.name + "-Tube", save=True, parentName=self.geo.terminals.inner.name
3170 )
3172 # Transition layers:
3173 self.innerTransitionNotchVolume = dimTags(
3174 name="innerTransitionNotch",
3175 save=True,
3176 )
3177 self.outerTransitionNotchVolume = dimTags(
3178 name="outerTransitionNotch",
3179 save=True,
3180 )
3182 # Inner air cylinder volume:
3183 self.centerAirCylinderVolume = dimTags(
3184 name=self.geo.air.name + "-InnerCylinder",
3185 save=True,
3186 parentName=self.geo.air.name,
3187 )
3189 # Top and bottom parts of the air volume:
3190 self.topAirPancakeWindingExtursionVolume = dimTags(
3191 name=self.geo.air.name + "-TopPancakeWindingExtursion",
3192 save=True,
3193 parentName=self.geo.air.name,
3194 )
3195 self.topAirPancakeContactLayerExtursionVolume = dimTags(
3196 name=self.geo.air.name + "-TopPancakeContactLayerExtursion",
3197 save=True,
3198 parentName=self.geo.air.name,
3199 )
3200 self.topAirTerminalsExtrusionVolume = dimTags(
3201 name=self.geo.air.name + "-TopTerminalsExtrusion",
3202 save=True,
3203 parentName=self.geo.air.name,
3204 )
3205 self.topAirTubeTerminalsExtrusionVolume = dimTags(
3206 name=self.geo.air.name + "-TopTubeTerminalsExtrusion",
3207 save=True,
3208 parentName=self.geo.air.name,
3209 )
3211 self.bottomAirPancakeWindingExtursionVolume = dimTags(
3212 name=self.geo.air.name + "-BottomPancakeWindingExtursion",
3213 save=True,
3214 parentName=self.geo.air.name,
3215 )
3216 self.bottomAirPancakeContactLayerExtursionVolume = dimTags(
3217 name=self.geo.air.name + "-BottomPancakeContactLayerExtursion",
3218 save=True,
3219 parentName=self.geo.air.name,
3220 )
3221 self.bottomAirTerminalsExtrusionVolume = dimTags(
3222 name=self.geo.air.name + "-BottomTerminalsExtrusion",
3223 save=True,
3224 parentName=self.geo.air.name,
3225 )
3226 self.bottomAirTubeTerminalsExtrusionVolume = dimTags(
3227 name=self.geo.air.name + "-BottomTubeTerminalsExtrusion",
3228 save=True,
3229 parentName=self.geo.air.name,
3230 )
3232 # Gap air:
3233 self.gapAirVolume = dimTags(
3234 name=self.geo.air.name + "-Gap", save=True, parentName=self.geo.air.name
3235 )
3237 # Create additional/optional volume storages (they might be used in the meshing
3238 # process):
3239 self.firstTerminalVolume = dimTags(name=self.geo.terminals.firstName, save=True)
3240 self.lastTerminalVolume = dimTags(name=self.geo.terminals.lastName, save=True)
3242 # ==============================================================================
3243 # CREATING VOLUME STORAGES ENDS ================================================
3244 # ==============================================================================
3246 # self.fundamentalSurfaces is a dictionary of surface dimTags tuples. The keys
3247 # are the dimTags objects of the corresponding volumes. The values are the
3248 # dimTags tuples of the surfaces that are used to extrude the volumes. It is
3249 # created in self.generateFundamentalSurfaces method.
3250 self.fundamentalSurfaces = {}
3252 # self.pancakeIndex stores the index of the current pancake coil.
3253 self.pancakeIndex = 0
3255 # self.contactSurfaces is a dictionary of surface dimTags tuples. The keys are
3256 # the dimTags objects of the corresponding volumes. The values are the dimTags
3257 # tuples of the surfaces that are obtained from the previous extrusion and used
3258 # for the next extrusion. The same surface is used for the next extrusion to
3259 # avoid surface duplication. It is created in self.extrudeGapPart and
3260 # self.extrudeWindingPart methods.
3261 self.contactSurfaces = {}
3263 # They will be lists of dimTags objects:
3264 self.individualWinding = []
3265 self.individualContactLayer = []
3267 self.gapAirSurfacesDimTags = []
3269 for i in range(self.geo.numberOfPancakes):
3270 # Itterate over the number of pancake coils:
3271 self.individualWinding.append(
3272 dimTags(
3273 name=self.geo.winding.name + str(self.pancakeIndex + 1),
3274 save=True,
3275 parentName=self.geo.winding.name,
3276 )
3277 )
3278 self.individualContactLayer.append(
3279 dimTags(
3280 name=self.geo.contactLayer.name + str(self.pancakeIndex + 1),
3281 save=True,
3282 parentName=self.geo.contactLayer.name,
3283 )
3284 )
3286 # Generate the fundamental surfaces:
3287 self.fundamentalSurfaces = self.generateFundamentalSurfaces()
3289 # Create gap air or collect the gap air surfaces:
3290 if i != 0:
3291 bottomSurfacesDimTags = []
3292 for key, value in topSurfaces.items():
3293 if (
3294 key is self.individualWinding[self.pancakeIndex - 1]
3295 or key is self.individualContactLayer[self.pancakeIndex - 1]
3296 or key is self.outerTerminalTouchingVolume
3297 or key is self.innerTerminalTouchingVolume
3298 or key is self.innerTransitionNotchVolume
3299 or key is self.outerTransitionNotchVolume
3300 ):
3301 bottomSurfacesDimTags.extend(value)
3303 topSurfacesDimTags = []
3304 for key, value in self.fundamentalSurfaces.items():
3305 if (
3306 key is self.individualWinding[self.pancakeIndex]
3307 or key is self.individualContactLayer[self.pancakeIndex]
3308 or key is self.outerTerminalTouchingVolume
3309 or key is self.innerTerminalTouchingVolume
3310 or key is self.innerTransitionNotchVolume
3311 or key is self.outerTransitionNotchVolume
3312 ):
3313 topSurfacesDimTags.extend(value)
3315 sideSurfacesDimTags = []
3316 if i % 2 == 1:
3317 # Touches it tube and air tube
3318 bottomSurfacesDimTags.extend(
3319 topSurfaces[self.outerTerminalTubeVolume]
3320 )
3321 topSurfacesDimTags.extend(
3322 self.fundamentalSurfaces[self.outerTerminalTubeVolume]
3323 )
3325 if self.mesh.terminals.structured:
3326 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(
3327 3
3328 )[-4:]
3329 else:
3330 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(
3331 3
3332 )[-1:]
3334 lastItTubeSurfsDimTags = gmsh.model.getBoundary(
3335 lastItTubeVolDimTags, oriented=False
3336 )
3337 lastItTubeSideSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3338 lastItTubeSurfsDimTags
3339 )
3340 sideSurfacesDimTags.extend(
3341 findOuterOnes(lastItTubeSideSurfsDimTags)
3342 )
3344 if self.mesh.air.structured:
3345 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[
3346 -4:
3347 ]
3348 else:
3349 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[
3350 -1:
3351 ]
3352 lastAirTubeSurfsDimTags = gmsh.model.getBoundary(
3353 lastAirTubeVolDimTags, oriented=False
3354 )
3355 lastAirTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3356 lastAirTubeSurfsDimTags
3357 )
3358 sideSurfacesDimTags.extend(
3359 findOuterOnes(lastAirTubeSurfsDimTags, findInnerOnes=True)
3360 )
3362 else:
3363 # Touches ot tube and air cylinder
3364 bottomSurfacesDimTags.extend(
3365 topSurfaces[self.innerTerminalTubeVolume]
3366 )
3367 topSurfacesDimTags.extend(
3368 self.fundamentalSurfaces[self.innerTerminalTubeVolume]
3369 )
3370 if self.mesh.terminals.structured:
3371 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(
3372 3
3373 )[-4:]
3374 else:
3375 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(
3376 3
3377 )[-1:]
3379 lastOtTubeSurfsDimTags = gmsh.model.getBoundary(
3380 lastOtTubeVolDimTags, oriented=False
3381 )
3382 lastOtTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3383 lastOtTubeSurfsDimTags
3384 )
3385 sideSurfacesDimTags.extend(
3386 findOuterOnes(lastOtTubeSurfsDimTags, findInnerOnes=True)
3387 )
3389 if self.mesh.air.structured:
3390 lastAirCylinderVolDimTags = (
3391 self.centerAirCylinderVolume.getDimTags(3)[-4:]
3392 )
3393 else:
3394 lastAirCylinderVolDimTags = (
3395 self.centerAirCylinderVolume.getDimTags(3)[-1:]
3396 )
3398 lastAirCylinderSurfsDimTags = gmsh.model.getBoundary(
3399 lastAirCylinderVolDimTags, oriented=False
3400 )
3401 lastAirCylinderSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3402 lastAirCylinderSurfsDimTags
3403 )
3404 sideSurfacesDimTags.extend(
3405 findOuterOnes(lastAirCylinderSurfsDimTags)
3406 )
3408 allGapAirSurfacesDimTags = (
3409 bottomSurfacesDimTags + topSurfacesDimTags + sideSurfacesDimTags
3410 )
3412 # Technically, since all the boundary surfaces of the gap air volumes
3413 # are found here, we should be able to create the gap air volumes with
3414 # addSurfaceLoop and addVolume functions. However, when those are used,
3415 # Geometry.remove_all_duplicates() will indicate some
3416 # duplicates/ill-shaped geometry entities. The indication is
3417 # gmsh.model.occ.remove_all_duplicates() will change the geometry
3418 # (delete some volumes and create new ones), and I have always thought
3419 # that means there are big errors in the geometry and that geometry
3420 # should not be used.
3422 # Alternatively, using these surface tags, the gap air can be created
3423 # with fragment operations as well. Geometry.remove_all_duplicates()
3424 # will tell everything is fine when the fragment operation is used.
3426 # However, I checked manually as well, the way I am using the
3427 # addSurfaceLoop and addVolume should definitely work (because the end
3428 # result is the same with fragments), and I think it is a gmsh/occ
3429 # related problem. In the end, I realized creating the gap air with
3430 # addSurfaceLoop and addVolume won't even affect the mesh, and
3431 # everything seems conformal and nice. Since the fragment operation
3432 # is also very slow, I decided to use addSurfaceLoop and addVolume them.
3433 # However, I keep it as an option so that if the user feels something
3434 # funny about the geometry, the gap air can be created with fragment
3435 # operations as well.
3437 if not self.geo.air.generateGapAirWithFragment:
3438 allGapAirSurfacesTags = [
3439 dimTag[1] for dimTag in allGapAirSurfacesDimTags
3440 ]
3441 surfaceLoop = gmsh.model.occ.addSurfaceLoop(allGapAirSurfacesTags)
3442 volume = gmsh.model.occ.addVolume([surfaceLoop])
3443 self.gapAirVolume.add([(3, volume)])
3445 else:
3446 # Save the surface tags for a fast fragment operation:
3447 self.gapAirSurfacesDimTags.append(allGapAirSurfacesDimTags)
3449 # self.extrudeSurfaces uses self.fundamentalSurfaces for extrusion and adds
3450 # the new volumes to the dimTags objects and returns the dictionary of the
3451 # new top surfaces. The new top surfaces then will be used in extrudeGapPart
3452 # method.
3453 topSurfaces = self.extrudeWindingPart()
3455 if i == 0:
3456 # If it is the first pancake coil, fundemental surfaces are extruded
3457 # downwards to create the bottom air volume and terminal volume.
3458 _ = self.extrudeGapPart(
3459 self.fundamentalSurfaces,
3460 -self.geo.air.axialMargin,
3461 terminalDimTagsObject=self.outerTerminalTubeVolume,
3462 firstTerminal=True,
3463 )
3465 if not i == self.geo.numberOfPancakes - 1:
3466 # If it is not the last pancake coil, extrude the terminal surface to
3467 # create the next contactTerminalSurface and store the new volume in the
3468 # corresponding dimTags object.
3469 self.contactSurfaces = self.extrudeGapPart(topSurfaces)
3471 else:
3472 # If it is the last pancake coil, extrude the terminal surface all the
3473 # way up to the top and store the new volume in the corresponding
3474 # dimTags object.
3475 _ = self.extrudeGapPart(
3476 topSurfaces,
3477 self.geo.air.axialMargin,
3478 lastTerminal=True,
3479 )
3480 self.pancakeIndex = self.pancakeIndex + 1
3482 # Create the gap air volume:
3483 if self.geo.air.generateGapAirWithFragment and self.geo.numberOfPancakes > 1:
3484 self.generateGapAirWithFragment(self.gapAirSurfacesDimTags)
3486 logger.info(
3487 "Generating pancake coils has been finished in"
3488 f" {timeit.default_timer() - start_time:.2f} s."
3489 )
3491 def generateFundamentalSurfaces(self):
3492 """
3493 Generates the inner air, outer air, winding, contact layer, and terminal surfaces
3494 of the current pancake coil and returns them. It finds the z coordinate of the
3495 surfaces and the direction of the pancake coil, depending on the pancake index.
3497 :return: list of dimTags that contains fundamental surfaces
3498 :rtype: list[tuple[int, int]]
3499 """
3500 fundamentalSurfaces = {}
3502 # Select the direction of the spiral:
3503 if self.pancakeIndex % 2 == 0:
3504 spiralDirection = direction.ccw
3505 else:
3506 spiralDirection = direction.cw
3508 # Calculate the z coordinate of the surfaces:
3509 z = (
3510 -self.geo.air.h / 2
3511 + self.geo.air.axialMargin
3512 + self.pancakeIndex * (self.geo.winding.height + self.geo.gapBetweenPancakes)
3513 )
3515 # Create the winding and contact layer surface:
3516 surface = spiralSurface(
3517 self.geo.winding.innerRadius,
3518 self.geo.winding.thickness,
3519 self.geo.contactLayer.thickness,
3520 self.geo.winding.numberOfTurns,
3521 z,
3522 self.geo.winding.theta_i,
3523 self.geo.terminals.transitionNotchAngle,
3524 spiralDirection,
3525 thinShellApproximation=self.geo.contactLayer.thinShellApproximation,
3526 )
3527 #print("this is the surface")
3528 #print(surface)
3529 # Save the surface tags (if TSA, contactLayerSurfaceTags will be empty):
3530 fundamentalSurfaces[self.individualWinding[self.pancakeIndex]] = [
3531 (2, tag) for tag in surface.surfaceTags
3532 ]
3533 fundamentalSurfaces[self.individualContactLayer[self.pancakeIndex]] = [
3534 (2, tag) for tag in surface.contactLayerSurfaceTags
3535 ]
3537 if self.geo.air.type == "cylinder":
3538 outerAirSurf = outerAirSurface(
3539 self.geo.air.radius,
3540 self.geo.terminals.outer.r,
3541 type="cylinder",
3542 divideIntoFourParts=self.mesh.air.structured,
3543 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,
3544 )
3545 elif self.geo.air.type == "cuboid":
3546 outerAirSurf = outerAirSurface(
3547 self.geo.air.sideLength / 2,
3548 self.geo.terminals.outer.r,
3549 type="cuboid",
3550 divideIntoFourParts=self.mesh.air.structured,
3551 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,
3552 )
3554 outerTerminalSurf = outerTerminalSurface(
3555 self.geo.terminals.outer.r,
3556 self.geo.terminals.outer.thickness,
3557 divideIntoFourParts=self.mesh.terminals.structured,
3558 )
3559 innerTerminalSurf = innerTerminalSurface(
3560 self.geo.terminals.inner.r,
3561 self.geo.terminals.inner.thickness,
3562 divideIntoFourParts=self.mesh.terminals.structured,
3563 )
3564 innerAirSurf = innerAirSurface(
3565 self.geo.terminals.inner.r,
3566 divideIntoFourParts=self.mesh.air.structured,
3567 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,
3568 )
3570 if self.contactSurfaces:
3571 # If self.contactSurfaces is not empty, it means that it is not the
3572 # first pancake coil. In that case, contactSurfaces should be used to
3573 # avoid surface duplication.
3575 # Create outer air:
3576 outerAirPSDimTags = self.contactSurfaces[self.outerAirTubeVolume]
3577 outerAirPSTags = [dimTag[1] for dimTag in outerAirPSDimTags]
3578 if self.geo.air.shellTransformation:
3579 if self.geo.air.type == "cuboid":
3580 cuboidShellDimTags1 = self.contactSurfaces[self.airShellVolumePart1]
3581 cuboidShellTags1 = [dimTag[1] for dimTag in cuboidShellDimTags1]
3582 cuboidShellDimTags2 = self.contactSurfaces[self.airShellVolumePart2]
3583 cuboidShellTags2 = [dimTag[1] for dimTag in cuboidShellDimTags2]
3584 cuboidShellDimTags3 = self.contactSurfaces[self.airShellVolumePart3]
3585 cuboidShellTags3 = [dimTag[1] for dimTag in cuboidShellDimTags3]
3586 cuboidShellDimTags4 = self.contactSurfaces[self.airShellVolumePart4]
3587 cuboidShellTags4 = [dimTag[1] for dimTag in cuboidShellDimTags4]
3588 outerAirSurf.setPrecreatedSurfaceTags(
3589 outerAirPSTags,
3590 cuboidShellTags1=cuboidShellTags1,
3591 cuboidShellTags2=cuboidShellTags2,
3592 cuboidShellTags3=cuboidShellTags3,
3593 cuboidShellTags4=cuboidShellTags4,
3594 )
3595 elif self.geo.air.type == "cylinder":
3596 cylinderShellDimTags = self.contactSurfaces[self.airShellVolume]
3597 cylinderShellTags = [dimTag[1] for dimTag in cylinderShellDimTags]
3598 outerAirSurf.setPrecreatedSurfaceTags(
3599 outerAirPSTags,
3600 cylinderShellTags=cylinderShellTags,
3601 )
3602 else:
3603 outerAirSurf.setPrecreatedSurfaceTags(outerAirPSTags)
3605 # Create inner air:
3606 innerAirPSDimTags = self.contactSurfaces[self.centerAirCylinderVolume]
3607 innerAirPSTags = [dimTag[1] for dimTag in innerAirPSDimTags]
3608 innerAirSurf.setPrecreatedSurfaceTags(innerAirPSTags)
3610 if self.pancakeIndex % 2 == 0:
3611 # In this case, we should create all the surfaces for the inner terminal
3612 # but not for outer terminal. Because it is a pancake coil with an even
3613 # index (self.pancakeIndex%2==0) which means that it is connected to the
3614 # previous pancake coil with outer terminal and the outer terminal
3615 # surface is ready (extruded before).
3617 # Create outer terminal:
3618 outerTerminalTubePSDimTags = self.contactSurfaces[
3619 self.outerTerminalTubeVolume
3620 ]
3621 outerTerminalTubePSTags = [
3622 dimTag[1] for dimTag in outerTerminalTubePSDimTags
3623 ]
3624 outerTerminalSurf.createWithWindingAndTubeTags(
3625 surface, outerTerminalTubePSTags, self.pancakeIndex
3626 )
3628 # Create inner terminal:
3629 innerTerminalSurf.createWithInnerAirAndWinding(
3630 innerAirSurf, surface, self.pancakeIndex
3631 )
3633 else:
3634 # In this case, we should create all the surfaces for the outer terminal
3635 # but not for inner terminal. Because it is a pancake coil with an odd
3636 # index (self.pancakeIndex%2==1) which means that it is connected to the
3637 # previous pancake coil with inner terminal and the inner terminal
3638 # surface is ready (extruded before).
3640 # Create outer terminal:
3641 outerTerminalSurf.createWithOuterAirAndWinding(
3642 outerAirSurf, surface, self.pancakeIndex
3643 )
3645 # Create inner terminal:
3646 innerTerminalTubePSDimTags = self.contactSurfaces[
3647 self.innerTerminalTubeVolume
3648 ]
3649 innerTerminalTubePSTags = [
3650 dimTag[1] for dimTag in innerTerminalTubePSDimTags
3651 ]
3652 innerTerminalSurf.createWithWindingAndTubeTags(
3653 surface, innerTerminalTubePSTags, self.pancakeIndex
3654 )
3656 else:
3657 # If self.contactSurfaces is empty, it means that it is the first pancake
3658 # coil. In that case, the surfaces should be created from scratch.
3660 if self.geo.air.shellTransformation:
3661 if self.geo.air.type == "cuboid":
3662 outerAirSurf.createFromScratch(
3663 z,
3664 shellTransformation=True,
3665 shellRadius=self.geo.air.shellSideLength / 2,
3666 )
3667 else:
3668 outerAirSurf.createFromScratch(
3669 z,
3670 shellTransformation=True,
3671 shellRadius=self.geo.air.shellOuterRadius,
3672 )
3673 else:
3674 outerAirSurf.createFromScratch(z)
3676 innerAirSurf.createFromScratch(z)
3677 outerTerminalSurf.createWithOuterAirAndWinding(
3678 outerAirSurf, surface, self.pancakeIndex
3679 )
3680 innerTerminalSurf.createWithInnerAirAndWinding(
3681 innerAirSurf, surface, self.pancakeIndex
3682 )
3684 # Save the surface tags:
3685 fundamentalSurfaces[self.outerAirTubeVolume] = [
3686 (2, tag) for tag in outerAirSurf.surfaceTags
3687 ]
3689 fundamentalSurfaces[self.centerAirCylinderVolume] = [
3690 (2, tag) for tag in innerAirSurf.surfaceTags
3691 ]
3693 fundamentalSurfaces[self.outerTerminalTubeVolume] = [
3694 (2, tag) for tag in outerTerminalSurf.tubeSurfaceTags
3695 ]
3696 fundamentalSurfaces[self.outerTerminalTouchingVolume] = [
3697 (2, tag) for tag in outerTerminalSurf.nontubeSurfaceTags
3698 ]
3700 fundamentalSurfaces[self.innerTerminalTubeVolume] = [
3701 (2, tag) for tag in innerTerminalSurf.tubeSurfaceTags
3702 ]
3703 fundamentalSurfaces[self.innerTerminalTouchingVolume] = [
3704 (2, tag) for tag in innerTerminalSurf.nontubeSurfaceTags
3705 ]
3706 fundamentalSurfaces[self.innerTransitionNotchVolume] = [
3707 (2, tag) for tag in surface.innerNotchSurfaceTags
3708 ]
3709 fundamentalSurfaces[self.outerTransitionNotchVolume] = [
3710 (2, tag) for tag in surface.outerNotchSurfaceTags
3711 ]
3713 if self.geo.air.shellTransformation:
3714 if self.geo.air.type == "cuboid":
3715 fundamentalSurfaces[self.airShellVolumePart1] = [
3716 (2, tag) for tag in outerAirSurf.shellTagsPart1
3717 ]
3718 fundamentalSurfaces[self.airShellVolumePart2] = [
3719 (2, tag) for tag in outerAirSurf.shellTagsPart2
3720 ]
3721 fundamentalSurfaces[self.airShellVolumePart3] = [
3722 (2, tag) for tag in outerAirSurf.shellTagsPart3
3723 ]
3724 fundamentalSurfaces[self.airShellVolumePart4] = [
3725 (2, tag) for tag in outerAirSurf.shellTagsPart4
3726 ]
3727 elif self.geo.air.type == "cylinder":
3728 fundamentalSurfaces[self.airShellVolume] = [
3729 (2, tag) for tag in outerAirSurf.shellTags
3730 ]
3731 # windingSurfaces = {}
3732 # # Save only the winding surface tags:
3733 # windingSurfaces[self.individualWinding[self.pancakeIndex]] = [
3734 # (2, tag) for tag in surface.surfaceTags
3735 # ]
3736 # fundamentalSurfaces = windingSurfaces
3737 # print(fundamentalSurfaces)
3739 return fundamentalSurfaces
3741 def extrudeGapPart(
3742 self,
3743 surfacesDict,
3744 tZ: float = None,
3745 terminalDimTagsObject: dimTags = None,
3746 firstTerminal=False,
3747 lastTerminal=False,
3748 ):
3749 """
3750 Extrudes the given surfaces dimTags dictionary to a given height (tZ) and adds
3751 the created volumes to the corresponding dictionary keys (dimTags objects). It
3752 returns the extrusion's top surfaces as a dictionary again, where the keys are
3753 the corresponding dimTagsObjects and the values are the dimTags of the surfaces.
3755 If tZ is not given, then it is set to the gap height (self.geo.gapBetweenPancakes). This is the
3756 default value used for connecting the pancake coils. Only for the creation of
3757 the first and the last pancake coils different tZ values are used.
3759 If terminalDimTagsObject is not given, then the created volume is added
3760 automatically to the innerTerminalVolume or outerTerminalVolume dimTags object,
3761 depending on the value of self.pancakeIndex. However, giving
3762 terminalDimTagsObject is necessary for creating the first and the last terminal.
3763 Otherwise, finding out the correct terminal dimTagsObject would be very
3764 challenging.
3766 :param surfaces: the surface dimTag dictionary to be extruded. The keys are the
3767 dimTags objects and the values are the dimTags of the surfaces. The
3768 keys are used to easily add the corresponding volumes to the correct
3769 dimTags objects
3770 :type surfaces: dict[dimTags, list[tuple[int, int]]]
3771 :param tZ: the height of the extrusion
3772 :type tZ: float, optional
3773 :param terminalDimTagsObject: the dimTags object of the terminal to be extruded
3774 :type terminalDimTagsObject: dimTags, optional
3775 :return: top surfaces of the extrusion as a dictionary where the keys are the
3776 dimTags objects and the values are the dimTags of the surfaces
3777 :rtype: dict[dimTags, list[tuple[int, int]]]
3778 """
3779 bottomPart = False
3780 topPart = False
3781 if tZ is None:
3782 tZ = self.geo.gapBetweenPancakes
3783 elif tZ < 0:
3784 bottomPart = True
3785 elif tZ > 0:
3786 topPart = True
3788 if terminalDimTagsObject is None:
3789 # terminalDimTagsObject needs to be given for the first terminal that is
3790 # extruded downwards.
3791 if self.pancakeIndex % 2 == 0:
3792 terminalDimTagsObject = self.innerTerminalTubeVolume
3793 else:
3794 terminalDimTagsObject = self.outerTerminalTubeVolume
3796 # if terminalDimTagsObject is self.innerTerminalVolume:
3797 # otherTerminal = self.outerTerminalVolume
3798 # else:
3799 # otherTerminal = self.innerTerminalVolume
3801 # Create the list of surfaces to be extruded:
3802 listOfDimTags = []
3803 listOfDimTagsObjects = []
3804 listOfDimTagsForTopSurfaces = []
3805 if topPart:
3806 # Then in this case, most of the surfaces should be added to the air volumes
3807 # instead of the terminal, winding, and contact layer volumes.
3808 for key, dimTagsList in surfacesDict.items():
3809 if key is self.individualWinding[self.pancakeIndex]:
3810 dimTagsObjects = [self.topAirPancakeWindingExtursionVolume] * len(
3811 dimTagsList
3812 )
3813 elif key is self.individualContactLayer[self.pancakeIndex]:
3814 dimTagsObjects = [
3815 self.topAirPancakeContactLayerExtursionVolume
3816 ] * len(dimTagsList)
3817 elif (
3818 key is terminalDimTagsObject
3819 or key is self.airShellVolume
3820 or key is self.airShellVolumePart1
3821 or key is self.airShellVolumePart2
3822 or key is self.airShellVolumePart3
3823 or key is self.airShellVolumePart4
3824 or key is self.outerAirTubeVolume
3825 or key is self.centerAirCylinderVolume
3826 ):
3827 dimTagsObjects = [key] * len(dimTagsList)
3828 else:
3829 # key is self.outerTerminalTouchingVolume
3830 # or key is self.innerTerminalTouchingVolume
3831 # or key is (other terminal's tube volume)
3832 dimTagsObjects = [self.topAirTerminalsExtrusionVolume] * len(
3833 dimTagsList
3834 )
3835 if (
3836 key is self.innerTerminalTubeVolume
3837 or key is self.outerTerminalTubeVolume
3838 ):
3839 dimTagsObjects = [
3840 self.topAirTubeTerminalsExtrusionVolume
3841 ] * len(dimTagsList)
3843 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(
3844 dimTagsList
3845 )
3846 listOfDimTags = listOfDimTags + dimTagsList
3847 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3848 elif bottomPart:
3849 # Then in this case, most of the surfaces should be added to the air volumes
3850 # instead of the terminal, winding, and contact layer volumes.
3851 for key, dimTagsList in surfacesDict.items():
3852 if key is self.individualWinding[self.pancakeIndex]:
3853 dimTagsObjects = [
3854 self.bottomAirPancakeWindingExtursionVolume
3855 ] * len(dimTagsList)
3856 elif key is self.individualContactLayer[self.pancakeIndex]:
3857 dimTagsObjects = [
3858 self.bottomAirPancakeContactLayerExtursionVolume
3859 ] * len(dimTagsList)
3860 elif (
3861 key is terminalDimTagsObject
3862 or key is self.airShellVolume
3863 or key is self.airShellVolumePart1
3864 or key is self.airShellVolumePart2
3865 or key is self.airShellVolumePart3
3866 or key is self.airShellVolumePart4
3867 or key is self.outerAirTubeVolume
3868 or key is self.centerAirCylinderVolume
3869 ):
3870 dimTagsObjects = [key] * len(dimTagsList)
3871 else:
3872 # key is self.outerTerminalTouchingVolume
3873 # or key is self.innerTerminalTouchingVolume
3874 # or key is (other terminal's tube volume)
3875 dimTagsObjects = [self.bottomAirTerminalsExtrusionVolume] * len(
3876 dimTagsList
3877 )
3878 if (
3879 key is self.innerTerminalTubeVolume
3880 or key is self.outerTerminalTubeVolume
3881 ):
3882 dimTagsObjects = [
3883 self.bottomAirTubeTerminalsExtrusionVolume
3884 ] * len(dimTagsList)
3886 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(
3887 dimTagsList
3888 )
3889 listOfDimTags = listOfDimTags + dimTagsList
3890 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3891 else:
3892 for key, dimTagsList in surfacesDict.items():
3893 if (
3894 key is self.outerAirTubeVolume
3895 or key is self.centerAirCylinderVolume
3896 or key is self.airShellVolume
3897 or key is self.airShellVolumePart1
3898 or key is self.airShellVolumePart2
3899 or key is self.airShellVolumePart3
3900 or key is self.airShellVolumePart4
3901 or key is terminalDimTagsObject
3902 ):
3903 dimTagsObjects = [key] * len(dimTagsList)
3905 listOfDimTags = listOfDimTags + dimTagsList
3906 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3908 listOfDimTagsForTopSurfaces = listOfDimTagsObjects
3910 extrusionResult = dimTags()
3911 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, tZ))
3913 # Add the created volumes to the corresponding dimTags objects:
3914 volumeDimTags = extrusionResult.getDimTags(3)
3915 for i, volumeDimTag in enumerate(volumeDimTags):
3916 listOfDimTagsObjects[i].add(volumeDimTag)
3918 if firstTerminal:
3919 self.firstTerminalVolume.add(terminalDimTagsObject.getDimTags(3))
3920 elif lastTerminal:
3921 self.lastTerminalVolume.add(terminalDimTagsObject.getDimTags(3))
3923 topSurfacesDimTags = extrusionResult.getExtrusionTop()
3924 topSurfacesDict = {}
3925 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):
3926 if listOfDimTagsObjects[i] in topSurfacesDict:
3927 topSurfacesDict[listOfDimTagsForTopSurfaces[i]].append(topSurfaceDimTag)
3928 else:
3929 topSurfacesDict[listOfDimTagsForTopSurfaces[i]] = [topSurfaceDimTag]
3931 return topSurfacesDict
3933 def extrudeWindingPart(self):
3934 """
3935 Extrudes all the fundamental surfaces of the pancake coil by self.geo.winding.height and
3936 returns the next connection terminal's top surface dimTag, and other air dimTags
3937 in a dictionary so that they can be further extruded.
3939 :return: dictionary of top surfaces where the keys are the dimTags objects and
3940 the values are the dimTags of the surfaces
3941 :rtype: dict[dimTags, list[tuple[int, int]]]
3942 """
3943 # Create the list of surfaces to be extruded:
3944 listOfDimTags = []
3945 listOfDimTagsObjects = []
3946 for key, dimTagsList in self.fundamentalSurfaces.items():
3947 dimTagsObjects = [key] * len(dimTagsList)
3949 listOfDimTags = listOfDimTags + dimTagsList
3950 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3952 # Extrude the fundamental surfaces:
3953 extrusionResult = dimTags()
3954 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, self.geo.winding.height))
3956 # Add the created volumes to the corresponding dimTags objects:
3957 volumes = extrusionResult.getDimTags(3)
3958 for i, volumeDimTag in enumerate(volumes):
3959 listOfDimTagsObjects[i].add(volumeDimTag)
3961 if self.pancakeIndex == 0:
3962 # Note the first pancake (sometimes useful for creating regions in the
3963 # meshing part):
3964 for i, volumeDimTag in enumerate(volumes):
3965 if listOfDimTagsObjects[i].parentName == self.geo.terminals.outer.name:
3966 self.firstTerminalVolume.add(volumeDimTag)
3968 # Not elif! Because the first pancake coil is also the last pancake coil if
3969 # there is only one pancake coil.
3970 if self.pancakeIndex == self.geo.numberOfPancakes - 1:
3971 # Note the last pancake (sometimes useful for creating regions in the
3972 # meshing part):
3973 for i, volumeDimTag in enumerate(volumes):
3974 if (
3975 self.pancakeIndex % 2 == 1
3976 and listOfDimTagsObjects[i].parentName == self.geo.terminals.outer.name
3977 ):
3978 self.lastTerminalVolume.add(volumeDimTag)
3979 elif (
3980 self.pancakeIndex % 2 == 0
3981 and listOfDimTagsObjects[i].parentName == self.geo.terminals.inner.name
3982 ):
3983 self.lastTerminalVolume.add(volumeDimTag)
3985 # Return the top surfaces:
3986 # Add the created top surfaces to a new dictionary:
3987 topSurfacesDimTags = extrusionResult.getExtrusionTop()
3988 topSurfaces = {}
3989 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):
3990 if listOfDimTagsObjects[i] in topSurfaces:
3991 topSurfaces[listOfDimTagsObjects[i]].append(topSurfaceDimTag)
3992 else:
3993 topSurfaces[listOfDimTagsObjects[i]] = [topSurfaceDimTag]
3995 return topSurfaces
3997 def generateGapAirWithFragment(
3998 self, gapAirSurfacesDimTags: List[List[Tuple[int, int]]]
3999 ):
4000 """
4001 A class to fill the gap between the multiple pancake coils with air. First, it
4002 creates a dummy cylinder with the same radius as the outer terminal's outer
4003 radius. Then using gapAirSurfacesDimTags, gmsh.model.occ.fragment() operation is
4004 applied to the dummy cylinder volume in a for loop to create the gap air
4005 volumes. After each fragment operation, one of the volumes created is removed
4006 because it is the solid volume which is the combination of windings,
4007 contact layers, and terminals. In the end, dummy cylinder is removed as well.
4010 WARNING:
4011 Currently, this method doesn't work.
4013 :param geometry: geometry information
4014 :param pancakeCoils: pancakeCoilsWithAir object
4015 :type pancakeCoils: pancakeCoilsWithAir
4016 """
4017 logger.info("Generating gap air has been started.")
4018 start_time = timeit.default_timer()
4020 # Create the dummy air volume:
4021 dummyAir = gmsh.model.occ.addCylinder(
4022 0,
4023 0,
4024 -self.geo.air.h / 2,
4025 0,
4026 0,
4027 self.geo.air.h,
4028 self.geo.terminals.outer.r,
4029 )
4031 toBeDeletedDimTags = []
4032 gapAirVolumesCurrentDimTags = []
4033 for i in range(len(gapAirSurfacesDimTags)):
4034 # Get the outer surfaces of the pancake coils for cutting the pancake coils
4035 # from the dummy air. The outer surfaces are used instead of pancake volumes
4036 # to reduce the amount of work for gmsh. It makes it significantly faster.
4037 # if len(gapAirSurfacesDimTags[i]) !=12:
4038 fragmentResults = gmsh.model.occ.fragment(
4039 [(3, dummyAir)],
4040 gapAirSurfacesDimTags[i],
4041 removeObject=False,
4042 removeTool=False,
4043 )
4044 fragmentVolumeResultsDimTags = fragmentResults[1][0]
4045 toBeDeletedDimTags.append(fragmentVolumeResultsDimTags[0])
4046 gapAirVolumesCurrentDimTags.append(fragmentVolumeResultsDimTags[1])
4048 toBeDeletedDimTags.append((3, dummyAir))
4049 # Fragmnet operation both creates the air volume and solid pancake coils volume
4050 # because the surfaces are used for cutting. Therefore, the solid pancake coils
4051 # volume should be removed from the fragment results:
4052 gmsh.model.occ.remove(toBeDeletedDimTags)
4054 # Add results to the air volume storage. After the geometry is saves as a .brep
4055 # file, and loaded back, the gaps between the tags are avoided by moving the
4056 # the other tags. Therefore, this is how the tags are stored:
4057 toBeDeletedTags = [dimTag[1] for dimTag in toBeDeletedDimTags]
4058 volumeTagsStart = min(toBeDeletedTags)
4059 numberOfGapAirVolumes = len(gapAirVolumesCurrentDimTags)
4060 gapAirVolumesToBeSaved = [
4061 (3, volumeTagsStart + i) for i in range(numberOfGapAirVolumes)
4062 ]
4064 # For debugging purposes, physical groups are being created in the geometry
4065 # generation process as well. Normally, it us done during meshing because BREP
4066 # files cannot store physical groups. Since the tags above (airVolumes) will be
4067 # valid after only saving the geometry as a BREP file and loading it back, the
4068 # current tags are given to the airVolume.add() method as well. This is done to
4069 # be able to create the correct physical group.
4070 self.gapAirVolume.add(
4071 dimTagsList=gapAirVolumesToBeSaved,
4072 dimTagsListForPG=gapAirVolumesCurrentDimTags,
4073 )
4075 logger.info(
4076 "Generating gap air has been finished in"
4077 f" {timeit.default_timer() - start_time:.2f} s."
4078 )
4081class Geometry(Base):
4082 """
4083 Main geometry class for Pancake3D.
4085 :param fdm: FiQuS data model
4086 :param geom_folder: folder where the geometry files are saved
4087 :type geom_folder: str
4088 :param mesh_folder: folder where the mesh files are saved
4089 :type mesh_folder: str
4090 :param solution_folder: folder where the solution files are saved
4091 :type solution_folder: str
4092 """
4094 def __init__(
4095 self,
4096 fdm,
4097 geom_folder,
4098 mesh_folder,
4099 solution_folder,
4100 ) -> None:
4101 super().__init__(fdm, geom_folder, mesh_folder, solution_folder)
4103 # Clear if there is any existing dimTags storage:
4104 dimTagsStorage.clear()
4106 # Start GMSH:
4107 self.gu = GmshUtils(self.geom_folder)
4108 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh)
4110 # To speed up the GUI:
4111 gmsh.option.setNumber("Geometry.NumSubEdges", 10)
4113 # To see the surfaces in a better way in GUI:
4114 gmsh.option.setNumber("Geometry.SurfaceType", 1)
4116 # To avoid any unwanted modifications to the geometry, the automatic fixing of
4117 # the geometry is disabled:
4118 gmsh.option.setNumber("Geometry.OCCAutoFix", 0)
4120 # Set the tolerance:
4121 if self.geo.dimensionTolerance < gmsh.option.getNumber("Geometry.Tolerance"):
4122 gmsh.option.setNumber("Geometry.Tolerance", self.geo.dimensionTolerance)
4124 gmsh.option.setNumber("Geometry.ToleranceBoolean", self.geo.dimensionTolerance)
4126 spiralCurve.sectionsPerTurn = self.geo.winding.spt
4127 spiralCurve.curvesPerTurn = self.geo.winding.numberOfVolumesPerTurn
4129 def generate_geometry(self):
4130 """
4131 Generates geometry and saves it as a .brep file.
4134 """
4135 logger.info(
4136 f"Generating Pancake3D geometry ({self.brep_file}) has been started."
4137 )
4138 start_time = timeit.default_timer()
4140 self.pancakeCoil = pancakeCoilsWithAir(self.geo, self.mesh)
4142 gmsh.model.occ.synchronize()
4143 gmsh.write(self.brep_file)
4145 logger.info(
4146 f"Generating Pancake3D geometry ({self.brep_file}) has been finished in"
4147 f" {timeit.default_timer() - start_time:.2f} s."
4148 )
4150 def load_geometry(self):
4151 """
4152 Loads geometry from .brep file.
4153 """
4154 logger.info("Loading Pancake3D geometry has been started.")
4155 start_time = timeit.default_timer()
4157 previousGeo = FilesAndFolders.read_data_from_yaml(
4158 self.geometry_data_file, Pancake3DGeometry
4159 )
4161 if previousGeo.dict() != self.geo.dict():
4162 raise ValueError(
4163 "Geometry data has been changed. Please regenerate the geometry or load"
4164 " the previous geometry data."
4165 )
4167 gmsh.clear()
4168 gmsh.model.occ.importShapes(self.brep_file, format="brep")
4169 gmsh.model.occ.synchronize()
4171 logger.info(
4172 "Loading Pancake3D geometry has been finished in"
4173 f" {timeit.default_timer() - start_time:.2f} s."
4174 )
4176 def generate_vi_file(self):
4177 """
4178 Generates volume information file. Volume information file stores dimTags of all
4179 the stored volumes in the geometry. Without this file, regions couldn't be
4180 created, meaning that finite element simulation cannot be done.
4182 The file extension is custom because users are not supposed to edit or open this
4183 file, and it makes it intuitively clear that it is a volume information file.
4184 """
4185 logger.info(
4186 f"Generating volume information file ({self.vi_file}) has been started."
4187 )
4188 start_time = timeit.default_timer()
4190 dimTagsDict = dimTagsStorage.getDimTagsDict()
4191 json.dump(
4192 dimTagsDict,
4193 open(self.vi_file, "w"),
4194 )
4196 logger.info(
4197 "Generating volume information file has been finished in"
4198 f" {timeit.default_timer() - start_time:.2f} s."
4199 )
4201 if self.geo_gui:
4202 self.generate_physical_groups()
4203 self.gu.launch_interactive_GUI()
4204 else:
4205 gmsh.clear()
4206 gmsh.finalize()
4208 @staticmethod
4209 def generate_physical_groups():
4210 """
4211 Generates physical groups. Physical groups are not saved in the BREP file but
4212 it can be useful for debugging purposes.
4215 """
4216 gmsh.model.occ.synchronize()
4218 dimTagsDict = dimTagsStorage.getDimTagsDict(forPhysicalGroups=True)
4220 for key, value in dimTagsDict.items():
4221 tags = [dimTag[1] for dimTag in value]
4222 gmsh.model.addPhysicalGroup(
4223 3,
4224 tags,
4225 name=key,
4226 )
4228 @staticmethod
4229 def remove_all_duplicates():
4230 """
4231 Removes all the duplicates and then prints the entities that are created or
4232 removed during the operation. It prints the line number where the function is
4233 called as well. This function is helpful for debugging. Finding duplicates means
4234 there is a problem in geometry creation logic, and the meshes will not be
4235 conformal. It shouldn't be used in the final version of the code since removing
4236 duplicates is computationally expensive, and there shouldn't be duplicates at
4237 all.
4239 WARNING:
4240 This function currently does not work properly. It is not recommended to use
4241 right now. It finds duplicates even if there are no duplicates (topology
4242 problems).
4243 """
4245 logger.info(f"Removing all the duplicates has been started.")
4246 start_time = timeit.default_timer()
4248 gmsh.model.occ.synchronize()
4249 oldEntities = []
4250 oldEntities.extend(gmsh.model.getEntities(3))
4251 oldEntities.extend(gmsh.model.getEntities(2))
4252 oldEntities.extend(gmsh.model.getEntities(1))
4253 oldEntities.extend(gmsh.model.getEntities(0))
4254 oldEntities = set(oldEntities)
4256 gmsh.model.occ.removeAllDuplicates()
4258 gmsh.model.occ.synchronize()
4259 newEntities = []
4260 newEntities.extend(gmsh.model.getEntities(3))
4261 newEntities.extend(gmsh.model.getEntities(2))
4262 newEntities.extend(gmsh.model.getEntities(1))
4263 newEntities.extend(gmsh.model.getEntities(0))
4264 newEntities = set(newEntities)
4265 NewlyCreated = newEntities - oldEntities
4266 Removed = oldEntities - newEntities
4268 frameinfo = getframeinfo(currentframe().f_back)
4270 if len(NewlyCreated) > 0 or len(Removed) > 0:
4271 logger.warning(f"Duplicates found! Line: {frameinfo.lineno}")
4272 logger.warning(f"{len(NewlyCreated)}NewlyCreated = {list(NewlyCreated)}")
4273 logger.warning(f"{len(Removed)}Removed = {list(Removed)}")
4274 else:
4275 logger.info(f"No duplicates found! Line: {frameinfo.lineno}")
4277 logger.info(
4278 "Removing all the duplicates has been finished in"
4279 f" {timeit.default_timer() - start_time:.2f} s."
4280 )