Merge pull request #80 from simbilod/fix_resolutionspec

basic framework for class-based refinement

meshwell/model.py

@@ -482,7 +482,6 @@ def mesh(
                     refinement_field_indices,
                     refinement_max_index,
                     default_characteristic_length,
-                    final_entity_list,
                 )
 
             # Use the smallest element size overall

meshwell/resolution.py

@@ -1,103 +1,192 @@
 import numpy as np
 import copy
+from typing import Literal, Any
 from pydantic import BaseModel
 
 
 class ResolutionSpec(BaseModel):
+    """A ResolutionSpec is attached to a pre-CAD entity.
+
+    It sets a mesh size field (see child classes) to the resulting post-CAD volumes, surfaces, curves, or points.
+
+    The volumes, surfaces, curves can be filtered based on their mass (volume, area, length). Points can be filtered based on the length of the curve they belong to.
     """
-    Object holding resolution information for an entity and its boundaries.
 
-    # FIXME: make a better ResolutionSpec class that handles the entity/boundary distinction
+    # Eventually we can add flags here to also consider proximity to other specific physicals (e.g. shared interfaces)
+    apply_to: Literal["volumes", "surfaces", "curves", "points"]
+    min_mass: float = 0
+    max_mass: float = np.inf
 
-    Arguments:
-        # Volume
-        volume_resolution (float): resolution of the volume (3D). No effect if the entity is 2D.
-            defaults to inf --> default resolution, since the local resolution is always the minimum of all size fields
+    @property
+    def entity_str(self):
+        """Convenience wrapper."""
+        if self.apply_to == "volumes":
+            return "RegionsList"
+        elif self.apply_to == "surfaces":
+            return "SurfacesList"
+        elif self.apply_to == "curves":
+            return "CurvesList"
+        elif self.apply_to == "points":
+            return "PointsList"
 
-        # Surface
-        surface_resolution (float): resolution of the surface (2D) or of all the surfaces touching the volume (3D)
-            defaults to inf --> default resolution (2D) or volume_resolution (3D), since the local resolution is always the minimum of all size fields
+    @property
+    def target_dimension(self):
+        """Convenience wrapper."""
+        if self.apply_to == "volumes":
+            return 3
+        elif self.apply_to == "surfaces":
+            return 2
+        elif self.apply_to == "curves":
+            return 1
+        elif self.apply_to == "points":
+            return 0
 
-        # Curves
-        curve_resolution (float): resolution of curves constituting the volumes' surfaces (3D) or surfaces (2D)
-            defaults to inf --> surface_resolution, since the local resolution is always the minimum of all size fields
 
