Coverage for fiqus/geom_generators/GeometryPancake3D.py: 82%

1import math

2import logging

3from enum import Enum

4from inspect import currentframe, getframeinfo

5from typing import List, Tuple, Dict

6import operator

8import os

9import json

10import timeit

11import numpy as np

12import gmsh

14from fiqus.utils.Utils import GmshUtils

15from fiqus.utils.Utils import FilesAndFolders

16from fiqus.mains.MainPancake3D import Base

17from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry

19logger = logging.getLogger(__name__)

21# coordinateList = []

23def findSurfacesWithNormalsOnXYPlane(dimTags):

24 result = []

25 for dimTag in dimTags:

26 surfaceNormal = gmsh.model.getNormal(dimTag[1], [0.5, 0.5])

27 if abs(surfaceNormal[2]) < 1e-6:

28 result.append(dimTag)

30 return result

33def findOuterOnes(dimTags, findInnerOnes=False):

34 """

35 Finds the outermost surface/curve/point in a list of dimTags. The outermost means

36 the furthest from the origin.

37 """

38 dim = dimTags[0][0]

39 if dim == 2:

40 distances = []

41 for dimTag in dimTags:

42 _, curves = gmsh.model.occ.getCurveLoops(dimTag[1])

43 for curve in curves:

44 curve = list(curve)

45 pointTags = gmsh.model.getBoundary(

46 [(1, curveTag) for curveTag in curve],

47 oriented=False,

48 combined=False,

49 )

50 # Get the positions of the points:

51 points = []

52 for dimTag in pointTags:

53 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]

54 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]

55 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))

56 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))

58 distances.append(

59 max([point[0] ** 2 + point[1] ** 2 for point in points])

60 )

61 elif dim == 1:

62 distances = []

63 for dimTag in dimTags:

64 pointTags = gmsh.model.getBoundary(

65 [dimTag],

66 oriented=False,

67 combined=False,

68 )

69 # Get the positions of the points:

70 points = []

71 for dimTag in pointTags:

72 boundingbox1 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[:3]

73 boundingbox2 = gmsh.model.occ.getBoundingBox(0, dimTag[1])[3:]

74 boundingbox = list(map(operator.add, boundingbox1, boundingbox2))

75 points.append(list(map(operator.truediv, boundingbox, (2, 2, 2))))

77 distances.append(max([point[0] ** 2 + point[1] ** 2 for point in points]))

79 if findInnerOnes:

80 goalDistance = min(distances)

81 else:

82 goalDistance = max(distances)

84 result = []

85 for distance, dimTag in zip(distances, dimTags):

86 # Return all the dimTags with the hoal distance:

87 if math.isclose(distance, goalDistance, abs_tol=1e-6):

88 result.append(dimTag)

90 return result

93class dimTags:

94 """

95 This is a class for (dim, tag) tuple lists. DimTags are heavily used in GMSH, and

96 dimTags class makes it easier to store and manipulate them.

98 Every dimTags instance with the save option True will be stored in the

99 dimTagsStorage class. dimTags instances with the save option False will not be

100 stored. dimTagsStorage class stores all the dimTags with the corresponding names so

101 that later the dimTagsStorage class can be used to create the volume information

102 file. The volume information file is required in the meshing stage to create the

103 appropriate regions.

104

105 If parentName is specified, during the volume information file generation, another

106 key that is equal to parentName is created, and all the dimTags are also added

107 there. For example, air consists of many parts such as gap, outer tube, inner

108 cylinder parts, and their dimTags object all have the parentName of self.geo.air.name

109 (user input of the air name) so that in the volume information file, they are all

110 combined under the air name as well.

111

112 :param name: name of the dimTags object (default: None)

113 :type name: str, optional

114 :param parentName: name of the parent dimTags object (default: None)

115 :type parentName: str, optional

116 :param save: True if the instance to be stored in dimTagsStorage class, False

117 otherwise (default: False)

118 :type save: bool, optional

119 """

120

121 point = 0

122 curve = 1

123 surface = 2

124 volume = 3

125

126 def storageUpdateRequired(func):

127 """

128 A decorator for dimTags class. It will update the dimTagsStorage class if the

129 save option is True and the decorated function is called. Use this decorator

130 for every function that changes the dimTags instance so that the storage gets

131 updated.

132

133 :param func: function to be decorated

134 :type func: function

135 :return: decorated function

136 :rtype: function

137 """

138

139 def wrapper(self, *args, **kwargs):

140 func(self, *args, **kwargs)

141

142 if self.save:

143 # Update the dimTagsStorage:

144 dimTagsStorage.updateDimTags(self)

145

146 return wrapper

147

148 def __init__(

149 self,

150 name: str = None,

151 parentName: str = None,

152 save: bool = False,

153 ):

154 self.name = name

155

156 dimTagsObjects = dimTagsStorage.getDimTagsObject(name)

157 if dimTagsObjects != []:

158 dimTagsObject = dimTagsObjects[0]

159 self.physicalTag = dimTagsObject.physicalTag

160 # To store points, curves, surfaces, and volumes separately:

161 self.dimTags = dimTagsObject.dimTags

162 self.dimTagsForPG = dimTagsObject.dimTagsForPG

163 self.allDimTags = dimTagsObject.allDimTags

164 self.parentName = dimTagsObject.parentName

165 self.save = dimTagsObject.save

166 else:

167 self.physicalTag = None

168 # To store points, curves, surfaces, and volumes separately:

169 self.dimTags = [[] for _ in range(4)]

170 self.dimTagsForPG = [[] for _ in range(4)] # dim tags for physical groups

171 self.allDimTags = []

172 self.parentName = parentName

173 self.save = save

174 if self.save:

175 dimTagsStorage.updateDimTags(self)

176

177 @storageUpdateRequired

178 def add(

179 self,

180 dimTagsList: List[Tuple[int, int]],

181 dimTagsListForPG: List[Tuple[int, int]] = None,

182 ):

183 """

184 Adds a list of (dim, tag) tuples to the dimTags object.

185

186 dimTagsListForPG is also accepted as an argument because sometimes, the stored

187 dimTags and the current dimTags in the geometry generation can be different. For

188 example, if volume 61 is deleted, the other volume tags (62, 63, ...) won't

189 shift back. However, after saving the geometry as a BREP file and rereading it,

190 the volume tags will be shifted back. In this case, the stored dimTags should be

191 shifted as well. But to create the correct physical region in the

192 geometry-creating process (which won't be saved in the BREP file, just used for

193 debugging purposes), another optional dimTagsListForPG argument is accepted.

194

195 :param dimTagsList: list of (dim, tag) tuples

196 :type dimTagsList: list[tuple[int, int]]

197 :param dimTagsListForPG: list of (dim, tag) tuples for physical groups

198 (default: None). If dimTagsListForPG is None, dimTagsList will be used for

199 physical groups as well.

200 :type dimTagsListForPG: list[tuple[int, int]], optional

201

202 """

203 if not isinstance(dimTagsList, list):

204 dimTagsList = [dimTagsList]

205

206 if not all(isinstance(element, tuple) for element in dimTagsList):

207 raise TypeError("Dim tags must be a list of tuples!")

208

209 # Sometimes, negative entities can be added for topology.

210 for i, v in enumerate(dimTagsList):

211 if v[1] < 0:

212 dimTagsList[i] = (v[0], -v[1])

213

214 # Add dim tags if they are not already added:

215 for v in dimTagsList:

216 if v not in self.allDimTags:

217 self.dimTags[v[0]].append(v)

218 if dimTagsListForPG is None:

219 self.dimTagsForPG[v[0]].append(v)

220 self.allDimTags.append(v)

221

222 if dimTagsListForPG is not None:

223 if not isinstance(dimTagsListForPG, list):

224 dimTagsListForPG = [dimTagsListForPG]

225

226 if not all(isinstance(element, tuple) for element in dimTagsListForPG):

227 raise TypeError("Dim tags must be a list of tuples!")

228

229 for i, v in enumerate(dimTagsListForPG):

230 if v[1] < 0:

231 dimTagsListForPG[i] = (v[0], -v[1])

232 for v in dimTagsListForPG:

233 if v not in self.dimTagsForPG:

234 self.dimTagsForPG[v[0]].append(v)

235

236 def addWithTags(self, dim: int, tags: List[int], tagsForPG: List[int] = None):

237 """

238 Adds a list of tags with a specific dimension to the dimTags object. The

239 explanation of the tagsForPG argument is given in the add() method.

240

241 :param dim: dimension of the tags

242 :type dim: int

243 :param tags: list of tags

244 :type tags: list[tuple[int, int]]

245 :param tagsForPG: list of tags for physical groups (default: None)

246 :type tagsForPG: list[tuple[int, int]], optional

247

248 """

249 if not isinstance(tags, list):

250 tags = [tags]

251

252 if not isinstance(dim, int):

253 raise TypeError("Dimension must be an integer!")

254

255 if not all(isinstance(element, int) for element in tags):

256 raise TypeError("Tags must be a list of integers!")

257

258 dims = [dim] * len(tags)

259 dimTagsList = list(zip(dims, tags))

260

261 if tagsForPG is not None:

262 if not isinstance(tagsForPG, list):

263 tagsForPG = [tagsForPG]

264

265 if not all(isinstance(element, int) for element in tagsForPG):

266 raise TypeError("Tags must be a list of integers!")

267

268 dimTagsListForPG = list(zip(dims, tagsForPG))

269 self.add(dimTagsList, dimTagsListForPG)

270 else:

271 self.add(dimTagsList)

272

273 def getDimTags(self, dim: int = None) -> List[Tuple[int, int]]:

274 """

275 Returns the stored list of (dim, tag) tuples with a specific dimension. If dim

276 is not specified, returns all the (dim, tag) tuples.

277

278 :param dim: dimension of the tags to be returned (default: None)

279 :type dim: int, optional

280 :return: list of (dim, tag) tuples

281 :rtype: list[tuple[int, int]]

282 """

283 if dim is None:

284 return self.dimTags[0] + self.dimTags[1] + self.dimTags[2] + self.dimTags[3]

285 else:

286 return self.dimTags[dim]

287

288 def getDimTagsForPG(self, dim: int = None) -> List[Tuple[int, int]]:

289 """

290 Returns the stored list of (dim, tag) tuples for physical groups with a specific

291 dimension. If dim is not specified, returns all the (dim, tag) tuples for

292 physical groups.

293

294 :param dim: dimension of the tags to be returned (default: None)

295 :type dim: int, optional

296 :return: list of (dim, tag) tuples for physical groups

297 :rtype: list[tuple[int, int]]

298 """

299 if dim is None:

300 return (

301 self.dimTagsForPG[0]

302 + self.dimTagsForPG[1]

303 + self.dimTagsForPG[2]

304 + self.dimTagsForPG[3]

305 )

306 else:

307 return self.dimTagsForPG[dim]

308

309 def getTags(self, dim: int, forPhysicalGroup=False) -> List[int]:

310 """

311 Returns the stored list of tags with a specific dimension.

312

313 :param dim: dimension of the tags to be returned

314 :type dim: int

315 :return: list of tags

316 :rtype: list[int]

317 """

318 if forPhysicalGroup:

319 return [v[1] for v in self.dimTagsForPG[dim]]

320 else:

321 return [v[1] for v in self.dimTags[dim]]

322

323 def getExtrusionTop(self, dim=3):

324 """

325 Returns the top surfaces, lines, or points of an extrusion operation if the

326 dimTags object contains the tags of an extrusion operation.

327 gmsh.model.occ.extrusion() function returns all the entities that are created

328 by the extrusion as a dim tags list. The first element is always the top

329 surface, the second is the volume. However, when more than one surface is

330 extruded, extracting the top surfaces is not trivial. This function returns

331 the top surfaces of an extrusion operation. It does that by finding the entities

332 right before the volume entities. If dim is 2, the function will return the top

333 curves of an extrusion operation. If dim is 1, the function will return the top

334 points of an extrusion.

335

336 :param dim: dimension of the entity that is being created (default: 3)

337 :type dim: int, optional

338 :return: list of (dim, tag) tuples of the top entities of an extrusion operation

339 :rtype: list[tuple[int, int]]

340 """

341

342 topSurfaces = []

343 for index, dimTag in enumerate(self.allDimTags):

344 if dimTag[0] == dim:

345 topSurfaces.append(self.allDimTags[index - 1])

346

347 return topSurfaces

348

349 def getExtrusionSide(self, dim=3):

350 """

351 Returns the side surfaces, lines, or points of an extrusion operation if the

352 dimTags object contains the tags of an extrusion operation.

353 gmsh.model.occ.extrusion() function returns all the entities that are created

354 by the extrusion as a dim tags list. The first element is always the top

355 surface, the second is the volume. The other elements are the side surfaces.

356 However, when more than one surface is extruded, extracting the side surfaces

357 is not trivial. This function returns the side surfaces of an extrusion

358 operation. It does that by finding returning all the entities except the top

359 surface and the volume. If dim is 2, the function will return the side curves of

360 an extrusion operation.

361

362 :param dim: dimension of the entity that is being created (default: 3)

363 :type dim: int, optional

364 :return: list of (dim, tag) tuples of the side entities of an extrusion operation

365 :rtype: list[tuple[int, int]]

366 """

367 sideSurfaces = []

368 sideSurfaceStartIndex = None

369 for index, dimTag in enumerate(self.allDimTags):

370 if dimTag[0] == dim:

371 if sideSurfaceStartIndex is not None:

372 sideSurfaces.append(

373 self.allDimTags[sideSurfaceStartIndex : index - 1]

374 )

375 sideSurfaceStartIndex = index + 1

376 else:

377 sideSurfaceStartIndex = index + 1

378

379 sideSurfaces.append(self.allDimTags[sideSurfaceStartIndex:])

380

381 return sideSurfaces

382

383 def __add__(self, other):

384 """

385 Adds two dimTags objects and returns a new dimTags object with the same save and

386 name attirbues of the first dimTags object.

387

388 It might cause bugs because of the recursive behavior of the

389 @storageUpdateRequired decorator. Use with caution. Currently only used by

390 dimTagsStorage.updateDimTags method.

391

392 :param other: dimTags object to be added

393 :type other: dimTags

394 :return: dimTags object with the sum of the two dimTags objects

395 :rtype: dimTags

396 """

397 result = dimTags()

398 result.name = self.name

399 result.parentName = self.parentName

400 result.physicalTag = self.physicalTag

401 result.dimTags = self.dimTags

402 result.dimTagsForPG = self.dimTagsForPG

403 result.allDimTags = self.allDimTags

404

405 result.add(other.allDimTags)

406 result.save = self.save

407 return result

408

409 def __repr__(self):

410 """

411 Returns the string representation of the dimTags object. If the dimTags object

412 is saved, it will return "SAVED: name". If the dimTags object is not saved, it

413 will return "NOT SAVED: name".

414

415 dimTags objects are used as dictionary keys throughout the code. This

416 representation makes debugging easier.

417

418 :return: string representation of the dimTags object

419 :rtype: str

420 """

421 if self.save:

422 return "SAVED: " + self.name

423 else:

424 return "NOT SAVED: " + self.name

425

426

427class dimTagsStorage:

428 """

429 This is a global class to store the dimTags of important entities in the model.

430 Every dimTags instance with self.save = True will be stored in this class. Later,

431 the storage will be used to generate the volume information (*.vi) file. *.vi file

432 will be used for generating the physical regions in the meshing part.

433

434 Users should not use this class directly. Instead, they should use the dimTags

435 class. If they assign save = True to the dimTags instance, it will be stored in this

436 class.

437

438 This class is a singleton class. It means that there will be only one instance of

439 this class in the whole module. This is done to be able to use the same storage

440 throughout this module. See the singleton design pattern for more information.

441 """

442

443 __instance = None

444 __dimTagsDict = {} # Dictionary with the names of the dimTags objects as keys and

445 # dimTags objects as values

446

447 def __new__(cls):

448 if cls.__instance is None:

449 cls.__instance = super().__new__(cls)

450

451 return cls.__instance

452

453 @classmethod

454 def updateDimTags(cls, dimTagsObject: dimTags):

455 """

456 Either adds or updates the dimTags object in the storage.

457

458 :param dimTags: dimTags object to be added or updated.

459 :type dimTags: dimTags

460

461 """

462 if dimTagsObject.name in cls.__dimTagsDict:

463 newDimTags = dimTagsObject + cls.__dimTagsDict[dimTagsObject.name]

464 cls.__dimTagsDict[dimTagsObject.name] = newDimTags

465 else:

466 cls.__dimTagsDict[dimTagsObject.name] = dimTagsObject

467

468 # @classmethod

469 # def updateDimTagsFromDict(cls, name: str, dimTagsList: List[Tuple[int, int]]):

470 # """

471 # Updates or adds dimTags from a list of (dim, tag) tuples.

472 #

473 # :param name: Name of the dimTags entry to update or add.

474 # :type name: str

475 # :param dimTagsList: List of (dim, tag) tuples to be associated with the name.

476 # :type dimTagsList: List[Tuple[int, int]]

477 # """

478 # # Check if the entry exists; if so, update it, otherwise add a new entry

479 # if name in cls.__dimTagsDict:

480 # existingDimTags = cls.__dimTagsDict[

481 # name].dimTags # Assuming dimTags object has a dimTags attribute storing the list of tuples

482 # updatedDimTags = existingDimTags + dimTagsList

483 # cls.__dimTagsDict[name].dimTags = updatedDimTags # Update the list of tuples

484 # else:

485 # # Create a new dimTags object (this step depends on how your dimTags class is structured)

486 # newDimTagsObject = dimTags(name=name, dimTags=dimTagsList) # Assuming such a constructor exists

487 # cls.__dimTagsDict[name] = newDimTagsObject

488

489 @classmethod

490 def getDimTagsObject(cls, names: List[str]):

491 """

492 Returns the dimTags object with the given names.

493

494 :param names: names of the dimTags objects.

495 :type names: list[str]

496 :return: dimTags objects with the given name.

497 :rtype: list[dimTags]

498 """

499 if not isinstance(names, list):

500 names = [names]

501

502 dimTagsObjects = []

503 for name in names:

504 if name in cls.__dimTagsDict:

505 dimTagsObjects.append(cls.__dimTagsDict[name])

506

507 return dimTagsObjects

508

509 @classmethod

510 def getDimTags(cls, names: List[str], dim: int = None) -> List[Tuple[int, int]]:

511 """

512 Returns the stored list of (dim, tag) tuples with dimension a specific dimenions

513 and names. If dim is not specified, all the stored (dim, tag) tuples under the

514 given names will be returned.

515

516 :param names: names of the dimTags object that will be returned

517 :type names: list[str]

518 :param dim: dimension of the (dim, tag) tuples to be returned (default: None).

519 If dim is None, all the stored (dim, tag) tuples under the name names will

520 be returned.

521 :type dim: int, optional

522 :return: list of (dim, tag) tuples

523 """

524 if not isinstance(names, list):

525 names = [names]

526

527 dimTagsResult = []

528 for name in names:

529 dimTagsResult.extend(cls.__dimTagsDict[name].getDimTags(dim))

530

531 return dimTagsResult

532

533 @classmethod

534 def getTags(cls, names: List[str], dim: int) -> List[int]:

535 """

536 Returns the stored list of tags with dimension a specific dimension and names.

537

538 :param names: names of the dimTags objects

539 :type names: list[str]

540 :param dim: dimension of the tags to be returned

541 :type dim: int

542 :return: list of tags

543 :rtype: list[int]

544 """

545 dimTags = cls.getDimTags(names, dim)

546 tags = [dimTag[1] for dimTag in dimTags]

547

548 return tags

549

550 @classmethod

551 def getDimTagsDict(

552 cls, forPhysicalGroups=False

553 ) -> Dict[str, List[Tuple[int, int]]]:

554 """

555 Returns a dictionary with the names of the dimTags objects as keys and the

556 stored list of (dim, tag) tuples as values. This method is used to generate the

557 .vi file. If forPhysicalGroups is True, the dimTags for physical groups will be

558 returned instead of the dimTags.

559

560 :param forPhysicalGroups: True if the dimTags for physical groups should be

561 returned, False otherwise (default: False)

562 :type forPhysicalGroups: bool, optional

563 :return: dictionary with the names of the dimTags objects as keys and the

564 stored list of (dim, tag) tuples as values

565 :rtype: dict[str, list[tuple[int, int]]]

566 """

567 dictionary = {}

568 for name, dimTagsObject in cls.__dimTagsDict.items():

569 if dimTagsObject.parentName is not None:

570 if dimTagsObject.parentName in dictionary:

571 if forPhysicalGroups:

572 dictionary[dimTagsObject.parentName].extend(

573 dimTagsObject.getDimTagsForPG()

574 )

575 else:

576 dictionary[dimTagsObject.parentName].extend(

577 dimTagsObject.getDimTags()

578 )

579 else:

580 if forPhysicalGroups:

