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

1473 statements  

« prev     ^ index     » next       coverage.py v7.4.4, created at 2025-03-14 02:31 +0100

1import math 

2import logging 

3from enum import Enum 

4from inspect import currentframe, getframeinfo 

5from typing import List, Tuple, Dict 

6import operator 

7 

8import os 

9import json 

10import timeit 

11import numpy as np 

12import gmsh 

13 

14from fiqus.utils.Utils import GmshUtils 

15from fiqus.utils.Utils import FilesAndFolders 

16from fiqus.mains.MainPancake3D import Base 

17from fiqus.data.DataFiQuSPancake3D import Pancake3DGeometry 

18 

19logger = logging.getLogger(__name__) 

20 

21# coordinateList = [] 

22 

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) 

29 

30 return result 

31 

32 

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)))) 

57 

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)))) 

76 

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

78 

79 if findInnerOnes: 

80 goalDistance = min(distances) 

81 else: 

82 goalDistance = max(distances) 

83 

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) 

89 

90 return result 

91 

92 

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. 

97 

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 )