Merge pull request #72 from simbilod/multiple_physicals

Multiple physicals

meshwell/gmsh_entity.py

@@ -1,7 +1,7 @@
 import gmsh
 
 from typing import Any, Dict, Optional, List
-from pydantic import BaseModel
+from pydantic import BaseModel, validator
 from meshwell.resolution import ResolutionSpec
 
 
@@ -23,12 +23,19 @@ class GMSH_entity(BaseModel):
     gmsh_function_kwargs: Dict[str, Any]
     dimension: int
     model: Any
-    physical_name: Optional[str] = None
+    physical_name: Optional[str | tuple[str]] = None
     mesh_order: float | None = None
     mesh_bool: bool = True
     additive: bool = False
     resolutions: List[ResolutionSpec] | None = None
 
+    @validator("physical_name", pre=True, always=True)
+    def _set_physical_name(cls, physical_name: Optional[str | tuple[str]]):
+        if isinstance(physical_name, str):
+            return [physical_name]
+        else:
+            return physical_name
+
     def instanciate(self):
         """Returns dim tag from entity."""
         entity = self.gmsh_function(**self.gmsh_function_kwargs)

meshwell/labeledentity.py

@@ -10,7 +10,7 @@ class LabeledEntities(BaseModel):
     index: int
     model: Any
     dimtags: List[Tuple[int, int]]
-    physical_name: str
+    physical_name: str | tuple[str, ...]
     resolutions: List[ResolutionSpec] | None = None
     keep: bool
     boundaries: List[int] = []

meshwell/model.py

