Merge pull request #74 from simbilod/sampling

Sampling

docs/changelog.md

@@ -1,5 +1,9 @@
 # [Changelog](https://keepachangelog.com/en/1.0.0/)
 
+## 1.2.3
+
+- expose Sampling parameter [PR#74](https://github.com/simbilod/meshwell/pull/74)
+
 ## 1.2.2
 
 - sigmoid resolution tapering [PR#73](https://github.com/simbilod/meshwell/pull/73)

meshwell/labeledentity.py

@@ -55,6 +55,8 @@ def add_refinement_fields_to_model(
     ):
         """
         Adds refinement fields to the model based on base_resolution and resolution info.
+
+        # FIXME: make a better ResolutionSpec class that handles the entity/boundary distinction
         """
         n = refinement_max_index
 
@@ -98,6 +100,12 @@ def add_refinement_fields_to_model(
                     boundary_sizemax = resolutionspec.sizemax_surfaces
                     boundary_distmax = resolutionspec.distmax_surfaces
                     boundary_sigmoid = resolutionspec.surface_sigmoid
+                    boundaries_samplings = {
+                        boundary: resolutionspec.calculate_sampling_surface(
+                            self.model.occ.getMass(2, boundary)
+                        )
+                        for boundary in boundaries
+                    }
 
                     # Check condition on surface curves
                     boundary_lines = [
@@ -116,6 +124,12 @@ def add_refinement_fields_to_model(
                     boundary_lines_sizemax = resolutionspec.sizemax_curves
                     boundary_lines_distmax = resolutionspec.distmax_curves
                     boundary_line_sigmoid = resolutionspec.curve_sigmoid
+                    boundary_lines_samplings = {
+                        boundary_line: resolutionspec.calculate_sampling_curve(
+                            self.model.occ.getMass(1, boundary_line)
+                        )
+                        for boundary_line in boundary_lines
+                    }
 
                 elif self.get_dim() == 2:
                     entity_str = "SurfacesList"
@@ -150,6 +164,12 @@ def add_refinement_fields_to_model(
                     boundary_sizemax = resolutionspec.sizemax_curves
                     boundary_distmax = resolutionspec.distmax_curves
                     boundary_sigmoid = resolutionspec.curve_sigmoid
+                    boundaries_samplings = {
+                        boundary: resolutionspec.calculate_sampling_curve(
+                            self.model.occ.getMass(1, boundary)
+                        )
+                        for boundary in boundaries
+                    }
 
                 elif self.get_dim() == 1:
                     entity_str = "CurvesList"
@@ -171,10 +191,12 @@ def add_refinement_fields_to_model(
                     refinement_field_indices.extend((n + 1,))
                     n += 2
 
-                if boundaries:
+                for boundary in boundaries:
                     self.model.mesh.field.add("Distance", n)
-                    self.model.mesh.field.setNumbers(n, boundary_str, boundaries)
-                    self.model.mesh.field.setNumber(n, "Sampling", 100)
+                    self.model.mesh.field.setNumbers(n, boundary_str, [boundary])
+                    self.model.mesh.field.setNumber(
+                        n, "Sampling", boundaries_samplings[boundary]
+                    )
                     self.model.mesh.field.add("Threshold", n + 1)
                     self.model.mesh.field.setNumber(n + 1, "InField", n)
                     self.model.mesh.field.setNumber(
@@ -195,10 +217,12 @@ def add_refinement_fields_to_model(
                     refinement_field_indices.extend((n + 1,))
                     n += 2
 
-                if boundary_lines:
+                for boundary_line in boundary_lines:
                     self.model.mesh.field.add("Distance", n)
-                    self.model.mesh.field.setNumbers(n, curves_str, boundary_lines)
-                    self.model.mesh.field.setNumber(n, "Sampling", 100)
+                    self.model.mesh.field.setNumbers(n, curves_str, [boundary_line])
+                    self.model.mesh.field.setNumber(
+                        n, "Sampling", boundary_lines_samplings[boundary_line]
+                    )
                     self.model.mesh.field.add("Threshold", n + 1)
                     self.model.mesh.field.setNumber(n + 1, "InField", n)
                     self.model.mesh.field.setNumber(

meshwell/resolution.py

@@ -7,6 +7,8 @@ class ResolutionSpec(BaseModel):
     """
     Object holding resolution information for an entity and its boundaries.
 
+    # FIXME: make a better ResolutionSpec class that handles the entity/boundary distinction
+
     Arguments:
         # Volume
         volume_resolution (float): resolution of the volume (3D). No effect if the entity is 2D.
@@ -36,13 +38,17 @@ class ResolutionSpec(BaseModel):
     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 = np.inf
     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 = np.inf
     distmax_points: float | None = None
@@ -62,3 +68,29 @@ def refine(self, resolution_factor: float):
         if result.sizemax_points is not None:
             result.sizemax_points *= resolution_factor
         return result
+
+    def calculate_sampling(self, mass_per_sampling, mass, max_sampling):
+        if mass_per_sampling == np.inf:
+            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
+            )

pyproject.toml

@@ -14,7 +14,7 @@ classifiers = [
 	"Programming Language :: Python :: 3.11",
 	"Operating System :: OS Independent",
 ]
-version="1.2.2"
+version="1.2.3"
 authors = [
     {name = "Simon Bilodeau", email = "[email protected]"},
 ]