581 dictionary[dimTagsObject.parentName] = (

582 dimTagsObject.getDimTagsForPG()

583 )

584 else:

585 dictionary[dimTagsObject.parentName] = (

586 dimTagsObject.getDimTags()

587 )

588 if forPhysicalGroups:

589 dictionary[name] = dimTagsObject.getDimTagsForPG()

590 else:

591 dictionary[name] = dimTagsObject.getDimTags()

592

593 return dictionary

594

595 @classmethod

596 def getAllStoredDimTags(cls) -> List[Tuple[int, int]]:

597 """

598 Returns a list of all the stored (dim, tag) tuples, regardless of the name of

599 the dimTags object (i.e. all the dimTags objects are merged into one list).

600

601 :return: list of all the stored (dim, tag) tuples.

602 :rtype: list[tuple[int, int]]

603 """

604 AllStoredDimTags = []

605 for name, dimTagsObject in cls.__dimTagsDict.items():

606 AllStoredDimTags.extend(dimTagsObject.getDimTags())

607

608 return AllStoredDimTags

609

610 @classmethod

611 def clear(cls):

612 """

613 Clears the dimTagsStorage class.

614

615

616 """

617 cls.__instance = None

618 cls.__dimTagsDict = (

619 {}

620 ) # Dictionary with the names of the dimTags objects as keys and

621 # dimTags objects as values

622

623

624class coordinate(Enum):

625 """

626 A class to specify coordinate types easily.

627 """

628

629 rectangular = 0

630 cylindrical = 1

631 spherical = 2

632

633

634class direction(Enum):

635 """

636 A class to specify direction easily.

637 """

638

639 ccw = 0

640 cw = 1

641

642

643class point:

644 """

645 This is a class for creating points in GMSH. It supports rectangular and cylindrical

646 coordinates. Moreover, vector operations are supported.

647

648 :param r0: x, r, or r (default: 0.0)

649 :type r0: float, optional

650 :param r1: y, theta, or theta (default: 0.0)

651 :type r1: float, optional

652 :param r2: z, z, or phi (default: 0.0)

653 :type r2: float, optional

654 :param type: coordinate type (default: coordinate.rectangular)

655 :type type: coordinate, optional

656 """

657

658 def __init__(self, r0=0.0, r1=0.0, r2=0.0, type=coordinate.rectangular) -> None:

659

660 self.type = type # Store 'type' as an instance attribute

661

662 if type is coordinate.rectangular:

663 self.x = r0

664 self.y = r1

665 self.z = r2

666

667 self.r = math.sqrt(self.x**2 + self.y**2)

668 self.theta = math.atan2(self.y, self.x)

669 elif type is coordinate.cylindrical:

670 self.r = r0

671 self.theta = r1

672 self.x = self.r * math.cos(self.theta)

673 self.y = self.r * math.sin(self.theta)

674 self.z = r2

675 elif type is coordinate.spherical:

676 raise ValueError("Spherical coordinates are not supported yet!")

677 else:

678 raise ValueError("Improper coordinate type value!")

679

680 self.tag = gmsh.model.occ.addPoint(self.x, self.y, self.z)

681

682 def __repr__(self):

683 """

684 Returns the string representation of the point.

685

686 :return: string representation of the point

687 :rtype: str

688 """

689 return "point(%r, %r, %r, %r)" % (self.x, self.y, self.z, self.type)

690

691 def __abs__(self):

692 """

693 Returns the magnitude of the point vector.

694

695 :return: the magnitude of the point vector

696 :rtype: float

697 """

698 return math.hypot(self.x, self.y, self.z)

699

700 def __add__(self, other):

701 """

702 Returns the summation of two point vectors.

703

704 :param other: point vector to be added

705 :type other: point

706 :return: the summation of two point vectors

707 :rtype: point

708 """

709 x = self.x + other.x

710 y = self.y + other.y

711 z = self.z + other.z

712 return point(x, y, z, coordinate.rectangular)

713

714 def __mul__(self, scalar):

715 """

716 Returns the product of a point vector and a scalar.

717

718 :param scalar: a scalar value

719 :type scalar: float

720 :return: point

721 :rtype: point

722 """

723 return point(

724 self.x * scalar,

725 self.y * scalar,

726 self.z * scalar,

727 coordinate.rectangular,

728 )

729

730

731class spiralCurve():

732 """

733 A class to create a spiral curves parallel to XY plane in GMSH. The curve is defined

734 by a spline and it is divided into sub-curves. Sub-curves are used because it makes

735 the geometry creation process easier.

736

737 :param innerRadius: inner radius

738 :type innerRadius: float

739 :param gap: gap after each turn

740 :type gap: float

741 :param turns: number of turns

742 :type turns: float

743 :param z: z coordinate

744 :type z: float

745 :param initialTheta: initial theta angle in radians

746 :type initialTheta: float

747 :param direction: direction of the spiral (default: direction.ccw)

748 :type direction: direction, optional

749 :param cutPlaneNormal: normal vector of the plane that will cut the spiral curve

750 (default: None)

751 :type cutPlaneNormal: tuple[float, float, float], optional

752 """

753

754 # If the number of points used per turn is n, then the number of sections per turn

755 # is n-1. They set the resolution of the spiral curve. It sets the limit of the

756 # precision of the float number of turns that can be used to create the spiral

757 # curve. The value below might be modified in Geometry.__init__ method.

758 sectionsPerTurn = 16

759

760 # There will be curvesPerTurn curve entities per turn. It will be effectively the

761 # number of volumes per turn in the end. The value below might be modified in

762 # Geometry.__init__ method.

763 curvesPerTurn = 2

764

765 def __init__(

766 self,

767 innerRadius,

768 gap,

769 turns,

770 z,

771 initialTheta,

772 transitionNotchAngle,

773 direction=direction.ccw,

774 cutPlaneNormal=Tuple[float, float, float],

775 ) -> None:

776 spt = self.sectionsPerTurn # just to make the code shorter

777 self.turnRes = 1 / spt # turn resolution

778 cpt = self.curvesPerTurn # just to make the code shorter

779 self.turns = turns

780

781 # =============================================================================

782 # GENERATING POINTS STARTS ====================================================

783 # =============================================================================

784

785 # Calculate the coordinates of the points that define the spiral curve:

786 if direction is direction.ccw:

787 # If the spiral is counter-clockwise, the initial theta angle decreases,

788 # and r increases as the theta angle decreases.

789 multiplier = 1

790 elif direction is direction.cw:

791 # If the spiral is clockwise, the initial theta angle increases, and r

792 # increases as the theta angle increases.

793 multiplier = -1

794

795 NofPointsPerTurn = int(spt + 1)

796 thetaArrays = []

797 for turn in range(1, int(self.turns) + 1):

798 thetaArrays.append(

799 np.linspace(

800 initialTheta + (turn - 1) * 2 * math.pi * multiplier,

801 initialTheta + (turn) * 2 * math.pi * multiplier,

802 NofPointsPerTurn,

803 )

804 )

805

806 thetaArrays.append(

807 np.linspace(

808 initialTheta + (turn) * 2 * math.pi * multiplier,

809 initialTheta + (self.turns) * 2 * math.pi * multiplier,

810 round(spt * (self.turns - turn) + 1),

811 )

812 )

813

814 if cutPlaneNormal is not None:

815 # If the cutPlaneNormal is specified, the spiral curve will be cut by a

816 # plane that is normal to the cutPlaneNormal vector and passes through the

817 # origin.

818

819 alpha = math.atan2(cutPlaneNormal[1], cutPlaneNormal[0]) - math.pi / 2

820 alpha2 = alpha + math.pi

821

822 listOfBreakPoints = []

823 for turn in range(1, int(self.turns) + 2):

824 breakPoint1 = alpha + (turn - 1) * 2 * math.pi * multiplier

825 breakPoint2 = alpha2 + (turn - 1) * 2 * math.pi * multiplier

826 if (

827 breakPoint1 > initialTheta

828 and breakPoint1 < initialTheta + 2 * math.pi * self.turns

829 ):

830 listOfBreakPoints.append(breakPoint1)

831 if (

832 breakPoint2 > initialTheta

833 and breakPoint2 < initialTheta + 2 * math.pi * self.turns

834 ):

835 listOfBreakPoints.append(breakPoint2)

836

837 thetaArrays.append(np.array(listOfBreakPoints))

838

839 theta = np.concatenate(thetaArrays)

840 theta = np.round(theta, 10)

841 theta = np.unique(theta)

842 theta = np.sort(theta)

843 theta = theta[::multiplier]

844

845 r = innerRadius + (theta - initialTheta) / (2 * math.pi) * (gap) * multiplier

846 z = np.ones(theta.shape) * z

847

848 # Create the points and store their tags:

849 points = [] # point objects

850 pointTags = [] # point tags

851 breakPointObjectsDueToCutPlane = [] # only used if cutPlaneNormal is not None

852 breakPointTagsDueToCutPlane = [] # only used if cutPlaneNormal is not None

853 pointObjectsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None

854 pointTagsWithoutBreakPoints = [] # only used if cutPlaneNormal is not None

855 breakPointObjectsDueToTransition = []

856 breakPointTagsDueToTransition = []

857 coordinateList = []

858

859 for j in range(len(theta)):

860 pointObject = point(r[j], theta[j], z[j], coordinate.cylindrical)

861 [x_c, y_c, z_c] = [r[j], theta[j], z[j]]

862 #print([x_c, y_c, z_c])

863 coordinateList.append([x_c, y_c, z_c])

864 points.append(pointObject)

865 pointTags.append(pointObject.tag)

866 if cutPlaneNormal is not None:

867 if theta[j] in listOfBreakPoints:

868 breakPointObjectsDueToCutPlane.append(pointObject)

869 breakPointTagsDueToCutPlane.append(pointObject.tag)

870 else:

871 pointObjectsWithoutBreakPoints.append(pointObject)

872 pointTagsWithoutBreakPoints.append(pointObject.tag)

873

874 # identify if the point is a break point due to the layer transition:

875 angle1 = initialTheta + (2 * math.pi - transitionNotchAngle) * multiplier

876 angle2 = (

877 initialTheta

878 + ((self.turns % 1) * 2 * math.pi + transitionNotchAngle) * multiplier

879 )

880 if math.isclose(

881 math.fmod(theta[j], 2 * math.pi), angle1, abs_tol=1e-6

882 ) or math.isclose(math.fmod(theta[j], 2 * math.pi), angle2, abs_tol=1e-6):

883 breakPointObjectsDueToTransition.append(pointObject)

884 breakPointTagsDueToTransition.append(pointObject.tag)

885

886 # Plotter

887 # x_coords = [coord[0] for coord in coordinateList]

888 # y_coords = [coord[1] for coord in coordinateList]

889 # z_coords = [coord[2] for coord in coordinateList]

890

891 # print(f'number of divisions {self.Pancake3DMeshWinding.ane}')

892 # print(self.wi.ane)

893 # Creating the 3D plot

894 # fig = plt.figure()

895 # ax = fig.add_subplot(111, projection='3d')

896

897 # Plotting the coordinates

898 # ax.scatter(x_coords, y_coords, z_coords, c='r', marker='o')

899 #

900 # # Setting labels

901 # ax.set_xlabel('X Label')

902 # ax.set_ylabel('Y Label')

903 # ax.set_zlabel('Z Label')

904 #

905 # plt.show()

906

907 # Logic to break the points up into relevant geom coordinates

908 # Brick points structure (for X-Y plane only for now):

909 # [[[x1, y1, z1], [x2, y2, z2], [x3, y3, z3], [x4, y4, z4]], ...]

910 # Theoretically, very easy to extend to 8 points knowing the height of the

911

912 # Defining the coordinate lists to which the points are to be added

913

914 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number

915 winding_1 = []

916 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number

917 winding_2 = []

918 # winding one covers the list of points in the domain of theta [k, pi*k], where k is an integer number

919 winding_3 = []

920 # winding one covers the list of points in the domain of theta [pi*k, 2pi*k], where k is an integer number

921 winding_4 = []

922 #print(theta[10])

923 # heightPancake = self.geo.winding.height

924 # print(heightPancake)

925 for i in range(len(theta)-1): # range is reduced as no brick can be created starting at the last point

926 # Assuming theta is a numpy array and you're looking for the index of a value close to pi

927 value_to_find = theta[i]+np.pi

928 tolerance = 1e-10 # Define a small tolerance

929 # Find indices where the condition is true

930 indices = np.where(np.abs(theta - value_to_find) < tolerance)[0]

931 if len(indices) > 0:

932 windingUpIndex = indices[0] # Take the first index if there are multiple matches

933 try:

934 x_1 = r[i] * np.cos(theta[i])

935 y_1 = r[i] * np.sin(theta[i])

936 z_g = z[i]

937 x_2 = r[i+1] * np.cos(theta[i+1])

938 y_2 = r[i+1] * np.sin(theta[i+1])

939 x_3 = r[windingUpIndex] * np.cos(theta[windingUpIndex])

940 y_3 = r[windingUpIndex] * np.sin(theta[windingUpIndex])

941 x_4 = r[windingUpIndex+1] * np.cos(theta[windingUpIndex+1])

942 y_4 = r[windingUpIndex+1] * np.sin(theta[windingUpIndex+1])

943 addPoints = [[x_1, y_1, z_g], [x_2, y_2, z_g], [x_3, y_3, z_g], [x_4, y_4, z_g]]

944 k = theta[i]//(2*np.pi)

945 if (theta[i] <= np.pi*(k+1)):

946 # print('winding 1 or 3')

947 if (k%2 == 0):

948 # print('winding 1')

949 winding_1.append(addPoints)

950 else:

951 # print('winding 3')

952 winding_3.append(addPoints)

953

954 if (theta[i] >= np.pi*(k+1)):

955 # print('winding 2 or 4')

956 if (k%2 == 0):

957 # print('winding 2')

958 winding_2.append(addPoints)

959 else:

960 # print('winding 4')

961 winding_4.append(addPoints)

962 except IndexError:

963 print('All of the winding conductor points have been found')

964

965 # print(winding_1)

966 # print(winding_2)

967 # print(winding_3)

968 # print(winding_4)

969

970 # Plotter

971 # x_coords = [coord[0] for coord in winding_1]

972 # y_coords = [coord[1] for coord in winding_1]

973 # z_coords = [coord[2] for coord in winding_1]

974 #

975 # # Creating the 3D plot

976 # fig = plt.figure()

977 # ax = fig.add_subplot(111, projection='3d')

978 #

979 # # Plotting the coordinates

980 # ax.scatter(x_coords, y_coords, z_coords, c='r', marker='o')

981 #

982 # # Setting labels

983 # ax.set_xlabel('X Label')

984 # ax.set_ylabel('Y Label')

985 # ax.set_zlabel('Z Label')

986 #

987 # plt.show()

988

989 # if True:

990 # indexPoint = 1

991 # rangeUpdated = 0

992 # dict_cond = {0: {'SHAPE': 'BR8', 'XCENTRE': '0.0', 'YCENTRE': '0.0', 'ZCENTRE': '0.0', 'PHI1': '0.0', 'THETA1': '0.0', 'PSI1': '0.0', 'XCEN2': '0.0', 'YCEN2': '0.0', 'ZCEN2': '0.0', 'THETA2': '0.0', 'PHI2': '0.0', 'PSI2': '0.0', 'XP1': '-0.879570', 'YP1': '-0.002940', 'ZP1': '-1.131209', 'XP2': '-0.879570', 'YP2': '0.002940', 'ZP2': '-1.131209', 'XP3': '-0.881381', 'YP3': '0.002940', 'ZP3': '-1.114205', 'XP4': '-0.881381', 'YP4': '-0.002940', 'ZP4': '-1.114205', 'XP5': '-0.861227', 'YP5': '-0.002972', 'ZP5': '-1.129183', 'XP6': '-0.861208', 'YP6': '0.002908', 'ZP6': '-1.129182', 'XP7': '-0.863294', 'YP7': '0.002912', 'ZP7': '-1.112210', 'XP8': '-0.863313', 'YP8': '-0.002968', 'ZP8': '-1.112211', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}, 1: {'SHAPE': 'BR8', 'XCENTRE': '0.0', 'YCENTRE': '0.0', 'ZCENTRE': '0.0', 'PHI1': '0.0', 'THETA1': '0.0', 'PSI1': '0.0', 'XCEN2': '0.0', 'YCEN2': '0.0', 'ZCEN2': '0.0', 'THETA2': '0.0', 'PHI2': '0.0', 'PSI2': '0.0', 'XP1': '-0.861227', 'YP1': '-0.002972', 'ZP1': '-1.129183', 'XP2': '-0.861208', 'YP2': '0.002908', 'ZP2': '-1.129182', 'XP3': '-0.863294', 'YP3': '0.002912', 'ZP3': '-1.112210', 'XP4': '-0.863313', 'YP4': '-0.002968', 'ZP4': '-1.112211', 'XP5': '-0.842917', 'YP5': '-0.003066', 'ZP5': '-1.126858', 'XP6': '-0.842880', 'YP6': '0.002814', 'ZP6': '-1.126858', 'XP7': '-0.845242', 'YP7': '0.002830', 'ZP7': '-1.109922', 'XP8': '-0.845278', 'YP8': '-0.003050', 'ZP8': '-1.109922', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}, 2: {'SHAPE': 'BR8', 'XCENTRE': '0.0', 'YCENTRE': '0.0', 'ZCENTRE': '0.0', 'PHI1': '0.0', 'THETA1': '0.0', 'PSI1': '0.0', 'XCEN2': '0.0', 'YCEN2': '0.0', 'ZCEN2': '0.0', 'THETA2': '0.0', 'PHI2': '0.0', 'PSI2': '0.0', 'XP1': '-0.842917', 'YP1': '-0.003066', 'ZP1': '-1.126858', 'XP2': '-0.842880', 'YP2': '0.002814', 'ZP2': '-1.126858', 'XP3': '-0.845242', 'YP3': '0.002830', 'ZP3': '-1.109922', 'XP4': '-0.845278', 'YP4': '-0.003050', 'ZP4': '-1.109922', 'XP5': '-0.824646', 'YP5': '-0.003216', 'ZP5': '-1.124235', 'XP6': '-0.824593', 'YP6': '0.002664', 'ZP6': '-1.124239', 'XP7': '-0.827229', 'YP7': '0.002698', 'ZP7': '-1.107343', 'XP8': '-0.827282', 'YP8': '-0.003181', 'ZP8': '-1.107339', 'CURD': '201264967.975494', 'SYMMETRY': '1', 'DRIVELABEL': 'drive 0', 'IRXY': '0', 'IRYZ': '0', 'IRZX': '0', 'TOLERANCE': '1e-6'}}

993 # # print(dict_cond)

994 # for brick in dict_cond:

995 # for pointIndex in range (rangeUpdated, rangeUpdated+7):

996 # dict_cond[brick][f'XP{indexPoint}'] = str(coordinateList[pointIndex][0])

997 # dict_cond[brick][f'YP{indexPoint}'] = str(coordinateList[pointIndex][1])

998 # dict_cond[brick][f'ZP{indexPoint}'] = str(coordinateList[pointIndex][2])

999 # indexPoint+=1

1000 # indexPoint = 1

1001 # rangeUpdated = rangeUpdated + 8

1002

1003 # writing COND.json file

1004 # Define the path for the JSON file, one directory up from the current script

1005 # json_file_path = os.path.join(os.path.dirname(os.getcwd()), "BR8.json")

1006

1007 # Function to print the contents of a JSON file

1008 # def print_json_contents(path):

1009 # try:

1010 # with open(path, 'r') as file:

1011 # data = json.load(file)

1012 # print(json.dumps(data, indent=4))

1013 # except FileNotFoundError:

1014 # print("File not found.")

1015 # except json.JSONDecodeError:

1016 # print("File is empty or contains non-JSON conforming data.")

1017 #

1018 # # Print current contents

1019 # print("Current contents of BR8.json:")

1020 # print_json_contents(json_file_path)

1021

1022 # Overwrite the JSON file

1023 # with open(json_file_path, 'w') as file:

1024 # json.dump(dict_cond, file, indent=4)

1025 #

1026 # print("\nContents of BR8.json after overwriting:")

1027 # print_json_contents(json_file_path)

1028 # writing the .cond file

1029

1030 # p = tuh.Paths('tests/parsers', '')

1031 # FilesAndFolders.prep_folder(p.model_folder)

1032

1033 # Specify the target directory relative to the current working directory

1034 target_dir = os.path.join(os.getcwd(), 'tests', '_outputs', 'parsers')

1035

1036 # Ensure the target directory exists

1037 os.makedirs(target_dir, exist_ok=True)

1038

1039 # Define the output file path

1040 # out_file = os.path.join(target_dir, 'BR8.cond')

1041 # list_of_shapes = ['BR8']

1042 # for shape in list_of_shapes:

1043 # pc = ParserCOND()

1044 # input_dict = dict_cond

1045 # pc.write_cond(input_dict, out_file)