-        # Points
-        point_resolution (float): resolution of points constituting the volumes' surfaces' curves (3D) or surfaces' curves (2D)
-            defaults to inf --> curve_resolution, since the local resolution is always the minimum of all size fields
-            can be filtered by the length of the associated curves
-    """
+class ConstantInField(ResolutionSpec):
+    """Constant resolution within the entities."""
+
+    resolution: float
+
+    def apply(
+        self,
+        model: Any,
+        current_field_index: int,
+        entities_mass_dict,
+        refinement_field_indices,
+    ) -> int:
+        new_field_indices = []
 
-    # Volume
-    resolution_volumes: float | None = None
-    min_volumes: float = 0
-    max_volumes: float = np.inf
-    # Surface
-    resolution_surfaces: float | None = None
-    min_area_surfaces: float = 0
-    max_area_surfaces: float = np.inf
-    distmax_surfaces: float | None = None
-    sizemax_surfaces: float | None = None
-    surface_sigmoid: bool = False
-    surface_per_sampling_surfaces: float | None = None
-    sampling_surface_max: int = 100
-    # Curve
-    resolution_curves: float | None = None
-    min_length_curves: float = 0
-    max_length_curves: float = np.inf
-    distmax_curves: float | None = None
-    sizemax_curves: float | None = None
-    curve_sigmoid: bool = False
-    length_per_sampling_curves: float | None = None
-    sampling_curve_max: int = 100
-    # Point
-    resolution_points: float | None = None
-    min_length_curves_for_points: float = 0
-    max_length_curves_for_points: float = np.inf
-    distmax_points: float | None = None
-    sizemax_points: float | None = None
-    point_sigmoid: bool = False
+        model.mesh.field.add("MathEval", current_field_index)
+        model.mesh.field.setString(current_field_index, "F", f"{self.resolution}")
+        model.mesh.field.add("Restrict", current_field_index + 1)
+        model.mesh.field.setNumber(
+            current_field_index + 1, "InField", current_field_index
+        )
+        model.mesh.field.setNumbers(
+            current_field_index + 1,
+            self.entity_str,
+            list(entities_mass_dict.keys()),
+        )
+        new_field_indices = (current_field_index + 1,)
+        current_field_index += 2
+
+        return new_field_indices, current_field_index
 
     def refine(self, resolution_factor: float):
         result = copy.copy(self)
-        if result.resolution_volumes is not None:
-            result.resolution_volumes *= resolution_factor
-        if result.resolution_surfaces is not None:
-            result.resolution_surfaces *= resolution_factor
-        if result.sizemax_surfaces is not None:
-            result.sizemax_surfaces *= resolution_factor
-        if result.resolution_curves is not None:
-            result.resolution_curves *= resolution_factor
-        if result.sizemax_curves is not None:
-            result.sizemax_curves *= resolution_factor
-        if result.resolution_points is not None:
-            result.resolution_points *= resolution_factor
-        if result.sizemax_points is not None:
-            result.sizemax_points *= resolution_factor
+        if result.resolution is not None:
+            result.resolution *= resolution_factor
+
         return result
 
-    def calculate_sampling(self, mass_per_sampling, mass, max_sampling):
-        if mass_per_sampling is None:
-            return 2  # avoid int(inf) error
-        else:
-            return min(max(2, int(mass / mass_per_sampling)), max_sampling)
-
-    def calculate_sampling_surface(self, area):
-        if self.surface_per_sampling_surfaces:
-            return self.calculate_sampling(
-                self.surface_per_sampling_surfaces, area, self.sampling_surface_max
-            )
-        else:
-            return self.calculate_sampling(
-                0.5 * self.resolution_surfaces, area, self.sampling_surface_max
-            )
-
-    def calculate_sampling_curve(self, length):
-        if self.length_per_sampling_curves:
-            return self.calculate_sampling(
-                self.length_per_sampling_curves, length, self.sampling_curve_max
-            )
-        else:
-            return self.calculate_sampling(
-                0.5 * self.resolution_curves, length, self.sampling_curve_max
-            )
+
+class SampledField(ResolutionSpec):
+    """Shared functionality for size fields that require sampling the entities at points."""
+
+    mass_per_sampling: float | None = None
+    max_sampling: int = 100
+    sizemin: float
+
+    def calculate_samplings(self, entities_mass_dict):
+        """Calculates a more optimal sampling from the masses"""
+
+        if self.mass_per_sampling is None:
+            # Default sampling is half the minimum resolution
+            mass_per_sampling = 0.5 * self.sizemin
+
+        return {
+            tag: min(max(2, int(mass / mass_per_sampling)), self.max_sampling)
+            for tag, mass in entities_mass_dict.items()
+        }
+
+    def refine(self, resolution_factor: float):
+        result = copy.copy(self)
+        if result.sizemax is not None:
+            result.sizemax *= resolution_factor
+        if result.sizemin is not None:
+            result.sizemin *= resolution_factor
+
+
+class ThresholdField(SampledField):
+    """Linear growth of the resolution away from the entity"""
+
+    sizemax: float
+    sizemin: float
+    distmin: float = 0
+    distmax: float
+
+    def apply(
+        self,
+        model: Any,
+        current_field_index: int,
+        entities_mass_dict,
+        refinement_field_indices,
+    ) -> int:
+        new_field_indices = []
+
+        # Compute samplings
+        samplings_dict = self.calculate_samplings(entities_mass_dict)
+
+        # FIXME: It is computationally cheaper to have a large sampling on all the curves rather than one field per curve; but there is probably an optimum somewhere.
+        # FOr instance, the distribution should be very skewed (tiny vertical curves, tiny curves in bends, vs long horizontal ones), so there may be benefits for a small number of optimized fields.
+        samplings = max(samplings_dict.values())
+        entities = list(entities_mass_dict.keys())
+
+        model.mesh.field.add("Distance", current_field_index)
+        model.mesh.field.setNumbers(current_field_index, self.entity_str, entities)
+        model.mesh.field.setNumber(current_field_index, "Sampling", samplings)
+        model.mesh.field.add("Threshold", current_field_index + 1)
+        model.mesh.field.setNumber(
+            current_field_index + 1, "InField", current_field_index
+        )
+        model.mesh.field.setNumber(current_field_index + 1, "SizeMin", self.sizemin)
+        model.mesh.field.setNumber(current_field_index + 1, "DistMin", self.distmin)
+        if self.sizemax and self.distmax:
+            model.mesh.field.setNumber(current_field_index + 1, "SizeMax", self.sizemax)
+            model.mesh.field.setNumber(current_field_index + 1, "DistMax", self.distmax)
+        model.mesh.field.setNumber(current_field_index + 1, "StopAtDistMax", 1)
+        new_field_indices = (current_field_index + 1,)
+        current_field_index += 2
+
+        return new_field_indices, current_field_index
+
+
+class ExponentialField(SampledField):
+    """Exponential growth of the characteristic length away from the entity"""
+
+    growth_factor: float
+
+    def apply(
+        self,
+        model: Any,
+        current_field_index: int,
+        entities_mass_dict,
+        refinement_field_indices,
+    ) -> int:
+        new_field_indices = []
+
+        # Compute samplings
+        samplings_dict = self.calculate_samplings(entities_mass_dict)
+
+        # FIXME: It is computationally cheaper to have a large sampling on all the curves rather than one field per curve; but there is probably an optimum somewhere.
+        # FOr instance, the distribution should be very skewed (tiny vertical curves, tiny curves in bends, vs long horizontal ones), so there may be benefits for a small number of optimized fields.
+        samplings = max(samplings_dict.values())
+        entities = list(entities_mass_dict.keys())
+
+        # Sampled distance field
+        model.mesh.field.add("Distance", current_field_index)
+        model.mesh.field.setNumbers(current_field_index, self.entity_str, entities)
+        model.mesh.field.setNumber(current_field_index, "Sampling", samplings)
+
+        # Math field
+        model.mesh.field.add("MathEval", current_field_index + 1)
+        model.mesh.field.setString(
+            current_field_index + 1,
+            "F",
+            f"({self.growth_factor}^F{current_field_index} - 1) + {self.sizemin}",
+        )
+
+        new_field_indices = (current_field_index + 1,)
+        current_field_index += 2
+
+        return new_field_indices, current_field_index

meshwell/tag.py

@@ -13,11 +13,11 @@ def tag_entities(entity_list: List):
         3: {},
     }
     for entities in entity_list:
-        dim = entities.get_dim()
+        dim = entities.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())
+            names_to_tags[dim][physical_name].extend(entities.tags)
 
     for dim in names_to_tags.keys():
         for physical_name, tags in names_to_tags[dim].items():
@@ -34,12 +34,12 @@ def tag_interfaces(entity_list: List, max_dim: int, boundary_delimiter: str):
     for entity1, entity2 in combinations(entity_list, 2):
         if entity1.physical_name == entity2.physical_name:
             continue
-        elif entity1.get_dim() != entity2.get_dim():
+        elif entity1.dim != entity2.dim:
             continue
-        elif entity1.get_dim() != max_dim:
+        elif entity1.dim != max_dim:
             continue
         else:
-            dim = entity1.get_dim() - 1
+            dim = entity1.dim - 1
             common_interfaces = list(
                 set(entity1.boundaries).intersection(entity2.boundaries)
             )
@@ -73,9 +73,9 @@ def tag_boundaries(
         2: {},
     }
     for entity in entity_list:
-        if entity.get_dim() != max_dim:
+        if entity.dim != max_dim:
             continue
-        dim = entity.get_dim() - 1
+        dim = entity.dim - 1
         boundaries = list(set(entity.boundaries) - set(entity.interfaces))
         for entity_physical_name in entity.physical_name:
             boundary_name = (