Coverage for fiqus/mesh_generators/MeshPancake3D.py: 92%
1023 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-05-04 03:30 +0200
1import timeit
2import json
3import logging
4import math
5from enum import Enum
6import operator
7import itertools
8import re
10import gmsh
11import scipy.integrate
12import numpy as np
14from fiqus.data import RegionsModelFiQuS
15from fiqus.utils.Utils import GmshUtils
16from fiqus.data.RegionsModelFiQuS import RegionsModel
17from fiqus.mains.MainPancake3D import Base
18from fiqus.utils.Utils import FilesAndFolders
19from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry, Pancake3DMesh
22logger = logging.getLogger(__name__)
25class regionType(str, Enum):
26 """
27 A class to specify region type easily in the regions class.
28 """
30 powered = "powered"
31 insulator = "insulator"
32 air = "air"
33 air_far_field = "air_far_field"
36class entityType(str, Enum):
37 """
38 A class to specify entity type easily in the regions class.
39 """
41 vol = "vol"
42 vol_in = "vol_in"
43 vol_out = "vol_out"
44 surf = "surf"
45 surf_th = "surf_th"
46 surf_in = "surf_in"
47 surf_out = "surf_out"
48 surf_ext = "surf_ext"
49 cochain = "cochain"
50 curve = "curve"
51 point = "point"
54class regions:
55 """
56 A class to generate physical groups in GMSH and create the corresponding regions
57 file. The regions file is the file where the region tags are stored in the FiQuS
58 regions data model convention. The file is used to template the *.pro file (GetDP
59 input file).
60 """
62 def __init__(self):
63 # Main types of entities:
64 # The keys are the FiQuS region categories, and the values are the corresponding
65 # GMSH entity type.
66 self.entityMainType = {
67 "vol": "vol",
68 "vol_in": "vol",
69 "vol_out": "vol",
70 "surf": "surf",
71 "surf_th": "surf",
72 "surf_in": "surf",
73 "surf_out": "surf",
74 "surf_ext": "surf",
75 "cochain": "curve",
76 "curve": "curve",
77 "point": "point",
78 }
80 # Dimensions of entity types:
81 self.entityDim = {"vol": 3, "surf": 2, "curve": 1, "point": 0}
83 # Keys for regions file. The keys are appended to the name of the regions
84 # accordingly.
85 self.entityKey = {
86 "vol": "",
87 "vol_in": "",
88 "vol_out": "",
89 "surf": "_bd",
90 "surf_th": "_bd",
91 "surf_in": "_in",
92 "surf_out": "_out",
93 "surf_ext": "_ext",
94 "cochain": "_cut",
95 "curve": "_curve",
96 "point": "_point",
97 }
99 # FiQuS convetion for region numbers:
100 self.regionTags = {
101 "vol": 1000000, # volume region tag start
102 "surf": 2000000, # surface region tag start
103 "curve": 3000000, # curve region tag start
104 "point": 4000000, # point region tag start
105 }
107 # Initialize the regions model:
108 self.rm = RegionsModelFiQuS.RegionsModel()
110 # This is used because self.rm.powered["Pancake3D"] is not initialized in
111 # RegionsModelFiQuS.RegionsModel. It should be fixed in the future.
112 self.rm.powered["Pancake3D"] = RegionsModelFiQuS.Powered()
114 # Initializing the required variables (air and powered.vol_in and
115 # powered.vol_out are not initialized because they are not lists but numbers):
116 self.rm.powered["Pancake3D"].vol.names = []
117 self.rm.powered["Pancake3D"].vol.numbers = []
118 self.rm.powered["Pancake3D"].surf.names = []
119 self.rm.powered["Pancake3D"].surf.numbers = []
120 self.rm.powered["Pancake3D"].surf_th.names = []
121 self.rm.powered["Pancake3D"].surf_th.numbers = []
122 self.rm.powered["Pancake3D"].surf_in.names = []
123 self.rm.powered["Pancake3D"].surf_in.numbers = []
124 self.rm.powered["Pancake3D"].surf_out.names = []
125 self.rm.powered["Pancake3D"].surf_out.numbers = []
126 self.rm.powered["Pancake3D"].curve.names = []
127 self.rm.powered["Pancake3D"].curve.numbers = []
129 self.rm.insulator.vol.names = []
130 self.rm.insulator.vol.numbers = []
131 self.rm.insulator.surf.names = []
132 self.rm.insulator.surf.numbers = []
133 self.rm.insulator.curve.names = []
134 self.rm.insulator.curve.numbers = []
136 self.rm.air_far_field.vol.names = []
137 self.rm.air_far_field.vol.numbers = []
139 self.rm.air.cochain.names = []
140 self.rm.air.cochain.numbers = []
141 self.rm.air.point.names = []
142 self.rm.air.point.numbers = []
144 def addEntities(
145 self, name, entityTags, regionType: regionType, entityType: entityType
146 ):
147 """
148 Add entities as a physical group in GMSH and add the corresponding region to the
149 regions file data.
151 :param name: Name of the region (entityKey will be appended).
152 :type name: str
153 :param entityTags: Tags of the entities to be added as a physical group.
154 :type entityTags: list of integers (tags)
155 :param regionType: Type of the region. regionType class should be used.
156 :type regionType: regionType
157 :param entityType: Type of the entity. entityType class should be used.
158 :type entityType: entityType
159 """
160 if not isinstance(entityTags, list):
161 entityTags = [entityTags]
163 name = name + self.entityKey[entityType]
164 mainType = self.entityMainType[entityType]
165 dim = self.entityDim[mainType]
166 regionTag = self.regionTags[mainType]
168 if regionType is regionType.powered:
169 if entityType is entityType.vol_in or entityType is entityType.vol_out:
170 getattr(self.rm.powered["Pancake3D"], entityType).name = name
171 getattr(self.rm.powered["Pancake3D"], entityType).number = regionTag
173 else:
174 getattr(self.rm.powered["Pancake3D"], entityType).names.append(name)
175 getattr(self.rm.powered["Pancake3D"], entityType).numbers.append(
176 regionTag
177 )
178 elif regionType is regionType.insulator:
179 getattr(self.rm.insulator, entityType).names.append(name)
180 getattr(self.rm.insulator, entityType).numbers.append(regionTag)
181 elif regionType is regionType.air:
182 if entityType is entityType.cochain or entityType is entityType.point:
183 getattr(self.rm.air, entityType).names.append(name)
184 getattr(self.rm.air, entityType).numbers.append(regionTag)
185 else:
186 getattr(self.rm.air, entityType).name = name
187 getattr(self.rm.air, entityType).number = regionTag
188 elif regionType is regionType.air_far_field:
189 getattr(self.rm.air_far_field, entityType).names.append(name)
190 getattr(self.rm.air_far_field, entityType).numbers.append(regionTag)
192 gmsh.model.addPhysicalGroup(dim=dim, tags=entityTags, tag=regionTag, name=name)
193 self.regionTags[mainType] = self.regionTags[mainType] + 1
195 def generateRegionsFile(self, filename):
196 """
197 Generate the regions file from the final data.
199 :param filename: Name of the regions file (with extension).
200 :type filename: str
201 """
202 FilesAndFolders.write_data_to_yaml(filename, self.rm.model_dump())
205class curveType(Enum):
206 """
207 A class to specify curve type easily in the windingCurve class.
208 """
210 axial = 0
211 horizontal = 1
212 spiralArc = 2
213 circle = 3
216class curve:
217 """
218 Even though volume tags can be stored in a volume information file and can be used
219 after reading the BREP (geometry) file, surface tags and line tags cannot be stored
220 because their tags will be changed. However, we need to know which line is which to
221 create a proper mesh. For example, we would like to know which lines are on the XY
222 plane, which lines are straight, which lines are spirals, etc.
224 This class is created for recognizing lines of winding, contact layer, and top/bottom
225 air volumes (because they are extrusions of winding and contact layer). Line tags of
226 the volumes can be easily accessed with gmsh.model.getBoundary() function. Then a
227 line tag can be used to create an instance of this object. The class will analyze
228 the line's start and end points and decide if it's a spiral, axial, or horizontal
229 curve. Then, it calculates the length of the line. This information is required to
230 create a structured mesh for winding, contact layer, and top/bottom air volumes.
232 Every windingCurve object is a line that stores the line's type and length.
234 :param tag: Tag of the line.
235 :type tag: int
236 :param geometryData: Geometry data object.
237 """
239 def __init__(self, tag, geometryData):
240 self.geo = geometryData
242 self.tag = tag
244 pointDimTags = gmsh.model.getBoundary(
245 [(1, self.tag)], oriented=False, combined=True
246 )
247 self.pointTags = [dimTag[1] for dimTag in pointDimTags]
249 # Get the positions of the points:
250 self.points = []
251 for tag in self.pointTags:
252 boundingbox1 = gmsh.model.occ.getBoundingBox(0, tag)[:3]
253 boundingbox2 = gmsh.model.occ.getBoundingBox(0, tag)[3:]
254 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))
255 self.points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))
257 # Round the point positions to the nearest multiple of self.geo.dimensionTolerance to avoid
258 # numerical errors:
259 for i in range(len(self.points)):
260 for coord in range(3):
261 self.points[i][coord] = (
262 round(self.points[i][coord] / self.geo.dimensionTolerance) * self.geo.dimensionTolerance
263 )
265 if self.isCircle():
266 self.type = curveType.circle
267 # The length of the circle curves are not used.
269 elif self.isAxial():
270 self.type = curveType.axial
272 self.length = abs(self.points[0][2] - self.points[1][2])
274 elif self.isHorizontal():
275 self.type = curveType.horizontal
277 self.length = math.sqrt(
278 (self.points[0][0] - self.points[1][0]) ** 2
279 + (self.points[0][1] - self.points[1][1]) ** 2
280 )
282 else:
283 # If the curve is not axial or horizontal, it is a spiral curve:
284 self.type = curveType.spiralArc
286 # First point:
287 r1 = math.sqrt(self.points[0][0] ** 2 + self.points[0][1] ** 2)
288 theta1 = math.atan2(self.points[0][1], self.points[0][0])
290 # Second point:
291 r2 = math.sqrt(self.points[1][0] ** 2 + self.points[1][1] ** 2)
292 theta2 = math.atan2(self.points[1][1], self.points[1][0])
294 # Calculate the length of the spiral curve with numerical integration:
295 self.length = curve.calculateSpiralArcLength(r1, r2, theta1, theta2)
297 # Calculate starting turn number (n1, float) and ending turn number (n2,
298 # float): (note that they are float modulos of 1, and not the exact turn
299 # numbers)
300 self.n1 = (theta1 - self.geo.winding.theta_i) / 2 / math.pi
301 self.n1 = round(self.n1 / self.geo.winding.turnTol) * self.geo.winding.turnTol
303 self.n2 = (theta2 - self.geo.winding.theta_i) / 2 / math.pi
304 self.n2 = round(self.n2 / self.geo.winding.turnTol) * self.geo.winding.turnTol
306 def isAxial(self):
307 """
308 Checks if the curve is an axial curve. It does so by comparing the z-coordinates
309 of its starting and end points.
311 :return: True if the curve is axial, False otherwise.
312 :rtype: bool
313 """
314 return not math.isclose(
315 self.points[0][2], self.points[1][2], abs_tol=self.geo.dimensionTolerance
316 )
318 def isHorizontal(self):
319 """
320 Checks if the curve is a horizontal curve. It does so by comparing the center of
321 mass of the line and the average of the points' x and y coordinates. Having an
322 equal z-coordinate for both starting point and ending point is not enough since
323 spiral curves are on the horizontal plane as well.
325 :return: True if the curve is horizontal, False otherwise.
326 :rtype: bool
327 """
328 cm = gmsh.model.occ.getCenterOfMass(1, self.tag)
329 xcm = (self.points[0][0] + self.points[1][0]) / 2
330 ycm = (self.points[0][1] + self.points[1][1]) / 2
332 return math.isclose(cm[0], xcm, abs_tol=self.geo.dimensionTolerance) and math.isclose(
333 cm[1], ycm, abs_tol=self.geo.dimensionTolerance
334 )
336 def isCircle(self):
337 """
338 Checks if the curve is a circle. Since combined is set to True in
339 gmsh.model.getBoundary() function, the function won't return any points for
340 circle curves.
341 """
342 if len(self.points) == 0:
343 return True
344 else:
345 return False
347 @staticmethod
348 def calculateSpiralArcLength(r_1, r_2, theta_1, theta_2):
349 r"""
350 Numerically integrates the speed function of the spiral arc to calculate the
351 length of the arc.
353 In pancake coil design, spirals are cylindrical curves where the radius is
354 linearly increasing with theta. The parametric equation of a spiral sitting on
355 an XY plane can be given as follows:
357 $$
358 \\theta = t
359 $$
361 $$
362 r = a t + b
363 $$
365 $$
366 z = c
367 $$
369 where $a$, $b$, and $c$ are constants and $t$ is any real number on a given set.
371 How to calculate arc length?
373 The same spiral curve can be specified with a position vector in cylindrical
374 coordinates:
376 $$
377 \\text{position vector} = \\vec{r} = r \\vec{u}_r
378 $$
380 where $\\vec{u}_r$ is a unit vector that points towards the point.
382 Taking the derivative of the $\\vec{r}$ with respect to $t$ would give the
383 $\\text{velocity vector}$ ($\\vec{v}$) (note that both $\\vec{u}_r$ and
384 $\\vec{r}$ change with time, product rule needs to be used):
386 $$
387 \\text{velocity vector} = \\vec{\\dot{r}} = \\dot{r} \\vec{u}_r + (r \\dot{\\theta}) \\vec{u}_\\theta
388 $$
390 where $\\vec{\\dot{r}}$ and $\\dot{\\theta}$ are the derivatives of $r$ and
391 $\\theta$ with respect to $t$, and $\\vec{u}_\\theta$ is a unit vector that is
392 vertical to $\\vec{u}_r$ and points to the positive angle side.
394 The magnitude of the $\\vec{\\dot{r}}$ would result in speed. Speed's
395 integration with respect to time gives the arc length. The $\\theta$ and $r$ are
396 already specified above with the parametric equations. The only part left is
397 finding the $a$ and $b$ constants used in the parametric equations. Because TSA
398 and non-TSA spirals are a little different, the easiest way would be to
399 calculate them with a given two points on spirals, which are end and starts
400 points. The rest of the code is self-explanatory.
402 :param r_1: radial position of the starting point
403 :type r_1: float
404 :param r_2: radial position of the ending point
405 :type r_2: float
406 :param theta_1: angular position of the starting point
407 :type theta_1: float
408 :param theta_2: angular position of the ending point
409 :type theta_2: float
410 :return: length of the spiral arc
411 :rtype: float
412 """
413 # The same angle can be subtracted from both theta_1 and theta_2 to simplify the
414 # calculations:
415 theta2 = theta_2 - theta_1
416 theta1 = 0
418 # Since r = a * theta + b, r_1 = b since theta_1 = 0:
419 b = r_1
421 # Since r = a * theta + b, r_2 = a * theta2 + b:
422 a = (r_2 - b) / theta2
424 def integrand(t):
425 return math.sqrt(a**2 + (a * t + b) ** 2)
427 return abs(scipy.integrate.quad(integrand, theta1, theta2)[0])
430alreadyMeshedSurfaceTags = []
433class Mesh(Base):
434 """
435 Main mesh class for Pancake3D.
437 :param fdm: FiQuS data model
438 :param geom_folder: folder where the geometry files are saved
439 :type geom_folder: str
440 :param mesh_folder: folder where the mesh files are saved
441 :type mesh_folder: str
442 :param solution_folder: folder where the solution files are saved
443 :type solution_folder: str
444 """
446 def __init__(
447 self,
448 fdm,
449 geom_folder,
450 mesh_folder,
451 solution_folder,
452 ) -> None:
453 super().__init__(fdm, geom_folder, mesh_folder, solution_folder)
455 # Read volume information file:
456 with open(self.vi_file, "r") as f:
457 self.dimTags = json.load(f)
459 for key, value in self.dimTags.items():
460 self.dimTags[key] = [tuple(dimTag) for dimTag in value]
462 # Start GMSH:
463 self.gu = GmshUtils(self.mesh_folder)
464 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh)
466 self.contactLayerAndWindingRadialLines = [] # Store for strucured terminals
468 def generate_mesh(self):
469 """
470 Sets the mesh settings and generates the mesh.
473 """
474 logger.info("Generating Pancake3D mesh has been started.")
476 start_time = timeit.default_timer()
478 # =============================================================================
479 # MESHING WINDING AND CONTACT LAYER STARTS =======================================
480 # =============================================================================
481 allWindingAndCLSurfaceTags = []
482 allWindingAndCLLineTags = []
483 for i in range(self.geo.numberOfPancakes):
484 # Get the volume tags:
485 windingVolumeDimTags = self.dimTags[self.geo.winding.name + str(i + 1)]
486 windingVolumeTags = [dimTag[1] for dimTag in windingVolumeDimTags]
488 contactLayerVolumeDimTags = self.dimTags[self.geo.contactLayer.name + str(i + 1)]
489 contactLayerVolumeTags = [dimTag[1] for dimTag in contactLayerVolumeDimTags]
491 # Get the surface and line tags:
492 windingSurfaceTags, windingLineTags = self.get_boundaries(
493 windingVolumeDimTags, returnTags=True
494 )
495 allWindingAndCLSurfaceTags.extend(windingSurfaceTags)
496 allWindingAndCLLineTags.extend(windingLineTags)
497 contactLayerSurfaceTags, contactLayerLineTags = self.get_boundaries(
498 contactLayerVolumeDimTags, returnTags=True
499 )
500 allWindingAndCLSurfaceTags.extend(contactLayerSurfaceTags)
501 allWindingAndCLLineTags.extend(contactLayerLineTags)
503 self.structure_mesh(
504 windingVolumeTags,
505 windingSurfaceTags,
506 windingLineTags,
507 contactLayerVolumeTags,
508 contactLayerSurfaceTags,
509 contactLayerLineTags,
510 meshSettingIndex=i,
511 )
513 notchVolumesDimTags = (
514 self.dimTags["innerTransitionNotch"] + self.dimTags["outerTransitionNotch"]
515 )
516 notchVolumeTags = [dimTag[1] for dimTag in notchVolumesDimTags]
518 notchSurfaceTags, notchLineTags = self.get_boundaries(
519 notchVolumesDimTags, returnTags=True
520 )
522 for lineTag in notchLineTags:
523 if lineTag not in allWindingAndCLLineTags:
524 gmsh.model.mesh.setTransfiniteCurve(lineTag, 1)
526 recombine = self.mesh.winding.elementType[0] in ["hexahedron", "prism"]
527 for surfaceTag in notchSurfaceTags:
528 if surfaceTag not in allWindingAndCLSurfaceTags:
529 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
530 if recombine:
531 normal = gmsh.model.getNormal(surfaceTag, [0.5, 0.5])
532 if abs(normal[2]) > 1e-4:
533 pass
534 else:
535 gmsh.model.mesh.setRecombine(2, surfaceTag)
538 for volumeTag in notchVolumeTags:
539 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
541 # =============================================================================
542 # MESHING WINDING AND CONTACT LAYER ENDS =========================================
543 # =============================================================================
545 # =============================================================================
546 # MESHING AIR STARTS ==========================================================
547 # =============================================================================
548 # Winding and contact layer extrusions of the air:
549 # Get the volume tags:
550 airTopWindingExtrusionVolumeDimTags = self.dimTags[
551 self.geo.air.name + "-TopPancakeWindingExtursion"
552 ]
554 airTopContactLayerExtrusionVolumeDimTags = self.dimTags[
555 self.geo.air.name + "-TopPancakeContactLayerExtursion"
556 ]
558 airTopTerminalsExtrusionVolumeDimTags = self.dimTags[
559 self.geo.air.name + "-TopTerminalsExtrusion"
560 ]
562 airBottomWindingExtrusionVolumeDimTags = self.dimTags[
563 self.geo.air.name + "-BottomPancakeWindingExtursion"
564 ]
566 airBottomContactLayerExtrusionVolumeDimTags = self.dimTags[
567 self.geo.air.name + "-BottomPancakeContactLayerExtursion"
568 ]
570 airBottomTerminalsExtrusionVolumeDimTags = self.dimTags[
571 self.geo.air.name + "-BottomTerminalsExtrusion"
572 ]
574 removedAirVolumeDimTags = []
575 newAirVolumeDimTags = []
576 if self.mesh.air.structured:
577 # Then it means air type is cuboid!
578 airTopWindingExtrusionVolumeTags = [
579 dimTag[1] for dimTag in airTopWindingExtrusionVolumeDimTags
580 ]
581 airTopContactLayerExtrusionVolumeTags = [
582 dimTag[1] for dimTag in airTopContactLayerExtrusionVolumeDimTags
583 ]
584 airBottomWindingExtrusionVolumeTags = [
585 dimTag[1] for dimTag in airBottomWindingExtrusionVolumeDimTags
586 ]
587 airBottomContactLayerExtrusionVolumeTags = [
588 dimTag[1] for dimTag in airBottomContactLayerExtrusionVolumeDimTags
589 ]
591 # Calcualte axial number of elements for air:
592 axialElementsPerLengthForWinding = min(self.mesh.winding.axialNumberOfElements) / self.geo.winding.height
593 axneForAir = round(
594 axialElementsPerLengthForWinding * self.geo.air.axialMargin + 1e-6
595 )
597 # Get the surface and line tags:
598 (
599 airTopWindingExtrusionSurfaceTags,
600 airTopWindingExtrusionLineTags,
601 ) = self.get_boundaries(
602 airTopWindingExtrusionVolumeDimTags, returnTags=True
603 )
604 (
605 airTopContactLayerExtrusionSurfaceTags,
606 airTopContactLayerExtrusionLineTags,
607 ) = self.get_boundaries(
608 airTopContactLayerExtrusionVolumeDimTags, returnTags=True
609 )
611 self.structure_mesh(
612 airTopWindingExtrusionVolumeTags,
613 airTopWindingExtrusionSurfaceTags,
614 airTopWindingExtrusionLineTags,
615 airTopContactLayerExtrusionVolumeTags,
616 airTopContactLayerExtrusionSurfaceTags,
617 airTopContactLayerExtrusionLineTags,
618 meshSettingIndex=self.geo.numberOfPancakes - 1, # The last pancake coil
619 axialNumberOfElements=axneForAir,
620 bumpCoefficient=1,
621 )
623 # Get the surface and line tags:
624 (
625 airBottomWindingExtrusionSurfaceTags,
626 airBottomWindingExtrusionLineTags,
627 ) = self.get_boundaries(
628 airBottomWindingExtrusionVolumeDimTags, returnTags=True
629 )
630 (
631 airBottomContactLayerExtrusionSurfaceTags,
632 airBottomContactLayerExtrusionLineTags,
633 ) = self.get_boundaries(
634 airBottomContactLayerExtrusionVolumeDimTags, returnTags=True
635 )
637 self.structure_mesh(
638 airBottomWindingExtrusionVolumeTags,
639 airBottomWindingExtrusionSurfaceTags,
640 airBottomWindingExtrusionLineTags,
641 airBottomContactLayerExtrusionVolumeTags,
642 airBottomContactLayerExtrusionSurfaceTags,
643 airBottomContactLayerExtrusionLineTags,
644 meshSettingIndex=0, # The first pancake coil
645 axialNumberOfElements=axneForAir,
646 bumpCoefficient=1,
647 )
649 # Structure tubes of the air:
650 airOuterTubeVolumeDimTags = self.dimTags[self.geo.air.name + "-OuterTube"]
651 airOuterTubeVolumeTags = [dimTag[1] for dimTag in airOuterTubeVolumeDimTags]
653 airTopTubeTerminalsVolumeDimTags = self.dimTags[
654 self.geo.air.name + "-TopTubeTerminalsExtrusion"
655 ]
656 airTopTubeTerminalsVolumeTags = [
657 dimTag[1] for dimTag in airTopTubeTerminalsVolumeDimTags
658 ]
660 airBottomTubeTerminalsVolumeDimTags = self.dimTags[
661 self.geo.air.name + "-BottomTubeTerminalsExtrusion"
662 ]
663 airBottomTubeTerminalsVolumeTags = [
664 dimTag[1] for dimTag in airBottomTubeTerminalsVolumeDimTags
665 ]
667 # Structure inner cylinder of the air:
668 airInnerCylinderVolumeDimTags = self.dimTags[
669 self.geo.air.name + "-InnerCylinder"
670 ]
671 airInnerCylinderVolumeTags = [
672 dimTag[1] for dimTag in airInnerCylinderVolumeDimTags
673 ]
675 airTubesAndCylinders = airOuterTubeVolumeTags + airInnerCylinderVolumeTags
677 if self.geo.air.shellTransformation:
678 shellVolumes = self.dimTags[self.geo.air.shellVolumeName]
679 shellVolumeTags = [dimTag[1] for dimTag in shellVolumes]
680 airTubesAndCylinders.extend(shellVolumeTags)
682 airRadialElementMultiplier = 1 / self.mesh.air.radialElementSize
683 self.structure_tubes_and_cylinders(
684 airTubesAndCylinders,
685 radialElementMultiplier=airRadialElementMultiplier,
686 )
688 if self.mesh.terminals.structured:
689 terminalsRadialElementMultiplier = 1 / self.mesh.terminals.radialElementSize
691 self.structure_tubes_and_cylinders(
692 airTopTubeTerminalsVolumeTags + airBottomTubeTerminalsVolumeTags,
693 radialElementMultiplier=terminalsRadialElementMultiplier,
694 )
696 airTopTouchingTerminalsVolumeDimTags = list(
697 set(airTopTerminalsExtrusionVolumeDimTags)
698 - set(airTopTubeTerminalsVolumeDimTags)
699 )
700 airTopTouchingTerminalsVolumeTags = [
701 dimTag[1] for dimTag in airTopTouchingTerminalsVolumeDimTags
702 ]
704 airBottomTouchingTerminalsVolumeDimTags = list(
705 set(airBottomTerminalsExtrusionVolumeDimTags)
706 - set(airBottomTubeTerminalsVolumeDimTags)
707 )
708 airBottomTouchingTerminalsVolumeTags = [
709 dimTag[1] for dimTag in airBottomTouchingTerminalsVolumeDimTags
710 ]
712 self.structure_tubes_and_cylinders(
713 airTopTouchingTerminalsVolumeTags
714 + airBottomTouchingTerminalsVolumeTags,
715 terminalNonTubeParts=True,
716 radialElementMultiplier=terminalsRadialElementMultiplier,
717 )
719 else:
720 # Fuse top volumes:
721 airTopVolumeDimTags = (
722 airTopWindingExtrusionVolumeDimTags
723 + airTopContactLayerExtrusionVolumeDimTags
724 + airTopTerminalsExtrusionVolumeDimTags
725 )
726 airTopVolumeDimTag = Mesh.fuse_volumes(
727 airTopVolumeDimTags,
728 fuseSurfaces=True,
729 fusedSurfacesArePlane=True,
730 )
731 newAirVolumeDimTags.append(airTopVolumeDimTag)
732 removedAirVolumeDimTags.extend(airTopVolumeDimTags)
734 # Fuse bottom volumes:
735 airBottomVolumeDimTags = (
736 airBottomWindingExtrusionVolumeDimTags
737 + airBottomContactLayerExtrusionVolumeDimTags
738 + airBottomTerminalsExtrusionVolumeDimTags
739 )
740 airBottomVolumeDimTag = Mesh.fuse_volumes(
741 airBottomVolumeDimTags,
742 fuseSurfaces=True,
743 fusedSurfacesArePlane=True,
744 )
745 newAirVolumeDimTags.append(airBottomVolumeDimTag)
746 removedAirVolumeDimTags.extend(airBottomVolumeDimTags)
748 # Fuse inner cylinder and outer tube part of air:
749 airInnerCylinderVolumeDimTags = self.dimTags[
750 self.geo.air.name + "-InnerCylinder"
751 ]
752 if self.geo.numberOfPancakes > 1:
753 # Fuse the first two and the last two volumes separately (due to cuts):
754 firstTwoVolumes = airInnerCylinderVolumeDimTags[0:2]
755 lastTwoVolumes = airInnerCylinderVolumeDimTags[-2:]
756 airInnerCylinderVolumeDimTags = airInnerCylinderVolumeDimTags[2:-2]
757 airInnerCylinderVolumeDimTag = Mesh.fuse_volumes(
758 airInnerCylinderVolumeDimTags, fuseSurfaces=False
759 )
760 airInnerCylinderVolumeDimTagFirst = Mesh.fuse_volumes(
761 firstTwoVolumes,
762 fuseSurfaces=False,
763 )
764 airInnerCylinderVolumeDimTagLast = Mesh.fuse_volumes(
765 lastTwoVolumes,
766 fuseSurfaces=False,
767 )
768 newAirVolumeDimTags.append(airInnerCylinderVolumeDimTag)
769 newAirVolumeDimTags.append(airInnerCylinderVolumeDimTagFirst)
770 newAirVolumeDimTags.append(airInnerCylinderVolumeDimTagLast)
771 removedAirVolumeDimTags.extend(
772 airInnerCylinderVolumeDimTags + firstTwoVolumes + lastTwoVolumes
773 )
774 self.dimTags[self.geo.air.name + "-InnerCylinder"] = [
775 airInnerCylinderVolumeDimTag,
776 airInnerCylinderVolumeDimTagFirst,
777 airInnerCylinderVolumeDimTagLast,
778 ]
779 else:
780 pass
781 # self.dimTags[self.geo.air.name + "-InnerCylinder"] = [
782 # self.dimTags[self.geo.air.name + "-InnerCylinder"][1],
783 # self.dimTags[self.geo.air.name + "-InnerCylinder"][0],
784 # self.dimTags[self.geo.air.name + "-InnerCylinder"][2],
785 # ]
787 airOuterTubeVolumeDimTags = self.dimTags[self.geo.air.name + "-OuterTube"]
788 airOuterTubeVolumeDimTag = Mesh.fuse_volumes(
789 airOuterTubeVolumeDimTags,
790 fuseSurfaces=True,
791 fusedSurfacesArePlane=False,
792 )
793 newAirOuterTubeVolumeDimTag = airOuterTubeVolumeDimTag
794 removedAirVolumeDimTags.extend(airOuterTubeVolumeDimTags)
795 self.dimTags[self.geo.air.name + "-OuterTube"] = [newAirOuterTubeVolumeDimTag]
797 if self.geo.air.shellTransformation:
798 # Fuse air shell volumes:
799 if self.geo.air.type == "cylinder":
800 removedShellVolumeDimTags = []
801 shellVolumeDimTags = self.dimTags[self.geo.air.shellVolumeName]
802 shellVolumeDimTag = Mesh.fuse_volumes(
803 shellVolumeDimTags,
804 fuseSurfaces=True,
805 fusedSurfacesArePlane=False,
806 )
807 removedShellVolumeDimTags.extend(shellVolumeDimTags)
808 newShellVolumeDimTags = [shellVolumeDimTag]
809 for removedDimTag in removedShellVolumeDimTags:
810 self.dimTags[self.geo.air.shellVolumeName].remove(removedDimTag)
811 elif self.geo.air.type == "cuboid":
812 # Unfortunately, surfaces cannot be combined for the cuboid type of air.
813 # However, it doesn't affect the mesh quality that much.
814 newShellVolumeDimTags = []
816 shellPart1VolumeDimTag = Mesh.fuse_volumes(
817 self.dimTags[self.geo.air.shellVolumeName + "-Part1"],
818 fuseSurfaces=False,
819 )
820 self.dimTags[self.geo.air.shellVolumeName + "-Part1"] = [
821 shellPart1VolumeDimTag
822 ]
824 shellPart2VolumeDimTag = Mesh.fuse_volumes(
825 self.dimTags[self.geo.air.shellVolumeName + "-Part2"],
826 fuseSurfaces=False,
827 )
828 self.dimTags[self.geo.air.shellVolumeName + "-Part2"] = [
829 shellPart2VolumeDimTag
830 ]
832 shellPart3VolumeDimTag = Mesh.fuse_volumes(
833 self.dimTags[self.geo.air.shellVolumeName + "-Part3"],
834 fuseSurfaces=False,
835 )
836 self.dimTags[self.geo.air.shellVolumeName + "-Part3"] = [
837 shellPart3VolumeDimTag
838 ]
840 shellPart4VolumeDimTag = Mesh.fuse_volumes(
841 self.dimTags[self.geo.air.shellVolumeName + "-Part4"],
842 fuseSurfaces=False,
843 )
844 self.dimTags[self.geo.air.shellVolumeName + "-Part4"] = [
845 shellPart4VolumeDimTag
846 ]
848 # The problem is, shell volume and outer air tube volume has a common
849 # surface and that surface should be combined as well for high quality mesh.
850 # However, it can be only done for cylinder type of air for now.
851 # Get the combined boundary surfaces:
852 if self.geo.air.type == "cylinder":
853 (
854 newAirOuterTubeVolumeDimTag,
855 newShellVolumeDimTag,
856 ) = Mesh.fuse_common_surfaces_of_two_volumes(
857 [airOuterTubeVolumeDimTag],
858 newShellVolumeDimTags,
859 fuseOtherSurfaces=False,
860 surfacesArePlane=False,
861 )
862 self.dimTags[self.geo.air.name + "-OuterTube"] = [newAirOuterTubeVolumeDimTag]
864 airOuterTubeVolumeDimTag = newAirOuterTubeVolumeDimTag
865 self.dimTags[self.geo.air.shellVolumeName].append(
866 newShellVolumeDimTag
867 )
869 newAirVolumeDimTags.append(newAirOuterTubeVolumeDimTag)
871 # Update volume tags dictionary of air:
872 self.dimTags[self.geo.air.name] = list(
873 (
874 set(self.dimTags[self.geo.air.name]) - set(removedAirVolumeDimTags)
875 ).union(set(newAirVolumeDimTags))
876 )
878 # ==============================================================================
879 # MESHING AIR ENDS =============================================================
880 # ==============================================================================
882 # ==============================================================================
883 # MESHING TERMINALS STARTS =====================================================
884 # ==============================================================================
885 if self.mesh.terminals.structured:
886 # Structure tubes of the terminals:
887 terminalOuterTubeVolumeDimTags = self.dimTags[self.geo.terminals.outer.name + "-Tube"]
888 terminalOuterTubeVolumeTags = [
889 dimTag[1] for dimTag in terminalOuterTubeVolumeDimTags
890 ]
891 terminalInnerTubeVolumeDimTags = self.dimTags[self.geo.terminals.inner.name + "-Tube"]
892 terminalInnerTubeVolumeTags = [
893 dimTag[1] for dimTag in terminalInnerTubeVolumeDimTags
894 ]
896 terminalsRadialElementMultiplier = 1 / self.mesh.terminals.radialElementSize
897 self.structure_tubes_and_cylinders(
898 terminalOuterTubeVolumeTags + terminalInnerTubeVolumeTags,
899 radialElementMultiplier=terminalsRadialElementMultiplier,
900 )
902 # Structure nontube parts of the terminals:
903 terminalOuterNonTubeVolumeDimTags = self.dimTags[
904 self.geo.terminals.outer.name + "-Touching"
905 ]
906 terminalOuterNonTubeVolumeTags = [
907 dimTag[1] for dimTag in terminalOuterNonTubeVolumeDimTags
908 ]
909 terminalInnerNonTubeVolumeDimTags = self.dimTags[
910 self.geo.terminals.inner.name + "-Touching"
911 ]
912 terminalInnerNonTubeVolumeTags = [
913 dimTag[1] for dimTag in terminalInnerNonTubeVolumeDimTags
914 ]
916 self.structure_tubes_and_cylinders(
917 terminalInnerNonTubeVolumeTags + terminalOuterNonTubeVolumeTags,
918 terminalNonTubeParts=True,
919 radialElementMultiplier=terminalsRadialElementMultiplier,
920 )
921 # ==============================================================================
922 # MESHING TERMINALS ENDS =======================================================
923 # ==============================================================================
925 # ==============================================================================
926 # FIELD SETTINGS STARTS ========================================================
927 # ==============================================================================
929 # Mesh fields for the air:
930 # Meshes will grow as they get further from the field surfaces:
931 fieldSurfacesDimTags = gmsh.model.getBoundary(
932 self.dimTags[self.geo.winding.name], oriented=False, combined=True
933 )
934 fieldSurfacesTags = [dimTag[1] for dimTag in fieldSurfacesDimTags]
936 distanceField = gmsh.model.mesh.field.add("Distance")
938 gmsh.model.mesh.field.setNumbers(
939 distanceField,
940 "SurfacesList",
941 fieldSurfacesTags,
942 )
944 thresholdField = gmsh.model.mesh.field.add("Threshold")
945 gmsh.model.mesh.field.setNumber(thresholdField, "InField", distanceField)
946 gmsh.model.mesh.field.setNumber(thresholdField, "SizeMin", self.mesh.sizeMin)
947 gmsh.model.mesh.field.setNumber(thresholdField, "SizeMax", self.mesh.sizeMax)
948 gmsh.model.mesh.field.setNumber(
949 thresholdField, "DistMin", self.mesh.startGrowingDistance
950 )
952 gmsh.model.mesh.field.setNumber(
953 thresholdField, "DistMax", self.mesh.stopGrowingDistance
954 )
956 gmsh.model.mesh.field.setAsBackgroundMesh(thresholdField)
958 # ==============================================================================
959 # FIELD SETTINGS ENDS ==========================================================
960 # ==============================================================================
962 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
963 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
964 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
966 try:
967 # Only print warnings and errors:
968 # Don't print on terminal, because we will use logger:
969 gmsh.option.setNumber("General.Terminal", 0)
970 # Start logger:
971 gmsh.logger.start()
973 gmsh.option.setNumber("Mesh.Algorithm", 6)
974 gmsh.option.setNumber("Mesh.Algorithm3D", 1)
976 # Mesh:
977 gmsh.model.mesh.generate()
978 gmsh.model.mesh.optimize("Netgen")
980 # Print the log:
981 log = gmsh.logger.get()
982 for line in log:
983 if line.startswith("Info"):
984 logger.info(re.sub(r"Info:\s+", "", line))
985 elif line.startswith("Warning"):
986 logger.warning(re.sub(r"Warning:\s+", "", line))
988 gmsh.logger.stop()
989 except:
990 # Print the log:
991 log = gmsh.logger.get()
992 for line in log:
993 if line.startswith("Info"):
994 logger.info(re.sub(r"Info:\s+", "", line))
995 elif line.startswith("Warning"):
996 logger.warning(re.sub(r"Warning:\s+", "", line))
997 elif line.startswith("Error"):
998 logger.error(re.sub(r"Error:\s+", "", line))
1000 gmsh.logger.stop()
1002 self.generate_regions()
1004 logger.error(
1005 "Meshing Pancake3D magnet has failed. Try to change"
1006 " minimumElementSize and maximumElementSize parameters."
1007 )
1008 raise
1010 # SICN not implemented in 1D!
1011 allElementsDim2 = list(gmsh.model.mesh.getElements(dim=2)[1][0])
1012 allElementsDim3 = list(gmsh.model.mesh.getElements(dim=3)[1][0])
1013 allElements = allElementsDim2 + allElementsDim3
1014 elementQualities = gmsh.model.mesh.getElementQualities(allElements)
1015 lowestQuality = min(elementQualities)
1016 averageQuality = sum(elementQualities) / len(elementQualities)
1017 NofLowQualityElements = len(
1018 [quality for quality in elementQualities if quality < 0.01]
1019 )
1020 NofIllElemets = len(
1021 [quality for quality in elementQualities if quality < 0.001]
1022 )
1024 logger.info(
1025 f"The lowest quality among the elements is {lowestQuality:.4f} (SICN). The"
1026 " number of elements with quality lower than 0.01 is"
1027 f" {NofLowQualityElements}."
1028 )
1030 if NofIllElemets > 0:
1031 logger.warning(
1032 f"There are {NofIllElemets} elements with quality lower than 0.001. Try"
1033 " to change minimumElementSize and maximumElementSize parameters."
1034 )
1036 # Create cuts:
1037 # This is required to make the air a simply connected domain. This is required
1038 # for the solution part. You can read more about Homology in GMSH documentation.
1039 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name]]
1041 if self.geo.air.shellTransformation:
1042 shellTags = [
1043 dimTag[1] for dimTag in self.dimTags[self.geo.air.shellVolumeName]
1044 ]
1045 airTags.extend(shellTags)
1047 dummyAirRegion = gmsh.model.addPhysicalGroup(dim=3, tags=airTags)
1048 dummyAirRegionDimTag = (3, dummyAirRegion)
1050 innerCylinderTags = [self.dimTags[self.geo.air.name + "-InnerCylinder"][0][1]]
1051 gapTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name + "-Gap"]]
1052 # Only remove every second gap:
1053 gapTags = gapTags[1::2]
1055 dummyAirRegionWithoutInnerCylinder = gmsh.model.addPhysicalGroup(
1056 dim=3, tags=list(set(airTags) - set(innerCylinderTags) - set(gapTags))
1057 )
1058 dummyAirRegionWithoutInnerCylinderDimTag = (
1059 3,
1060 dummyAirRegionWithoutInnerCylinder,
1061 )
1063 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.winding.name]]
1064 dummyWindingRegion = gmsh.model.addPhysicalGroup(dim=3, tags=windingTags)
1065 dummyWindingRegionDimTag = (3, dummyWindingRegion)
1067 if self.geo.contactLayer.thinShellApproximation:
1068 # Find all the contact layer surfaces:
1069 allWindingDimTags = []
1070 for i in range(self.geo.numberOfPancakes):
1071 windingDimTags = self.dimTags[self.geo.winding.name + str(i + 1)]
1072 allWindingDimTags.extend(windingDimTags)
1074 windingBoundarySurfaces = gmsh.model.getBoundary(
1075 allWindingDimTags, combined=True, oriented=False
1076 )
1077 allWindingSurfaces = gmsh.model.getBoundary(
1078 allWindingDimTags, combined=False, oriented=False
1079 )
1081 contactLayerSurfacesDimTags = list(
1082 set(allWindingSurfaces) - set(windingBoundarySurfaces)
1083 )
1084 contactLayerTags = [dimTag[1] for dimTag in contactLayerSurfacesDimTags]
1086 # Get rid of non-contactLayer surfaces:
1087 realContactLayerTags = []
1088 for contactLayerTag in contactLayerTags:
1089 surfaceNormal = list(gmsh.model.getNormal(contactLayerTag, [0.5, 0.5]))
1090 centerOfMass = gmsh.model.occ.getCenterOfMass(2, contactLayerTag)
1092 if (
1093 abs(
1094 surfaceNormal[0] * centerOfMass[0]
1095 + surfaceNormal[1] * centerOfMass[1]
1096 )
1097 > 1e-6
1098 ):
1099 realContactLayerTags.append(contactLayerTag)
1101 # Get rid of surfaces that touch terminals:
1102 terminalSurfaces = gmsh.model.getBoundary(
1103 self.dimTags[self.geo.terminals.outer.name] + self.dimTags[self.geo.terminals.inner.name],
1104 combined=False,
1105 oriented=False,
1106 )
1107 terminalSurfaces = [dimTag[1] for dimTag in terminalSurfaces]
1108 finalContactLayerTags = [
1109 tag for tag in realContactLayerTags if tag not in terminalSurfaces
1110 ]
1112 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1113 dim=2, tags=finalContactLayerTags
1114 )
1115 dummyContactLayerRegionDimTag = (2, dummyContactLayerRegion)
1117 else:
1118 contactLayerTags = [dimTag[1] for dimTag in self.dimTags[self.geo.contactLayer.name]]
1120 # get rid of volumes that touch terminals:
1121 terminalSurfaces = gmsh.model.getBoundary(
1122 self.dimTags[self.geo.terminals.outer.name] + self.dimTags[self.geo.terminals.inner.name],
1123 combined=False,
1124 oriented=False,
1125 )
1126 finalContactLayerTags = []
1127 for contactLayerTag in contactLayerTags:
1128 insulatorSurfaces = gmsh.model.getBoundary(
1129 [(3, contactLayerTag)], combined=False, oriented=False
1130 )
1131 itTouchesTerminals = False
1132 for insulatorSurface in insulatorSurfaces:
1133 if insulatorSurface in terminalSurfaces:
1134 itTouchesTerminals = True
1135 break
1137 if not itTouchesTerminals:
1138 finalContactLayerTags.append(contactLayerTag)
1140 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1141 dim=3, tags=finalContactLayerTags
1142 )
1143 dummyContactLayerRegionDimTag = (3, dummyContactLayerRegion)
1145 # First cohomology request (normal cut for NI coils):
1146 gmsh.model.mesh.addHomologyRequest(
1147 "Cohomology",
1148 domainTags=[dummyAirRegion],
1149 subdomainTags=[],
1150 dims=[1],
1151 )
1153 if self.mesh.computeCohomologyForInsulating:
1154 # Second cohomology request (insulated cut for insulated coils):
1155 if self.geo.numberOfPancakes > 1:
1156 gmsh.model.mesh.addHomologyRequest(
1157 "Cohomology",
1158 domainTags=[
1159 dummyAirRegionWithoutInnerCylinder,
1160 dummyContactLayerRegion,
1161 ],
1162 subdomainTags=[],
1163 dims=[1],
1164 )
1165 else:
1166 gmsh.model.mesh.addHomologyRequest(
1167 "Cohomology",
1168 domainTags=[
1169 dummyAirRegion,
1170 dummyContactLayerRegion,
1171 ],
1172 subdomainTags=[],
1173 dims=[1],
1174 )
1176 # Third cohomology request (for cuts between pancake coils):
1177 gmsh.model.mesh.addHomologyRequest(
1178 "Cohomology",
1179 domainTags=[
1180 dummyAirRegion,
1181 dummyContactLayerRegion,
1182 dummyWindingRegion,
1183 ],
1184 subdomainTags=[],
1185 dims=[1],
1186 )
1188 # Start logger:
1189 gmsh.logger.start()
1191 cuts = gmsh.model.mesh.computeHomology()
1193 # Print the log:
1194 log = gmsh.logger.get()
1195 for line in log:
1196 if line.startswith("Info"):
1197 logger.info(re.sub(r"Info:\s+", "", line))
1198 elif line.startswith("Warning"):
1199 logger.warning(re.sub(r"Warning:\s+", "", line))
1200 gmsh.logger.stop()
1202 if self.geo.numberOfPancakes > 1:
1203 cutsDictionary = {
1204 "H^1{1}": self.geo.air.cutName,
1205 "H^1{1,4,3}": "CutsBetweenPancakes",
1206 "H^1{2,4}": "CutsForPerfectInsulation",
1207 }
1208 else:
1209 cutsDictionary = {
1210 "H^1{1}": self.geo.air.cutName,
1211 "H^1{1,4,3}": "CutsBetweenPancakes",
1212 "H^1{1,4}": "CutsForPerfectInsulation",
1213 }
1214 cutTags = [dimTag[1] for dimTag in cuts]
1215 cutEntities = []
1216 for tag in cutTags:
1217 name = gmsh.model.getPhysicalName(1, tag)
1218 cutEntities = list(gmsh.model.getEntitiesForPhysicalGroup(1, tag))
1219 cutEntitiesDimTags = [(1, cutEntity) for cutEntity in cutEntities]
1220 for key in cutsDictionary:
1221 if key in name:
1222 if cutsDictionary[key] in self.dimTags:
1223 self.dimTags[cutsDictionary[key]].extend(cutEntitiesDimTags)
1224 else:
1225 self.dimTags[cutsDictionary[key]] = cutEntitiesDimTags
1227 # Remove newly created physical groups because they will be created again in
1228 # generate_regions method.
1229 gmsh.model.removePhysicalGroups(
1230 [dummyContactLayerRegionDimTag]
1231 + [dummyAirRegionDimTag]
1232 + [dummyAirRegionWithoutInnerCylinderDimTag]
1233 + [dummyWindingRegionDimTag]
1234 + cuts
1235 )
1237 logger.info(
1238 "Generating Pancake3D mesh has been finished in"
1239 f" {timeit.default_timer() - start_time:.2f} s."
1240 )
1242 def structure_mesh(
1243 self,
1244 windingVolumeTags,
1245 windingSurfaceTags,
1246 windingLineTags,
1247 contactLayerVolumeTags,
1248 contactLayerSurfaceTags,
1249 contactLayerLineTags,
1250 meshSettingIndex,
1251 axialNumberOfElements=None,
1252 bumpCoefficient=None,
1253 ):
1254 """
1255 Structures the winding and contact layer meshed depending on the user inputs. If
1256 the bottom and top part of the air is to be structured, the same method is used.
1258 :param windingVolumeTags: tags of the winding volumes
1259 :type windingVolumeTags: list[int]
1260 :param windingSurfaceTags: tags of the winding surfaces
1261 :type windingSurfaceTags: list[int]
1262 :param windingLineTags: tags of the winding lines
1263 :type windingLineTags: list[int]
1264 :param contactLayerVolumeTags: tags of the contact layer volumes
1265 :type contactLayerVolumeTags: list[int]
1266 :param contactLayerSurfaceTags: tags of the contact layer surfaces
1267 :type contactLayerSurfaceTags: list[int]
1268 :param contactLayerLineTags: tags of the contact layer lines
1269 :type contactLayerLineTags: list[int]
1270 :param meshSettingIndex: index of the mesh setting
1271 :type meshSettingIndex: int
1272 :param axialNumberOfElements: number of axial elements
1273 :type axialNumberOfElements: int, optional
1274 :param bumpCoefficient: bump coefficient for axial meshing
1275 :type bumpCoefficient: float, optional
1277 """
1278 # Transfinite settings:
1279 # Arc lenght of the innermost one turn of spiral:
1280 if self.geo.contactLayer.thinShellApproximation:
1281 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1282 self.geo.winding.innerRadius,
1283 self.geo.winding.innerRadius
1284 + self.geo.winding.thickness
1285 + self.geo.contactLayer.thickness * (self.geo.numberOfPancakes - 1) / self.geo.numberOfPancakes,
1286 0,
1287 2 * math.pi,
1288 )
1289 else:
1290 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1291 self.geo.winding.innerRadius,
1292 self.geo.winding.innerRadius + self.geo.winding.thickness,
1293 0,
1294 2 * math.pi,
1295 )
1297 # Arc length of one element:
1298 arcElementLength = oneTurnSpiralLength / self.mesh.winding.azimuthalNumberOfElementsPerTurn[meshSettingIndex]
1300 # Number of azimuthal elements per turn:
1301 arcNumElementsPerTurn = round(oneTurnSpiralLength / arcElementLength)
1303 # Make all the lines transfinite:
1304 for j, lineTags in enumerate([windingLineTags, contactLayerLineTags]):
1305 for lineTag in lineTags:
1306 lineObject = curve(lineTag, self.geo)
1308 if lineObject.type is curveType.horizontal:
1309 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1310 # used.
1311 if self.geo.contactLayer.thinShellApproximation:
1312 numNodes = self.mesh.winding.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1314 else:
1315 if j == 0:
1316 # This line is the winding's horizontal line:
1317 numNodes = self.mesh.winding.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1319 else:
1320 # This line is the contact layer's horizontal line:
1321 numNodes = self.mesh.contactLayer.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1323 # Set transfinite curve:
1324 self.contactLayerAndWindingRadialLines.append(lineTag)
1325 gmsh.model.mesh.setTransfiniteCurve(lineTag, numNodes)
1327 elif lineObject.type is curveType.axial:
1328 # The curve is axial, so axialNumberOfElements entry is used.
1329 if axialNumberOfElements is None:
1330 numNodes = self.mesh.winding.axialNumberOfElements[meshSettingIndex] + 1
1331 else:
1332 numNodes = axialNumberOfElements + 1
1334 if bumpCoefficient is None:
1335 bumpCoefficient = self.mesh.winding.axialDistributionCoefficient[meshSettingIndex]
1336 gmsh.model.mesh.setTransfiniteCurve(
1337 lineTag, numNodes, meshType="Bump", coef=bumpCoefficient
1338 )
1340 else:
1341 # The line is an arc, so the previously calculated arcNumElementsPerTurn
1342 # is used. All the number of elements per turn must be the same
1343 # independent of radial position. Otherwise, transfinite meshing cannot
1344 # be performed. However, to support the float number of turns, number
1345 # of nodes are being calculated depending on the start and end turns of
1346 # the arc.d
1347 lengthInTurns = abs(lineObject.n2 - lineObject.n1)
1348 if lengthInTurns > 0.5:
1349 # The arc can never be longer than half a turn.
1350 lengthInTurns = 1 - lengthInTurns
1352 lengthInTurns = (
1353 round(lengthInTurns / self.geo.winding.turnTol) * self.geo.winding.turnTol
1354 )
1356 arcNumEl = round(arcNumElementsPerTurn * lengthInTurns)
1358 arcNumNodes = int(arcNumEl + 1)
1360 # Set transfinite curve:
1361 gmsh.model.mesh.setTransfiniteCurve(lineTag, arcNumNodes)
1363 for j, surfTags in enumerate([windingSurfaceTags, contactLayerSurfaceTags]):
1364 for surfTag in surfTags:
1365 # Make all the surfaces transfinite:
1366 gmsh.model.mesh.setTransfiniteSurface(surfTag)
1368 if self.mesh.winding.elementType[meshSettingIndex] == "hexahedron":
1369 # If the element type is hexahedron, recombine all the surfaces:
1370 gmsh.model.mesh.setRecombine(2, surfTag)
1371 elif self.mesh.winding.elementType[meshSettingIndex] == "prism":
1372 # If the element type is prism, recombine only the side surfaces:
1373 surfaceNormal = list(gmsh.model.getNormal(surfTag, [0.5, 0.5]))
1374 if abs(surfaceNormal[2]) < 1e-6:
1375 gmsh.model.mesh.setRecombine(2, surfTag)
1377 # If the element type is tetrahedron, do not recombine any surface.
1379 for volTag in windingVolumeTags + contactLayerVolumeTags:
1380 # Make all the volumes transfinite:
1381 gmsh.model.mesh.setTransfiniteVolume(volTag)
1383 def structure_tubes_and_cylinders(
1384 self, volumeTags, terminalNonTubeParts=False, radialElementMultiplier=1
1385 ):
1386 # Number of azimuthal elements per quarter:
1387 arcNumElementsPerQuarter = int(self.mesh.winding.azimuthalNumberOfElementsPerTurn[0] / 4)
1388 radialNumberOfElementsPerLength = (
1389 self.mesh.winding.radialNumberOfElementsPerTurn[0] / self.geo.winding.thickness * radialElementMultiplier
1390 )
1392 surfacesDimTags = gmsh.model.getBoundary(
1393 [(3, tag) for tag in volumeTags], combined=False, oriented=False
1394 )
1395 surfacesTags = [dimTag[1] for dimTag in surfacesDimTags]
1396 surfacesTags = list(set(surfacesTags))
1398 curvesDimTags = gmsh.model.getBoundary(
1399 surfacesDimTags, combined=False, oriented=False
1400 )
1401 curvesTags = [dimTag[1] for dimTag in curvesDimTags]
1403 # Make all the lines transfinite:
1404 for curveTag in curvesTags:
1405 curveObject = curve(curveTag, self.geo)
1407 if curveObject.type is curveType.horizontal:
1408 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1409 # used.
1411 # But, the curve might be a part of the transitionNotch.
1412 isTransitionNotch = False
1413 point2 = curveObject.points[1]
1414 point1 = curveObject.points[0]
1415 if (
1416 abs(point2[0] - point1[0]) > 1e-5
1417 and abs(point2[1] - point1[1]) > 1e-5
1418 ):
1419 isTransitionNotch = True
1421 if isTransitionNotch:
1422 gmsh.model.mesh.setTransfiniteCurve(curveTag, 2)
1423 else:
1424 if terminalNonTubeParts:
1425 if curveTag not in self.contactLayerAndWindingRadialLines:
1426 numNodes = (
1427 round(radialNumberOfElementsPerLength * self.geo.terminals.inner.thickness)
1428 + 1
1429 )
1430 # Set transfinite curve:
1431 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1432 else:
1433 numNodes = (
1434 round(radialNumberOfElementsPerLength * curveObject.length)
1435 + 1
1436 )
1437 # Set transfinite curve:
1438 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1440 elif curveObject.type is curveType.axial:
1441 # The curve is axial, so axialNumberOfElements entry is used.
1442 if math.isclose(curveObject.length, self.geo.winding.height, rel_tol=1e-7):
1443 numNodes = min(self.mesh.winding.axialNumberOfElements) + 1
1444 else:
1445 axialElementsPerLength = min(self.mesh.winding.axialNumberOfElements) / self.geo.winding.height
1446 numNodes = (
1447 round(axialElementsPerLength * curveObject.length + 1e-6) + 1
1448 )
1450 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1452 else:
1453 # The line is an arc
1454 lengthInTurns = abs(curveObject.n2 - curveObject.n1)
1455 if lengthInTurns > 0.5:
1456 # The arc can never be longer than half a turn.
1457 lengthInTurns = 1 - lengthInTurns
1459 lengthInTurns = (
1460 round(lengthInTurns / self.geo.winding.turnTol) * self.geo.winding.turnTol
1461 )
1463 arcNumEl = round(arcNumElementsPerQuarter * 4 * lengthInTurns)
1465 arcNumNodes = int(arcNumEl + 1)
1467 # Set transfinite curve:
1469 gmsh.model.mesh.setTransfiniteCurve(curveTag, arcNumNodes)
1471 for surfaceTag in surfacesTags:
1472 # Make all the surfaces transfinite:
1474 if self.mesh.winding.elementType[0] == "hexahedron":
1475 # If the element type is hexahedron, recombine all the surfaces:
1476 gmsh.model.mesh.setRecombine(2, surfaceTag)
1477 elif self.mesh.winding.elementType[0] == "prism":
1478 # If the element type is prism, recombine only the side surfaces:
1479 surfaceNormal = list(gmsh.model.getNormal(surfaceTag, [0.5, 0.5]))
1480 if abs(surfaceNormal[2]) < 1e-5:
1481 gmsh.model.mesh.setRecombine(2, surfaceTag)
1483 curves = gmsh.model.getBoundary(
1484 [(2, surfaceTag)], combined=False, oriented=False
1485 )
1486 numberOfCurves = len(curves)
1487 if terminalNonTubeParts:
1488 if numberOfCurves == 4:
1489 numberOfHorizontalCurves = 0
1490 for curveTag in curves:
1491 curveObject = curve(curveTag[1], self.geo)
1492 if curveObject.type is curveType.horizontal:
1493 numberOfHorizontalCurves += 1
1495 if numberOfHorizontalCurves == 3:
1496 pass
1497 else:
1498 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1500 elif numberOfCurves == 3:
1501 pass
1502 else:
1503 curves = gmsh.model.getBoundary(
1504 [(2, surfaceTag)], combined=False, oriented=False
1505 )
1506 curveObjects = [curve(line[1], self.geo) for line in curves]
1508 divisionCurves = []
1509 for curveObject in curveObjects:
1510 if curveObject.type is curveType.horizontal:
1511 point1 = curveObject.points[0]
1512 point2 = curveObject.points[1]
1513 if not (
1514 abs(point2[0] - point1[0]) > 1e-5
1515 and abs(point2[1] - point1[1]) > 1e-5
1516 ):
1517 divisionCurves.append(curveObject)
1519 cornerPoints = (
1520 divisionCurves[0].pointTags + divisionCurves[1].pointTags
1521 )
1523 if surfaceTag not in alreadyMeshedSurfaceTags:
1524 alreadyMeshedSurfaceTags.append(surfaceTag)
1525 gmsh.model.mesh.setTransfiniteSurface(
1526 surfaceTag, cornerTags=cornerPoints
1527 )
1528 else:
1529 if numberOfCurves == 3:
1530 # Then it is a pie, corner points should be adjusted:
1531 originPointTag = None
1532 curveObject1 = curve(curves[0][1], self.geo)
1533 for point, tag in zip(curveObject1.points, curveObject1.pointTags):
1534 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1535 originPointTag = tag
1537 if originPointTag is None:
1538 curveObject2 = curve(curves[1][1], self.geo)
1539 for point, tag in zip(
1540 curveObject2.points, curveObject2.pointTags
1541 ):
1542 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1543 originPointTag = tag
1545 otherPointDimTags = gmsh.model.getBoundary(
1546 [(2, surfaceTag)],
1547 combined=False,
1548 oriented=False,
1549 recursive=True,
1550 )
1551 otherPointTags = [dimTag[1] for dimTag in otherPointDimTags]
1552 otherPointTags.remove(originPointTag)
1553 cornerTags = [originPointTag] + otherPointTags
1554 gmsh.model.mesh.setTransfiniteSurface(
1555 surfaceTag, cornerTags=cornerTags
1556 )
1557 else:
1558 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1560 for volumeTag in volumeTags:
1561 if terminalNonTubeParts:
1562 surfaces = gmsh.model.getBoundary(
1563 [(3, volumeTag)], combined=False, oriented=False
1564 )
1565 curves = gmsh.model.getBoundary(
1566 surfaces, combined=False, oriented=False
1567 )
1568 curves = list(set(curves))
1570 if len(curves) == 12:
1571 numberOfArcs = 0
1572 for curveTag in curves:
1573 curveObject = curve(curveTag[1], self.geo)
1574 if curveObject.type is curveType.spiralArc:
1575 numberOfArcs += 1
1576 if numberOfArcs == 2:
1577 pass
1578 else:
1579 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1580 # elif len(curves) == 15:
1581 # divisionCurves = []
1582 # for curveTag in curves:
1583 # curveObject = curve(curveTag[1], self.geo)
1584 # if curveObject.type is curveType.horizontal:
1585 # point1 = curveObject.points[0]
1586 # point2 = curveObject.points[1]
1587 # if not (
1588 # abs(point2[0] - point1[0]) > 1e-5
1589 # and abs(point2[1] - point1[1]) > 1e-5
1590 # ):
1591 # divisionCurves.append(curveObject)
1593 # cornerPoints = (
1594 # divisionCurves[0].pointTags
1595 # + divisionCurves[1].pointTags
1596 # + divisionCurves[2].pointTags
1597 # + divisionCurves[3].pointTags
1598 # )
1599 # gmsh.model.mesh.setTransfiniteVolume(
1600 # volumeTag, cornerTags=cornerPoints
1601 # )
1602 else:
1603 # Make all the volumes transfinite:
1604 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1606 @staticmethod
1607 def get_boundaries(volumeDimTags, returnTags=False):
1608 """
1609 Returns all the surface and line dimTags or tags of a given list of volume
1610 dimTags.
1612 :param volumeDimTags: dimTags of the volumes
1613 :type volumeDimTags: list[tuple[int, int]]
1614 :param returnTags: if True, returns tags instead of dimTags
1615 :type returnTags: bool, optional
1616 :return: surface and line dimTags or tags
1617 :rtype: tuple[list[tuple[int, int]], list[tuple[int, int]]] or
1618 tuple[list[int], list[int]]
1619 """
1620 # Get the surface tags:
1621 surfaceDimTags = list(
1622 set(
1623 gmsh.model.getBoundary(
1624 volumeDimTags,
1625 combined=False,
1626 oriented=False,
1627 recursive=False,
1628 )
1629 )
1630 )
1632 # Get the line tags:
1633 lineDimTags = list(
1634 set(
1635 gmsh.model.getBoundary(
1636 surfaceDimTags,
1637 combined=False,
1638 oriented=False,
1639 recursive=False,
1640 )
1641 )
1642 )
1644 if returnTags:
1645 surfaceTags = [dimTag[1] for dimTag in surfaceDimTags]
1646 lineTags = [dimTag[1] for dimTag in lineDimTags]
1647 return surfaceTags, lineTags
1648 else:
1649 return surfaceDimTags, lineDimTags
1651 @staticmethod
1652 def fuse_volumes(volumeDimTags, fuseSurfaces=True, fusedSurfacesArePlane=False):
1653 """
1654 Fuses all the volumes in a given list of volume dimTags, removes old volumes,
1655 and returns the new volume dimTag. Also, if compundSurfacces is True, it fuses
1656 the surfaces that only belong to the volume. fusedSurfacesArePlane can be
1657 used to change the behavior of the fuse_surfaces method.
1659 :param volumeDimTags: dimTags of the volumes
1660 :type volumeDimTags: list[tuple[int, int]]
1661 :param fuseSurfaces: if True, fuses the surfaces that only belong to the
1662 volume
1663 :type fuseSurfaces: bool, optional
1664 :param fusedSurfacesArePlane: if True, fused surfaces are assumed to be
1665 plane, and fusion is performed accordingly
1666 :return: new volume's dimTag
1667 :rtype: tuple[int, int]
1668 """
1670 # Get the combined boundary surfaces:
1671 boundarySurfacesDimTags = gmsh.model.getBoundary(
1672 volumeDimTags,
1673 combined=True,
1674 oriented=False,
1675 recursive=False,
1676 )
1677 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1679 # Get all the boundary surfaces:
1680 allBoundarySurfacesDimTags = gmsh.model.getBoundary(
1681 volumeDimTags,
1682 combined=False,
1683 oriented=False,
1684 recursive=False,
1685 )
1687 # Find internal (common) surfaces:
1688 internalSurfacesDimTags = list(
1689 set(allBoundarySurfacesDimTags) - set(boundarySurfacesDimTags)
1690 )
1692 # Get the combined boundary lines:
1693 boundaryLinesDimTags = gmsh.model.getBoundary(
1694 allBoundarySurfacesDimTags,
1695 combined=True,
1696 oriented=False,
1697 recursive=False,
1698 )
1699 boundarLinesTags = [dimTag[1] for dimTag in boundaryLinesDimTags]
1701 # Get all the boundary lines:
1702 allBoundaryLinesDimTags = gmsh.model.getBoundary(
1703 allBoundarySurfacesDimTags,
1704 combined=False,
1705 oriented=False,
1706 recursive=False,
1707 )
1709 # Find internal (common) lines:
1710 internalLinesDimTags = list(
1711 set(allBoundaryLinesDimTags) - set(boundarLinesTags)
1712 )
1714 # Remove the old volumes:
1715 removedVolumeDimTags = volumeDimTags
1716 gmsh.model.occ.remove(removedVolumeDimTags, recursive=False)
1718 # Remove the internal surfaces:
1719 gmsh.model.occ.remove(internalSurfacesDimTags, recursive=False)
1721 # Remove the internal lines:
1722 gmsh.model.occ.remove(internalLinesDimTags, recursive=False)
1724 # Create a new single volume (even thought we don't use the new volume tag
1725 # directly, it is required for finding the surfaces that only belong to the
1726 # volume):
1727 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1728 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1729 newVolumeDimTag = (3, newVolumeTag)
1730 gmsh.model.occ.synchronize()
1732 if fuseSurfaces:
1733 newVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1734 (3, newVolumeTag), surfacesArePlane=fusedSurfacesArePlane
1735 )
1737 return newVolumeDimTag
1739 @staticmethod
1740 def fuse_common_surfaces_of_two_volumes(
1741 volume1DimTags, volume2DimTags, fuseOtherSurfaces=False, surfacesArePlane=False
1742 ):
1743 """
1744 Fuses common surfaces of two volumes. Volumes are given as a list of dimTags,
1745 but they are assumed to form a single volume, and this function fuses those
1746 multiple volumes into a single volume as well. If fuseOtherSurfaces is set to
1747 True, it tries to fuse surfaces that only belong to one volume too; however,
1748 that feature is not used in Pancake3D currently.
1750 :param volume1DimTags: dimTags of the first volume
1751 :type volume1DimTags: list[tuple[int, int]]
1752 :param volume2DimTags: dimTags of the second volume
1753 :type volume2DimTags: list[tuple[int, int]]
1754 :param fuseOtherSurfaces: if True, fuses the surfaces that only belong to one
1755 volume
1756 :type fuseOtherSurfaces: bool, optional
1757 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1758 fusion is performed accordingly
1759 :type surfacesArePlane: bool, optional
1760 :return: new volumes dimTags
1761 :rtype: tuple[tuple[int, int], tuple[int, int]]
1762 """
1763 vol1BoundarySurfacesDimTags = gmsh.model.getBoundary(
1764 volume1DimTags,
1765 combined=True,
1766 oriented=False,
1767 recursive=False,
1768 )
1770 vol2BoundarySurfacesDimTags = gmsh.model.getBoundary(
1771 volume2DimTags,
1772 combined=True,
1773 oriented=False,
1774 recursive=False,
1775 )
1777 # Remove the old volumes:
1778 gmsh.model.occ.remove(volume1DimTags + volume2DimTags, recursive=False)
1780 # Find common surfaces:
1781 commonSurfacesDimTags = list(
1782 set(vol2BoundarySurfacesDimTags).intersection(
1783 set(vol1BoundarySurfacesDimTags)
1784 )
1785 )
1787 # Fuse common surfaces:
1788 fusedCommonSurfaceDimTags = Mesh.fuse_surfaces(
1789 commonSurfacesDimTags, surfacesArePlane=surfacesArePlane
1790 )
1792 # Create the new volumes:
1793 for commonSurfaceDimTag in commonSurfacesDimTags:
1794 vol1BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1795 vol2BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1797 vol1BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1798 vol1BoundarySurfaceTags = [dimTag[1] for dimTag in vol1BoundarySurfacesDimTags]
1799 vol2BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1800 vol2BoundarySurfaceTags = [dimTag[1] for dimTag in vol2BoundarySurfacesDimTags]
1802 vol1SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1803 vol1BoundarySurfaceTags, sewing=True
1804 )
1805 vol1NewVolumeDimTag = (3, gmsh.model.occ.addVolume([vol1SurfaceLoop]))
1807 vol2SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1808 vol2BoundarySurfaceTags, sewing=True
1809 )
1810 vol2NewVolumeDimTag = (
1811 3,
1812 gmsh.model.occ.addVolume([vol2SurfaceLoop]),
1813 )
1815 gmsh.model.occ.synchronize()
1817 if fuseOtherSurfaces:
1818 vol1NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1819 vol1NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1820 )
1821 vol2NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1822 vol2NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1823 )
1825 return vol1NewVolumeDimTag, vol2NewVolumeDimTag
1827 @staticmethod
1828 def fuse_possible_surfaces_of_a_volume(volumeDimTag, surfacesArePlane=False):
1829 """
1830 Fuses surfaces that only belong to the volumeDimTag.
1832 :param volumeDimTag: dimTag of the volume
1833 :type volumeDimTag: tuple[int, int]
1834 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1835 fusion is performed accordingly
1836 :type surfacesArePlane: bool, optional
1837 :return: new volume dimTag
1838 :rtype: tuple[int, int]
1839 """
1840 boundarySurfacesDimTags = gmsh.model.getBoundary(
1841 [volumeDimTag],
1842 combined=True,
1843 oriented=False,
1844 recursive=False,
1845 )
1846 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1848 # Combine surfaces that only belong to the volume:
1849 toBeFusedSurfacesDimTags = []
1850 surfacesNormals = []
1851 for surfaceDimTag in boundarySurfacesDimTags:
1852 upward, _ = gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1854 if len(list(upward)) == 1:
1855 toBeFusedSurfacesDimTags.append(surfaceDimTag)
1856 # Get the normal of the surface:
1857 surfacesNormals.append(
1858 list(gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5]))
1859 )
1861 # Remove the old volume (it is not required anymore):
1862 gmsh.model.occ.remove([volumeDimTag], recursive=False)
1863 gmsh.model.occ.synchronize()
1865 # Categorize surfaces based on their normals so that they can be combined
1866 # correctly. Without these, perpendicular surfaces will cause problems.
1868 # Define a threshold to determine if two surface normals are similar or not
1869 threshold = 1e-6
1871 # Initialize an empty list to store the sets of surfaces
1872 setsOfSurfaces = []
1874 # Calculate the Euclidean distance between each pair of objects
1875 for i in range(len(toBeFusedSurfacesDimTags)):
1876 surfaceDimTag = toBeFusedSurfacesDimTags[i]
1877 surfaceTouchingVolumeTags, _ = list(
1878 gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1879 )
1880 surfaceNormal = surfacesNormals[i]
1881 assignedToASet = False
1883 for surfaceSet in setsOfSurfaces:
1884 representativeSurfaceDimTag = surfaceSet[0]
1885 representativeSurfaceTouchingVolumeTags, _ = list(
1886 gmsh.model.getAdjacencies(
1887 representativeSurfaceDimTag[0],
1888 representativeSurfaceDimTag[1],
1889 )
1890 )
1891 representativeNormal = list(
1892 gmsh.model.getNormal(representativeSurfaceDimTag[1], [0.5, 0.5])
1893 )
1895 # Calculate the difference between surfaceNormal and
1896 # representativeNormal:
1897 difference = math.sqrt(
1898 sum(
1899 (x - y) ** 2
1900 for x, y in zip(surfaceNormal, representativeNormal)
1901 )
1902 )
1904 # Check if the distance is below the threshold
1905 if difference < threshold and set(surfaceTouchingVolumeTags) == set(
1906 representativeSurfaceTouchingVolumeTags
1907 ):
1908 # Add the object to an existing category
1909 surfaceSet.append(surfaceDimTag)
1910 assignedToASet = True
1911 break
1913 if not assignedToASet:
1914 # Create a new category with the current object if none of the
1915 # existing sets match
1916 setsOfSurfaces.append([surfaceDimTag])
1918 for surfaceSet in setsOfSurfaces:
1919 if len(surfaceSet) > 1:
1920 oldSurfaceDimTags = surfaceSet
1921 newSurfaceDimTags = Mesh.fuse_surfaces(
1922 oldSurfaceDimTags, surfacesArePlane=surfacesArePlane
1923 )
1924 newSurfaceTags = [dimTag[1] for dimTag in newSurfaceDimTags]
1926 oldSurfaceTags = [dimTag[1] for dimTag in oldSurfaceDimTags]
1927 boundarSurfacesTags = [
1928 tag for tag in boundarSurfacesTags if tag not in oldSurfaceTags
1929 ]
1930 boundarSurfacesTags.extend(newSurfaceTags)
1932 # Create a new single volume:
1933 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1934 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1935 gmsh.model.occ.synchronize()
1937 return (3, newVolumeTag)
1939 @staticmethod
1940 def fuse_surfaces(surfaceDimTags, surfacesArePlane=False, categorizeSurfaces=False):
1941 """
1942 Fuses all the surfaces in a given list of surface dimTags, removes the old
1943 surfaces, and returns the new dimTags. If surfacesArePlane is True, the surfaces
1944 are assumed to be plane, and fusing will be done without gmsh.model.occ.fuse
1945 method, which is faster.
1947 :param surfaceDimTags: dimTags of the surfaces
1948 :type surfaceDimTags: list[tuple[int, int]]
1949 :param surfacesArePlane: if True, surfaces are assumed to be plane
1950 :type surfacesArePlane: bool, optional
1951 :return: newly created surface dimTags
1952 :rtype: list[tuple[int, int]]
1953 """
1954 oldSurfaceDimTags = surfaceDimTags
1956 if surfacesArePlane:
1957 # Get the combined boundary curves:
1958 boundaryCurvesDimTags = gmsh.model.getBoundary(
1959 oldSurfaceDimTags,
1960 combined=True,
1961 oriented=False,
1962 recursive=False,
1963 )
1965 # Get all the boundary curves:
1966 allCurvesDimTags = gmsh.model.getBoundary(
1967 oldSurfaceDimTags,
1968 combined=False,
1969 oriented=False,
1970 recursive=False,
1971 )
1973 # Find internal (common) curves:
1974 internalCurvesDimTags = list(
1975 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
1976 )
1978 # Remove the old surfaces:
1979 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
1981 # Remove the internal curves:
1982 gmsh.model.occ.remove(internalCurvesDimTags, recursive=True)
1984 # Create a new single surface:
1985 def findOuterOnes(dimTags, findInnerOnes=False):
1986 """
1987 Finds the outermost surface/curve/point in a list of dimTags. The outermost means
1988 the furthest from the origin.
1989 """
1990 dim = dimTags[0][0]
1992 if dim == 2:
1993 distances = []
1994 for dimTag in dimTags:
1995 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])
1996 for curve in curves:
1997 curve = list(curve)
1998 gmsh.model.occ.synchronize()
1999 pointTags = gmsh.model.getBoundary(
2000 [(1, curveTag) for curveTag in curve],
2001 oriented=False,
2002 combined=False,
2003 )
2004 # Get the positions of the points:
2005 points = []
2006 for dimTag in pointTags:
2007 boundingbox1 = gmsh.model.occ.getBoundingBox(
2008 0, dimTag[1]
2009 )[:3]
2010 boundingbox2 = gmsh.model.occ.getBoundingBox(
2011 0, dimTag[1]
2012 )[3:]
2013 boundingbox = list(
2014 map(operator.add, boundingbox1, boundingbox2)
2015 )
2016 points.append(
2017 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2018 )
2020 distances.append(
2021 max([point[0] ** 2 + point[1] ** 2 for point in points])
2022 )
2023 elif dim == 1:
2024 distances = []
2025 for dimTag in dimTags:
2026 gmsh.model.occ.synchronize()
2027 pointTags = gmsh.model.getBoundary(
2028 [dimTag],
2029 oriented=False,
2030 combined=False,
2031 )
2032 # Get the positions of the points:
2033 points = []
2034 for dimTag in pointTags:
2035 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2036 :3
2037 ]
2038 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2039 3:
2040 ]
2041 boundingbox = list(
2042 map(operator.add, boundingbox1, boundingbox2)
2043 )
2044 points.append(
2045 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2046 )
2048 distances.append(
2049 max([point[0] ** 2 + point[1] ** 2 for point in points])
2050 )
2052 if findInnerOnes:
2053 goalDistance = min(distances)
2054 else:
2055 goalDistance = max(distances)
2057 result = []
2058 for distance, dimTag in zip(distances, dimTags):
2059 # Return all the dimTags with the hoal distance:
2060 if math.isclose(distance, goalDistance, abs_tol=1e-6):
2061 result.append(dimTag)
2063 return result
2065 outerCurvesDimTags = findOuterOnes(boundaryCurvesDimTags)
2066 outerCurvesTags = [dimTag[1] for dimTag in outerCurvesDimTags]
2067 curveLoopOuter = gmsh.model.occ.addCurveLoop(outerCurvesTags)
2069 innerCurvesDimTags = findOuterOnes(
2070 boundaryCurvesDimTags, findInnerOnes=True
2071 )
2072 innerCurvesTags = [dimTag[1] for dimTag in innerCurvesDimTags]
2073 curveLoopInner = gmsh.model.occ.addCurveLoop(innerCurvesTags)
2075 newSurfaceTag = gmsh.model.occ.addPlaneSurface(
2076 [curveLoopOuter, curveLoopInner]
2077 )
2079 gmsh.model.occ.synchronize()
2081 return [(2, newSurfaceTag)]
2082 else:
2083 # Create a new single surface:
2084 # The order of tags seems to be important for the fuse method to work
2085 # and not crash with a segmentation fault.
2086 try:
2087 fuseResults = gmsh.model.occ.fuse(
2088 [oldSurfaceDimTags[0]],
2089 oldSurfaceDimTags[1:],
2090 removeObject=False,
2091 removeTool=False,
2092 )
2093 newSurfaceDimTags = fuseResults[0]
2094 except:
2095 return oldSurfaceDimTags
2097 # Get the combined boundary curves:
2098 gmsh.model.occ.synchronize()
2099 boundaryCurvesDimTags = gmsh.model.getBoundary(
2100 newSurfaceDimTags,
2101 combined=True,
2102 oriented=False,
2103 recursive=False,
2104 )
2106 # Get all the boundary curves:
2107 allCurvesDimTags = gmsh.model.getBoundary(
2108 oldSurfaceDimTags,
2109 combined=False,
2110 oriented=False,
2111 recursive=False,
2112 )
2114 # Find internal (common) curves:
2115 internalCurvesDimTags = list(
2116 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
2117 )
2119 # Remove the old surfaces:
2120 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
2122 # Remove the internal curves:
2123 gmsh.model.occ.remove(internalCurvesDimTags, recursive=False)
2125 gmsh.model.occ.synchronize()
2127 return newSurfaceDimTags
2129 def generate_regions(self):
2130 """
2131 Generates physical groups and the regions file. Physical groups are generated in
2132 GMSH, and their tags and names are saved in the regions file. FiQuS use the
2133 regions file to create the corresponding .pro file.
2135 .vi file sends the information about geometry from geometry class to mesh class.
2136 .regions file sends the information about the physical groups formed out of
2137 elementary entities from the mesh class to the solution class.
2139 The file extension for the regions file is custom because users should not edit
2140 or even see this file.
2142 Regions are generated in the meshing part because BREP files cannot store
2143 regions.
2144 """
2145 logger.info("Generating physical groups and regions file has been started.")
2146 start_time = timeit.default_timer()
2148 # Create regions instance to both generate regions file and physical groups:
2149 self.regions = regions()
2151 # ==============================================================================
2152 # WINDING AND CONTACT LAYER REGIONS START =========================================
2153 # ==============================================================================
2154 if not self.geo.contactLayer.thinShellApproximation:
2155 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.winding.name]]
2156 self.regions.addEntities(
2157 self.geo.winding.name, windingTags, regionType.powered, entityType.vol
2158 )
2160 insulatorTags = [dimTag[1] for dimTag in self.dimTags[self.geo.contactLayer.name]]
2162 terminalDimTags = (
2163 self.dimTags[self.geo.terminals.inner.name] + self.dimTags[self.geo.terminals.outer.name]
2164 )
2165 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2166 terminalDimTags, combined=False, oriented=False
2167 )
2168 transitionNotchSurfaces = gmsh.model.getBoundary(
2169 self.dimTags["innerTransitionNotch"]
2170 + self.dimTags["outerTransitionNotch"],
2171 combined=False,
2172 oriented=False,
2173 )
2175 contactLayer = []
2176 contactLayerBetweenTerminalsAndWinding = []
2177 for insulatorTag in insulatorTags:
2178 insulatorSurfaces = gmsh.model.getBoundary(
2179 [(3, insulatorTag)], combined=False, oriented=False
2180 )
2181 itTouchesTerminals = False
2182 for insulatorSurface in insulatorSurfaces:
2183 if (
2184 insulatorSurface
2185 in terminalAndNotchSurfaces + transitionNotchSurfaces
2186 ):
2187 contactLayerBetweenTerminalsAndWinding.append(insulatorTag)
2188 itTouchesTerminals = True
2189 break
2191 if not itTouchesTerminals:
2192 contactLayer.append(insulatorTag)
2194 self.regions.addEntities(
2195 self.geo.contactLayer.name, contactLayer, regionType.insulator, entityType.vol
2196 )
2198 self.regions.addEntities(
2199 "WindingAndTerminalContactLayer",
2200 contactLayerBetweenTerminalsAndWinding,
2201 regionType.insulator,
2202 entityType.vol,
2203 )
2204 else:
2205 # Calculate the number of stacks for each individual winding. Number of
2206 # stacks is the number of volumes per turn. It affects how the regions
2207 # are created because of the TSA's pro file formulation.
2209 # find the smallest prime number that divides NofVolumes:
2210 windingDimTags = self.dimTags[self.geo.winding.name + "1"]
2211 windingTags = [dimTag[1] for dimTag in windingDimTags]
2212 NofVolumes = self.geo.winding.numberOfVolumesPerTurn
2213 smallest_prime_divisor = 2
2214 while NofVolumes % smallest_prime_divisor != 0:
2215 smallest_prime_divisor += 1
2217 # the number of stacks is the region divison per turn:
2218 NofStacks = smallest_prime_divisor
2220 # the number of sets are the total number of regions for all windings and
2221 # contact layers:
2222 NofSets = 2 * NofStacks
2224 allInnerTerminalSurfaces = gmsh.model.getBoundary(
2225 self.dimTags[self.geo.terminals.inner.name] + self.dimTags["innerTransitionNotch"],
2226 combined=False,
2227 oriented=False,
2228 )
2229 allInnerTerminalContactLayerSurfaces = []
2230 for innerTerminalSurface in allInnerTerminalSurfaces:
2231 normal = gmsh.model.getNormal(innerTerminalSurface[1], [0.5, 0.5])
2232 if abs(normal[2]) < 1e-5:
2233 curves = gmsh.model.getBoundary(
2234 [innerTerminalSurface], combined=False, oriented=False
2235 )
2236 curveTags = [dimTag[1] for dimTag in curves]
2237 for curveTag in curveTags:
2238 curveObject = curve(curveTag, self.geo)
2239 if curveObject.type is curveType.spiralArc:
2240 allInnerTerminalContactLayerSurfaces.append(
2241 innerTerminalSurface[1]
2242 )
2244 finalWindingSets = []
2245 finalContactLayerSets = []
2246 for i in range(NofSets):
2247 finalWindingSets.append([])
2248 finalContactLayerSets.append([])
2250 for i in range(self.geo.numberOfPancakes):
2251 windingDimTags = self.dimTags[self.geo.winding.name + str(i + 1)]
2252 windingTags = [dimTag[1] for dimTag in windingDimTags]
2254 NofVolumes = len(windingDimTags)
2256 windings = []
2257 for windingTag in windingTags:
2258 surfaces = gmsh.model.getBoundary(
2259 [(3, windingTag)], combined=False, oriented=False
2260 )
2261 curves = gmsh.model.getBoundary(
2262 surfaces, combined=False, oriented=False
2263 )
2264 curveTags = list(set([dimTag[1] for dimTag in curves]))
2265 for curveTag in curveTags:
2266 curveObject = curve(curveTag, self.geo)
2267 if curveObject.type is curveType.spiralArc:
2268 windingVolumeLengthInTurns = abs(
2269 curveObject.n2 - curveObject.n1
2270 )
2271 if windingVolumeLengthInTurns > 0.5:
2272 # The arc can never be longer than half a turn.
2273 windingVolumeLengthInTurns = (
2274 1 - windingVolumeLengthInTurns
2275 )
2277 windings.append((windingTag, windingVolumeLengthInTurns))
2279 windingStacks = []
2280 while len(windings) > 0:
2281 stack = []
2282 stackLength = 0
2283 for windingTag, windingVolumeLengthInTurns in windings:
2284 if stackLength < 1 / NofStacks - 1e-6:
2285 stack.append(windingTag)
2286 stackLength += windingVolumeLengthInTurns
2287 else:
2288 break
2289 # remove all the windings that are already added to the stack:
2290 windings = [
2291 (windingTag, windingVolumeLengthInTurns)
2292 for windingTag, windingVolumeLengthInTurns in windings
2293 if windingTag not in stack
2294 ]
2296 # find spiral surfaces of the stack:
2297 stackDimTags = [(3, windingTag) for windingTag in stack]
2298 stackSurfacesDimTags = gmsh.model.getBoundary(
2299 stackDimTags, combined=True, oriented=False
2300 )
2301 stackCurvesDimTags = gmsh.model.getBoundary(
2302 stackSurfacesDimTags, combined=False, oriented=False
2303 )
2304 # find the curve furthest from the origin:
2305 curveObjects = []
2306 for curveDimTag in stackCurvesDimTags:
2307 curveObject = curve(curveDimTag[1], self.geo)
2308 if curveObject.type is curveType.spiralArc:
2309 curveObjectDistanceFromOrigin = math.sqrt(
2310 curveObject.points[0][0] ** 2
2311 + curveObject.points[0][1] ** 2
2312 )
2313 curveObjects.append(
2314 (curveObject, curveObjectDistanceFromOrigin)
2315 )
2317 # sort the curves based on their distance from the origin (furthest first)
2318 curveObjects.sort(key=lambda x: x[1], reverse=True)
2320 curveTags = [curveObject[0].tag for curveObject in curveObjects]
2322 # only keep half of the curveTags:
2323 furthestCurveTags = curveTags[: len(curveTags) // 2]
2325 stackSpiralSurfaces = []
2326 for surfaceDimTag in stackSurfacesDimTags:
2327 normal = gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5])
2328 if abs(normal[2]) < 1e-5:
2329 curves = gmsh.model.getBoundary(
2330 [surfaceDimTag], combined=False, oriented=False
2331 )
2332 curveTags = [dimTag[1] for dimTag in curves]
2333 for curveTag in curveTags:
2334 if curveTag in furthestCurveTags:
2335 stackSpiralSurfaces.append(surfaceDimTag[1])
2336 break
2338 # add inner terminal surfaces too:
2339 if len(windingStacks) >= NofStacks:
2340 correspondingWindingStack = windingStacks[
2341 len(windingStacks) - NofStacks
2342 ]
2343 correspondingWindings = correspondingWindingStack[0]
2344 correspondingSurfaces = gmsh.model.getBoundary(
2345 [(3, windingTag) for windingTag in correspondingWindings],
2346 combined=True,
2347 oriented=False,
2348 )
2349 correspondingSurfaceTags = [
2350 dimTag[1] for dimTag in correspondingSurfaces
2351 ]
2352 for surface in allInnerTerminalContactLayerSurfaces:
2353 if surface in correspondingSurfaceTags:
2354 stackSpiralSurfaces.append(surface)
2356 windingStacks.append((stack, stackSpiralSurfaces))
2358 windingSets = []
2359 contactLayerSets = []
2360 for j in range(NofSets):
2361 windingTags = [
2362 windingTags for windingTags, _ in windingStacks[j::NofSets]
2363 ]
2364 windingTags = list(itertools.chain.from_iterable(windingTags))
2366 surfaceTags = [
2367 surfaceTags for _, surfaceTags in windingStacks[j::NofSets]
2368 ]
2369 surfaceTags = list(itertools.chain.from_iterable(surfaceTags))
2371 windingSets.append(windingTags)
2372 contactLayerSets.append(surfaceTags)
2374 # windingSets is a list with a length of NofSets.
2375 # finalWindingSets is also a list with a length of NofSets.
2376 for j, (windingSet, contactLayerSet) in enumerate(
2377 zip(windingSets, contactLayerSets)
2378 ):
2379 finalWindingSets[j].extend(windingSet)
2380 finalContactLayerSets[j].extend(contactLayerSet)
2382 # Seperate transition layer:
2383 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2384 self.dimTags[self.geo.terminals.inner.name]
2385 + self.dimTags[self.geo.terminals.outer.name]
2386 + self.dimTags["innerTransitionNotch"]
2387 + self.dimTags["outerTransitionNotch"],
2388 combined=False,
2389 oriented=False,
2390 )
2391 terminalAndNotchSurfaceTags = set(
2392 [dimTag[1] for dimTag in terminalAndNotchSurfaces]
2393 )
2395 contactLayerSets = []
2396 terminalWindingContactLayerSets = []
2397 for j in range(NofSets):
2398 contactLayerSets.append([])
2399 terminalWindingContactLayerSets.append([])
2401 for j in range(NofSets):
2402 allContactLayersInTheSet = finalContactLayerSets[j]
2404 insulatorList = []
2405 windingTerminalInsulatorList = []
2406 for contactLayer in allContactLayersInTheSet:
2407 if contactLayer in terminalAndNotchSurfaceTags:
2408 windingTerminalInsulatorList.append(contactLayer)
2409 else:
2410 insulatorList.append(contactLayer)
2412 contactLayerSets[j].extend(set(insulatorList))
2413 terminalWindingContactLayerSets[j].extend(set(windingTerminalInsulatorList))
2415 allContactLayerSurfacesForAllPancakes = []
2416 for j in range(NofSets):
2417 # Add winding volumes:
2418 self.regions.addEntities(
2419 self.geo.winding.name + "-" + str(j + 1),
2420 finalWindingSets[j],
2421 regionType.powered,
2422 entityType.vol,
2423 )
2425 # Add insulator surfaces:
2426 self.regions.addEntities(
2427 self.geo.contactLayer.name + "-" + str(j + 1),
2428 contactLayerSets[j],
2429 regionType.insulator,
2430 entityType.surf,
2431 )
2432 allContactLayerSurfacesForAllPancakes.extend(contactLayerSets[j])
2434 # Add terminal and winding contact layer:
2435 allContactLayerSurfacesForAllPancakes.extend(
2436 terminalWindingContactLayerSets[j]
2437 )
2439 # Add intersection of transition notch boundary and WindingAndTerminalContactLayer:
2440 transitionNotchSurfaces = gmsh.model.getBoundary(
2441 self.dimTags["innerTransitionNotch"]
2442 + self.dimTags["outerTransitionNotch"],
2443 combined=False,
2444 oriented=False,
2445 recursive=False
2446 )
2448 terminalContactLayerMinusNotch = set(terminalWindingContactLayerSets[j]).difference([tag for (dim, tag) in transitionNotchSurfaces])
2450 self.regions.addEntities(
2451 "WindingAndTerminalContactLayerWithoutNotch" + "-" + str(j + 1),
2452 list(terminalContactLayerMinusNotch),
2453 regionType.insulator,
2454 entityType.surf,
2455 )
2457 notchAndTerminalContactLayerIntersection = set([tag for (dim, tag) in transitionNotchSurfaces]).intersection(terminalWindingContactLayerSets[j])
2459 self.regions.addEntities(
2460 "WindingAndTerminalContactLayerOnlyNotch" + "-" + str(j + 1),
2461 list(notchAndTerminalContactLayerIntersection),
2462 regionType.insulator,
2463 entityType.surf,
2464 )
2466 allContactLayerSurfacesForAllPancakes = list(
2467 set(allContactLayerSurfacesForAllPancakes)
2468 )
2469 # Get insulator's boundary line that touches the air (required for the
2470 # pro file formulation):
2471 allContactLayerSurfacesForAllPancakesDimTags = [
2472 (2, surfaceTag) for surfaceTag in allContactLayerSurfacesForAllPancakes
2473 ]
2474 insulatorBoundary = gmsh.model.getBoundary(
2475 allContactLayerSurfacesForAllPancakesDimTags,
2476 combined=True,
2477 oriented=False,
2478 )
2479 insulatorBoundaryTags = [dimTag[1] for dimTag in insulatorBoundary]
2481 # Add insulator boundary lines:
2482 # Vertical lines should be removed from the insulator boundary because
2483 # they touch the terminals, not the air:
2484 verticalInsulatorBoundaryTags = []
2485 insulatorBoundaryTagsCopy = insulatorBoundaryTags.copy()
2486 for lineTag in insulatorBoundaryTagsCopy:
2487 lineObject = curve(lineTag, self.geo)
2488 if lineObject.type is curveType.axial:
2489 verticalInsulatorBoundaryTags.append(lineTag)
2490 insulatorBoundaryTags.remove(lineTag)
2492 # Create regions:
2493 self.regions.addEntities(
2494 self.geo.contactLayerBoundaryName,
2495 insulatorBoundaryTags,
2496 regionType.insulator,
2497 entityType.curve,
2498 )
2499 self.regions.addEntities(
2500 self.geo.contactLayerBoundaryName + "-TouchingTerminal",
2501 verticalInsulatorBoundaryTags,
2502 regionType.insulator,
2503 entityType.curve,
2504 )
2506 innerTransitionNotchTags = [
2507 dimTag[1] for dimTag in self.dimTags["innerTransitionNotch"]
2508 ]
2509 outerTransitionNotchTags = [
2510 dimTag[1] for dimTag in self.dimTags["outerTransitionNotch"]
2511 ]
2512 self.regions.addEntities(
2513 "innerTransitionNotch",
2514 innerTransitionNotchTags,
2515 regionType.powered,
2516 entityType.vol,
2517 )
2518 self.regions.addEntities(
2519 "outerTransitionNotch",
2520 outerTransitionNotchTags,
2521 regionType.powered,
2522 entityType.vol,
2523 )
2524 # ==============================================================================
2525 # WINDING AND CONTACT LAYER REGIONS ENDS =======================================
2526 # ==============================================================================
2528 # ==============================================================================
2529 # TERMINAL REGIONS START =======================================================
2530 # ==============================================================================
2532 innerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.terminals.inner.name]]
2533 self.regions.addEntities(
2534 self.geo.terminals.inner.name, innerTerminalTags, regionType.powered, entityType.vol_in
2535 )
2536 outerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.terminals.outer.name]]
2537 self.regions.addEntities(
2538 self.geo.terminals.outer.name,
2539 outerTerminalTags,
2540 regionType.powered,
2541 entityType.vol_out,
2542 )
2544 # Top and bottom terminal surfaces:
2545 firstTerminalDimTags = self.dimTags[self.geo.terminals.firstName]
2546 lastTerminalDimTags = self.dimTags[self.geo.terminals.lastName]
2548 if self.mesh.terminals.structured:
2549 topSurfaceDimTags = []
2550 for i in [1, 2, 3, 4]:
2551 lastTerminalSurfaces = gmsh.model.getBoundary(
2552 [lastTerminalDimTags[-i]], combined=False, oriented=False
2553 )
2554 topSurfaceDimTags.append(lastTerminalSurfaces[-1])
2555 else:
2556 lastTerminalSurfaces = gmsh.model.getBoundary(
2557 [lastTerminalDimTags[-1]], combined=False, oriented=False
2558 )
2559 topSurfaceDimTags = [lastTerminalSurfaces[-1]]
2560 topSurfaceTags = [dimTag[1] for dimTag in topSurfaceDimTags]
2562 if self.mesh.terminals.structured:
2563 bottomSurfaceDimTags = []
2564 for i in [1, 2, 3, 4]:
2565 firstTerminalSurfaces = gmsh.model.getBoundary(
2566 [firstTerminalDimTags[-i]], combined=False, oriented=False
2567 )
2568 bottomSurfaceDimTags.append(firstTerminalSurfaces[-1])
2569 else:
2570 firstTerminalSurfaces = gmsh.model.getBoundary(
2571 [firstTerminalDimTags[-1]], combined=False, oriented=False
2572 )
2573 bottomSurfaceDimTags = [firstTerminalSurfaces[-1]]
2574 bottomSurfaceTags = [dimTag[1] for dimTag in bottomSurfaceDimTags]
2576 self.regions.addEntities(
2577 "TopSurface",
2578 topSurfaceTags,
2579 regionType.powered,
2580 entityType.surf_out,
2581 )
2582 self.regions.addEntities(
2583 "BottomSurface",
2584 bottomSurfaceTags,
2585 regionType.powered,
2586 entityType.surf_in,
2587 )
2589 # if self.geo.contactLayer.tsa:
2590 # outerTerminalSurfaces = gmsh.model.getBoundary(
2591 # self.dimTags[self.geo.terminals.o.name], combined=True, oriented=False
2592 # )
2593 # outerTerminalSurfaces = [dimTag[1] for dimTag in outerTerminalSurfaces]
2594 # innerTerminalSurfaces = gmsh.model.getBoundary(
2595 # self.dimTags[self.geo.terminals.i.name], combined=True, oriented=False
2596 # )
2597 # innerTerminalSurfaces = [dimTag[1] for dimTag in innerTerminalSurfaces]
2598 # windingSurfaces = gmsh.model.getBoundary(
2599 # self.dimTags[self.geo.winding.name] + self.dimTags[self.geo.contactLayer.name],
2600 # combined=True,
2601 # oriented=False,
2602 # )
2603 # windingSurfaces = [dimTag[1] for dimTag in windingSurfaces]
2605 # windingAndOuterTerminalCommonSurfaces = list(
2606 # set(windingSurfaces).intersection(set(outerTerminalSurfaces))
2607 # )
2608 # windingAndInnerTerminalCommonSurfaces = list(
2609 # set(windingSurfaces).intersection(set(innerTerminalSurfaces))
2610 # )
2612 # self.regions.addEntities(
2613 # "WindingAndTerminalContactLayer",
2614 # windingAndOuterTerminalCommonSurfaces
2615 # + windingAndInnerTerminalCommonSurfaces,
2616 # regionType.insulator,
2617 # entityType.surf,
2618 # )
2620 # ==============================================================================
2621 # TERMINAL REGIONS ENDS ========================================================
2622 # ==============================================================================
2624 # ==============================================================================
2625 # AIR AND AIR SHELL REGIONS STARTS =============================================
2626 # ==============================================================================
2627 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name]]
2628 self.regions.addEntities(
2629 self.geo.air.name, airTags, regionType.air, entityType.vol
2630 )
2632 # Create a region with two points on air to be used in the pro file formulation:
2633 # To those points, Phi=0 boundary condition will be applied to set the gauge.
2634 outerAirSurfaces = gmsh.model.getBoundary(
2635 self.dimTags[self.geo.air.name + "-OuterTube"], combined=True, oriented=False
2636 )
2637 outerAirSurface = outerAirSurfaces[-1]
2638 outerAirCurves = gmsh.model.getBoundary(
2639 [outerAirSurface], combined=True, oriented=False
2640 )
2641 outerAirCurve = outerAirCurves[-1]
2642 outerAirPoint = gmsh.model.getBoundary(
2643 [outerAirCurve], combined=False, oriented=False
2644 )
2645 outerAirPointTag = outerAirPoint[0][1]
2646 self.regions.addEntities(
2647 "OuterAirPoint",
2648 [outerAirPointTag],
2649 regionType.air,
2650 entityType.point,
2651 )
2653 innerAirSurfaces = gmsh.model.getBoundary(
2654 self.dimTags[self.geo.air.name + "-InnerCylinder"],
2655 combined=True,
2656 oriented=False,
2657 )
2658 innerAirSurface = innerAirSurfaces[0]
2659 innerAirCurves = gmsh.model.getBoundary(
2660 [innerAirSurface], combined=True, oriented=False
2661 )
2662 innerAirCurve = innerAirCurves[-1]
2663 innerAirPoint = gmsh.model.getBoundary(
2664 [innerAirCurve], combined=False, oriented=False
2665 )
2666 innerAirPointTag = innerAirPoint[0][1]
2667 self.regions.addEntities(
2668 "InnerAirPoint",
2669 [innerAirPointTag],
2670 regionType.air,
2671 entityType.point,
2672 )
2674 if self.geo.air.shellTransformation:
2675 if self.geo.air.type == "cylinder":
2676 airShellTags = [
2677 dimTag[1] for dimTag in self.dimTags[self.geo.air.shellVolumeName]
2678 ]
2679 self.regions.addEntities(
2680 self.geo.air.shellVolumeName,
2681 airShellTags,
2682 regionType.air_far_field,
2683 entityType.vol,
2684 )
2685 elif self.geo.air.type == "cuboid":
2686 airShell1Tags = [
2687 dimTag[1]
2688 for dimTag in self.dimTags[self.geo.air.shellVolumeName + "-Part1"]
2689 + self.dimTags[self.geo.air.shellVolumeName + "-Part3"]
2690 ]
2691 airShell2Tags = [
2692 dimTag[1]
2693 for dimTag in self.dimTags[self.geo.air.shellVolumeName + "-Part2"]
2694 + self.dimTags[self.geo.air.shellVolumeName + "-Part4"]
2695 ]
2696 self.regions.addEntities(
2697 self.geo.air.shellVolumeName + "-PartX",
2698 airShell1Tags,
2699 regionType.air_far_field,
2700 entityType.vol,
2701 )
2702 self.regions.addEntities(
2703 self.geo.air.shellVolumeName + "-PartY",
2704 airShell2Tags,
2705 regionType.air_far_field,
2706 entityType.vol,
2707 )
2708 # ==============================================================================
2709 # AIR AND AIR SHELL REGIONS ENDS ===============================================
2710 # ==============================================================================
2712 # ==============================================================================
2713 # CUTS STARTS ==================================================================
2714 # ==============================================================================
2715 if self.geo.air.cutName in self.dimTags:
2716 cutTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.cutName]]
2717 self.regions.addEntities(
2718 self.geo.air.cutName, cutTags, regionType.air, entityType.cochain
2719 )
2721 if "CutsForPerfectInsulation" in self.dimTags:
2722 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsForPerfectInsulation"]]
2723 self.regions.addEntities(
2724 "CutsForPerfectInsulation", cutTags, regionType.air, entityType.cochain
2725 )
2727 if "CutsBetweenPancakes" in self.dimTags:
2728 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsBetweenPancakes"]]
2729 self.regions.addEntities(
2730 "CutsBetweenPancakes", cutTags, regionType.air, entityType.cochain
2731 )
2732 # ==============================================================================
2733 # CUTS ENDS ====================================================================
2734 # ==============================================================================
2736 # ==============================================================================
2737 # PANCAKE BOUNDARY SURFACE STARTS ==============================================
2738 # ==============================================================================
2739 # Pancake3D Boundary Surface:
2740 allPancakeVolumes = (
2741 self.dimTags[self.geo.winding.name]
2742 + self.dimTags[self.geo.terminals.inner.name]
2743 + self.dimTags[self.geo.terminals.outer.name]
2744 + self.dimTags[self.geo.contactLayer.name]
2745 + self.dimTags["innerTransitionNotch"]
2746 + self.dimTags["outerTransitionNotch"]
2747 )
2748 Pancake3DAllBoundary = gmsh.model.getBoundary(
2749 allPancakeVolumes, combined=True, oriented=False
2750 )
2751 Pancake3DBoundaryDimTags = list(
2752 set(Pancake3DAllBoundary)
2753 - set(topSurfaceDimTags)
2754 - set(bottomSurfaceDimTags)
2755 )
2756 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2757 self.regions.addEntities(
2758 self.geo.pancakeBoundaryName,
2759 pancake3DBoundaryTags,
2760 regionType.powered,
2761 entityType.surf,
2762 )
2764 if self.geo.contactLayer.thinShellApproximation:
2765 # add non-winding boundary for convective cooling
2766 windingBoundaryDimTags = gmsh.model.getBoundary(
2767 [(3, tag) for tag in itertools.chain(*finalWindingSets)],
2768 combined=True,
2769 oriented=False,
2770 )
2772 inner_terminal_and_transition_notch_all_boundaries = gmsh.model.getBoundary(
2773 self.dimTags[self.geo.terminals.inner.name] + self.dimTags["innerTransitionNotch"],
2774 combined=True,
2775 oriented=False
2776 )
2778 inner_terminal_and_transition_notch_boundary_dim_tags = set(Pancake3DBoundaryDimTags).intersection(inner_terminal_and_transition_notch_all_boundaries)
2779 inner_terminal_and_transition_notch_boundary_tags = [dimTag[1] for dimTag in inner_terminal_and_transition_notch_boundary_dim_tags]
2780 self.regions.addEntities(
2781 f"{self.geo.pancakeBoundaryName}-InnerTerminalAndTransitionNotch",
2782 inner_terminal_and_transition_notch_boundary_tags,
2783 regionType.powered,
2784 entityType.surf_th,
2785 )
2787 outer_terminal_and_transition_notch_all_boundaries = gmsh.model.getBoundary(
2788 self.dimTags[self.geo.terminals.outer.name] + self.dimTags["outerTransitionNotch"],
2789 combined=True,
2790 oriented=False
2791 )
2793 outer_terminal_and_transition_notch_boundary_dim_tags = set(Pancake3DBoundaryDimTags).intersection(outer_terminal_and_transition_notch_all_boundaries)
2794 outer_terminal_and_transition_notch_boundary_tags = [dimTag[1] for dimTag in outer_terminal_and_transition_notch_boundary_dim_tags]
2795 self.regions.addEntities(
2796 f"{self.geo.pancakeBoundaryName}-outerTerminalAndTransitionNotch",
2797 outer_terminal_and_transition_notch_boundary_tags,
2798 regionType.powered,
2799 entityType.surf_th,
2800 )
2802 # add pancake boundary for convective cooling following the winding numbering logic
2803 # computes the intersection between pancake boundary and the boundary of each winding group
2804 for j in range(NofSets):
2806 windingBoundaryDimTags = gmsh.model.getBoundary(
2807 [(3, tag) for tag in finalWindingSets[j]],
2808 combined=True,
2809 oriented=False,
2810 )
2812 windingBoundaryDimTags = set(windingBoundaryDimTags).intersection(Pancake3DBoundaryDimTags)
2814 windingBoundaryTags = [dimTag[1] for dimTag in windingBoundaryDimTags]
2815 self.regions.addEntities(
2816 f"{self.geo.pancakeBoundaryName}-Winding{j+1}",
2817 windingBoundaryTags,
2818 regionType.powered,
2819 entityType.surf_th,
2820 )
2822 if not self.geo.contactLayer.thinShellApproximation:
2823 # Pancake3D Boundary Surface with only winding and terminals:
2824 allPancakeVolumes = (
2825 self.dimTags[self.geo.winding.name]
2826 + self.dimTags[self.geo.terminals.inner.name]
2827 + self.dimTags[self.geo.terminals.outer.name]
2828 + self.dimTags["innerTransitionNotch"]
2829 + self.dimTags["outerTransitionNotch"]
2830 + [(3, tag) for tag in contactLayerBetweenTerminalsAndWinding]
2831 )
2832 Pancake3DAllBoundary = gmsh.model.getBoundary(
2833 allPancakeVolumes, combined=True, oriented=False
2834 )
2835 Pancake3DBoundaryDimTags = list(
2836 set(Pancake3DAllBoundary)
2837 - set(topSurfaceDimTags)
2838 - set(bottomSurfaceDimTags)
2839 )
2840 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2841 self.regions.addEntities(
2842 self.geo.pancakeBoundaryName + "-OnlyWindingAndTerminals",
2843 pancake3DBoundaryTags,
2844 regionType.powered,
2845 entityType.surf,
2846 )
2848 # ==============================================================================
2849 # PANCAKE BOUNDARY SURFACE ENDS ================================================
2850 # ==============================================================================
2852 # Generate regions file:
2853 self.regions.generateRegionsFile(self.regions_file)
2854 self.rm = FilesAndFolders.read_data_from_yaml(self.regions_file, RegionsModel)
2856 logger.info(
2857 "Generating physical groups and regions file has been finished in"
2858 f" {timeit.default_timer() - start_time:.2f} s."
2859 )
2861 def generate_mesh_file(self):
2862 """
2863 Saves mesh file to disk.
2866 """
2867 logger.info(
2868 f"Generating Pancake3D mesh file ({self.mesh_file}) has been started."
2869 )
2870 start_time = timeit.default_timer()
2872 gmsh.write(self.mesh_file)
2874 logger.info(
2875 f"Generating Pancake3D mesh file ({self.mesh_file}) has been finished"
2876 f" in {timeit.default_timer() - start_time:.2f} s."
2877 )
2879 if self.mesh_gui:
2880 gmsh.option.setNumber("Geometry.Volumes", 0)
2881 gmsh.option.setNumber("Geometry.Surfaces", 0)
2882 gmsh.option.setNumber("Geometry.Curves", 0)
2883 gmsh.option.setNumber("Geometry.Points", 0)
2884 self.gu.launch_interactive_GUI()
2885 else:
2886 gmsh.clear()
2887 gmsh.finalize()
2889 def load_mesh(self):
2890 """
2891 Loads mesh from .msh file.
2894 """
2895 logger.info("Loading Pancake3D mesh has been started.")
2896 start_time = timeit.default_timer()
2898 previousGeo = FilesAndFolders.read_data_from_yaml(
2899 self.geometry_data_file, Pancake3DGeometry
2900 )
2901 previousMesh = FilesAndFolders.read_data_from_yaml(
2902 self.mesh_data_file, Pancake3DMesh
2903 )
2905 if previousGeo.model_dump() != self.geo.model_dump():
2906 raise ValueError(
2907 "Geometry data has been changed. Please regenerate the geometry or load"
2908 " the previous geometry data."
2909 )
2910 elif previousMesh.model_dump() != self.mesh.model_dump():
2911 raise ValueError(
2912 "Mesh data has been changed. Please regenerate the mesh or load the"
2913 " previous mesh data."
2914 )
2916 gmsh.clear()
2917 gmsh.open(self.mesh_file)
2919 logger.info(
2920 "Loading Pancake3D mesh has been finished in"
2921 f" {timeit.default_timer() - start_time:.2f} s."
2922 )