1046 # print('path')

1047 # print(out_file)

1048 # print('hello world')

1049 # =============================================================================

1050 # GENERATING POINTS ENDS ======================================================

1051 # =============================================================================

1052

1053 # =============================================================================

1054 # GENERATING SPLINES STARTS ===================================================

1055 # =============================================================================

1056

1057 # Create the spline with the points:

1058 spline = gmsh.model.occ.addSpline(pointTags)

1059

1060 # Split the spline into sub-curves:

1061 sectionsPerCurve = int(spt / cpt)

1062

1063 # Create a list of point tags that will split the spline:

1064 # Basically, they are the points to divide the spirals into sectionsPerCurve

1065 # turns. However, some other points are also included to support the float

1066 # number of turns. It is best to visually look at the divisions with

1067 # gmsh.fltk.run() to understand why the split points are chosen the way they are

1068 # selected.

1069

1070 if cutPlaneNormal is None:

1071 pointObjectsWithoutBreakPoints = points

1072 pointTagsWithoutBreakPoints = pointTags

1073

1074 splitPointTags = list(

1075 set(pointTagsWithoutBreakPoints[:-1:sectionsPerCurve])

1076 | set(pointTagsWithoutBreakPoints[-spt - 1 :: -spt])

1077 | set(breakPointTagsDueToCutPlane)

1078 | set(breakPointTagsDueToTransition)

1079 )

1080 splitPointTags = sorted(splitPointTags)

1081 # Remove the first element of the list (starting point):

1082 _, *splitPointTags = splitPointTags

1083

1084 # Also create a list of corresponding point objects:

1085 splitPoints = list(

1086 set(pointObjectsWithoutBreakPoints[:-1:sectionsPerCurve])

1087 | set(pointObjectsWithoutBreakPoints[-spt - 1 :: -spt])

1088 | set(breakPointObjectsDueToCutPlane)

1089 | set(breakPointObjectsDueToTransition)

1090 )

1091 splitPoints = sorted(splitPoints, key=lambda x: x.tag)

1092 # Remove the first element of the list (starting point):

1093 _, *splitPoints = splitPoints

1094

1095 # Split the spline:

1096 dims = [0] * len(splitPointTags)

1097 _, splines = gmsh.model.occ.fragment(

1098 [(1, spline)],

1099 list(zip(dims, splitPointTags)),

1100 removeObject=True,

1101 removeTool=True,

1102 )

1103 splines = splines[0]

1104 self.splineTags = [j for _, j in splines]

1105

1106 # Note the turn number of each spline. This will be used in getSplineTag and

1107 # getSplineTags methods.

1108 self.splineTurns = []

1109 for i in range(len(self.splineTags)):

1110 if i == 0:

1111 startPoint = points[0]

1112 endPoint = splitPoints[0]

1113 elif i == len(self.splineTags) - 1:

1114 startPoint = splitPoints[-1]

1115 endPoint = points[-1]

1116 else:

1117 startPoint = splitPoints[i - 1]

1118 endPoint = splitPoints[i]

1119

1120 startTurn = (startPoint.theta - initialTheta) / (2 * math.pi)

1121 startTurn = round(startTurn / self.turnRes) * self.turnRes

1122 endTurn = (endPoint.theta - initialTheta) / (2 * math.pi)

1123 endTurn = round(endTurn / self.turnRes) * self.turnRes

1124

1125 if direction is direction.ccw:

1126 self.splineTurns.append((startTurn, endTurn))

1127 else:

1128 self.splineTurns.append((-startTurn, -endTurn))

1129

1130 # Check if splineTurn tuples starts with the small turn number:

1131 for i in range(len(self.splineTurns)):

1132 self.splineTurns[i] = sorted(self.splineTurns[i])

1133

1134 # =============================================================================

1135 # GENERATING SPLINES ENDS =====================================================

1136 # =============================================================================

1137

1138 # Find start and end points of the spiral curve:

1139 gmsh.model.occ.synchronize() # synchronize the model to make getBoundary work

1140 self.startPointTag = gmsh.model.getBoundary([(1, self.getSplineTag(0))])[1][1]

1141 self.endPointTag = gmsh.model.getBoundary(

1142 [(1, self.getSplineTag(self.turns, endPoint=True))]

1143 )[1][1]

1144

1145 def getSplineTag(self, turn, endPoint=False):

1146 """

1147 Returns the spline tag at a specific turn. It returns the spline tag of the

1148 section that is on the turn except its end point.

1149

1150 :param turn: turn number (it can be a float)

1151 :type turn: float

1152 :param endPoint: if True, return the spline tag of the section that is on the

1153 turn including its end point but not its start point (default: False)

1154 :type endPoint: bool, optional

1155 :return: spline tag

1156 """

1157 if endPoint:

1158 for i, r in enumerate(self.splineTurns):

1159 if r[0] + (self.turnRes / 2) < turn <= r[1] + (self.turnRes / 2):

1160 return self.splineTags[i]

1161 else:

1162 for i, r in enumerate(self.splineTurns):

1163 if r[0] - (self.turnRes / 2) <= turn < r[1] - (self.turnRes / 2):

1164 return self.splineTags[i]

1165

1166 def getPointTag(self, turn, endPoint=False):

1167 """

1168 Returns the point object at a specific turn.

1169

1170 :param turn: turn number (it can be a float)

1171 :type turn: float

1172 :return: point object

1173 :rtype: point

1174 """

1175 if turn < 0 or turn > self.turns:

1176 raise ValueError("Turn number is out of range!")

1177

1178 if turn == 0:

1179 return self.startPointTag

1180 elif turn == self.turns:

1181 return self.endPointTag

1182 else:

1183 curveTag = self.getSplineTag(turn, endPoint=endPoint)

1184 if endPoint:

1185 points = gmsh.model.getBoundary([(1, curveTag)])

1186 return points[1][1]

1187 else:

1188 points = gmsh.model.getBoundary([(1, curveTag)])

1189 return points[0][1]

1190

1191 def getSplineTags(self, turnStart=None, turnEnd=None):

1192 """

1193 Get the spline tags from a specific turn to another specific turn. If turnStart

1194 and turnEnd are not specified, it returns all the spline tags.

1195

1196 :param turnStart: start turn number (it can be a float) (default: None)

1197 :type turnStart: float, optional

1198 :param turnEnd: end turn number (it can be a float) (default: None)

1199 :type turnEnd: float, optional

1200 :return: spline tags

1201 :rtype: list[int]

1202 """

1203 if turnStart is None and turnEnd is None:

1204 return self.splineTags

1205 elif turnStart is None or turnEnd is None:

1206 raise ValueError(

1207 "turnStart and turnEnd must be both specified or both not specified."

1208 " You specified only one of them."

1209 )

1210 else:

1211 start = self.splineTags.index(self.getSplineTag(turnStart, False))

1212 end = self.splineTags.index(self.getSplineTag(turnEnd, True)) + 1

1213 return self.splineTags[start:end]

1214

1215

1216class spiralSurface:

1217 """

1218 This is a class to create a spiral surface parallel to the XY plane in GMSH. If

1219 thinShellApproximation is set to False, it creates two spiral surfaces parallel to

1220 the XY plane, and their inner and outer curve loops in GMSH. One of the surfaces is

1221 the main surface specified, which is the winding surface, and the other is the

1222 contact layer surface (the gap between the winding surface). If thinShellApproximation

1223 is set to True, it creates only one spiral surface that touches each other

1224 (conformal).

1225

1226 Note that surfaces are subdivided depending on the spiral curve divisions because

1227 this is required for the thin-shell approximation. Otherwise, when

1228 thinShellApproximation is set to True, it would be a disc rather than a spiral since

1229 it touches each other. However, this can be avoided by dividing the surfaces into

1230 small parts and making them conformal. Dividing the surfaces is not necessary when

1231 thinShellApproximation is set to false, but in order to use the same logic with TSA,

1232 it is divided anyway.

1233

1234 :param innerRadius: inner radius

1235 :type innerRadius: float

1236 :param thickness: thickness

1237 :type thickness: float

1238 :param contactLayerThickness: contact layer thickness

1239 :type contactLayerThickness: float

1240 :param turns: number of turns

1241 :type turns: float

1242 :param z: z coordinate

1243 :type z: float

1244 :param initialTheta: initial theta angle in radians

1245 :type initialTheta: float

1246 :param spiralDirection: direction of the spiral (default: direction.ccw)

1247 :type spiralDirection: direction, optional

1248 :param thinShellApproximation: if True, the thin shell approximation is used

1249 (default: False)

1250 :type thinShellApproximation: bool, optional

1251 :param cutPlaneNormal: normal vector of the plane that will cut the spiral surface

1252 (default: None)

1253 :type cutPlaneNormal: tuple[float, float, float], optional

1254 """

1255

1256 def __init__(

1257 self,

1258 innerRadius,

1259 thickness,

1260 contactLayerThickness,

1261 turns,

1262 z,

1263 initialTheta,

1264 transitionNotchAngle,

1265 spiralDirection=direction.ccw,

1266 thinShellApproximation=False,

1267 cutPlaneNormal=None,

1268 ) -> None:

1269 r_i = innerRadius

1270 t = thickness

1271 theta_i = initialTheta

1272 self.theta_i = theta_i

1273

1274 self.surfaceTags = []

1275 self.contactLayerSurfaceTags = []

1276

1277 self.direction = spiralDirection

1278 self.tsa = thinShellApproximation

1279 self.transitionNotchAngle = transitionNotchAngle

1280 # =============================================================================

1281 # GENERATING SPIRAL CURVES STARTS =============================================

1282 # =============================================================================

1283 if thinShellApproximation:

1284 # Create only one spiral curve because the thin shell approximation is used:

1285 # Winding thickness is increased slightly to ensure that the outer radius

1286 # would be the same without the thin shell approximation.

1287 turns = (

1288 turns + 1

1289 ) # for TSA, spiral has (n+1) turns but spiral surface has n

1290 spiral = spiralCurve(

1291 r_i,

1292 t + contactLayerThickness * (turns - 1) / (turns),

1293 turns,

1294 z,

1295 theta_i,

1296 transitionNotchAngle,

1297 spiralDirection,

1298 cutPlaneNormal=cutPlaneNormal,

1299 )

1300

1301 # These are created to be able to use the same code with TSA and without TSA:

1302 innerSpiral = spiral

1303 outerSpiral = spiral

1304 else:

1305 # Create two spiral curves because the thin shell approximation is not used:

1306 innerSpiral = spiralCurve(

1307 r_i - contactLayerThickness,

1308 t + contactLayerThickness,

1309 turns + 1,

1310 z,

1311 theta_i,

1312 transitionNotchAngle,

1313 spiralDirection,

1314 cutPlaneNormal=cutPlaneNormal,

1315 )

1316 outerSpiral = spiralCurve(

1317 r_i,

1318 t + contactLayerThickness,

1319 turns + 1,

1320 z,

1321 theta_i,

1322 transitionNotchAngle,

1323 spiralDirection,

1324 cutPlaneNormal=cutPlaneNormal,

1325 )

1326

1327 self.innerSpiral = innerSpiral

1328 self.outerSpiral = outerSpiral

1329 self.turns = turns

1330 # =============================================================================

1331 # GENERATING SPIRAL CURVES ENDS ===============================================

1332 # =============================================================================

1333

1334 # =============================================================================

1335 # GENERATING SURFACES STARTS ==================================================

1336 # =============================================================================

1337 endLines = []

1338 endInsLines = []

1339

1340 # This is used to check if all the contact layers are finished:

1341 allContactLayersAreFinished = False

1342

1343 # Itterate over the spline tags:

1344 for i in range(len(innerSpiral.splineTags)):

1345 if thinShellApproximation:

1346 # The current spline will be the inner spline:

1347 innerSplineTag = spiral.splineTags[i]

1348

1349 # Find the spline tag of the outer spline by finding the spline tag of

1350 # the spline that is exactly on the next turn:

1351

1352 # Note the turn number of the current spline's start point:

1353 startTurn = spiral.splineTurns[i][0]

1354

1355 if startTurn + 1 + 1e-4 > turns:

1356 # If the current spline is on the outer surface, break the loop,

1357 # because the whole surface is finished:

1358 break

1359

1360 # Find the outer spline tag:

1361 isItBroken = True

1362 for j, turnTuple in enumerate(spiral.splineTurns):

1363 # Equality can not be checked with == because of the floating point

1364 # errors:

1365 if abs(turnTuple[0] - 1 - startTurn) < 1e-4:

1366 outerSplineTag = spiral.splineTags[j]

1367 isItBroken = False

1368 break

1369

1370 if isItBroken:

1371 raise RuntimeError(

1372 "Something went wrong while creating the spiral surface. Outer"

1373 f" spline tag of {innerSplineTag} could not be found for TSA."

1374 )

1375

1376 else:

1377 # Store the tags of the current splines:

1378 innerSplineTag = innerSpiral.splineTags[i]

1379 outerSplineTag = outerSpiral.splineTags[i]

1380

1381 # The current outer spline will be the inner spline of the

1382 # contact layer:

1383 innerInsSplineTag = outerSpiral.splineTags[i]

1384

1385 # Find the spline tag of the contact layer's outer spline by finding the

1386 # spline tag of the spline that is exactly on the next turn:

1387

1388 # Note the turn number of the current spline's start point:

1389 startTurn = outerSpiral.splineTurns[i][0]

1390

1391 if startTurn + 1 + 1e-4 > turns + 1:

1392 # If the current spline is on the outer surface, note that all the

1393 # contact layers are finished:

1394 allContactLayersAreFinished = True

1395

1396 # Find the contact layer's outer spline tag:

1397 for j, turnTuple in enumerate(innerSpiral.splineTurns):

1398 if math.isclose(turnTuple[0], 1 + startTurn, abs_tol=1e-6):

1399 outerInsSplineTag = innerSpiral.splineTags[j]

1400 break

1401

1402 # Create the lines to connect the two splines so that a surface can be

1403 # created:

1404

1405 # Create start line:

1406 if i == 0:

1407 # If it is the first spline, start line should be created.

1408

1409 # Create points:

1410 isStartPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]

1411 if thinShellApproximation:

1412 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]

1413 else:

1414 osStartPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]

1415

1416 # Create the line:

1417 startLine = gmsh.model.occ.addLine(osStartPoint, isStartPoint)

1418 firstStartLine = startLine

1419

1420 # Create lines for the contact layer if the thin shell approximation is not

1421 # used:

1422 if not thinShellApproximation and not allContactLayersAreFinished:

1423 isInsStartPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[

1424 1

1425 ][1]

1426 osInsStartPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[

1427 0

1428 ][1]

1429

1430 # Create the line:

1431 startInsLine = gmsh.model.occ.addLine(

1432 osInsStartPoint, isInsStartPoint

1433 )

1434 firstInsStartLine = startInsLine

1435

1436 else:

1437 # If it is not the first spline, the start line is the end line of the

1438 # previous surface. This guarantees that the surfaces are connected

1439 # (conformality).

1440 startLine = endLines[i - 1]

1441

1442 # Do the same for the contact layer if the thin shell approximation is not

1443 # used:

1444 if not thinShellApproximation and not allContactLayersAreFinished:

1445 startInsLine = endInsLines[i - 1]

1446

1447 # Create end line:

1448

1449 # Create points:

1450 # The ifs are used because getBoundary is not consistent somehow.

1451 if i == 0:

1452 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[0][1]

1453 else:

1454 isEndPoint = gmsh.model.getBoundary([(1, innerSplineTag)])[1][1]

1455

1456 if (not i == 0) or thinShellApproximation:

1457 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[1][1]

1458 else:

1459 osEndPoint = gmsh.model.getBoundary([(1, outerSplineTag)])[0][1]

1460

1461 # Create the line:

1462 endLine = gmsh.model.occ.addLine(isEndPoint, osEndPoint)

1463 endLines.append(endLine)

1464

1465 # Create lines for the contact layer if the thin shell approximation is not

1466 # used:

1467 if not thinShellApproximation and not allContactLayersAreFinished:

1468 if i == 0:

1469 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[0][

1470 1

1471 ]

1472 else:

1473 isInsEndPoint = gmsh.model.getBoundary([(1, innerInsSplineTag)])[1][

1474 1

1475 ]

1476

1477 osInsEndPoint = gmsh.model.getBoundary([(1, outerInsSplineTag)])[1][1]

1478

1479 # Create the line:

1480 endInsLine = gmsh.model.occ.addLine(isInsEndPoint, osInsEndPoint)

1481 endInsLines.append(endInsLine)

1482

1483 # Create the surface:

1484 curveLoop = gmsh.model.occ.addCurveLoop(

1485 [startLine, innerSplineTag, endLine, outerSplineTag]

1486 )

1487 self.surfaceTags.append(gmsh.model.occ.addPlaneSurface([curveLoop]))

1488

1489 # Create the surface for the contact layer if the thin shell approximation is

1490 # not used:

1491 if not thinShellApproximation and not allContactLayersAreFinished:

1492 curveLoop = gmsh.model.occ.addCurveLoop(

1493 [startInsLine, innerInsSplineTag, endInsLine, outerInsSplineTag]

1494 )

1495 self.contactLayerSurfaceTags.append(

1496 gmsh.model.occ.addPlaneSurface([curveLoop])

1497 )

1498

1499 # =============================================================================

1500 # GENERATING SURFACES ENDS ====================================================

1501 # =============================================================================

1502

1503 # =============================================================================

1504 # GENERATING CURVE LOOPS STARTS ===============================================

1505 # =============================================================================

1506

1507 # Create the inner and outer curve loops (for both TSA == True and TSA == False):

1508

1509 # VERY IMPORTANT NOTES ABOUT THE DIRECTION OF THE CURVE LOOPS

1510 # 1- GMSH doesn't like duplicates. Or the user doesn't like duplicates if they

1511 # want conformality. Actually, it's a positive thing about debugging because

1512 # you can always use `Geometry.remove_all_duplicates()` to see if there are

1513 # any duplicates. If there are, the problem is found. Solve it.

1514 # 2- The problem arises when one uses surface loops or curve loops. Because even

1515 # if you think there are no duplicates, GMSH/OCC might create some during

1516 # addCurveLoops and addSurfaceLoops operations. Even though

1517 # `geometry.remove_all_duplicates()` tells that there are duplicates, the

1518 # user doesn't suspect about addCurveLoop and addSurfaceLoop at first,

1519 # because it doesn't make sense.

1520 # 3- How you put the curves in the curve loops is very important! The same curve

1521 # loop with the same lines might cause problems if the user puts them in a

1522 # different order. For example, to create a plane surface with two curve

1523 # loops, the direction of the curve loops should be the same. That's why the

1524 # code has both innerCurveLoopTag and innerOppositeCurveLoopTag (the same

1525 # thing for the outer curve loop).

1526

1527 # create the transition layer (notch):

1528 # Inner curve loop:

1529 notchStartPoint = innerSpiral.getPointTag(

1530 1 - transitionNotchAngle / (2 * math.pi)

1531 )

1532 notchLeftPoint = innerSpiral.getPointTag(0)

1533 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)

1534 notchRightLine = innerSpiral.getSplineTag(

1535 1 - transitionNotchAngle / (2 * math.pi)

1536 )

1537

1538 if thinShellApproximation:

1539 innerStartCurves = [firstStartLine]

1540 else:

1541 innerStartCurves = [firstInsStartLine, firstStartLine]

1542

1543 if thinShellApproximation:

1544

1545 notchCurveLoop = gmsh.model.occ.addCurveLoop(

1546 [notchLeftLine, notchRightLine] + innerStartCurves

1547 )

1548 self.innerNotchSurfaceTags = [

1549 gmsh.model.occ.addPlaneSurface([notchCurveLoop])

1550 ]

1551 else:

1552 notchMiddlePoint = outerSpiral.getPointTag(0)

1553 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)

1554

1555 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(

1556 [notchLeftLine, notchMiddleLine, firstStartLine]

1557 )

1558 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(

1559 [notchMiddleLine, notchRightLine, firstInsStartLine]

1560 )

1561 self.innerNotchSurfaceTags = [

1562 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),

1563 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),

1564 ]

1565

1566 lines = innerSpiral.getSplineTags(

1567 0, 1 - transitionNotchAngle / (2 * math.pi)

1568 ) # The first turn of the spline

1569 innerCurves = lines + [notchLeftLine]

1570 self.innerNotchLeftLine = notchLeftLine

1571

1572 self.innerStartCurves = innerStartCurves

1573

1574 self.innerCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves)

1575 self.innerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(innerCurves[::-1])

1576

1577 # Outer curve loop:

1578 # The last turn of the spline:

1579 if thinShellApproximation:

1580 notchStartPoint = innerSpiral.getPointTag(

1581 self.turns + transitionNotchAngle / (2 * math.pi) - 1

1582 )

1583 notchLeftPoint = innerSpiral.getPointTag(self.turns)

1584 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)

