Coverage for fiqus/mesh_generators/MeshPancake3D.py: 88%
999 statements
« prev ^ index » next coverage.py v6.4.4, created at 2024-05-20 03:24 +0200
1import timeit
2import json
3import logging
4import math
5from enum import Enum
6import operator
7import itertools
8import re
10import gmsh
11import scipy.integrate
12import numpy as np
14from fiqus.data import RegionsModelFiQuS
15from fiqus.utils.Utils import GmshUtils
16from fiqus.data.RegionsModelFiQuS import RegionsModel
17from fiqus.mains.MainPancake3D import Base
18from fiqus.utils.Utils import FilesAndFolders
19from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry, Pancake3DMesh
22logger = logging.getLogger(__name__)
25class regionType(str, Enum):
26 """
27 A class to specify region type easily in the regions class.
28 """
30 powered = "powered"
31 insulator = "insulator"
32 air = "air"
33 air_far_field = "air_far_field"
36class entityType(str, Enum):
37 """
38 A class to specify entity type easily in the regions class.
39 """
41 vol = "vol"
42 vol_in = "vol_in"
43 vol_out = "vol_out"
44 surf = "surf"
45 surf_th = "surf_th"
46 surf_in = "surf_in"
47 surf_out = "surf_out"
48 surf_ext = "surf_ext"
49 cochain = "cochain"
50 curve = "curve"
51 point = "point"
54class regions:
55 """
56 A class to generate physical groups in GMSH and create the corresponding regions
57 file. The regions file is the file where the region tags are stored in the FiQuS
58 regions data model convention. The file is used to template the *.pro file (GetDP
59 input file).
60 """
62 def __init__(self):
63 # Main types of entities:
64 # The keys are the FiQuS region categories, and the values are the corresponding
65 # GMSH entity type.
66 self.entityMainType = {
67 "vol": "vol",
68 "vol_in": "vol",
69 "vol_out": "vol",
70 "surf": "surf",
71 "surf_th": "surf",
72 "surf_in": "surf",
73 "surf_out": "surf",
74 "surf_ext": "surf",
75 "cochain": "curve",
76 "curve": "curve",
77 "point": "point",
78 }
80 # Dimensions of entity types:
81 self.entityDim = {"vol": 3, "surf": 2, "curve": 1, "point": 0}
83 # Keys for regions file. The keys are appended to the name of the regions
84 # accordingly.
85 self.entityKey = {
86 "vol": "",
87 "vol_in": "",
88 "vol_out": "",
89 "surf": "_bd",
90 "surf_th": "_bd",
91 "surf_in": "_in",
92 "surf_out": "_out",
93 "surf_ext": "_ext",
94 "cochain": "_cut",
95 "curve": "_curve",
96 "point": "_point",
97 }
99 # FiQuS convetion for region numbers:
100 self.regionTags = {
101 "vol": 1000000, # volume region tag start
102 "surf": 2000000, # surface region tag start
103 "curve": 3000000, # curve region tag start
104 "point": 4000000, # point region tag start
105 }
107 # Initialize the regions model:
108 self.rm = RegionsModelFiQuS.RegionsModel()
110 # This is used because self.rm.powered is not initialized in
111 # RegionsModelFiQuS.RegionsModel. It should be fixed in the future.
112 self.rm.powered = 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.vol.names = []
117 self.rm.powered.vol.numbers = []
118 self.rm.powered.surf.names = []
119 self.rm.powered.surf.numbers = []
120 self.rm.powered.surf_th.names = []
121 self.rm.powered.surf_th.numbers = []
122 self.rm.powered.surf_in.names = []
123 self.rm.powered.surf_in.numbers = []
124 self.rm.powered.surf_out.names = []
125 self.rm.powered.surf_out.numbers = []
126 self.rm.powered.curve.names = []
127 self.rm.powered.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, entityType).name = name
171 getattr(self.rm.powered, entityType).number = regionTag
173 else:
174 getattr(self.rm.powered, entityType).names.append(name)
175 getattr(self.rm.powered, 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.dimTol 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.dimTol) * self.geo.dimTol
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.wi.theta_i) / 2 / math.pi
301 self.n1 = round(self.n1 / self.geo.wi.turnTol) * self.geo.wi.turnTol
303 self.n2 = (theta2 - self.geo.wi.theta_i) / 2 / math.pi
304 self.n2 = round(self.n2 / self.geo.wi.turnTol) * self.geo.wi.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.dimTol
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.dimTol) and math.isclose(
333 cm[1], ycm, abs_tol=self.geo.dimTol
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()
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.N):
484 # Get the volume tags:
485 windingVolumeDimTags = self.dimTags[self.geo.wi.name + str(i + 1)]
486 windingVolumeTags = [dimTag[1] for dimTag in windingVolumeDimTags]
488 contactLayerVolumeDimTags = self.dimTags[self.geo.ii.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.wi.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.ai.name + "-TopPancakeWindingExtursion"
552 ]
554 airTopContactLayerExtrusionVolumeDimTags = self.dimTags[
555 self.geo.ai.name + "-TopPancakeContactLayerExtursion"
556 ]
558 airTopTerminalsExtrusionVolumeDimTags = self.dimTags[
559 self.geo.ai.name + "-TopTerminalsExtrusion"
560 ]
562 airBottomWindingExtrusionVolumeDimTags = self.dimTags[
563 self.geo.ai.name + "-BottomPancakeWindingExtursion"
564 ]
566 airBottomContactLayerExtrusionVolumeDimTags = self.dimTags[
567 self.geo.ai.name + "-BottomPancakeContactLayerExtursion"
568 ]
570 airBottomTerminalsExtrusionVolumeDimTags = self.dimTags[
571 self.geo.ai.name + "-BottomTerminalsExtrusion"
572 ]
574 removedAirVolumeDimTags = []
575 newAirVolumeDimTags = []
576 if self.mesh.ai.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.wi.axne) / self.geo.wi.h
593 axneForAir = round(
594 axialElementsPerLengthForWinding * self.geo.ai.margin + 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.N - 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.ai.name + "-OuterTube"]
651 airOuterTubeVolumeTags = [dimTag[1] for dimTag in airOuterTubeVolumeDimTags]
653 airTopTubeTerminalsVolumeDimTags = self.dimTags[
654 self.geo.ai.name + "-TopTubeTerminalsExtrusion"
655 ]
656 airTopTubeTerminalsVolumeTags = [
657 dimTag[1] for dimTag in airTopTubeTerminalsVolumeDimTags
658 ]
660 airBottomTubeTerminalsVolumeDimTags = self.dimTags[
661 self.geo.ai.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.ai.name + "-InnerCylinder"
670 ]
671 airInnerCylinderVolumeTags = [
672 dimTag[1] for dimTag in airInnerCylinderVolumeDimTags
673 ]
675 airTubesAndCylinders = airOuterTubeVolumeTags + airInnerCylinderVolumeTags
677 if self.geo.ai.shellTransformation:
678 shellVolumes = self.dimTags[self.geo.ai.shellVolumeName]
679 shellVolumeTags = [dimTag[1] for dimTag in shellVolumes]
680 airTubesAndCylinders.extend(shellVolumeTags)
682 airRadialElementMultiplier = 1 / self.mesh.ai.radialElementSize
683 self.structure_tubes_and_cylinders(
684 airTubesAndCylinders,
685 radialElementMultiplier=airRadialElementMultiplier,
686 )
688 if self.mesh.ti.structured:
689 terminalsRadialElementMultiplier = 1 / self.mesh.ti.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.ai.name + "-InnerCylinder"
751 ]
752 if self.geo.N > 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.ai.name + "-InnerCylinder"] = [
775 airInnerCylinderVolumeDimTag,
776 airInnerCylinderVolumeDimTagFirst,
777 airInnerCylinderVolumeDimTagLast,
778 ]
779 else:
780 pass
781 # self.dimTags[self.geo.ai.name + "-InnerCylinder"] = [
782 # self.dimTags[self.geo.ai.name + "-InnerCylinder"][1],
783 # self.dimTags[self.geo.ai.name + "-InnerCylinder"][0],
784 # self.dimTags[self.geo.ai.name + "-InnerCylinder"][2],
785 # ]
787 airOuterTubeVolumeDimTags = self.dimTags[self.geo.ai.name + "-OuterTube"]
788 airOuterTubeVolumeDimTag = Mesh.fuse_volumes(
789 airOuterTubeVolumeDimTags,
790 fuseSurfaces=True,
791 fusedSurfacesArePlane=False,
792 )
793 newAirOuterTubeVolumeDimTag = airOuterTubeVolumeDimTag
794 removedAirVolumeDimTags.extend(airOuterTubeVolumeDimTags)
796 if self.geo.ai.shellTransformation:
797 # Fuse air shell volumes:
798 if self.geo.ai.type == "cylinder":
799 removedShellVolumeDimTags = []
800 shellVolumeDimTags = self.dimTags[self.geo.ai.shellVolumeName]
801 shellVolumeDimTag = Mesh.fuse_volumes(
802 shellVolumeDimTags,
803 fuseSurfaces=True,
804 fusedSurfacesArePlane=False,
805 )
806 removedShellVolumeDimTags.extend(shellVolumeDimTags)
807 newShellVolumeDimTags = [shellVolumeDimTag]
808 for removedDimTag in removedShellVolumeDimTags:
809 self.dimTags[self.geo.ai.shellVolumeName].remove(removedDimTag)
810 elif self.geo.ai.type == "cuboid":
811 # Unfortunately, surfaces cannot be combined for the cuboid type of air.
812 # However, it doesn't affect the mesh quality that much.
813 newShellVolumeDimTags = []
815 shellPart1VolumeDimTag = Mesh.fuse_volumes(
816 self.dimTags[self.geo.ai.shellVolumeName + "-Part1"],
817 fuseSurfaces=False,
818 )
819 self.dimTags[self.geo.ai.shellVolumeName + "-Part1"] = [
820 shellPart1VolumeDimTag
821 ]
823 shellPart2VolumeDimTag = Mesh.fuse_volumes(
824 self.dimTags[self.geo.ai.shellVolumeName + "-Part2"],
825 fuseSurfaces=False,
826 )
827 self.dimTags[self.geo.ai.shellVolumeName + "-Part2"] = [
828 shellPart2VolumeDimTag
829 ]
831 shellPart3VolumeDimTag = Mesh.fuse_volumes(
832 self.dimTags[self.geo.ai.shellVolumeName + "-Part3"],
833 fuseSurfaces=False,
834 )
835 self.dimTags[self.geo.ai.shellVolumeName + "-Part3"] = [
836 shellPart3VolumeDimTag
837 ]
839 shellPart4VolumeDimTag = Mesh.fuse_volumes(
840 self.dimTags[self.geo.ai.shellVolumeName + "-Part4"],
841 fuseSurfaces=False,
842 )
843 self.dimTags[self.geo.ai.shellVolumeName + "-Part4"] = [
844 shellPart4VolumeDimTag
845 ]
847 # The problem is, shell volume and outer air tube volume has a common
848 # surface and that surface should be combined as well for high quality mesh.
849 # However, it can be only done for cylinder type of air for now.
850 # Get the combined boundary surfaces:
851 if self.geo.ai.type == "cylinder":
852 (
853 newAirOuterTubeVolumeDimTag,
854 newShellVolumeDimTag,
855 ) = Mesh.fuse_common_surfaces_of_two_volumes(
856 [airOuterTubeVolumeDimTag],
857 newShellVolumeDimTags,
858 fuseOtherSurfaces=False,
859 surfacesArePlane=False,
860 )
861 airOuterTubeVolumeDimTag = newAirOuterTubeVolumeDimTag
862 self.dimTags[self.geo.ai.shellVolumeName].append(
863 newShellVolumeDimTag
864 )
866 newAirVolumeDimTags.append(newAirOuterTubeVolumeDimTag)
868 # Update volume tags dictionary of air:
869 self.dimTags[self.geo.ai.name] = list(
870 (
871 set(self.dimTags[self.geo.ai.name]) - set(removedAirVolumeDimTags)
872 ).union(set(newAirVolumeDimTags))
873 )
875 # ==============================================================================
876 # MESHING AIR ENDS =============================================================
877 # ==============================================================================
879 # ==============================================================================
880 # MESHING TERMINALS STARTS =====================================================
881 # ==============================================================================
882 if self.mesh.ti.structured:
883 # Structure tubes of the terminals:
884 terminalOuterTubeVolumeDimTags = self.dimTags[self.geo.ti.o.name + "-Tube"]
885 terminalOuterTubeVolumeTags = [
886 dimTag[1] for dimTag in terminalOuterTubeVolumeDimTags
887 ]
888 terminalInnerTubeVolumeDimTags = self.dimTags[self.geo.ti.i.name + "-Tube"]
889 terminalInnerTubeVolumeTags = [
890 dimTag[1] for dimTag in terminalInnerTubeVolumeDimTags
891 ]
893 terminalsRadialElementMultiplier = 1 / self.mesh.ti.radialElementSize
894 self.structure_tubes_and_cylinders(
895 terminalOuterTubeVolumeTags + terminalInnerTubeVolumeTags,
896 radialElementMultiplier=terminalsRadialElementMultiplier,
897 )
899 # Structure nontube parts of the terminals:
900 terminalOuterNonTubeVolumeDimTags = self.dimTags[
901 self.geo.ti.o.name + "-Touching"
902 ]
903 terminalOuterNonTubeVolumeTags = [
904 dimTag[1] for dimTag in terminalOuterNonTubeVolumeDimTags
905 ]
906 terminalInnerNonTubeVolumeDimTags = self.dimTags[
907 self.geo.ti.i.name + "-Touching"
908 ]
909 terminalInnerNonTubeVolumeTags = [
910 dimTag[1] for dimTag in terminalInnerNonTubeVolumeDimTags
911 ]
913 self.structure_tubes_and_cylinders(
914 terminalInnerNonTubeVolumeTags + terminalOuterNonTubeVolumeTags,
915 terminalNonTubeParts=True,
916 radialElementMultiplier=terminalsRadialElementMultiplier,
917 )
918 # ==============================================================================
919 # MESHING TERMINALS ENDS =======================================================
920 # ==============================================================================
922 # ==============================================================================
923 # FIELD SETTINGS STARTS ========================================================
924 # ==============================================================================
926 # Mesh fields for the air:
927 # Meshes will grow as they get further from the field surfaces:
928 fieldSurfacesDimTags = gmsh.model.getBoundary(
929 self.dimTags[self.geo.wi.name], oriented=False, combined=True
930 )
931 fieldSurfacesTags = [dimTag[1] for dimTag in fieldSurfacesDimTags]
933 distanceField = gmsh.model.mesh.field.add("Distance")
935 gmsh.model.mesh.field.setNumbers(
936 distanceField,
937 "SurfacesList",
938 fieldSurfacesTags,
939 )
941 thresholdField = gmsh.model.mesh.field.add("Threshold")
942 gmsh.model.mesh.field.setNumber(thresholdField, "InField", distanceField)
943 gmsh.model.mesh.field.setNumber(thresholdField, "SizeMin", self.mesh.sizeMin)
944 gmsh.model.mesh.field.setNumber(thresholdField, "SizeMax", self.mesh.sizeMax)
945 gmsh.model.mesh.field.setNumber(
946 thresholdField, "DistMin", self.mesh.startGrowingDistance
947 )
949 gmsh.model.mesh.field.setNumber(
950 thresholdField, "DistMax", self.mesh.stopGrowingDistance
951 )
953 gmsh.model.mesh.field.setAsBackgroundMesh(thresholdField)
955 # ==============================================================================
956 # FIELD SETTINGS ENDS ==========================================================
957 # ==============================================================================
959 gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
960 gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
961 gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)
963 try:
964 # Only print warnings and errors:
965 # Don't print on terminal, because we will use logger:
966 gmsh.option.setNumber("General.Terminal", 0)
967 # Start logger:
968 gmsh.logger.start()
970 # Mesh:
971 gmsh.model.mesh.generate()
972 gmsh.model.mesh.optimize()
974 # Print the log:
975 log = gmsh.logger.get()
976 for line in log:
977 if line.startswith("Info"):
978 logger.info(re.sub(r"Info:\s+", "", line))
979 elif line.startswith("Warning"):
980 logger.warning(re.sub(r"Warning:\s+", "", line))
982 gmsh.logger.stop()
983 except:
984 # Print the log:
985 log = gmsh.logger.get()
986 for line in log:
987 if line.startswith("Info"):
988 logger.info(re.sub(r"Info:\s+", "", line))
989 elif line.startswith("Warning"):
990 logger.warning(re.sub(r"Warning:\s+", "", line))
991 elif line.startswith("Error"):
992 logger.error(re.sub(r"Error:\s+", "", line))
994 gmsh.logger.stop()
996 self.generate_regions()
998 logger.error(
999 "Meshing Pancake3D magnet has failed. Try to change"
1000 " minimumElementSize and maximumElementSize parameters."
1001 )
1002 raise
1004 # SICN not implemented in 1D!
1005 allElementsDim2 = list(gmsh.model.mesh.getElements(dim=2)[1][0])
1006 allElementsDim3 = list(gmsh.model.mesh.getElements(dim=3)[1][0])
1007 allElements = allElementsDim2 + allElementsDim3
1008 elementQualities = gmsh.model.mesh.getElementQualities(allElements)
1009 lowestQuality = min(elementQualities)
1010 averageQuality = sum(elementQualities) / len(elementQualities)
1011 NofLowQualityElements = len(
1012 [quality for quality in elementQualities if quality < 0.01]
1013 )
1014 NofIllElemets = len(
1015 [quality for quality in elementQualities if quality < 0.001]
1016 )
1018 logger.info(
1019 f"The lowest quality among the elements is {lowestQuality:.4f} (SICN). The"
1020 " number of elements with quality lower than 0.01 is"
1021 f" {NofLowQualityElements}."
1022 )
1024 if NofIllElemets > 0:
1025 logger.warning(
1026 f"There are {NofIllElemets} elements with quality lower than 0.001. Try"
1027 " to change minimumElementSize and maximumElementSize parameters."
1028 )
1030 # Create cuts:
1031 # This is required to make the air a simply connected domain. This is required
1032 # for the solution part. You can read more about Homology in GMSH documentation.
1033 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ai.name]]
1035 if self.geo.ai.shellTransformation:
1036 shellTags = [
1037 dimTag[1] for dimTag in self.dimTags[self.geo.ai.shellVolumeName]
1038 ]
1039 airTags.extend(shellTags)
1041 dummyAirRegion = gmsh.model.addPhysicalGroup(dim=3, tags=airTags)
1042 dummyAirRegionDimTag = (3, dummyAirRegion)
1044 innerCylinderTags = [self.dimTags[self.geo.ai.name + "-InnerCylinder"][0][1]]
1045 gapTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ai.name + "-Gap"]]
1046 # Only remove every second gap:
1047 gapTags = gapTags[1::2]
1049 dummyAirRegionWithoutInnerCylinder = gmsh.model.addPhysicalGroup(
1050 dim=3, tags=list(set(airTags) - set(innerCylinderTags) - set(gapTags))
1051 )
1052 dummyAirRegionWithoutInnerCylinderDimTag = (
1053 3,
1054 dummyAirRegionWithoutInnerCylinder,
1055 )
1057 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.wi.name]]
1058 dummyWindingRegion = gmsh.model.addPhysicalGroup(dim=3, tags=windingTags)
1059 dummyWindingRegionDimTag = (3, dummyWindingRegion)
1061 if self.geo.ii.tsa:
1062 # Find all the contact layer surfaces:
1063 allWindingDimTags = []
1064 for i in range(self.geo.N):
1065 windingDimTags = self.dimTags[self.geo.wi.name + str(i + 1)]
1066 allWindingDimTags.extend(windingDimTags)
1068 windingBoundarySurfaces = gmsh.model.getBoundary(
1069 allWindingDimTags, combined=True, oriented=False
1070 )
1071 allWindingSurfaces = gmsh.model.getBoundary(
1072 allWindingDimTags, combined=False, oriented=False
1073 )
1075 contactLayerSurfacesDimTags = list(
1076 set(allWindingSurfaces) - set(windingBoundarySurfaces)
1077 )
1078 contactLayerTags = [dimTag[1] for dimTag in contactLayerSurfacesDimTags]
1080 # Get rid of non-contactLayer surfaces:
1081 realContactLayerTags = []
1082 for contactLayerTag in contactLayerTags:
1083 surfaceNormal = list(gmsh.model.getNormal(contactLayerTag, [0.5, 0.5]))
1084 centerOfMass = gmsh.model.occ.getCenterOfMass(2, contactLayerTag)
1086 if (
1087 abs(
1088 surfaceNormal[0] * centerOfMass[0]
1089 + surfaceNormal[1] * centerOfMass[1]
1090 )
1091 > 1e-6
1092 ):
1093 realContactLayerTags.append(contactLayerTag)
1095 # Get rid of surfaces that touch terminals:
1096 terminalSurfaces = gmsh.model.getBoundary(
1097 self.dimTags[self.geo.ti.o.name] + self.dimTags[self.geo.ti.i.name],
1098 combined=False,
1099 oriented=False,
1100 )
1101 terminalSurfaces = [dimTag[1] for dimTag in terminalSurfaces]
1102 finalContactLayerTags = [
1103 tag for tag in realContactLayerTags if tag not in terminalSurfaces
1104 ]
1106 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1107 dim=2, tags=finalContactLayerTags
1108 )
1109 dummyContactLayerRegionDimTag = (2, dummyContactLayerRegion)
1111 else:
1112 contactLayerTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ii.name]]
1114 # get rid of volumes that touch terminals:
1115 terminalSurfaces = gmsh.model.getBoundary(
1116 self.dimTags[self.geo.ti.o.name] + self.dimTags[self.geo.ti.i.name],
1117 combined=False,
1118 oriented=False,
1119 )
1120 finalContactLayerTags = []
1121 for contactLayerTag in contactLayerTags:
1122 insulatorSurfaces = gmsh.model.getBoundary(
1123 [(3, contactLayerTag)], combined=False, oriented=False
1124 )
1125 itTouchesTerminals = False
1126 for insulatorSurface in insulatorSurfaces:
1127 if insulatorSurface in terminalSurfaces:
1128 itTouchesTerminals = True
1129 break
1131 if not itTouchesTerminals:
1132 finalContactLayerTags.append(contactLayerTag)
1134 dummyContactLayerRegion = gmsh.model.addPhysicalGroup(
1135 dim=3, tags=finalContactLayerTags
1136 )
1137 dummyContactLayerRegionDimTag = (3, dummyContactLayerRegion)
1139 # First cohomology request (normal cut for NI coils):
1140 gmsh.model.mesh.addHomologyRequest(
1141 "Cohomology",
1142 domainTags=[dummyAirRegion],
1143 subdomainTags=[],
1144 dims=[1],
1145 )
1147 # Second cohomology request (insulated cut for insulated coils):
1148 if self.geo.N > 1:
1149 gmsh.model.mesh.addHomologyRequest(
1150 "Cohomology",
1151 domainTags=[
1152 dummyAirRegionWithoutInnerCylinder,
1153 dummyContactLayerRegion,
1154 ],
1155 subdomainTags=[],
1156 dims=[1],
1157 )
1158 else:
1159 gmsh.model.mesh.addHomologyRequest(
1160 "Cohomology",
1161 domainTags=[
1162 dummyAirRegion,
1163 dummyContactLayerRegion,
1164 ],
1165 subdomainTags=[],
1166 dims=[1],
1167 )
1169 # Third cohomology request (for cuts between pancake coils):
1170 gmsh.model.mesh.addHomologyRequest(
1171 "Cohomology",
1172 domainTags=[
1173 dummyAirRegion,
1174 dummyContactLayerRegion,
1175 dummyWindingRegion,
1176 ],
1177 subdomainTags=[],
1178 dims=[1],
1179 )
1181 # Start logger:
1182 gmsh.logger.start()
1184 cuts = gmsh.model.mesh.computeHomology()
1186 # Print the log:
1187 log = gmsh.logger.get()
1188 for line in log:
1189 if line.startswith("Info"):
1190 logger.info(re.sub(r"Info:\s+", "", line))
1191 elif line.startswith("Warning"):
1192 logger.warning(re.sub(r"Warning:\s+", "", line))
1193 gmsh.logger.stop()
1195 if self.geo.N > 1:
1196 cutsDictionary = {
1197 "H^1{1}": self.geo.ai.cutName,
1198 "H^1{1,4,3}": "CutsBetweenPancakes",
1199 "H^1{2,4}": "CutsForPerfectInsulation",
1200 }
1201 else:
1202 cutsDictionary = {
1203 "H^1{1}": self.geo.ai.cutName,
1204 "H^1{1,4,3}": "CutsBetweenPancakes",
1205 "H^1{1,4}": "CutsForPerfectInsulation",
1206 }
1207 cutTags = [dimTag[1] for dimTag in cuts]
1208 cutEntities = []
1209 for tag in cutTags:
1210 name = gmsh.model.getPhysicalName(1, tag)
1211 cutEntities = list(gmsh.model.getEntitiesForPhysicalGroup(1, tag))
1212 cutEntitiesDimTags = [(1, cutEntity) for cutEntity in cutEntities]
1213 for key in cutsDictionary:
1214 if key in name:
1215 if cutsDictionary[key] in self.dimTags:
1216 self.dimTags[cutsDictionary[key]].extend(cutEntitiesDimTags)
1217 else:
1218 self.dimTags[cutsDictionary[key]] = cutEntitiesDimTags
1220 # Remove newly created physical groups because they will be created again in
1221 # generate_regions method.
1222 gmsh.model.removePhysicalGroups(
1223 [dummyContactLayerRegionDimTag]
1224 + [dummyAirRegionDimTag]
1225 + [dummyAirRegionWithoutInnerCylinderDimTag]
1226 + [dummyWindingRegionDimTag]
1227 + cuts
1228 )
1230 logger.info(
1231 "Generating Pancake3D mesh has been finished in"
1232 f" {timeit.default_timer() - start_time:.2f} s."
1233 )
1235 def structure_mesh(
1236 self,
1237 windingVolumeTags,
1238 windingSurfaceTags,
1239 windingLineTags,
1240 contactLayerVolumeTags,
1241 contactLayerSurfaceTags,
1242 contactLayerLineTags,
1243 meshSettingIndex,
1244 axialNumberOfElements=None,
1245 bumpCoefficient=None,
1246 ):
1247 """
1248 Structures the winding and contact layer meshed depending on the user inputs. If
1249 the bottom and top part of the air is to be structured, the same method is used.
1251 :param windingVolumeTags: tags of the winding volumes
1252 :type windingVolumeTags: list[int]
1253 :param windingSurfaceTags: tags of the winding surfaces
1254 :type windingSurfaceTags: list[int]
1255 :param windingLineTags: tags of the winding lines
1256 :type windingLineTags: list[int]
1257 :param contactLayerVolumeTags: tags of the contact layer volumes
1258 :type contactLayerVolumeTags: list[int]
1259 :param contactLayerSurfaceTags: tags of the contact layer surfaces
1260 :type contactLayerSurfaceTags: list[int]
1261 :param contactLayerLineTags: tags of the contact layer lines
1262 :type contactLayerLineTags: list[int]
1263 :param meshSettingIndex: index of the mesh setting
1264 :type meshSettingIndex: int
1265 :param axialNumberOfElements: number of axial elements
1266 :type axialNumberOfElements: int, optional
1267 :param bumpCoefficient: bump coefficient for axial meshing
1268 :type bumpCoefficient: float, optional
1270 """
1271 # Transfinite settings:
1272 # Arc lenght of the innermost one turn of spiral:
1273 if self.geo.ii.tsa:
1274 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1275 self.geo.wi.r_i,
1276 self.geo.wi.r_i
1277 + self.geo.wi.t
1278 + self.geo.ii.t * (self.geo.N - 1) / self.geo.N,
1279 0,
1280 2 * math.pi,
1281 )
1282 else:
1283 oneTurnSpiralLength = curve.calculateSpiralArcLength(
1284 self.geo.wi.r_i,
1285 self.geo.wi.r_i + self.geo.wi.t,
1286 0,
1287 2 * math.pi,
1288 )
1290 # Arc length of one element:
1291 arcElementLength = oneTurnSpiralLength / self.mesh.wi.ane[meshSettingIndex]
1293 # Number of azimuthal elements per turn:
1294 arcNumElementsPerTurn = round(oneTurnSpiralLength / arcElementLength)
1296 # Make all the lines transfinite:
1297 for j, lineTags in enumerate([windingLineTags, contactLayerLineTags]):
1298 for lineTag in lineTags:
1299 lineObject = curve(lineTag, self.geo)
1301 if lineObject.type is curveType.horizontal:
1302 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1303 # used.
1304 if self.geo.ii.tsa:
1305 numNodes = self.mesh.wi.rne[meshSettingIndex] + 1
1307 else:
1308 if j == 0:
1309 # This line is the winding's horizontal line:
1310 numNodes = self.mesh.wi.rne[meshSettingIndex] + 1
1312 else:
1313 # This line is the contact layer's horizontal line:
1314 numNodes = self.mesh.ii.rne[meshSettingIndex] + 1
1316 # Set transfinite curve:
1317 self.contactLayerAndWindingRadialLines.append(lineTag)
1318 gmsh.model.mesh.setTransfiniteCurve(lineTag, numNodes)
1320 elif lineObject.type is curveType.axial:
1321 # The curve is axial, so axialNumberOfElements entry is used.
1322 if axialNumberOfElements is None:
1323 numNodes = self.mesh.wi.axne[meshSettingIndex] + 1
1324 else:
1325 numNodes = axialNumberOfElements + 1
1327 if bumpCoefficient is None:
1328 bumpCoefficient = self.mesh.wi.axbc[meshSettingIndex]
1329 gmsh.model.mesh.setTransfiniteCurve(
1330 lineTag, numNodes, meshType="Bump", coef=bumpCoefficient
1331 )
1333 else:
1334 # The line is an arc, so the previously calculated arcNumElementsPerTurn
1335 # is used. All the number of elements per turn must be the same
1336 # independent of radial position. Otherwise, transfinite meshing cannot
1337 # be performed. However, to support the float number of turns, number
1338 # of nodes are being calculated depending on the start and end turns of
1339 # the arc.d
1340 lengthInTurns = abs(lineObject.n2 - lineObject.n1)
1341 if lengthInTurns > 0.5:
1342 # The arc can never be longer than half a turn.
1343 lengthInTurns = 1 - lengthInTurns
1345 lengthInTurns = (
1346 round(lengthInTurns / self.geo.wi.turnTol) * self.geo.wi.turnTol
1347 )
1349 arcNumEl = round(arcNumElementsPerTurn * lengthInTurns)
1351 arcNumNodes = int(arcNumEl + 1)
1353 # Set transfinite curve:
1354 gmsh.model.mesh.setTransfiniteCurve(lineTag, arcNumNodes)
1356 for j, surfTags in enumerate([windingSurfaceTags, contactLayerSurfaceTags]):
1357 for surfTag in surfTags:
1358 # Make all the surfaces transfinite:
1359 gmsh.model.mesh.setTransfiniteSurface(surfTag)
1361 if self.mesh.wi.elementType[meshSettingIndex] == "hexahedron":
1362 # If the element type is hexahedron, recombine all the surfaces:
1363 gmsh.model.mesh.setRecombine(2, surfTag)
1364 elif self.mesh.wi.elementType[meshSettingIndex] == "prism":
1365 # If the element type is prism, recombine only the side surfaces:
1366 surfaceNormal = list(gmsh.model.getNormal(surfTag, [0.5, 0.5]))
1367 if abs(surfaceNormal[2]) < 1e-6:
1368 gmsh.model.mesh.setRecombine(2, surfTag)
1370 # If the element type is tetrahedron, do not recombine any surface.
1372 for volTag in windingVolumeTags + contactLayerVolumeTags:
1373 # Make all the volumes transfinite:
1374 gmsh.model.mesh.setTransfiniteVolume(volTag)
1376 def structure_tubes_and_cylinders(
1377 self, volumeTags, terminalNonTubeParts=False, radialElementMultiplier=1
1378 ):
1379 # Number of azimuthal elements per quarter:
1380 arcNumElementsPerQuarter = int(self.mesh.wi.ane[0] / 4)
1381 radialNumberOfElementsPerLength = (
1382 self.mesh.wi.rne[0] / self.geo.wi.t * radialElementMultiplier
1383 )
1385 surfacesDimTags = gmsh.model.getBoundary(
1386 [(3, tag) for tag in volumeTags], combined=False, oriented=False
1387 )
1388 surfacesTags = [dimTag[1] for dimTag in surfacesDimTags]
1389 surfacesTags = list(set(surfacesTags))
1391 curvesDimTags = gmsh.model.getBoundary(
1392 surfacesDimTags, combined=False, oriented=False
1393 )
1394 curvesTags = [dimTag[1] for dimTag in curvesDimTags]
1396 # Make all the lines transfinite:
1397 for curveTag in curvesTags:
1398 curveObject = curve(curveTag, self.geo)
1400 if curveObject.type is curveType.horizontal:
1401 # The curve is horizontal, so radialNumberOfElementsPerTurn entry is
1402 # used.
1404 # But, the curve might be a part of the transitionNotch.
1405 isTransitionNotch = False
1406 point2 = curveObject.points[1]
1407 point1 = curveObject.points[0]
1408 if (
1409 abs(point2[0] - point1[0]) > 1e-5
1410 and abs(point2[1] - point1[1]) > 1e-5
1411 ):
1412 isTransitionNotch = True
1414 if isTransitionNotch:
1415 gmsh.model.mesh.setTransfiniteCurve(curveTag, 2)
1416 else:
1417 if terminalNonTubeParts:
1418 if curveTag not in self.contactLayerAndWindingRadialLines:
1419 numNodes = (
1420 round(radialNumberOfElementsPerLength * self.geo.ti.i.t)
1421 + 1
1422 )
1423 # Set transfinite curve:
1424 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1425 else:
1426 numNodes = (
1427 round(radialNumberOfElementsPerLength * curveObject.length)
1428 + 1
1429 )
1430 # Set transfinite curve:
1431 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1433 elif curveObject.type is curveType.axial:
1434 # The curve is axial, so axialNumberOfElements entry is used.
1435 if math.isclose(curveObject.length, self.geo.wi.h, rel_tol=1e-7):
1436 numNodes = min(self.mesh.wi.axne) + 1
1437 else:
1438 axialElementsPerLength = min(self.mesh.wi.axne) / self.geo.wi.h
1439 numNodes = (
1440 round(axialElementsPerLength * curveObject.length + 1e-6) + 1
1441 )
1443 gmsh.model.mesh.setTransfiniteCurve(curveTag, numNodes)
1445 else:
1446 # The line is an arc
1447 lengthInTurns = abs(curveObject.n2 - curveObject.n1)
1448 if lengthInTurns > 0.5:
1449 # The arc can never be longer than half a turn.
1450 lengthInTurns = 1 - lengthInTurns
1452 lengthInTurns = (
1453 round(lengthInTurns / self.geo.wi.turnTol) * self.geo.wi.turnTol
1454 )
1456 arcNumEl = round(arcNumElementsPerQuarter * 4 * lengthInTurns)
1458 arcNumNodes = int(arcNumEl + 1)
1460 # Set transfinite curve:
1462 gmsh.model.mesh.setTransfiniteCurve(curveTag, arcNumNodes)
1464 for surfaceTag in surfacesTags:
1465 # Make all the surfaces transfinite:
1467 if self.mesh.wi.elementType[0] == "hexahedron":
1468 # If the element type is hexahedron, recombine all the surfaces:
1469 gmsh.model.mesh.setRecombine(2, surfaceTag)
1470 elif self.mesh.wi.elementType[0] == "prism":
1471 # If the element type is prism, recombine only the side surfaces:
1472 surfaceNormal = list(gmsh.model.getNormal(surfaceTag, [0.5, 0.5]))
1473 if abs(surfaceNormal[2]) < 1e-5:
1474 gmsh.model.mesh.setRecombine(2, surfaceTag)
1476 curves = gmsh.model.getBoundary(
1477 [(2, surfaceTag)], combined=False, oriented=False
1478 )
1479 numberOfCurves = len(curves)
1480 if terminalNonTubeParts:
1481 if numberOfCurves == 4:
1482 numberOfHorizontalCurves = 0
1483 for curveTag in curves:
1484 curveObject = curve(curveTag[1], self.geo)
1485 if curveObject.type is curveType.horizontal:
1486 numberOfHorizontalCurves += 1
1488 if numberOfHorizontalCurves == 3:
1489 pass
1490 else:
1491 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1493 elif numberOfCurves == 3:
1494 pass
1495 else:
1496 curves = gmsh.model.getBoundary(
1497 [(2, surfaceTag)], combined=False, oriented=False
1498 )
1499 curveObjects = [curve(line[1], self.geo) for line in curves]
1501 divisionCurves = []
1502 for curveObject in curveObjects:
1503 if curveObject.type is curveType.horizontal:
1504 point1 = curveObject.points[0]
1505 point2 = curveObject.points[1]
1506 if not (
1507 abs(point2[0] - point1[0]) > 1e-5
1508 and abs(point2[1] - point1[1]) > 1e-5
1509 ):
1510 divisionCurves.append(curveObject)
1512 cornerPoints = (
1513 divisionCurves[0].pointTags + divisionCurves[1].pointTags
1514 )
1516 if surfaceTag not in alreadyMeshedSurfaceTags:
1517 alreadyMeshedSurfaceTags.append(surfaceTag)
1518 gmsh.model.mesh.setTransfiniteSurface(
1519 surfaceTag, cornerTags=cornerPoints
1520 )
1521 else:
1522 if numberOfCurves == 3:
1523 # Then it is a pie, corner points should be adjusted:
1524 originPointTag = None
1525 curveObject1 = curve(curves[0][1], self.geo)
1526 for point, tag in zip(curveObject1.points, curveObject1.pointTags):
1527 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1528 originPointTag = tag
1530 if originPointTag is None:
1531 curveObject2 = curve(curves[1][1], self.geo)
1532 for point, tag in zip(
1533 curveObject2.points, curveObject2.pointTags
1534 ):
1535 if math.sqrt(point[0] ** 2 + point[1] ** 2) < 1e-6:
1536 originPointTag = tag
1538 otherPointDimTags = gmsh.model.getBoundary(
1539 [(2, surfaceTag)],
1540 combined=False,
1541 oriented=False,
1542 recursive=True,
1543 )
1544 otherPointTags = [dimTag[1] for dimTag in otherPointDimTags]
1545 otherPointTags.remove(originPointTag)
1546 cornerTags = [originPointTag] + otherPointTags
1547 gmsh.model.mesh.setTransfiniteSurface(
1548 surfaceTag, cornerTags=cornerTags
1549 )
1550 else:
1551 gmsh.model.mesh.setTransfiniteSurface(surfaceTag)
1553 for volumeTag in volumeTags:
1554 if terminalNonTubeParts:
1555 surfaces = gmsh.model.getBoundary(
1556 [(3, volumeTag)], combined=False, oriented=False
1557 )
1558 curves = gmsh.model.getBoundary(
1559 surfaces, combined=False, oriented=False
1560 )
1561 curves = list(set(curves))
1563 if len(curves) == 12:
1564 numberOfArcs = 0
1565 for curveTag in curves:
1566 curveObject = curve(curveTag[1], self.geo)
1567 if curveObject.type is curveType.spiralArc:
1568 numberOfArcs += 1
1569 if numberOfArcs == 2:
1570 pass
1571 else:
1572 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1573 # elif len(curves) == 15:
1574 # divisionCurves = []
1575 # for curveTag in curves:
1576 # curveObject = curve(curveTag[1], self.geo)
1577 # if curveObject.type is curveType.horizontal:
1578 # point1 = curveObject.points[0]
1579 # point2 = curveObject.points[1]
1580 # if not (
1581 # abs(point2[0] - point1[0]) > 1e-5
1582 # and abs(point2[1] - point1[1]) > 1e-5
1583 # ):
1584 # divisionCurves.append(curveObject)
1586 # cornerPoints = (
1587 # divisionCurves[0].pointTags
1588 # + divisionCurves[1].pointTags
1589 # + divisionCurves[2].pointTags
1590 # + divisionCurves[3].pointTags
1591 # )
1592 # gmsh.model.mesh.setTransfiniteVolume(
1593 # volumeTag, cornerTags=cornerPoints
1594 # )
1595 else:
1596 # Make all the volumes transfinite:
1597 gmsh.model.mesh.setTransfiniteVolume(volumeTag)
1599 @staticmethod
1600 def get_boundaries(volumeDimTags, returnTags=False):
1601 """
1602 Returns all the surface and line dimTags or tags of a given list of volume
1603 dimTags.
1605 :param volumeDimTags: dimTags of the volumes
1606 :type volumeDimTags: list[tuple[int, int]]
1607 :param returnTags: if True, returns tags instead of dimTags
1608 :type returnTags: bool, optional
1609 :return: surface and line dimTags or tags
1610 :rtype: tuple[list[tuple[int, int]], list[tuple[int, int]]] or
1611 tuple[list[int], list[int]]
1612 """
1613 # Get the surface tags:
1614 surfaceDimTags = list(
1615 set(
1616 gmsh.model.getBoundary(
1617 volumeDimTags,
1618 combined=False,
1619 oriented=False,
1620 recursive=False,
1621 )
1622 )
1623 )
1625 # Get the line tags:
1626 lineDimTags = list(
1627 set(
1628 gmsh.model.getBoundary(
1629 surfaceDimTags,
1630 combined=False,
1631 oriented=False,
1632 recursive=False,
1633 )
1634 )
1635 )
1637 if returnTags:
1638 surfaceTags = [dimTag[1] for dimTag in surfaceDimTags]
1639 lineTags = [dimTag[1] for dimTag in lineDimTags]
1640 return surfaceTags, lineTags
1641 else:
1642 return surfaceDimTags, lineDimTags
1644 @staticmethod
1645 def fuse_volumes(volumeDimTags, fuseSurfaces=True, fusedSurfacesArePlane=False):
1646 """
1647 Fuses all the volumes in a given list of volume dimTags, removes old volumes,
1648 and returns the new volume dimTag. Also, if compundSurfacces is True, it fuses
1649 the surfaces that only belong to the volume. fusedSurfacesArePlane can be
1650 used to change the behavior of the fuse_surfaces method.
1652 :param volumeDimTags: dimTags of the volumes
1653 :type volumeDimTags: list[tuple[int, int]]
1654 :param fuseSurfaces: if True, fuses the surfaces that only belong to the
1655 volume
1656 :type fuseSurfaces: bool, optional
1657 :param fusedSurfacesArePlane: if True, fused surfaces are assumed to be
1658 plane, and fusion is performed accordingly
1659 :return: new volume's dimTag
1660 :rtype: tuple[int, int]
1661 """
1663 # Get the combined boundary surfaces:
1664 boundarySurfacesDimTags = gmsh.model.getBoundary(
1665 volumeDimTags,
1666 combined=True,
1667 oriented=False,
1668 recursive=False,
1669 )
1670 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1672 # Get all the boundary surfaces:
1673 allBoundarySurfacesDimTags = gmsh.model.getBoundary(
1674 volumeDimTags,
1675 combined=False,
1676 oriented=False,
1677 recursive=False,
1678 )
1680 # Find internal (common) surfaces:
1681 internalSurfacesDimTags = list(
1682 set(allBoundarySurfacesDimTags) - set(boundarySurfacesDimTags)
1683 )
1685 # Get the combined boundary lines:
1686 boundaryLinesDimTags = gmsh.model.getBoundary(
1687 allBoundarySurfacesDimTags,
1688 combined=True,
1689 oriented=False,
1690 recursive=False,
1691 )
1692 boundarLinesTags = [dimTag[1] for dimTag in boundaryLinesDimTags]
1694 # Get all the boundary lines:
1695 allBoundaryLinesDimTags = gmsh.model.getBoundary(
1696 allBoundarySurfacesDimTags,
1697 combined=False,
1698 oriented=False,
1699 recursive=False,
1700 )
1702 # Find internal (common) lines:
1703 internalLinesDimTags = list(
1704 set(allBoundaryLinesDimTags) - set(boundarLinesTags)
1705 )
1707 # Remove the old volumes:
1708 removedVolumeDimTags = volumeDimTags
1709 gmsh.model.occ.remove(removedVolumeDimTags, recursive=False)
1711 # Remove the internal surfaces:
1712 gmsh.model.occ.remove(internalSurfacesDimTags, recursive=False)
1714 # Remove the internal lines:
1715 gmsh.model.occ.remove(internalLinesDimTags, recursive=False)
1717 # Create a new single volume (even thought we don't use the new volume tag
1718 # directly, it is required for finding the surfaces that only belong to the
1719 # volume):
1720 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1721 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1722 newVolumeDimTag = (3, newVolumeTag)
1723 gmsh.model.occ.synchronize()
1725 if fuseSurfaces:
1726 newVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1727 (3, newVolumeTag), surfacesArePlane=fusedSurfacesArePlane
1728 )
1730 return newVolumeDimTag
1732 @staticmethod
1733 def fuse_common_surfaces_of_two_volumes(
1734 volume1DimTags, volume2DimTags, fuseOtherSurfaces=False, surfacesArePlane=False
1735 ):
1736 """
1737 Fuses common surfaces of two volumes. Volumes are given as a list of dimTags,
1738 but they are assumed to form a single volume, and this function fuses those
1739 multiple volumes into a single volume as well. If fuseOtherSurfaces is set to
1740 True, it tries to fuse surfaces that only belong to one volume too; however,
1741 that feature is not used in Pancake3D currently.
1743 :param volume1DimTags: dimTags of the first volume
1744 :type volume1DimTags: list[tuple[int, int]]
1745 :param volume2DimTags: dimTags of the second volume
1746 :type volume2DimTags: list[tuple[int, int]]
1747 :param fuseOtherSurfaces: if True, fuses the surfaces that only belong to one
1748 volume
1749 :type fuseOtherSurfaces: bool, optional
1750 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1751 fusion is performed accordingly
1752 :type surfacesArePlane: bool, optional
1753 :return: new volumes dimTags
1754 :rtype: tuple[tuple[int, int], tuple[int, int]]
1755 """
1756 vol1BoundarySurfacesDimTags = gmsh.model.getBoundary(
1757 volume1DimTags,
1758 combined=True,
1759 oriented=False,
1760 recursive=False,
1761 )
1763 vol2BoundarySurfacesDimTags = gmsh.model.getBoundary(
1764 volume2DimTags,
1765 combined=True,
1766 oriented=False,
1767 recursive=False,
1768 )
1770 # Remove the old volumes:
1771 gmsh.model.occ.remove(volume1DimTags + volume2DimTags, recursive=False)
1773 # Find common surfaces:
1774 commonSurfacesDimTags = list(
1775 set(vol2BoundarySurfacesDimTags).intersection(
1776 set(vol1BoundarySurfacesDimTags)
1777 )
1778 )
1780 # Fuse common surfaces:
1781 fusedCommonSurfaceDimTags = Mesh.fuse_surfaces(
1782 commonSurfacesDimTags, surfacesArePlane=surfacesArePlane
1783 )
1785 # Create the new volumes:
1786 for commonSurfaceDimTag in commonSurfacesDimTags:
1787 vol1BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1788 vol2BoundarySurfacesDimTags.remove(commonSurfaceDimTag)
1790 vol1BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1791 vol1BoundarySurfaceTags = [dimTag[1] for dimTag in vol1BoundarySurfacesDimTags]
1792 vol2BoundarySurfacesDimTags.extend(fusedCommonSurfaceDimTags)
1793 vol2BoundarySurfaceTags = [dimTag[1] for dimTag in vol2BoundarySurfacesDimTags]
1795 vol1SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1796 vol1BoundarySurfaceTags, sewing=True
1797 )
1798 vol1NewVolumeDimTag = (3, gmsh.model.occ.addVolume([vol1SurfaceLoop]))
1800 vol2SurfaceLoop = gmsh.model.occ.addSurfaceLoop(
1801 vol2BoundarySurfaceTags, sewing=True
1802 )
1803 vol2NewVolumeDimTag = (
1804 3,
1805 gmsh.model.occ.addVolume([vol2SurfaceLoop]),
1806 )
1808 gmsh.model.occ.synchronize()
1810 if fuseOtherSurfaces:
1811 vol1NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1812 vol1NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1813 )
1814 vol2NewVolumeDimTag = Mesh.fuse_possible_surfaces_of_a_volume(
1815 vol2NewVolumeDimTag, surfacesArePlane=surfacesArePlane
1816 )
1818 return vol1NewVolumeDimTag, vol2NewVolumeDimTag
1820 @staticmethod
1821 def fuse_possible_surfaces_of_a_volume(volumeDimTag, surfacesArePlane=False):
1822 """
1823 Fuses surfaces that only belong to the volumeDimTag.
1825 :param volumeDimTag: dimTag of the volume
1826 :type volumeDimTag: tuple[int, int]
1827 :param surfacesArePlane: if True, fused surfaces are assumed to be plane, and
1828 fusion is performed accordingly
1829 :type surfacesArePlane: bool, optional
1830 :return: new volume dimTag
1831 :rtype: tuple[int, int]
1832 """
1833 boundarySurfacesDimTags = gmsh.model.getBoundary(
1834 [volumeDimTag],
1835 combined=True,
1836 oriented=False,
1837 recursive=False,
1838 )
1839 boundarSurfacesTags = [dimTag[1] for dimTag in boundarySurfacesDimTags]
1841 # Combine surfaces that only belong to the volume:
1842 toBeFusedSurfacesDimTags = []
1843 surfacesNormals = []
1844 for surfaceDimTag in boundarySurfacesDimTags:
1845 upward, _ = gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1847 if len(list(upward)) == 1:
1848 toBeFusedSurfacesDimTags.append(surfaceDimTag)
1849 # Get the normal of the surface:
1850 surfacesNormals.append(
1851 list(gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5]))
1852 )
1854 # Remove the old volume (it is not required anymore):
1855 gmsh.model.occ.remove([volumeDimTag], recursive=False)
1856 gmsh.model.occ.synchronize()
1858 # Categorize surfaces based on their normals so that they can be combined
1859 # correctly. Without these, perpendicular surfaces will cause problems.
1861 # Define a threshold to determine if two surface normals are similar or not
1862 threshold = 1e-6
1864 # Initialize an empty list to store the sets of surfaces
1865 setsOfSurfaces = []
1867 # Calculate the Euclidean distance between each pair of objects
1868 for i in range(len(toBeFusedSurfacesDimTags)):
1869 surfaceDimTag = toBeFusedSurfacesDimTags[i]
1870 surfaceTouchingVolumeTags, _ = list(
1871 gmsh.model.getAdjacencies(surfaceDimTag[0], surfaceDimTag[1])
1872 )
1873 surfaceNormal = surfacesNormals[i]
1874 assignedToASet = False
1876 for surfaceSet in setsOfSurfaces:
1877 representativeSurfaceDimTag = surfaceSet[0]
1878 representativeSurfaceTouchingVolumeTags, _ = list(
1879 gmsh.model.getAdjacencies(
1880 representativeSurfaceDimTag[0],
1881 representativeSurfaceDimTag[1],
1882 )
1883 )
1884 representativeNormal = list(
1885 gmsh.model.getNormal(representativeSurfaceDimTag[1], [0.5, 0.5])
1886 )
1888 # Calculate the difference between surfaceNormal and
1889 # representativeNormal:
1890 difference = math.sqrt(
1891 sum(
1892 (x - y) ** 2
1893 for x, y in zip(surfaceNormal, representativeNormal)
1894 )
1895 )
1897 # Check if the distance is below the threshold
1898 if difference < threshold and set(surfaceTouchingVolumeTags) == set(
1899 representativeSurfaceTouchingVolumeTags
1900 ):
1901 # Add the object to an existing category
1902 surfaceSet.append(surfaceDimTag)
1903 assignedToASet = True
1904 break
1906 if not assignedToASet:
1907 # Create a new category with the current object if none of the
1908 # existing sets match
1909 setsOfSurfaces.append([surfaceDimTag])
1911 for surfaceSet in setsOfSurfaces:
1912 if len(surfaceSet) > 1:
1913 oldSurfaceDimTags = surfaceSet
1914 newSurfaceDimTags = Mesh.fuse_surfaces(
1915 oldSurfaceDimTags, surfacesArePlane=surfacesArePlane
1916 )
1917 newSurfaceTags = [dimTag[1] for dimTag in newSurfaceDimTags]
1919 oldSurfaceTags = [dimTag[1] for dimTag in oldSurfaceDimTags]
1920 boundarSurfacesTags = [
1921 tag for tag in boundarSurfacesTags if tag not in oldSurfaceTags
1922 ]
1923 boundarSurfacesTags.extend(newSurfaceTags)
1925 # Create a new single volume:
1926 surfaceLoop = gmsh.model.occ.addSurfaceLoop(boundarSurfacesTags, sewing=True)
1927 newVolumeTag = gmsh.model.occ.addVolume([surfaceLoop])
1928 gmsh.model.occ.synchronize()
1930 return (3, newVolumeTag)
1932 @staticmethod
1933 def fuse_surfaces(surfaceDimTags, surfacesArePlane=False, categorizeSurfaces=False):
1934 """
1935 Fuses all the surfaces in a given list of surface dimTags, removes the old
1936 surfaces, and returns the new dimTags. If surfacesArePlane is True, the surfaces
1937 are assumed to be plane, and fusing will be done without gmsh.model.occ.fuse
1938 method, which is faster.
1940 :param surfaceDimTags: dimTags of the surfaces
1941 :type surfaceDimTags: list[tuple[int, int]]
1942 :param surfacesArePlane: if True, surfaces are assumed to be plane
1943 :type surfacesArePlane: bool, optional
1944 :return: newly created surface dimTags
1945 :rtype: list[tuple[int, int]]
1946 """
1947 oldSurfaceDimTags = surfaceDimTags
1949 if surfacesArePlane:
1950 # Get the combined boundary curves:
1951 boundaryCurvesDimTags = gmsh.model.getBoundary(
1952 oldSurfaceDimTags,
1953 combined=True,
1954 oriented=False,
1955 recursive=False,
1956 )
1958 # Get all the boundary curves:
1959 allCurvesDimTags = gmsh.model.getBoundary(
1960 oldSurfaceDimTags,
1961 combined=False,
1962 oriented=False,
1963 recursive=False,
1964 )
1966 # Find internal (common) curves:
1967 internalCurvesDimTags = list(
1968 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
1969 )
1971 # Remove the old surfaces:
1972 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
1974 # Remove the internal curves:
1975 gmsh.model.occ.remove(internalCurvesDimTags, recursive=True)
1977 # Create a new single surface:
1978 def findOuterOnes(dimTags, findInnerOnes=False):
1979 """
1980 Finds the outermost surface/curve/point in a list of dimTags. The outermost means
1981 the furthest from the origin.
1982 """
1983 dim = dimTags[0][0]
1985 if dim == 2:
1986 distances = []
1987 for dimTag in dimTags:
1988 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])
1989 for curve in curves:
1990 curve = list(curve)
1991 gmsh.model.occ.synchronize()
1992 pointTags = gmsh.model.getBoundary(
1993 [(1, curveTag) for curveTag in curve],
1994 oriented=False,
1995 combined=False,
1996 )
1997 # Get the positions of the points:
1998 points = []
1999 for dimTag in pointTags:
2000 boundingbox1 = gmsh.model.occ.getBoundingBox(
2001 0, dimTag[1]
2002 )[:3]
2003 boundingbox2 = gmsh.model.occ.getBoundingBox(
2004 0, dimTag[1]
2005 )[3:]
2006 boundingbox = list(
2007 map(operator.add, boundingbox1, boundingbox2)
2008 )
2009 points.append(
2010 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2011 )
2013 distances.append(
2014 max([point[0] ** 2 + point[1] ** 2 for point in points])
2015 )
2016 elif dim == 1:
2017 distances = []
2018 for dimTag in dimTags:
2019 gmsh.model.occ.synchronize()
2020 pointTags = gmsh.model.getBoundary(
2021 [dimTag],
2022 oriented=False,
2023 combined=False,
2024 )
2025 # Get the positions of the points:
2026 points = []
2027 for dimTag in pointTags:
2028 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2029 :3
2030 ]
2031 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[
2032 3:
2033 ]
2034 boundingbox = list(
2035 map(operator.add, boundingbox1, boundingbox2)
2036 )
2037 points.append(
2038 list(map(operator.truediv, boundingbox, (2, 2, 2)))
2039 )
2041 distances.append(
2042 max([point[0] ** 2 + point[1] ** 2 for point in points])
2043 )
2045 if findInnerOnes:
2046 goalDistance = min(distances)
2047 else:
2048 goalDistance = max(distances)
2050 result = []
2051 for distance, dimTag in zip(distances, dimTags):
2052 # Return all the dimTags with the hoal distance:
2053 if math.isclose(distance, goalDistance, abs_tol=1e-6):
2054 result.append(dimTag)
2056 return result
2058 outerCurvesDimTags = findOuterOnes(boundaryCurvesDimTags)
2059 outerCurvesTags = [dimTag[1] for dimTag in outerCurvesDimTags]
2060 curveLoopOuter = gmsh.model.occ.addCurveLoop(outerCurvesTags)
2062 innerCurvesDimTags = findOuterOnes(
2063 boundaryCurvesDimTags, findInnerOnes=True
2064 )
2065 innerCurvesTags = [dimTag[1] for dimTag in innerCurvesDimTags]
2066 curveLoopInner = gmsh.model.occ.addCurveLoop(innerCurvesTags)
2068 newSurfaceTag = gmsh.model.occ.addPlaneSurface(
2069 [curveLoopOuter, curveLoopInner]
2070 )
2072 gmsh.model.occ.synchronize()
2074 return [(2, newSurfaceTag)]
2075 else:
2076 # Create a new single surface:
2077 try:
2078 fuseResults = gmsh.model.occ.fuse(
2079 [oldSurfaceDimTags[-1]],
2080 oldSurfaceDimTags[0:-1],
2081 removeObject=False,
2082 removeTool=False,
2083 )
2084 newSurfaceDimTags = fuseResults[0]
2085 except:
2086 return oldSurfaceDimTags
2088 # Get the combined boundary curves:
2089 gmsh.model.occ.synchronize()
2090 boundaryCurvesDimTags = gmsh.model.getBoundary(
2091 newSurfaceDimTags,
2092 combined=True,
2093 oriented=False,
2094 recursive=False,
2095 )
2097 # Get all the boundary curves:
2098 allCurvesDimTags = gmsh.model.getBoundary(
2099 oldSurfaceDimTags,
2100 combined=False,
2101 oriented=False,
2102 recursive=False,
2103 )
2105 # Find internal (common) curves:
2106 internalCurvesDimTags = list(
2107 set(allCurvesDimTags) - set(boundaryCurvesDimTags)
2108 )
2110 # Remove the old surfaces:
2111 gmsh.model.occ.remove(oldSurfaceDimTags, recursive=False)
2113 # Remove the internal curves:
2114 gmsh.model.occ.remove(internalCurvesDimTags, recursive=False)
2116 gmsh.model.occ.synchronize()
2118 return newSurfaceDimTags
2120 def generate_regions(self):
2121 """
2122 Generates physical groups and the regions file. Physical groups are generated in
2123 GMSH, and their tags and names are saved in the regions file. FiQuS use the
2124 regions file to create the corresponding .pro file.
2126 .vi file sends the information about geometry from geometry class to mesh class.
2127 .regions file sends the information about the physical groups formed out of
2128 elementary entities from the mesh class to the solution class.
2130 The file extension for the regions file is custom because users should not edit
2131 or even see this file.
2133 Regions are generated in the meshing part because BREP files cannot store
2134 regions.
2135 """
2136 logger.info("Generating physical groups and regions file has been started.")
2137 start_time = timeit.default_timer()
2139 # Create regions instance to both generate regions file and physical groups:
2140 self.regions = regions()
2142 # ==============================================================================
2143 # WINDING AND CONTACT LAYER REGIONS START =========================================
2144 # ==============================================================================
2145 if not self.geo.ii.tsa:
2146 windingTags = [dimTag[1] for dimTag in self.dimTags[self.geo.wi.name]]
2147 self.regions.addEntities(
2148 self.geo.wi.name, windingTags, regionType.powered, entityType.vol
2149 )
2151 insulatorTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ii.name]]
2153 terminalDimTags = (
2154 self.dimTags[self.geo.ti.i.name] + self.dimTags[self.geo.ti.o.name]
2155 )
2156 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2157 terminalDimTags, combined=False, oriented=False
2158 )
2159 transitionNotchSurfaces = gmsh.model.getBoundary(
2160 self.dimTags["innerTransitionNotch"]
2161 + self.dimTags["outerTransitionNotch"],
2162 combined=False,
2163 oriented=False,
2164 )
2166 contactLayer = []
2167 contactLayerBetweenTerminalsAndWinding = []
2168 for insulatorTag in insulatorTags:
2169 insulatorSurfaces = gmsh.model.getBoundary(
2170 [(3, insulatorTag)], combined=False, oriented=False
2171 )
2172 itTouchesTerminals = False
2173 for insulatorSurface in insulatorSurfaces:
2174 if (
2175 insulatorSurface
2176 in terminalAndNotchSurfaces + transitionNotchSurfaces
2177 ):
2178 contactLayerBetweenTerminalsAndWinding.append(insulatorTag)
2179 itTouchesTerminals = True
2180 break
2182 if not itTouchesTerminals:
2183 contactLayer.append(insulatorTag)
2185 self.regions.addEntities(
2186 self.geo.ii.name, contactLayer, regionType.insulator, entityType.vol
2187 )
2189 self.regions.addEntities(
2190 "WindingAndTerminalContactLayer",
2191 contactLayerBetweenTerminalsAndWinding,
2192 regionType.insulator,
2193 entityType.vol,
2194 )
2195 else:
2196 # Calculate the number of stacks for each individual winding. Number of
2197 # stacks is the number of volumes per turn. It affects how the regions
2198 # are created because of the TSA's pro file formulation.
2200 # find the smallest prime number that divides NofVolumes:
2201 windingDimTags = self.dimTags[self.geo.wi.name + "1"]
2202 windingTags = [dimTag[1] for dimTag in windingDimTags]
2203 NofVolumes = self.geo.wi.NofVolPerTurn
2204 smallest_prime_divisor = 2
2205 while NofVolumes % smallest_prime_divisor != 0:
2206 smallest_prime_divisor += 1
2208 # the number of stacks is the region divison per turn:
2209 NofStacks = smallest_prime_divisor
2211 # the number of sets are the total number of regions for all windings and
2212 # contact layers:
2213 NofSets = 2 * NofStacks
2215 allInnerTerminalSurfaces = gmsh.model.getBoundary(
2216 self.dimTags[self.geo.ti.i.name] + self.dimTags["innerTransitionNotch"],
2217 combined=False,
2218 oriented=False,
2219 )
2220 allInnerTerminalContactLayerSurfaces = []
2221 for innerTerminalSurface in allInnerTerminalSurfaces:
2222 normal = gmsh.model.getNormal(innerTerminalSurface[1], [0.5, 0.5])
2223 if abs(normal[2]) < 1e-5:
2224 curves = gmsh.model.getBoundary(
2225 [innerTerminalSurface], combined=False, oriented=False
2226 )
2227 curveTags = [dimTag[1] for dimTag in curves]
2228 for curveTag in curveTags:
2229 curveObject = curve(curveTag, self.geo)
2230 if curveObject.type is curveType.spiralArc:
2231 allInnerTerminalContactLayerSurfaces.append(
2232 innerTerminalSurface[1]
2233 )
2235 finalWindingSets = []
2236 finalContactLayerSets = []
2237 for i in range(NofSets):
2238 finalWindingSets.append([])
2239 finalContactLayerSets.append([])
2241 for i in range(self.geo.N):
2242 windingDimTags = self.dimTags[self.geo.wi.name + str(i + 1)]
2243 windingTags = [dimTag[1] for dimTag in windingDimTags]
2245 NofVolumes = len(windingDimTags)
2247 windings = []
2248 for windingTag in windingTags:
2249 surfaces = gmsh.model.getBoundary(
2250 [(3, windingTag)], combined=False, oriented=False
2251 )
2252 curves = gmsh.model.getBoundary(
2253 surfaces, combined=False, oriented=False
2254 )
2255 curveTags = list(set([dimTag[1] for dimTag in curves]))
2256 for curveTag in curveTags:
2257 curveObject = curve(curveTag, self.geo)
2258 if curveObject.type is curveType.spiralArc:
2259 windingVolumeLengthInTurns = abs(
2260 curveObject.n2 - curveObject.n1
2261 )
2262 if windingVolumeLengthInTurns > 0.5:
2263 # The arc can never be longer than half a turn.
2264 windingVolumeLengthInTurns = (
2265 1 - windingVolumeLengthInTurns
2266 )
2268 windings.append((windingTag, windingVolumeLengthInTurns))
2270 windingStacks = []
2271 while len(windings) > 0:
2272 stack = []
2273 stackLength = 0
2274 for windingTag, windingVolumeLengthInTurns in windings:
2275 if stackLength < 1 / NofStacks - 1e-6:
2276 stack.append(windingTag)
2277 stackLength += windingVolumeLengthInTurns
2278 else:
2279 break
2280 # remove all the windings that are already added to the stack:
2281 windings = [
2282 (windingTag, windingVolumeLengthInTurns)
2283 for windingTag, windingVolumeLengthInTurns in windings
2284 if windingTag not in stack
2285 ]
2287 # find spiral surfaces of the stack:
2288 stackDimTags = [(3, windingTag) for windingTag in stack]
2289 stackSurfacesDimTags = gmsh.model.getBoundary(
2290 stackDimTags, combined=True, oriented=False
2291 )
2292 stackCurvesDimTags = gmsh.model.getBoundary(
2293 stackSurfacesDimTags, combined=False, oriented=False
2294 )
2295 # find the curve furthest from the origin:
2296 curveObjects = []
2297 for curveDimTag in stackCurvesDimTags:
2298 curveObject = curve(curveDimTag[1], self.geo)
2299 if curveObject.type is curveType.spiralArc:
2300 curveObjectDistanceFromOrigin = math.sqrt(
2301 curveObject.points[0][0] ** 2
2302 + curveObject.points[0][1] ** 2
2303 )
2304 curveObjects.append(
2305 (curveObject, curveObjectDistanceFromOrigin)
2306 )
2308 # sort the curves based on their distance from the origin (furthest first)
2309 curveObjects.sort(key=lambda x: x[1], reverse=True)
2311 curveTags = [curveObject[0].tag for curveObject in curveObjects]
2313 # only keep half of the curveTags:
2314 furthestCurveTags = curveTags[: len(curveTags) // 2]
2316 stackSpiralSurfaces = []
2317 for surfaceDimTag in stackSurfacesDimTags:
2318 normal = gmsh.model.getNormal(surfaceDimTag[1], [0.5, 0.5])
2319 if abs(normal[2]) < 1e-5:
2320 curves = gmsh.model.getBoundary(
2321 [surfaceDimTag], combined=False, oriented=False
2322 )
2323 curveTags = [dimTag[1] for dimTag in curves]
2324 for curveTag in curveTags:
2325 if curveTag in furthestCurveTags:
2326 stackSpiralSurfaces.append(surfaceDimTag[1])
2327 break
2329 # add inner terminal surfaces too:
2330 if len(windingStacks) >= NofStacks:
2331 correspondingWindingStack = windingStacks[
2332 len(windingStacks) - NofStacks
2333 ]
2334 correspondingWindings = correspondingWindingStack[0]
2335 correspondingSurfaces = gmsh.model.getBoundary(
2336 [(3, windingTag) for windingTag in correspondingWindings],
2337 combined=True,
2338 oriented=False,
2339 )
2340 correspondingSurfaceTags = [
2341 dimTag[1] for dimTag in correspondingSurfaces
2342 ]
2343 for surface in allInnerTerminalContactLayerSurfaces:
2344 if surface in correspondingSurfaceTags:
2345 stackSpiralSurfaces.append(surface)
2347 windingStacks.append((stack, stackSpiralSurfaces))
2349 windingSets = []
2350 contactLayerSets = []
2351 for j in range(NofSets):
2352 windingTags = [
2353 windingTags for windingTags, _ in windingStacks[j::NofSets]
2354 ]
2355 windingTags = list(itertools.chain.from_iterable(windingTags))
2357 surfaceTags = [
2358 surfaceTags for _, surfaceTags in windingStacks[j::NofSets]
2359 ]
2360 surfaceTags = list(itertools.chain.from_iterable(surfaceTags))
2362 windingSets.append(windingTags)
2363 contactLayerSets.append(surfaceTags)
2365 # windingSets is a list with a length of NofSets.
2366 # finalWindingSets is also a list with a length of NofSets.
2367 for j, (windingSet, contactLayerSet) in enumerate(
2368 zip(windingSets, contactLayerSets)
2369 ):
2370 finalWindingSets[j].extend(windingSet)
2371 finalContactLayerSets[j].extend(contactLayerSet)
2373 # Seperate transition layer:
2374 terminalAndNotchSurfaces = gmsh.model.getBoundary(
2375 self.dimTags[self.geo.ti.i.name]
2376 + self.dimTags[self.geo.ti.o.name]
2377 + self.dimTags["innerTransitionNotch"]
2378 + self.dimTags["outerTransitionNotch"],
2379 combined=False,
2380 oriented=False,
2381 )
2382 terminalAndNotchSurfaceTags = set(
2383 [dimTag[1] for dimTag in terminalAndNotchSurfaces]
2384 )
2386 contactLayerSets = []
2387 terminalWindingContactLayerSets = []
2388 for j in range(NofSets):
2389 contactLayerSets.append([])
2390 terminalWindingContactLayerSets.append([])
2392 for j in range(NofSets):
2393 allContactLayersInTheSet = finalContactLayerSets[j]
2395 insulatorList = []
2396 windingTerminalInsulatorList = []
2397 for contactLayer in allContactLayersInTheSet:
2398 if contactLayer in terminalAndNotchSurfaceTags:
2399 windingTerminalInsulatorList.append(contactLayer)
2400 else:
2401 insulatorList.append(contactLayer)
2403 contactLayerSets[j].extend(set(insulatorList))
2404 terminalWindingContactLayerSets[j].extend(set(windingTerminalInsulatorList))
2406 allContactLayerSurfacesForAllPancakes = []
2407 for j in range(NofSets):
2408 # Add winding volumes:
2409 self.regions.addEntities(
2410 self.geo.wi.name + "-" + str(j + 1),
2411 finalWindingSets[j],
2412 regionType.powered,
2413 entityType.vol,
2414 )
2416 # Add insulator surfaces:
2417 self.regions.addEntities(
2418 self.geo.ii.name + "-" + str(j + 1),
2419 contactLayerSets[j],
2420 regionType.insulator,
2421 entityType.surf,
2422 )
2423 allContactLayerSurfacesForAllPancakes.extend(contactLayerSets[j])
2425 # Add terminal and winding contact layer:
2426 self.regions.addEntities(
2427 "WindingAndTerminalContactLayer" + "-" + str(j + 1),
2428 terminalWindingContactLayerSets[j],
2429 regionType.insulator,
2430 entityType.surf,
2431 )
2432 allContactLayerSurfacesForAllPancakes.extend(
2433 terminalWindingContactLayerSets[j]
2434 )
2436 allContactLayerSurfacesForAllPancakes = list(
2437 set(allContactLayerSurfacesForAllPancakes)
2438 )
2439 # Get insulator's boundary line that touches the air (required for the
2440 # pro file formulation):
2441 allContactLayerSurfacesForAllPancakesDimTags = [
2442 (2, surfaceTag) for surfaceTag in allContactLayerSurfacesForAllPancakes
2443 ]
2444 insulatorBoundary = gmsh.model.getBoundary(
2445 allContactLayerSurfacesForAllPancakesDimTags,
2446 combined=True,
2447 oriented=False,
2448 )
2449 insulatorBoundaryTags = [dimTag[1] for dimTag in insulatorBoundary]
2451 # Add insulator boundary lines:
2452 # Vertical lines should be removed from the insulator boundary because
2453 # they touch the terminals, not the air:
2454 verticalInsulatorBoundaryTags = []
2455 insulatorBoundaryTagsCopy = insulatorBoundaryTags.copy()
2456 for lineTag in insulatorBoundaryTagsCopy:
2457 lineObject = curve(lineTag, self.geo)
2458 if lineObject.type is curveType.axial:
2459 verticalInsulatorBoundaryTags.append(lineTag)
2460 insulatorBoundaryTags.remove(lineTag)
2462 # Create regions:
2463 self.regions.addEntities(
2464 self.geo.contactLayerBoundaryName,
2465 insulatorBoundaryTags,
2466 regionType.insulator,
2467 entityType.curve,
2468 )
2469 self.regions.addEntities(
2470 self.geo.contactLayerBoundaryName + "-TouchingTerminal",
2471 verticalInsulatorBoundaryTags,
2472 regionType.insulator,
2473 entityType.curve,
2474 )
2476 innerTransitionNotchTags = [
2477 dimTag[1] for dimTag in self.dimTags["innerTransitionNotch"]
2478 ]
2479 outerTransitionNotchTags = [
2480 dimTag[1] for dimTag in self.dimTags["outerTransitionNotch"]
2481 ]
2482 self.regions.addEntities(
2483 "innerTransitionNotch",
2484 innerTransitionNotchTags,
2485 regionType.powered,
2486 entityType.vol,
2487 )
2488 self.regions.addEntities(
2489 "outerTransitionNotch",
2490 outerTransitionNotchTags,
2491 regionType.powered,
2492 entityType.vol,
2493 )
2494 # ==============================================================================
2495 # WINDING AND CONTACT LAYER REGIONS ENDS =======================================
2496 # ==============================================================================
2498 # ==============================================================================
2499 # TERMINAL REGIONS START =======================================================
2500 # ==============================================================================
2502 innerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ti.i.name]]
2503 self.regions.addEntities(
2504 self.geo.ti.i.name, innerTerminalTags, regionType.powered, entityType.vol_in
2505 )
2506 outerTerminalTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ti.o.name]]
2507 self.regions.addEntities(
2508 self.geo.ti.o.name,
2509 outerTerminalTags,
2510 regionType.powered,
2511 entityType.vol_out,
2512 )
2514 # Top and bottom terminal surfaces:
2515 firstTerminalDimTags = self.dimTags[self.geo.ti.firstName]
2516 lastTerminalDimTags = self.dimTags[self.geo.ti.lastName]
2518 if self.mesh.ti.structured:
2519 topSurfaceDimTags = []
2520 for i in [1, 2, 3, 4]:
2521 lastTerminalSurfaces = gmsh.model.getBoundary(
2522 [lastTerminalDimTags[-i]], combined=False, oriented=False
2523 )
2524 topSurfaceDimTags.append(lastTerminalSurfaces[-1])
2525 else:
2526 lastTerminalSurfaces = gmsh.model.getBoundary(
2527 [lastTerminalDimTags[-1]], combined=False, oriented=False
2528 )
2529 topSurfaceDimTags = [lastTerminalSurfaces[-1]]
2530 topSurfaceTags = [dimTag[1] for dimTag in topSurfaceDimTags]
2532 if self.mesh.ti.structured:
2533 bottomSurfaceDimTags = []
2534 for i in [1, 2, 3, 4]:
2535 firstTerminalSurfaces = gmsh.model.getBoundary(
2536 [firstTerminalDimTags[-i]], combined=False, oriented=False
2537 )
2538 bottomSurfaceDimTags.append(firstTerminalSurfaces[-1])
2539 else:
2540 firstTerminalSurfaces = gmsh.model.getBoundary(
2541 [firstTerminalDimTags[-1]], combined=False, oriented=False
2542 )
2543 bottomSurfaceDimTags = [firstTerminalSurfaces[-1]]
2544 bottomSurfaceTags = [dimTag[1] for dimTag in bottomSurfaceDimTags]
2546 self.regions.addEntities(
2547 "TopSurface",
2548 topSurfaceTags,
2549 regionType.powered,
2550 entityType.surf_out,
2551 )
2552 self.regions.addEntities(
2553 "BottomSurface",
2554 bottomSurfaceTags,
2555 regionType.powered,
2556 entityType.surf_in,
2557 )
2559 # if self.geo.ii.tsa:
2560 # outerTerminalSurfaces = gmsh.model.getBoundary(
2561 # self.dimTags[self.geo.ti.o.name], combined=True, oriented=False
2562 # )
2563 # outerTerminalSurfaces = [dimTag[1] for dimTag in outerTerminalSurfaces]
2564 # innerTerminalSurfaces = gmsh.model.getBoundary(
2565 # self.dimTags[self.geo.ti.i.name], combined=True, oriented=False
2566 # )
2567 # innerTerminalSurfaces = [dimTag[1] for dimTag in innerTerminalSurfaces]
2568 # windingSurfaces = gmsh.model.getBoundary(
2569 # self.dimTags[self.geo.wi.name] + self.dimTags[self.geo.ii.name],
2570 # combined=True,
2571 # oriented=False,
2572 # )
2573 # windingSurfaces = [dimTag[1] for dimTag in windingSurfaces]
2575 # windingAndOuterTerminalCommonSurfaces = list(
2576 # set(windingSurfaces).intersection(set(outerTerminalSurfaces))
2577 # )
2578 # windingAndInnerTerminalCommonSurfaces = list(
2579 # set(windingSurfaces).intersection(set(innerTerminalSurfaces))
2580 # )
2582 # self.regions.addEntities(
2583 # "WindingAndTerminalContactLayer",
2584 # windingAndOuterTerminalCommonSurfaces
2585 # + windingAndInnerTerminalCommonSurfaces,
2586 # regionType.insulator,
2587 # entityType.surf,
2588 # )
2590 # ==============================================================================
2591 # TERMINAL REGIONS ENDS ========================================================
2592 # ==============================================================================
2594 # ==============================================================================
2595 # AIR AND AIR SHELL REGIONS STARTS =============================================
2596 # ==============================================================================
2597 airTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ai.name]]
2598 self.regions.addEntities(
2599 self.geo.ai.name, airTags, regionType.air, entityType.vol
2600 )
2602 # Create a region with two points on air to be used in the pro file formulation:
2603 # To those points, Phi=0 boundary condition will be applied to set the gauge.
2604 outerAirSurfaces = gmsh.model.getBoundary(
2605 self.dimTags[self.geo.ai.name], combined=True, oriented=False
2606 )
2607 outerAirSurface = outerAirSurfaces[-1]
2608 outerAirCurves = gmsh.model.getBoundary(
2609 [outerAirSurface], combined=True, oriented=False
2610 )
2611 outerAirCurve = outerAirCurves[-1]
2612 outerAirPoint = gmsh.model.getBoundary(
2613 [outerAirCurve], combined=False, oriented=False
2614 )
2615 outerAirPointTag = outerAirPoint[0][1]
2616 self.regions.addEntities(
2617 "OuterAirPoint",
2618 [outerAirPointTag],
2619 regionType.air,
2620 entityType.point,
2621 )
2623 innerAirSurfaces = gmsh.model.getBoundary(
2624 self.dimTags[self.geo.ai.name + "-InnerCylinder"],
2625 combined=True,
2626 oriented=False,
2627 )
2628 innerAirSurface = innerAirSurfaces[0]
2629 innerAirCurves = gmsh.model.getBoundary(
2630 [innerAirSurface], combined=True, oriented=False
2631 )
2632 innerAirCurve = innerAirCurves[-1]
2633 innerAirPoint = gmsh.model.getBoundary(
2634 [innerAirCurve], combined=False, oriented=False
2635 )
2636 innerAirPointTag = innerAirPoint[0][1]
2637 self.regions.addEntities(
2638 "InnerAirPoint",
2639 [innerAirPointTag],
2640 regionType.air,
2641 entityType.point,
2642 )
2644 if self.geo.ai.shellTransformation:
2645 if self.geo.ai.type == "cylinder":
2646 airShellTags = [
2647 dimTag[1] for dimTag in self.dimTags[self.geo.ai.shellVolumeName]
2648 ]
2649 self.regions.addEntities(
2650 self.geo.ai.shellVolumeName,
2651 airShellTags,
2652 regionType.air_far_field,
2653 entityType.vol,
2654 )
2655 elif self.geo.ai.type == "cuboid":
2656 airShell1Tags = [
2657 dimTag[1]
2658 for dimTag in self.dimTags[self.geo.ai.shellVolumeName + "-Part1"]
2659 + self.dimTags[self.geo.ai.shellVolumeName + "-Part3"]
2660 ]
2661 airShell2Tags = [
2662 dimTag[1]
2663 for dimTag in self.dimTags[self.geo.ai.shellVolumeName + "-Part2"]
2664 + self.dimTags[self.geo.ai.shellVolumeName + "-Part4"]
2665 ]
2666 self.regions.addEntities(
2667 self.geo.ai.shellVolumeName + "-PartX",
2668 airShell1Tags,
2669 regionType.air_far_field,
2670 entityType.vol,
2671 )
2672 self.regions.addEntities(
2673 self.geo.ai.shellVolumeName + "-PartY",
2674 airShell2Tags,
2675 regionType.air_far_field,
2676 entityType.vol,
2677 )
2678 # ==============================================================================
2679 # AIR AND AIR SHELL REGIONS ENDS ===============================================
2680 # ==============================================================================
2682 # ==============================================================================
2683 # CUTS STARTS ==================================================================
2684 # ==============================================================================
2685 if self.geo.ai.cutName in self.dimTags:
2686 cutTags = [dimTag[1] for dimTag in self.dimTags[self.geo.ai.cutName]]
2687 self.regions.addEntities(
2688 self.geo.ai.cutName, cutTags, regionType.air, entityType.cochain
2689 )
2691 if "CutsForPerfectInsulation" in self.dimTags:
2692 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsForPerfectInsulation"]]
2693 self.regions.addEntities(
2694 "CutsForPerfectInsulation", cutTags, regionType.air, entityType.cochain
2695 )
2697 if "CutsBetweenPancakes" in self.dimTags:
2698 cutTags = [dimTag[1] for dimTag in self.dimTags["CutsBetweenPancakes"]]
2699 self.regions.addEntities(
2700 "CutsBetweenPancakes", cutTags, regionType.air, entityType.cochain
2701 )
2702 # ==============================================================================
2703 # CUTS ENDS ====================================================================
2704 # ==============================================================================
2706 # ==============================================================================
2707 # PANCAKE BOUNDARY SURFACE STARTS ==============================================
2708 # ==============================================================================
2709 # Pancake3D Boundary Surface:
2710 allPancakeVolumes = (
2711 self.dimTags[self.geo.wi.name]
2712 + self.dimTags[self.geo.ti.i.name]
2713 + self.dimTags[self.geo.ti.o.name]
2714 + self.dimTags[self.geo.ii.name]
2715 + self.dimTags["innerTransitionNotch"]
2716 + self.dimTags["outerTransitionNotch"]
2717 )
2718 Pancake3DAllBoundary = gmsh.model.getBoundary(
2719 allPancakeVolumes, combined=True, oriented=False
2720 )
2721 Pancake3DBoundaryDimTags = list(
2722 set(Pancake3DAllBoundary)
2723 - set(topSurfaceDimTags)
2724 - set(bottomSurfaceDimTags)
2725 )
2726 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2727 self.regions.addEntities(
2728 self.geo.pancakeBoundaryName,
2729 pancake3DBoundaryTags,
2730 regionType.powered,
2731 entityType.surf,
2732 )
2734 if not self.geo.ii.tsa:
2735 # Pancake3D Boundary Surface with only winding and terminals:
2736 allPancakeVolumes = (
2737 self.dimTags[self.geo.wi.name]
2738 + self.dimTags[self.geo.ti.i.name]
2739 + self.dimTags[self.geo.ti.o.name]
2740 + self.dimTags["innerTransitionNotch"]
2741 + self.dimTags["outerTransitionNotch"]
2742 + [(3, tag) for tag in contactLayerBetweenTerminalsAndWinding]
2743 )
2744 Pancake3DAllBoundary = gmsh.model.getBoundary(
2745 allPancakeVolumes, combined=True, oriented=False
2746 )
2747 Pancake3DBoundaryDimTags = list(
2748 set(Pancake3DAllBoundary)
2749 - set(topSurfaceDimTags)
2750 - set(bottomSurfaceDimTags)
2751 )
2752 pancake3DBoundaryTags = [dimTag[1] for dimTag in Pancake3DBoundaryDimTags]
2753 self.regions.addEntities(
2754 self.geo.pancakeBoundaryName + "-OnlyWindingAndTerminals",
2755 pancake3DBoundaryTags,
2756 regionType.powered,
2757 entityType.surf,
2758 )
2760 # ==============================================================================
2761 # PANCAKE BOUNDARY SURFACE ENDS ================================================
2762 # ==============================================================================
2764 # Generate regions file:
2765 self.regions.generateRegionsFile(self.regions_file)
2766 self.rm = FilesAndFolders.read_data_from_yaml(self.regions_file, RegionsModel)
2768 logger.info(
2769 "Generating physical groups and regions file has been finished in"
2770 f" {timeit.default_timer() - start_time:.2f} s."
2771 )
2773 def generate_mesh_file(self):
2774 """
2775 Saves mesh file to disk.
2778 """
2779 logger.info(
2780 f"Generating Pancake3D mesh file ({self.mesh_file}) has been started."
2781 )
2782 start_time = timeit.default_timer()
2784 gmsh.write(self.mesh_file)
2786 logger.info(
2787 f"Generating Pancake3D mesh file ({self.mesh_file}) has been finished"
2788 f" in {timeit.default_timer() - start_time:.2f} s."
2789 )
2791 if self.mesh_gui:
2792 gmsh.option.setNumber("Geometry.Volumes", 0)
2793 gmsh.option.setNumber("Geometry.Surfaces", 0)
2794 gmsh.option.setNumber("Geometry.Curves", 0)
2795 gmsh.option.setNumber("Geometry.Points", 0)
2796 self.gu.launch_interactive_GUI()
2797 else:
2798 gmsh.clear()
2799 gmsh.finalize()
2801 def load_mesh(self):
2802 """
2803 Loads mesh from .msh file.
2806 """
2807 logger.info("Loading Pancake3D mesh has been started.")
2808 start_time = timeit.default_timer()
2810 previousGeo = FilesAndFolders.read_data_from_yaml(
2811 self.geometry_data_file, Pancake3DGeometry
2812 )
2813 previousMesh = FilesAndFolders.read_data_from_yaml(
2814 self.mesh_data_file, Pancake3DMesh
2815 )
2817 if previousGeo.dict() != self.geo.dict():
2818 raise ValueError(
2819 "Geometry data has been changed. Please regenerate the geometry or load"
2820 " the previous geometry data."
2821 )
2822 elif previousMesh.dict() != self.mesh.dict():
2823 raise ValueError(
2824 "Mesh data has been changed. Please regenerate the mesh or load the"
2825 " previous mesh data."
2826 )
2828 gmsh.clear()
2829 gmsh.open(self.mesh_file)
2831 logger.info(
2832 "Loading Pancake3D mesh has been finished in"
2833 f" {timeit.default_timer() - start_time:.2f} s."
2834 )