Merge pull request #13 from OpenSEMBA/feature/unshielded-multiwires

Changes near vacuum to a rectangle

Author: lmdiazangulo

src/BoundingBox.py

@@ -4,28 +4,42 @@
 
 class BoundingBox():
     edges: Dict[str, float]
-    def __init__(self, listOfCoordinates:Tuple[int, int, int, int, int, int]):
+    def __init__(self, listOfCoordinates:Tuple[float,float,float,float,float,float]):
         self.edges = {
-                'xmin': listOfCoordinates[0],
-                'ymin': listOfCoordinates[1],
-                'zmin': listOfCoordinates[2],
-                'xmax': listOfCoordinates[3],
-                'ymax': listOfCoordinates[4],
-                'zmax': listOfCoordinates[5],
+                'XMin': listOfCoordinates[0],
+                'YMin': listOfCoordinates[1],
+                'ZMin': listOfCoordinates[2],
+                'XMax': listOfCoordinates[3],
+                'YMax': listOfCoordinates[4],
+                'ZMax': listOfCoordinates[5],
             }
 
+    def getOrigin(self) -> Tuple[float,float,float]:
+        return (
+            self.edges['XMin'],
+            self.edges['YMin'],
+            self.edges['ZMin']
+        )
+
     def getCenter(self) -> Tuple[float,float,float]:
         return (
-            (self.edges['xmax'] + self.edges['xmin']) / 2,
-            (self.edges['ymax'] + self.edges['ymin']) / 2,
-            (self.edges['zmax'] + self.edges['zmin']) / 2
+            (self.edges['XMax'] + self.edges['XMin']) / 2,
+            (self.edges['YMax'] + self.edges['YMin']) / 2,
+            (self.edges['ZMax'] + self.edges['ZMin']) / 2
         )
     def getDiagonal(self) -> float:
-        dx = self.edges['xmax'] - self.edges['xmin']
-        dy = self.edges['ymax'] - self.edges['ymin']
-        dz = self.edges['zmax'] - self.edges['zmin']
+        dx = self.edges['XMax'] - self.edges['XMin']
+        dy = self.edges['YMax'] - self.edges['YMin']
+        dz = self.edges['ZMax'] - self.edges['ZMin']
         return (dx**2 + dy**2 + dz**2) ** 0.5
 
+    def getLengths(self) -> Tuple[float, float, float]:
+        return (
+            self.edges['XMax'] - self.edges['XMin'],
+            self.edges['YMax'] - self.edges['YMin'],
+            self.edges['ZMax'] - self.edges['ZMin']
+        )
+    
     @staticmethod
     def _getBoundingBox(element:Tuple[int,int]) -> 'BoundingBox':
         boundingBox:BoundingBox = BoundingBox(gmsh.model.occ.get_bounding_box(*element))
@@ -39,25 +53,25 @@ def getBoundingBoxFromGroup(elements:List[Tuple[int,int]]) -> 'BoundingBox':
 
         if len(boundingBoxs) != 0:
             edges: Dict[str, List[float]] = {
-                'xmin': [],
-                'ymin': [],
-                'zmin': [],
-                'xmax': [],
-                'ymax': [],
-                'zmax': [],
+                'XMin': [],
+                'YMin': [],
+                'ZMin': [],
+                'XMax': [],
+                'YMax': [],
+                'ZMax': [],
             }
             for boundingBox in boundingBoxs:
                 for key in edges.keys():
                     edges[key].append(boundingBox.edges[key])
 
             return BoundingBox(
                 (
-                    utils.getMinFromList(edges['xmin']),
-                    utils.getMinFromList(edges['ymin']),
-                    utils.getMinFromList(edges['zmin']),
-                    utils.getMaxFromList(edges['xmax']),
-                    utils.getMaxFromList(edges['ymax']),
-                    utils.getMaxFromList(edges['zmax'])
+                    utils.getMinFromList(edges['XMin']),
+                    utils.getMinFromList(edges['YMin']),
+                    utils.getMinFromList(edges['ZMin']),
+                    utils.getMaxFromList(edges['XMax']),
+                    utils.getMaxFromList(edges['YMax']),
+                    utils.getMaxFromList(edges['ZMax'])
                 )
             )
         else:

src/ShapesClassification.py

@@ -3,6 +3,7 @@
 import gmsh
 from .BoundingBox import BoundingBox
 from itertools import chain