1585 notchRightLine = innerSpiral.getSplineTag(

1586 self.turns - 1 + transitionNotchAngle / (2 * math.pi), endPoint=True

1587 )

1588 else:

1589 notchStartPoint = outerSpiral.getPointTag(

1590 self.turns + transitionNotchAngle / (2 * math.pi)

1591 )

1592 notchMiddlePoint = innerSpiral.getPointTag(self.turns + 1)

1593 notchLeftPoint = outerSpiral.getPointTag(self.turns + 1)

1594 notchLeftLine = gmsh.model.occ.addLine(notchStartPoint, notchLeftPoint)

1595 notchMiddleLine = gmsh.model.occ.addLine(notchStartPoint, notchMiddlePoint)

1596 notchRightLine = outerSpiral.getSplineTag(

1597 self.turns + transitionNotchAngle / (2 * math.pi), self.turns

1598 )

1599 if thinShellApproximation:

1600 lines = outerSpiral.getSplineTags(turns - 1, turns)

1601 else:

1602 lines = outerSpiral.getSplineTags(turns, turns + 1)

1603

1604 if thinShellApproximation:

1605 outerEndCurves = [endLines[-1]]

1606 else:

1607 outerEndCurves = [endInsLines[-1], endLines[-1]]

1608

1609 if thinShellApproximation:

1610 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(

1611 [notchLeftLine, notchRightLine, endLines[-1]]

1612 )

1613 self.outerNotchSurfaceTags = [

1614 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),

1615 ]

1616 else:

1617 notchCurveLoop1 = gmsh.model.occ.addCurveLoop(

1618 [notchLeftLine, notchMiddleLine, endLines[-1]]

1619 )

1620 notchCurveLoop2 = gmsh.model.occ.addCurveLoop(

1621 [notchMiddleLine, notchRightLine, endInsLines[-1]]

1622 )

1623 self.outerNotchSurfaceTags = [

1624 gmsh.model.occ.addPlaneSurface([notchCurveLoop1]),

1625 gmsh.model.occ.addPlaneSurface([notchCurveLoop2]),

1626 ]

1627

1628 if thinShellApproximation:

1629 lines = innerSpiral.getSplineTags(

1630 self.turns - 1 + transitionNotchAngle / (2 * math.pi), self.turns

1631 ) # The first turn of the spline

1632 else:

1633 lines = outerSpiral.getSplineTags(

1634 self.turns + transitionNotchAngle / (2 * math.pi), self.turns + 1

1635 )

1636 outerCurves = lines + [notchLeftLine]

1637 self.outerNotchLeftLine = notchLeftLine

1638

1639 self.outerEndCurves = outerEndCurves

1640

1641 self.outerCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves)

1642 self.outerOppositeCurveLoopTag = gmsh.model.occ.addCurveLoop(outerCurves[::-1])

1643 # =============================================================================

1644 # GENERATING CURVE LOOPS ENDS =================================================

1645 # =============================================================================

1646

1647 if not thinShellApproximation:

1648 surfaceTags = self.surfaceTags

1649 self.surfaceTags = self.contactLayerSurfaceTags

1650 self.contactLayerSurfaceTags = surfaceTags

1651

1652 def getInnerRightPointTag(self):

1653 """

1654 Returns the point tag of the inner left point.

1655

1656 :return: point tag

1657 :rtype: int

1658 """

1659 return self.innerSpiral.getPointTag(0)

1660

1661 def getInnerUpperPointTag(self):

1662 """

1663 Returns the point tag of the inner right point.

1664

1665 :return: point tag

1666 :rtype: int

1667 """

1668 if self.direction is direction.ccw:

1669 return self.innerSpiral.getPointTag(0.25)

1670 else:

1671 return self.innerSpiral.getPointTag(0.75)

1672

1673 def getInnerLeftPointTag(self):

1674 """

1675 Returns the point tag of the inner upper point.

1676

1677 :return: point tag

1678 :rtype: int

1679 """

1680 return self.innerSpiral.getPointTag(0.5)

1681

1682 def getInnerLowerPointTag(self):

1683 """

1684 Returns the point tag of the inner lower point.

1685

1686 :return: point tag

1687 :rtype: int

1688 """

1689 if self.direction is direction.ccw:

1690 return self.innerSpiral.getPointTag(0.75)

1691 else:

1692 return self.innerSpiral.getPointTag(0.25)

1693

1694 def getOuterRightPointTag(self):

1695 """

1696 Returns the point tag of the outer left point.

1697

1698 :return: point tag

1699 :rtype: int

1700 """

1701 if self.tsa:

1702 turns = self.turns

1703 else:

1704 turns = self.turns + 1

1705 return self.outerSpiral.getPointTag(turns, endPoint=False)

1706

1707 def getOuterLowerPointTag(self):

1708 """

1709 Returns the point tag of the outer right point.

1710

1711 :return: point tag

1712 :rtype: int

1713 """

1714 if self.tsa:

1715 turns = self.turns

1716 else:

1717 turns = self.turns + 1

1718 if self.direction is direction.ccw:

1719 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)

1720 else:

1721 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)

1722

1723 def getOuterLeftPointTag(self):

1724 """

1725 Returns the point tag of the outer upper point.

1726

1727 :return: point tag

1728 :rtype: int

1729 """

1730 if self.tsa:

1731 turns = self.turns

1732 else:

1733 turns = self.turns + 1

1734 return self.outerSpiral.getPointTag(turns - 0.5, endPoint=False)

1735

1736 def getOuterUpperPointTag(self):

1737 """

1738 Returns the point tag of the outer lower point.

1739

1740 :return: point tag

1741 :rtype: int

1742 """

1743 if self.tsa:

1744 turns = self.turns

1745 else:

1746 turns = self.turns + 1

1747 if self.direction is direction.ccw:

1748 return self.outerSpiral.getPointTag(turns - 0.75, endPoint=False)

1749 else:

1750 return self.outerSpiral.getPointTag(turns - 0.25, endPoint=False)

1751

1752 def getInnerUpperRightCurves(self):

1753 """

1754 Returns the curve tags of the upper right curves.

1755

1756 :return: curve tags

1757 :rtype: list[int]

1758 """

1759 if self.direction is direction.ccw:

1760 curves = self.innerSpiral.getSplineTags(0, 0.25)

1761 else:

1762 lines = self.innerSpiral.getSplineTags(

1763 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)

1764 ) # The first turn of the spline

1765 lines = lines + [self.innerNotchLeftLine]

1766

1767 return lines

1768

1769 return curves

1770

1771 def getInnerUpperLeftCurves(self):

1772 """

1773 Returns the curve tags of the upper left curves.

1774

1775 :return: curve tags

1776 :rtype: list[int]

1777 """

1778 if self.direction is direction.ccw:

1779 curves = self.innerSpiral.getSplineTags(0.25, 0.5)

1780 else:

1781 curves = self.innerSpiral.getSplineTags(0.5, 0.75)

1782

1783 return curves

1784

1785 def getInnerLowerLeftCurves(self):

1786 """

1787 Returns the curve tags of the lower left curves.

1788

1789 :return: curve tags

1790 :rtype: list[int]

1791 """

1792 if self.direction is direction.ccw:

1793 curves = self.innerSpiral.getSplineTags(0.5, 0.75)

1794 else:

1795 curves = self.innerSpiral.getSplineTags(0.25, 0.5)

1796

1797 return curves

1798

1799 def getInnerLowerRightCurves(self):

1800 """

1801 Returns the curve tags of the lower right curves.

1802

1803 :return: curve tags

1804 :rtype: list[int]

1805 """

1806 if self.direction is direction.ccw:

1807 lines = self.innerSpiral.getSplineTags(

1808 0.75, 1 - self.transitionNotchAngle / (2 * math.pi)

1809 ) # The first turn of the spline

1810 lines = lines + [self.innerNotchLeftLine]

1811

1812 return lines

1813 else:

1814 curves = self.innerSpiral.getSplineTags(0, 0.25)

1815

1816 return curves

1817

1818 def getOuterUpperRightCurves(self):

1819 """

1820 Returns the curve tags of the upper right curves.

1821

1822 :return: curve tags

1823 :rtype: list[int]

1824 """

1825 if self.tsa:

1826 turns = self.turns

1827 else:

1828 turns = self.turns + 1

1829

1830 if self.direction is direction.ccw:

1831 if self.tsa:

1832 lines = self.innerSpiral.getSplineTags(

1833 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),

1834 self.turns - 0.75,

1835 ) # The first turn of the spline

1836 else:

1837 lines = self.outerSpiral.getSplineTags(

1838 self.turns + self.transitionNotchAngle / (2 * math.pi),

1839 self.turns + 1 - 0.75,

1840 )

1841 lines = lines + [self.outerNotchLeftLine]

1842

1843 return lines

1844 else:

1845 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)

1846

1847 return curves

1848

1849 def getOuterUpperLeftCurves(self):

1850 """

1851 Returns the curve tags of the lower right curves.

1852

1853 :return: curve tags

1854 :rtype: list[int]

1855 """

1856 if self.tsa:

1857 turns = self.turns

1858 else:

1859 turns = self.turns + 1

1860 if self.direction is direction.ccw:

1861 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)

1862 else:

1863 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)

1864

1865 return curves

1866

1867 def getOuterLowerLeftCurves(self):

1868 """

1869 Returns the curve tags of the lower left curves.

1870

1871 :return: curve tags

1872 :rtype: list[int]

1873 """

1874 if self.tsa:

1875 turns = self.turns

1876 else:

1877 turns = self.turns + 1

1878 if self.direction is direction.ccw:

1879 curves = self.outerSpiral.getSplineTags(turns - 0.5, turns - 0.25)

1880 else:

1881 curves = self.outerSpiral.getSplineTags(turns - 0.75, turns - 0.5)

1882

1883 return curves

1884

1885 def getOuterLowerRightCurves(self):

1886 """

1887 Returns the curve tags of the upper left curves.

1888

1889 :return: curve tags

1890 :rtype: list[int]

1891 """

1892 if self.tsa:

1893 turns = self.turns

1894 else:

1895 turns = self.turns + 1

1896 if self.direction is direction.ccw:

1897 curves = self.outerSpiral.getSplineTags(turns - 0.25, turns)

1898 else:

1899 if self.tsa:

1900 lines = self.innerSpiral.getSplineTags(

1901 self.turns - 1 + self.transitionNotchAngle / (2 * math.pi),

1902 self.turns - 0.75,

1903 ) # The first turn of the spline

1904 else:

1905 lines = self.outerSpiral.getSplineTags(

1906 self.turns + self.transitionNotchAngle / (2 * math.pi),

1907 self.turns + 1 - 0.75,

1908 )

1909 lines = lines + [self.outerNotchLeftLine]

1910

1911 return lines

1912

1913 return curves

1914

1915 def getInnerStartCurves(self):

1916 """

1917 Returns the curve tags of the start curves.

1918

1919 :return: curve tags

1920 :rtype: list[int]

1921 """

1922 return self.innerStartCurves

1923

1924 def getOuterEndCurves(self):

1925 """

1926 Returns the curve tags of the end curves.

1927

1928 :return: curve tags

1929 :rtype: list[int]

1930 """

1931 return self.outerEndCurves

1932

1933

1934class circleWithFourCurves:

1935 def __init__(

1936 self,

1937 x,

1938 y,

1939 z,

1940 radius,

1941 upperRightTag=None,

1942 upperLeftTag=None,

1943 lowerLeftTag=None,

1944 lowerRightTag=None,

1945 leftPointTag=None,

1946 rightPointTag=None,

1947 upperPointTag=None,

1948 lowerPointTag=None,

1949 ):

1950 self.x = x

1951 self.y = y

1952 self.z = z

1953 self.r = radius

1954 if upperRightTag is None:

1955 dummyCircle = gmsh.model.occ.addCircle(self.x, self.y, self.z, self.r)

1956 self.leftPointTag = gmsh.model.occ.addPoint(self.x - self.r, self.y, self.z)

1957 self.rightPointTag = gmsh.model.occ.addPoint(

1958 self.x + self.r, self.y, self.z

1959 )

1960 self.upperPointTag = gmsh.model.occ.addPoint(

1961 self.x, self.y + self.r, self.z

1962 )

1963 self.lowerPointTag = gmsh.model.occ.addPoint(

1964 self.x, self.y - self.r, self.z

1965 )

1966

1967 fragmentResults = gmsh.model.occ.fragment(

1968 [(1, dummyCircle)],

1969 [

1970 (0, self.leftPointTag),

1971 (0, self.rightPointTag),

1972 (0, self.upperPointTag),

1973 (0, self.lowerPointTag),

1974 ],

1975 )[0]

1976 linesDimTags = [dimTag for dimTag in fragmentResults if dimTag[0] == 1]

1977

1978 self.upperRightTag = linesDimTags[0][1]

1979 self.upperLeftTag = linesDimTags[1][1]

1980 self.lowerLeftTag = linesDimTags[2][1]

1981 self.lowerRightTag = linesDimTags[3][1]

1982 else:

1983 self.upperRightTag = upperRightTag

1984 self.upperLeftTag = upperLeftTag

1985 self.lowerLeftTag = lowerLeftTag

1986 self.lowerRightTag = lowerRightTag

1987

1988 self.leftPointTag = leftPointTag

1989 self.rightPointTag = rightPointTag

1990 self.upperPointTag = upperPointTag

1991 self.lowerPointTag = lowerPointTag

1992

1993

1994class outerAirSurface:

1995 def __init__(

1996 self,

1997 outerRadius,

1998 innerRadius,

1999 type="cylinder",

2000 divideIntoFourParts=False,

2001 divideTerminalPartIntoFourParts=False,

2002 ):

2003 self.surfaceTags = []

2004

2005 self.divideIntoFourParts = divideIntoFourParts

2006 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts

2007

2008 # for cylinder:

2009 self.shellTags = []

2010

2011 # for cuboid:

2012 self.shellTagsPart1 = []

2013 self.shellTagsPart2 = []

2014 self.shellTagsPart3 = []

2015 self.shellTagsPart4 = []

2016

2017 self.type = type

2018 self.outerRadius = outerRadius

2019 self.innerRadius = innerRadius

2020

2021 def createFromScratch(self, z, shellTransformation=False, shellRadius=None):

2022 self.z = z

2023

2024 if self.divideIntoFourParts:

2025 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)

2026 else:

2027 if self.divideTerminalPartIntoFourParts:

2028 self.innerCircle = circleWithFourCurves(0, 0, z, self.innerRadius)

2029 innerCL = gmsh.model.occ.addCurveLoop(

2030 [

2031 self.innerCircle.upperRightTag,

2032 self.innerCircle.upperLeftTag,

2033 self.innerCircle.lowerLeftTag,

2034 self.innerCircle.lowerRightTag,

2035 ]

2036 )

2037 else:

2038 innerCL = gmsh.model.occ.addCircle(0, 0, z, self.innerRadius)

2039 innerCL = gmsh.model.occ.addCurveLoop([innerCL])

2040

2041 if self.type == "cylinder" and self.divideIntoFourParts:

2042 outerCircle = circleWithFourCurves(0, 0, z, self.outerRadius)

2043

2044 leftLineTag = gmsh.model.occ.addLine(

2045 outerCircle.leftPointTag, self.innerCircle.leftPointTag

2046 )

2047 rightLineTag = gmsh.model.occ.addLine(

2048 outerCircle.rightPointTag, self.innerCircle.rightPointTag

2049 )

2050 upperLineTag = gmsh.model.occ.addLine(

2051 outerCircle.upperPointTag, self.innerCircle.upperPointTag

2052 )

2053 lowerLineTag = gmsh.model.occ.addLine(

2054 outerCircle.lowerPointTag, self.innerCircle.lowerPointTag

2055 )

2056

2057 # Create surfaces:

2058 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2059 [

2060 outerCircle.upperRightTag,

2061 rightLineTag,

2062 self.innerCircle.upperRightTag,

2063 -upperLineTag,

2064 ]

2065 )

2066 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2067 self.surfaceTags.append(self.upperRightTag)

2068

2069 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2070 [

2071 outerCircle.upperLeftTag,

2072 leftLineTag,

2073 self.innerCircle.upperLeftTag,

2074 -upperLineTag,

2075 ]

2076 )

2077 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2078 self.surfaceTags.append(self.upperLeftTag)

2079

2080 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2081 [

2082 outerCircle.lowerLeftTag,

2083 leftLineTag,

2084 self.innerCircle.lowerLeftTag,

2085 -lowerLineTag,

2086 ]

2087 )

2088 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2089 self.surfaceTags.append(self.lowerLeftTag)

2090

2091 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2092 [

2093 outerCircle.lowerRightTag,

2094 rightLineTag,

2095 self.innerCircle.lowerRightTag,

2096 -lowerLineTag,

2097 ]

2098 )

2099 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2100 self.surfaceTags.append(self.lowerRightTag)

2101

2102 if shellTransformation:

2103 shellOuterCircle = circleWithFourCurves(0, 0, z, shellRadius)

2104 shellLeftLineTag = gmsh.model.occ.addLine(

2105 shellOuterCircle.leftPointTag, outerCircle.leftPointTag

2106 )

2107 shellRightLineTag = gmsh.model.occ.addLine(

2108 shellOuterCircle.rightPointTag, outerCircle.rightPointTag

2109 )

2110 shellUpperLineTag = gmsh.model.occ.addLine(

2111 shellOuterCircle.upperPointTag, outerCircle.upperPointTag

2112 )

2113 shellLowerLineTag = gmsh.model.occ.addLine(

2114 shellOuterCircle.lowerPointTag, outerCircle.lowerPointTag

2115 )

2116

2117 # Create surfaces:

2118 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2119 [

2120 shellOuterCircle.upperRightTag,

2121 shellRightLineTag,

2122 outerCircle.upperRightTag,

2123 -shellUpperLineTag,

2124 ]

2125 )

2126 self.upperRightTag = gmsh.model.occ.addPlaneSurface(

2127 [upperRightCurveLoop]

2128 )

2129 self.shellTags.append(self.upperRightTag)

2130

2131 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2132 [

2133 shellOuterCircle.upperLeftTag,

2134 shellLeftLineTag,

2135 outerCircle.upperLeftTag,

2136 -shellUpperLineTag,

2137 ]

2138 )

2139 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2140 self.shellTags.append(self.upperLeftTag)

2141

2142 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2143 [

2144 shellOuterCircle.lowerLeftTag,

2145 shellLeftLineTag,

2146 outerCircle.lowerLeftTag,

2147 -shellLowerLineTag,

2148 ]

2149 )

2150 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2151 self.shellTags.append(self.lowerLeftTag)

2152

2153 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2154 [

2155 shellOuterCircle.lowerRightTag,

2156 shellRightLineTag,

2157 outerCircle.lowerRightTag,

2158 -shellLowerLineTag,

2159 ]

2160 )

2161 self.lowerRightTag = gmsh.model.occ.addPlaneSurface(

2162 [lowerRightCurveLoop]

2163 )

2164 self.shellTags.append(self.lowerRightTag)

2165

2166 elif self.type == "cylinder" and not self.divideIntoFourParts:

2167 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.outerRadius)

2168 outerCL = gmsh.model.occ.addCurveLoop([outerCL])

2169

2170 if shellTransformation:

2171 shellOuterCL = gmsh.model.occ.addCircle(0, 0, z, shellRadius)

2172 shellOuterCL = gmsh.model.occ.addCurveLoop([shellOuterCL])

2173

2174 shellSurfaceTag = gmsh.model.occ.addPlaneSurface(

2175 [shellOuterCL, outerCL]

2176 )

2177 self.shellTags.append(shellSurfaceTag)

2178

2179 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])

2180 self.surfaceTags.append(surfaceTag)

2181

2182 elif self.type == "cuboid":

2183 # LL: lower left

2184 # LR: lower right

2185 # UR: upper right

2186 # UL: upper left

2187 airLLpointTag = point(-self.outerRadius, -self.outerRadius, z).tag

2188 airLRpointTag = point(self.outerRadius, -self.outerRadius, z).tag

2189 airURpointTag = point(self.outerRadius, self.outerRadius, z).tag

2190 airULpointTag = point(-self.outerRadius, self.outerRadius, z).tag

2191

2192 # LH: lower horizontal

2193 # UH: upper horizontal

2194 # LV: left vertical

2195 # RV: right vertical

2196 airLHlineTag = gmsh.model.occ.addLine(airLLpointTag, airLRpointTag)

2197 airRVLineTag = gmsh.model.occ.addLine(airLRpointTag, airURpointTag)

2198 airUHLineTag = gmsh.model.occ.addLine(airURpointTag, airULpointTag)

2199 airLVLineTag = gmsh.model.occ.addLine(airULpointTag, airLLpointTag)

2200

2201 outerCL = gmsh.model.occ.addCurveLoop(

2202 [airLHlineTag, airRVLineTag, airUHLineTag, airLVLineTag]

2203 )

2204

2205 if self.divideIntoFourParts:

