Merge pull request #14 from OpenSEMBA/dev

Dev

Author: lmdiazangulo

.github/workflows/tests.yml

@@ -7,9 +7,12 @@ on:
   push:
     branches:
       - master
+      - dev
+
   pull_request:
     branches:
       - master
+      - dev
 
 env:
   STEP2GMSH_TOP_DIR: step2gmsh

.gitignore

@@ -2,3 +2,7 @@
 *.pyc
 *.msh
 .vscode/
+tmpFolder/
+venv/
+gmshDoc/
+*.vtk

README.md

@@ -26,7 +26,7 @@ The tested input step files have been generated with [FreeCAD](https://www.freec
 
 - A layer named `Conductor_N` with `N` being an integer represents a perfect conductor. `Conductor_0` is a special case of which represents the ground and defines the global domain. For layers named `Conductor_N` with `N` different to zero their areas will be substracted from the computational domain and removed.
 - Layers named as `Dielectric_N` are used to identify regions which will have a material assigned.
-- Open and semi-open problems can be defined using a single layer called `OpenRegion`.
+- Open and semi-open problems can be defined using a single layer called `OpenBoundary`.
 
 Below is shown an example of a closed case with 6 conductors and 5 dielectrics, the external boundary corresponds to `Conductor_0`. The case is modeled with FreeCAD and can be found in the `testData/five_wires` folder together with the exported as a step file. The resulting mesh after applying `step2gmsh` is shown below.
 

__init__.py

@@ -1,2 +1,3 @@
 gmsh >= 4.11
-pytest >= 7.3
+pytest >= 7.3
+numpy

requirements.txt

@@ -0,0 +1,37 @@
+import json
+import gmsh
+from typing import Dict, List
+import numpy as np
+class AreaExporterService:
+    _EMPTY_NAME_CASE = ""
+    computedAreas:Dict[str,List]
+    geometry: Dict
+    def __init__(self):
+        self.computedAreas = {
+            "geometries": []
+        }
+    
+    def addComputedArea(self, geometry:str, area:float):
+        geometry:Dict ={
+            "geometry": geometry,
+            "area": area,
+        }
+        self.computedAreas['geometries'].append(geometry)
+
+    def addPhysicalModelOfDimension(self, dimension=2):
+        physicalGroups = gmsh.model.getPhysicalGroups(dimension)
+        for physicalGroup in physicalGroups:
+            entityTags = gmsh.model.getEntitiesForPhysicalGroup(*physicalGroup)
+            geometryName = gmsh.model.getPhysicalName(*physicalGroup)
+            for tag in entityTags:
+                if dimension == 1:
+                    rad = gmsh.model.occ.getMass(dimension, tag) / (2*np.pi)
+                    area = rad*rad*np.pi
+                if dimension == 2:
+                    area = gmsh.model.occ.getMass(dimension, tag)
+                if geometryName != AreaExporterService._EMPTY_NAME_CASE:
+                    self.addComputedArea(geometryName, area)
+
+    def exportToJson(self, exportFileName:str):
+        with open(exportFileName + ".areas.json", 'w') as f:
+            json.dump(self.computedAreas, f, indent=3)

src/AreaExporterService.py

@@ -0,0 +1,78 @@
+from typing import Tuple, List, Dict
+import gmsh
+from . import utils
+
+class BoundingBox():
+    edges: Dict[str, float]
+    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],
+            }
+        
+    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
+        )
+    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']
+        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))
+        return boundingBox
+
+    @staticmethod
+    def getBoundingBoxFromGroup(elements:List[Tuple[int,int]]) -> 'BoundingBox':
+        boundingBoxs:List[BoundingBox] = []
+        for element in elements:
+            boundingBoxs.append(BoundingBox._getBoundingBox(element))
+        
+        if len(boundingBoxs) != 0:
+            edges: Dict[str, List[float]] = {
+                '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'])
+                )
+            )
+        else:
+            return BoundingBox((0,0,0,0,0,0))

src/BoundingBox.py

