Coverage for fiqus/mesh_generators/MeshPancake3D.py: 92%
1021 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-12-25 02:54 +0100
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.dict())
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 # Mesh:
974 gmsh.model.mesh.generate()
975 gmsh.model.mesh.optimize()
977 # Print the log:
978 log = gmsh.logger.get()
979 for line in log:
980 if line.startswith("Info"):
981 logger.info(re.sub(r"Info:\s+", "", line))
982 elif line.startswith("Warning"):
983 logger.warning(re.sub(r"Warning:\s+", "", line))
985 gmsh.logger.stop()
986 except:
987 # Print the log:
988 log = gmsh.logger.get()
989 for line in log:
990 if line.startswith("Info"):
991 logger.info(re.sub(r"Info:\s+", "", line))
992 elif line.startswith("Warning"):
993 logger.warning(re.sub(r"Warning:\s+", "", line))
994 elif line.startswith("Error"):
995 logger.error(re.sub(r"Error:\s+", "", line))
997 gmsh.logger.stop()
999 self.generate_regions()
1001 logger.error(
1002 "Meshing Pancake3D magnet has failed. Try to change"
1003 " minimumElementSize and maximumElementSize parameters."
1004 )
1005 raise
1007 # SICN not implemented in 1D!
1008 allElementsDim2 = list(gmsh.model.mesh.getElements(dim=2)[1][0])
1009 allElementsDim3 = list(gmsh.model.mesh.getElements(dim=3)[1][0])
1010 allElements = allElementsDim2 + allElementsDim3
1011 elementQualities = gmsh.model.mesh.getElementQualities(allElements)
1012 lowestQuality = min(elementQualities)
1013 averageQuality = sum(elementQualities) / len(elementQualities)
1014 NofLowQualityElements = len(
1015 [quality for quality in elementQualities if quality < 0.01]
1016 )
1017 NofIllElemets = len(
1018 [quality for quality in elementQualities if quality < 0.001]
1019 )
1021 logger.info(
1022 f"The lowest quality among the elements is {lowestQuality:.4f} (SICN). The"
1023 " number of elements with quality lower than 0.01 is"
1024 f" {NofLowQualityElements}."
1025 )
1027 if NofIllElemets > 0:
1028 logger.warning(
1029 f"There are {NofIllElemets} elements with quality lower than 0.001. Try"
1030 " to change minimumElementSize and maximumElementSize parameters."
1031 )
1033 # Create cuts:
1034 # This is required to make the air a simply connected domain. This is required
1035 # for the solution part. You can read more about Homology in GMSH documentation.
1036 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name]]
1038 if self.geo.air.shellTransformation:
1039 shellTags = [
1040 dimTag[1] for dimTag in self.dimTags[self.geo.air.shellVolumeName]
1041 ]
1042 airTags.extend(shellTags)
1044 dummyAirRegion = gmsh.model.addPhysicalGroup(dim=3, tags=airTags)
1045 dummyAirRegionDimTag = (3, dummyAirRegion)
1047 innerCylinderTags = [self.dimTags[self.geo.air.name + "-InnerCylinder"][0][1]]
1048 gapTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name + "-Gap"]]
1049 # Only remove every second gap:
1050 gapTags = gapTags[1::2]
1052 dummyAirRegionWithoutInnerCylinder = gmsh.model.addPhysicalGroup(
1053 dim=3, tags=list(set(airTags) - set(innerCylinderTags) - set(gapTags))
1054 )
1055 dummyAirRegionWithoutInnerCylinderDimTag = (
1056 3,
1057 dummyAirRegionWithoutInnerCylinder,
1058 )
1060 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.winding.name]]
1061 dummyWindingRegion = gmsh.model.addPhysicalGroup(dim=3, tags=windingTags)
1062 dummyWindingRegionDimTag = (3, dummyWindingRegion)
1064 if self.geo.contactLayer.thinShellApproximation:
1065 # Find all the contact layer surfaces:
1066 allWindingDimTags = []
1067 for i in range(self.geo.numberOfPancakes):
1068 windingDimTags = self.dimTags[self.geo.winding.name + str(i + 1)]
1069 allWindingDimTags.extend(windingDimTags)
1071 windingBoundarySurfaces = gmsh.model.getBoundary(
1072 allWindingDimTags, combined=True, oriented=False
1073 )
1074 allWindingSurfaces = gmsh.model.getBoundary(
1075 allWindingDimTags, combined=False, oriented=False
1076 )
1078 contactLayerSurfacesDimTags = list(
1079 set(allWindingSurfaces) - set(windingBoundarySurfaces)
1080 )
1081 contactLayerTags = [dimTag[1] for dimTag in contactLayerSurfacesDimTags]
1083 # Get rid of non-contactLayer surfaces:
1084 realContactLayerTags = []
1085 for contactLayerTag in contactLayerTags:
1086 surfaceNormal = list(gmsh.model.getNormal(contactLayerTag, [0.5, 0.5]))
1087 centerOfMass = gmsh.model.occ.getCenterOfMass(2, contactLayerTag)
1089 if (
1090 abs(
1091 surfaceNormal[0] * centerOfMass[0]
1092 + surfaceNormal[1] * centerOfMass[1]
1093 )
1094 > 1e-6
1095 ):
1096 realContactLayerTags.append(contactLayerTag)
1098 # Get rid of surfaces that touch terminals:
1099 terminalSurfaces = gmsh.model.getBoundary(
1100 self.dimTags[self.geo.terminals.outer.name] + self.dimTags[self.geo.terminals.inner.name],
1101 combined=False,
1102 oriented=False,
1103 )
1104 terminalSurfaces = [dimTag[1] for dimTag in terminalSurfaces]
1105 finalContactLayerTags = [
1106 tag for tag in realContactLayerTags if tag not in terminalSurfaces
1107 ]
1109 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1110 dim=2, tags=finalContactLayerTags
1111 )
1112 dummyContactLayerRegionDimTag = (2, dummyContactLayerRegion)
1114 else:
1115 contactLayerTags = [dimTag[1] for dimTag in self.dimTags[self.geo.contactLayer.name]]
1117 # get rid of volumes that touch terminals:
1118 terminalSurfaces = gmsh.model.getBoundary(
1119 self.dimTags[self.geo.terminals.outer.name] + self.dimTags[self.geo.terminals.inner.name],
1120 combined=False,
1121 oriented=False,
1122 )
1123 finalContactLayerTags = []
1124 for contactLayerTag in contactLayerTags:
1125 insulatorSurfaces = gmsh.model.getBoundary(
1126 [(3, contactLayerTag)], combined=False, oriented=False
1127 )
1128 itTouchesTerminals = False
1129 for insulatorSurface in insulatorSurfaces:
1130 if insulatorSurface in terminalSurfaces:
1131 itTouchesTerminals = True
1132 break
1134 if not itTouchesTerminals:
1135 finalContactLayerTags.append(contactLayerTag)
1137 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1138 dim=3, tags=finalContactLayerTags
1139 )
1140 dummyContactLayerRegionDimTag = (3, dummyContactLayerRegion)
1142 # First cohomology request (normal cut for NI coils):
1143 gmsh.model.mesh.addHomologyRequest(
1144 "Cohomology",
1145 domainTags=[dummyAirRegion],
1146 subdomainTags=[],
1147 dims=[1],
1148 )
1150 if self.mesh.computeCohomologyForInsulating:
1151 # Second cohomology request (insulated cut for insulated coils):
1152 if self.geo.numberOfPancakes > 1:
1153 gmsh.model.mesh.addHomologyRequest(
1154 "Cohomology",
1155 domainTags=[
1156 dummyAirRegionWithoutInnerCylinder,
1157 dummyContactLayerRegion,
1158 ],
1159 subdomainTags=[],
1160 dims=[1],
1161 )
1162 else:
1163 gmsh.model.mesh.addHomologyRequest(
1164 "Cohomology",
1165 domainTags=[
1166 dummyAirRegion,
1167 dummyContactLayerRegion,
1168 ],
1169 subdomainTags=[],
1170 dims=[1],
1171 )
1173 # Third cohomology request (for cuts between pancake coils):
1174 gmsh.model.mesh.addHomologyRequest(
1175 "Cohomology",
1176 domainTags=[
1177 dummyAirRegion,
1178 dummyContactLayerRegion,
1179 dummyWindingRegion,
1180 ],
1181 subdomainTags=[],
1182 dims=[1],
1183 )
1185 # Start logger:
1186 gmsh.logger.start()
1188 cuts = gmsh.model.mesh.computeHomology()
1190 # Print the log:
1191 log = gmsh.logger.get()
1192 for line in log:
1193 if line.startswith("Info"):
1194 logger.info(re.sub(r"Info:\s+", "", line))
1195 elif line.startswith("Warning"):
1196 logger.warning(re.sub(r"Warning:\s+", "", line))
1197 gmsh.logger.stop()
1199 if self.geo.numberOfPancakes > 1:
1200 cutsDictionary = {
1201 "H^1{1}": self.geo.air.cutName,
1202 "H^1{1,4,3}": "CutsBetweenPancakes",
1203 "H^1{2,4}": "CutsForPerfectInsulation",
1204 }
1205 else:
1206 cutsDictionary = {
1207 "H^1{1}": self.geo.air.cutName,
1208 "H^1{1,4,3}": "CutsBetweenPancakes",
1209 "H^1{1,4}": "CutsForPerfectInsulation",
1210 }
1211 cutTags = [dimTag[1] for dimTag in cuts]
1212 cutEntities = []
1213 for tag in cutTags:
1214 name = gmsh.model.getPhysicalName(1, tag)
1215 cutEntities = list(gmsh.model.getEntitiesForPhysicalGroup(1, tag))
1216 cutEntitiesDimTags = [(1, cutEntity) for cutEntity in cutEntities]
1217 for key in cutsDictionary:
1218 if key in name:
1219 if cutsDictionary[key] in self.dimTags:
1220 self.dimTags[cutsDictionary[key]].extend(cutEntitiesDimTags)
1221 else:
1222 self.dimTags[cutsDictionary[key]] = cutEntitiesDimTags
1224 # Remove newly created physical groups because they will be created again in
1225 # generate_regions method.
1226 gmsh.model.removePhysicalGroups(
1227 [dummyContactLayerRegionDimTag]
1228 + [dummyAirRegionDimTag]
1229 + [dummyAirRegionWithoutInnerCylinderDimTag]
1230 + [dummyWindingRegionDimTag]
1231 + cuts
1232 )
1234 logger.info(
1235 "Generating Pancake3D mesh has been finished in"
1236 f" {timeit.default_timer() - start_time:.2f} s."
1237 )
1239 def structure_mesh(
1240 self,
1241 windingVolumeTags,
1242 windingSurfaceTags,
1243 windingLineTags,
1244 contactLayerVolumeTags,
1245 contactLayerSurfaceTags,
1246 contactLayerLineTags,
1247 meshSettingIndex,
1248 axialNumberOfElements=None,
1249 bumpCoefficient=None,
1250 ):
1251 """
1252 Structures the winding and contact layer meshed depending on the user inputs. If
1253 the bottom and top part of the air is to be structured, the same method is used.
1255 :param windingVolumeTags: tags of the winding volumes
1256 :type windingVolumeTags: list[int]
1257 :param windingSurfaceTags: tags of the winding surfaces
1258 :type windingSurfaceTags: list[int]
1259 :param windingLineTags: tags of the winding lines
1260 :type windingLineTags: list[int]
1261 :param contactLayerVolumeTags: tags of the contact layer volumes
1262 :type contactLayerVolumeTags: list[int]
1263 :param contactLayerSurfaceTags: tags of the contact layer surfaces
1264 :type contactLayerSurfaceTags: list[int]
1265 :param contactLayerLineTags: tags of the contact layer lines
1266 :type contactLayerLineTags: list[int]
1267 :param meshSettingIndex: index of the mesh setting
1268 :type meshSettingIndex: int
1269 :param axialNumberOfElements: number of axial elements
1270 :type axialNumberOfElements: int, optional
1271 :param bumpCoefficient: bump coefficient for axial meshing
1272 :type bumpCoefficient: float, optional
1274 """
1275 # Transfinite settings:
1276 # Arc lenght of the innermost one turn of spiral:
1277 if self.geo.contactLayer.thinShellApproximation:
1278 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1279 self.geo.winding.innerRadius,
1280 self.geo.winding.innerRadius
1281 + self.geo.winding.thickness
1282 + self.geo.contactLayer.thickness * (self.geo.numberOfPancakes - 1) / self.geo.numberOfPancakes,
1283 0,
1284 2 * math.pi,
1285 )
1286 else:
1287 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1288 self.geo.winding.innerRadius,
1289 self.geo.winding.innerRadius + self.geo.winding.thickness,
1290 0,
1291 2 * math.pi,
1292 )
1294 # Arc length of one element:
1295 arcElementLength = oneTurnSpiralLength / self.mesh.winding.azimuthalNumberOfElementsPerTurn[meshSettingIndex]
1297 # Number of azimuthal elements per turn:
1298 arcNumElementsPerTurn = round(oneTurnSpiralLength / arcElementLength)
1300 # Make all the lines transfinite:
1301 for j, lineTags in enumerate([windingLineTags, contactLayerLineTags]):
1302 for lineTag in lineTags:
1303 lineObject = curve(lineTag, self.geo)
1305 if lineObject.type is curveType.horizontal:
1306 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1307 # used.
1308 if self.geo.contactLayer.thinShellApproximation:
1309 numNodes = self.mesh.winding.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1311 else:
1312 if j == 0:
1313 # This line is the winding's horizontal line:
1314 numNodes = self.mesh.winding.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1316 else:
1317 # This line is the contact layer's horizontal line:
1318 numNodes = self.mesh.contactLayer.radialNumberOfElementsPerTurn[meshSettingIndex] + 1
1320 # Set transfinite curve:
1321 self.contactLayerAndWindingRadialLines.append(lineTag)
1322 gmsh.model.mesh.setTransfiniteCurve(lineTag, numNodes)
1324 elif lineObject.type is curveType.axial:
1325 # The curve is axial, so axialNumberOfElements entry is used.
1326 if axialNumberOfElements is None:
1327 numNodes = self.mesh.winding.axialNumberOfElements[meshSettingIndex] + 1
1328 else:
1329 numNodes = axialNumberOfElements + 1
1331 if bumpCoefficient is None:
1332 bumpCoefficient = self.mesh.winding.axialDistributionCoefficient[meshSettingIndex]
1333 gmsh.model.mesh.setTransfiniteCurve(
1334 lineTag, numNodes, meshType="Bump", coef=bumpCoefficient
1335 )
1337 else:
1338 # The line is an arc, so the previously calculated arcNumElementsPerTurn
1339 # is used. All the number of elements per turn must be the same
1340 # independent of radial position. Otherwise, transfinite meshing cannot
1341 # be performed. However, to support the float number of turns, number
1342 # of nodes are being calculated depending on the start and end turns of
1343 # the arc.d
1344 lengthInTurns = abs(lineObject.n2 - lineObject.n1)
1345 if lengthInTurns > 0.5:
1346 # The arc can never be longer than half a turn.
1347 lengthInTurns = 1 - lengthInTurns
1349 lengthInTurns = (
1350 round(lengthInTurns / self.geo.winding.turnTol) * self.geo.winding.turnTol
1351 )
1353 arcNumEl = round(arcNumElementsPerTurn * lengthInTurns)
1355 arcNumNodes = int(arcNumEl + 1)
1357 # Set transfinite curve:
1358 gmsh.model.mesh.setTransfiniteCurve(lineTag, arcNumNodes)
1360 for j, surfTags in enumerate([windingSurfaceTags, contactLayerSurfaceTags]):
1361 for surfTag in surfTags:
1362 # Make all the surfaces transfinite:
1363 gmsh.model.mesh.setTransfiniteSurface(surfTag)
1365 if self.mesh.winding.elementType[meshSettingIndex] == "hexahedron":
1366 # If the element type is hexahedron, recombine all the surfaces:
1367 gmsh.model.mesh.setRecombine(2, surfTag)
1368 elif self.mesh.winding.elementType[meshSettingIndex] == "prism":
1369 # If the element type is prism, recombine only the side surfaces:
1370 surfaceNormal = list(gmsh.model.getNormal(surfTag, [0.5, 0.5]))
1371 if abs(surfaceNormal[2]) < 1e-6:
1372 gmsh.model.mesh.setRecombine(2, surfTag)
1374 # If the element type is tetrahedron, do not recombine any surface.
1376 for volTag in windingVolumeTags + contactLayerVolumeTags:
1377 # Make all the volumes transfinite:
1378 gmsh.model.mesh.setTransfiniteVolume(volTag)
1380 def structure_tubes_and_cylinders(
1381 self, volumeTags, terminalNonTubeParts=False, radialElementMultiplier=1
1382 ):
1383 # Number of azimuthal elements per quarter:
1384 arcNumElementsPerQuarter = int(self.mesh.winding.azimuthalNumberOfElementsPerTurn[0] / 4)
1385 radialNumberOfElementsPerLength = (
1386 self.mesh.winding.radialNumberOfElementsPerTurn[0] / self.geo.winding.thickness * radialElementMultiplier
1387 )
1389 surfacesDimTags = gmsh.model.getBoundary(
1390 [(3, tag) for tag in volumeTags], combined=False, oriented=False
1391 )
1392 surfacesTags = [dimTag[1] for dimTag in surfacesDimTags]
1393 surfacesTags = list(set(surfacesTags))
1395 curvesDimTags = gmsh.model.getBoundary(
1396 surfacesDimTags, combined=False, oriented=False
1397 )
1398 curvesTags = [dimTag[1] for dimTag in curvesDimTags]
1400 # Make all the lines transfinite:
1401 for curveTag in curvesTags:
1402 curveObject = curve(curveTag, self.geo)
1404 if curveObject.type is curveType.horizontal:
1405 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1406 # used.
1408 # But, the curve might be a part of the transitionNotch.
1409 isTransitionNotch = False
1410 point2 = curveObject.points[1]
1411 point1 = curveObject.points[0]
1412 if (
1413 abs(point2[0] - point1[0]) > 1e-5
1414 and abs(point2[1] - point1[1]) > 1e-5
1415 ):
1416 isTransitionNotch = True
1418 if isTransitionNotch:
1419 gmsh.model.mesh.setTransfiniteCurve(curveTag, 2)
1420 else:
1421 if terminalNonTubeParts:
1422 if curveTag not in self.contactLayerAndWindingRadialLines:
1423 numNodes = (
1424 round(radialNumberOfElementsPerLength * self.geo.terminals.inner.thickness)
1425 + 1
1426 )
1427 # Set transfinite curve:
1428 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1429 else:
1430 numNodes = (
1431 round(radialNumberOfElementsPerLength * curveObject.length)
1432 + 1
1433 )
1434 # Set transfinite curve:
1435 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1437 elif curveObject.type is curveType.axial:
1438 # The curve is axial, so axialNumberOfElements entry is used.
1439 if math.isclose(curveObject.length, self.geo.winding.height, rel_tol=1e-7):
1440 numNodes = min(self.mesh.winding.axialNumberOfElements) + 1
1441 else:
1442 axialElementsPerLength = min(self.mesh.winding.axialNumberOfElements) / self.geo.winding.height
1443 numNodes = (
1444 round(axialElementsPerLength * curveObject.length + 1e-6) + 1
1445 )
1447 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1449 else:
1450 # The line is an arc
1451 lengthInTurns = abs(curveObject.n2 - curveObject.n1)
1452 if lengthInTurns > 0.5:
1453 # The arc can never be longer than half a turn.
1454 lengthInTurns = 1 - lengthInTurns
1456 lengthInTurns = (
1457 round(lengthInTurns / self.geo.winding.turnTol) * self.geo.winding.turnTol
1458 )
1460 arcNumEl = round(arcNumElementsPerQuarter * 4 * lengthInTurns)
1462 arcNumNodes = int(arcNumEl + 1)
1464 # Set transfinite curve:
1466 gmsh.model.mesh.setTransfiniteCurve(curveTag, arcNumNodes)
1468 for surfaceTag in surfacesTags:
1469 # Make all the surfaces transfinite:
1471 if self.mesh.winding.elementType[0] == "hexahedron":
1472 # If the element type is hexahedron, recombine all the surfaces:
1473 gmsh.model.mesh.setRecombine(2, surfaceTag)
1474 elif self.mesh.winding.elementType[0] == "prism":
1475 # If the element type is prism, recombine only the side surfaces:
1476 surfaceNormal = list(gmsh.model.getNormal(surfaceTag, [0.5, 0.5]))
1477 if abs(surfaceNormal[2]) < 1e-5:
1478 gmsh.model.mesh.setRecombine(2, surfaceTag)
1480 curves = gmsh.model.getBoundary(
1481 [(2, surfaceTag)], combined=False, oriented=False
1482 )
1483 numberOfCurves = len(curves)
1484 if terminalNonTubeParts:
1485 if numberOfCurves == 4:
1486 numberOfHorizontalCurves = 0
1487 for curveTag in curves:
1488 curveObject = curve(curveTag[1], self.geo)
1489 if curveObject.type is curveType.horizontal:
1490 numberOfHorizontalCurves += 1
1492 if numberOfHorizontalCurves == 3:
1493 pass
1494 else:
1495 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1497 elif numberOfCurves == 3:
1498 pass
1499 else:
1500 curves = gmsh.model.getBoundary(
1501 [(2, surfaceTag)], combined=False, oriented=False
1502 )
1503 curveObjects = [curve(line[1], self.geo) for line in curves]
1505 divisionCurves = []
1506 for curveObject in curveObjects:
1507 if curveObject.type is curveType.horizontal:
1508 point1 = curveObject.points[0]
1509 point2 = curveObject.points[1]
1510 if not (
1511 abs(point2[0] - point1[0]) > 1e-5
1512 and abs(point2[1] - point1[1]) > 1e-5
1513 ):
1514 divisionCurves.append(curveObject)
1516 cornerPoints = (
1517 divisionCurves[0].pointTags + divisionCurves[1].pointTags
1518 )
1520 if surfaceTag not in alreadyMeshedSurfaceTags:
1521 alreadyMeshedSurfaceTags.append(surfaceTag)
1522 gmsh.model.mesh.setTransfiniteSurface(
1523 surfaceTag, cornerTags=cornerPoints
1524 )
1525 else:
1526 if numberOfCurves == 3:
1527 # Then it is a pie, corner points should be adjusted:
1528 originPointTag = None
1529 curveObject1 = curve(curves[0][1], self.geo)
1530 for point, tag in zip(curveObject1.points, curveObject1.pointTags):
1531 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1532 originPointTag = tag
1534 if originPointTag is None:
1535 curveObject2 = curve(curves[1][1], self.geo)
1536 for point, tag in zip(
1537 curveObject2.points, curveObject2.pointTags
1538 ):
1539 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1540 originPointTag = tag
1542 otherPointDimTags = gmsh.model.getBoundary(
1543 [(2, surfaceTag)],
1544 combined=False,
1545 oriented=False,
1546 recursive=True,
1547 )
1548 otherPointTags = [dimTag[1] for dimTag in otherPointDimTags]
1549 otherPointTags.remove(originPointTag)
1550 cornerTags = [originPointTag] + otherPointTags
1551 gmsh.model.mesh.setTransfiniteSurface(
1552 surfaceTag, cornerTags=cornerTags
1553 )
1554 else:
1555 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1557 for volumeTag in volumeTags:
1558 if terminalNonTubeParts:
1559 surfaces = gmsh.model.getBoundary(
1560 [(3, volumeTag)], combined=False, oriented=False
1561 )
1562 curves = gmsh.model.getBoundary(
1563 surfaces, combined=False, oriented=False
1564 )
1565 curves = list(set(curves))
1567 if len(curves) == 12:
1568 numberOfArcs = 0
1569 for curveTag in curves:
1570 curveObject = curve(curveTag[1], self.geo)
1571 if curveObject.type is curveType.spiralArc:
1572 numberOfArcs += 1
1573 if numberOfArcs == 2:
1574 pass
1575 else:
1576 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1577 # elif len(curves) == 15:
1578 # divisionCurves = []
1579 # for curveTag in curves:
1580 # curveObject = curve(curveTag[1], self.geo)
1581 # if curveObject.type is curveType.horizontal:
1582 # point1 = curveObject.points[0]
1583 # point2 = curveObject.points[1]
1584 # if not (
1585 # abs(point2[0] - point1[0]) > 1e-5
1586 # and abs(point2[1] - point1[1]) > 1e-5
1587 # ):
1588 # divisionCurves.append(curveObject)
1590 # cornerPoints = (
1591 # divisionCurves[0].pointTags
1592 # + divisionCurves[1].pointTags
1593 # + divisionCurves[2].pointTags
1594 # + divisionCurves[3].pointTags
1595 # )
1596 # gmsh.model.mesh.setTransfiniteVolume(
1597 # volumeTag, cornerTags=cornerPoints
1598 # )
1599 else:
1600 # Make all the volumes transfinite:
1601 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1603 @staticmethod
1604 def get_boundaries(volumeDimTags, returnTags=False):
1605 """
1606 Returns all the surface and line dimTags or tags of a given list of volume
1607 dimTags.
1609 :param volumeDimTags: dimTags of the volumes
1610 :type volumeDimTags: list[tuple[int, int]]
1611 :param returnTags: if True, returns tags instead of dimTags
1612 :type returnTags: bool, optional
1613 :return: surface and line dimTags or tags
1614 :rtype: tuple[list[tuple[int, int]], list[tuple[int, int]]] or
1615 tuple[list[int], list[int]]
1616 """
1617 # Get the surface tags:
1618 surfaceDimTags = list(
1619 set(
1620 gmsh.model.getBoundary(
1621 volumeDimTags,
1622 combined=False,
1623 oriented=False,
1624 recursive=False,
1625 )
1626 )
1627 )
1629 # Get the line tags:
1630 lineDimTags = list(
1631 set(
1632 gmsh.model.getBoundary(
1633 surfaceDimTags,
1634 combined=False,
1635 oriented=False,
1636 recursive=False,
1637 )
1638 )
1639 )
1641 if returnTags:
1642 surfaceTags = [dimTag[1] for dimTag in surfaceDimTags]
1643 lineTags = [dimTag[1] for dimTag in lineDimTags]
1644 return surfaceTags, lineTags
1645 else:
1646 return surfaceDimTags, lineDimTags
1648 @staticmethod
1649 def fuse_volumes(volumeDimTags, fuseSurfaces=True, fusedSurfacesArePlane=False):
1650 """
1651 Fuses all the volumes in a given list of volume dimTags, removes old volumes,
1652 and returns the new volume dimTag. Also, if compundSurfacces is True, it fuses
1653 the surfaces that only belong to the volume. fusedSurfacesArePlane can be
1654 used to change the behavior of the fuse_surfaces method.
1656 :param volumeDimTags: dimTags of the volumes
1657 :type volumeDimTags: list[tuple[int, int]]
1658 :param fuseSurfaces: if True, fuses the surfaces that only belong to the
1659 volume
1660 :type fuseSurfaces: bool, optional
1661 :param fusedSurfacesArePlane: if True, fused surfaces are assumed to be
1662 plane, and fusion is performed accordingly
1663 :return: new volume's dimTag
1664 :rtype: tuple[int, int]
1665 """
1667 # Get the combined boundary surfaces:
1668 boundarySurfacesDimTags = gmsh.model.getBoundary(
1669 volumeDimTags,
1670 combined=True,
1671 oriented=False,
1672 recursive=False,
1673 )
1674 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1676 # Get all the boundary surfaces:
1677 allBoundarySurfacesDimTags = gmsh.model.getBoundary(
1678 volumeDimTags,
1679 combined=False,
1680 oriented=False,
1681 recursive=False,
1682 )
1684 # Find internal (common) surfaces:
1685 internalSurfacesDimTags = list(
1686 set(allBoundarySurfacesDimTags) - set(boundarySurfacesDimTags)
1687 )
1689 # Get the combined boundary lines:
1690 boundaryLinesDimTags = gmsh.model.getBoundary(
1691 allBoundarySurfacesDimTags,
1692 combined=True,
1693 oriented=False,
1694 recursive=False,
1695 )
1696 boundarLinesTags = [dimTag[1] for dimTag in boundaryLinesDimTags]
1698 # Get all the boundary lines:
1699 allBoundaryLinesDimTags = gmsh.model.getBoundary(
1700 allBoundarySurfacesDimTags,
1701 combined=False,
1702 oriented=False,
1703 recursive=False,
1704 )
1706 # Find internal (common) lines:
1707 internalLinesDimTags = list(
1708 set(allBoundaryLinesDimTags) - set(boundarLinesTags)
1709 )
1711 # Remove the old volumes:
1712 removedVolumeDimTags = volumeDimTags
1713 gmsh.model.occ.remove(removedVolumeDimTags, recursive=False)
1715 # Remove the internal surfaces:
1716 gmsh.model.occ.remove(internalSurfacesDimTags, recursive=False)
1718 # Remove the internal lines:
1719 gmsh.model.occ.remove(internalLinesDimTags, recursive=False)
1721 # Create a new single volume (even thought we don't use the new volume tag
1722 # directly, it is required for finding the surfaces that only belong to the
1723 # volume):
1724 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1725 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1726 newVolumeDimTag = (3, newVolumeTag)
1727 gmsh.model.occ.synchronize()
1729 if fuseSurfaces:
1730 newVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1731 (3, newVolumeTag), surfacesArePlane=fusedSurfacesArePlane
1732 )
1734 return newVolumeDimTag
1736 @staticmethod
1737 def fuse_common_surfaces_of_two_volumes(
1738 volume1DimTags, volume2DimTags, fuseOtherSurfaces=False, surfacesArePlane=False
1739 ):
1740 """
1741 Fuses common surfaces of two volumes. Volumes are given as a list of dimTags,
1742 but they are assumed to form a single volume, and this function fuses those
1743 multiple volumes into a single volume as well. If fuseOtherSurfaces is set to
1744 True, it tries to fuse surfaces that only belong to one volume too; however,
1745 that feature is not used in Pancake3D currently.
1747 :param volume1DimTags: dimTags of the first volume
1748 :type volume1DimTags: list[tuple[int, int]]
1749 :param volume2DimTags: dimTags of the second volume
1750 :type volume2DimTags: list[tuple[int, int]]
1751 :param fuseOtherSurfaces: if True, fuses the surfaces that only belong to one
1752 volume
1753 :type fuseOtherSurfaces: bool, optional
1754 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1755 fusion is performed accordingly
1756 :type surfacesArePlane: bool, optional
1757 :return: new volumes dimTags
1758 :rtype: tuple[tuple[int, int], tuple[int, int]]
1759 """
1760 vol1BoundarySurfacesDimTags = gmsh.model.getBoundary(
1761 volume1DimTags,
1762 combined=True,
1763 oriented=False,
1764 recursive=False,
1765 )
1767 vol2BoundarySurfacesDimTags = gmsh.model.getBoundary(
1768 volume2DimTags,
1769 combined=True,
1770 oriented=False,
1771 recursive=False,
1772 )
1774 # Remove the old volumes:
1775 gmsh.model.occ.remove(volume1DimTags + volume2DimTags, recursive=False)
1777 # Find common surfaces:
1778 commonSurfacesDimTags = list(
1779 set(vol2BoundarySurfacesDimTags).intersection(
1780 set(vol1BoundarySurfacesDimTags)
1781 )
1782 )
1784 # Fuse common surfaces:
1785 fusedCommonSurfaceDimTags = Mesh.fuse_surfaces(
1786 commonSurfacesDimTags, surfacesArePlane=surfacesArePlane
1787 )
1789 # Create the new volumes:
1790 for commonSurfaceDimTag in commonSurfacesDimTags:
1791 vol1BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1792 vol2BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1794 vol1BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1795 vol1BoundarySurfaceTags = [dimTag[1] for dimTag in vol1BoundarySurfacesDimTags]
1796 vol2BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1797 vol2BoundarySurfaceTags = [dimTag[1] for dimTag in vol2BoundarySurfacesDimTags]
1799 vol1SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1800 vol1BoundarySurfaceTags, sewing=True
1801 )
1802 vol1NewVolumeDimTag = (3, gmsh.model.occ.addVolume([vol1SurfaceLoop]))
1804 vol2SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1805 vol2BoundarySurfaceTags, sewing=True
1806 )
1807 vol2NewVolumeDimTag = (
1808 3,
1809 gmsh.model.occ.addVolume([vol2SurfaceLoop]),
1810 )
1812 gmsh.model.occ.synchronize()
1814 if fuseOtherSurfaces:
1815 vol1NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1816 vol1NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1817 )
1818 vol2NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1819 vol2NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1820 )
1822 return vol1NewVolumeDimTag, vol2NewVolumeDimTag
1824 @staticmethod
1825 def fuse_possible_surfaces_of_a_volume(volumeDimTag, surfacesArePlane=False):
1826 """
1827 Fuses surfaces that only belong to the volumeDimTag.
1829 :param volumeDimTag: dimTag of the volume
1830 :type volumeDimTag: tuple[int, int]
1831 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1832 fusion is performed accordingly
1833 :type surfacesArePlane: bool, optional
1834 :return: new volume dimTag
1835 :rtype: tuple[int, int]
1836 """
1837 boundarySurfacesDimTags = gmsh.model.getBoundary(
1838 [volumeDimTag],
1839 combined=True,
1840 oriented=False,
1841 recursive=False,
1842 )
1843 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1845 # Combine surfaces that only belong to the volume:
1846 toBeFusedSurfacesDimTags = []
1847 surfacesNormals = []
1848 for surfaceDimTag in boundarySurfacesDimTags:
1849 upward, _ = gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1851 if len(list(upward)) == 1:
1852 toBeFusedSurfacesDimTags.append(surfaceDimTag)
1853 # Get the normal of the surface:
1854 surfacesNormals.append(
1855 list(gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5]))
1856 )
1858 # Remove the old volume (it is not required anymore):
1859 gmsh.model.occ.remove([volumeDimTag], recursive=False)
1860 gmsh.model.occ.synchronize()
1862 # Categorize surfaces based on their normals so that they can be combined
1863 # correctly. Without these, perpendicular surfaces will cause problems.
1865 # Define a threshold to determine if two surface normals are similar or not
1866 threshold = 1e-6
1868 # Initialize an empty list to store the sets of surfaces
1869 setsOfSurfaces = []
1871 # Calculate the Euclidean distance between each pair of objects
1872 for i in range(len(toBeFusedSurfacesDimTags)):
1873 surfaceDimTag = toBeFusedSurfacesDimTags[i]
1874 surfaceTouchingVolumeTags, _ = list(
1875 gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1876 )
1877 surfaceNormal = surfacesNormals[i]
1878 assignedToASet = False
1880 for surfaceSet in setsOfSurfaces:
1881 representativeSurfaceDimTag = surfaceSet[0]
1882 representativeSurfaceTouchingVolumeTags, _ = list(
1883 gmsh.model.getAdjacencies(
1884 representativeSurfaceDimTag[0],
1885 representativeSurfaceDimTag[1],
1886 )
1887 )
1888 representativeNormal = list(
1889 gmsh.model.getNormal(representativeSurfaceDimTag[1], [0.5, 0.5])
1890 )
1892 # Calculate the difference between surfaceNormal and
1893 # representativeNormal:
1894 difference = math.sqrt(
1895 sum(
1896 (x - y) ** 2
1897 for x, y in zip(surfaceNormal, representativeNormal)
1898 )
1899 )
1901 # Check if the distance is below the threshold
1902 if difference < threshold and set(surfaceTouchingVolumeTags) == set(
1903 representativeSurfaceTouchingVolumeTags
1904 ):
1905 # Add the object to an existing category
1906 surfaceSet.append(surfaceDimTag)
1907 assignedToASet = True
1908 break
1910 if not assignedToASet:
1911 # Create a new category with the current object if none of the
1912 # existing sets match
1913 setsOfSurfaces.append([surfaceDimTag])
1915 for surfaceSet in setsOfSurfaces:
1916 if len(surfaceSet) > 1:
1917 oldSurfaceDimTags = surfaceSet
1918 newSurfaceDimTags = Mesh.fuse_surfaces(
1919 oldSurfaceDimTags, surfacesArePlane=surfacesArePlane
1920 )
1921 newSurfaceTags = [dimTag[1] for dimTag in newSurfaceDimTags]
1923 oldSurfaceTags = [dimTag[1] for dimTag in oldSurfaceDimTags]
1924 boundarSurfacesTags = [
1925 tag for tag in boundarSurfacesTags if tag not in oldSurfaceTags
1926 ]
1927 boundarSurfacesTags.extend(newSurfaceTags)
1929 # Create a new single volume:
1930 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1931 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1932 gmsh.model.occ.synchronize()
1934 return (3, newVolumeTag)
1936 @staticmethod
1937 def fuse_surfaces(surfaceDimTags, surfacesArePlane=False, categorizeSurfaces=False):
1938 """
1939 Fuses all the surfaces in a given list of surface dimTags, removes the old
1940 surfaces, and returns the new dimTags. If surfacesArePlane is True, the surfaces
1941 are assumed to be plane, and fusing will be done without gmsh.model.occ.fuse
1942 method, which is faster.
1944 :param surfaceDimTags: dimTags of the surfaces
1945 :type surfaceDimTags: list[tuple[int, int]]
1946 :param surfacesArePlane: if True, surfaces are assumed to be plane
1947 :type surfacesArePlane: bool, optional
1948 :return: newly created surface dimTags
1949 :rtype: list[tuple[int, int]]
1950 """
1951 oldSurfaceDimTags = surfaceDimTags
1953 if surfacesArePlane:
1954 # Get the combined boundary curves:
1955 boundaryCurvesDimTags = gmsh.model.getBoundary(
1956 oldSurfaceDimTags,
1957 combined=True,
1958 oriented=False,
1959 recursive=False,
1960 )
1962 # Get all the boundary curves:
1963 allCurvesDimTags = gmsh.model.getBoundary(
1964 oldSurfaceDimTags,
1965 combined=False,
1966 oriented=False,
1967 recursive=False,
1968 )
1970 # Find internal (common) curves:
1971 internalCurvesDimTags = list(
1972 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
1973 )
1975 # Remove the old surfaces:
1976 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
1978 # Remove the internal curves:
1979 gmsh.model.occ.remove(internalCurvesDimTags, recursive=True)
1981 # Create a new single surface:
1982 def findOuterOnes(dimTags, findInnerOnes=False):
1983 """
1984 Finds the outermost surface/curve/point in a list of dimTags. The outermost means
1985 the furthest from the origin.
1986 """
1987 dim = dimTags[0][0]
1989 if dim == 2:
1990 distances = []
1991 for dimTag in dimTags:
1992 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])
1993 for curve in curves:
1994 curve = list(curve)
1995 gmsh.model.occ.synchronize()
1996 pointTags = gmsh.model.getBoundary(
1997 [(1, curveTag) for curveTag in curve],
1998 oriented=False,
1999 combined=False,
2000 )
2001 # Get the positions of the points:
2002 points = []
2003 for dimTag in pointTags:
2004 boundingbox1 = gmsh.model.occ.getBoundingBox(
2005 0, dimTag[1]
2006 )[:3]
2007 boundingbox2 = gmsh.model.occ.getBoundingBox(
2008 0, dimTag[1]
2009 )[3:]
2010 boundingbox = list(
2011 map(operator.add, boundingbox1, boundingbox2)
2012 )
2013 points.append(
2014 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2015 )
2017 distances.append(
2018 max([point[0] ** 2 + point[1] ** 2 for point in points])
2019 )
2020 elif dim == 1:
2021 distances = []
2022 for dimTag in dimTags:
2023 gmsh.model.occ.synchronize()
2024 pointTags = gmsh.model.getBoundary(
2025 [dimTag],
2026 oriented=False,
2027 combined=False,
2028 )
2029 # Get the positions of the points:
2030 points = []
2031 for dimTag in pointTags:
2032 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2033 :3
2034 ]
2035 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2036 3:
2037 ]
2038 boundingbox = list(
2039 map(operator.add, boundingbox1, boundingbox2)
2040 )
2041 points.append(
2042 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2043 )
2045 distances.append(
2046 max([point[0] ** 2 + point[1] ** 2 for point in points])
2047 )
2049 if findInnerOnes:
2050 goalDistance = min(distances)
2051 else:
2052 goalDistance = max(distances)
2054 result = []
2055 for distance, dimTag in zip(distances, dimTags):
2056 # Return all the dimTags with the hoal distance:
2057 if math.isclose(distance, goalDistance, abs_tol=1e-6):
2058 result.append(dimTag)
2060 return result
2062 outerCurvesDimTags = findOuterOnes(boundaryCurvesDimTags)
2063 outerCurvesTags = [dimTag[1] for dimTag in outerCurvesDimTags]
2064 curveLoopOuter = gmsh.model.occ.addCurveLoop(outerCurvesTags)
2066 innerCurvesDimTags = findOuterOnes(
2067 boundaryCurvesDimTags, findInnerOnes=True
2068 )
2069 innerCurvesTags = [dimTag[1] for dimTag in innerCurvesDimTags]
2070 curveLoopInner = gmsh.model.occ.addCurveLoop(innerCurvesTags)
2072 newSurfaceTag = gmsh.model.occ.addPlaneSurface(
2073 [curveLoopOuter, curveLoopInner]
2074 )
2076 gmsh.model.occ.synchronize()
2078 return [(2, newSurfaceTag)]
2079 else:
2080 # Create a new single surface:
2081 # The order of tags seems to be important for the fuse method to work
2082 # and not crash with a segmentation fault.
2083 try:
2084 fuseResults = gmsh.model.occ.fuse(
2085 [oldSurfaceDimTags[0]],
2086 oldSurfaceDimTags[1:],
2087 removeObject=False,
2088 removeTool=False,
2089 )
2090 newSurfaceDimTags = fuseResults[0]
2091 except:
2092 return oldSurfaceDimTags
2094 # Get the combined boundary curves:
2095 gmsh.model.occ.synchronize()
2096 boundaryCurvesDimTags = gmsh.model.getBoundary(
2097 newSurfaceDimTags,
2098 combined=True,
2099 oriented=False,
2100 recursive=False,
2101 )
2103 # Get all the boundary curves:
2104 allCurvesDimTags = gmsh.model.getBoundary(
2105 oldSurfaceDimTags,
2106 combined=False,
2107 oriented=False,
2108 recursive=False,
2109 )
2111 # Find internal (common) curves:
2112 internalCurvesDimTags = list(
2113 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
2114 )
2116 # Remove the old surfaces:
2117 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
2119 # Remove the internal curves:
2120 gmsh.model.occ.remove(internalCurvesDimTags, recursive=False)
2122 gmsh.model.occ.synchronize()
2124 return newSurfaceDimTags
2126 def generate_regions(self):
2127 """
2128 Generates physical groups and the regions file. Physical groups are generated in
2129 GMSH, and their tags and names are saved in the regions file. FiQuS use the
2130 regions file to create the corresponding .pro file.
2132 .vi file sends the information about geometry from geometry class to mesh class.
2133 .regions file sends the information about the physical groups formed out of
2134 elementary entities from the mesh class to the solution class.
2136 The file extension for the regions file is custom because users should not edit
2137 or even see this file.
2139 Regions are generated in the meshing part because BREP files cannot store
2140 regions.
2141 """
2142 logger.info("Generating physical groups and regions file has been started.")
2143 start_time = timeit.default_timer()
2145 # Create regions instance to both generate regions file and physical groups:
2146 self.regions = regions()
2148 # ==============================================================================
2149 # WINDING AND CONTACT LAYER REGIONS START =========================================
2150 # ==============================================================================
2151 if not self.geo.contactLayer.thinShellApproximation:
2152 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.winding.name]]
2153 self.regions.addEntities(
2154 self.geo.winding.name, windingTags, regionType.powered, entityType.vol
2155 )
2157 insulatorTags = [dimTag[1] for dimTag in self.dimTags[self.geo.contactLayer.name]]
2159 terminalDimTags = (
2160 self.dimTags[self.geo.terminals.inner.name] + self.dimTags[self.geo.terminals.outer.name]
2161 )
2162 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2163 terminalDimTags, combined=False, oriented=False
2164 )
2165 transitionNotchSurfaces = gmsh.model.getBoundary(
2166 self.dimTags["innerTransitionNotch"]
2167 + self.dimTags["outerTransitionNotch"],
2168 combined=False,
2169 oriented=False,
2170 )
2172 contactLayer = []
2173 contactLayerBetweenTerminalsAndWinding = []
2174 for insulatorTag in insulatorTags:
2175 insulatorSurfaces = gmsh.model.getBoundary(
2176 [(3, insulatorTag)], combined=False, oriented=False
2177 )
2178 itTouchesTerminals = False
2179 for insulatorSurface in insulatorSurfaces:
2180 if (
2181 insulatorSurface
2182 in terminalAndNotchSurfaces + transitionNotchSurfaces
2183 ):
2184 contactLayerBetweenTerminalsAndWinding.append(insulatorTag)
2185 itTouchesTerminals = True
2186 break
2188 if not itTouchesTerminals:
2189 contactLayer.append(insulatorTag)
2191 self.regions.addEntities(
2192 self.geo.contactLayer.name, contactLayer, regionType.insulator, entityType.vol
2193 )
2195 self.regions.addEntities(
2196 "WindingAndTerminalContactLayer",
2197 contactLayerBetweenTerminalsAndWinding,
2198 regionType.insulator,
2199 entityType.vol,
2200 )
2201 else:
2202 # Calculate the number of stacks for each individual winding. Number of
2203 # stacks is the number of volumes per turn. It affects how the regions
2204 # are created because of the TSA's pro file formulation.
2206 # find the smallest prime number that divides NofVolumes:
2207 windingDimTags = self.dimTags[self.geo.winding.name + "1"]
2208 windingTags = [dimTag[1] for dimTag in windingDimTags]
2209 NofVolumes = self.geo.winding.numberOfVolumesPerTurn
2210 smallest_prime_divisor = 2
2211 while NofVolumes % smallest_prime_divisor != 0:
2212 smallest_prime_divisor += 1
2214 # the number of stacks is the region divison per turn:
2215 NofStacks = smallest_prime_divisor
2217 # the number of sets are the total number of regions for all windings and
2218 # contact layers:
2219 NofSets = 2 * NofStacks
2221 allInnerTerminalSurfaces = gmsh.model.getBoundary(
2222 self.dimTags[self.geo.terminals.inner.name] + self.dimTags["innerTransitionNotch"],
2223 combined=False,
2224 oriented=False,
2225 )
2226 allInnerTerminalContactLayerSurfaces = []
2227 for innerTerminalSurface in allInnerTerminalSurfaces:
2228 normal = gmsh.model.getNormal(innerTerminalSurface[1], [0.5, 0.5])
2229 if abs(normal[2]) < 1e-5:
2230 curves = gmsh.model.getBoundary(
2231 [innerTerminalSurface], combined=False, oriented=False
2232 )
2233 curveTags = [dimTag[1] for dimTag in curves]
2234 for curveTag in curveTags:
2235 curveObject = curve(curveTag, self.geo)
2236 if curveObject.type is curveType.spiralArc:
2237 allInnerTerminalContactLayerSurfaces.append(
2238 innerTerminalSurface[1]
2239 )
2241 finalWindingSets = []
2242 finalContactLayerSets = []
2243 for i in range(NofSets):
2244 finalWindingSets.append([])
2245 finalContactLayerSets.append([])
2247 for i in range(self.geo.numberOfPancakes):
2248 windingDimTags = self.dimTags[self.geo.winding.name + str(i + 1)]
2249 windingTags = [dimTag[1] for dimTag in windingDimTags]
2251 NofVolumes = len(windingDimTags)
2253 windings = []
2254 for windingTag in windingTags:
2255 surfaces = gmsh.model.getBoundary(
2256 [(3, windingTag)], combined=False, oriented=False
2257 )
2258 curves = gmsh.model.getBoundary(
2259 surfaces, combined=False, oriented=False
2260 )
2261 curveTags = list(set([dimTag[1] for dimTag in curves]))
2262 for curveTag in curveTags:
2263 curveObject = curve(curveTag, self.geo)
2264 if curveObject.type is curveType.spiralArc:
2265 windingVolumeLengthInTurns = abs(
2266 curveObject.n2 - curveObject.n1
2267 )
2268 if windingVolumeLengthInTurns > 0.5:
2269 # The arc can never be longer than half a turn.
2270 windingVolumeLengthInTurns = (
2271 1 - windingVolumeLengthInTurns
2272 )
2274 windings.append((windingTag, windingVolumeLengthInTurns))
2276 windingStacks = []
2277 while len(windings) > 0:
2278 stack = []
2279 stackLength = 0
2280 for windingTag, windingVolumeLengthInTurns in windings:
2281 if stackLength < 1 / NofStacks - 1e-6:
2282 stack.append(windingTag)
2283 stackLength += windingVolumeLengthInTurns
2284 else:
2285 break
2286 # remove all the windings that are already added to the stack:
2287 windings = [
2288 (windingTag, windingVolumeLengthInTurns)
2289 for windingTag, windingVolumeLengthInTurns in windings
2290 if windingTag not in stack
2291 ]
2293 # find spiral surfaces of the stack:
2294 stackDimTags = [(3, windingTag) for windingTag in stack]
2295 stackSurfacesDimTags = gmsh.model.getBoundary(
2296 stackDimTags, combined=True, oriented=False
2297 )
2298 stackCurvesDimTags = gmsh.model.getBoundary(
2299 stackSurfacesDimTags, combined=False, oriented=False
2300 )
2301 # find the curve furthest from the origin:
2302 curveObjects = []
2303 for curveDimTag in stackCurvesDimTags:
2304 curveObject = curve(curveDimTag[1], self.geo)
2305 if curveObject.type is curveType.spiralArc:
2306 curveObjectDistanceFromOrigin = math.sqrt(
2307 curveObject.points[0][0] ** 2
2308 + curveObject.points[0][1] ** 2
2309 )
2310 curveObjects.append(
2311 (curveObject, curveObjectDistanceFromOrigin)
2312 )
2314 # sort the curves based on their distance from the origin (furthest first)
2315 curveObjects.sort(key=lambda x: x[1], reverse=True)
2317 curveTags = [curveObject[0].tag for curveObject in curveObjects]
2319 # only keep half of the curveTags:
2320 furthestCurveTags = curveTags[: len(curveTags) // 2]
2322 stackSpiralSurfaces = []
2323 for surfaceDimTag in stackSurfacesDimTags:
2324 normal = gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5])
2325 if abs(normal[2]) < 1e-5:
2326 curves = gmsh.model.getBoundary(
2327 [surfaceDimTag], combined=False, oriented=False
2328 )
2329 curveTags = [dimTag[1] for dimTag in curves]
2330 for curveTag in curveTags:
2331 if curveTag in furthestCurveTags:
2332 stackSpiralSurfaces.append(surfaceDimTag[1])
2333 break
2335 # add inner terminal surfaces too:
2336 if len(windingStacks) >= NofStacks:
2337 correspondingWindingStack = windingStacks[
2338 len(windingStacks) - NofStacks
2339 ]
2340 correspondingWindings = correspondingWindingStack[0]
2341 correspondingSurfaces = gmsh.model.getBoundary(
2342 [(3, windingTag) for windingTag in correspondingWindings],
2343 combined=True,
2344 oriented=False,
2345 )
2346 correspondingSurfaceTags = [
2347 dimTag[1] for dimTag in correspondingSurfaces
2348 ]
2349 for surface in allInnerTerminalContactLayerSurfaces:
2350 if surface in correspondingSurfaceTags:
2351 stackSpiralSurfaces.append(surface)
2353 windingStacks.append((stack, stackSpiralSurfaces))
2355 windingSets = []
2356 contactLayerSets = []
2357 for j in range(NofSets):
2358 windingTags = [
2359 windingTags for windingTags, _ in windingStacks[j::NofSets]
2360 ]
2361 windingTags = list(itertools.chain.from_iterable(windingTags))
2363 surfaceTags = [
2364 surfaceTags for _, surfaceTags in windingStacks[j::NofSets]
2365 ]
2366 surfaceTags = list(itertools.chain.from_iterable(surfaceTags))
2368 windingSets.append(windingTags)
2369 contactLayerSets.append(surfaceTags)
2371 # windingSets is a list with a length of NofSets.
2372 # finalWindingSets is also a list with a length of NofSets.
2373 for j, (windingSet, contactLayerSet) in enumerate(
2374 zip(windingSets, contactLayerSets)
2375 ):
2376 finalWindingSets[j].extend(windingSet)
2377 finalContactLayerSets[j].extend(contactLayerSet)
2379 # Seperate transition layer:
2380 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2381 self.dimTags[self.geo.terminals.inner.name]
2382 + self.dimTags[self.geo.terminals.outer.name]
2383 + self.dimTags["innerTransitionNotch"]
2384 + self.dimTags["outerTransitionNotch"],
2385 combined=False,
2386 oriented=False,
2387 )
2388 terminalAndNotchSurfaceTags = set(
2389 [dimTag[1] for dimTag in terminalAndNotchSurfaces]
2390 )
2392 contactLayerSets = []
2393 terminalWindingContactLayerSets = []
2394 for j in range(NofSets):
2395 contactLayerSets.append([])
2396 terminalWindingContactLayerSets.append([])
2398 for j in range(NofSets):
2399 allContactLayersInTheSet = finalContactLayerSets[j]
2401 insulatorList = []
2402 windingTerminalInsulatorList = []
2403 for contactLayer in allContactLayersInTheSet:
2404 if contactLayer in terminalAndNotchSurfaceTags:
2405 windingTerminalInsulatorList.append(contactLayer)
2406 else:
2407 insulatorList.append(contactLayer)
2409 contactLayerSets[j].extend(set(insulatorList))
2410 terminalWindingContactLayerSets[j].extend(set(windingTerminalInsulatorList))
2412 allContactLayerSurfacesForAllPancakes = []
2413 for j in range(NofSets):
2414 # Add winding volumes:
2415 self.regions.addEntities(
2416 self.geo.winding.name + "-" + str(j + 1),
2417 finalWindingSets[j],
2418 regionType.powered,
2419 entityType.vol,
2420 )
2422 # Add insulator surfaces:
2423 self.regions.addEntities(
2424 self.geo.contactLayer.name + "-" + str(j + 1),
2425 contactLayerSets[j],
2426 regionType.insulator,
2427 entityType.surf,
2428 )
2429 allContactLayerSurfacesForAllPancakes.extend(contactLayerSets[j])
2431 # Add terminal and winding contact layer:
2432 allContactLayerSurfacesForAllPancakes.extend(
2433 terminalWindingContactLayerSets[j]
2434 )
2436 # Add intersection of transition notch boundary and WindingAndTerminalContactLayer:
2437 transitionNotchSurfaces = gmsh.model.getBoundary(
2438 self.dimTags["innerTransitionNotch"]
2439 + self.dimTags["outerTransitionNotch"],
2440 combined=False,
2441 oriented=False,
2442 recursive=False
2443 )
2445 terminalContactLayerMinusNotch = set(terminalWindingContactLayerSets[j]).difference([tag for (dim, tag) in transitionNotchSurfaces])
2447 self.regions.addEntities(
2448 "WindingAndTerminalContactLayerWithoutNotch" + "-" + str(j + 1),
2449 list(terminalContactLayerMinusNotch),
2450 regionType.insulator,
2451 entityType.surf,
2452 )
2454 notchAndTerminalContactLayerIntersection = set([tag for (dim, tag) in transitionNotchSurfaces]).intersection(terminalWindingContactLayerSets[j])
2456 self.regions.addEntities(
2457 "WindingAndTerminalContactLayerOnlyNotch" + "-" + str(j + 1),
2458 list(notchAndTerminalContactLayerIntersection),
2459 regionType.insulator,
2460 entityType.surf,
2461 )
2463 allContactLayerSurfacesForAllPancakes = list(
2464 set(allContactLayerSurfacesForAllPancakes)
2465 )
2466 # Get insulator's boundary line that touches the air (required for the
2467 # pro file formulation):
2468 allContactLayerSurfacesForAllPancakesDimTags = [
2469 (2, surfaceTag) for surfaceTag in allContactLayerSurfacesForAllPancakes
2470 ]
2471 insulatorBoundary = gmsh.model.getBoundary(
2472 allContactLayerSurfacesForAllPancakesDimTags,
2473 combined=True,
2474 oriented=False,
2475 )
2476 insulatorBoundaryTags = [dimTag[1] for dimTag in insulatorBoundary]
2478 # Add insulator boundary lines:
2479 # Vertical lines should be removed from the insulator boundary because
2480 # they touch the terminals, not the air:
2481 verticalInsulatorBoundaryTags = []
2482 insulatorBoundaryTagsCopy = insulatorBoundaryTags.copy()
2483 for lineTag in insulatorBoundaryTagsCopy:
2484 lineObject = curve(lineTag, self.geo)
2485 if lineObject.type is curveType.axial:
2486 verticalInsulatorBoundaryTags.append(lineTag)
2487 insulatorBoundaryTags.remove(lineTag)
2489 # Create regions:
2490 self.regions.addEntities(
2491 self.geo.contactLayerBoundaryName,
2492 insulatorBoundaryTags,
2493 regionType.insulator,
2494 entityType.curve,
2495 )
2496 self.regions.addEntities(
2497 self.geo.contactLayerBoundaryName + "-TouchingTerminal",
2498 verticalInsulatorBoundaryTags,
2499 regionType.insulator,
2500 entityType.curve,
2501 )
2503 innerTransitionNotchTags = [
2504 dimTag[1] for dimTag in self.dimTags["innerTransitionNotch"]
2505 ]
2506 outerTransitionNotchTags = [
2507 dimTag[1] for dimTag in self.dimTags["outerTransitionNotch"]
2508 ]
2509 self.regions.addEntities(
2510 "innerTransitionNotch",
2511 innerTransitionNotchTags,
2512 regionType.powered,
2513 entityType.vol,
2514 )
2515 self.regions.addEntities(
2516 "outerTransitionNotch",
2517 outerTransitionNotchTags,
2518 regionType.powered,
2519 entityType.vol,
2520 )
2521 # ==============================================================================
2522 # WINDING AND CONTACT LAYER REGIONS ENDS =======================================
2523 # ==============================================================================
2525 # ==============================================================================
2526 # TERMINAL REGIONS START =======================================================
2527 # ==============================================================================
2529 innerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.terminals.inner.name]]
2530 self.regions.addEntities(
2531 self.geo.terminals.inner.name, innerTerminalTags, regionType.powered, entityType.vol_in
2532 )
2533 outerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.terminals.outer.name]]
2534 self.regions.addEntities(
2535 self.geo.terminals.outer.name,
2536 outerTerminalTags,
2537 regionType.powered,
2538 entityType.vol_out,
2539 )
2541 # Top and bottom terminal surfaces:
2542 firstTerminalDimTags = self.dimTags[self.geo.terminals.firstName]
2543 lastTerminalDimTags = self.dimTags[self.geo.terminals.lastName]
2545 if self.mesh.terminals.structured:
2546 topSurfaceDimTags = []
2547 for i in [1, 2, 3, 4]:
2548 lastTerminalSurfaces = gmsh.model.getBoundary(
2549 [lastTerminalDimTags[-i]], combined=False, oriented=False
2550 )
2551 topSurfaceDimTags.append(lastTerminalSurfaces[-1])
2552 else:
2553 lastTerminalSurfaces = gmsh.model.getBoundary(
2554 [lastTerminalDimTags[-1]], combined=False, oriented=False
2555 )
2556 topSurfaceDimTags = [lastTerminalSurfaces[-1]]
2557 topSurfaceTags = [dimTag[1] for dimTag in topSurfaceDimTags]
2559 if self.mesh.terminals.structured:
2560 bottomSurfaceDimTags = []
2561 for i in [1, 2, 3, 4]:
2562 firstTerminalSurfaces = gmsh.model.getBoundary(
2563 [firstTerminalDimTags[-i]], combined=False, oriented=False
2564 )
2565 bottomSurfaceDimTags.append(firstTerminalSurfaces[-1])
2566 else:
2567 firstTerminalSurfaces = gmsh.model.getBoundary(
2568 [firstTerminalDimTags[-1]], combined=False, oriented=False
2569 )
2570 bottomSurfaceDimTags = [firstTerminalSurfaces[-1]]
2571 bottomSurfaceTags = [dimTag[1] for dimTag in bottomSurfaceDimTags]
2573 self.regions.addEntities(
2574 "TopSurface",
2575 topSurfaceTags,
2576 regionType.powered,
2577 entityType.surf_out,
2578 )
2579 self.regions.addEntities(
2580 "BottomSurface",
2581 bottomSurfaceTags,
2582 regionType.powered,
2583 entityType.surf_in,
2584 )
2586 # if self.geo.contactLayer.tsa:
2587 # outerTerminalSurfaces = gmsh.model.getBoundary(
2588 # self.dimTags[self.geo.terminals.o.name], combined=True, oriented=False
2589 # )
2590 # outerTerminalSurfaces = [dimTag[1] for dimTag in outerTerminalSurfaces]
2591 # innerTerminalSurfaces = gmsh.model.getBoundary(
2592 # self.dimTags[self.geo.terminals.i.name], combined=True, oriented=False
2593 # )
2594 # innerTerminalSurfaces = [dimTag[1] for dimTag in innerTerminalSurfaces]
2595 # windingSurfaces = gmsh.model.getBoundary(
2596 # self.dimTags[self.geo.winding.name] + self.dimTags[self.geo.contactLayer.name],
2597 # combined=True,
2598 # oriented=False,
2599 # )
2600 # windingSurfaces = [dimTag[1] for dimTag in windingSurfaces]
2602 # windingAndOuterTerminalCommonSurfaces = list(
2603 # set(windingSurfaces).intersection(set(outerTerminalSurfaces))
2604 # )
2605 # windingAndInnerTerminalCommonSurfaces = list(
2606 # set(windingSurfaces).intersection(set(innerTerminalSurfaces))
2607 # )
2609 # self.regions.addEntities(
2610 # "WindingAndTerminalContactLayer",
2611 # windingAndOuterTerminalCommonSurfaces
2612 # + windingAndInnerTerminalCommonSurfaces,
2613 # regionType.insulator,
2614 # entityType.surf,
2615 # )
2617 # ==============================================================================
2618 # TERMINAL REGIONS ENDS ========================================================
2619 # ==============================================================================
2621 # ==============================================================================
2622 # AIR AND AIR SHELL REGIONS STARTS =============================================
2623 # ==============================================================================
2624 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.name]]
2625 self.regions.addEntities(
2626 self.geo.air.name, airTags, regionType.air, entityType.vol
2627 )
2629 # Create a region with two points on air to be used in the pro file formulation:
2630 # To those points, Phi=0 boundary condition will be applied to set the gauge.
2631 outerAirSurfaces = gmsh.model.getBoundary(
2632 self.dimTags[self.geo.air.name + "-OuterTube"], combined=True, oriented=False
2633 )
2634 outerAirSurface = outerAirSurfaces[-1]
2635 outerAirCurves = gmsh.model.getBoundary(
2636 [outerAirSurface], combined=True, oriented=False
2637 )
2638 outerAirCurve = outerAirCurves[-1]
2639 outerAirPoint = gmsh.model.getBoundary(
2640 [outerAirCurve], combined=False, oriented=False
2641 )
2642 outerAirPointTag = outerAirPoint[0][1]
2643 self.regions.addEntities(
2644 "OuterAirPoint",
2645 [outerAirPointTag],
2646 regionType.air,
2647 entityType.point,
2648 )
2650 innerAirSurfaces = gmsh.model.getBoundary(
2651 self.dimTags[self.geo.air.name + "-InnerCylinder"],
2652 combined=True,
2653 oriented=False,
2654 )
2655 innerAirSurface = innerAirSurfaces[0]
2656 innerAirCurves = gmsh.model.getBoundary(
2657 [innerAirSurface], combined=True, oriented=False
2658 )
2659 innerAirCurve = innerAirCurves[-1]
2660 innerAirPoint = gmsh.model.getBoundary(
2661 [innerAirCurve], combined=False, oriented=False
2662 )
2663 innerAirPointTag = innerAirPoint[0][1]
2664 self.regions.addEntities(
2665 "InnerAirPoint",
2666 [innerAirPointTag],
2667 regionType.air,
2668 entityType.point,
2669 )
2671 if self.geo.air.shellTransformation:
2672 if self.geo.air.type == "cylinder":
2673 airShellTags = [
2674 dimTag[1] for dimTag in self.dimTags[self.geo.air.shellVolumeName]
2675 ]
2676 self.regions.addEntities(
2677 self.geo.air.shellVolumeName,
2678 airShellTags,
2679 regionType.air_far_field,
2680 entityType.vol,
2681 )
2682 elif self.geo.air.type == "cuboid":
2683 airShell1Tags = [
2684 dimTag[1]
2685 for dimTag in self.dimTags[self.geo.air.shellVolumeName + "-Part1"]
2686 + self.dimTags[self.geo.air.shellVolumeName + "-Part3"]
2687 ]
2688 airShell2Tags = [
2689 dimTag[1]
2690 for dimTag in self.dimTags[self.geo.air.shellVolumeName + "-Part2"]
2691 + self.dimTags[self.geo.air.shellVolumeName + "-Part4"]
2692 ]
2693 self.regions.addEntities(
2694 self.geo.air.shellVolumeName + "-PartX",
2695 airShell1Tags,
2696 regionType.air_far_field,
2697 entityType.vol,
2698 )
2699 self.regions.addEntities(
2700 self.geo.air.shellVolumeName + "-PartY",
2701 airShell2Tags,
2702 regionType.air_far_field,
2703 entityType.vol,
2704 )
2705 # ==============================================================================
2706 # AIR AND AIR SHELL REGIONS ENDS ===============================================
2707 # ==============================================================================
2709 # ==============================================================================
2710 # CUTS STARTS ==================================================================
2711 # ==============================================================================
2712 if self.geo.air.cutName in self.dimTags:
2713 cutTags = [dimTag[1] for dimTag in self.dimTags[self.geo.air.cutName]]
2714 self.regions.addEntities(
2715 self.geo.air.cutName, cutTags, regionType.air, entityType.cochain
2716 )
2718 if "CutsForPerfectInsulation" in self.dimTags:
2719 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsForPerfectInsulation"]]
2720 self.regions.addEntities(
2721 "CutsForPerfectInsulation", cutTags, regionType.air, entityType.cochain
2722 )
2724 if "CutsBetweenPancakes" in self.dimTags:
2725 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsBetweenPancakes"]]
2726 self.regions.addEntities(
2727 "CutsBetweenPancakes", cutTags, regionType.air, entityType.cochain
2728 )
2729 # ==============================================================================
2730 # CUTS ENDS ====================================================================
2731 # ==============================================================================
2733 # ==============================================================================
2734 # PANCAKE BOUNDARY SURFACE STARTS ==============================================
2735 # ==============================================================================
2736 # Pancake3D Boundary Surface:
2737 allPancakeVolumes = (
2738 self.dimTags[self.geo.winding.name]
2739 + self.dimTags[self.geo.terminals.inner.name]
2740 + self.dimTags[self.geo.terminals.outer.name]
2741 + self.dimTags[self.geo.contactLayer.name]
2742 + self.dimTags["innerTransitionNotch"]
2743 + self.dimTags["outerTransitionNotch"]
2744 )
2745 Pancake3DAllBoundary = gmsh.model.getBoundary(
2746 allPancakeVolumes, combined=True, oriented=False
2747 )
2748 Pancake3DBoundaryDimTags = list(
2749 set(Pancake3DAllBoundary)
2750 - set(topSurfaceDimTags)
2751 - set(bottomSurfaceDimTags)
2752 )
2753 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2754 self.regions.addEntities(
2755 self.geo.pancakeBoundaryName,
2756 pancake3DBoundaryTags,
2757 regionType.powered,
2758 entityType.surf,
2759 )
2761 if self.geo.contactLayer.thinShellApproximation:
2762 # add non-winding boundary for convective cooling
2763 windingBoundaryDimTags = gmsh.model.getBoundary(
2764 [(3, tag) for tag in itertools.chain(*finalWindingSets)],
2765 combined=True,
2766 oriented=False,
2767 )
2769 inner_terminal_and_transition_notch_all_boundaries = gmsh.model.getBoundary(
2770 self.dimTags[self.geo.terminals.inner.name] + self.dimTags["innerTransitionNotch"],
2771 combined=True,
2772 oriented=False
2773 )
2775 inner_terminal_and_transition_notch_boundary_dim_tags = set(Pancake3DBoundaryDimTags).intersection(inner_terminal_and_transition_notch_all_boundaries)
2776 inner_terminal_and_transition_notch_boundary_tags = [dimTag[1] for dimTag in inner_terminal_and_transition_notch_boundary_dim_tags]
2777 self.regions.addEntities(
2778 f"{self.geo.pancakeBoundaryName}-InnerTerminalAndTransitionNotch",
2779 inner_terminal_and_transition_notch_boundary_tags,
2780 regionType.powered,
2781 entityType.surf_th,
2782 )
2784 outer_terminal_and_transition_notch_all_boundaries = gmsh.model.getBoundary(
2785 self.dimTags[self.geo.terminals.outer.name] + self.dimTags["outerTransitionNotch"],
2786 combined=True,
2787 oriented=False
2788 )
2790 outer_terminal_and_transition_notch_boundary_dim_tags = set(Pancake3DBoundaryDimTags).intersection(outer_terminal_and_transition_notch_all_boundaries)
2791 outer_terminal_and_transition_notch_boundary_tags = [dimTag[1] for dimTag in outer_terminal_and_transition_notch_boundary_dim_tags]
2792 self.regions.addEntities(
2793 f"{self.geo.pancakeBoundaryName}-outerTerminalAndTransitionNotch",
2794 outer_terminal_and_transition_notch_boundary_tags,
2795 regionType.powered,
2796 entityType.surf_th,
2797 )
2799 # add pancake boundary for convective cooling following the winding numbering logic
2800 # computes the intersection between pancake boundary and the boundary of each winding group
2801 for j in range(NofSets):
2803 windingBoundaryDimTags = gmsh.model.getBoundary(
2804 [(3, tag) for tag in finalWindingSets[j]],
2805 combined=True,
2806 oriented=False,
2807 )
2809 windingBoundaryDimTags = set(windingBoundaryDimTags).intersection(Pancake3DBoundaryDimTags)
2811 windingBoundaryTags = [dimTag[1] for dimTag in windingBoundaryDimTags]
2812 self.regions.addEntities(
2813 f"{self.geo.pancakeBoundaryName}-Winding{j+1}",
2814 windingBoundaryTags,
2815 regionType.powered,
2816 entityType.surf_th,
2817 )
2819 if not self.geo.contactLayer.thinShellApproximation:
2820 # Pancake3D Boundary Surface with only winding and terminals:
2821 allPancakeVolumes = (
2822 self.dimTags[self.geo.winding.name]
2823 + self.dimTags[self.geo.terminals.inner.name]
2824 + self.dimTags[self.geo.terminals.outer.name]
2825 + self.dimTags["innerTransitionNotch"]
2826 + self.dimTags["outerTransitionNotch"]
2827 + [(3, tag) for tag in contactLayerBetweenTerminalsAndWinding]
2828 )
2829 Pancake3DAllBoundary = gmsh.model.getBoundary(
2830 allPancakeVolumes, combined=True, oriented=False
2831 )
2832 Pancake3DBoundaryDimTags = list(
2833 set(Pancake3DAllBoundary)
2834 - set(topSurfaceDimTags)
2835 - set(bottomSurfaceDimTags)
2836 )
2837 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2838 self.regions.addEntities(
2839 self.geo.pancakeBoundaryName + "-OnlyWindingAndTerminals",
2840 pancake3DBoundaryTags,
2841 regionType.powered,
2842 entityType.surf,
2843 )
2845 # ==============================================================================
2846 # PANCAKE BOUNDARY SURFACE ENDS ================================================
2847 # ==============================================================================
2849 # Generate regions file:
2850 self.regions.generateRegionsFile(self.regions_file)
2851 self.rm = FilesAndFolders.read_data_from_yaml(self.regions_file, RegionsModel)
2853 logger.info(
2854 "Generating physical groups and regions file has been finished in"
2855 f" {timeit.default_timer() - start_time:.2f} s."
2856 )
2858 def generate_mesh_file(self):
2859 """
2860 Saves mesh file to disk.
2863 """
2864 logger.info(
2865 f"Generating Pancake3D mesh file ({self.mesh_file}) has been started."
2866 )
2867 start_time = timeit.default_timer()
2869 gmsh.write(self.mesh_file)
2871 logger.info(
2872 f"Generating Pancake3D mesh file ({self.mesh_file}) has been finished"
2873 f" in {timeit.default_timer() - start_time:.2f} s."
2874 )
2876 if self.mesh_gui:
2877 gmsh.option.setNumber("Geometry.Volumes", 0)
2878 gmsh.option.setNumber("Geometry.Surfaces", 0)
2879 gmsh.option.setNumber("Geometry.Curves", 0)
2880 gmsh.option.setNumber("Geometry.Points", 0)
2881 self.gu.launch_interactive_GUI()
2882 else:
2883 gmsh.clear()
2884 gmsh.finalize()
2886 def load_mesh(self):
2887 """
2888 Loads mesh from .msh file.
2891 """
2892 logger.info("Loading Pancake3D mesh has been started.")
2893 start_time = timeit.default_timer()
2895 previousGeo = FilesAndFolders.read_data_from_yaml(
2896 self.geometry_data_file, Pancake3DGeometry
2897 )
2898 previousMesh = FilesAndFolders.read_data_from_yaml(
2899 self.mesh_data_file, Pancake3DMesh
2900 )
2902 if previousGeo.dict() != self.geo.dict():
2903 raise ValueError(
2904 "Geometry data has been changed. Please regenerate the geometry or load"
2905 " the previous geometry data."
2906 )
2907 elif previousMesh.dict() != self.mesh.dict():
2908 raise ValueError(
2909 "Mesh data has been changed. Please regenerate the mesh or load the"
2910 " previous mesh data."
2911 )
2913 gmsh.clear()
2914 gmsh.open(self.mesh_file)
2916 logger.info(
2917 "Loading Pancake3D mesh has been finished in"
2918 f" {timeit.default_timer() - start_time:.2f} s."
2919 )