2206 innerCL = gmsh.model.occ.addCurveLoop(

2207 [

2208 self.innerCircle.upperRightTag,

2209 self.innerCircle.lowerRightTag,

2210 self.innerCircle.lowerLeftTag,

2211 self.innerCircle.upperLeftTag,

2212 ]

2213 )

2214

2215 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL, innerCL])

2216 self.surfaceTags.append(surfaceTag)

2217

2218 if shellTransformation:

2219 # LL: lower left

2220 # LR: lower right

2221 # UR: upper right

2222 # UL: upper left

2223 shellLLpointTag = point(

2224 -shellRadius,

2225 -shellRadius,

2226 z,

2227 ).tag

2228 shellLRpointTag = point(

2229 shellRadius,

2230 -shellRadius,

2231 z,

2232 ).tag

2233 shellURpointTag = point(

2234 shellRadius,

2235 shellRadius,

2236 z,

2237 ).tag

2238 shellULpointTag = point(

2239 -shellRadius,

2240 shellRadius,

2241 z,

2242 ).tag

2243

2244 # LH: lower horizontal

2245 # UH: upper horizontal

2246 # LV: left vertical

2247 # RV: right vertical

2248 shellLHlineTag = gmsh.model.occ.addLine(

2249 shellLLpointTag, shellLRpointTag

2250 )

2251 shellRVLineTag = gmsh.model.occ.addLine(

2252 shellLRpointTag, shellURpointTag

2253 )

2254 shellUHLineTag = gmsh.model.occ.addLine(

2255 shellURpointTag, shellULpointTag

2256 )

2257 shellLVLineTag = gmsh.model.occ.addLine(

2258 shellULpointTag, shellLLpointTag

2259 )

2260

2261 shellLowerLeftLineTag = gmsh.model.occ.addLine(

2262 shellLLpointTag, airLLpointTag

2263 )

2264 shellLowerRightLineTag = gmsh.model.occ.addLine(

2265 shellLRpointTag, airLRpointTag

2266 )

2267 shellUpperLeftLineTag = gmsh.model.occ.addLine(

2268 shellULpointTag, airULpointTag

2269 )

2270 shellUpperRightLineTag = gmsh.model.occ.addLine(

2271 shellURpointTag, airURpointTag

2272 )

2273

2274 # Shell lower surface:

2275 shellLowerPSTag = gmsh.model.occ.addCurveLoop(

2276 [

2277 shellLowerLeftLineTag,

2278 airLHlineTag,

2279 shellLowerRightLineTag,

2280 shellLHlineTag,

2281 ]

2282 )

2283 shellLowerPSTag = gmsh.model.occ.addPlaneSurface([shellLowerPSTag])

2284 self.shellTagsPart1.append(shellLowerPSTag)

2285

2286 # Shell right surface:

2287 shellRightPSTag = gmsh.model.occ.addCurveLoop(

2288 [

2289 shellLowerRightLineTag,

2290 airRVLineTag,

2291 shellUpperRightLineTag,

2292 shellRVLineTag,

2293 ]

2294 )

2295 shellRightPSTag = gmsh.model.occ.addPlaneSurface([shellRightPSTag])

2296 self.shellTagsPart2.append(shellRightPSTag)

2297

2298 # Shell upper surface:

2299 shellUpperPSTag = gmsh.model.occ.addCurveLoop(

2300 [

2301 shellUpperLeftLineTag,

2302 airUHLineTag,

2303 shellUpperRightLineTag,

2304 shellUHLineTag,

2305 ]

2306 )

2307 shellUpperPSTag = gmsh.model.occ.addPlaneSurface([shellUpperPSTag])

2308 self.shellTagsPart3.append(shellUpperPSTag)

2309

2310 # Shell left surface:

2311 shellLeftPSTag = gmsh.model.occ.addCurveLoop(

2312 [

2313 shellLowerLeftLineTag,

2314 airLVLineTag,

2315 shellUpperLeftLineTag,

2316 shellLVLineTag,

2317 ]

2318 )

2319 shellLeftPSTag = gmsh.model.occ.addPlaneSurface([shellLeftPSTag])

2320 self.shellTagsPart4.append(shellLeftPSTag)

2321

2322 def setPrecreatedSurfaceTags(

2323 self,

2324 surfaceTags,

2325 cylinderShellTags=None,

2326 cuboidShellTags1=None,

2327 cuboidShellTags2=None,

2328 cuboidShellTags3=None,

2329 cuboidShellTags4=None,

2330 ):

2331 if not isinstance(surfaceTags, list):

2332 raise TypeError("surfaceTags must be a list.")

2333

2334 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2]

2335 self.surfaceTags.extend(surfaceTags)

2336

2337 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:

2338 # Create innerCircle object from the tags:

2339 curves = gmsh.model.getBoundary(

2340 [(2, tag) for tag in surfaceTags], oriented=False

2341 )

2342 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)

2343 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]

2344 innerCurveTags.sort()

2345 upperRightCurve = innerCurveTags[0]

2346 upperLeftCurve = innerCurveTags[1]

2347 lowerLeftCurve = innerCurveTags[2]

2348 lowerRightCurve = innerCurveTags[3]

2349

2350 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)

2351 pointTags = [dimTag[1] for dimTag in points]

2352 pointTags.sort()

2353 upperPointTag = pointTags[0]

2354 leftPointTag = pointTags[1]

2355

2356 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)

2357 pointTags = [dimTag[1] for dimTag in points]

2358 pointTags.sort()

2359 rightPointTag = pointTags[0]

2360 lowerPointTag = pointTags[1]

2361

2362 self.innerCircle = circleWithFourCurves(

2363 0,

2364 0,

2365 self.z,

2366 self.outerRadius,

2367 upperRightTag=upperRightCurve,

2368 upperLeftTag=upperLeftCurve,

2369 lowerLeftTag=lowerLeftCurve,

2370 lowerRightTag=lowerRightCurve,

2371 leftPointTag=leftPointTag,

2372 rightPointTag=rightPointTag,

2373 upperPointTag=upperPointTag,

2374 lowerPointTag=lowerPointTag,

2375 )

2376

2377 if cylinderShellTags is not None:

2378 self.shellTags.extend(cylinderShellTags)

2379

2380 if cuboidShellTags1 is not None:

2381 self.shellTagsPart1.extend(cuboidShellTags1)

2382 self.shellTagsPart2.extend(cuboidShellTags2)

2383 self.shellTagsPart3.extend(cuboidShellTags3)

2384 self.shellTagsPart4.extend(cuboidShellTags4)

2385

2386 def getInnerCL(self):

2387 # checked!

2388 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])

2389 # curves = list(curves)

2390 # curves = list(curves[1])

2391

2392 # innerCL = gmsh.model.occ.addCurveLoop(curves)

2393

2394 gmsh.model.occ.synchronize() # don't delete this line

2395 curves = gmsh.model.getBoundary(

2396 [(2, tag) for tag in self.surfaceTags], oriented=False

2397 )

2398 innerCurveDimTags = findOuterOnes(curves, findInnerOnes=True)

2399 innerCurveTags = [dimTag[1] for dimTag in innerCurveDimTags]

2400

2401 innerCL = gmsh.model.occ.addCurveLoop(innerCurveTags)

2402 return innerCL

2403

2404

2405class outerTerminalSurface:

2406 def __init__(

2407 self,

2408 outerRadius,

2409 tubeThickness,

2410 divideIntoFourParts=False,

2411 ):

2412 self.tubeSurfaceTags = []

2413 self.nontubeSurfaceTags = []

2414

2415 self.divideIntoFourParts = divideIntoFourParts

2416

2417 self.outerRadius = outerRadius

2418 self.tubeThickness = tubeThickness

2419

2420 def createNontubePartWithMiddleCircleAndWinding(

2421 self, middleCircle: circleWithFourCurves, winding: spiralSurface

2422 ):

2423 leftLineTag = gmsh.model.occ.addLine(

2424 middleCircle.leftPointTag, winding.getOuterLeftPointTag()

2425 )

2426 rightLineTag = gmsh.model.occ.addLine(

2427 middleCircle.rightPointTag, winding.getOuterRightPointTag()

2428 )

2429 upperLineTag = gmsh.model.occ.addLine(

2430 middleCircle.upperPointTag, winding.getOuterUpperPointTag()

2431 )

2432 lowerLineTag = gmsh.model.occ.addLine(

2433 middleCircle.lowerPointTag, winding.getOuterLowerPointTag()

2434 )

2435

2436 # Create surfaces for the nontube part:

2437 if winding.direction is direction.ccw:

2438 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2439 winding.getOuterUpperRightCurves()

2440 # + winding.getOuterEndCurves()

2441 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]

2442 )

2443 else:

2444 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2445 winding.getOuterUpperRightCurves()

2446 + [rightLineTag, middleCircle.upperRightTag, upperLineTag]

2447 )

2448 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2449 self.nontubeSurfaceTags.append(self.upperRightTag)

2450

2451 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2452 winding.getOuterUpperLeftCurves()

2453 + [leftLineTag, middleCircle.upperLeftTag, upperLineTag]

2454 )

2455 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2456 self.nontubeSurfaceTags.append(self.upperLeftTag)

2457

2458 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2459 winding.getOuterLowerLeftCurves()

2460 + [leftLineTag, middleCircle.lowerLeftTag, lowerLineTag]

2461 )

2462 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2463 self.nontubeSurfaceTags.append(self.lowerLeftTag)

2464

2465 if winding.direction is direction.ccw:

2466 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2467 winding.getOuterLowerRightCurves()

2468 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]

2469 )

2470 else:

2471 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2472 winding.getOuterLowerRightCurves()

2473 # + winding.getOuterEndCurves()

2474 + [rightLineTag, middleCircle.lowerRightTag, lowerLineTag]

2475 )

2476 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2477 self.nontubeSurfaceTags.append(self.lowerRightTag)

2478

2479 def createWithOuterAirAndWinding(

2480 self, outerAir: outerAirSurface, winding: spiralSurface, pancakeIndex

2481 ):

2482 # Tube part:

2483 z = outerAir.z

2484

2485 if self.divideIntoFourParts:

2486 outerCircle = outerAir.innerCircle

2487 middleCircle = circleWithFourCurves(

2488 0, 0, z, self.outerRadius - self.tubeThickness

2489 )

2490

2491 leftLineTag = gmsh.model.occ.addLine(

2492 outerCircle.leftPointTag, middleCircle.leftPointTag

2493 )

2494 rightLineTag = gmsh.model.occ.addLine(

2495 outerCircle.rightPointTag, middleCircle.rightPointTag

2496 )

2497 upperLineTag = gmsh.model.occ.addLine(

2498 outerCircle.upperPointTag, middleCircle.upperPointTag

2499 )

2500 lowerLineTag = gmsh.model.occ.addLine(

2501 outerCircle.lowerPointTag, middleCircle.lowerPointTag

2502 )

2503

2504 # Create surfaces for the tube part:

2505 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2506 [

2507 outerCircle.upperRightTag,

2508 rightLineTag,

2509 middleCircle.upperRightTag,

2510 -upperLineTag,

2511 ]

2512 )

2513 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2514 self.tubeSurfaceTags.append(self.upperRightTag)

2515

2516 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2517 [

2518 outerCircle.upperLeftTag,

2519 leftLineTag,

2520 middleCircle.upperLeftTag,

2521 -upperLineTag,

2522 ]

2523 )

2524 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2525 self.tubeSurfaceTags.append(self.upperLeftTag)

2526

2527 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2528 [

2529 outerCircle.lowerLeftTag,

2530 leftLineTag,

2531 middleCircle.lowerLeftTag,

2532 -lowerLineTag,

2533 ]

2534 )

2535 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2536 self.tubeSurfaceTags.append(self.lowerLeftTag)

2537

2538 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2539 [

2540 outerCircle.lowerRightTag,

2541 rightLineTag,

2542 middleCircle.lowerRightTag,

2543 -lowerLineTag,

2544 ]

2545 )

2546 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2547 self.tubeSurfaceTags.append(self.lowerRightTag)

2548

2549 else:

2550 outerCL = outerAir.getInnerCL()

2551

2552 middleCL = gmsh.model.occ.addCircle(

2553 0, 0, z, self.outerRadius - self.tubeThickness

2554 )

2555 middleCL = gmsh.model.occ.addCurveLoop([middleCL])

2556

2557 tubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])

2558 self.tubeSurfaceTags.append(tubeSurface)

2559

2560 # Nontube part:

2561 if self.divideIntoFourParts:

2562 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)

2563

2564 else:

2565 # Inner and outer curve loops. Sometimes the opposite curve loops are used

2566 # because the other one comes from the self.contactTerminalSurface. To create

2567 # a valid surface, the topology of the curve loops should be consistent. See the

2568 # note in the spiralSurface class.

2569 if pancakeIndex % 2 == 0:

2570 innerCL = winding.outerCurveLoopTag

2571 elif pancakeIndex % 2 == 1:

2572 innerCL = winding.outerOppositeCurveLoopTag

2573

2574 # potential bug (curve order might be wrong)

2575 if self.divideIntoFourParts:

2576 middleCL = gmsh.model.occ.addCurveLoop(

2577 [

2578 middleCircle.upperRightTag,

2579 middleCircle.upperLeftTag,

2580 middleCircle.lowerLeftTag,

2581 middleCircle.lowerRightTag,

2582 ]

2583 )

2584

2585 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])

2586 self.nontubeSurfaceTags.append(nontubeSurface)

2587

2588 def createWithWindingAndTubeTags(

2589 self, winding: spiralSurface, tubeTags, pancakeIndex

2590 ):

2591 if not isinstance(tubeTags, list):

2592 raise TypeError("tubeTags must be a list.")

2593

2594 self.tubeSurfaceTags.extend(tubeTags)

2595

2596 middleCurves = gmsh.model.getBoundary(

2597 [(2, tag) for tag in tubeTags], oriented=False

2598 )

2599 middleCurveDimTags = findOuterOnes(middleCurves, findInnerOnes=True)

2600 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]

2601

2602 if self.divideIntoFourParts:

2603 # Create middleCircle object from the tags:

2604 middleCurveTags.sort()

2605 upperRightCurve = middleCurveTags[0]

2606 upperLeftCurve = middleCurveTags[1]

2607 lowerLeftCurve = middleCurveTags[2]

2608 lowerRightCurve = middleCurveTags[3]

2609

2610 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)

2611 pointTags = [dimTag[1] for dimTag in points]

2612 pointTags.sort()

2613 upperPointTag = pointTags[0]

2614 leftPointTag = pointTags[1]

2615

2616 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)

2617 pointTags = [dimTag[1] for dimTag in points]

2618 pointTags.sort()

2619 rightPointTag = pointTags[0]

2620 lowerPointTag = pointTags[1]

2621

2622 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]

2623 middleCircle = circleWithFourCurves(

2624 0,

2625 0,

2626 z,

2627 self.outerRadius - self.tubeThickness,

2628 upperRightTag=upperRightCurve,

2629 upperLeftTag=upperLeftCurve,

2630 lowerLeftTag=lowerLeftCurve,

2631 lowerRightTag=lowerRightCurve,

2632 leftPointTag=leftPointTag,

2633 rightPointTag=rightPointTag,

2634 upperPointTag=upperPointTag,

2635 lowerPointTag=lowerPointTag,

2636 )

2637

2638 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)

2639 else:

2640 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)

2641

2642 if pancakeIndex % 2 == 0:

2643 innerCL = winding.outerCurveLoopTag

2644 elif pancakeIndex % 2 == 1:

2645 innerCL = winding.outerOppositeCurveLoopTag

2646

2647 nontubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])

2648 self.nontubeSurfaceTags.append(nontubeSurface)

2649

2650

2651class innerAirSurface:

2652 def __init__(

2653 self, radius, divideIntoFourParts=False, divideTerminalPartIntoFourParts=False

2654 ):

2655 self.surfaceTags = []

2656

2657 self.divideIntoFourParts = divideIntoFourParts

2658 self.divideTerminalPartIntoFourParts = divideTerminalPartIntoFourParts

2659

2660 self.radius = radius

2661

2662 def createFromScratch(self, z):

2663 self.z = z

2664 if self.divideIntoFourParts:

2665 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)

2666

2667 originTag = point(0, 0, z).tag

2668

2669 leftLineTag = gmsh.model.occ.addLine(

2670 self.outerCircle.leftPointTag, originTag

2671 )

2672 rightLineTag = gmsh.model.occ.addLine(

2673 self.outerCircle.rightPointTag, originTag

2674 )

2675 upperLineTag = gmsh.model.occ.addLine(

2676 self.outerCircle.upperPointTag, originTag

2677 )

2678 lowerLineTag = gmsh.model.occ.addLine(

2679 self.outerCircle.lowerPointTag, originTag

2680 )

2681

2682 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2683 [self.outerCircle.upperRightTag, rightLineTag, -upperLineTag]

2684 )

2685 upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2686 self.surfaceTags.append(upperRightTag)

2687

2688 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2689 [self.outerCircle.upperLeftTag, leftLineTag, -upperLineTag]

2690 )

2691 upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2692 self.surfaceTags.append(upperLeftTag)

2693

2694 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2695 [self.outerCircle.lowerLeftTag, leftLineTag, -lowerLineTag]

2696 )

2697 lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2698 self.surfaceTags.append(lowerLeftTag)

2699

2700 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2701 [self.outerCircle.lowerRightTag, rightLineTag, -lowerLineTag]

2702 )

2703 lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2704 self.surfaceTags.append(lowerRightTag)

2705

2706 else:

2707 if self.divideTerminalPartIntoFourParts:

2708 self.outerCircle = circleWithFourCurves(0, 0, z, self.radius)

2709 outerCL = gmsh.model.occ.addCurveLoop(

2710 [

2711 self.outerCircle.upperRightTag,

2712 self.outerCircle.upperLeftTag,

2713 self.outerCircle.lowerLeftTag,

2714 self.outerCircle.lowerRightTag,

2715 ]

2716 )

2717 else:

2718 outerCL = gmsh.model.occ.addCircle(0, 0, z, self.radius)

2719 outerCL = gmsh.model.occ.addCurveLoop([outerCL])

2720

2721 surfaceTag = gmsh.model.occ.addPlaneSurface([outerCL])

2722 self.surfaceTags.append(surfaceTag)

2723

2724 def setPrecreatedSurfaceTags(self, surfaceTags):

2725 if not isinstance(surfaceTags, list):

2726 raise TypeError("surfaceTags must be a list.")

2727

2728 self.z = gmsh.model.occ.getCenterOfMass(2, surfaceTags[0])[2] # potential bug

2729 self.surfaceTags = []

2730 self.surfaceTags.extend(surfaceTags)

2731

2732 if self.divideIntoFourParts or self.divideTerminalPartIntoFourParts:

2733 # Create outerCirle object from the tags:

2734 curves = gmsh.model.getBoundary(

2735 [(2, tag) for tag in surfaceTags], oriented=False

2736 )

2737 outerCurveDimTags = findOuterOnes(curves)

2738 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]

2739 outerCurveTags.sort()

2740 upperRightCurve = outerCurveTags[0]

2741 upperLeftCurve = outerCurveTags[1]

2742 lowerLeftCurve = outerCurveTags[2]

2743 lowerRightCurve = outerCurveTags[3]

2744

2745 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)

2746 pointTags = [dimTag[1] for dimTag in points]

2747 pointTags.sort()

2748 upperPointTag = pointTags[0]

2749 leftPointTag = pointTags[1]

2750

2751 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)

2752 pointTags = [dimTag[1] for dimTag in points]

2753 pointTags.sort()

2754 rightPointTag = pointTags[0]

2755 lowerPointTag = pointTags[1]

2756

2757 self.outerCircle = circleWithFourCurves(

2758 0,

2759 0,

2760 self.z,

2761 self.radius,

2762 upperRightTag=upperRightCurve,

2763 upperLeftTag=upperLeftCurve,

2764 lowerLeftTag=lowerLeftCurve,

2765 lowerRightTag=lowerRightCurve,

2766 leftPointTag=leftPointTag,

2767 rightPointTag=rightPointTag,

2768 upperPointTag=upperPointTag,

2769 lowerPointTag=lowerPointTag,

2770 )

2771

2772 def getOuterCL(self):

2773 # checked!

2774 # _, curves = gmsh.model.occ.getCurveLoops(self.surfaceTags[0])

2775 # curves = list(curves)

2776 # curves = [int(curves[0])]

2777

2778 # outerCL = gmsh.model.occ.addCurveLoop([curves[0]])

2779

2780 # return outerCL

2781

2782 curves = gmsh.model.getBoundary(

2783 [(2, tag) for tag in self.surfaceTags], oriented=False

2784 )

2785 outerCurveDimTags = findOuterOnes(curves)

2786 outerCurveTags = [dimTag[1] for dimTag in outerCurveDimTags]

2787

2788 outerCL = gmsh.model.occ.addCurveLoop(outerCurveTags)