+import numpy as np
 
 class ShapesClassification:
     isOpenCase:bool
@@ -87,12 +88,12 @@ def buildVacuumDomain(self):
         if self.isOpenCase and len(self.open) == 0:
             self.vacuum = self._buildDefaultVacuumDomain()
         elif self.isOpenCase and len(self.open) > 0:
-            self.vacuum = self._buildOpenVacuumDomain()
+            self.vacuum = self._buildVacuumDomainFromOpenBoundary()
         else:
             self.vacuum = self._buildClosedVacuumDomain()
         return self.vacuum
 
-    def _buildOpenVacuumDomain(self) -> Dict[int, List[int]]:
+    def _buildVacuumDomainFromOpenBoundary(self) -> Dict[int, List[int]]:
         dom = self.open[0]
 
         surfsToRemove = []
@@ -131,14 +132,18 @@ def _buildDefaultVacuumDomain(self):
             nonVacuumSurfaces.extend(surf)
 
         boundingBox = BoundingBox.getBoundingBoxFromGroup(nonVacuumSurfaces)
-        boundingBoxCenter = boundingBox.getCenter()
-        boundingBoxDiagonal = boundingBox.getDiagonal()
 
-        nearVacuumDiameter = boundingBoxDiagonal
-        nearVacuum = [(2, gmsh.model.occ.addDisk(*boundingBoxCenter, nearVacuumDiameter, nearVacuumDiameter))]
-
-        farVacuumDiameter = 4*boundingBoxDiagonal
-        farVacuum = [(2, gmsh.model.occ.addDisk(*boundingBoxCenter, farVacuumDiameter, farVacuumDiameter))]
+    
+        bbMaxLength = np.max(boundingBox.getLengths())
+        nVOrigin = tuple(
+            np.subtract(boundingBox.getCenter(), 
+                        (bbMaxLength, bbMaxLength, 0.0)))
+        nearVacuum = [(2, gmsh.model.occ.addRectangle(*nVOrigin, *(bbMaxLength*2.0,)*2))]
+
+        farVacuumDiameter = 4.0 * boundingBox.getDiagonal()
+        farVacuum = [(2, gmsh.model.occ.addDisk(
+            *boundingBox.getCenter(), 
+            farVacuumDiameter, farVacuumDiameter))]
 
         gmsh.model.occ.synchronize()
 
@@ -148,8 +153,16 @@ def _buildDefaultVacuumDomain(self):
         nearVacuum = gmsh.model.occ.cut(
             nearVacuum, nonVacuumSurfaces, removeObject=True, removeTool=False)[0]
 
-        self.open = dict([[0, gmsh.model.getBoundary(farVacuum)]])
+        # -- Set mesh size for near vacuum region
+        bb = BoundingBox(
+            gmsh.model.getBoundingBox(2, nearVacuum[0][1]))
+        minSide = np.min(np.array([bb.getLengths()[0], bb.getLengths()[1]]))
 
+        innerRegion = gmsh.model.getBoundary(nearVacuum, recursive=True)
+        gmsh.model.mesh.setSize(innerRegion, minSide / 20)
+        
+        
+        self.open = dict([[0, gmsh.model.getBoundary(farVacuum)]])
         gmsh.model.occ.synchronize()
 
         return dict([[0, nearVacuum], [1, farVacuum]])

src/mesher.py

@@ -5,16 +5,18 @@
 
 from src.AreaExporterService import AreaExporterService
 from .ShapesClassification import ShapesClassification