@@ -0,0 +1,170 @@
+from typing import Any, Tuple, List, Dict
+
+import gmsh
+from .BoundingBox import BoundingBox
+from itertools import chain
+import numpy as np
+
+class ShapesClassification:
+    isOpenCase:bool
+
+
+    def __init__(self, shapes):
+        gmsh.model.occ.synchronize()
+
+        self.allShapes = shapes
+        self.pecs = self.get_surfaces_with_label(shapes, "Conductor_")
+        self.dielectrics = self.get_surfaces_with_label(shapes, "Dielectric_")
+        self.open = self.get_surfaces_with_label(shapes, "OpenBoundary_")
+        self.vacuum = dict()
+
+        self.isOpenCase = self.isOpenProblem()
+
+        if len(self.open) > 1:
+            raise ValueError("Only one open region is allowed.")
+
+    @staticmethod
+    def getNumberFromName(entity_name: str, label: str):
+        ini = entity_name.rindex(label) + len(label)
+        num = int(entity_name[ini:])
+        return num
+
+    @staticmethod
+    def get_surfaces_with_label(entity_tags, label: str):
+        surfaces = dict()
+        for s in entity_tags:
+            name = gmsh.model.get_entity_name(*s)
+            if s[0] != 2 or label not in name:
+                continue
+            num = ShapesClassification.getNumberFromName(name, label)
+            surfaces[num] = [s]
+
+        return surfaces
+
+    def isOpenProblem(self):
+        elements = list(chain(self.pecs.values()))
+        for idx, element in enumerate(elements):
+            for otheridx, otherElement in enumerate(elements[idx+1:]):
+                if element != otherElement:
+                    intersect = gmsh.model.occ.intersect(
+                        element, 
+                        otherElement,
+                        removeObject=False,
+                        tag=300+otheridx,
+                        removeTool=False
+                    )[0]
+                    if intersect:
+                        return False   
+        return True
+    
+    def removeConductorsFromDielectrics(self):
+        for num, diel in self.dielectrics.items():
+            pec_surfs = []
+            for num2, pec_surf in self.pecs.items():
+                if num2 == 0 and not self.isOpenCase:
+                    continue
+                pec_surfs.extend(pec_surf)
+            self.dielectrics[num] = gmsh.model.occ.cut(diel, pec_surfs, removeTool=False)[0]
+
+        gmsh.model.occ.synchronize()
+
+    def ensureDielectricsDoNotOverlap(self):
+        for n1, diel1 in self.dielectrics.items():
+            others = list(
+                chain(
+                    *[x[1] for x in self.dielectrics.items() if x[0] != n1]
+                )
+            )
+
+            if len(others) == 0:
+                continue
+
+            self.dielectrics[n1] = gmsh.model.occ.cut(
+                self.dielectrics[n1], others, removeObject=True, removeTool=False)[0]
+
+        gmsh.model.occ.synchronize()
+
+    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._buildVacuumDomainFromOpenBoundary()
+        else:
+            self.vacuum = self._buildClosedVacuumDomain()
+        return self.vacuum
+    
+    def _buildVacuumDomainFromOpenBoundary(self) -> Dict[int, List[int]]:
+        dom = self.open[0]
+        
+        surfsToRemove = []
+        for num, surf in self.pecs.items():
+            surfsToRemove.extend(surf)
+
+        for _, surf in self.dielectrics.items():
+            surfsToRemove.extend(surf)
+
+        dom = gmsh.model.occ.cut(
+            dom, surfsToRemove, removeObject=False, removeTool=False)[0]
+        gmsh.model.occ.synchronize()
+
+        return dict([[0, dom]])
+    
+    def _buildClosedVacuumDomain(self) -> Tuple[int, int]:
+        dom = self.pecs[0]
+        surfsToRemove = []
+        for num, surf in self.pecs.items():
+            if num == 0:
+                continue
+            surfsToRemove.extend(surf)
+
+        for _, surf in self.dielectrics.items():
+            surfsToRemove.extend(surf)
+        dom = gmsh.model.occ.cut(
+            dom, surfsToRemove, removeObject=False, removeTool=False)[0]
+        gmsh.model.occ.synchronize()
+        return dict([[0, dom]])
+    
+    def _buildDefaultVacuumDomain(self):
+        nonVacuumSurfaces = []
+        for _, surf in self.pecs.items():
+            nonVacuumSurfaces.extend(surf)
+        for _, surf in self.dielectrics.items():
+            nonVacuumSurfaces.extend(surf)
+            
+        boundingBox = BoundingBox.getBoundingBoxFromGroup(nonVacuumSurfaces)
+
+    
+        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()
+
+        farVacuum = gmsh.model.occ.cut(
+            farVacuum, nearVacuum, removeObject=True, removeTool=False)[0]
+
+        nearVacuum = gmsh.model.occ.cut(
+            nearVacuum, nonVacuumSurfaces, removeObject=True, removeTool=False)[0]
+        
+        # -- 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]])
+    
+