2789 return outerCL

2790

2791

2792class innerTerminalSurface:

2793 def __init__(self, innerRadius, tubeThickness, divideIntoFourParts=False):

2794 self.tubeSurfaceTags = []

2795 self.nontubeSurfaceTags = []

2796

2797 self.divideIntoFourParts = divideIntoFourParts

2798

2799 self.innerRadius = innerRadius

2800 self.tubeThickness = tubeThickness

2801

2802 def createNontubePartWithMiddleCircleAndWinding(

2803 self, middleCircle: circleWithFourCurves, winding: spiralSurface

2804 ):

2805 leftLineTag = gmsh.model.occ.addLine(

2806 winding.getInnerLeftPointTag(), middleCircle.leftPointTag

2807 )

2808 rightLineTag = gmsh.model.occ.addLine(

2809 winding.getInnerRightPointTag(), middleCircle.rightPointTag

2810 )

2811 upperLineTag = gmsh.model.occ.addLine(

2812 winding.getInnerUpperPointTag(), middleCircle.upperPointTag

2813 )

2814 lowerLineTag = gmsh.model.occ.addLine(

2815 winding.getInnerLowerPointTag(), middleCircle.lowerPointTag

2816 )

2817

2818 # Create surfaces for the nontube part:

2819 if winding.direction is direction.ccw:

2820 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2821 winding.getInnerUpperRightCurves()

2822 + [

2823 upperLineTag,

2824 middleCircle.upperRightTag,

2825 rightLineTag,

2826 ]

2827 )

2828 else:

2829 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2830 winding.getInnerUpperRightCurves()

2831 + [

2832 upperLineTag,

2833 middleCircle.upperRightTag,

2834 rightLineTag,

2835 ]

2836 # + winding.getInnerStartCurves()

2837 )

2838 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2839 self.nontubeSurfaceTags.append(self.upperRightTag)

2840

2841 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2842 winding.getInnerUpperLeftCurves()

2843 + [

2844 upperLineTag,

2845 middleCircle.upperLeftTag,

2846 leftLineTag,

2847 ]

2848 )

2849 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2850 self.nontubeSurfaceTags.append(self.upperLeftTag)

2851

2852 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2853 winding.getInnerLowerLeftCurves()

2854 + [

2855 lowerLineTag,

2856 middleCircle.lowerLeftTag,

2857 leftLineTag,

2858 ]

2859 )

2860 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2861 self.nontubeSurfaceTags.append(self.lowerLeftTag)

2862

2863 if winding.direction is direction.ccw:

2864 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2865 winding.getInnerLowerRightCurves()

2866 + [

2867 lowerLineTag,

2868 middleCircle.lowerRightTag,

2869 rightLineTag,

2870 ]

2871 # + winding.getInnerStartCurves()

2872 )

2873 else:

2874 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2875 winding.getInnerLowerRightCurves()

2876 + [

2877 lowerLineTag,

2878 middleCircle.lowerRightTag,

2879 rightLineTag,

2880 ]

2881 )

2882 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2883 self.nontubeSurfaceTags.append(self.lowerRightTag)

2884

2885 def createWithInnerAirAndWinding(

2886 self, innerAir: innerAirSurface, winding: spiralSurface, pancakeIndex

2887 ):

2888 z = innerAir.z

2889

2890 # Tube part:

2891 if self.divideIntoFourParts:

2892 innerCircle = innerAir.outerCircle

2893 middleCircle = circleWithFourCurves(

2894 0, 0, z, self.innerRadius + self.tubeThickness

2895 )

2896

2897 leftLineTag = gmsh.model.occ.addLine(

2898 middleCircle.leftPointTag, innerCircle.leftPointTag

2899 )

2900 rightLineTag = gmsh.model.occ.addLine(

2901 middleCircle.rightPointTag, innerCircle.rightPointTag

2902 )

2903 upperLineTag = gmsh.model.occ.addLine(

2904 middleCircle.upperPointTag, innerCircle.upperPointTag

2905 )

2906 lowerLineTag = gmsh.model.occ.addLine(

2907 middleCircle.lowerPointTag, innerCircle.lowerPointTag

2908 )

2909

2910 # Create surfaces for the tube part:

2911 upperRightCurveLoop = gmsh.model.occ.addCurveLoop(

2912 [

2913 middleCircle.upperRightTag,

2914 rightLineTag,

2915 innerCircle.upperRightTag,

2916 -upperLineTag,

2917 ]

2918 )

2919 self.upperRightTag = gmsh.model.occ.addPlaneSurface([upperRightCurveLoop])

2920 self.tubeSurfaceTags.append(self.upperRightTag)

2921

2922 upperLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2923 [

2924 middleCircle.upperLeftTag,

2925 leftLineTag,

2926 innerCircle.upperLeftTag,

2927 -upperLineTag,

2928 ]

2929 )

2930 self.upperLeftTag = gmsh.model.occ.addPlaneSurface([upperLeftCurveLoop])

2931 self.tubeSurfaceTags.append(self.upperLeftTag)

2932

2933 lowerLeftCurveLoop = gmsh.model.occ.addCurveLoop(

2934 [

2935 middleCircle.lowerLeftTag,

2936 leftLineTag,

2937 innerCircle.lowerLeftTag,

2938 -lowerLineTag,

2939 ]

2940 )

2941 self.lowerLeftTag = gmsh.model.occ.addPlaneSurface([lowerLeftCurveLoop])

2942 self.tubeSurfaceTags.append(self.lowerLeftTag)

2943

2944 lowerRightCurveLoop = gmsh.model.occ.addCurveLoop(

2945 [

2946 middleCircle.lowerRightTag,

2947 rightLineTag,

2948 innerCircle.lowerRightTag,

2949 -lowerLineTag,

2950 ]

2951 )

2952 self.lowerRightTag = gmsh.model.occ.addPlaneSurface([lowerRightCurveLoop])

2953 self.tubeSurfaceTags.append(self.lowerRightTag)

2954

2955 else:

2956 innerCL = innerAir.getOuterCL()

2957

2958 middleCL = gmsh.model.occ.addCircle(

2959 0, 0, z, self.innerRadius + self.tubeThickness

2960 )

2961 middleCL = gmsh.model.occ.addCurveLoop([middleCL])

2962

2963 tubeSurface = gmsh.model.occ.addPlaneSurface([middleCL, innerCL])

2964 self.tubeSurfaceTags.append(tubeSurface)

2965

2966 # Nontube part:

2967 if self.divideIntoFourParts:

2968 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)

2969

2970 else:

2971 # Inner and outer curve loops. Sometimes the opposite curve loops are used

2972 # because the other one comes from the self.contactTerminalSurface. To create

2973 # a valid surface, the topology of the curve loops should be consistent. See the

2974 # note in the spiralSurface class.

2975 if pancakeIndex == 0:

2976 outerCL = winding.innerCurveLoopTag

2977

2978 elif pancakeIndex % 2 == 0:

2979 outerCL = winding.innerOppositeCurveLoopTag

2980

2981 elif pancakeIndex % 2 == 1:

2982 outerCL = winding.innerCurveLoopTag

2983

2984 # potential bug (curve order might be wrong)

2985 if self.divideIntoFourParts:

2986 middleCL = gmsh.model.occ.addCurveLoop(

2987 [

2988 middleCircle.upperRightTag,

2989 middleCircle.upperLeftTag,

2990 middleCircle.lowerLeftTag,

2991 middleCircle.lowerRightTag,

2992 ]

2993 )

2994

2995 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])

2996 self.nontubeSurfaceTags.append(nontubeSurface)

2997

2998 def createWithWindingAndTubeTags(

2999 self, winding: spiralSurface, tubeTags, pancakeIndex

3000 ):

3001 if not isinstance(tubeTags, list):

3002 raise TypeError("tubeTags must be a list.")

3003

3004 self.tubeSurfaceTags.extend(tubeTags)

3005

3006 middleCurves = gmsh.model.getBoundary(

3007 [(2, tag) for tag in tubeTags], oriented=False

3008 )

3009 middleCurveDimTags = findOuterOnes(middleCurves)

3010 middleCurveTags = [dimTag[1] for dimTag in middleCurveDimTags]

3011

3012 if self.divideIntoFourParts:

3013 # Create middleCircle object from the tags:

3014 middleCurveTags.sort()

3015 upperRightCurve = middleCurveTags[0]

3016 upperLeftCurve = middleCurveTags[1]

3017 lowerLeftCurve = middleCurveTags[2]

3018 lowerRightCurve = middleCurveTags[3]

3019

3020 points = gmsh.model.getBoundary([(1, upperLeftCurve)], oriented=False)

3021 pointTags = [dimTag[1] for dimTag in points]

3022 pointTags.sort()

3023 upperPointTag = pointTags[0]

3024 leftPointTag = pointTags[1]

3025

3026 points = gmsh.model.getBoundary([(1, lowerRightCurve)], oriented=False)

3027 pointTags = [dimTag[1] for dimTag in points]

3028 pointTags.sort()

3029 rightPointTag = pointTags[0]

3030 lowerPointTag = pointTags[1]

3031

3032 z = gmsh.model.occ.getCenterOfMass(1, upperRightCurve)[2]

3033 middleCircle = circleWithFourCurves(

3034 0,

3035 0,

3036 z,

3037 self.innerRadius + self.tubeThickness,

3038 upperRightTag=upperRightCurve,

3039 upperLeftTag=upperLeftCurve,

3040 lowerLeftTag=lowerLeftCurve,

3041 lowerRightTag=lowerRightCurve,

3042 leftPointTag=leftPointTag,

3043 rightPointTag=rightPointTag,

3044 upperPointTag=upperPointTag,

3045 lowerPointTag=lowerPointTag,

3046 )

3047

3048 self.createNontubePartWithMiddleCircleAndWinding(middleCircle, winding)

3049 else:

3050 middleCL = gmsh.model.occ.addCurveLoop(middleCurveTags)

3051

3052 if pancakeIndex == 0:

3053 outerCL = winding.innerCurveLoopTag

3054

3055 elif pancakeIndex % 2 == 0:

3056 outerCL = winding.innerOppositeCurveLoopTag

3057

3058 elif pancakeIndex % 2 == 1:

3059 outerCL = winding.innerCurveLoopTag

3060

3061 nontubeSurface = gmsh.model.occ.addPlaneSurface([outerCL, middleCL])

3062 self.nontubeSurfaceTags.append(nontubeSurface)

3063

3064

3065class pancakeCoilsWithAir:

3066 """

3067 A class to create Pancake3D coil. With this class, any number of pancake coils stack

3068 can be created in GMSH. Moreover, the class also creates some parts of the air

3069 volume as well. It creates the inner cylinder air volume and the outer tube air

3070 volume. However, the air between the pancake coils is not created. It is created in

3071 the gapAir class.

3072

3073 self.fundamentalSurfaces are the surfaces at the bottom of each pancake coil. They

3074 are created one by one with self.generateFundamentalSurfaces() method. The first

3075 created self.fundamentalSurfaces are the first pancake coil's bottom surfaces. Those

3076 surfaces include the outer air shell surface, outer air tube surface, outer

3077 terminal's outer tube part, outer terminal's touching part, winding surfaces,

3078 contact layer surfaces, inner terminal's touching part, inner terminal's inner tube

3079 part, and inner air disc. Terminals are divided into two because they can only be

3080 connected with the perfect tubes since each pancake coil is rotated in a different

3081 direction.

3082

3083 For the first pancake coil, self.fundamentalSurfaces are extruded downwards to

3084 connect the terminal to the end of the geometry at the bottom.

3085

3086 Then self.fundamentalSurfaces are extruded upwards to create the first pancake coil

3087 with self.extrudeWindingPart method. The method returns the extrusion's top

3088 surfaces, which are saved in the topSurfaces variable.

3089

3090 Then those topSurfaces variable is given to another method, self.extrudeGapPart, and

3091 they are further extruded upwards up to the bottom of the next pancake coil.

3092 However, only the air tube, air cylinder, and connection terminal (the perfect inner

3093 terminal tube or outer terminal tube) are extruded. Otherwise, conformality would be

3094 impossible. The gaps are filled with air in gapAir class with fragment operation

3095 later. Then the top surfaces are returned by the method and saved in

3096 self.contactSurfaces variable.

3097

3098 Then using the self.contactSurfaces, self.generateFundamentalSurfaces method creates

3099 the new fundamental surfaces. All the surfaces from the self.contactSurfaces are

3100 used in the new self.fundamentalSurfaces variable to avoid surface duplication.

3101

3102 The logic goes until the last topSurfaces are extruded upwards to connect the last

3103 terminal to the top of the geometry.

3104

3105 Every pancake coil's rotation direction is different each time. Otherwise, their

3106 magnetic fields would neutralize each other.

3107

3108 The first and second pancake coils are connected with the inner terminal. Then the

3109 second and the third pancake coils are connected with the outer terminal. And so on.

3110

3111 :param geometryData: geometry information

3112 """

3113

3114 def __init__(self, geometryData, meshData) -> None:

3115 logger.info("Generating pancake coils has been started.")

3116 start_time = timeit.default_timer()

3117

3118 # Data:

3119 self.geo = geometryData

3120 self.mesh = meshData

3121

3122 # ==============================================================================

3123 # CREATING VOLUME STORAGES STARTS ==============================================

3124 # ==============================================================================

3125 # Air shell (they will be empty if shellTransformation == False):

3126 # For cylinder type:

3127 self.airShellVolume = dimTags(name=self.geo.air.shellVolumeName, save=True)

3128

3129 # For cuboid type:

3130 self.airShellVolumePart1 = dimTags(

3131 name=self.geo.air.shellVolumeName + "-Part1", save=True

3132 )

3133 self.airShellVolumePart2 = dimTags(

3134 name=self.geo.air.shellVolumeName + "-Part2", save=True

3135 )

3136 self.airShellVolumePart3 = dimTags(

3137 name=self.geo.air.shellVolumeName + "-Part3", save=True

3138 )

3139 self.airShellVolumePart4 = dimTags(

3140 name=self.geo.air.shellVolumeName + "-Part4", save=True

3141 )

3142

3143 # Outer air tube volume (actually it is not a tube if the air type is cuboid):

3144 self.outerAirTubeVolume = dimTags(

3145 name=self.geo.air.name + "-OuterTube", save=True, parentName=self.geo.air.name

3146 )

3147

3148 # Outer terminal's outer tube part:

3149 self.outerTerminalTubeVolume = dimTags(

3150 name=self.geo.terminals.outer.name + "-Tube", save=True, parentName=self.geo.terminals.outer.name

3151 )

3152

3153 # Outer terminal's volume that touches the winding:

3154 self.outerTerminalTouchingVolume = dimTags(

3155 name=self.geo.terminals.outer.name + "-Touching",

3156 save=True,

3157 parentName=self.geo.terminals.outer.name,

3158 )

3159

3160 # Inner terminal's volume that touches the winding:

3161 self.innerTerminalTouchingVolume = dimTags(

3162 name=self.geo.terminals.inner.name + "-Touching",

3163 save=True,

3164 parentName=self.geo.terminals.inner.name,

3165 )

3166

3167 # Inner terminal's inner tube part:

3168 self.innerTerminalTubeVolume = dimTags(

3169 name=self.geo.terminals.inner.name + "-Tube", save=True, parentName=self.geo.terminals.inner.name

3170 )

3171

3172 # Transition layers:

3173 self.innerTransitionNotchVolume = dimTags(

3174 name="innerTransitionNotch",

3175 save=True,

3176 )

3177 self.outerTransitionNotchVolume = dimTags(

3178 name="outerTransitionNotch",

3179 save=True,

3180 )

3181

3182 # Inner air cylinder volume:

3183 self.centerAirCylinderVolume = dimTags(

3184 name=self.geo.air.name + "-InnerCylinder",

3185 save=True,

3186 parentName=self.geo.air.name,

3187 )

3188

3189 # Top and bottom parts of the air volume:

3190 self.topAirPancakeWindingExtursionVolume = dimTags(

3191 name=self.geo.air.name + "-TopPancakeWindingExtursion",

3192 save=True,

3193 parentName=self.geo.air.name,

3194 )

3195 self.topAirPancakeContactLayerExtursionVolume = dimTags(

3196 name=self.geo.air.name + "-TopPancakeContactLayerExtursion",

3197 save=True,

3198 parentName=self.geo.air.name,

3199 )

3200 self.topAirTerminalsExtrusionVolume = dimTags(

3201 name=self.geo.air.name + "-TopTerminalsExtrusion",

3202 save=True,

3203 parentName=self.geo.air.name,

3204 )

3205 self.topAirTubeTerminalsExtrusionVolume = dimTags(

3206 name=self.geo.air.name + "-TopTubeTerminalsExtrusion",

3207 save=True,

3208 parentName=self.geo.air.name,

3209 )

3210

3211 self.bottomAirPancakeWindingExtursionVolume = dimTags(

3212 name=self.geo.air.name + "-BottomPancakeWindingExtursion",

3213 save=True,

3214 parentName=self.geo.air.name,

3215 )

3216 self.bottomAirPancakeContactLayerExtursionVolume = dimTags(

3217 name=self.geo.air.name + "-BottomPancakeContactLayerExtursion",

3218 save=True,

3219 parentName=self.geo.air.name,

3220 )

3221 self.bottomAirTerminalsExtrusionVolume = dimTags(

3222 name=self.geo.air.name + "-BottomTerminalsExtrusion",

3223 save=True,

3224 parentName=self.geo.air.name,

3225 )

3226 self.bottomAirTubeTerminalsExtrusionVolume = dimTags(

3227 name=self.geo.air.name + "-BottomTubeTerminalsExtrusion",

3228 save=True,

3229 parentName=self.geo.air.name,

3230 )

3231

3232 # Gap air:

3233 self.gapAirVolume = dimTags(

3234 name=self.geo.air.name + "-Gap", save=True, parentName=self.geo.air.name

3235 )

3236

3237 # Create additional/optional volume storages (they might be used in the meshing

3238 # process):

3239 self.firstTerminalVolume = dimTags(name=self.geo.terminals.firstName, save=True)

3240 self.lastTerminalVolume = dimTags(name=self.geo.terminals.lastName, save=True)

3241

3242 # ==============================================================================

3243 # CREATING VOLUME STORAGES ENDS ================================================

3244 # ==============================================================================

3245

3246 # self.fundamentalSurfaces is a dictionary of surface dimTags tuples. The keys

3247 # are the dimTags objects of the corresponding volumes. The values are the

3248 # dimTags tuples of the surfaces that are used to extrude the volumes. It is

3249 # created in self.generateFundamentalSurfaces method.

3250 self.fundamentalSurfaces = {}

3251

3252 # self.pancakeIndex stores the index of the current pancake coil.

3253 self.pancakeIndex = 0

3254

3255 # self.contactSurfaces is a dictionary of surface dimTags tuples. The keys are

3256 # the dimTags objects of the corresponding volumes. The values are the dimTags

3257 # tuples of the surfaces that are obtained from the previous extrusion and used

3258 # for the next extrusion. The same surface is used for the next extrusion to

3259 # avoid surface duplication. It is created in self.extrudeGapPart and

3260 # self.extrudeWindingPart methods.

3261 self.contactSurfaces = {}

3262

3263 # They will be lists of dimTags objects:

3264 self.individualWinding = []

3265 self.individualContactLayer = []

3266

3267 self.gapAirSurfacesDimTags = []

3268

3269 for i in range(self.geo.numberOfPancakes):

3270 # Itterate over the number of pancake coils:

3271 self.individualWinding.append(

3272 dimTags(

3273 name=self.geo.winding.name + str(self.pancakeIndex + 1),

3274 save=True,

3275 parentName=self.geo.winding.name,

3276 )

3277 )

3278 self.individualContactLayer.append(

3279 dimTags(

3280 name=self.geo.contactLayer.name + str(self.pancakeIndex + 1),

3281 save=True,

3282 parentName=self.geo.contactLayer.name,

3283 )

3284 )

3285

3286 # Generate the fundamental surfaces:

3287 self.fundamentalSurfaces = self.generateFundamentalSurfaces()

3288

3289 # Create gap air or collect the gap air surfaces:

3290 if i != 0:

3291 bottomSurfacesDimTags = []

3292 for key, value in topSurfaces.items():

3293 if (

3294 key is self.individualWinding[self.pancakeIndex - 1]

3295 or key is self.individualContactLayer[self.pancakeIndex - 1]

3296 or key is self.outerTerminalTouchingVolume

3297 or key is self.innerTerminalTouchingVolume

3298 or key is self.innerTransitionNotchVolume

3299 or key is self.outerTransitionNotchVolume

3300 ):

3301 bottomSurfacesDimTags.extend(value)

3302

3303 topSurfacesDimTags = []

3304 for key, value in self.fundamentalSurfaces.items():