+from .BoundingBox import BoundingBox
+import numpy as np
 
 class Mesher():
     DEFAULT_MESHING_OPTIONS = {
 
-        "Mesh.MshFileVersion": 2.2,   # Mandatory for MFEM compatibility
+        "Mesh.MshFileVersion": 2.2,   # Required for MFEM compatibility
         "Mesh.MeshSizeFromCurvature": 50,
         "Mesh.ElementOrder": 3,
         "Mesh.ScalingFactor": 1e-3,
         "Mesh.SurfaceFaces": 1,
-        "Mesh.MeshSizeMax": 50,
+        "Mesh.MeshSizeMax": 40,
 
         "General.DrawBoundingBoxes": 1,
         "General.Axes": 1,
@@ -32,7 +34,7 @@ def runFromInput(self, inputFile, runGui=False):
         self.meshFromStep(inputFile, caseName, self.DEFAULT_MESHING_OPTIONS)
         self.exportGeometryAreas(caseName)
         gmsh.write(caseName + '.msh')
-        gmsh.write(caseName + '.vtk')
+        gmsh.write(caseName + '.vtk') # vtk export is just for debugging. 
         if runGui:
             gmsh.fltk.run()
 
@@ -64,6 +66,8 @@ def meshFromStep(self, inputFile: str, caseName: str, meshingOptions=None):
         for [opt, val] in meshingOptions.items():
             gmsh.option.setNumber(opt, val)
 
+        # --- Mesh generation ---
+        
         gmsh.model.mesh.generate(2)
 
     def exportGeometryAreas(self, caseName:str):

test/test_BoundingBox.py

@@ -10,24 +10,24 @@ class testBoundingBox(unittest.TestCase):
     def testCanDefineBoundingBox(self):
         inputExample = (-1, 2.0, 3, 5, 6, 7)
         expectedBoundingBoxEdges: Dict[str, float] = {
-                    'xmin': -1,
-                    'ymin': 2.0,
-                    'zmin': 3,
-                    'xmax': 5,
-                    'ymax': 6,
-                    'zmax': 7,
+                    'XMin': -1,
+                    'YMin': 2.0,
+                    'ZMin': 3,
+                    'XMax': 5,
+                    'YMax': 6,
+                    'ZMax': 7,
                 }
         boundingBox = BoundingBox(inputExample)
         self.assertDictEqual(boundingBox.edges, expectedBoundingBoxEdges)
 
     def testCanGetBoundingBoxFromSurface(self):
         expectedBoundingBoxEdges: Dict[str, float] = {
-            'xmin': -10,
-            'ymin': -10,
-            'zmin': 0,
-            'xmax': 10,
-            'ymax': 10,
-            'zmax': 0,
+            'XMin': -10,
+            'YMin': -10,
+            'ZMin': 0,
+            'XMax': 10,
+            'YMax': 10,
+            'ZMax': 0,
         }
 
         gmsh.initialize()
@@ -39,12 +39,12 @@ def testCanGetBoundingBoxFromSurface(self):
 
     def testCanGetBoundingBoxFromListOfSurfaces(self):
         expectedBoundingBoxEdges: Dict[str, float] = {
-            'xmin': -10,
-            'ymin': -10,
-            'zmin': 0,
-            'xmax': 35,
-            'ymax': 40,
-            'zmax': 0,
+            'XMin': -10,
+            'YMin': -10,
+            'ZMin': 0,
+            'XMax': 35,
+            'YMax': 40,
+            'ZMax': 0,
         }
 
         gmsh.initialize()

test/test_mesher.py

@@ -223,6 +223,8 @@ def test_unshielded_multiwire(self):
         caseName = 'unshielded_multiwire'
         Mesher().meshFromStep(self.inputFileFromCaseName(caseName), caseName)
 
+        gmsh.write(caseName + '.vtk')
+
         pGs = gmsh.model.getPhysicalGroups()
         pGNames = [gmsh.model.getPhysicalName(*pG) for pG in pGs]
         expectedNames = ['Conductor_0', 'Conductor_1',  'Dielectric_1', 
@@ -236,5 +238,27 @@ def test_unshielded_multiwire(self):
         for idx, name in enumerate(expectedNames):
             self.assertEqual(self.countEntitiesInPhysicalGroupWithName(name), expectedEntities[idx], name)
 
+    def test_lansink2024_single_wire_multipolar(self):
+        caseName = 'lansink2024_single_wire_multipolar'
+        Mesher().meshFromStep(self.inputFileFromCaseName(caseName), caseName)
+        
+        gmsh.write(caseName + '.msh')
+        gmsh.write(caseName + '.vtk')
+
+        pGs = gmsh.model.getPhysicalGroups()
+        pGNames = [gmsh.model.getPhysicalName(*pG) for pG in pGs]
+        expectedNames = ['Conductor_0', 
+                         'Dielectric_0', 
+                         'OpenBoundary_0',
+                         'Vacuum_0', 'Vacuum_1']
+        expectedEntities = [1, 1, 1,
+                            1,
+                            1, 1]
+        self.assertEqual(sorted(pGNames), sorted(expectedNames))
+
+        for idx, name in enumerate(expectedNames):
+            self.assertEqual(self.countEntitiesInPhysicalGroupWithName(name), expectedEntities[idx], name)
+
+
 if __name__ == '__main__':
     unittest.main()