@@ -160,6 +160,9 @@ def mesh(
         interface_delimiter: str = "___",
         boundary_delimiter: str = "None",
         addition_delimiter: str = "+",
+        addition_intersection_physicals: bool = True,
+        addition_addition_physicals: bool = True,
+        addition_structural_physicals: bool = True,
         gmsh_version: Optional[float] = None,
         finalize: bool = True,
         reinitialize: bool = True,
@@ -180,6 +183,10 @@ def mesh(
             verbosity: GMSH stdout while meshing (True or False)
             interface_delimiter: string characters to use when naming interfaces between entities
             boundary_delimiter: string characters to use when defining an interface between an entity and nothing (simulation boundary)
+            addition_delimiter: string characters to use with addition_intersection_physicals
+            addition_addition_physicals: if True, the intersecting additive entities also carry the additive physical name
+            addition_structural_physicals: if True, the intersecting additive entities also carry the underlying entity names
+            addition_intersection_physicals: if True, the intersecting additive entities also carry a unique {underlying}{addition_delimiter}{additive} name
             gmsh_version: Gmsh mesh format version. For example, Palace requires an older version of 2.2,
                 see https://mfem.org/mesh-formats/#gmsh-mesh-formats.
             finalize: if True (default), finalizes the GMSH model after execution
@@ -243,14 +250,16 @@ def mesh(
             if progress_bars:
                 if physical_name:
                     enumerator.set_description(
-                        f"{physical_name:<30} - {'instanciate':<15}"
+                        f"{str(physical_name):<30} - {'instanciate':<15}"
                     )
             # First create the shape
             dimtags_out = entity_obj.instanciate()
 
             if progress_bars:
                 if physical_name:
-                    enumerator.set_description(f"{physical_name:<30} - {'dimtags':<15}")
+                    enumerator.set_description(
+                        f"{str(physical_name):<30} - {'dimtags':<15}"
+                    )
             # Parse dimension
             dim = validate_dimtags(dimtags_out)
             max_dim = max(dim, max_dim)
@@ -259,7 +268,7 @@ def mesh(
             if progress_bars:
                 if physical_name:
                     enumerator.set_description(
-                        f"{physical_name:<30} - {'entities':<15}"
+                        f"{str(physical_name):<30} - {'entities':<15}"
                     )
             # Assemble with other shapes
             current_entities = LabeledEntities(
@@ -272,7 +281,9 @@ def mesh(
             )
             if progress_bars:
                 if physical_name:
-                    enumerator.set_description(f"{physical_name:<30} - {'boolean':<15}")
+                    enumerator.set_description(
+                        f"{str(physical_name):<30} - {'boolean':<15}"
+                    )
             if index != 0:
                 cut = self.occ.cut(
                     current_entities.dimtags,
@@ -288,13 +299,13 @@ def mesh(
                 if progress_bars:
                     if physical_name:
                         enumerator.set_description(
-                            f"{physical_name:<30} - {'duplicates':<15}"
+                            f"{str(physical_name):<30} - {'duplicates':<15}"
                         )
                 self.occ.removeAllDuplicates()
                 if progress_bars:
                     if physical_name:
                         enumerator.set_description(
-                            f"{physical_name:<30} - {'sync':<15}"
+                            f"{str(physical_name):<30} - {'sync':<15}"
                         )
                 self.sync_model()
                 current_entities.dimtags = list(set(cut[0]))
@@ -332,17 +343,33 @@ def mesh(
                     else additive_entity.resolutions
                 )
                 updated_entities = []
+                # Parse additive physical names
                 for index, structural_entity in enumerate(structural_entity_list):
                     removeTool = False if index + 1 < len(structural_entities) else True
-                    # Create new composite physical_name
-                    additive_name = f"{structural_entity.physical_name}{addition_delimiter}{additive_entity.physical_name}"
+
+                    # Add all physical names
+                    additive_names = []
+                    if addition_addition_physicals:
+                        additive_names.extend(additive_entity.physical_name)
+                    if addition_structural_physicals:
+                        additive_names.extend(structural_entity.physical_name)
+                    if addition_intersection_physicals:
+                        additive_names.extend(
+                            [
+                                f"{x}{addition_delimiter}{y}"
+                                for x in structural_entity.physical_name
+                                for y in additive_entity.physical_name
+                            ]
+                        )
+                    additive_names = tuple(additive_names)
+
+                    # Add the intersection and complement as new LabeledEntities
                     intersection = self.occ.intersect(
                         structural_entity.dimtags,
                         additive_dimtags,
                         removeObject=False,  # Don't remove yet
                         removeTool=removeTool,  # Only remove at the end
                     )
-                    # Add the intersection and complement as new LabeledEntities
                     # No overlap case (no intersection)
                     if intersection[0] is []:
                         self.occ.synchronize()
@@ -380,7 +407,7 @@ def mesh(
                             LabeledEntities(
                                 index=structural_entity.index,
                                 dimtags=intersection[0],
-                                physical_name=additive_name,
+                                physical_name=additive_names,
                                 keep=structural_entity.keep,
                                 model=self.model,
                                 resolutions=structural_entity_resolutions

meshwell/polysurface.py

@@ -19,7 +19,7 @@ class PolySurface(BaseModel):
         ...
     )
     model: Any
-    physical_name: Optional[str] = Field(None)
+    physical_name: Optional[str | tuple[str, ...]] = Field(None)
     mesh_order: float | None = None
     mesh_bool: bool = Field(True)
     additive: bool = Field(False)
@@ -32,7 +32,7 @@ def __init__(
         self,
         polygons: Union[Polygon, List[Polygon], MultiPolygon, List[MultiPolygon]],
         model: Any,
-        physical_name: Optional[str] = None,
+        physical_name: Optional[str | tuple[str, ...]] = None,
         mesh_order: float | None = None,
         mesh_bool: bool = True,
         additive: bool = False,
@@ -57,7 +57,10 @@ def __init__(
         self.model = model
 
         self.mesh_order = mesh_order
-        self.physical_name = physical_name
+        if isinstance(physical_name, str):
+            self.physical_name = [physical_name]
+        else:
+            self.physical_name = physical_name
         self.mesh_bool = mesh_bool
         self.dimension = 2
         self.resolutions = resolutions

meshwell/prism.py

@@ -21,7 +21,7 @@ class Prism(BaseModel):
     )
     buffers: Dict[float, float] = Field(...)
     model: Any
-    physical_name: Optional[str] = Field(None)
+    physical_name: Optional[str | tuple[str, ...]] = Field(None)
     mesh_order: float | None = None
     mesh_bool: bool = Field(True)
     additive: bool = Field(False)
@@ -39,7 +39,7 @@ def __init__(
         polygons: Union[Polygon, List[Polygon], MultiPolygon, List[MultiPolygon]],
         buffers: Dict[float, float],
         model: Any,
-        physical_name: Optional[str] = None,
+        physical_name: Optional[str | tuple[str, ...]] = None,
         mesh_order: float | None = None,
         mesh_bool: bool = True,
         additive: bool = False,
@@ -78,7 +78,12 @@ def __init__(
 
         # Mesh order and name
         self.mesh_order = mesh_order
-        self.physical_name = physical_name
+
+        # Format physical name
+        if isinstance(physical_name, str):
+            self.physical_name = [physical_name]
+        else:
+            self.physical_name = physical_name
         self.mesh_bool = mesh_bool
         self.dimension = 3
         self.resolutions = resolutions

meshwell/tag.py

@@ -1,19 +1,36 @@
 from typing import List
 import gmsh
-from itertools import combinations
+from itertools import combinations, product
 
 
 def tag_entities(entity_list: List):
     """Adds physical labels to the entities in the model."""
+    # One pass to get the global name --> dimtags mapping
+    names_to_tags = {
+        0: {},
+        1: {},
+        2: {},
+        3: {},
+    }
     for entities in entity_list:
-        if entities.physical_name:
-            gmsh.model.addPhysicalGroup(
-                entities.get_dim(), entities.get_tags(), name=entities.physical_name
-            )
+        dim = entities.get_dim()
+        for physical_name in entities.physical_name:
+            if physical_name not in names_to_tags[dim]:
+                names_to_tags[dim][physical_name] = []
+            names_to_tags[dim][physical_name].extend(entities.get_tags())
+
+    for dim in names_to_tags.keys():
+        for physical_name, tags in names_to_tags[dim].items():
+            gmsh.model.addPhysicalGroup(dim, tags, name=physical_name)
 
 
 def tag_interfaces(entity_list: List, max_dim: int, boundary_delimiter: str):
     """Adds physical labels to the interfaces between entities in entity_list."""
+    names_to_tags = {
+        0: {},
+        1: {},
+        2: {},
+    }
     for entity1, entity2 in combinations(entity_list, 2):
         if entity1.physical_name == entity2.physical_name:
             continue
@@ -22,18 +39,26 @@ def tag_interfaces(entity_list: List, max_dim: int, boundary_delimiter: str):
         elif entity1.get_dim() != max_dim:
             continue
         else:
+            dim = entity1.get_dim() - 1
             common_interfaces = list(
                 set(entity1.boundaries).intersection(entity2.boundaries)
             )
             if common_interfaces:
-                # Remember which boundaries were interfaces with another entity
+                # Update entity interface logs
                 entity1.interfaces.extend(common_interfaces)
                 entity2.interfaces.extend(common_interfaces)
-                gmsh.model.addPhysicalGroup(
-                    max_dim - 1,
-                    common_interfaces,
-                    name=f"{entity1.physical_name}{boundary_delimiter}{entity2.physical_name}",
-                )
+                # Prepare physical tags
+                for entity1_physical_name, entity2_physical_name in product(
+                    entity1.physical_name, entity2.physical_name
+                ):
+                    interface_name = f"{entity1_physical_name}{boundary_delimiter}{entity2_physical_name}"
+                    if interface_name not in names_to_tags[dim]:
+                        names_to_tags[dim][interface_name] = []
+                    names_to_tags[dim][interface_name].extend(common_interfaces)
+
+    for dim in names_to_tags.keys():
+        for physical_name, tags in names_to_tags[dim].items():
+            gmsh.model.addPhysicalGroup(dim, tags, name=physical_name)
 
     return entity_list
 
@@ -42,12 +67,24 @@ def tag_boundaries(
     entity_list: List, max_dim: int, boundary_delimiter: str, mesh_edge_name: str
 ):
     """Adds physical labels to the boundaries of the entities in entity_list."""
+    names_to_tags = {
+        0: {},
+        1: {},
+        2: {},
+    }
     for entity in entity_list:
         if entity.get_dim() != max_dim:
             continue
+        dim = entity.get_dim() - 1
         boundaries = list(set(entity.boundaries) - set(entity.interfaces))
-        gmsh.model.addPhysicalGroup(
-            max_dim - 1,
-            boundaries,
-            name=f"{entity.physical_name}{boundary_delimiter}{mesh_edge_name}",
-        )
+        for entity_physical_name in entity.physical_name:
+            boundary_name = (
+                f"{entity_physical_name}{boundary_delimiter}{mesh_edge_name}"
+            )
+            if boundary_name not in names_to_tags[dim]:
+                names_to_tags[dim][boundary_name] = []
+            names_to_tags[dim][boundary_name].extend(boundaries)
+
+    for dim in names_to_tags.keys():
+        for physical_name, tags in names_to_tags[dim].items():
+            gmsh.model.addPhysicalGroup(dim, tags, name=physical_name)