3305 if (

3306 key is self.individualWinding[self.pancakeIndex]

3307 or key is self.individualContactLayer[self.pancakeIndex]

3308 or key is self.outerTerminalTouchingVolume

3309 or key is self.innerTerminalTouchingVolume

3310 or key is self.innerTransitionNotchVolume

3311 or key is self.outerTransitionNotchVolume

3312 ):

3313 topSurfacesDimTags.extend(value)

3314

3315 sideSurfacesDimTags = []

3316 if i % 2 == 1:

3317 # Touches it tube and air tube

3318 bottomSurfacesDimTags.extend(

3319 topSurfaces[self.outerTerminalTubeVolume]

3320 )

3321 topSurfacesDimTags.extend(

3322 self.fundamentalSurfaces[self.outerTerminalTubeVolume]

3323 )

3324

3325 if self.mesh.terminals.structured:

3326 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(

3327 3

3328 )[-4:]

3329 else:

3330 lastItTubeVolDimTags = self.innerTerminalTubeVolume.getDimTags(

3331 3

3332 )[-1:]

3333

3334 lastItTubeSurfsDimTags = gmsh.model.getBoundary(

3335 lastItTubeVolDimTags, oriented=False

3336 )

3337 lastItTubeSideSurfsDimTags = findSurfacesWithNormalsOnXYPlane(

3338 lastItTubeSurfsDimTags

3339 )

3340 sideSurfacesDimTags.extend(

3341 findOuterOnes(lastItTubeSideSurfsDimTags)

3342 )

3343

3344 if self.mesh.air.structured:

3345 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[

3346 -4:

3347 ]

3348 else:

3349 lastAirTubeVolDimTags = self.outerAirTubeVolume.getDimTags(3)[

3350 -1:

3351 ]

3352 lastAirTubeSurfsDimTags = gmsh.model.getBoundary(

3353 lastAirTubeVolDimTags, oriented=False

3354 )

3355 lastAirTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(

3356 lastAirTubeSurfsDimTags

3357 )

3358 sideSurfacesDimTags.extend(

3359 findOuterOnes(lastAirTubeSurfsDimTags, findInnerOnes=True)

3360 )

3361

3362 else:

3363 # Touches ot tube and air cylinder

3364 bottomSurfacesDimTags.extend(

3365 topSurfaces[self.innerTerminalTubeVolume]

3366 )

3367 topSurfacesDimTags.extend(

3368 self.fundamentalSurfaces[self.innerTerminalTubeVolume]

3369 )

3370 if self.mesh.terminals.structured:

3371 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(

3372 3

3373 )[-4:]

3374 else:

3375 lastOtTubeVolDimTags = self.outerTerminalTubeVolume.getDimTags(

3376 3

3377 )[-1:]

3378

3379 lastOtTubeSurfsDimTags = gmsh.model.getBoundary(

3380 lastOtTubeVolDimTags, oriented=False

3381 )

3382 lastOtTubeSurfsDimTags = findSurfacesWithNormalsOnXYPlane(

3383 lastOtTubeSurfsDimTags

3384 )

3385 sideSurfacesDimTags.extend(

3386 findOuterOnes(lastOtTubeSurfsDimTags, findInnerOnes=True)

3387 )

3388

3389 if self.mesh.air.structured:

3390 lastAirCylinderVolDimTags = (

3391 self.centerAirCylinderVolume.getDimTags(3)[-4:]

3392 )

3393 else:

3394 lastAirCylinderVolDimTags = (

3395 self.centerAirCylinderVolume.getDimTags(3)[-1:]

3396 )

3397

3398 lastAirCylinderSurfsDimTags = gmsh.model.getBoundary(

3399 lastAirCylinderVolDimTags, oriented=False

3400 )

3401 lastAirCylinderSurfsDimTags = findSurfacesWithNormalsOnXYPlane(

3402 lastAirCylinderSurfsDimTags

3403 )

3404 sideSurfacesDimTags.extend(

3405 findOuterOnes(lastAirCylinderSurfsDimTags)

3406 )

3407

3408 allGapAirSurfacesDimTags = (

3409 bottomSurfacesDimTags + topSurfacesDimTags + sideSurfacesDimTags

3410 )

3411

3412 # Technically, since all the boundary surfaces of the gap air volumes

3413 # are found here, we should be able to create the gap air volumes with

3414 # addSurfaceLoop and addVolume functions. However, when those are used,

3415 # Geometry.remove_all_duplicates() will indicate some

3416 # duplicates/ill-shaped geometry entities. The indication is

3417 # gmsh.model.occ.remove_all_duplicates() will change the geometry

3418 # (delete some volumes and create new ones), and I have always thought

3419 # that means there are big errors in the geometry and that geometry

3420 # should not be used.

3421

3422 # Alternatively, using these surface tags, the gap air can be created

3423 # with fragment operations as well. Geometry.remove_all_duplicates()

3424 # will tell everything is fine when the fragment operation is used.

3425

3426 # However, I checked manually as well, the way I am using the

3427 # addSurfaceLoop and addVolume should definitely work (because the end

3428 # result is the same with fragments), and I think it is a gmsh/occ

3429 # related problem. In the end, I realized creating the gap air with

3430 # addSurfaceLoop and addVolume won't even affect the mesh, and

3431 # everything seems conformal and nice. Since the fragment operation

3432 # is also very slow, I decided to use addSurfaceLoop and addVolume them.

3433 # However, I keep it as an option so that if the user feels something

3434 # funny about the geometry, the gap air can be created with fragment

3435 # operations as well.

3436

3437 if not self.geo.air.generateGapAirWithFragment:

3438 allGapAirSurfacesTags = [

3439 dimTag[1] for dimTag in allGapAirSurfacesDimTags

3440 ]

3441 surfaceLoop = gmsh.model.occ.addSurfaceLoop(allGapAirSurfacesTags)

3442 volume = gmsh.model.occ.addVolume([surfaceLoop])

3443 self.gapAirVolume.add([(3, volume)])

3444

3445 else:

3446 # Save the surface tags for a fast fragment operation:

3447 self.gapAirSurfacesDimTags.append(allGapAirSurfacesDimTags)

3448

3449 # self.extrudeSurfaces uses self.fundamentalSurfaces for extrusion and adds

3450 # the new volumes to the dimTags objects and returns the dictionary of the

3451 # new top surfaces. The new top surfaces then will be used in extrudeGapPart

3452 # method.

3453 topSurfaces = self.extrudeWindingPart()

3454

3455 if i == 0:

3456 # If it is the first pancake coil, fundemental surfaces are extruded

3457 # downwards to create the bottom air volume and terminal volume.

3458 _ = self.extrudeGapPart(

3459 self.fundamentalSurfaces,

3460 -self.geo.air.axialMargin,

3461 terminalDimTagsObject=self.outerTerminalTubeVolume,

3462 firstTerminal=True,

3463 )

3464

3465 if not i == self.geo.numberOfPancakes - 1:

3466 # If it is not the last pancake coil, extrude the terminal surface to

3467 # create the next contactTerminalSurface and store the new volume in the

3468 # corresponding dimTags object.

3469 self.contactSurfaces = self.extrudeGapPart(topSurfaces)

3470

3471 else:

3472 # If it is the last pancake coil, extrude the terminal surface all the

3473 # way up to the top and store the new volume in the corresponding

3474 # dimTags object.

3475 _ = self.extrudeGapPart(

3476 topSurfaces,

3477 self.geo.air.axialMargin,

3478 lastTerminal=True,

3479 )

3480 self.pancakeIndex = self.pancakeIndex + 1

3481

3482 # Create the gap air volume:

3483 if self.geo.air.generateGapAirWithFragment and self.geo.numberOfPancakes > 1:

3484 self.generateGapAirWithFragment(self.gapAirSurfacesDimTags)

3485

3486 logger.info(

3487 "Generating pancake coils has been finished in"

3488 f" {timeit.default_timer() - start_time:.2f} s."

3489 )

3490

3491 def generateFundamentalSurfaces(self):

3492 """

3493 Generates the inner air, outer air, winding, contact layer, and terminal surfaces

3494 of the current pancake coil and returns them. It finds the z coordinate of the

3495 surfaces and the direction of the pancake coil, depending on the pancake index.

3496

3497 :return: list of dimTags that contains fundamental surfaces

3498 :rtype: list[tuple[int, int]]

3499 """

3500 fundamentalSurfaces = {}

3501

3502 # Select the direction of the spiral:

3503 if self.pancakeIndex % 2 == 0:

3504 spiralDirection = direction.ccw

3505 else:

3506 spiralDirection = direction.cw

3507

3508 # Calculate the z coordinate of the surfaces:

3509 z = (

3510 -self.geo.air.h / 2

3511 + self.geo.air.axialMargin

3512 + self.pancakeIndex * (self.geo.winding.height + self.geo.gapBetweenPancakes)

3513 )

3514

3515 # Create the winding and contact layer surface:

3516 surface = spiralSurface(

3517 self.geo.winding.innerRadius,

3518 self.geo.winding.thickness,

3519 self.geo.contactLayer.thickness,

3520 self.geo.winding.numberOfTurns,

3521 z,

3522 self.geo.winding.theta_i,

3523 self.geo.terminals.transitionNotchAngle,

3524 spiralDirection,

3525 thinShellApproximation=self.geo.contactLayer.thinShellApproximation,

3526 )

3527 #print("this is the surface")

3528 #print(surface)

3529 # Save the surface tags (if TSA, contactLayerSurfaceTags will be empty):

3530 fundamentalSurfaces[self.individualWinding[self.pancakeIndex]] = [

3531 (2, tag) for tag in surface.surfaceTags

3532 ]

3533 fundamentalSurfaces[self.individualContactLayer[self.pancakeIndex]] = [

3534 (2, tag) for tag in surface.contactLayerSurfaceTags

3535 ]

3536

3537 if self.geo.air.type == "cylinder":

3538 outerAirSurf = outerAirSurface(

3539 self.geo.air.radius,

3540 self.geo.terminals.outer.r,

3541 type="cylinder",

3542 divideIntoFourParts=self.mesh.air.structured,

3543 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,

3544 )

3545 elif self.geo.air.type == "cuboid":

3546 outerAirSurf = outerAirSurface(

3547 self.geo.air.sideLength / 2,

3548 self.geo.terminals.outer.r,

3549 type="cuboid",

3550 divideIntoFourParts=self.mesh.air.structured,

3551 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,

3552 )

3553

3554 outerTerminalSurf = outerTerminalSurface(

3555 self.geo.terminals.outer.r,

3556 self.geo.terminals.outer.thickness,

3557 divideIntoFourParts=self.mesh.terminals.structured,

3558 )

3559 innerTerminalSurf = innerTerminalSurface(

3560 self.geo.terminals.inner.r,

3561 self.geo.terminals.inner.thickness,

3562 divideIntoFourParts=self.mesh.terminals.structured,

3563 )

3564 innerAirSurf = innerAirSurface(

3565 self.geo.terminals.inner.r,

3566 divideIntoFourParts=self.mesh.air.structured,

3567 divideTerminalPartIntoFourParts=self.mesh.terminals.structured,

3568 )

3569

3570 if self.contactSurfaces:

3571 # If self.contactSurfaces is not empty, it means that it is not the

3572 # first pancake coil. In that case, contactSurfaces should be used to

3573 # avoid surface duplication.

3574

3575 # Create outer air:

3576 outerAirPSDimTags = self.contactSurfaces[self.outerAirTubeVolume]

3577 outerAirPSTags = [dimTag[1] for dimTag in outerAirPSDimTags]

3578 if self.geo.air.shellTransformation:

3579 if self.geo.air.type == "cuboid":

3580 cuboidShellDimTags1 = self.contactSurfaces[self.airShellVolumePart1]

3581 cuboidShellTags1 = [dimTag[1] for dimTag in cuboidShellDimTags1]

3582 cuboidShellDimTags2 = self.contactSurfaces[self.airShellVolumePart2]

3583 cuboidShellTags2 = [dimTag[1] for dimTag in cuboidShellDimTags2]

3584 cuboidShellDimTags3 = self.contactSurfaces[self.airShellVolumePart3]

3585 cuboidShellTags3 = [dimTag[1] for dimTag in cuboidShellDimTags3]

3586 cuboidShellDimTags4 = self.contactSurfaces[self.airShellVolumePart4]

3587 cuboidShellTags4 = [dimTag[1] for dimTag in cuboidShellDimTags4]

3588 outerAirSurf.setPrecreatedSurfaceTags(

3589 outerAirPSTags,

3590 cuboidShellTags1=cuboidShellTags1,

3591 cuboidShellTags2=cuboidShellTags2,

3592 cuboidShellTags3=cuboidShellTags3,

3593 cuboidShellTags4=cuboidShellTags4,

3594 )

3595 elif self.geo.air.type == "cylinder":

3596 cylinderShellDimTags = self.contactSurfaces[self.airShellVolume]

3597 cylinderShellTags = [dimTag[1] for dimTag in cylinderShellDimTags]

3598 outerAirSurf.setPrecreatedSurfaceTags(

3599 outerAirPSTags,

3600 cylinderShellTags=cylinderShellTags,

3601 )

3602 else:

3603 outerAirSurf.setPrecreatedSurfaceTags(outerAirPSTags)

3604

3605 # Create inner air:

3606 innerAirPSDimTags = self.contactSurfaces[self.centerAirCylinderVolume]

3607 innerAirPSTags = [dimTag[1] for dimTag in innerAirPSDimTags]

3608 innerAirSurf.setPrecreatedSurfaceTags(innerAirPSTags)

3609

3610 if self.pancakeIndex % 2 == 0:

3611 # In this case, we should create all the surfaces for the inner terminal

3612 # but not for outer terminal. Because it is a pancake coil with an even

3613 # index (self.pancakeIndex%2==0) which means that it is connected to the

3614 # previous pancake coil with outer terminal and the outer terminal

3615 # surface is ready (extruded before).

3616

3617 # Create outer terminal:

3618 outerTerminalTubePSDimTags = self.contactSurfaces[

3619 self.outerTerminalTubeVolume

3620 ]

3621 outerTerminalTubePSTags = [

3622 dimTag[1] for dimTag in outerTerminalTubePSDimTags

3623 ]

3624 outerTerminalSurf.createWithWindingAndTubeTags(

3625 surface, outerTerminalTubePSTags, self.pancakeIndex

3626 )

3627

3628 # Create inner terminal:

3629 innerTerminalSurf.createWithInnerAirAndWinding(

3630 innerAirSurf, surface, self.pancakeIndex

3631 )

3632

3633 else:

3634 # In this case, we should create all the surfaces for the outer terminal

3635 # but not for inner terminal. Because it is a pancake coil with an odd

3636 # index (self.pancakeIndex%2==1) which means that it is connected to the

3637 # previous pancake coil with inner terminal and the inner terminal

3638 # surface is ready (extruded before).

3639

3640 # Create outer terminal:

3641 outerTerminalSurf.createWithOuterAirAndWinding(

3642 outerAirSurf, surface, self.pancakeIndex

3643 )

3644

3645 # Create inner terminal:

3646 innerTerminalTubePSDimTags = self.contactSurfaces[

3647 self.innerTerminalTubeVolume

3648 ]

3649 innerTerminalTubePSTags = [

3650 dimTag[1] for dimTag in innerTerminalTubePSDimTags

3651 ]

3652 innerTerminalSurf.createWithWindingAndTubeTags(

3653 surface, innerTerminalTubePSTags, self.pancakeIndex

3654 )

3655

3656 else:

3657 # If self.contactSurfaces is empty, it means that it is the first pancake

3658 # coil. In that case, the surfaces should be created from scratch.

3659

3660 if self.geo.air.shellTransformation:

3661 if self.geo.air.type == "cuboid":

3662 outerAirSurf.createFromScratch(

3663 z,

3664 shellTransformation=True,

3665 shellRadius=self.geo.air.shellSideLength / 2,

3666 )

3667 else:

3668 outerAirSurf.createFromScratch(

3669 z,

3670 shellTransformation=True,

3671 shellRadius=self.geo.air.shellOuterRadius,

3672 )

3673 else:

3674 outerAirSurf.createFromScratch(z)

3675

3676 innerAirSurf.createFromScratch(z)

3677 outerTerminalSurf.createWithOuterAirAndWinding(

3678 outerAirSurf, surface, self.pancakeIndex

3679 )

3680 innerTerminalSurf.createWithInnerAirAndWinding(

3681 innerAirSurf, surface, self.pancakeIndex

3682 )

3683

3684 # Save the surface tags:

3685 fundamentalSurfaces[self.outerAirTubeVolume] = [

3686 (2, tag) for tag in outerAirSurf.surfaceTags

3687 ]

3688

3689 fundamentalSurfaces[self.centerAirCylinderVolume] = [

3690 (2, tag) for tag in innerAirSurf.surfaceTags

3691 ]

3692

3693 fundamentalSurfaces[self.outerTerminalTubeVolume] = [

3694 (2, tag) for tag in outerTerminalSurf.tubeSurfaceTags

3695 ]

3696 fundamentalSurfaces[self.outerTerminalTouchingVolume] = [

3697 (2, tag) for tag in outerTerminalSurf.nontubeSurfaceTags

3698 ]

3699

3700 fundamentalSurfaces[self.innerTerminalTubeVolume] = [

3701 (2, tag) for tag in innerTerminalSurf.tubeSurfaceTags

3702 ]

3703 fundamentalSurfaces[self.innerTerminalTouchingVolume] = [

3704 (2, tag) for tag in innerTerminalSurf.nontubeSurfaceTags

3705 ]

3706 fundamentalSurfaces[self.innerTransitionNotchVolume] = [

3707 (2, tag) for tag in surface.innerNotchSurfaceTags

3708 ]

3709 fundamentalSurfaces[self.outerTransitionNotchVolume] = [

3710 (2, tag) for tag in surface.outerNotchSurfaceTags

3711 ]

3712

3713 if self.geo.air.shellTransformation:

3714 if self.geo.air.type == "cuboid":

3715 fundamentalSurfaces[self.airShellVolumePart1] = [

3716 (2, tag) for tag in outerAirSurf.shellTagsPart1

3717 ]

3718 fundamentalSurfaces[self.airShellVolumePart2] = [

3719 (2, tag) for tag in outerAirSurf.shellTagsPart2

3720 ]

3721 fundamentalSurfaces[self.airShellVolumePart3] = [

3722 (2, tag) for tag in outerAirSurf.shellTagsPart3

3723 ]

3724 fundamentalSurfaces[self.airShellVolumePart4] = [

3725 (2, tag) for tag in outerAirSurf.shellTagsPart4

3726 ]

3727 elif self.geo.air.type == "cylinder":

3728 fundamentalSurfaces[self.airShellVolume] = [

3729 (2, tag) for tag in outerAirSurf.shellTags

3730 ]

3731 # windingSurfaces = {}

3732 # # Save only the winding surface tags:

3733 # windingSurfaces[self.individualWinding[self.pancakeIndex]] = [

3734 # (2, tag) for tag in surface.surfaceTags

3735 # ]

3736 # fundamentalSurfaces = windingSurfaces

3737 # print(fundamentalSurfaces)

3738

3739 return fundamentalSurfaces

3740

3741 def extrudeGapPart(

3742 self,

3743 surfacesDict,

3744 tZ: float = None,

3745 terminalDimTagsObject: dimTags = None,

3746 firstTerminal=False,

3747 lastTerminal=False,

3748 ):

3749 """

3750 Extrudes the given surfaces dimTags dictionary to a given height (tZ) and adds

3751 the created volumes to the corresponding dictionary keys (dimTags objects). It

3752 returns the extrusion's top surfaces as a dictionary again, where the keys are

3753 the corresponding dimTagsObjects and the values are the dimTags of the surfaces.

3754

3755 If tZ is not given, then it is set to the gap height (self.geo.gapBetweenPancakes). This is the

3756 default value used for connecting the pancake coils. Only for the creation of

3757 the first and the last pancake coils different tZ values are used.

3758

3759 If terminalDimTagsObject is not given, then the created volume is added

3760 automatically to the innerTerminalVolume or outerTerminalVolume dimTags object,

3761 depending on the value of self.pancakeIndex. However, giving

3762 terminalDimTagsObject is necessary for creating the first and the last terminal.

3763 Otherwise, finding out the correct terminal dimTagsObject would be very

3764 challenging.

3765

3766 :param surfaces: the surface dimTag dictionary to be extruded. The keys are the

3767 dimTags objects and the values are the dimTags of the surfaces. The

3768 keys are used to easily add the corresponding volumes to the correct

3769 dimTags objects

3770 :type surfaces: dict[dimTags, list[tuple[int, int]]]

3771 :param tZ: the height of the extrusion

3772 :type tZ: float, optional

3773 :param terminalDimTagsObject: the dimTags object of the terminal to be extruded

3774 :type terminalDimTagsObject: dimTags, optional

3775 :return: top surfaces of the extrusion as a dictionary where the keys are the

3776 dimTags objects and the values are the dimTags of the surfaces

3777 :rtype: dict[dimTags, list[tuple[int, int]]]

3778 """

