Coverage for fiqus/geom_generators/GeometryPancake3D.py: 81%
1435 statements
« prev ^ index » next coverage.py v6.4.4, created at 2024-05-20 03:24 +0200
1import math
2import logging
3from enum import Enum
4from inspect import currentframe, getframeinfo
5from typing import List, Tuple, Dict
6import operator
8import json
9import timeit
10import numpy as np
11import gmsh
13from fiqus.utils.Utils import GmshUtils
14from fiqus.utils.Utils import FilesAndFolders
15from fiqus.mains.MainPancake3D import Base
16from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry
18logger = logging.getLogger(__name__)
21def findSurfacesWithNormalsOnXYPlane(dimTags):
22 result = []
23 for dimTag in dimTags:
24 surfaceNormal = gmsh.model.getNormal(dimTag[1], [0.5, 0.5])
25 if abs(surfaceNormal[2]) < 1e-6:
26 result.append(dimTag)
28 return result
31def findOuterOnes(dimTags, findInnerOnes=False):
32 """
33 Finds the outermost surface/curve/point in a list of dimTags. The outermost means
34 the furthest from the origin.
35 """
36 dim = dimTags[0][0]
37 if dim == 2:
38 distances = []
39 for dimTag in dimTags:
40 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])
41 for curve in curves:
42 curve = list(curve)
43 pointTags = gmsh.model.getBoundary(
44 [(1, curveTag) for curveTag in curve],
45 oriented=False,
46 combined=False,
47 )
48 # Get the positions of the points:
49 points = []
50 for dimTag in pointTags:
51 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]
52 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]
53 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))
54 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))
56 distances.append(
57 max([point[0] ** 2 + point[1] ** 2 for point in points])
58 )
59 elif dim == 1:
60 distances = []
61 for dimTag in dimTags:
62 pointTags = gmsh.model.getBoundary(
63 [dimTag],
64 oriented=False,
65 combined=False,
66 )
67 # Get the positions of the points:
68 points = []
69 for dimTag in pointTags:
70 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]
71 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]
72 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))
73 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))
75 distances.append(max([point[0] ** 2 + point[1] ** 2 for point in points]))
77 if findInnerOnes:
78 goalDistance = min(distances)
79 else:
80 goalDistance = max(distances)
82 result = []
83 for distance, dimTag in zip(distances, dimTags):
84 # Return all the dimTags with the hoal distance:
85 if math.isclose(distance, goalDistance, abs_tol=1e-6):
86 result.append(dimTag)
88 return result
91class dimTags:
92 """
93 This is a class for (dim, tag) tuple lists. DimTags are heavily used in GMSH, and
94 dimTags class makes it easier to store and manipulate them.
96 Every dimTags instance with the save option True will be stored in the
97 dimTagsStorage class. dimTags instances with the save option False will not be
98 stored. dimTagsStorage class stores all the dimTags with the corresponding names so
99 that later the dimTagsStorage class can be used to create the volume information
100 file. The volume information file is required in the meshing stage to create the
101 appropriate regions.
103 If parentName is specified, during the volume information file generation, another
104 key that is equal to parentName is created, and all the dimTags are also added
105 there. For example, air consists of many parts such as gap, outer tube, inner
106 cylinder parts, and their dimTags object all have the parentName of self.geo.ai.name
107 (user input of the air name) so that in the volume information file, they are all
108 combined under the air name as well.
110 :param name: name of the dimTags object (default: None)
111 :type name: str, optional
112 :param parentName: name of the parent dimTags object (default: None)
113 :type parentName: str, optional
114 :param save: True if the instance to be stored in dimTagsStorage class, False
115 otherwise (default: False)
116 :type save: bool, optional
117 """
119 point = 0
120 curve = 1
121 surface = 2
122 volume = 3
124 def storageUpdateRequired(func):
125 """
126 A decorator for dimTags class. It will update the dimTagsStorage class if the
127 save option is True and the decorated function is called. Use this decorator
128 for every function that changes the dimTags instance so that the storage gets
129 updated.
131 :param func: function to be decorated
132 :type func: function
133 :return: decorated function
134 :rtype: function
135 """
137 def wrapper(self, *args, **kwargs):
138 func(self, *args, **kwargs)
140 if self.save:
141 # Update the dimTagsStorage:
142 dimTagsStorage.updateDimTags(self)
144 return wrapper
146 def __init__(
147 self,
148 name: str = None,
149 parentName: str = None,
150 save: bool = False,
151 ):
152 self.name = name
154 dimTagsObjects = dimTagsStorage.getDimTagsObject(name)
155 if dimTagsObjects != []:
156 dimTagsObject = dimTagsObjects[0]
157 self.physicalTag = dimTagsObject.physicalTag
158 # To store points, curves, surfaces, and volumes separately:
159 self.dimTags = dimTagsObject.dimTags
160 self.dimTagsForPG = dimTagsObject.dimTagsForPG
161 self.allDimTags = dimTagsObject.allDimTags
162 self.parentName = dimTagsObject.parentName
163 self.save = dimTagsObject.save
164 else:
165 self.physicalTag = None
166 # To store points, curves, surfaces, and volumes separately:
167 self.dimTags = [[] for _ in range(4)]
168 self.dimTagsForPG = [[] for _ in range(4)] # dim tags for physical groups
169 self.allDimTags = []
170 self.parentName = parentName
171 self.save = save
172 if self.save:
173 dimTagsStorage.updateDimTags(self)
175 @storageUpdateRequired
176 def add(
177 self,
178 dimTagsList: List[Tuple[int, int]],
179 dimTagsListForPG: List[Tuple[int, int]] = None,
180 ):
181 """
182 Adds a list of (dim, tag) tuples to the dimTags object.
184 dimTagsListForPG is also accepted as an argument because sometimes, the stored
185 dimTags and the current dimTags in the geometry generation can be different. For
186 example, if volume 61 is deleted, the other volume tags (62, 63, ...) won't
187 shift back. However, after saving the geometry as a BREP file and rereading it,
188 the volume tags will be shifted back. In this case, the stored dimTags should be
189 shifted as well. But to create the correct physical region in the
190 geometry-creating process (which won't be saved in the BREP file, just used for
191 debugging purposes), another optional dimTagsListForPG argument is accepted.
193 :param dimTagsList: list of (dim, tag) tuples
194 :type dimTagsList: list[tuple[int, int]]
195 :param dimTagsListForPG: list of (dim, tag) tuples for physical groups
196 (default: None). If dimTagsListForPG is None, dimTagsList will be used for
197 physical groups as well.
198 :type dimTagsListForPG: list[tuple[int, int]], optional
200 """
201 if not isinstance(dimTagsList, list):
202 dimTagsList = [dimTagsList]
204 if not all(isinstance(element, tuple) for element in dimTagsList):
205 raise TypeError("Dim tags must be a list of tuples!")
207 # Sometimes, negative entities can be added for topology.
208 for i, v in enumerate(dimTagsList):
209 if v[1] < 0:
210 dimTagsList[i] = (v[0], -v[1])
212 # Add dim tags if they are not already added:
213 for v in dimTagsList:
214 if v not in self.allDimTags:
215 self.dimTags[v[0]].append(v)
216 if dimTagsListForPG is None:
217 self.dimTagsForPG[v[0]].append(v)
218 self.allDimTags.append(v)
220 if dimTagsListForPG is not None:
221 if not isinstance(dimTagsListForPG, list):
222 dimTagsListForPG = [dimTagsListForPG]
224 if not all(isinstance(element, tuple) for element in dimTagsListForPG):
225 raise TypeError("Dim tags must be a list of tuples!")
227 for i, v in enumerate(dimTagsListForPG):
228 if v[1] < 0:
229 dimTagsListForPG[i] = (v[0], -v[1])
230 for v in dimTagsListForPG:
231 if v not in self.dimTagsForPG:
232 self.dimTagsForPG[v[0]].append(v)
234 def addWithTags(self, dim: int, tags: List[int], tagsForPG: List[int] = None):
235 """
236 Adds a list of tags with a specific dimension to the dimTags object. The
237 explanation of the tagsForPG argument is given in the add() method.
239 :param dim: dimension of the tags
240 :type dim: int
241 :param tags: list of tags
242 :type tags: list[tuple[int, int]]
243 :param tagsForPG: list of tags for physical groups (default: None)
244 :type tagsForPG: list[tuple[int, int]], optional
246 """
247 if not isinstance(tags, list):
248 tags = [tags]
250 if not isinstance(dim, int):
251 raise TypeError("Dimension must be an integer!")
253 if not all(isinstance(element, int) for element in tags):
254 raise TypeError("Tags must be a list of integers!")
256 dims = [dim] * len(tags)
257 dimTagsList = list(zip(dims, tags))
259 if tagsForPG is not None:
260 if not isinstance(tagsForPG, list):
261 tagsForPG = [tagsForPG]
263 if not all(isinstance(element, int) for element in tagsForPG):
264 raise TypeError("Tags must be a list of integers!")
266 dimTagsListForPG = list(zip(dims, tagsForPG))
267 self.add(dimTagsList, dimTagsListForPG)
268 else:
269 self.add(dimTagsList)
271 def getDimTags(self, dim: int = None) -> List[Tuple[int, int]]:
272 """
273 Returns the stored list of (dim, tag) tuples with a specific dimension. If dim
274 is not specified, returns all the (dim, tag) tuples.
276 :param dim: dimension of the tags to be returned (default: None)
277 :type dim: int, optional
278 :return: list of (dim, tag) tuples
279 :rtype: list[tuple[int, int]]
280 """
281 if dim is None:
282 return self.dimTags[0] + self.dimTags[1] + self.dimTags[2] + self.dimTags[3]
283 else:
284 return self.dimTags[dim]
286 def getDimTagsForPG(self, dim: int = None) -> List[Tuple[int, int]]:
287 """
288 Returns the stored list of (dim, tag) tuples for physical groups with a specific
289 dimension. If dim is not specified, returns all the (dim, tag) tuples for
290 physical groups.
292 :param dim: dimension of the tags to be returned (default: None)
293 :type dim: int, optional
294 :return: list of (dim, tag) tuples for physical groups
295 :rtype: list[tuple[int, int]]
296 """
297 if dim is None:
298 return (
299 self.dimTagsForPG[0]
300 + self.dimTagsForPG[1]
301 + self.dimTagsForPG[2]
302 + self.dimTagsForPG[3]
303 )
304 else:
305 return self.dimTagsForPG[dim]
307 def getTags(self, dim: int, forPhysicalGroup=False) -> List[int]:
308 """
309 Returns the stored list of tags with a specific dimension.
311 :param dim: dimension of the tags to be returned
312 :type dim: int
313 :return: list of tags
314 :rtype: list[int]
315 """
316 if forPhysicalGroup:
317 return [v[1] for v in self.dimTagsForPG[dim]]
318 else:
319 return [v[1] for v in self.dimTags[dim]]
321 def getExtrusionTop(self, dim=3):
322 """
323 Returns the top surfaces, lines, or points of an extrusion operation if the
324 dimTags object contains the tags of an extrusion operation.
325 gmsh.model.occ.extrusion() function returns all the entities that are created
326 by the extrusion as a dim tags list. The first element is always the top
327 surface, the second is the volume. However, when more than one surface is
328 extruded, extracting the top surfaces is not trivial. This function returns
329 the top surfaces of an extrusion operation. It does that by finding the entities
330 right before the volume entities. If dim is 2, the function will return the top
331 curves of an extrusion operation. If dim is 1, the function will return the top
332 points of an extrusion.
334 :param dim: dimension of the entity that is being created (default: 3)
335 :type dim: int, optional
336 :return: list of (dim, tag) tuples of the top entities of an extrusion operation
337 :rtype: list[tuple[int, int]]
338 """
340 topSurfaces = []
341 for index, dimTag in enumerate(self.allDimTags):
342 if dimTag[0] == dim:
343 topSurfaces.append(self.allDimTags[index - 1])
345 return topSurfaces
347 def getExtrusionSide(self, dim=3):
348 """
349 Returns the side surfaces, lines, or points of an extrusion operation if the
350 dimTags object contains the tags of an extrusion operation.
351 gmsh.model.occ.extrusion() function returns all the entities that are created
352 by the extrusion as a dim tags list. The first element is always the top
353 surface, the second is the volume. The other elements are the side surfaces.
354 However, when more than one surface is extruded, extracting the side surfaces
355 is not trivial. This function returns the side surfaces of an extrusion
356 operation. It does that by finding returning all the entities except the top
357 surface and the volume. If dim is 2, the function will return the side curves of
358 an extrusion operation.
360 :param dim: dimension of the entity that is being created (default: 3)
361 :type dim: int, optional
362 :return: list of (dim, tag) tuples of the side entities of an extrusion operation
363 :rtype: list[tuple[int, int]]
364 """
365 sideSurfaces = []
366 sideSurfaceStartIndex = None
367 for index, dimTag in enumerate(self.allDimTags):
368 if dimTag[0] == dim:
369 if sideSurfaceStartIndex is not None:
370 sideSurfaces.append(
371 self.allDimTags[sideSurfaceStartIndex : index - 1]
372 )
373 sideSurfaceStartIndex = index + 1
374 else:
375 sideSurfaceStartIndex = index + 1
377 sideSurfaces.append(self.allDimTags[sideSurfaceStartIndex:])
379 return sideSurfaces
381 def __add__(self, other):
382 """
383 Adds two dimTags objects and returns a new dimTags object with the same save and
384 name attirbues of the first dimTags object.
386 It might cause bugs because of the recursive behavior of the
387 @storageUpdateRequired decorator. Use with caution. Currently only used by
388 dimTagsStorage.updateDimTags method.
390 :param other: dimTags object to be added
391 :type other: dimTags
392 :return: dimTags object with the sum of the two dimTags objects
393 :rtype: dimTags
394 """
395 result = dimTags()
396 result.name = self.name
397 result.parentName = self.parentName
398 result.physicalTag = self.physicalTag
399 result.dimTags = self.dimTags
400 result.dimTagsForPG = self.dimTagsForPG
401 result.allDimTags = self.allDimTags
403 result.add(other.allDimTags)
404 result.save = self.save
405 return result
407 def __repr__(self):
408 """
409 Returns the string representation of the dimTags object. If the dimTags object
410 is saved, it will return "SAVED: name". If the dimTags object is not saved, it
411 will return "NOT SAVED: name".
413 dimTags objects are used as dictionary keys throughout the code. This
414 representation makes debugging easier.
416 :return: string representation of the dimTags object
417 :rtype: str
418 """
419 if self.save:
420 return "SAVED: " + self.name
421 else:
422 return "NOT SAVED: " + self.name
425class dimTagsStorage:
426 """
427 This is a global class to store the dimTags of important entities in the model.
428 Every dimTags instance with self.save = True will be stored in this class. Later,
429 the storage will be used to generate the volume information (*.vi) file. *.vi file
430 will be used for generating the physical regions in the meshing part.
432 Users should not use this class directly. Instead, they should use the dimTags
433 class. If they assign save = True to the dimTags instance, it will be stored in this
434 class.
436 This class is a singleton class. It means that there will be only one instance of
437 this class in the whole module. This is done to be able to use the same storage
438 throughout this module. See the singleton design pattern for more information.
439 """
441 __instance = None
442 __dimTagsDict = {} # Dictionary with the names of the dimTags objects as keys and
443 # dimTags objects as values
445 def __new__(cls):
446 if cls.__instance is None:
447 cls.__instance = super().__new__(cls)
449 return cls.__instance
451 @classmethod
452 def updateDimTags(cls, dimTagsObject: dimTags):
453 """
454 Either adds or updates the dimTags object in the storage.
456 :param dimTags: dimTags object to be added or updated.
457 :type dimTags: dimTags
459 """
460 if dimTagsObject.name in cls.__dimTagsDict:
461 newDimTags = dimTagsObject + cls.__dimTagsDict[dimTagsObject.name]
462 cls.__dimTagsDict[dimTagsObject.name] = newDimTags
463 else:
464 cls.__dimTagsDict[dimTagsObject.name] = dimTagsObject
466 @classmethod
467 def getDimTagsObject(cls, names: List[str]):
468 """
469 Returns the dimTags object with the given names.
471 :param names: names of the dimTags objects.
472 :type names: list[str]
473 :return: dimTags objects with the given name.
474 :rtype: list[dimTags]
475 """
476 if not isinstance(names, list):
477 names = [names]
479 dimTagsObjects = []
480 for name in names:
481 if name in cls.__dimTagsDict:
482 dimTagsObjects.append(cls.__dimTagsDict[name])
484 return dimTagsObjects
486 @classmethod
487 def getDimTags(cls, names: List[str], dim: int = None) -> List[Tuple[int, int]]:
488 """
489 Returns the stored list of (dim, tag) tuples with dimension a specific dimenions
490 and names. If dim is not specified, all the stored (dim, tag) tuples under the
491 given names will be returned.
493 :param names: names of the dimTags object that will be returned
494 :type names: list[str]
495 :param dim: dimension of the (dim, tag) tuples to be returned (default: None).
496 If dim is None, all the stored (dim, tag) tuples under the name names will
497 be returned.
498 :type dim: int, optional
499 :return: list of (dim, tag) tuples
500 """
501 if not isinstance(names, list):
502 names = [names]
504 dimTagsResult = []
505 for name in names:
506 dimTagsResult.extend(cls.__dimTagsDict[name].getDimTags(dim))
508 return dimTagsResult
510 @classmethod
511 def getTags(cls, names: List[str], dim: int) -> List[int]:
512 """
513 Returns the stored list of tags with dimension a specific dimension and names.
515 :param names: names of the dimTags objects
516 :type names: list[str]
517 :param dim: dimension of the tags to be returned
518 :type dim: int
519 :return: list of tags
520 :rtype: list[int]
521 """
522 dimTags = cls.getDimTags(names, dim)
523 tags = [dimTag[1] for dimTag in dimTags]
525 return tags
527 @classmethod
528 def getDimTagsDict(
529 cls, forPhysicalGroups=False
530 ) -> Dict[str, List[Tuple[int, int]]]:
531 """
532 Returns a dictionary with the names of the dimTags objects as keys and the
533 stored list of (dim, tag) tuples as values. This method is used to generate the
534 .vi file. If forPhysicalGroups is True, the dimTags for physical groups will be
535 returned instead of the dimTags.
537 :param forPhysicalGroups: True if the dimTags for physical groups should be
538 returned, False otherwise (default: False)
539 :type forPhysicalGroups: bool, optional
540 :return: dictionary with the names of the dimTags objects as keys and the
541 stored list of (dim, tag) tuples as values
542 :rtype: dict[str, list[tuple[int, int]]]
543 """
544 dictionary = {}
545 for name, dimTagsObject in cls.__dimTagsDict.items():
546 if dimTagsObject.parentName is not None:
547 if dimTagsObject.parentName in dictionary:
548 if forPhysicalGroups:
549 dictionary[dimTagsObject.parentName].extend(
550 dimTagsObject.getDimTagsForPG()
551 )
552 else:
553 dictionary[dimTagsObject.parentName].extend(
554 dimTagsObject.getDimTags()
555 )
556 else:
557 if forPhysicalGroups:
558 dictionary[dimTagsObject.parentName] = (
559 dimTagsObject.getDimTagsForPG()
560 )
561 else:
562 dictionary[dimTagsObject.parentName] = (
563 dimTagsObject.getDimTags()
564 )
565 if forPhysicalGroups:
566 dictionary[name] = dimTagsObject.getDimTagsForPG()
567 else:
568 dictionary[name] = dimTagsObject.getDimTags()
570 return dictionary
572 @classmethod
573 def getAllStoredDimTags(cls) -> List[Tuple[int, int]]:
574 """
575 Returns a list of all the stored (dim, tag) tuples, regardless of the name of
576 the dimTags object (i.e. all the dimTags objects are merged into one list).
578 :return: list of all the stored (dim, tag) tuples.
579 :rtype: list[tuple[int, int]]
580 """
581 AllStoredDimTags = []
582 for name, dimTagsObject in cls.__dimTagsDict.items():
583 AllStoredDimTags.extend(dimTagsObject.getDimTags())
585 return AllStoredDimTags
587 @classmethod
588 def clear(cls):
589 """
590 Clears the dimTagsStorage class.
593 """
594 cls.__instance = None
595 cls.__dimTagsDict = (
596 {}
597 ) # Dictionary with the names of the dimTags objects as keys and
598 # dimTags objects as values
601class coordinate(Enum):
602 """
603 A class to specify coordinate types easily.
604 """
606 rectangular = 0
607 cylindrical = 1
608 spherical = 2
611class direction(Enum):
612 """
613 A class to specify direction easily.
614 """
616 ccw = 0
617 cw = 1
620class point:
621 """
622 This is a class for creating points in GMSH. It supports rectangular and cylindrical
623 coordinates. Moreover, vector operations are supported.
625 :param r0: x, r, or r (default: 0.0)
626 :type r0: float, optional
627 :param r1: y, theta, or theta (default: 0.0)
628 :type r1: float, optional
629 :param r2: z, z, or phi (default: 0.0)
630 :type r2: float, optional
631 :param type: coordinate type (default: coordinate.rectangular)
632 :type type: coordinate, optional
633 """
635 def __init__(self, r0=0.0, r1=0.0, r2=0.0, type=coordinate.rectangular) -> None:
636 if type is coordinate.rectangular:
637 self.x = r0
638 self.y = r1
639 self.z = r2
641 self.r = math.sqrt(self.x**2 + self.y**2)
642 self.theta = math.atan2(self.y, self.x)
643 elif type is coordinate.cylindrical:
644 self.r = r0
645 self.theta = r1
646 self.x = self.r * math.cos(self.theta)
647 self.y = self.r * math.sin(self.theta)
648 self.z = r2
649 elif type is coordinate.spherical:
650 raise ValueError("Spherical coordinates are not supported yet!")
651 else:
652 raise ValueError("Improper coordinate type value!")
654 self.tag = gmsh.model.occ.addPoint(self.x, self.y, self.z)
656 def __repr__(self):
657 """
658 Returns the string representation of the point.
660 :return: string representation of the point
661 :rtype: str
662 """
663 return "point(%r, %r, %r, %r)" % (self.x, self.y, self.z, self.type)
665 def __abs__(self):
666 """
667 Returns the magnitude of the point vector.
669 :return: the magnitude of the point vector
670 :rtype: float
671 """
672 return math.hypot(self.x, self.y, self.z)
674 def __add__(self, other):
675 """
676 Returns the summation of two point vectors.
678 :param other: point vector to be added
679 :type other: point
680 :return: the summation of two point vectors
681 :rtype: point
682 """
683 x = self.x + other.x
684 y = self.y + other.y
685 z = self.z + other.z
686 return point(x, y, z, coordinate.rectangular)
688 def __mul__(self, scalar):
689 """
690 Returns the product of a point vector and a scalar.
692 :param scalar: a scalar value
693 :type scalar: float
694 :return: point
695 :rtype: point
696 """
697 return point(
698 self.x * scalar,
699 self.y * scalar,
700 self.z * scalar,
701 coordinate.rectangular,
702 )
705class spiralCurve:
706 """
707 A class to create a spiral curves parallel to XY plane in GMSH. The curve is defined
708 by a spline and it is divided into sub-curves. Sub-curves are used because it makes
709 the geometry creation process easier.
711 :param innerRadius: inner radius
712 :type innerRadius: float
713 :param gap: gap after each turn
714 :type gap: float
715 :param turns: number of turns
716 :type turns: float
717 :param z: z coordinate
718 :type z: float
719 :param initialTheta: initial theta angle in radians
720 :type initialTheta: float
721 :param direction: direction of the spiral (default: direction.ccw)
722 :type direction: direction, optional
723 :param cutPlaneNormal: normal vector of the plane that will cut the spiral curve
724 (default: None)
725 :type cutPlaneNormal: tuple[float, float, float], optional
726 """
728 # If the number of points used per turn is n, then the number of sections per turn
729 # is n-1. They set the resolution of the spiral curve. It sets the limit of the
730 # precision of the float number of turns that can be used to create the spiral
731 # curve. The value below might be modified in Geometry.__init__ method.
732 sectionsPerTurn = 16
734 # There will be curvesPerTurn curve entities per turn. It will be effectively the
735 # number of volumes per turn in the end. The value below might be modified in
736 # Geometry.__init__ method.
737 curvesPerTurn = 2
739 def __init__(
740 self,
741 innerRadius,
742 gap,
743 turns,
744 z,
745 initialTheta,
746 transitionNotchAngle,
747 direction=direction.ccw,
748 cutPlaneNormal=Tuple[float, float, float],
749 ) -> None:
750 spt = self.sectionsPerTurn # just to make the code shorter
751 self.turnRes = 1 / spt # turn resolution
752 cpt = self.curvesPerTurn # just to make the code shorter
753 self.turns = turns
755 # =============================================================================
756 # GENERATING POINTS STARTS ====================================================
757 # =============================================================================
759 # Calculate the coordinates of the points that define the spiral curve:
760 if direction is direction.ccw:
761 # If the spiral is counter-clockwise, the initial theta angle decreases,
762 # and r increases as the theta angle decreases.
763 multiplier = 1
764 elif direction is direction.cw:
765 # If the spiral is clockwise, the initial theta angle increases, and r
766 # increases as the theta angle increases.
767 multiplier = -1
769 NofPointsPerTurn = int(spt + 1)
770 thetaArrays = []
771 for turn in range(1, int(self.turns) + 1):
772 thetaArrays.append(
773 np.linspace(
774 initialTheta + (turn - 1) * 2 * math.pi * multiplier,
775 initialTheta + (turn) * 2 * math.pi * multiplier,
776 NofPointsPerTurn,
777 )
778 )
780 thetaArrays.append(
781 np.linspace(
782 initialTheta + (turn) * 2 * math.pi * multiplier,
783 initialTheta + (self.turns) * 2 * math.pi * multiplier,
784 round(spt * (self.turns - turn) + 1),
785 )
786 )
788 if cutPlaneNormal is not None:
789 # If the cutPlaneNormal is specified, the spiral curve will be cut by a
790 # plane that is normal to the cutPlaneNormal vector and passes through the
791 # origin.
793 alpha = math.atan2(cutPlaneNormal[1], cutPlaneNormal[0]) - math.pi / 2
794 alpha2 = alpha + math.pi
796 listOfBreakPoints = []
797 for turn in range(1, int(self.turns) + 2):
798 breakPoint1 = alpha + (turn - 1) * 2 * math.pi * multiplier
799 breakPoint2 = alpha2 + (turn - 1) * 2 * math.pi * multiplier
800 if (
801 breakPoint1 > initialTheta
802 and breakPoint1 < initialTheta + 2 * math.pi * self.turns
803 ):
804 listOfBreakPoints.append(breakPoint1)
805 if (
806 breakPoint2 > initialTheta
807 and breakPoint2 < initialTheta + 2 * math.pi * self.turns
808 ):
809 listOfBreakPoints.append(breakPoint2)
811 thetaArrays.append(np.array(listOfBreakPoints))
813 theta = np.concatenate(thetaArrays)
814 theta = np.round(theta, 10)
815 theta = np.unique(theta)
816 theta = np.sort(theta)
817 theta = theta[::multiplier]
819 r = innerRadius + (theta - initialTheta) / (2 * math.pi) * (gap) * multiplier
820 z = np.ones(theta.shape) * z
822 # Create the points and store their tags:
823 points = [] # point objects
824 pointTags = [] # point tags
825 breakPointObjectsDueToCutPlane = [] # only used if cutPlaneNormal is not None
826 breakPointTagsDueToCutPlane = [] # only used if cutPlaneNormal is not None
827 pointObjectsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None
828 pointTagsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None
829 breakPointObjectsDueToTransition = []
830 breakPointTagsDueToTransition = []
832 for j in range(len(theta)):
833 pointObject = point(r[j], theta[j], z[j], coordinate.cylindrical)
834 points.append(pointObject)
835 pointTags.append(pointObject.tag)
836 if cutPlaneNormal is not None:
837 if theta[j] in listOfBreakPoints:
838 breakPointObjectsDueToCutPlane.append(pointObject)
839 breakPointTagsDueToCutPlane.append(pointObject.tag)
840 else:
841 pointObjectsWithoutBreakPoints.append(pointObject)
842 pointTagsWithoutBreakPoints.append(pointObject.tag)
844 # identify if the point is a break point due to the layer transition:
845 angle1 = initialTheta + (2 * math.pi - transitionNotchAngle) * multiplier
846 angle2 = (
847 initialTheta
848 + ((self.turns % 1) * 2 * math.pi + transitionNotchAngle) * multiplier
849 )
850 if math.isclose(
851 math.fmod(theta[j], 2 * math.pi), angle1, abs_tol=1e-6
852 ) or math.isclose(math.fmod(theta[j], 2 * math.pi), angle2, abs_tol=1e-6):
853 breakPointObjectsDueToTransition.append(pointObject)
854 breakPointTagsDueToTransition.append(pointObject.tag)
856 # =============================================================================
857 # GENERATING POINTS ENDS ======================================================
858 # =============================================================================
860 # =============================================================================
861 # GENERATING SPLINES STARTS ===================================================
862 # =============================================================================
864 # Create the spline with the points:
865 spline = gmsh.model.occ.addSpline(pointTags)
867 # Split the spline into sub-curves:
868 sectionsPerCurve = int(spt / cpt)
870 # Create a list of point tags that will split the spline:
871 # Basically, they are the points to divide the spirals into sectionsPerCurve
872 # turns. However, some other points are also included to support the float
873 # number of turns. It is best to visually look at the divisions with
874 # gmsh.fltk.run() to understand why the split points are chosen the way they are
875 # selected.
876 if cutPlaneNormal is None:
877 pointObjectsWithoutBreakPoints = points
878 pointTagsWithoutBreakPoints = pointTags
880 splitPointTags = list(
881 set(pointTagsWithoutBreakPoints[:-1:sectionsPerCurve])
882 | set(pointTagsWithoutBreakPoints[-spt - 1 :: -spt])
883 | set(breakPointTagsDueToCutPlane)
884 | set(breakPointTagsDueToTransition)
885 )
886 splitPointTags = sorted(splitPointTags)
887 # Remove the first element of the list (starting point):
888 _, *splitPointTags = splitPointTags
890 # Also create a list of corresponding point objects:
891 splitPoints = list(
892 set(pointObjectsWithoutBreakPoints[:-1:sectionsPerCurve])
893 | set(pointObjectsWithoutBreakPoints[-spt - 1 :: -spt])
894 | set(breakPointObjectsDueToCutPlane)
895 | set(breakPointObjectsDueToTransition)
896 )
897 splitPoints = sorted(splitPoints, key=lambda x: x.tag)
898 # Remove the first element of the list (starting point):
899 _, *splitPoints = splitPoints
901 # Split the spline:
902 dims = [0] * len(splitPointTags)
903 _, splines = gmsh.model.occ.fragment(
904 [(1, spline)],
905 list(zip(dims, splitPointTags)),
906 removeObject=True,
907 removeTool=True,
908 )
909 splines = splines[0]
910 self.splineTags = [j for _, j in splines]
912 # Note the turn number of each spline. This will be used in getSplineTag and
913 # getSplineTags methods.
914 self.splineTurns = []
915 for i in range(len(self.splineTags)):
916 if i == 0:
917 startPoint = points[0]
918 endPoint = splitPoints[0]
919 elif i == len(self.splineTags) - 1:
920 startPoint = splitPoints[-1]
921 endPoint = points[-1]
922 else:
923 startPoint = splitPoints[i - 1]
924 endPoint = splitPoints[i]
926 startTurn = (startPoint.theta - initialTheta) / (2 * math.pi)
927 startTurn = round(startTurn / self.turnRes) * self.turnRes
928 endTurn = (endPoint.theta - initialTheta) / (2 * math.pi)
929 endTurn = round(endTurn / self.turnRes) * self.turnRes
931 if direction is direction.ccw:
932 self.splineTurns.append((startTurn, endTurn))
933 else:
934 self.splineTurns.append((-startTurn, -endTurn))
936 # Check if splineTurn tuples starts with the small turn number:
937 for i in range(len(self.splineTurns)):
938 self.splineTurns[i] = sorted(self.splineTurns[i])
940 # =============================================================================
941 # GENERATING SPLINES ENDS =====================================================
942 # =============================================================================
944 # Find start and end points of the spiral curve:
945 gmsh.model.occ.synchronize() # synchronize the model to make getBoundary work
946 self.startPointTag = gmsh.model.getBoundary([(1, self.getSplineTag(0))])[1][1]
947 self.endPointTag = gmsh.model.getBoundary(
948 [(1, self.getSplineTag(self.turns, endPoint=True))]
949 )[1][1]
951 def getSplineTag(self, turn, endPoint=False):
952 """
953 Returns the spline tag at a specific turn. It returns the spline tag of the
954 section that is on the turn except its end point.
956 :param turn: turn number (it can be a float)
957 :type turn: float
958 :param endPoint: if True, return the spline tag of the section that is on the
959 turn including its end point but not its start point (default: False)
960 :type endPoint: bool, optional
961 :return: spline tag
962 """
963 if endPoint:
964 for i, r in enumerate(self.splineTurns):
965 if r[0] + (self.turnRes / 2) < turn <= r[1] + (self.turnRes / 2):
966 return self.splineTags[i]
967 else:
968 for i, r in enumerate(self.splineTurns):
969 if r[0] - (self.turnRes / 2) <= turn < r[1] - (self.turnRes / 2):
970 return self.splineTags[i]
972 def getPointTag(self, turn, endPoint=False):
973 """
974 Returns the point object at a specific turn.
976 :param turn: turn number (it can be a float)
977 :type turn: float
978 :return: point object
979 :rtype: point
980 """
981 if turn < 0 or turn > self.turns:
982 raise ValueError("Turn number is out of range!")
984 if turn == 0:
985 return self.startPointTag
986 elif turn == self.turns:
987 return self.endPointTag
988 else:
989 curveTag = self.getSplineTag(turn, endPoint=endPoint)
990 if endPoint:
991 points = gmsh.model.getBoundary([(1, curveTag)])
992 return points[1][1]
993 else:
994 points = gmsh.model.getBoundary([(1, curveTag)])
995 return points[0][1]
997 def getSplineTags(self, turnStart=None, turnEnd=None):
998 """
999 Get the spline tags from a specific turn to another specific turn. If turnStart
1000 and turnEnd are not specified, it returns all the spline tags.
1002 :param turnStart: start turn number (it can be a float) (default: None)
1003 :type turnStart: float, optional
1004 :param turnEnd: end turn number (it can be a float) (default: None)
1005 :type turnEnd: float, optional
1006 :return: spline tags
1007 :rtype: list[int]
1008 """
1009 if turnStart is None and turnEnd is None:
1010 return self.splineTags
1011 elif turnStart is None or turnEnd is None:
1012 raise ValueError(
1013 "turnStart and turnEnd must be both specified or both not specified."
1014 " You specified only one of them."
1015 )
1016 else:
1017 start = self.splineTags.index(self.getSplineTag(turnStart, False))
1018 end = self.splineTags.index(self.getSplineTag(turnEnd, True)) + 1
1019 return self.splineTags[start:end]
1022class spiralSurface:
1023 """
1024 This is a class to create a spiral surface parallel to the XY plane in GMSH. If
1025 thinShellApproximation is set to False, it creates two spiral surfaces parallel to
1026 the XY plane, and their inner and outer curve loops in GMSH. One of the surfaces is
1027 the main surface specified, which is the winding surface, and the other is the
1028 contact layer surface (the gap between the winding surface). If thinShellApproximation
1029 is set to True, it creates only one spiral surface that touches each other
1030 (conformal).
1032 Note that surfaces are subdivided depending on the spiral curve divisions because
1033 this is required for the thin-shell approximation. Otherwise, when
1034 thinShellApproximation is set to True, it would be a disc rather than a spiral since
1035 it touches each other. However, this can be avoided by dividing the surfaces into
1036 small parts and making them conformal. Dividing the surfaces is not necessary when
1037 thinShellApproximation is set to false, but in order to use the same logic with TSA,
1038 it is divided anyway.
1040 :param innerRadius: inner radius
1041 :type innerRadius: float
1042 :param thickness: thickness
1043 :type thickness: float
1044 :param contactLayerThickness: contact layer thickness
1045 :type contactLayerThickness: float
1046 :param turns: number of turns
1047 :type turns: float
1048 :param z: z coordinate
1049 :type z: float
1050 :param initialTheta: initial theta angle in radians
1051 :type initialTheta: float
1052 :param spiralDirection: direction of the spiral (default: direction.ccw)
1053 :type spiralDirection: direction, optional
1054 :param thinShellApproximation: if True, the thin shell approximation is used
1055 (default: False)
1056 :type thinShellApproximation: bool, optional
1057 :param cutPlaneNormal: normal vector of the plane that will cut the spiral surface
1058 (default: None)
1059 :type cutPlaneNormal: tuple[float, float, float], optional
1060 """
1062 def __init__(
1063 self,
1064 innerRadius,
1065 thickness,
1066 contactLayerThickness,
1067 turns,
1068 z,
1069 initialTheta,
1070 transitionNotchAngle,
1071 spiralDirection=direction.ccw,
1072 thinShellApproximation=False,
1073 cutPlaneNormal=None,
1074 ) -> None:
1075 r_i = innerRadius
1076 t = thickness
1077 theta_i = initialTheta
1078 self.theta_i = theta_i
1080 self.surfaceTags = []
1081 self.contactLayerSurfaceTags = []
1083 self.direction = spiralDirection
1084 self.tsa = thinShellApproximation
1085 self.transitionNotchAngle = transitionNotchAngle
1086 # =============================================================================
1087 # GENERATING SPIRAL CURVES STARTS =============================================
1088 # =============================================================================
1089 if thinShellApproximation:
1090 # Create only one spiral curve because the thin shell approximation is used:
1091 # Winding thickness is increased slightly to ensure that the outer radius
1092 # would be the same without the thin shell approximation.
1093 turns = (
1094 turns + 1
1095 ) # for TSA, spiral has (n+1) turns but spiral surface has n
1096 spiral = spiralCurve(
1097 r_i,
1098 t + contactLayerThickness * (turns - 1) / (turns),
1099 turns,
1100 z,
1101 theta_i,
1102 transitionNotchAngle,
1103 spiralDirection,
1104 cutPlaneNormal=cutPlaneNormal,
1105 )
1107 # These are created to be able to use the same code with TSA and without TSA:
1108 innerSpiral = spiral
1109 outerSpiral = spiral
1110 else:
1111 # Create two spiral curves because the thin shell approximation is not used:
1112 innerSpiral = spiralCurve(
1113 r_i - contactLayerThickness,
1114 t + contactLayerThickness,
1115 turns + 1,
1116 z,
1117 theta_i,
1118 transitionNotchAngle,
1119 spiralDirection,
1120 cutPlaneNormal=cutPlaneNormal,
1121 )
1122 outerSpiral = spiralCurve(
1123 r_i,
1124 t + contactLayerThickness,
1125 turns + 1,
1126 z,
1127 theta_i,
1128 transitionNotchAngle,
1129 spiralDirection,
1130 cutPlaneNormal=cutPlaneNormal,
1131 )
1133 self.innerSpiral = innerSpiral
1134 self.outerSpiral = outerSpiral
1135 self.turns = turns
1136 # =============================================================================
1137 # GENERATING SPIRAL CURVES ENDS ===============================================
1138 # =============================================================================
1140 # =============================================================================
1141 # GENERATING SURFACES STARTS ==================================================
1142 # =============================================================================
1143 endLines = []
1144 endInsLines = []
1146 # This is used to check if all the contact layers are finished:
1147 allContactLayersAreFinished = False
1149 # Itterate over the spline tags:
1150 for i in range(len(innerSpiral.splineTags)):
1151 if thinShellApproximation:
1152 # The current spline will be the inner spline:
1153 innerSplineTag = spiral.splineTags[i]
1155 # Find the spline tag of the outer spline by finding the spline tag of
1156 # the spline that is exactly on the next turn:
1158 # Note the turn number of the current spline's start point:
1159 startTurn = spiral.splineTurns[i][0]
1161 if startTurn + 1 + 1e-4 > turns:
1162 # If the current spline is on the outer surface, break the loop,
1163 # because the whole surface is finished:
1164 break
1166 # Find the outer spline tag:
1167 isItBroken = True
1168 for j, turnTuple in enumerate(spiral.splineTurns):
1169 # Equality can not be checked with == because of the floating point
1170 # errors:
1171 if abs(turnTuple[0] - 1 - startTurn) < 1e-4:
1172 outerSplineTag = spiral.splineTags[j]
1173 isItBroken = False
1174 break
1176 if isItBroken:
1177 raise RuntimeError(
1178 "Something went wrong while creating the spiral surface. Outer"
1179 f" spline tag of {innerSplineTag} could not be found for TSA."
1180 )
1182 else:
1183 # Store the tags of the current splines:
1184 innerSplineTag = innerSpiral.splineTags[i]
1185 outerSplineTag = outerSpiral.splineTags[i]
1187 # The current outer spline will be the inner spline of the
1188 # contact layer:
1189 innerInsSplineTag = outerSpiral.splineTags[i]
1191 # Find the spline tag of the contact layer's outer spline by finding the
1192 # spline tag of the spline that is exactly on the next turn:
1194 # Note the turn number of the current spline's start point:
1195 startTurn = outerSpiral.splineTurns[i][0]
1197 if startTurn + 1 + 1e-4 > turns + 1:
1198 # If the current spline is on the outer surface, note that all the
1199 # contact layers are finished:
1200 allContactLayersAreFinished = True
1202 # Find the contact layer's outer spline tag:
1203 for j, turnTuple in enumerate(innerSpiral.splineTurns):
1204 if math.isclose(turnTuple[0], 1 + startTurn, abs_tol=1e-6):
1205 outerInsSplineTag = innerSpiral.splineTags[j]
1206 break
1208 # Create the lines to connect the two splines so that a surface can be
1209 # created:
1211 # Create start line:
1212 if i == 0:
1213 # If it is the first spline, start line should be created.
1215 # Create points:
1216 isStartPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]
1217 if thinShellApproximation:
1218 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]
1219 else:
1220 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]
1222 # Create the line:
1223 startLine = gmsh.model.occ.addLine(osStartPoint, isStartPoint)
1224 firstStartLine = startLine
1226 # Create lines for the contact layer if the thin shell approximation is not
1227 # used:
1228 if not thinShellApproximation and not allContactLayersAreFinished:
1229 isInsStartPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[
1230 1
1231 ][1]
1232 osInsStartPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[
1233 0
1234 ][1]
1236 # Create the line:
1237 startInsLine = gmsh.model.occ.addLine(
1238 osInsStartPoint, isInsStartPoint
1239 )
1240 firstInsStartLine = startInsLine
1242 else:
1243 # If it is not the first spline, the start line is the end line of the
1244 # previous surface. This guarantees that the surfaces are connected
1245 # (conformality).
1246 startLine = endLines[i - 1]
1248 # Do the same for the contact layer if the thin shell approximation is not
1249 # used:
1250 if not thinShellApproximation and not allContactLayersAreFinished:
1251 startInsLine = endInsLines[i - 1]
1253 # Create end line:
1255 # Create points:
1256 # The ifs are used because getBoundary is not consistent somehow.
1257 if i == 0:
1258 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[0][1]
1259 else:
1260 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]
1262 if (not i == 0) or thinShellApproximation:
1263 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]
1264 else:
1265 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]
1267 # Create the line:
1268 endLine = gmsh.model.occ.addLine(isEndPoint, osEndPoint)
1269 endLines.append(endLine)
1271 # Create lines for the contact layer if the thin shell approximation is not
1272 # used:
1273 if not thinShellApproximation and not allContactLayersAreFinished:
1274 if i == 0:
1275 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[0][
1276 1
1277 ]
1278 else:
1279 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[1][
1280 1
1281 ]
1283 osInsEndPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[1][1]
1285 # Create the line:
1286 endInsLine = gmsh.model.occ.addLine(isInsEndPoint, osInsEndPoint)
1287 endInsLines.append(endInsLine)
1289 # Create the surface:
1290 curveLoop = gmsh.model.occ.addCurveLoop(
1291 [startLine, innerSplineTag, endLine, outerSplineTag]
1292 )
1293 self.surfaceTags.append(gmsh.model.occ.addPlaneSurface([curveLoop]))
1295 # Create the surface for the contact layer if the thin shell approximation is
1296 # not used:
1297 if not thinShellApproximation and not allContactLayersAreFinished:
1298 curveLoop = gmsh.model.occ.addCurveLoop(
1299 [startInsLine, innerInsSplineTag, endInsLine, outerInsSplineTag]
1300 )
1301 self.contactLayerSurfaceTags.append(
1302 gmsh.model.occ.addPlaneSurface([curveLoop])
1303 )
1305 # =============================================================================
1306 # GENERATING SURFACES ENDS ====================================================
1307 # =============================================================================
1309 # =============================================================================
1310 # GENERATING CURVE LOOPS STARTS ===============================================
1311 # =============================================================================
1313 # Create the inner and outer curve loops (for both TSA == True and TSA == False):
1315 # VERY IMPORTANT NOTES ABOUT THE DIRECTION OF THE CURVE LOOPS
1316 # 1- GMSH doesn't like duplicates. Or the user doesn't like duplicates if they
1317 # want conformality. Actually, it's a positive thing about debugging because
1318 # you can always use `Geometry.remove_all_duplicates()` to see if there are
1319 # any duplicates. If there are, the problem is found. Solve it.
1320 # 2- The problem arises when one uses surface loops or curve loops. Because even
1321 # if you think there are no duplicates, GMSH/OCC might create some during
1322 # addCurveLoops and addSurfaceLoops operations. Even though
1323 # `geometry.remove_all_duplicates()` tells that there are duplicates, the
1324 # user doesn't suspect about addCurveLoop and addSurfaceLoop at first,
1325 # because it doesn't make sense.
1326 # 3- How you put the curves in the curve loops is very important! The same curve
1327 # loop with the same lines might cause problems if the user puts them in a
1328 # different order. For example, to create a plane surface with two curve
1329 # loops, the direction of the curve loops should be the same. That's why the
1330 # code has both innerCurveLoopTag and innerOppositeCurveLoopTag (the same
1331 # thing for the outer curve loop).
1333 # create the transition layer (notch):
1334 # Inner curve loop:
1335 notchStartPoint = innerSpiral.getPointTag(
1336 1 - transitionNotchAngle / (2 * math.pi)
1337 )
1338 notchLeftPoint = innerSpiral.getPointTag(0)
1339 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1340 notchRightLine = innerSpiral.getSplineTag(
1341 1 - transitionNotchAngle / (2 * math.pi)
1342 )
1344 if thinShellApproximation:
1345 innerStartCurves = [firstStartLine]
1346 else:
1347 innerStartCurves = [firstInsStartLine, firstStartLine]
1349 if thinShellApproximation:
1351 notchCurveLoop = gmsh.model.occ.addCurveLoop(
1352 [notchLeftLine, notchRightLine] + innerStartCurves
1353 )
1354 self.innerNotchSurfaceTags = [
1355 gmsh.model.occ.addPlaneSurface([notchCurveLoop])
1356 ]
1357 else:
1358 notchMiddlePoint = outerSpiral.getPointTag(0)
1359 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)
1361 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1362 [notchLeftLine, notchMiddleLine, firstStartLine]
1363 )
1364 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(
1365 [notchMiddleLine, notchRightLine, firstInsStartLine]
1366 )
1367 self.innerNotchSurfaceTags = [
1368 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1369 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),
1370 ]
1372 lines = innerSpiral.getSplineTags(
1373 0, 1 - transitionNotchAngle / (2 * math.pi)
1374 ) # The first turn of the spline
1375 innerCurves = lines + [notchLeftLine]
1376 self.innerNotchLeftLine = notchLeftLine
1378 self.innerStartCurves = innerStartCurves
1380 self.innerCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves)
1381 self.innerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves[::-1])
1383 # Outer curve loop:
1384 # The last turn of the spline:
1385 if thinShellApproximation:
1386 notchStartPoint = innerSpiral.getPointTag(
1387 self.turns + transitionNotchAngle / (2 * math.pi) - 1
1388 )
1389 notchLeftPoint = innerSpiral.getPointTag(self.turns)
1390 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1391 notchRightLine = innerSpiral.getSplineTag(
1392 self.turns - 1 + transitionNotchAngle / (2 * math.pi), endPoint=True
1393 )
1394 else:
1395 notchStartPoint = outerSpiral.getPointTag(
1396 self.turns + transitionNotchAngle / (2 * math.pi)
1397 )
1398 notchMiddlePoint = innerSpiral.getPointTag(self.turns + 1)
1399 notchLeftPoint = outerSpiral.getPointTag(self.turns + 1)
1400 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)
1401 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)
1402 notchRightLine = outerSpiral.getSplineTag(
1403 self.turns + transitionNotchAngle / (2 * math.pi), self.turns
1404 )
1405 if thinShellApproximation:
1406 lines = outerSpiral.getSplineTags(turns - 1, turns)
1407 else:
1408 lines = outerSpiral.getSplineTags(turns, turns + 1)
1410 if thinShellApproximation:
1411 outerEndCurves = [endLines[-1]]
1412 else:
1413 outerEndCurves = [endInsLines[-1], endLines[-1]]
1415 if thinShellApproximation:
1416 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1417 [notchLeftLine, notchRightLine, endLines[-1]]
1418 )
1419 self.outerNotchSurfaceTags = [
1420 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1421 ]
1422 else:
1423 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(
1424 [notchLeftLine, notchMiddleLine, endLines[-1]]
1425 )
1426 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(
1427 [notchMiddleLine, notchRightLine, endInsLines[-1]]
1428 )
1429 self.outerNotchSurfaceTags = [
1430 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),
1431 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),
1432 ]
1434 if thinShellApproximation:
1435 lines = innerSpiral.getSplineTags(
1436 self.turns - 1 + transitionNotchAngle / (2 * math.pi), self.turns
1437 ) # The first turn of the spline
1438 else:
1439 lines = outerSpiral.getSplineTags(
1440 self.turns + transitionNotchAngle / (2 * math.pi), self.turns + 1
1441 )
1442 outerCurves = lines + [notchLeftLine]
1443 self.outerNotchLeftLine = notchLeftLine
1445 self.outerEndCurves = outerEndCurves
1447 self.outerCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves)
1448 self.outerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves[::-1])
1449 # =============================================================================
1450 # GENERATING CURVE LOOPS ENDS =================================================
1451 # =============================================================================
1453 if not thinShellApproximation:
1454 surfaceTags = self.surfaceTags
1455 self.surfaceTags = self.contactLayerSurfaceTags
1456 self.contactLayerSurfaceTags = surfaceTags
1458 def getInnerRightPointTag(self):
1459 """
1460 Returns the point tag of the inner left point.
1462 :return: point tag
1463 :rtype: int
1464 """
1465 return self.innerSpiral.getPointTag(0)
1467 def getInnerUpperPointTag(self):
1468 """
1469 Returns the point tag of the inner right point.
1471 :return: point tag
1472 :rtype: int
1473 """
1474 if self.direction is direction.ccw:
1475 return self.innerSpiral.getPointTag(0.25)
1476 else:
1477 return self.innerSpiral.getPointTag(0.75)
1479 def getInnerLeftPointTag(self):
1480 """
1481 Returns the point tag of the inner upper point.
1483 :return: point tag
1484 :rtype: int
1485 """
1486 return self.innerSpiral.getPointTag(0.5)
1488 def getInnerLowerPointTag(self):
1489 """
1490 Returns the point tag of the inner lower point.
1492 :return: point tag
1493 :rtype: int
1494 """
1495 if self.direction is direction.ccw:
1496 return self.innerSpiral.getPointTag(0.75)
1497 else:
1498 return self.innerSpiral.getPointTag(0.25)
1500 def getOuterRightPointTag(self):
1501 """
1502 Returns the point tag of the outer left point.
1504 :return: point tag
1505 :rtype: int
1506 """
1507 if self.tsa:
1508 turns = self.turns
1509 else:
1510 turns = self.turns + 1
1511 return self.outerSpiral.getPointTag(turns, endPoint=False)
1513 def getOuterLowerPointTag(self):
1514 """
1515 Returns the point tag of the outer right point.
1517 :return: point tag
1518 :rtype: int
1519 """
1520 if self.tsa:
1521 turns = self.turns
1522 else:
1523 turns = self.turns + 1
1524 if self.direction is direction.ccw:
1525 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)
1526 else:
1527 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)
1529 def getOuterLeftPointTag(self):
1530 """
1531 Returns the point tag of the outer upper point.
1533 :return: point tag
1534 :rtype: int
1535 """
1536 if self.tsa:
1537 turns = self.turns
1538 else:
1539 turns = self.turns + 1
1540 return self.outerSpiral.getPointTag(turns - 0.5, endPoint=False)
1542 def getOuterUpperPointTag(self):
1543 """
1544 Returns the point tag of the outer lower point.
1546 :return: point tag
1547 :rtype: int
1548 """
1549 if self.tsa:
1550 turns = self.turns
1551 else:
1552 turns = self.turns + 1
1553 if self.direction is direction.ccw:
1554 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)
1555 else:
1556 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)
1558 def getInnerUpperRightCurves(self):
1559 """
1560 Returns the curve tags of the upper right curves.
1562 :return: curve tags
1563 :rtype: list[int]
1564 """
1565 if self.direction is direction.ccw:
1566 curves = self.innerSpiral.getSplineTags(0, 0.25)
1567 else:
1568 lines = self.innerSpiral.getSplineTags(
1569 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)
1570 ) # The first turn of the spline
1571 lines = lines + [self.innerNotchLeftLine]
1573 return lines
1575 return curves
1577 def getInnerUpperLeftCurves(self):
1578 """
1579 Returns the curve tags of the upper left curves.
1581 :return: curve tags
1582 :rtype: list[int]
1583 """
1584 if self.direction is direction.ccw:
1585 curves = self.innerSpiral.getSplineTags(0.25, 0.5)
1586 else:
1587 curves = self.innerSpiral.getSplineTags(0.5, 0.75)
1589 return curves
1591 def getInnerLowerLeftCurves(self):
1592 """
1593 Returns the curve tags of the lower left curves.
1595 :return: curve tags
1596 :rtype: list[int]
1597 """
1598 if self.direction is direction.ccw:
1599 curves = self.innerSpiral.getSplineTags(0.5, 0.75)
1600 else:
1601 curves = self.innerSpiral.getSplineTags(0.25, 0.5)
1603 return curves
1605 def getInnerLowerRightCurves(self):
1606 """
1607 Returns the curve tags of the lower right curves.
1609 :return: curve tags
1610 :rtype: list[int]
1611 """
1612 if self.direction is direction.ccw:
1613 lines = self.innerSpiral.getSplineTags(
1614 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)
1615 ) # The first turn of the spline
1616 lines = lines + [self.innerNotchLeftLine]
1618 return lines
1619 else:
1620 curves = self.innerSpiral.getSplineTags(0, 0.25)
1622 return curves
1624 def getOuterUpperRightCurves(self):
1625 """
1626 Returns the curve tags of the upper right curves.
1628 :return: curve tags
1629 :rtype: list[int]
1630 """
1631 if self.tsa:
1632 turns = self.turns
1633 else:
1634 turns = self.turns + 1
1636 if self.direction is direction.ccw:
1637 if self.tsa:
1638 lines = self.innerSpiral.getSplineTags(
1639 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),
1640 self.turns - 0.75,
1641 ) # The first turn of the spline
1642 else:
1643 lines = self.outerSpiral.getSplineTags(
1644 self.turns + self.transitionNotchAngle / (2 * math.pi),
1645 self.turns + 1 - 0.75,
1646 )
1647 lines = lines + [self.outerNotchLeftLine]
1649 return lines
1650 else:
1651 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)
1653 return curves
1655 def getOuterUpperLeftCurves(self):
1656 """
1657 Returns the curve tags of the lower right curves.
1659 :return: curve tags
1660 :rtype: list[int]
1661 """
1662 if self.tsa:
1663 turns = self.turns
1664 else:
1665 turns = self.turns + 1
1666 if self.direction is direction.ccw:
1667 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)
1668 else:
1669 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)
1671 return curves
1673 def getOuterLowerLeftCurves(self):
1674 """
1675 Returns the curve tags of the lower left curves.
1677 :return: curve tags
1678 :rtype: list[int]
1679 """
1680 if self.tsa:
1681 turns = self.turns
1682 else:
1683 turns = self.turns + 1
1684 if self.direction is direction.ccw:
1685 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)
1686 else:
1687 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)
1689 return curves
1691 def getOuterLowerRightCurves(self):
1692 """
1693 Returns the curve tags of the upper left curves.
1695 :return: curve tags
1696 :rtype: list[int]
1697 """
1698 if self.tsa:
1699 turns = self.turns
1700 else:
1701 turns = self.turns + 1
1702 if self.direction is direction.ccw:
1703 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)
1704 else:
1705 if self.tsa:
1706 lines = self.innerSpiral.getSplineTags(
1707 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),
1708 self.turns - 0.75,
1709 ) # The first turn of the spline
1710 else:
1711 lines = self.outerSpiral.getSplineTags(
1712 self.turns + self.transitionNotchAngle / (2 * math.pi),
1713 self.turns + 1 - 0.75,
1714 )
1715 lines = lines + [self.outerNotchLeftLine]
1717 return lines
1719 return curves
1721 def getInnerStartCurves(self):
1722 """
1723 Returns the curve tags of the start curves.
1725 :return: curve tags
1726 :rtype: list[int]
1727 """
1728 return self.innerStartCurves
1730 def getOuterEndCurves(self):
1731 """
1732 Returns the curve tags of the end curves.
1734 :return: curve tags
1735 :rtype: list[int]
1736 """
1737 return self.outerEndCurves
1740class circleWithFourCurves:
1741 def __init__(
1742 self,
1743 x,
1744 y,
1745 z,
1746 radius,
1747 upperRightTag=None,
1748 upperLeftTag=None,
1749 lowerLeftTag=None,
1750 lowerRightTag=None,
1751 leftPointTag=None,
1752 rightPointTag=None,
1753 upperPointTag=None,
1754 lowerPointTag=None,
1755 ):
1756 self.x = x
1757 self.y = y
1758 self.z = z
1759 self.r = radius
1760 if upperRightTag is None:
1761 dummyCircle = gmsh.model.occ.addCircle(self.x, self.y, self.z, self.r)
1762 self.leftPointTag = gmsh.model.occ.addPoint(self.x - self.r, self.y, self.z)
1763 self.rightPointTag = gmsh.model.occ.addPoint(
1764 self.x + self.r, self.y, self.z
1765 )
1766 self.upperPointTag = gmsh.model.occ.addPoint(
1767 self.x, self.y + self.r, self.z
1768 )
1769 self.lowerPointTag = gmsh.model.occ.addPoint(
1770 self.x, self.y - self.r, self.z
1771 )
1773 fragmentResults = gmsh.model.occ.fragment(
1774 [(1, dummyCircle)],
1775 [
1776 (0, self.leftPointTag),
1777 (0, self.rightPointTag),
1778 (0, self.upperPointTag),
1779 (0, self.lowerPointTag),
1780 ],
1781 )[0]
1782 linesDimTags = [dimTag for dimTag in fragmentResults if dimTag[0] == 1]
1784 self.upperRightTag = linesDimTags[0][1]
1785 self.upperLeftTag = linesDimTags[1][1]
1786 self.lowerLeftTag = linesDimTags[2][1]
1787 self.lowerRightTag = linesDimTags[3][1]
1788 else:
1789 self.upperRightTag = upperRightTag
1790 self.upperLeftTag = upperLeftTag
1791 self.lowerLeftTag = lowerLeftTag
1792 self.lowerRightTag = lowerRightTag
1794 self.leftPointTag = leftPointTag
1795 self.rightPointTag = rightPointTag
1796 self.upperPointTag = upperPointTag
1797 self.lowerPointTag = lowerPointTag
1800class outerAirSurface:
1801 def __init__(
1802 self,
1803 outerRadius,
1804 innerRadius,
1805 type="cylinder",
1806 divideIntoFourParts=False,
1807 divideTerminalPartIntoFourParts=False,
1808 ):
1809 self.surfaceTags = []
1811 self.divideIntoFourParts = divideIntoFourParts
1812 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts
1814 # for cylinder:
1815 self.shellTags = []
1817 # for cuboid:
1818 self.shellTagsPart1 = []
1819 self.shellTagsPart2 = []
1820 self.shellTagsPart3 = []
1821 self.shellTagsPart4 = []
1823 self.type = type
1824 self.outerRadius = outerRadius
1825 self.innerRadius = innerRadius
1827 def createFromScratch(self, z, shellTransformation=False, shellRadius=None):
1828 self.z = z
1830 if self.divideIntoFourParts:
1831 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)
1832 else:
1833 if self.divideTerminalPartIntoFourParts:
1834 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)
1835 innerCL = gmsh.model.occ.addCurveLoop(
1836 [
1837 self.innerCircle.upperRightTag,
1838 self.innerCircle.upperLeftTag,
1839 self.innerCircle.lowerLeftTag,
1840 self.innerCircle.lowerRightTag,
1841 ]
1842 )
1843 else:
1844 innerCL = gmsh.model.occ.addCircle(0, 0, z, self.innerRadius)
1845 innerCL = gmsh.model.occ.addCurveLoop([innerCL])
1847 if self.type == "cylinder" and self.divideIntoFourParts:
1848 outerCircle = circleWithFourCurves(0, 0, z, self.outerRadius)
1850 leftLineTag = gmsh.model.occ.addLine(
1851 outerCircle.leftPointTag, self.innerCircle.leftPointTag
1852 )
1853 rightLineTag = gmsh.model.occ.addLine(
1854 outerCircle.rightPointTag, self.innerCircle.rightPointTag
1855 )
1856 upperLineTag = gmsh.model.occ.addLine(
1857 outerCircle.upperPointTag, self.innerCircle.upperPointTag
1858 )
1859 lowerLineTag = gmsh.model.occ.addLine(
1860 outerCircle.lowerPointTag, self.innerCircle.lowerPointTag
1861 )
1863 # Create surfaces:
1864 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
1865 [
1866 outerCircle.upperRightTag,
1867 rightLineTag,
1868 self.innerCircle.upperRightTag,
1869 -upperLineTag,
1870 ]
1871 )
1872 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
1873 self.surfaceTags.append(self.upperRightTag)
1875 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
1876 [
1877 outerCircle.upperLeftTag,
1878 leftLineTag,
1879 self.innerCircle.upperLeftTag,
1880 -upperLineTag,
1881 ]
1882 )
1883 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
1884 self.surfaceTags.append(self.upperLeftTag)
1886 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
1887 [
1888 outerCircle.lowerLeftTag,
1889 leftLineTag,
1890 self.innerCircle.lowerLeftTag,
1891 -lowerLineTag,
1892 ]
1893 )
1894 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
1895 self.surfaceTags.append(self.lowerLeftTag)
1897 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
1898 [
1899 outerCircle.lowerRightTag,
1900 rightLineTag,
1901 self.innerCircle.lowerRightTag,
1902 -lowerLineTag,
1903 ]
1904 )
1905 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
1906 self.surfaceTags.append(self.lowerRightTag)
1908 if shellTransformation:
1909 shellOuterCircle = circleWithFourCurves(0, 0, z, shellRadius)
1910 shellLeftLineTag = gmsh.model.occ.addLine(
1911 shellOuterCircle.leftPointTag, outerCircle.leftPointTag
1912 )
1913 shellRightLineTag = gmsh.model.occ.addLine(
1914 shellOuterCircle.rightPointTag, outerCircle.rightPointTag
1915 )
1916 shellUpperLineTag = gmsh.model.occ.addLine(
1917 shellOuterCircle.upperPointTag, outerCircle.upperPointTag
1918 )
1919 shellLowerLineTag = gmsh.model.occ.addLine(
1920 shellOuterCircle.lowerPointTag, outerCircle.lowerPointTag
1921 )
1923 # Create surfaces:
1924 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
1925 [
1926 shellOuterCircle.upperRightTag,
1927 shellRightLineTag,
1928 outerCircle.upperRightTag,
1929 -shellUpperLineTag,
1930 ]
1931 )
1932 self.upperRightTag = gmsh.model.occ.addPlaneSurface(
1933 [upperRightCurveLoop]
1934 )
1935 self.shellTags.append(self.upperRightTag)
1937 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
1938 [
1939 shellOuterCircle.upperLeftTag,
1940 shellLeftLineTag,
1941 outerCircle.upperLeftTag,
1942 -shellUpperLineTag,
1943 ]
1944 )
1945 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
1946 self.shellTags.append(self.upperLeftTag)
1948 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
1949 [
1950 shellOuterCircle.lowerLeftTag,
1951 shellLeftLineTag,
1952 outerCircle.lowerLeftTag,
1953 -shellLowerLineTag,
1954 ]
1955 )
1956 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
1957 self.shellTags.append(self.lowerLeftTag)
1959 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
1960 [
1961 shellOuterCircle.lowerRightTag,
1962 shellRightLineTag,
1963 outerCircle.lowerRightTag,
1964 -shellLowerLineTag,
1965 ]
1966 )
1967 self.lowerRightTag = gmsh.model.occ.addPlaneSurface(
1968 [lowerRightCurveLoop]
1969 )
1970 self.shellTags.append(self.lowerRightTag)
1972 elif self.type == "cylinder" and not self.divideIntoFourParts:
1973 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.outerRadius)
1974 outerCL = gmsh.model.occ.addCurveLoop([outerCL])
1976 if shellTransformation:
1977 shellOuterCL = gmsh.model.occ.addCircle(0, 0, z, shellRadius)
1978 shellOuterCL = gmsh.model.occ.addCurveLoop([shellOuterCL])
1980 shellSurfaceTag = gmsh.model.occ.addPlaneSurface(
1981 [shellOuterCL, outerCL]
1982 )
1983 self.shellTags.append(shellSurfaceTag)
1985 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])
1986 self.surfaceTags.append(surfaceTag)
1988 elif self.type == "cuboid":
1989 # LL: lower left
1990 # LR: lower right
1991 # UR: upper right
1992 # UL: upper left
1993 airLLpointTag = point(-self.outerRadius, -self.outerRadius, z).tag
1994 airLRpointTag = point(self.outerRadius, -self.outerRadius, z).tag
1995 airURpointTag = point(self.outerRadius, self.outerRadius, z).tag
1996 airULpointTag = point(-self.outerRadius, self.outerRadius, z).tag
1998 # LH: lower horizontal
1999 # UH: upper horizontal
2000 # LV: left vertical
2001 # RV: right vertical
2002 airLHlineTag = gmsh.model.occ.addLine(airLLpointTag, airLRpointTag)
2003 airRVLineTag = gmsh.model.occ.addLine(airLRpointTag, airURpointTag)
2004 airUHLineTag = gmsh.model.occ.addLine(airURpointTag, airULpointTag)
2005 airLVLineTag = gmsh.model.occ.addLine(airULpointTag, airLLpointTag)
2007 outerCL = gmsh.model.occ.addCurveLoop(
2008 [airLHlineTag, airRVLineTag, airUHLineTag, airLVLineTag]
2009 )
2011 if self.divideIntoFourParts:
2012 innerCL = gmsh.model.occ.addCurveLoop(
2013 [
2014 self.innerCircle.upperRightTag,
2015 self.innerCircle.lowerRightTag,
2016 self.innerCircle.lowerLeftTag,
2017 self.innerCircle.upperLeftTag,
2018 ]
2019 )
2021 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])
2022 self.surfaceTags.append(surfaceTag)
2024 if shellTransformation:
2025 # LL: lower left
2026 # LR: lower right
2027 # UR: upper right
2028 # UL: upper left
2029 shellLLpointTag = point(
2030 -shellRadius,
2031 -shellRadius,
2032 z,
2033 ).tag
2034 shellLRpointTag = point(
2035 shellRadius,
2036 -shellRadius,
2037 z,
2038 ).tag
2039 shellURpointTag = point(
2040 shellRadius,
2041 shellRadius,
2042 z,
2043 ).tag
2044 shellULpointTag = point(
2045 -shellRadius,
2046 shellRadius,
2047 z,
2048 ).tag
2050 # LH: lower horizontal
2051 # UH: upper horizontal
2052 # LV: left vertical
2053 # RV: right vertical
2054 shellLHlineTag = gmsh.model.occ.addLine(
2055 shellLLpointTag, shellLRpointTag
2056 )
2057 shellRVLineTag = gmsh.model.occ.addLine(
2058 shellLRpointTag, shellURpointTag
2059 )
2060 shellUHLineTag = gmsh.model.occ.addLine(
2061 shellURpointTag, shellULpointTag
2062 )
2063 shellLVLineTag = gmsh.model.occ.addLine(
2064 shellULpointTag, shellLLpointTag
2065 )
2067 shellLowerLeftLineTag = gmsh.model.occ.addLine(
2068 shellLLpointTag, airLLpointTag
2069 )
2070 shellLowerRightLineTag = gmsh.model.occ.addLine(
2071 shellLRpointTag, airLRpointTag
2072 )
2073 shellUpperLeftLineTag = gmsh.model.occ.addLine(
2074 shellULpointTag, airULpointTag
2075 )
2076 shellUpperRightLineTag = gmsh.model.occ.addLine(
2077 shellURpointTag, airURpointTag
2078 )
2080 # Shell lower surface:
2081 shellLowerPSTag = gmsh.model.occ.addCurveLoop(
2082 [
2083 shellLowerLeftLineTag,
2084 airLHlineTag,
2085 shellLowerRightLineTag,
2086 shellLHlineTag,
2087 ]
2088 )
2089 shellLowerPSTag = gmsh.model.occ.addPlaneSurface([shellLowerPSTag])
2090 self.shellTagsPart1.append(shellLowerPSTag)
2092 # Shell right surface:
2093 shellRightPSTag = gmsh.model.occ.addCurveLoop(
2094 [
2095 shellLowerRightLineTag,
2096 airRVLineTag,
2097 shellUpperRightLineTag,
2098 shellRVLineTag,
2099 ]
2100 )
2101 shellRightPSTag = gmsh.model.occ.addPlaneSurface([shellRightPSTag])
2102 self.shellTagsPart2.append(shellRightPSTag)
2104 # Shell upper surface:
2105 shellUpperPSTag = gmsh.model.occ.addCurveLoop(
2106 [
2107 shellUpperLeftLineTag,
2108 airUHLineTag,
2109 shellUpperRightLineTag,
2110 shellUHLineTag,
2111 ]
2112 )
2113 shellUpperPSTag = gmsh.model.occ.addPlaneSurface([shellUpperPSTag])
2114 self.shellTagsPart3.append(shellUpperPSTag)
2116 # Shell left surface:
2117 shellLeftPSTag = gmsh.model.occ.addCurveLoop(
2118 [
2119 shellLowerLeftLineTag,
2120 airLVLineTag,
2121 shellUpperLeftLineTag,
2122 shellLVLineTag,
2123 ]
2124 )
2125 shellLeftPSTag = gmsh.model.occ.addPlaneSurface([shellLeftPSTag])
2126 self.shellTagsPart4.append(shellLeftPSTag)
2128 def setPrecreatedSurfaceTags(
2129 self,
2130 surfaceTags,
2131 cylinderShellTags=None,
2132 cuboidShellTags1=None,
2133 cuboidShellTags2=None,
2134 cuboidShellTags3=None,
2135 cuboidShellTags4=None,
2136 ):
2137 if not isinstance(surfaceTags, list):
2138 raise TypeError("surfaceTags must be a list.")
2140 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2]
2141 self.surfaceTags.extend(surfaceTags)
2143 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:
2144 # Create innerCircle object from the tags:
2145 curves = gmsh.model.getBoundary(
2146 [(2, tag) for tag in surfaceTags], oriented=False
2147 )
2148 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)
2149 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]
2150 innerCurveTags.sort()
2151 upperRightCurve = innerCurveTags[0]
2152 upperLeftCurve = innerCurveTags[1]
2153 lowerLeftCurve = innerCurveTags[2]
2154 lowerRightCurve = innerCurveTags[3]
2156 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2157 pointTags = [dimTag[1] for dimTag in points]
2158 pointTags.sort()
2159 upperPointTag = pointTags[0]
2160 leftPointTag = pointTags[1]
2162 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2163 pointTags = [dimTag[1] for dimTag in points]
2164 pointTags.sort()
2165 rightPointTag = pointTags[0]
2166 lowerPointTag = pointTags[1]
2168 self.innerCircle = circleWithFourCurves(
2169 0,
2170 0,
2171 self.z,
2172 self.outerRadius,
2173 upperRightTag=upperRightCurve,
2174 upperLeftTag=upperLeftCurve,
2175 lowerLeftTag=lowerLeftCurve,
2176 lowerRightTag=lowerRightCurve,
2177 leftPointTag=leftPointTag,
2178 rightPointTag=rightPointTag,
2179 upperPointTag=upperPointTag,
2180 lowerPointTag=lowerPointTag,
2181 )
2183 if cylinderShellTags is not None:
2184 self.shellTags.extend(cylinderShellTags)
2186 if cuboidShellTags1 is not None:
2187 self.shellTagsPart1.extend(cuboidShellTags1)
2188 self.shellTagsPart2.extend(cuboidShellTags2)
2189 self.shellTagsPart3.extend(cuboidShellTags3)
2190 self.shellTagsPart4.extend(cuboidShellTags4)
2192 def getInnerCL(self):
2193 # checked!
2194 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])
2195 # curves = list(curves)
2196 # curves = list(curves[1])
2198 # innerCL = gmsh.model.occ.addCurveLoop(curves)
2200 gmsh.model.occ.synchronize() # don't delete this line
2201 curves = gmsh.model.getBoundary(
2202 [(2, tag) for tag in self.surfaceTags], oriented=False
2203 )
2204 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)
2205 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]
2207 innerCL = gmsh.model.occ.addCurveLoop(innerCurveTags)
2208 return innerCL
2211class outerTerminalSurface:
2212 def __init__(
2213 self,
2214 outerRadius,
2215 tubeThickness,
2216 divideIntoFourParts=False,
2217 ):
2218 self.tubeSurfaceTags = []
2219 self.nontubeSurfaceTags = []
2221 self.divideIntoFourParts = divideIntoFourParts
2223 self.outerRadius = outerRadius
2224 self.tubeThickness = tubeThickness
2226 def createNontubePartWithMiddleCircleAndWinding(
2227 self, middleCircle: circleWithFourCurves, winding: spiralSurface
2228 ):
2229 leftLineTag = gmsh.model.occ.addLine(
2230 middleCircle.leftPointTag, winding.getOuterLeftPointTag()
2231 )
2232 rightLineTag = gmsh.model.occ.addLine(
2233 middleCircle.rightPointTag, winding.getOuterRightPointTag()
2234 )
2235 upperLineTag = gmsh.model.occ.addLine(
2236 middleCircle.upperPointTag, winding.getOuterUpperPointTag()
2237 )
2238 lowerLineTag = gmsh.model.occ.addLine(
2239 middleCircle.lowerPointTag, winding.getOuterLowerPointTag()
2240 )
2242 # Create surfaces for the nontube part:
2243 if winding.direction is direction.ccw:
2244 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2245 winding.getOuterUpperRightCurves()
2246 # + winding.getOuterEndCurves()
2247 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]
2248 )
2249 else:
2250 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2251 winding.getOuterUpperRightCurves()
2252 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]
2253 )
2254 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2255 self.nontubeSurfaceTags.append(self.upperRightTag)
2257 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2258 winding.getOuterUpperLeftCurves()
2259 + [leftLineTag, middleCircle.upperLeftTag, upperLineTag]
2260 )
2261 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2262 self.nontubeSurfaceTags.append(self.upperLeftTag)
2264 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2265 winding.getOuterLowerLeftCurves()
2266 + [leftLineTag, middleCircle.lowerLeftTag, lowerLineTag]
2267 )
2268 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2269 self.nontubeSurfaceTags.append(self.lowerLeftTag)
2271 if winding.direction is direction.ccw:
2272 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2273 winding.getOuterLowerRightCurves()
2274 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]
2275 )
2276 else:
2277 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2278 winding.getOuterLowerRightCurves()
2279 # + winding.getOuterEndCurves()
2280 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]
2281 )
2282 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2283 self.nontubeSurfaceTags.append(self.lowerRightTag)
2285 def createWithOuterAirAndWinding(
2286 self, outerAir: outerAirSurface, winding: spiralSurface, pancakeIndex
2287 ):
2288 # Tube part:
2289 z = outerAir.z
2291 if self.divideIntoFourParts:
2292 outerCircle = outerAir.innerCircle
2293 middleCircle = circleWithFourCurves(
2294 0, 0, z, self.outerRadius - self.tubeThickness
2295 )
2297 leftLineTag = gmsh.model.occ.addLine(
2298 outerCircle.leftPointTag, middleCircle.leftPointTag
2299 )
2300 rightLineTag = gmsh.model.occ.addLine(
2301 outerCircle.rightPointTag, middleCircle.rightPointTag
2302 )
2303 upperLineTag = gmsh.model.occ.addLine(
2304 outerCircle.upperPointTag, middleCircle.upperPointTag
2305 )
2306 lowerLineTag = gmsh.model.occ.addLine(
2307 outerCircle.lowerPointTag, middleCircle.lowerPointTag
2308 )
2310 # Create surfaces for the tube part:
2311 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2312 [
2313 outerCircle.upperRightTag,
2314 rightLineTag,
2315 middleCircle.upperRightTag,
2316 -upperLineTag,
2317 ]
2318 )
2319 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2320 self.tubeSurfaceTags.append(self.upperRightTag)
2322 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2323 [
2324 outerCircle.upperLeftTag,
2325 leftLineTag,
2326 middleCircle.upperLeftTag,
2327 -upperLineTag,
2328 ]
2329 )
2330 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2331 self.tubeSurfaceTags.append(self.upperLeftTag)
2333 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2334 [
2335 outerCircle.lowerLeftTag,
2336 leftLineTag,
2337 middleCircle.lowerLeftTag,
2338 -lowerLineTag,
2339 ]
2340 )
2341 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2342 self.tubeSurfaceTags.append(self.lowerLeftTag)
2344 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2345 [
2346 outerCircle.lowerRightTag,
2347 rightLineTag,
2348 middleCircle.lowerRightTag,
2349 -lowerLineTag,
2350 ]
2351 )
2352 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2353 self.tubeSurfaceTags.append(self.lowerRightTag)
2355 else:
2356 outerCL = outerAir.getInnerCL()
2358 middleCL = gmsh.model.occ.addCircle(
2359 0, 0, z, self.outerRadius - self.tubeThickness
2360 )
2361 middleCL = gmsh.model.occ.addCurveLoop([middleCL])
2363 tubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
2364 self.tubeSurfaceTags.append(tubeSurface)
2366 # Nontube part:
2367 if self.divideIntoFourParts:
2368 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2370 else:
2371 # Inner and outer curve loops. Sometimes the opposite curve loops are used
2372 # because the other one comes from the self.contactTerminalSurface. To create
2373 # a valid surface, the topology of the curve loops should be consistent. See the
2374 # note in the spiralSurface class.
2375 if pancakeIndex % 2 == 0:
2376 innerCL = winding.outerCurveLoopTag
2377 elif pancakeIndex % 2 == 1:
2378 innerCL = winding.outerOppositeCurveLoopTag
2380 # potential bug (curve order might be wrong)
2381 if self.divideIntoFourParts:
2382 middleCL = gmsh.model.occ.addCurveLoop(
2383 [
2384 middleCircle.upperRightTag,
2385 middleCircle.upperLeftTag,
2386 middleCircle.lowerLeftTag,
2387 middleCircle.lowerRightTag,
2388 ]
2389 )
2391 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2392 self.nontubeSurfaceTags.append(nontubeSurface)
2394 def createWithWindingAndTubeTags(
2395 self, winding: spiralSurface, tubeTags, pancakeIndex
2396 ):
2397 if not isinstance(tubeTags, list):
2398 raise TypeError("tubeTags must be a list.")
2400 self.tubeSurfaceTags.extend(tubeTags)
2402 middleCurves = gmsh.model.getBoundary(
2403 [(2, tag) for tag in tubeTags], oriented=False
2404 )
2405 middleCurveDimTags = findOuterOnes(middleCurves, findInnerOnes=True)
2406 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]
2408 if self.divideIntoFourParts:
2409 # Create middleCircle object from the tags:
2410 middleCurveTags.sort()
2411 upperRightCurve = middleCurveTags[0]
2412 upperLeftCurve = middleCurveTags[1]
2413 lowerLeftCurve = middleCurveTags[2]
2414 lowerRightCurve = middleCurveTags[3]
2416 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2417 pointTags = [dimTag[1] for dimTag in points]
2418 pointTags.sort()
2419 upperPointTag = pointTags[0]
2420 leftPointTag = pointTags[1]
2422 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2423 pointTags = [dimTag[1] for dimTag in points]
2424 pointTags.sort()
2425 rightPointTag = pointTags[0]
2426 lowerPointTag = pointTags[1]
2428 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]
2429 middleCircle = circleWithFourCurves(
2430 0,
2431 0,
2432 z,
2433 self.outerRadius - self.tubeThickness,
2434 upperRightTag=upperRightCurve,
2435 upperLeftTag=upperLeftCurve,
2436 lowerLeftTag=lowerLeftCurve,
2437 lowerRightTag=lowerRightCurve,
2438 leftPointTag=leftPointTag,
2439 rightPointTag=rightPointTag,
2440 upperPointTag=upperPointTag,
2441 lowerPointTag=lowerPointTag,
2442 )
2444 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2445 else:
2446 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)
2448 if pancakeIndex % 2 == 0:
2449 innerCL = winding.outerCurveLoopTag
2450 elif pancakeIndex % 2 == 1:
2451 innerCL = winding.outerOppositeCurveLoopTag
2453 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2454 self.nontubeSurfaceTags.append(nontubeSurface)
2457class innerAirSurface:
2458 def __init__(
2459 self, radius, divideIntoFourParts=False, divideTerminalPartIntoFourParts=False
2460 ):
2461 self.surfaceTags = []
2463 self.divideIntoFourParts = divideIntoFourParts
2464 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts
2466 self.radius = radius
2468 def createFromScratch(self, z):
2469 self.z = z
2470 if self.divideIntoFourParts:
2471 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)
2473 originTag = point(0, 0, z).tag
2475 leftLineTag = gmsh.model.occ.addLine(
2476 self.outerCircle.leftPointTag, originTag
2477 )
2478 rightLineTag = gmsh.model.occ.addLine(
2479 self.outerCircle.rightPointTag, originTag
2480 )
2481 upperLineTag = gmsh.model.occ.addLine(
2482 self.outerCircle.upperPointTag, originTag
2483 )
2484 lowerLineTag = gmsh.model.occ.addLine(
2485 self.outerCircle.lowerPointTag, originTag
2486 )
2488 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2489 [self.outerCircle.upperRightTag, rightLineTag, -upperLineTag]
2490 )
2491 upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2492 self.surfaceTags.append(upperRightTag)
2494 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2495 [self.outerCircle.upperLeftTag, leftLineTag, -upperLineTag]
2496 )
2497 upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2498 self.surfaceTags.append(upperLeftTag)
2500 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2501 [self.outerCircle.lowerLeftTag, leftLineTag, -lowerLineTag]
2502 )
2503 lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2504 self.surfaceTags.append(lowerLeftTag)
2506 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2507 [self.outerCircle.lowerRightTag, rightLineTag, -lowerLineTag]
2508 )
2509 lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2510 self.surfaceTags.append(lowerRightTag)
2512 else:
2513 if self.divideTerminalPartIntoFourParts:
2514 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)
2515 outerCL = gmsh.model.occ.addCurveLoop(
2516 [
2517 self.outerCircle.upperRightTag,
2518 self.outerCircle.upperLeftTag,
2519 self.outerCircle.lowerLeftTag,
2520 self.outerCircle.lowerRightTag,
2521 ]
2522 )
2523 else:
2524 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.radius)
2525 outerCL = gmsh.model.occ.addCurveLoop([outerCL])
2527 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL])
2528 self.surfaceTags.append(surfaceTag)
2530 def setPrecreatedSurfaceTags(self, surfaceTags):
2531 if not isinstance(surfaceTags, list):
2532 raise TypeError("surfaceTags must be a list.")
2534 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2] # potential bug
2535 self.surfaceTags = []
2536 self.surfaceTags.extend(surfaceTags)
2538 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:
2539 # Create outerCirle object from the tags:
2540 curves = gmsh.model.getBoundary(
2541 [(2, tag) for tag in surfaceTags], oriented=False
2542 )
2543 outerCurveDimTags = findOuterOnes(curves)
2544 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]
2545 outerCurveTags.sort()
2546 upperRightCurve = outerCurveTags[0]
2547 upperLeftCurve = outerCurveTags[1]
2548 lowerLeftCurve = outerCurveTags[2]
2549 lowerRightCurve = outerCurveTags[3]
2551 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2552 pointTags = [dimTag[1] for dimTag in points]
2553 pointTags.sort()
2554 upperPointTag = pointTags[0]
2555 leftPointTag = pointTags[1]
2557 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2558 pointTags = [dimTag[1] for dimTag in points]
2559 pointTags.sort()
2560 rightPointTag = pointTags[0]
2561 lowerPointTag = pointTags[1]
2563 self.outerCircle = circleWithFourCurves(
2564 0,
2565 0,
2566 self.z,
2567 self.radius,
2568 upperRightTag=upperRightCurve,
2569 upperLeftTag=upperLeftCurve,
2570 lowerLeftTag=lowerLeftCurve,
2571 lowerRightTag=lowerRightCurve,
2572 leftPointTag=leftPointTag,
2573 rightPointTag=rightPointTag,
2574 upperPointTag=upperPointTag,
2575 lowerPointTag=lowerPointTag,
2576 )
2578 def getOuterCL(self):
2579 # checked!
2580 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])
2581 # curves = list(curves)
2582 # curves = [int(curves[0])]
2584 # outerCL = gmsh.model.occ.addCurveLoop([curves[0]])
2586 # return outerCL
2588 curves = gmsh.model.getBoundary(
2589 [(2, tag) for tag in self.surfaceTags], oriented=False
2590 )
2591 outerCurveDimTags = findOuterOnes(curves)
2592 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]
2594 outerCL = gmsh.model.occ.addCurveLoop(outerCurveTags)
2595 return outerCL
2598class innerTerminalSurface:
2599 def __init__(self, innerRadius, tubeThickness, divideIntoFourParts=False):
2600 self.tubeSurfaceTags = []
2601 self.nontubeSurfaceTags = []
2603 self.divideIntoFourParts = divideIntoFourParts
2605 self.innerRadius = innerRadius
2606 self.tubeThickness = tubeThickness
2608 def createNontubePartWithMiddleCircleAndWinding(
2609 self, middleCircle: circleWithFourCurves, winding: spiralSurface
2610 ):
2611 leftLineTag = gmsh.model.occ.addLine(
2612 winding.getInnerLeftPointTag(), middleCircle.leftPointTag
2613 )
2614 rightLineTag = gmsh.model.occ.addLine(
2615 winding.getInnerRightPointTag(), middleCircle.rightPointTag
2616 )
2617 upperLineTag = gmsh.model.occ.addLine(
2618 winding.getInnerUpperPointTag(), middleCircle.upperPointTag
2619 )
2620 lowerLineTag = gmsh.model.occ.addLine(
2621 winding.getInnerLowerPointTag(), middleCircle.lowerPointTag
2622 )
2624 # Create surfaces for the nontube part:
2625 if winding.direction is direction.ccw:
2626 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2627 winding.getInnerUpperRightCurves()
2628 + [
2629 upperLineTag,
2630 middleCircle.upperRightTag,
2631 rightLineTag,
2632 ]
2633 )
2634 else:
2635 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2636 winding.getInnerUpperRightCurves()
2637 + [
2638 upperLineTag,
2639 middleCircle.upperRightTag,
2640 rightLineTag,
2641 ]
2642 # + winding.getInnerStartCurves()
2643 )
2644 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2645 self.nontubeSurfaceTags.append(self.upperRightTag)
2647 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2648 winding.getInnerUpperLeftCurves()
2649 + [
2650 upperLineTag,
2651 middleCircle.upperLeftTag,
2652 leftLineTag,
2653 ]
2654 )
2655 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2656 self.nontubeSurfaceTags.append(self.upperLeftTag)
2658 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2659 winding.getInnerLowerLeftCurves()
2660 + [
2661 lowerLineTag,
2662 middleCircle.lowerLeftTag,
2663 leftLineTag,
2664 ]
2665 )
2666 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2667 self.nontubeSurfaceTags.append(self.lowerLeftTag)
2669 if winding.direction is direction.ccw:
2670 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2671 winding.getInnerLowerRightCurves()
2672 + [
2673 lowerLineTag,
2674 middleCircle.lowerRightTag,
2675 rightLineTag,
2676 ]
2677 # + winding.getInnerStartCurves()
2678 )
2679 else:
2680 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2681 winding.getInnerLowerRightCurves()
2682 + [
2683 lowerLineTag,
2684 middleCircle.lowerRightTag,
2685 rightLineTag,
2686 ]
2687 )
2688 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2689 self.nontubeSurfaceTags.append(self.lowerRightTag)
2691 def createWithInnerAirAndWinding(
2692 self, innerAir: innerAirSurface, winding: spiralSurface, pancakeIndex
2693 ):
2694 z = innerAir.z
2696 # Tube part:
2697 if self.divideIntoFourParts:
2698 innerCircle = innerAir.outerCircle
2699 middleCircle = circleWithFourCurves(
2700 0, 0, z, self.innerRadius + self.tubeThickness
2701 )
2703 leftLineTag = gmsh.model.occ.addLine(
2704 middleCircle.leftPointTag, innerCircle.leftPointTag
2705 )
2706 rightLineTag = gmsh.model.occ.addLine(
2707 middleCircle.rightPointTag, innerCircle.rightPointTag
2708 )
2709 upperLineTag = gmsh.model.occ.addLine(
2710 middleCircle.upperPointTag, innerCircle.upperPointTag
2711 )
2712 lowerLineTag = gmsh.model.occ.addLine(
2713 middleCircle.lowerPointTag, innerCircle.lowerPointTag
2714 )
2716 # Create surfaces for the tube part:
2717 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(
2718 [
2719 middleCircle.upperRightTag,
2720 rightLineTag,
2721 innerCircle.upperRightTag,
2722 -upperLineTag,
2723 ]
2724 )
2725 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])
2726 self.tubeSurfaceTags.append(self.upperRightTag)
2728 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2729 [
2730 middleCircle.upperLeftTag,
2731 leftLineTag,
2732 innerCircle.upperLeftTag,
2733 -upperLineTag,
2734 ]
2735 )
2736 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])
2737 self.tubeSurfaceTags.append(self.upperLeftTag)
2739 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(
2740 [
2741 middleCircle.lowerLeftTag,
2742 leftLineTag,
2743 innerCircle.lowerLeftTag,
2744 -lowerLineTag,
2745 ]
2746 )
2747 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])
2748 self.tubeSurfaceTags.append(self.lowerLeftTag)
2750 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(
2751 [
2752 middleCircle.lowerRightTag,
2753 rightLineTag,
2754 innerCircle.lowerRightTag,
2755 -lowerLineTag,
2756 ]
2757 )
2758 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])
2759 self.tubeSurfaceTags.append(self.lowerRightTag)
2761 else:
2762 innerCL = innerAir.getOuterCL()
2764 middleCL = gmsh.model.occ.addCircle(
2765 0, 0, z, self.innerRadius + self.tubeThickness
2766 )
2767 middleCL = gmsh.model.occ.addCurveLoop([middleCL])
2769 tubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])
2770 self.tubeSurfaceTags.append(tubeSurface)
2772 # Nontube part:
2773 if self.divideIntoFourParts:
2774 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2776 else:
2777 # Inner and outer curve loops. Sometimes the opposite curve loops are used
2778 # because the other one comes from the self.contactTerminalSurface. To create
2779 # a valid surface, the topology of the curve loops should be consistent. See the
2780 # note in the spiralSurface class.
2781 if pancakeIndex == 0:
2782 outerCL = winding.innerCurveLoopTag
2784 elif pancakeIndex % 2 == 0:
2785 outerCL = winding.innerOppositeCurveLoopTag
2787 elif pancakeIndex % 2 == 1:
2788 outerCL = winding.innerCurveLoopTag
2790 # potential bug (curve order might be wrong)
2791 if self.divideIntoFourParts:
2792 middleCL = gmsh.model.occ.addCurveLoop(
2793 [
2794 middleCircle.upperRightTag,
2795 middleCircle.upperLeftTag,
2796 middleCircle.lowerLeftTag,
2797 middleCircle.lowerRightTag,
2798 ]
2799 )
2801 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
2802 self.nontubeSurfaceTags.append(nontubeSurface)
2804 def createWithWindingAndTubeTags(
2805 self, winding: spiralSurface, tubeTags, pancakeIndex
2806 ):
2807 if not isinstance(tubeTags, list):
2808 raise TypeError("tubeTags must be a list.")
2810 self.tubeSurfaceTags.extend(tubeTags)
2812 middleCurves = gmsh.model.getBoundary(
2813 [(2, tag) for tag in tubeTags], oriented=False
2814 )
2815 middleCurveDimTags = findOuterOnes(middleCurves)
2816 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]
2818 if self.divideIntoFourParts:
2819 # Create middleCircle object from the tags:
2820 middleCurveTags.sort()
2821 upperRightCurve = middleCurveTags[0]
2822 upperLeftCurve = middleCurveTags[1]
2823 lowerLeftCurve = middleCurveTags[2]
2824 lowerRightCurve = middleCurveTags[3]
2826 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)
2827 pointTags = [dimTag[1] for dimTag in points]
2828 pointTags.sort()
2829 upperPointTag = pointTags[0]
2830 leftPointTag = pointTags[1]
2832 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)
2833 pointTags = [dimTag[1] for dimTag in points]
2834 pointTags.sort()
2835 rightPointTag = pointTags[0]
2836 lowerPointTag = pointTags[1]
2838 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]
2839 middleCircle = circleWithFourCurves(
2840 0,
2841 0,
2842 z,
2843 self.innerRadius + self.tubeThickness,
2844 upperRightTag=upperRightCurve,
2845 upperLeftTag=upperLeftCurve,
2846 lowerLeftTag=lowerLeftCurve,
2847 lowerRightTag=lowerRightCurve,
2848 leftPointTag=leftPointTag,
2849 rightPointTag=rightPointTag,
2850 upperPointTag=upperPointTag,
2851 lowerPointTag=lowerPointTag,
2852 )
2854 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)
2855 else:
2856 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)
2858 if pancakeIndex == 0:
2859 outerCL = winding.innerCurveLoopTag
2861 elif pancakeIndex % 2 == 0:
2862 outerCL = winding.innerOppositeCurveLoopTag
2864 elif pancakeIndex % 2 == 1:
2865 outerCL = winding.innerCurveLoopTag
2867 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])
2868 self.nontubeSurfaceTags.append(nontubeSurface)
2871class pancakeCoilsWithAir:
2872 """
2873 A class to create Pancake3D coil. With this class, any number of pancake coils stack
2874 can be created in GMSH. Moreover, the class also creates some parts of the air
2875 volume as well. It creates the inner cylinder air volume and the outer tube air
2876 volume. However, the air between the pancake coils is not created. It is created in
2877 the gapAir class.
2879 self.fundamentalSurfaces are the surfaces at the bottom of each pancake coil. They
2880 are created one by one with self.generateFundamentalSurfaces() method. The first
2881 created self.fundamentalSurfaces are the first pancake coil's bottom surfaces. Those
2882 surfaces include the outer air shell surface, outer air tube surface, outer
2883 terminal's outer tube part, outer terminal's touching part, winding surfaces,
2884 contact layer surfaces, inner terminal's touching part, inner terminal's inner tube
2885 part, and inner air disc. Terminals are divided into two because they can only be
2886 connected with the perfect tubes since each pancake coil is rotated in a different
2887 direction.
2889 For the first pancake coil, self.fundamentalSurfaces are extruded downwards to
2890 connect the terminal to the end of the geometry at the bottom.
2892 Then self.fundamentalSurfaces are extruded upwards to create the first pancake coil
2893 with self.extrudeWindingPart method. The method returns the extrusion's top
2894 surfaces, which are saved in the topSurfaces variable.
2896 Then those topSurfaces variable is given to another method, self.extrudeGapPart, and
2897 they are further extruded upwards up to the bottom of the next pancake coil.
2898 However, only the air tube, air cylinder, and connection terminal (the perfect inner
2899 terminal tube or outer terminal tube) are extruded. Otherwise, conformality would be
2900 impossible. The gaps are filled with air in gapAir class with fragment operation
2901 later. Then the top surfaces are returned by the method and saved in
2902 self.contactSurfaces variable.
2904 Then using the self.contactSurfaces, self.generateFundamentalSurfaces method creates
2905 the new fundamental surfaces. All the surfaces from the self.contactSurfaces are
2906 used in the new self.fundamentalSurfaces variable to avoid surface duplication.
2908 The logic goes until the last topSurfaces are extruded upwards to connect the last
2909 terminal to the top of the geometry.
2911 Every pancake coil's rotation direction is different each time. Otherwise, their
2912 magnetic fields would neutralize each other.
2914 The first and second pancake coils are connected with the inner terminal. Then the
2915 second and the third pancake coils are connected with the outer terminal. And so on.
2917 :param geometryData: geometry information
2918 """
2920 def __init__(self, geometryData, meshData) -> None:
2921 logger.info("Generating pancake coils has been started.")
2922 start_time = timeit.default_timer()
2924 # Data:
2925 self.geo = geometryData
2926 self.mesh = meshData
2928 # ==============================================================================
2929 # CREATING VOLUME STORAGES STARTS ==============================================
2930 # ==============================================================================
2931 # Air shell (they will be empty if shellTransformation == False):
2932 # For cylinder type:
2933 self.airShellVolume = dimTags(name=self.geo.ai.shellVolumeName, save=True)
2935 # For cuboid type:
2936 self.airShellVolumePart1 = dimTags(
2937 name=self.geo.ai.shellVolumeName + "-Part1", save=True
2938 )
2939 self.airShellVolumePart2 = dimTags(
2940 name=self.geo.ai.shellVolumeName + "-Part2", save=True
2941 )
2942 self.airShellVolumePart3 = dimTags(
2943 name=self.geo.ai.shellVolumeName + "-Part3", save=True
2944 )
2945 self.airShellVolumePart4 = dimTags(
2946 name=self.geo.ai.shellVolumeName + "-Part4", save=True
2947 )
2949 # Outer air tube volume (actually it is not a tube if the air type is cuboid):
2950 self.outerAirTubeVolume = dimTags(
2951 name=self.geo.ai.name + "-OuterTube", save=True, parentName=self.geo.ai.name
2952 )
2954 # Outer terminal's outer tube part:
2955 self.outerTerminalTubeVolume = dimTags(
2956 name=self.geo.ti.o.name + "-Tube", save=True, parentName=self.geo.ti.o.name
2957 )
2959 # Outer terminal's volume that touches the winding:
2960 self.outerTerminalTouchingVolume = dimTags(
2961 name=self.geo.ti.o.name + "-Touching",
2962 save=True,
2963 parentName=self.geo.ti.o.name,
2964 )
2966 # Inner terminal's volume that touches the winding:
2967 self.innerTerminalTouchingVolume = dimTags(
2968 name=self.geo.ti.i.name + "-Touching",
2969 save=True,
2970 parentName=self.geo.ti.i.name,
2971 )
2973 # Inner terminal's inner tube part:
2974 self.innerTerminalTubeVolume = dimTags(
2975 name=self.geo.ti.i.name + "-Tube", save=True, parentName=self.geo.ti.i.name
2976 )
2978 # Transition layers:
2979 self.innerTransitionNotchVolume = dimTags(
2980 name="innerTransitionNotch",
2981 save=True,
2982 )
2983 self.outerTransitionNotchVolume = dimTags(
2984 name="outerTransitionNotch",
2985 save=True,
2986 )
2988 # Inner air cylinder volume:
2989 self.centerAirCylinderVolume = dimTags(
2990 name=self.geo.ai.name + "-InnerCylinder",
2991 save=True,
2992 parentName=self.geo.ai.name,
2993 )
2995 # Top and bottom parts of the air volume:
2996 self.topAirPancakeWindingExtursionVolume = dimTags(
2997 name=self.geo.ai.name + "-TopPancakeWindingExtursion",
2998 save=True,
2999 parentName=self.geo.ai.name,
3000 )
3001 self.topAirPancakeContactLayerExtursionVolume = dimTags(
3002 name=self.geo.ai.name + "-TopPancakeContactLayerExtursion",
3003 save=True,
3004 parentName=self.geo.ai.name,
3005 )
3006 self.topAirTerminalsExtrusionVolume = dimTags(
3007 name=self.geo.ai.name + "-TopTerminalsExtrusion",
3008 save=True,
3009 parentName=self.geo.ai.name,
3010 )
3011 self.topAirTubeTerminalsExtrusionVolume = dimTags(
3012 name=self.geo.ai.name + "-TopTubeTerminalsExtrusion",
3013 save=True,
3014 parentName=self.geo.ai.name,
3015 )
3017 self.bottomAirPancakeWindingExtursionVolume = dimTags(
3018 name=self.geo.ai.name + "-BottomPancakeWindingExtursion",
3019 save=True,
3020 parentName=self.geo.ai.name,
3021 )
3022 self.bottomAirPancakeContactLayerExtursionVolume = dimTags(
3023 name=self.geo.ai.name + "-BottomPancakeContactLayerExtursion",
3024 save=True,
3025 parentName=self.geo.ai.name,
3026 )
3027 self.bottomAirTerminalsExtrusionVolume = dimTags(
3028 name=self.geo.ai.name + "-BottomTerminalsExtrusion",
3029 save=True,
3030 parentName=self.geo.ai.name,
3031 )
3032 self.bottomAirTubeTerminalsExtrusionVolume = dimTags(
3033 name=self.geo.ai.name + "-BottomTubeTerminalsExtrusion",
3034 save=True,
3035 parentName=self.geo.ai.name,
3036 )
3038 # Gap air:
3039 self.gapAirVolume = dimTags(
3040 name=self.geo.ai.name + "-Gap", save=True, parentName=self.geo.ai.name
3041 )
3043 # Create additional/optional volume storages (they might be used in the meshing
3044 # process):
3045 self.firstTerminalVolume = dimTags(name=self.geo.ti.firstName, save=True)
3046 self.lastTerminalVolume = dimTags(name=self.geo.ti.lastName, save=True)
3048 # ==============================================================================
3049 # CREATING VOLUME STORAGES ENDS ================================================
3050 # ==============================================================================
3052 # self.fundamentalSurfaces is a dictionary of surface dimTags tuples. The keys
3053 # are the dimTags objects of the corresponding volumes. The values are the
3054 # dimTags tuples of the surfaces that are used to extrude the volumes. It is
3055 # created in self.generateFundamentalSurfaces method.
3056 self.fundamentalSurfaces = {}
3058 # self.pancakeIndex stores the index of the current pancake coil.
3059 self.pancakeIndex = 0
3061 # self.contactSurfaces is a dictionary of surface dimTags tuples. The keys are
3062 # the dimTags objects of the corresponding volumes. The values are the dimTags
3063 # tuples of the surfaces that are obtained from the previous extrusion and used
3064 # for the next extrusion. The same surface is used for the next extrusion to
3065 # avoid surface duplication. It is created in self.extrudeGapPart and
3066 # self.extrudeWindingPart methods.
3067 self.contactSurfaces = {}
3069 # They will be lists of dimTags objects:
3070 self.individualWinding = []
3071 self.individualContactLayer = []
3073 self.gapAirSurfacesDimTags = []
3075 for i in range(self.geo.N):
3076 # Itterate over the number of pancake coils:
3077 self.individualWinding.append(
3078 dimTags(
3079 name=self.geo.wi.name + str(self.pancakeIndex + 1),
3080 save=True,
3081 parentName=self.geo.wi.name,
3082 )
3083 )
3084 self.individualContactLayer.append(
3085 dimTags(
3086 name=self.geo.ii.name + str(self.pancakeIndex + 1),
3087 save=True,
3088 parentName=self.geo.ii.name,
3089 )
3090 )
3092 # Generate the fundamental surfaces:
3093 self.fundamentalSurfaces = self.generateFundamentalSurfaces()
3095 # Create gap air or collect the gap air surfaces:
3096 if i != 0:
3097 bottomSurfacesDimTags = []
3098 for key, value in topSurfaces.items():
3099 if (
3100 key is self.individualWinding[self.pancakeIndex - 1]
3101 or key is self.individualContactLayer[self.pancakeIndex - 1]
3102 or key is self.outerTerminalTouchingVolume
3103 or key is self.innerTerminalTouchingVolume
3104 or key is self.innerTransitionNotchVolume
3105 or key is self.outerTransitionNotchVolume
3106 ):
3107 bottomSurfacesDimTags.extend(value)
3109 topSurfacesDimTags = []
3110 for key, value in self.fundamentalSurfaces.items():
3111 if (
3112 key is self.individualWinding[self.pancakeIndex]
3113 or key is self.individualContactLayer[self.pancakeIndex]
3114 or key is self.outerTerminalTouchingVolume
3115 or key is self.innerTerminalTouchingVolume
3116 or key is self.innerTransitionNotchVolume
3117 or key is self.outerTransitionNotchVolume
3118 ):
3119 topSurfacesDimTags.extend(value)
3121 sideSurfacesDimTags = []
3122 if i % 2 == 1:
3123 # Touches it tube and air tube
3124 bottomSurfacesDimTags.extend(
3125 topSurfaces[self.outerTerminalTubeVolume]
3126 )
3127 topSurfacesDimTags.extend(
3128 self.fundamentalSurfaces[self.outerTerminalTubeVolume]
3129 )
3131 if self.mesh.ti.structured:
3132 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(
3133 3
3134 )[-4:]
3135 else:
3136 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(
3137 3
3138 )[-1:]
3140 lastItTubeSurfsDimTags = gmsh.model.getBoundary(
3141 lastItTubeVolDimTags, oriented=False
3142 )
3143 lastItTubeSideSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3144 lastItTubeSurfsDimTags
3145 )
3146 sideSurfacesDimTags.extend(
3147 findOuterOnes(lastItTubeSideSurfsDimTags)
3148 )
3150 if self.mesh.ai.structured:
3151 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[
3152 -4:
3153 ]
3154 else:
3155 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[
3156 -1:
3157 ]
3158 lastAirTubeSurfsDimTags = gmsh.model.getBoundary(
3159 lastAirTubeVolDimTags, oriented=False
3160 )
3161 lastAirTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3162 lastAirTubeSurfsDimTags
3163 )
3164 sideSurfacesDimTags.extend(
3165 findOuterOnes(lastAirTubeSurfsDimTags, findInnerOnes=True)
3166 )
3168 else:
3169 # Touches ot tube and air cylinder
3170 bottomSurfacesDimTags.extend(
3171 topSurfaces[self.innerTerminalTubeVolume]
3172 )
3173 topSurfacesDimTags.extend(
3174 self.fundamentalSurfaces[self.innerTerminalTubeVolume]
3175 )
3176 if self.mesh.ti.structured:
3177 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(
3178 3
3179 )[-4:]
3180 else:
3181 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(
3182 3
3183 )[-1:]
3185 lastOtTubeSurfsDimTags = gmsh.model.getBoundary(
3186 lastOtTubeVolDimTags, oriented=False
3187 )
3188 lastOtTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3189 lastOtTubeSurfsDimTags
3190 )
3191 sideSurfacesDimTags.extend(
3192 findOuterOnes(lastOtTubeSurfsDimTags, findInnerOnes=True)
3193 )
3195 if self.mesh.ai.structured:
3196 lastAirCylinderVolDimTags = (
3197 self.centerAirCylinderVolume.getDimTags(3)[-4:]
3198 )
3199 else:
3200 lastAirCylinderVolDimTags = (
3201 self.centerAirCylinderVolume.getDimTags(3)[-1:]
3202 )
3204 lastAirCylinderSurfsDimTags = gmsh.model.getBoundary(
3205 lastAirCylinderVolDimTags, oriented=False
3206 )
3207 lastAirCylinderSurfsDimTags = findSurfacesWithNormalsOnXYPlane(
3208 lastAirCylinderSurfsDimTags
3209 )
3210 sideSurfacesDimTags.extend(
3211 findOuterOnes(lastAirCylinderSurfsDimTags)
3212 )
3214 allGapAirSurfacesDimTags = (
3215 bottomSurfacesDimTags + topSurfacesDimTags + sideSurfacesDimTags
3216 )
3218 # Technically, since all the boundary surfaces of the gap air volumes
3219 # are found here, we should be able to create the gap air volumes with
3220 # addSurfaceLoop and addVolume functions. However, when those are used,
3221 # Geometry.remove_all_duplicates() will indicate some
3222 # duplicates/ill-shaped geometry entities. The indication is
3223 # gmsh.model.occ.remove_all_duplicates() will change the geometry
3224 # (delete some volumes and create new ones), and I have always thought
3225 # that means there are big errors in the geometry and that geometry
3226 # should not be used.
3228 # Alternatively, using these surface tags, the gap air can be created
3229 # with fragment operations as well. Geometry.remove_all_duplicates()
3230 # will tell everything is fine when the fragment operation is used.
3232 # However, I checked manually as well, the way I am using the
3233 # addSurfaceLoop and addVolume should definitely work (because the end
3234 # result is the same with fragments), and I think it is a gmsh/occ
3235 # related problem. In the end, I realized creating the gap air with
3236 # addSurfaceLoop and addVolume won't even affect the mesh, and
3237 # everything seems conformal and nice. Since the fragment operation
3238 # is also very slow, I decided to use addSurfaceLoop and addVolume them.
3239 # However, I keep it as an option so that if the user feels something
3240 # funny about the geometry, the gap air can be created with fragment
3241 # operations as well.
3243 if not self.geo.ai.fragment:
3244 allGapAirSurfacesTags = [
3245 dimTag[1] for dimTag in allGapAirSurfacesDimTags
3246 ]
3247 surfaceLoop = gmsh.model.occ.addSurfaceLoop(allGapAirSurfacesTags)
3248 volume = gmsh.model.occ.addVolume([surfaceLoop])
3249 self.gapAirVolume.add([(3, volume)])
3251 else:
3252 # Save the surface tags for a fast fragment operation:
3253 self.gapAirSurfacesDimTags.append(allGapAirSurfacesDimTags)
3255 # self.extrudeSurfaces uses self.fundamentalSurfaces for extrusion and adds
3256 # the new volumes to the dimTags objects and returns the dictionary of the
3257 # new top surfaces. The new top surfaces then will be used in extrudeGapPart
3258 # method.
3259 topSurfaces = self.extrudeWindingPart()
3261 if i == 0:
3262 # If it is the first pancake coil, fundemental surfaces are extruded
3263 # downwards to create the bottom air volume and terminal volume.
3264 _ = self.extrudeGapPart(
3265 self.fundamentalSurfaces,
3266 -self.geo.ai.margin,
3267 terminalDimTagsObject=self.outerTerminalTubeVolume,
3268 firstTerminal=True,
3269 )
3271 if not i == self.geo.N - 1:
3272 # If it is not the last pancake coil, extrude the terminal surface to
3273 # create the next contactTerminalSurface and store the new volume in the
3274 # corresponding dimTags object.
3275 self.contactSurfaces = self.extrudeGapPart(topSurfaces)
3277 else:
3278 # If it is the last pancake coil, extrude the terminal surface all the
3279 # way up to the top and store the new volume in the corresponding
3280 # dimTags object.
3281 _ = self.extrudeGapPart(
3282 topSurfaces,
3283 self.geo.ai.margin,
3284 lastTerminal=True,
3285 )
3286 self.pancakeIndex = self.pancakeIndex + 1
3288 # Create the gap air volume:
3289 if self.geo.ai.fragment and self.geo.N > 1:
3290 self.generateGapAirWithFragment(self.gapAirSurfacesDimTags)
3292 logger.info(
3293 "Generating pancake coils has been finished in"
3294 f" {timeit.default_timer() - start_time:.2f} s."
3295 )
3297 def generateFundamentalSurfaces(self):
3298 """
3299 Generates the inner air, outer air, winding, contact layer, and terminal surfaces
3300 of the current pancake coil and returns them. It finds the z coordinate of the
3301 surfaces and the direction of the pancake coil, depending on the pancake index.
3303 :return: list of dimTags that contains fundamental surfaces
3304 :rtype: list[tuple[int, int]]
3305 """
3306 fundamentalSurfaces = {}
3308 # Select the direction of the spiral:
3309 if self.pancakeIndex % 2 == 0:
3310 spiralDirection = direction.ccw
3311 else:
3312 spiralDirection = direction.cw
3314 # Calculate the z coordinate of the surfaces:
3315 z = (
3316 -self.geo.ai.h / 2
3317 + self.geo.ai.margin
3318 + self.pancakeIndex * (self.geo.wi.h + self.geo.gap)
3319 )
3321 # Create the winding and contact layer surface:
3322 surface = spiralSurface(
3323 self.geo.wi.r_i,
3324 self.geo.wi.t,
3325 self.geo.ii.t,
3326 self.geo.wi.N,
3327 z,
3328 self.geo.wi.theta_i,
3329 self.geo.ti.transitionNotchAngle,
3330 spiralDirection,
3331 thinShellApproximation=self.geo.ii.tsa,
3332 )
3334 # Save the surface tags (if TSA, contactLayerSurfaceTags will be empty):
3335 fundamentalSurfaces[self.individualWinding[self.pancakeIndex]] = [
3336 (2, tag) for tag in surface.surfaceTags
3337 ]
3338 fundamentalSurfaces[self.individualContactLayer[self.pancakeIndex]] = [
3339 (2, tag) for tag in surface.contactLayerSurfaceTags
3340 ]
3342 if self.geo.ai.type == "cylinder":
3343 outerAirSurf = outerAirSurface(
3344 self.geo.ai.r,
3345 self.geo.ti.o.r,
3346 type="cylinder",
3347 divideIntoFourParts=self.mesh.ai.structured,
3348 divideTerminalPartIntoFourParts=self.mesh.ti.structured,
3349 )
3350 elif self.geo.ai.type == "cuboid":
3351 outerAirSurf = outerAirSurface(
3352 self.geo.ai.a / 2,
3353 self.geo.ti.o.r,
3354 type="cuboid",
3355 divideIntoFourParts=self.mesh.ai.structured,
3356 divideTerminalPartIntoFourParts=self.mesh.ti.structured,
3357 )
3359 outerTerminalSurf = outerTerminalSurface(
3360 self.geo.ti.o.r,
3361 self.geo.ti.o.t,
3362 divideIntoFourParts=self.mesh.ti.structured,
3363 )
3364 innerTerminalSurf = innerTerminalSurface(
3365 self.geo.ti.i.r,
3366 self.geo.ti.i.t,
3367 divideIntoFourParts=self.mesh.ti.structured,
3368 )
3369 innerAirSurf = innerAirSurface(
3370 self.geo.ti.i.r,
3371 divideIntoFourParts=self.mesh.ai.structured,
3372 divideTerminalPartIntoFourParts=self.mesh.ti.structured,
3373 )
3375 if self.contactSurfaces:
3376 # If self.contactSurfaces is not empty, it means that it is not the
3377 # first pancake coil. In that case, contactSurfaces should be used to
3378 # avoid surface duplication.
3380 # Create outer air:
3381 outerAirPSDimTags = self.contactSurfaces[self.outerAirTubeVolume]
3382 outerAirPSTags = [dimTag[1] for dimTag in outerAirPSDimTags]
3383 if self.geo.ai.shellTransformation:
3384 if self.geo.ai.type == "cuboid":
3385 cuboidShellDimTags1 = self.contactSurfaces[self.airShellVolumePart1]
3386 cuboidShellTags1 = [dimTag[1] for dimTag in cuboidShellDimTags1]
3387 cuboidShellDimTags2 = self.contactSurfaces[self.airShellVolumePart2]
3388 cuboidShellTags2 = [dimTag[1] for dimTag in cuboidShellDimTags2]
3389 cuboidShellDimTags3 = self.contactSurfaces[self.airShellVolumePart3]
3390 cuboidShellTags3 = [dimTag[1] for dimTag in cuboidShellDimTags3]
3391 cuboidShellDimTags4 = self.contactSurfaces[self.airShellVolumePart4]
3392 cuboidShellTags4 = [dimTag[1] for dimTag in cuboidShellDimTags4]
3393 outerAirSurf.setPrecreatedSurfaceTags(
3394 outerAirPSTags,
3395 cuboidShellTags1=cuboidShellTags1,
3396 cuboidShellTags2=cuboidShellTags2,
3397 cuboidShellTags3=cuboidShellTags3,
3398 cuboidShellTags4=cuboidShellTags4,
3399 )
3400 elif self.geo.ai.type == "cylinder":
3401 cylinderShellDimTags = self.contactSurfaces[self.airShellVolume]
3402 cylinderShellTags = [dimTag[1] for dimTag in cylinderShellDimTags]
3403 outerAirSurf.setPrecreatedSurfaceTags(
3404 outerAirPSTags,
3405 cylinderShellTags=cylinderShellTags,
3406 )
3407 else:
3408 outerAirSurf.setPrecreatedSurfaceTags(outerAirPSTags)
3410 # Create inner air:
3411 innerAirPSDimTags = self.contactSurfaces[self.centerAirCylinderVolume]
3412 innerAirPSTags = [dimTag[1] for dimTag in innerAirPSDimTags]
3413 innerAirSurf.setPrecreatedSurfaceTags(innerAirPSTags)
3415 if self.pancakeIndex % 2 == 0:
3416 # In this case, we should create all the surfaces for the inner terminal
3417 # but not for outer terminal. Because it is a pancake coil with an even
3418 # index (self.pancakeIndex%2==0) which means that it is connected to the
3419 # previous pancake coil with outer terminal and the outer terminal
3420 # surface is ready (extruded before).
3422 # Create outer terminal:
3423 outerTerminalTubePSDimTags = self.contactSurfaces[
3424 self.outerTerminalTubeVolume
3425 ]
3426 outerTerminalTubePSTags = [
3427 dimTag[1] for dimTag in outerTerminalTubePSDimTags
3428 ]
3429 outerTerminalSurf.createWithWindingAndTubeTags(
3430 surface, outerTerminalTubePSTags, self.pancakeIndex
3431 )
3433 # Create inner terminal:
3434 innerTerminalSurf.createWithInnerAirAndWinding(
3435 innerAirSurf, surface, self.pancakeIndex
3436 )
3438 else:
3439 # In this case, we should create all the surfaces for the outer terminal
3440 # but not for inner terminal. Because it is a pancake coil with an odd
3441 # index (self.pancakeIndex%2==1) which means that it is connected to the
3442 # previous pancake coil with inner terminal and the inner terminal
3443 # surface is ready (extruded before).
3445 # Create outer terminal:
3446 outerTerminalSurf.createWithOuterAirAndWinding(
3447 outerAirSurf, surface, self.pancakeIndex
3448 )
3450 # Create inner terminal:
3451 innerTerminalTubePSDimTags = self.contactSurfaces[
3452 self.innerTerminalTubeVolume
3453 ]
3454 innerTerminalTubePSTags = [
3455 dimTag[1] for dimTag in innerTerminalTubePSDimTags
3456 ]
3457 innerTerminalSurf.createWithWindingAndTubeTags(
3458 surface, innerTerminalTubePSTags, self.pancakeIndex
3459 )
3461 else:
3462 # If self.contactSurfaces is empty, it means that it is the first pancake
3463 # coil. In that case, the surfaces should be created from scratch.
3465 if self.geo.ai.shellTransformation:
3466 if self.geo.ai.type == "cuboid":
3467 outerAirSurf.createFromScratch(
3468 z,
3469 shellTransformation=True,
3470 shellRadius=self.geo.ai.shellSideLength / 2,
3471 )
3472 else:
3473 outerAirSurf.createFromScratch(
3474 z,
3475 shellTransformation=True,
3476 shellRadius=self.geo.ai.shellOuterRadius,
3477 )
3478 else:
3479 outerAirSurf.createFromScratch(z)
3481 innerAirSurf.createFromScratch(z)
3482 outerTerminalSurf.createWithOuterAirAndWinding(
3483 outerAirSurf, surface, self.pancakeIndex
3484 )
3485 innerTerminalSurf.createWithInnerAirAndWinding(
3486 innerAirSurf, surface, self.pancakeIndex
3487 )
3489 # Save the surface tags:
3490 fundamentalSurfaces[self.outerAirTubeVolume] = [
3491 (2, tag) for tag in outerAirSurf.surfaceTags
3492 ]
3494 fundamentalSurfaces[self.centerAirCylinderVolume] = [
3495 (2, tag) for tag in innerAirSurf.surfaceTags
3496 ]
3498 fundamentalSurfaces[self.outerTerminalTubeVolume] = [
3499 (2, tag) for tag in outerTerminalSurf.tubeSurfaceTags
3500 ]
3501 fundamentalSurfaces[self.outerTerminalTouchingVolume] = [
3502 (2, tag) for tag in outerTerminalSurf.nontubeSurfaceTags
3503 ]
3505 fundamentalSurfaces[self.innerTerminalTubeVolume] = [
3506 (2, tag) for tag in innerTerminalSurf.tubeSurfaceTags
3507 ]
3508 fundamentalSurfaces[self.innerTerminalTouchingVolume] = [
3509 (2, tag) for tag in innerTerminalSurf.nontubeSurfaceTags
3510 ]
3511 fundamentalSurfaces[self.innerTransitionNotchVolume] = [
3512 (2, tag) for tag in surface.innerNotchSurfaceTags
3513 ]
3514 fundamentalSurfaces[self.outerTransitionNotchVolume] = [
3515 (2, tag) for tag in surface.outerNotchSurfaceTags
3516 ]
3518 if self.geo.ai.shellTransformation:
3519 if self.geo.ai.type == "cuboid":
3520 fundamentalSurfaces[self.airShellVolumePart1] = [
3521 (2, tag) for tag in outerAirSurf.shellTagsPart1
3522 ]
3523 fundamentalSurfaces[self.airShellVolumePart2] = [
3524 (2, tag) for tag in outerAirSurf.shellTagsPart2
3525 ]
3526 fundamentalSurfaces[self.airShellVolumePart3] = [
3527 (2, tag) for tag in outerAirSurf.shellTagsPart3
3528 ]
3529 fundamentalSurfaces[self.airShellVolumePart4] = [
3530 (2, tag) for tag in outerAirSurf.shellTagsPart4
3531 ]
3532 elif self.geo.ai.type == "cylinder":
3533 fundamentalSurfaces[self.airShellVolume] = [
3534 (2, tag) for tag in outerAirSurf.shellTags
3535 ]
3537 return fundamentalSurfaces
3539 def extrudeGapPart(
3540 self,
3541 surfacesDict,
3542 tZ: float = None,
3543 terminalDimTagsObject: dimTags = None,
3544 firstTerminal=False,
3545 lastTerminal=False,
3546 ):
3547 """
3548 Extrudes the given surfaces dimTags dictionary to a given height (tZ) and adds
3549 the created volumes to the corresponding dictionary keys (dimTags objects). It
3550 returns the extrusion's top surfaces as a dictionary again, where the keys are
3551 the corresponding dimTagsObjects and the values are the dimTags of the surfaces.
3553 If tZ is not given, then it is set to the gap height (self.geo.gap). This is the
3554 default value used for connecting the pancake coils. Only for the creation of
3555 the first and the last pancake coils different tZ values are used.
3557 If terminalDimTagsObject is not given, then the created volume is added
3558 automatically to the innerTerminalVolume or outerTerminalVolume dimTags object,
3559 depending on the value of self.pancakeIndex. However, giving
3560 terminalDimTagsObject is necessary for creating the first and the last terminal.
3561 Otherwise, finding out the correct terminal dimTagsObject would be very
3562 challenging.
3564 :param surfaces: the surface dimTag dictionary to be extruded. The keys are the
3565 dimTags objects and the values are the dimTags of the surfaces. The
3566 keys are used to easily add the corresponding volumes to the correct
3567 dimTags objects
3568 :type surfaces: dict[dimTags, list[tuple[int, int]]]
3569 :param tZ: the height of the extrusion
3570 :type tZ: float, optional
3571 :param terminalDimTagsObject: the dimTags object of the terminal to be extruded
3572 :type terminalDimTagsObject: dimTags, optional
3573 :return: top surfaces of the extrusion as a dictionary where the keys are the
3574 dimTags objects and the values are the dimTags of the surfaces
3575 :rtype: dict[dimTags, list[tuple[int, int]]]
3576 """
3577 bottomPart = False
3578 topPart = False
3579 if tZ is None:
3580 tZ = self.geo.gap
3581 elif tZ < 0:
3582 bottomPart = True
3583 elif tZ > 0:
3584 topPart = True
3586 if terminalDimTagsObject is None:
3587 # terminalDimTagsObject needs to be given for the first terminal that is
3588 # extruded downwards.
3589 if self.pancakeIndex % 2 == 0:
3590 terminalDimTagsObject = self.innerTerminalTubeVolume
3591 else:
3592 terminalDimTagsObject = self.outerTerminalTubeVolume
3594 # if terminalDimTagsObject is self.innerTerminalVolume:
3595 # otherTerminal = self.outerTerminalVolume
3596 # else:
3597 # otherTerminal = self.innerTerminalVolume
3599 # Create the list of surfaces to be extruded:
3600 listOfDimTags = []
3601 listOfDimTagsObjects = []
3602 listOfDimTagsForTopSurfaces = []
3603 if topPart:
3604 # Then in this case, most of the surfaces should be added to the air volumes
3605 # instead of the terminal, winding, and contact layer volumes.
3606 for key, dimTagsList in surfacesDict.items():
3607 if key is self.individualWinding[self.pancakeIndex]:
3608 dimTagsObjects = [self.topAirPancakeWindingExtursionVolume] * len(
3609 dimTagsList
3610 )
3611 elif key is self.individualContactLayer[self.pancakeIndex]:
3612 dimTagsObjects = [
3613 self.topAirPancakeContactLayerExtursionVolume
3614 ] * len(dimTagsList)
3615 elif (
3616 key is terminalDimTagsObject
3617 or key is self.airShellVolume
3618 or key is self.airShellVolumePart1
3619 or key is self.airShellVolumePart2
3620 or key is self.airShellVolumePart3
3621 or key is self.airShellVolumePart4
3622 or key is self.outerAirTubeVolume
3623 or key is self.centerAirCylinderVolume
3624 ):
3625 dimTagsObjects = [key] * len(dimTagsList)
3626 else:
3627 # key is self.outerTerminalTouchingVolume
3628 # or key is self.innerTerminalTouchingVolume
3629 # or key is (other terminal's tube volume)
3630 dimTagsObjects = [self.topAirTerminalsExtrusionVolume] * len(
3631 dimTagsList
3632 )
3633 if (
3634 key is self.innerTerminalTubeVolume
3635 or key is self.outerTerminalTubeVolume
3636 ):
3637 dimTagsObjects = [
3638 self.topAirTubeTerminalsExtrusionVolume
3639 ] * len(dimTagsList)
3641 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(
3642 dimTagsList
3643 )
3644 listOfDimTags = listOfDimTags + dimTagsList
3645 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3646 elif bottomPart:
3647 # Then in this case, most of the surfaces should be added to the air volumes
3648 # instead of the terminal, winding, and contact layer volumes.
3649 for key, dimTagsList in surfacesDict.items():
3650 if key is self.individualWinding[self.pancakeIndex]:
3651 dimTagsObjects = [
3652 self.bottomAirPancakeWindingExtursionVolume
3653 ] * len(dimTagsList)
3654 elif key is self.individualContactLayer[self.pancakeIndex]:
3655 dimTagsObjects = [
3656 self.bottomAirPancakeContactLayerExtursionVolume
3657 ] * len(dimTagsList)
3658 elif (
3659 key is terminalDimTagsObject
3660 or key is self.airShellVolume
3661 or key is self.airShellVolumePart1
3662 or key is self.airShellVolumePart2
3663 or key is self.airShellVolumePart3
3664 or key is self.airShellVolumePart4
3665 or key is self.outerAirTubeVolume
3666 or key is self.centerAirCylinderVolume
3667 ):
3668 dimTagsObjects = [key] * len(dimTagsList)
3669 else:
3670 # key is self.outerTerminalTouchingVolume
3671 # or key is self.innerTerminalTouchingVolume
3672 # or key is (other terminal's tube volume)
3673 dimTagsObjects = [self.bottomAirTerminalsExtrusionVolume] * len(
3674 dimTagsList
3675 )
3676 if (
3677 key is self.innerTerminalTubeVolume
3678 or key is self.outerTerminalTubeVolume
3679 ):
3680 dimTagsObjects = [
3681 self.bottomAirTubeTerminalsExtrusionVolume
3682 ] * len(dimTagsList)
3684 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(
3685 dimTagsList
3686 )
3687 listOfDimTags = listOfDimTags + dimTagsList
3688 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3689 else:
3690 for key, dimTagsList in surfacesDict.items():
3691 if (
3692 key is self.outerAirTubeVolume
3693 or key is self.centerAirCylinderVolume
3694 or key is self.airShellVolume
3695 or key is self.airShellVolumePart1
3696 or key is self.airShellVolumePart2
3697 or key is self.airShellVolumePart3
3698 or key is self.airShellVolumePart4
3699 or key is terminalDimTagsObject
3700 ):
3701 dimTagsObjects = [key] * len(dimTagsList)
3703 listOfDimTags = listOfDimTags + dimTagsList
3704 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3706 listOfDimTagsForTopSurfaces = listOfDimTagsObjects
3708 extrusionResult = dimTags()
3709 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, tZ))
3711 # Add the created volumes to the corresponding dimTags objects:
3712 volumeDimTags = extrusionResult.getDimTags(3)
3713 for i, volumeDimTag in enumerate(volumeDimTags):
3714 listOfDimTagsObjects[i].add(volumeDimTag)
3716 if firstTerminal:
3717 self.firstTerminalVolume.add(terminalDimTagsObject.getDimTags(3))
3718 elif lastTerminal:
3719 self.lastTerminalVolume.add(terminalDimTagsObject.getDimTags(3))
3721 topSurfacesDimTags = extrusionResult.getExtrusionTop()
3722 topSurfacesDict = {}
3723 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):
3724 if listOfDimTagsObjects[i] in topSurfacesDict:
3725 topSurfacesDict[listOfDimTagsForTopSurfaces[i]].append(topSurfaceDimTag)
3726 else:
3727 topSurfacesDict[listOfDimTagsForTopSurfaces[i]] = [topSurfaceDimTag]
3729 return topSurfacesDict
3731 def extrudeWindingPart(self):
3732 """
3733 Extrudes all the fundamental surfaces of the pancake coil by self.geo.wi.h and
3734 returns the next connection terminal's top surface dimTag, and other air dimTags
3735 in a dictionary so that they can be further extruded.
3737 :return: dictionary of top surfaces where the keys are the dimTags objects and
3738 the values are the dimTags of the surfaces
3739 :rtype: dict[dimTags, list[tuple[int, int]]]
3740 """
3741 # Create the list of surfaces to be extruded:
3742 listOfDimTags = []
3743 listOfDimTagsObjects = []
3744 for key, dimTagsList in self.fundamentalSurfaces.items():
3745 dimTagsObjects = [key] * len(dimTagsList)
3747 listOfDimTags = listOfDimTags + dimTagsList
3748 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects
3750 # Extrude the fundamental surfaces:
3751 extrusionResult = dimTags()
3752 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, self.geo.wi.h))
3754 # Add the created volumes to the corresponding dimTags objects:
3755 volumes = extrusionResult.getDimTags(3)
3756 for i, volumeDimTag in enumerate(volumes):
3757 listOfDimTagsObjects[i].add(volumeDimTag)
3759 if self.pancakeIndex == 0:
3760 # Note the first pancake (sometimes useful for creating regions in the
3761 # meshing part):
3762 for i, volumeDimTag in enumerate(volumes):
3763 if listOfDimTagsObjects[i].parentName == self.geo.ti.o.name:
3764 self.firstTerminalVolume.add(volumeDimTag)
3766 # Not elif! Because the first pancake coil is also the last pancake coil if
3767 # there is only one pancake coil.
3768 if self.pancakeIndex == self.geo.N - 1:
3769 # Note the last pancake (sometimes useful for creating regions in the
3770 # meshing part):
3771 for i, volumeDimTag in enumerate(volumes):
3772 if (
3773 self.pancakeIndex % 2 == 1
3774 and listOfDimTagsObjects[i].parentName == self.geo.ti.o.name
3775 ):
3776 self.lastTerminalVolume.add(volumeDimTag)
3777 elif (
3778 self.pancakeIndex % 2 == 0
3779 and listOfDimTagsObjects[i].parentName == self.geo.ti.i.name
3780 ):
3781 self.lastTerminalVolume.add(volumeDimTag)
3783 # Return the top surfaces:
3784 # Add the created top surfaces to a new dictionary:
3785 topSurfacesDimTags = extrusionResult.getExtrusionTop()
3786 topSurfaces = {}
3787 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):
3788 if listOfDimTagsObjects[i] in topSurfaces:
3789 topSurfaces[listOfDimTagsObjects[i]].append(topSurfaceDimTag)
3790 else:
3791 topSurfaces[listOfDimTagsObjects[i]] = [topSurfaceDimTag]
3793 return topSurfaces
3795 def generateGapAirWithFragment(
3796 self, gapAirSurfacesDimTags: List[List[Tuple[int, int]]]
3797 ):
3798 """
3799 A class to fill the gap between the multiple pancake coils with air. First, it
3800 creates a dummy cylinder with the same radius as the outer terminal's outer
3801 radius. Then using gapAirSurfacesDimTags, gmsh.model.occ.fragment() operation is
3802 applied to the dummy cylinder volume in a for loop to create the gap air
3803 volumes. After each fragment operation, one of the volumes created is removed
3804 because it is the solid volume which is the combination of windings,
3805 contact layers, and terminals. In the end, dummy cylinder is removed as well.
3808 WARNING:
3809 Currently, this method doesn't work.
3811 :param geometry: geometry information
3812 :param pancakeCoils: pancakeCoilsWithAir object
3813 :type pancakeCoils: pancakeCoilsWithAir
3814 """
3815 logger.info("Generating gap air has been started.")
3816 start_time = timeit.default_timer()
3818 # Create the dummy air volume:
3819 dummyAir = gmsh.model.occ.addCylinder(
3820 0,
3821 0,
3822 -self.geo.ai.h / 2,
3823 0,
3824 0,
3825 self.geo.ai.h,
3826 self.geo.ti.o.r,
3827 )
3829 toBeDeletedDimTags = []
3830 gapAirVolumesCurrentDimTags = []
3831 for i in range(len(gapAirSurfacesDimTags)):
3832 # Get the outer surfaces of the pancake coils for cutting the pancake coils
3833 # from the dummy air. The outer surfaces are used instead of pancake volumes
3834 # to reduce the amount of work for gmsh. It makes it significantly faster.
3835 # if len(gapAirSurfacesDimTags[i]) !=12:
3836 fragmentResults = gmsh.model.occ.fragment(
3837 [(3, dummyAir)],
3838 gapAirSurfacesDimTags[i],
3839 removeObject=False,
3840 removeTool=False,
3841 )
3842 fragmentVolumeResultsDimTags = fragmentResults[1][0]
3843 toBeDeletedDimTags.append(fragmentVolumeResultsDimTags[0])
3844 gapAirVolumesCurrentDimTags.append(fragmentVolumeResultsDimTags[1])
3846 toBeDeletedDimTags.append((3, dummyAir))
3847 # Fragmnet operation both creates the air volume and solid pancake coils volume
3848 # because the surfaces are used for cutting. Therefore, the solid pancake coils
3849 # volume should be removed from the fragment results:
3850 gmsh.model.occ.remove(toBeDeletedDimTags)
3852 # Add results to the air volume storage. After the geometry is saves as a .brep
3853 # file, and loaded back, the gaps between the tags are avoided by moving the
3854 # the other tags. Therefore, this is how the tags are stored:
3855 toBeDeletedTags = [dimTag[1] for dimTag in toBeDeletedDimTags]
3856 volumeTagsStart = min(toBeDeletedTags)
3857 numberOfGapAirVolumes = len(gapAirVolumesCurrentDimTags)
3858 gapAirVolumesToBeSaved = [
3859 (3, volumeTagsStart + i) for i in range(numberOfGapAirVolumes)
3860 ]
3862 # For debugging purposes, physical groups are being created in the geometry
3863 # generation process as well. Normally, it us done during meshing because BREP
3864 # files cannot store physical groups. Since the tags above (airVolumes) will be
3865 # valid after only saving the geometry as a BREP file and loading it back, the
3866 # current tags are given to the airVolume.add() method as well. This is done to
3867 # be able to create the correct physical group.
3868 self.gapAirVolume.add(
3869 dimTagsList=gapAirVolumesToBeSaved,
3870 dimTagsListForPG=gapAirVolumesCurrentDimTags,
3871 )
3873 logger.info(
3874 "Generating gap air has been finished in"
3875 f" {timeit.default_timer() - start_time:.2f} s."
3876 )
3879class Geometry(Base):
3880 """
3881 Main geometry class for Pancake3D.
3883 :param fdm: FiQuS data model
3884 :param geom_folder: folder where the geometry files are saved
3885 :type geom_folder: str
3886 :param mesh_folder: folder where the mesh files are saved
3887 :type mesh_folder: str
3888 :param solution_folder: folder where the solution files are saved
3889 :type solution_folder: str
3890 """
3892 def __init__(
3893 self,
3894 fdm,
3895 geom_folder,
3896 mesh_folder,
3897 solution_folder,
3898 ) -> None:
3899 super().__init__(fdm, geom_folder, mesh_folder, solution_folder)
3901 # Clear if there is any existing dimTags storage:
3902 dimTagsStorage.clear()
3904 # Start GMSH:
3905 self.gu = GmshUtils(self.geom_folder)
3906 self.gu.initialize()
3908 # To speed up the GUI:
3909 gmsh.option.setNumber("Geometry.NumSubEdges", 10)
3911 # To see the surfaces in a better way in GUI:
3912 gmsh.option.setNumber("Geometry.SurfaceType", 1)
3914 # Makes the gmsh.mode.occ.cut function faster:
3915 gmsh.option.setNumber("Geometry.OCCParallel", 1)
3917 # To avoid any unwanted modifications to the geometry, the automatic fixing of
3918 # the geometry is disabled:
3919 gmsh.option.setNumber("Geometry.OCCAutoFix", 0)
3921 # Set the tolerance:
3922 if self.geo.dimTol < gmsh.option.getNumber("Geometry.Tolerance"):
3923 gmsh.option.setNumber("Geometry.Tolerance", self.geo.dimTol)
3925 gmsh.option.setNumber("Geometry.ToleranceBoolean", self.geo.dimTol)
3927 spiralCurve.sectionsPerTurn = self.geo.wi.spt
3928 spiralCurve.curvesPerTurn = self.geo.wi.NofVolPerTurn
3930 def generate_geometry(self):
3931 """
3932 Generates geometry and saves it as a .brep file.
3935 """
3936 logger.info(
3937 f"Generating Pancake3D geometry ({self.brep_file}) has been started."
3938 )
3939 start_time = timeit.default_timer()
3941 self.pancakeCoil = pancakeCoilsWithAir(self.geo, self.mesh)
3943 gmsh.model.occ.synchronize()
3944 gmsh.write(self.brep_file)
3946 logger.info(
3947 f"Generating Pancake3D geometry ({self.brep_file}) has been finished in"
3948 f" {timeit.default_timer() - start_time:.2f} s."
3949 )
3951 def load_geometry(self):
3952 """
3953 Loads geometry from .brep file.
3954 """
3955 logger.info("Loading Pancake3D geometry has been started.")
3956 start_time = timeit.default_timer()
3958 previousGeo = FilesAndFolders.read_data_from_yaml(
3959 self.geometry_data_file, Pancake3DGeometry
3960 )
3962 if previousGeo.dict() != self.geo.dict():
3963 raise ValueError(
3964 "Geometry data has been changed. Please regenerate the geometry or load"
3965 " the previous geometry data."
3966 )
3968 gmsh.clear()
3969 gmsh.model.occ.importShapes(self.brep_file, format="brep")
3970 gmsh.model.occ.synchronize()
3972 logger.info(
3973 "Loading Pancake3D geometry has been finished in"
3974 f" {timeit.default_timer() - start_time:.2f} s."
3975 )
3977 def generate_vi_file(self):
3978 """
3979 Generates volume information file. Volume information file stores dimTags of all
3980 the stored volumes in the geometry. Without this file, regions couldn't be
3981 created, meaning that finite element simulation cannot be done.
3983 The file extension is custom because users are not supposed to edit or open this
3984 file, and it makes it intuitively clear that it is a volume information file.
3985 """
3986 logger.info(
3987 f"Generating volume information file ({self.vi_file}) has been started."
3988 )
3989 start_time = timeit.default_timer()
3991 dimTagsDict = dimTagsStorage.getDimTagsDict()
3992 json.dump(
3993 dimTagsDict,
3994 open(self.vi_file, "w"),
3995 )
3997 logger.info(
3998 "Generating volume information file has been finished in"
3999 f" {timeit.default_timer() - start_time:.2f} s."
4000 )
4002 if self.geo_gui:
4003 self.generate_physical_groups()
4004 self.gu.launch_interactive_GUI()
4005 else:
4006 gmsh.clear()
4007 gmsh.finalize()
4009 @staticmethod
4010 def generate_physical_groups():
4011 """
4012 Generates physical groups. Physical groups are not saved in the BREP file but
4013 it can be useful for debugging purposes.
4016 """
4017 gmsh.model.occ.synchronize()
4019 dimTagsDict = dimTagsStorage.getDimTagsDict(forPhysicalGroups=True)
4021 for key, value in dimTagsDict.items():
4022 tags = [dimTag[1] for dimTag in value]
4023 gmsh.model.addPhysicalGroup(
4024 3,
4025 tags,
4026 name=key,
4027 )
4029 @staticmethod
4030 def remove_all_duplicates():
4031 """
4032 Removes all the duplicates and then prints the entities that are created or
4033 removed during the operation. It prints the line number where the function is
4034 called as well. This function is helpful for debugging. Finding duplicates means
4035 there is a problem in geometry creation logic, and the meshes will not be
4036 conformal. It shouldn't be used in the final version of the code since removing
4037 duplicates is computationally expensive, and there shouldn't be duplicates at
4038 all.
4040 WARNING:
4041 This function currently does not work properly. It is not recommended to use
4042 right now. It finds duplicates even if there are no duplicates (topology
4043 problems).
4044 """
4046 logger.info(f"Removing all the duplicates has been started.")
4047 start_time = timeit.default_timer()
4049 gmsh.model.occ.synchronize()
4050 oldEntities = []
4051 oldEntities.extend(gmsh.model.getEntities(3))
4052 oldEntities.extend(gmsh.model.getEntities(2))
4053 oldEntities.extend(gmsh.model.getEntities(1))
4054 oldEntities.extend(gmsh.model.getEntities(0))
4055 oldEntities = set(oldEntities)
4057 gmsh.model.occ.removeAllDuplicates()
4059 gmsh.model.occ.synchronize()
4060 newEntities = []
4061 newEntities.extend(gmsh.model.getEntities(3))
4062 newEntities.extend(gmsh.model.getEntities(2))
4063 newEntities.extend(gmsh.model.getEntities(1))
4064 newEntities.extend(gmsh.model.getEntities(0))
4065 newEntities = set(newEntities)
4066 NewlyCreated = newEntities - oldEntities
4067 Removed = oldEntities - newEntities
4069 frameinfo = getframeinfo(currentframe().f_back)
4071 if len(NewlyCreated) > 0 or len(Removed) > 0:
4072 logger.warning(f"Duplicates found! Line: {frameinfo.lineno}")
4073 logger.warning(f"{len(NewlyCreated)}NewlyCreated = {list(NewlyCreated)}")
4074 logger.warning(f"{len(Removed)}Removed = {list(Removed)}")
4075 else:
4076 logger.info(f"No duplicates found! Line: {frameinfo.lineno}")
4078 logger.info(
4079 "Removing all the duplicates has been finished in"
4080 f" {timeit.default_timer() - start_time:.2f} s."
4081 )