3779 bottomPart = False

3780 topPart = False

3781 if tZ is None:

3782 tZ = self.geo.gapBetweenPancakes

3783 elif tZ < 0:

3784 bottomPart = True

3785 elif tZ > 0:

3786 topPart = True

3787

3788 if terminalDimTagsObject is None:

3789 # terminalDimTagsObject needs to be given for the first terminal that is

3790 # extruded downwards.

3791 if self.pancakeIndex % 2 == 0:

3792 terminalDimTagsObject = self.innerTerminalTubeVolume

3793 else:

3794 terminalDimTagsObject = self.outerTerminalTubeVolume

3795

3796 # if terminalDimTagsObject is self.innerTerminalVolume:

3797 # otherTerminal = self.outerTerminalVolume

3798 # else:

3799 # otherTerminal = self.innerTerminalVolume

3800

3801 # Create the list of surfaces to be extruded:

3802 listOfDimTags = []

3803 listOfDimTagsObjects = []

3804 listOfDimTagsForTopSurfaces = []

3805 if topPart:

3806 # Then in this case, most of the surfaces should be added to the air volumes

3807 # instead of the terminal, winding, and contact layer volumes.

3808 for key, dimTagsList in surfacesDict.items():

3809 if key is self.individualWinding[self.pancakeIndex]:

3810 dimTagsObjects = [self.topAirPancakeWindingExtursionVolume] * len(

3811 dimTagsList

3812 )

3813 elif key is self.individualContactLayer[self.pancakeIndex]:

3814 dimTagsObjects = [

3815 self.topAirPancakeContactLayerExtursionVolume

3816 ] * len(dimTagsList)

3817 elif (

3818 key is terminalDimTagsObject

3819 or key is self.airShellVolume

3820 or key is self.airShellVolumePart1

3821 or key is self.airShellVolumePart2

3822 or key is self.airShellVolumePart3

3823 or key is self.airShellVolumePart4

3824 or key is self.outerAirTubeVolume

3825 or key is self.centerAirCylinderVolume

3826 ):

3827 dimTagsObjects = [key] * len(dimTagsList)

3828 else:

3829 # key is self.outerTerminalTouchingVolume

3830 # or key is self.innerTerminalTouchingVolume

3831 # or key is (other terminal's tube volume)

3832 dimTagsObjects = [self.topAirTerminalsExtrusionVolume] * len(

3833 dimTagsList

3834 )

3835 if (

3836 key is self.innerTerminalTubeVolume

3837 or key is self.outerTerminalTubeVolume

3838 ):

3839 dimTagsObjects = [

3840 self.topAirTubeTerminalsExtrusionVolume

3841 ] * len(dimTagsList)

3842

3843 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(

3844 dimTagsList

3845 )

3846 listOfDimTags = listOfDimTags + dimTagsList

3847 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects

3848 elif bottomPart:

3849 # Then in this case, most of the surfaces should be added to the air volumes

3850 # instead of the terminal, winding, and contact layer volumes.

3851 for key, dimTagsList in surfacesDict.items():

3852 if key is self.individualWinding[self.pancakeIndex]:

3853 dimTagsObjects = [

3854 self.bottomAirPancakeWindingExtursionVolume

3855 ] * len(dimTagsList)

3856 elif key is self.individualContactLayer[self.pancakeIndex]:

3857 dimTagsObjects = [

3858 self.bottomAirPancakeContactLayerExtursionVolume

3859 ] * len(dimTagsList)

3860 elif (

3861 key is terminalDimTagsObject

3862 or key is self.airShellVolume

3863 or key is self.airShellVolumePart1

3864 or key is self.airShellVolumePart2

3865 or key is self.airShellVolumePart3

3866 or key is self.airShellVolumePart4

3867 or key is self.outerAirTubeVolume

3868 or key is self.centerAirCylinderVolume

3869 ):

3870 dimTagsObjects = [key] * len(dimTagsList)

3871 else:

3872 # key is self.outerTerminalTouchingVolume

3873 # or key is self.innerTerminalTouchingVolume

3874 # or key is (other terminal's tube volume)

3875 dimTagsObjects = [self.bottomAirTerminalsExtrusionVolume] * len(

3876 dimTagsList

3877 )

3878 if (

3879 key is self.innerTerminalTubeVolume

3880 or key is self.outerTerminalTubeVolume

3881 ):

3882 dimTagsObjects = [

3883 self.bottomAirTubeTerminalsExtrusionVolume

3884 ] * len(dimTagsList)

3885

3886 listOfDimTagsForTopSurfaces = listOfDimTagsForTopSurfaces + [key] * len(

3887 dimTagsList

3888 )

3889 listOfDimTags = listOfDimTags + dimTagsList

3890 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects

3891 else:

3892 for key, dimTagsList in surfacesDict.items():

3893 if (

3894 key is self.outerAirTubeVolume

3895 or key is self.centerAirCylinderVolume

3896 or key is self.airShellVolume

3897 or key is self.airShellVolumePart1

3898 or key is self.airShellVolumePart2

3899 or key is self.airShellVolumePart3

3900 or key is self.airShellVolumePart4

3901 or key is terminalDimTagsObject

3902 ):

3903 dimTagsObjects = [key] * len(dimTagsList)

3904

3905 listOfDimTags = listOfDimTags + dimTagsList

3906 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects

3907

3908 listOfDimTagsForTopSurfaces = listOfDimTagsObjects

3909

3910 extrusionResult = dimTags()

3911 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, tZ))

3912

3913 # Add the created volumes to the corresponding dimTags objects:

3914 volumeDimTags = extrusionResult.getDimTags(3)

3915 for i, volumeDimTag in enumerate(volumeDimTags):

3916 listOfDimTagsObjects[i].add(volumeDimTag)

3917

3918 if firstTerminal:

3919 self.firstTerminalVolume.add(terminalDimTagsObject.getDimTags(3))

3920 elif lastTerminal:

3921 self.lastTerminalVolume.add(terminalDimTagsObject.getDimTags(3))

3922

3923 topSurfacesDimTags = extrusionResult.getExtrusionTop()

3924 topSurfacesDict = {}

3925 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):

3926 if listOfDimTagsObjects[i] in topSurfacesDict:

3927 topSurfacesDict[listOfDimTagsForTopSurfaces[i]].append(topSurfaceDimTag)

3928 else:

3929 topSurfacesDict[listOfDimTagsForTopSurfaces[i]] = [topSurfaceDimTag]

3930

3931 return topSurfacesDict

3932

3933 def extrudeWindingPart(self):

3934 """

3935 Extrudes all the fundamental surfaces of the pancake coil by self.geo.winding.height and

3936 returns the next connection terminal's top surface dimTag, and other air dimTags

3937 in a dictionary so that they can be further extruded.

3938

3939 :return: dictionary of top surfaces where the keys are the dimTags objects and

3940 the values are the dimTags of the surfaces

3941 :rtype: dict[dimTags, list[tuple[int, int]]]

3942 """

3943 # Create the list of surfaces to be extruded:

3944 listOfDimTags = []

3945 listOfDimTagsObjects = []

3946 for key, dimTagsList in self.fundamentalSurfaces.items():

3947 dimTagsObjects = [key] * len(dimTagsList)

3948

3949 listOfDimTags = listOfDimTags + dimTagsList

3950 listOfDimTagsObjects = listOfDimTagsObjects + dimTagsObjects

3951

3952 # Extrude the fundamental surfaces:

3953 extrusionResult = dimTags()

3954 extrusionResult.add(gmsh.model.occ.extrude(listOfDimTags, 0, 0, self.geo.winding.height))

3955

3956 # Add the created volumes to the corresponding dimTags objects:

3957 volumes = extrusionResult.getDimTags(3)

3958 for i, volumeDimTag in enumerate(volumes):

3959 listOfDimTagsObjects[i].add(volumeDimTag)

3960

3961 if self.pancakeIndex == 0:

3962 # Note the first pancake (sometimes useful for creating regions in the

3963 # meshing part):

3964 for i, volumeDimTag in enumerate(volumes):

3965 if listOfDimTagsObjects[i].parentName == self.geo.terminals.outer.name:

3966 self.firstTerminalVolume.add(volumeDimTag)

3967

3968 # Not elif! Because the first pancake coil is also the last pancake coil if

3969 # there is only one pancake coil.

3970 if self.pancakeIndex == self.geo.numberOfPancakes - 1:

3971 # Note the last pancake (sometimes useful for creating regions in the

3972 # meshing part):

3973 for i, volumeDimTag in enumerate(volumes):

3974 if (

3975 self.pancakeIndex % 2 == 1

3976 and listOfDimTagsObjects[i].parentName == self.geo.terminals.outer.name

3977 ):

3978 self.lastTerminalVolume.add(volumeDimTag)

3979 elif (

3980 self.pancakeIndex % 2 == 0

3981 and listOfDimTagsObjects[i].parentName == self.geo.terminals.inner.name

3982 ):

3983 self.lastTerminalVolume.add(volumeDimTag)

3984

3985 # Return the top surfaces:

3986 # Add the created top surfaces to a new dictionary:

3987 topSurfacesDimTags = extrusionResult.getExtrusionTop()

3988 topSurfaces = {}

3989 for i, topSurfaceDimTag in enumerate(topSurfacesDimTags):

3990 if listOfDimTagsObjects[i] in topSurfaces:

3991 topSurfaces[listOfDimTagsObjects[i]].append(topSurfaceDimTag)

3992 else:

3993 topSurfaces[listOfDimTagsObjects[i]] = [topSurfaceDimTag]

3994

3995 return topSurfaces

3996

3997 def generateGapAirWithFragment(

3998 self, gapAirSurfacesDimTags: List[List[Tuple[int, int]]]

3999 ):

4000 """

4001 A class to fill the gap between the multiple pancake coils with air. First, it

4002 creates a dummy cylinder with the same radius as the outer terminal's outer

4003 radius. Then using gapAirSurfacesDimTags, gmsh.model.occ.fragment() operation is

4004 applied to the dummy cylinder volume in a for loop to create the gap air

4005 volumes. After each fragment operation, one of the volumes created is removed

4006 because it is the solid volume which is the combination of windings,

4007 contact layers, and terminals. In the end, dummy cylinder is removed as well.

4008

4009

4010 WARNING:

4011 Currently, this method doesn't work.

4012

4013 :param geometry: geometry information

4014 :param pancakeCoils: pancakeCoilsWithAir object

4015 :type pancakeCoils: pancakeCoilsWithAir

4016 """

4017 logger.info("Generating gap air has been started.")

4018 start_time = timeit.default_timer()

4019

4020 # Create the dummy air volume:

4021 dummyAir = gmsh.model.occ.addCylinder(

4022 0,

4023 0,

4024 -self.geo.air.h / 2,

4025 0,

4026 0,

4027 self.geo.air.h,

4028 self.geo.terminals.outer.r,

4029 )

4030

4031 toBeDeletedDimTags = []

4032 gapAirVolumesCurrentDimTags = []

4033 for i in range(len(gapAirSurfacesDimTags)):

4034 # Get the outer surfaces of the pancake coils for cutting the pancake coils

4035 # from the dummy air. The outer surfaces are used instead of pancake volumes

4036 # to reduce the amount of work for gmsh. It makes it significantly faster.

4037 # if len(gapAirSurfacesDimTags[i]) !=12:

4038 fragmentResults = gmsh.model.occ.fragment(

4039 [(3, dummyAir)],

4040 gapAirSurfacesDimTags[i],

4041 removeObject=False,

4042 removeTool=False,

4043 )

4044 fragmentVolumeResultsDimTags = fragmentResults[1][0]

4045 toBeDeletedDimTags.append(fragmentVolumeResultsDimTags[0])

4046 gapAirVolumesCurrentDimTags.append(fragmentVolumeResultsDimTags[1])

4047

4048 toBeDeletedDimTags.append((3, dummyAir))

4049 # Fragmnet operation both creates the air volume and solid pancake coils volume

4050 # because the surfaces are used for cutting. Therefore, the solid pancake coils

4051 # volume should be removed from the fragment results:

4052 gmsh.model.occ.remove(toBeDeletedDimTags)

4053

4054 # Add results to the air volume storage. After the geometry is saves as a .brep

4055 # file, and loaded back, the gaps between the tags are avoided by moving the

4056 # the other tags. Therefore, this is how the tags are stored:

4057 toBeDeletedTags = [dimTag[1] for dimTag in toBeDeletedDimTags]

4058 volumeTagsStart = min(toBeDeletedTags)

4059 numberOfGapAirVolumes = len(gapAirVolumesCurrentDimTags)

4060 gapAirVolumesToBeSaved = [

4061 (3, volumeTagsStart + i) for i in range(numberOfGapAirVolumes)

4062 ]

4063

4064 # For debugging purposes, physical groups are being created in the geometry

4065 # generation process as well. Normally, it us done during meshing because BREP

4066 # files cannot store physical groups. Since the tags above (airVolumes) will be

4067 # valid after only saving the geometry as a BREP file and loading it back, the

4068 # current tags are given to the airVolume.add() method as well. This is done to

4069 # be able to create the correct physical group.

4070 self.gapAirVolume.add(

4071 dimTagsList=gapAirVolumesToBeSaved,

4072 dimTagsListForPG=gapAirVolumesCurrentDimTags,

4073 )

4074

4075 logger.info(

4076 "Generating gap air has been finished in"

4077 f" {timeit.default_timer() - start_time:.2f} s."

4078 )

4079

4080

4081class Geometry(Base):

4082 """

4083 Main geometry class for Pancake3D.

4084

4085 :param fdm: FiQuS data model

4086 :param geom_folder: folder where the geometry files are saved

4087 :type geom_folder: str

4088 :param mesh_folder: folder where the mesh files are saved

4089 :type mesh_folder: str

4090 :param solution_folder: folder where the solution files are saved

4091 :type solution_folder: str

4092 """

4093

4094 def __init__(

4095 self,

4096 fdm,

4097 geom_folder,

4098 mesh_folder,

4099 solution_folder,

4100 ) -> None:

4101 super().__init__(fdm, geom_folder, mesh_folder, solution_folder)

4102

4103 # Clear if there is any existing dimTags storage:

4104 dimTagsStorage.clear()

4105

4106 # Start GMSH:

4107 self.gu = GmshUtils(self.geom_folder)

4108 self.gu.initialize(verbosity_Gmsh=fdm.run.verbosity_Gmsh)

4109

4110 # To speed up the GUI:

4111 gmsh.option.setNumber("Geometry.NumSubEdges", 10)

4112

4113 # To see the surfaces in a better way in GUI:

4114 gmsh.option.setNumber("Geometry.SurfaceType", 1)

4115

4116 # To avoid any unwanted modifications to the geometry, the automatic fixing of

4117 # the geometry is disabled:

4118 gmsh.option.setNumber("Geometry.OCCAutoFix", 0)

4119

4120 # Set the tolerance:

4121 if self.geo.dimensionTolerance < gmsh.option.getNumber("Geometry.Tolerance"):

4122 gmsh.option.setNumber("Geometry.Tolerance", self.geo.dimensionTolerance)

4123

4124 gmsh.option.setNumber("Geometry.ToleranceBoolean", self.geo.dimensionTolerance)

4125

4126 spiralCurve.sectionsPerTurn = self.geo.winding.spt

4127 spiralCurve.curvesPerTurn = self.geo.winding.numberOfVolumesPerTurn

4128

4129 def generate_geometry(self):

4130 """

4131 Generates geometry and saves it as a .brep file.

4132

4133

4134 """

4135 logger.info(

4136 f"Generating Pancake3D geometry ({self.brep_file}) has been started."

4137 )

4138 start_time = timeit.default_timer()

4139

4140 self.pancakeCoil = pancakeCoilsWithAir(self.geo, self.mesh)

4141

4142 gmsh.model.occ.synchronize()

4143 gmsh.write(self.brep_file)

4144

4145 logger.info(

4146 f"Generating Pancake3D geometry ({self.brep_file}) has been finished in"

4147 f" {timeit.default_timer() - start_time:.2f} s."

4148 )

4149

4150 def load_geometry(self):

4151 """

4152 Loads geometry from .brep file.

4153 """

4154 logger.info("Loading Pancake3D geometry has been started.")

4155 start_time = timeit.default_timer()

4156

4157 previousGeo = FilesAndFolders.read_data_from_yaml(

4158 self.geometry_data_file, Pancake3DGeometry

4159 )

4160

4161 if previousGeo.dict() != self.geo.dict():

4162 raise ValueError(

4163 "Geometry data has been changed. Please regenerate the geometry or load"

4164 " the previous geometry data."

4165 )

4166

4167 gmsh.clear()

4168 gmsh.model.occ.importShapes(self.brep_file, format="brep")

4169 gmsh.model.occ.synchronize()

4170

4171 logger.info(

4172 "Loading Pancake3D geometry has been finished in"

4173 f" {timeit.default_timer() - start_time:.2f} s."

4174 )

4175

4176 def generate_vi_file(self):

4177 """

4178 Generates volume information file. Volume information file stores dimTags of all

4179 the stored volumes in the geometry. Without this file, regions couldn't be

4180 created, meaning that finite element simulation cannot be done.

4181

4182 The file extension is custom because users are not supposed to edit or open this

4183 file, and it makes it intuitively clear that it is a volume information file.

4184 """

4185 logger.info(

4186 f"Generating volume information file ({self.vi_file}) has been started."

4187 )

4188 start_time = timeit.default_timer()

4189

4190 dimTagsDict = dimTagsStorage.getDimTagsDict()

4191 json.dump(

4192 dimTagsDict,

4193 open(self.vi_file, "w"),

4194 )

4195

4196 logger.info(

4197 "Generating volume information file has been finished in"

4198 f" {timeit.default_timer() - start_time:.2f} s."

4199 )

4200

4201 if self.geo_gui:

4202 self.generate_physical_groups()

4203 self.gu.launch_interactive_GUI()

4204 else:

4205 gmsh.clear()

4206 gmsh.finalize()

4207

4208 @staticmethod

4209 def generate_physical_groups():

4210 """

4211 Generates physical groups. Physical groups are not saved in the BREP file but

4212 it can be useful for debugging purposes.

4213

4214

4215 """

4216 gmsh.model.occ.synchronize()

4217

4218 dimTagsDict = dimTagsStorage.getDimTagsDict(forPhysicalGroups=True)

4219

4220 for key, value in dimTagsDict.items():

4221 tags = [dimTag[1] for dimTag in value]

4222 gmsh.model.addPhysicalGroup(

4223 3,

4224 tags,

4225 name=key,

4226 )

4227

4228 @staticmethod

4229 def remove_all_duplicates():

4230 """

4231 Removes all the duplicates and then prints the entities that are created or

4232 removed during the operation. It prints the line number where the function is

4233 called as well. This function is helpful for debugging. Finding duplicates means

4234 there is a problem in geometry creation logic, and the meshes will not be

4235 conformal. It shouldn't be used in the final version of the code since removing

4236 duplicates is computationally expensive, and there shouldn't be duplicates at

4237 all.

4238

4239 WARNING:

4240 This function currently does not work properly. It is not recommended to use

4241 right now. It finds duplicates even if there are no duplicates (topology

4242 problems).

4243 """

4244

4245 logger.info(f"Removing all the duplicates has been started.")

4246 start_time = timeit.default_timer()

4247

4248 gmsh.model.occ.synchronize()

4249 oldEntities = []

4250 oldEntities.extend(gmsh.model.getEntities(3))

4251 oldEntities.extend(gmsh.model.getEntities(2))

4252 oldEntities.extend(gmsh.model.getEntities(1))

4253 oldEntities.extend(gmsh.model.getEntities(0))

4254 oldEntities = set(oldEntities)

4255

4256 gmsh.model.occ.removeAllDuplicates()

4257

4258 gmsh.model.occ.synchronize()

4259 newEntities = []

4260 newEntities.extend(gmsh.model.getEntities(3))

4261 newEntities.extend(gmsh.model.getEntities(2))

4262 newEntities.extend(gmsh.model.getEntities(1))

4263 newEntities.extend(gmsh.model.getEntities(0))

4264 newEntities = set(newEntities)

4265 NewlyCreated = newEntities - oldEntities

4266 Removed = oldEntities - newEntities

4267

4268 frameinfo = getframeinfo(currentframe().f_back)

4269

4270 if len(NewlyCreated) > 0 or len(Removed) > 0:

4271 logger.warning(f"Duplicates found! Line: {frameinfo.lineno}")

4272 logger.warning(f"{len(NewlyCreated)}NewlyCreated = {list(NewlyCreated)}")

4273 logger.warning(f"{len(Removed)}Removed = {list(Removed)}")

4274 else:

4275 logger.info(f"No duplicates found! Line: {frameinfo.lineno}")

4276

4277 logger.info(

4278 "Removing all the duplicates has been finished in"

4279 f" {timeit.default_timer() - start_time:.2f} s."

4280 )