diff --git a/.run/pytest test_field.run.xml b/.run/pytest test_field.run.xml
index 059c23d..f3e4989 100644
--- a/.run/pytest test_field.run.xml	
+++ b/.run/pytest test_field.run.xml	
@@ -1,13 +1,16 @@
 <component name="ProjectRunConfigurationManager">
   <configuration default="false" name="pytest test_field" type="tests" factoryName="py.test">
     <module name="bgem" />
+    <option name="ENV_FILES" value="" />
     <option name="INTERPRETER_OPTIONS" value="" />
     <option name="PARENT_ENVS" value="true" />
     <option name="SDK_HOME" value="/usr/bin/python3.8" />
     <option name="WORKING_DIRECTORY" value="" />
-    <option name="IS_MODULE_SDK" value="false" />
+    <option name="IS_MODULE_SDK" value="true" />
     <option name="ADD_CONTENT_ROOTS" value="true" />
     <option name="ADD_SOURCE_ROOTS" value="true" />
+    <option name="DEBUG_JUST_MY_CODE" value="true" />
+    <option name="RUN_TOOL" value="" />
     <option name="_new_keywords" value="&quot;&quot;" />
     <option name="_new_parameters" value="&quot;-s&quot;" />
     <option name="_new_additionalArguments" value="&quot;&quot;" />
diff --git a/.run/pytest test_gmsh.run.xml b/.run/pytest test_gmsh.run.xml
index d726504..d7dc2d1 100644
--- a/.run/pytest test_gmsh.run.xml	
+++ b/.run/pytest test_gmsh.run.xml	
@@ -1,13 +1,16 @@
 <component name="ProjectRunConfigurationManager">
   <configuration default="false" name="pytest test_gmsh" type="tests" factoryName="py.test">
     <module name="bgem" />
+    <option name="ENV_FILES" value="" />
     <option name="INTERPRETER_OPTIONS" value="" />
     <option name="PARENT_ENVS" value="true" />
     <option name="SDK_HOME" value="/usr/bin/python3.8" />
     <option name="WORKING_DIRECTORY" value="" />
-    <option name="IS_MODULE_SDK" value="false" />
+    <option name="IS_MODULE_SDK" value="true" />
     <option name="ADD_CONTENT_ROOTS" value="true" />
     <option name="ADD_SOURCE_ROOTS" value="true" />
+    <option name="DEBUG_JUST_MY_CODE" value="true" />
+    <option name="RUN_TOOL" value="" />
     <option name="_new_keywords" value="&quot;&quot;" />
     <option name="_new_parameters" value="&quot;&quot;" />
     <option name="_new_additionalArguments" value="&quot;&quot;" />
diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 0000000..b6f133b
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,96 @@
+# CODEX context for BGEM library
+
+## Project summary
+
+Bgem library is a intended as an extension or wrapper of GMSH python API.
+It helps to keep tracks of the geometry objects during geometrical operations;
+it introduces `ObjectSet` class keeping the tracks of dimtags and provides wrapper
+functions to GMSH API calls.
+
+Further it provides better interface to GMSH fields, reader/writer for GMSH meshes
+and a HealMesh class that can help fixing degenerated elements in a given mesh.
+
+Submodule called `Stochastic` provides generator for DFN (Discrete Fracture Network)
+based on distributions used in geology.
+
+## Project summary (short)
+
+Robust open source tool for creation of parametric geometries and computational meshes via. Python code. 
+Primary focus are hydrogeological applications with geometries including both random fractures and deterministic natural or anthropogenic features.
+Bgem library is a intended as an extension and/or wrapper of GMSH python API.
+
+## CODEX Ignore Folders
+- Not specified
+
+## CODEX Readonly Folders
+- Not specified
+
+## Project Structure
+- `src/bgem` - source files
+- `tests` - pytest based unit tests of individual source blocks
+- `tutorials` - integrated tests, currently insufficient
+
+## `src/bgem/` Project Source Structure
+
+- `bspline/` - 
+- `geometry/` - 
+- `gmsh/` - main functions providing interface for GMSH API
+- `polygons/` - 
+- `stochastic/` - code for DFN generation
+
+
+## CODEX Guidelines
+- treat keyword 'AGENT:' in comments as a source context dependent message for your further development
+- any comment containing `AGENT:` is an active developer instruction
+- NEVER remove, rewrite, or move an `AGENT:` comment unless you implement that instruction in the same change
+- if you make only a local fix around an `AGENT:` comment, leave the comment untouched
+- if an `AGENT:` instruction looks outdated or wrong, ask before removing it
+- never stage or commit changes yourself
+- Always review your changes before finishing for human review.
+- Do not touch anything outside the repository directory, unless directly asked!
+
+## Status tracking
+- `STATUS.md` is the handoff log for interrupted or multi-turn work. Update it when a session ends with unfinished relevant work, when the user asks for a status review, or when a fresh checkpoint would help the next session continue without re-discovery.
+- Keep newest entry first. Do not rewrite older entries except to fix clearly factual mistakes.
+- Start each entry with one line in this form:
+  `` `YYYY-MM-DD`: `<commit>` @ `<branch>` by `<author>` ``
+- Under each entry keep exactly these sections in order:
+  `## Goal`
+  short statement of the intended task or checkpoint scope
+  `## Changes summary`
+  flat bullets describing committed changes first, then important staged/unstaged/untracked changes if they are relevant to continuing the task
+  `## Verified`
+  flat bullets with commands actually run and the important observed result
+  `## Open items`
+  flat bullets for remaining risks, missing verification, known breakage, or next recommended step
+- Record only repo-relevant facts that help continuation. Skip conversational history, speculation, and incidental noise.
+- When the worktree is dirty, explicitly distinguish committed, staged, unstaged, and untracked changes, but mention only files relevant to the tracked task.
+- Prefer clickable file links for important files mentioned in `STATUS.md`.
+- If verification was partial, say so plainly. Do not imply a full test pass when only compile checks, a single testcase, or CLI smoke checks were run.
+- If a run failed, record the failing command and the actionable failure mode instead of hiding it.
+- Before finishing a task that changed the practical project state, review whether `STATUS.md` still matches the actual branch/worktree state and update it if needed.
+
+## Coding rules
+- Best code, is no code!
+- prefer functional style with poor functions;
+  ideally do not change objects after construction, all methods do calculations
+  only reading the data in the class
+- prefer high level code: numpy, pandas, xarray instead loops and native python structures (lists, dicts)
+- use logging
+- Use logging for debug outputs.
+- use pathlib
+- use attrs for dataclasses
+- use attrs staticmethod/classmethod technique to construct from other data then is stored in the dataclass
+- Be defensive, with strong checks, but only for the user input data.
+  That means error inputs must raise early. Therefore only check for existing keys in input dicts
+  if these will be required down in a long calculation. Otherwise just let KeyError do the job.
+- do not use "guess" default values, only obvious defaults
+- Do not use other config keys in the case of a KeyError, just throw early.
+- Do just basic asserts for consistency for function inputs.
+- Can add more asserts if needed during debugging.
+- NEVER resolve test errors by try blocks
+- NEVER add runtime fallbacks or import shims to compensate for a broken or incomplete environment.
+  If a declared dependency or tool is missing, report the environment problem plainly and fix the environment or tests around it, but do not implement code workarounds.
+- NEVER write "self explanatory" into comments
+- in comments indicate by ?? if you are not certain about intent of particular variable, function, parameter ...
+
diff --git a/requirements.txt b/requirements.txt
index 52590b4..f16d61c 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -3,3 +3,4 @@ numpy
 scipy
 bih
 plotly
+attrs
diff --git a/setup.py b/setup.py
index 99da665..436838f 100644
--- a/setup.py
+++ b/setup.py
@@ -56,7 +56,7 @@
     include_package_data=True,
     zip_safe=False,
     #install_requires=['numpy', 'scipy', 'bih', 'gmsh-sdk<=4.5.1'],
-    install_requires=['numpy>=1.13.4', 'pandas', 'scipy', 'bih', 'gmsh>=4.10.4'],
+    install_requires=['numpy>=1.13.4', 'pandas', 'scipy', 'bih', 'gmsh>=4.10.4', 'attrs'],
     # incompatible changes in SDK in release 4.6.0 to be changed in the new release of bgem
     python_requires='>=3',
     # extras_require={
diff --git a/src/bgem/gmsh/field.py b/src/bgem/gmsh/field.py
index 9190847..95b7218 100644
--- a/src/bgem/gmsh/field.py
+++ b/src/bgem/gmsh/field.py
@@ -347,6 +347,21 @@ def construct(self, model):
         # creata MathEval
 
 
+# def constant(value):
+#     """
+#     Define constant field.
+#     Unfortunately much slower then the Box field alternative.
+#     Seems that MathExpr evaluation time is proportional to the string length, i.e. it is parsed every time.
+#     :param value:
+#     :return:
+#     """
+#     return FieldExpr(str(value), [])
+
+def constant(value):
+    """
+    Make a field with constant value = value.
+    """
+    return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value)
 
 """
 Predefined coordinate fields.
@@ -414,11 +429,6 @@ def box(pt_min: Point, pt_max: Point, v_in:float, v_out:float=1e300) -> Field:
 
 
 
-def constant(value):
-    """
-    Make a field with constant value = value.
-    """
-    return box((-1, -1, -1), (1, 1, 1), v_in=value, v_out=value)
 
 
 
@@ -549,30 +559,41 @@ def minimum(*fields: Field) -> Field:
 # Structured
 
 
-def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
+
+
+def linear(x_field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
     """
-    Geometric interpolation between points a = (x0, y0) and b = (x1, y1).
+    Linear transformation of the scalar field given by two points.
+    Prescribed mapping a[0] to a[1], and b[0] to b[1].
     """
-    x = field
     x0, y0 = a
     x1, y1 = b
-    l0, l1 = math.log(y0), math.log(y1)
-    a = (l1 - l0) / (x1 - x0)
-    x_avg = (l1 * x0 - l0 * x1) / (l1 - l0)
-    return current_module.exp(a * (x - x_avg))
+    a = (y1 - y0) / (x1 - x0)
+    #x_avg = (y1 * x0 - y0 * x1) / (y1 - y0)
+    x_avg = (x0 + x1) / 2.0
+    y_avg = (y0 + y1) / 2.0
+
+    if a > 1e10:
+        # Following formula could provide more stable results
+        # but needs one more scalar addition
+        # that makes the test 3D meshing about 2% slower.
+        return a * (x_field - x_avg) + y_avg
+    else:
+        # Just two field operations.
+        add = - (a * x0 - y0)
+        return a * x_field + add
 
 
-def linear(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
+def geometric(field:Field, a:Tuple[float, float], b:Tuple[float, float]) -> Field:
     """
-    Linear transformation of the scalar field given by two points.
-    Prescribed mapping a[0] to a[1], and b[0] to b[1].
+    Geometric interpolation between points a = (x0, y0) and b = (x1, y1).
+
     """
     x = field
     x0, y0 = a
     x1, y1 = b
-    a = (y1 - y0) / (x1 - x0)
-    x_avg = (y1 * x0 - y0 * x1) / (y1 - y0)
-    return a * (x - x_avg)
+    l0, l1 = math.log(y0), math.log(y1)
+    return current_module.exp(linear(field, (x0, l0), (x1, l1)))
 
 def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: float = 1) -> Field:
     """
@@ -583,9 +604,7 @@ def polynomial(field:Field, a:Tuple[float, float], b:Tuple[float, float], q: flo
     x0, y0 = a
     x1, y1 = b
     l0, l1 = y0 ** (1/q), y1 ** (1/q)
-    a = (l1 - l0) / (x1 - x0)
-    x_avg = (l1 * x0 - l0 * x1) / (l1 - l0)
-    return (a * (x - x_avg)) ** q
+    return linear(field, (x0, l0), (x1, l1)) ** q
 
 def threshold(field:Field, lower_bound:Tuple[float, float],
               upper_bound:Tuple[float, float]=None, sigmoid:bool=False) -> Field:
diff --git a/src/bgem/gmsh/gmsh.py b/src/bgem/gmsh/gmsh.py
index 412891c..83de3d6 100644
--- a/src/bgem/gmsh/gmsh.py
+++ b/src/bgem/gmsh/gmsh.py
@@ -7,6 +7,8 @@
 import gmsh
 import re
 import warnings
+import inspect
+import pathlib
 
 from bgem.gmsh import gmsh_exceptions
 from bgem.gmsh import options as gmsh_options
@@ -485,6 +487,9 @@ def import_shapes(self, fileName, highestDimOnly=True):
         :param highestDimOnly:
 
         """
+        file_path = pathlib.Path(fileName)
+        if not file_path.exists():
+            raise FileNotFoundError(f"Error: The file '{file_path}' does not exist!")
         shapes = self.model.importShapes(fileName, highestDimOnly=highestDimOnly)
         self._need_synchronize = True
         return ObjectSet(self, shapes, [Region.default_region[dim] for dim, _ in shapes])
@@ -659,6 +664,80 @@ def make_mesh(self, objects: List['ObjectSet'], dim=3, eliminate=True) -> None:
         if bad_entities:
             print("Bad entities:", bad_entities)
 
+    def get_elements(self, obj_set: 'ObjectSet'):
+        """
+        Get a (dim, tag) dict -> `elementTypes', `elementTags', `nodeTags'
+        :param obj_set:
+        :return:
+        """
+        return {
+                (dim, tag):  gmsh.model.mesh.getElements(dim, tag)
+                for dim, tag in obj_set.dimtags
+            }
+
+    def mesh_new_region_from_elements(self, dim, region_name, el_tags, el_coords):
+        el_type = [15, 1, 2, 4][dim]
+        i_entity = self.model.addDiscreteEntity(dim, -1, [])
+        self.model.mesh.addElementsByType(i_entity,  el_type, el_tags, el_coords)
+
+    # def mesh_level_field_regions(self, obj_set:'ObjectSet', level_fn=None, bins:List[float] = None, region: Union[str, List[str]] = None):
+    #     """
+    #     Assign region to the elements of the given model entities (dimtags).
+    #     level_fn is given, the regions in form "<region>_<i>" are assigned to the elements
+    #     on which the level_fn is in the bin [x_i, x_{i+1}]. The level_fn is evaluated in nodes,
+    #     elements on the boundary are assigned to the bin for which the nodes have larger total weight.
+    #     Weight of node X is distance of level_fn(X) from the bin borders.
+    #     :param model_entities:
+    #     :param level_fn: a scalar level-set function on set of 3D points, Nx3 -> N
+    #     :param bins: [x_0, x_1, .... ,x_n] define n bins
+    #     default single bin [0,1]
+    #     :param region: str
+    #     default: region = "levelset"
+    #     for List of regions, the ragions assigned to individual bins
+    #     :return:
+    #
+    #     TODO: It is more tricky to implement as GMSH could assign a single PhysicalGroup to an entity.
+    #     See tutorial X2 from GMSH manual.
+    #     So one has to create a new discrete entity for each discrete level bin.
+    #     gmsh.model.addDiscreteEntity(dim, tag)
+    #     Boundary entities may be specified to allo GMSH perform correct topology operations.
+    #     No need in our case as we will just write down the mesh.
+    #     Then we could possibly assign elements:
+    #     gmsh.model.mesh.addElementsByType(tag, elementType, elementTags, nodeTags):
+    #     """
+    #     entity_dict = self.get_elements(obj_set)
+    #     node_dict = {}
+    #     for dim,tag in entity_dict.keys():
+    #         node_ids, node_coords = self.model.getNodes(dim, tag)
+    #         node_dict.update(dict(zip(node_ids, node_coords)))
+    #     unique_ids, unique_coords = zip((k, v) for k, v in node_dict.items())
+    #     unique_levels = level_fn(unique_coords)
+    #     node_levels = dict(unique_ids, unique_levels)
+    #
+    #     self.model.model.addPhysicalGroup
+    #
+    #
+    #     def mean_level(nodes):
+    #         return np.mean(node_levels[node] for node in nodes)
+    #
+    #     def bin(level):
+    #         return 1
+    #     def reg_el_tags(el_tags, el_nodes):
+    #         el_regions = [bin(mean_level(nodes)) for nodes in el_nodes]
+    #         zip(el_tags, el_regions)
+    #
+    #
+    #
+    #     el_data = zip(el_tuple for el_tuple in entity_dict.values())
+    #     concatenate = lambda x : list(itertools.chain.from_iterable(x))
+    #     el_ids, el_tags, el_nodes = [concatenate(x) for x in el_data]
+    #     self.mesh_new_region_from_elements(dim, region, el_tags, el_nodes)
+    #
+    #
+    #     unique_bins = bin(unique_levels)
+    #     {n_id, for n_id, coord in node_dict.items()
+
+
     def write_brep(self, filename=None):
         self.synchronize()
         if filename is None:
@@ -743,7 +822,6 @@ def group(self, *obj_list: Union['ObjectSet', List['ObjectSet']]) -> 'ObjectSet'
 
 
 
-
 class ObjectSet:
     default_mesh_step = 0
 
@@ -793,6 +871,18 @@ def group(*obj_list: 'ObjectSet') -> 'ObjectSet':
         g.mesh_step_size = mesh_step_size
         return g
 
+    @property
+    def attr(self):
+        """
+        A dynamic dataclass of attribute maps.
+        - single attribute : defined name and aggregation operation, own global dimtag map
+        - global map assigning attribute dicts to dimtags
+          or track dimtags through operations with associated attributes.
+          However region and mesh_Step tracking didn't quite work.
+        -
+        :return:
+        """
+        pass
 
     def __init__(self, factory: 'GeometryOCC', dim_tags: List[DimTag], regions: List[Region]) -> None:
         self.factory = factory
@@ -818,6 +908,10 @@ def tags(self):
     def size(self):
         return len(self.dim_tags)
 
+    def min_dim(self):
+        dims = [d for d, t in self.dim_tags]
+        return min(dims, default=1)
+
     def max_dim(self):
         dims = [d for d, t in self.dim_tags]
         return max(dims, default=1)
@@ -850,6 +944,75 @@ def modify_regions(self, format: str):
         self.regions = regions
         return self
 
+    def dt_equal(self, other: 'ObjectSet'):
+        """ Tests two ObjectSets equality over dimtags.
+            Dimtags does not need to be sorted.
+         """
+        assert other
+        return sorted(self.dim_tags) == sorted(other.dim_tags)
+
+    def dt_copy(self) -> 'ObjectSet':
+        """
+        Create a shallow copy of this ObjectSet with independent data lists.
+        :return: new ObjectSet.
+        """
+        result = ObjectSet(factory=self.factory, dim_tags=list(self.dim_tags), regions=list(self.regions))
+        result.mesh_step_size = list(self.mesh_step_size)
+        return result
+
+
+    def dt_intersection(self, *obj_list: 'ObjectSet') -> 'ObjectSet':
+        """
+        Create intersection self and given list of  ObjectSets (its dimtags) over dimtags.
+        :param obj_list: List of ObjectSets to be intersected.
+        :return: new ObjectSet
+        """
+        assert obj_list
+        dim_tags = set()
+        for item in obj_list:
+            if isinstance(item, ObjectSet):
+                dim_tags.update(item.dim_tags)
+                # for dt in item.dim_tags:
+                #     if dt in self.dim_tags:
+                #         dim_tags.add(dt)
+                #     # debug output:
+                #     else:
+                #         message = "Intersection skip dimtag: {} - not in self object.".format(dt)
+                #         print(message)
+            else:
+                raise Exception(f"group: Wrong argument of type {type(item)}, expecting ObjectSet..")
+
+        result = ObjectSet(factory=self.factory, dim_tags=[], regions=[])
+        for idx, dt in enumerate(self.dim_tags):
+            if dt in dim_tags:
+                result.dim_tags.append(dt)
+                result.regions.append(self.regions[idx])
+                result.mesh_step_size.append(self.mesh_step_size[idx])
+        return result
+
+    def dt_drop(self, *obj_list: 'ObjectSet') -> 'ObjectSet':
+        """
+        Drop any number of ObjectSets (its dimtags) from the self object.
+        If dimtags not found, it is skipped.
+        :param obj_list: List of ObjectSets which dimtags should be dropped.
+        :return: self
+        """
+        assert obj_list
+        for item in obj_list:
+            if isinstance(item, ObjectSet):
+                for dt in item.dim_tags:
+                    try:
+                        idx = self.dim_tags.index(dt)
+                        self.dim_tags.pop(idx)
+                        self.regions.pop(idx)
+                        self.mesh_step_size.pop(idx)
+                    except ValueError as err:
+                        message = "Drop skip dimtag: {} - not in self object.".format(dt)
+                        print(message)
+            else:
+                raise Exception(f"group: Wrong argument of type {type(item)}, expecting ObjectSet..")
+        return self
+
     def translate(self, vector):
         self.factory.model.translate(self.dim_tags, *vector)
         self.factory._need_synchronize = True
@@ -909,16 +1072,116 @@ def revolve(self, center, axis, angle, numElements=[], heights=[], recombine=Fal
         # split the Objectset by dimtags
         return all_obj.split_by_dimension()
 
+    def fillet(self, curve_obj, radii, removeVolume=True):
+        """
+        gmsh.model.occ.fillet(volumeTags, curveTags, radii, removeVolume=True)
+
+        Fillet the volumes `volumeTags' on the curves `curveTags' with `radii'.
+        It has effect on the edge curves of the volume rounding them with given radius.
+
+         The `radii' vector can either contain a single radius, as many
+        radii as `curveTags', or twice as many as `curveTags' (in which case
+        different radii are provided for the begin and end points of the curves).
+        Return the filleted entities in `outDimTags' as a vector of (dim, tag)
+        pairs. Remove the original volume if `removeVolume' is set.
+
+        Return `outDimTags'.
+
+        Types:
+        - `volumeTags': vector of integers
+        - `curveTags': vector of integers
+        - `radii': vector of doubles
+        - `outDimTags': vector of pairs of integers
+        - `removeVolume': boolean
+        """
+        assert self.min_dim() == self.max_dim() == 3
+        volume_tags = self.tags
+        assert curve_obj.min_dim() == curve_obj.max_dim() == 1
+        curve_tags = curve_obj.tags
+        radii = np.atleast_1d(radii)
+        try:
+            outDimTags = self.factory.model.fillet(volume_tags, curve_tags, radii, removeVolume)
+        except ValueError as err:
+            message = "\nFillet failed!\ndimtags: {}".format(str(self.dim_tags[:10]))
+            gerr = gmsh_exceptions.BoolOperationError(message)
+            self._raise_gmsh_exception(gerr, err)
+
+        regions = [Region.default_region[dim] for dim, tag in outDimTags]
+        all_obj = ObjectSet(self.factory, outDimTags, regions)
+
+        self.factory._need_synchronize = True
+        # split the Objectset by dimtags
+        return all_obj
+
+        def chamfer(volumeTags, curveTags, surfaceTags, distances, removeVolume=True):
+            """
+            gmsh.model.occ.chamfer(volumeTags, curveTags, surfaceTags, distances, removeVolume=True)
+
+            Chamfer the volumes `volumeTags' on the curves `curveTags' with distances
+            `distances' measured on surfaces `surfaceTags'. The `distances' vector can
+            either contain a single distance, as many distances as `curveTags' and
+            `surfaceTags', or twice as many as `curveTags' and `surfaceTags' (in which
+            case the first in each pair is measured on the corresponding surface in
+            `surfaceTags', the other on the other adjacent surface). Return the
+            chamfered entities in `outDimTags'. Remove the original volume if
+            `removeVolume' is set.
+
+            Return `outDimTags'.
+
+            Types:
+            - `volumeTags': vector of integers
+            - `curveTags': vector of integers
+            - `surfaceTags': vector of integers
+            - `distances': vector of doubles
+            - `outDimTags': vector of pairs of integers
+            - `removeVolume': boolean
+            """
+            pass
+
+    def get_mass(self):
+        """
+        Total mass of the object dim tags, no sense for mixed dimensions
+        :return:
+        """
+        return sum(self.factory.model.getMass(dim, tag) for dim, tag in self.dim_tags)
+
     def copy(self) -> 'ObjectSet':
+        dim_lists = self.split_by_dimension()
+        copy_dimensions = [
+            dim_obj._copy_simple()
+            for dim_obj in dim_lists
+        ]
+        return ObjectSet.group(*copy_dimensions)
+    def _copy_simple(self) -> 'ObjectSet':
         """
         Problem: gmsh.model.occ.copy fails to copy boundary dimtags.
+        Problem: it seems to copy only highest dimension tags even if the lower dim tags are not part of the boundary.
         """
         copy_tags = self.factory.model.copy(self.dim_tags)
         self.factory._need_synchronize = True
-        copy_obj = ObjectSet(self.factory, copy_tags, self.regions)
+        copy_regions = self.regions
+        copy_obj = ObjectSet(self.factory, copy_tags, copy_regions)
         copy_obj.mesh_step_size = self.mesh_step_size.copy()
         return copy_obj
 
+    def unique(self):
+        """
+        Unique dim_tag list.
+        :return:
+        """
+        seen = set()
+        dimtags = []
+        regions = []
+
+        # remove duplicate dimtags
+        for dt, r in zip(self.dim_tags, self.regions):
+            if dt not in seen:
+                seen.add(dt)
+                dimtags.append(dt)
+                regions.append(r)
+
+        return ObjectSet(self.factory, dimtags, regions)
+
     def get_boundary(self, combined=False):
         """
         Get the boundary of the model entities dimTags.
@@ -1228,6 +1491,40 @@ def remove_small_mass(self, mass_limit):
         self.regions = regions
         return self
 
+def _transfrom_op(op, obj_list, *args, **kwargs):
+    if isinstance(obj_list, (list, tuple, set)):
+        type_constructor = type(obj_list)
+        return type_constructor(
+            _transfrom_op(op, item, *args, **kwargs)
+            for item in obj_list)
+    elif isinstance(obj_list, dict):
+        return {
+            k: _transfrom_op(op, item, *args, **kwargs)
+            for k, item in obj_list
+        }
+    else:
+        return op(obj_list, *args, **kwargs)
+
+def _apply_recursive(fn):
+    def recursive_fn(data, *args, **kwargs):
+        return _transfrom_op(fn, data, *args, **kwargs)
+    return recursive_fn
+# Apply transformop to all ObjectSet methods
+
+for name, method in inspect.getmembers(ObjectSet, predicate=inspect.isfunction):
+    globals()[name] = _apply_recursive(method)
+
+
+# def translate(obj_list, vector):
+#     return _transfrom_op(ObjectSet.translate, obj_list, vector)
+#
+# def rotate(obj_list, axis, angle, center=[0, 0, 0]):
+#     return _transfrom_op(ObjectSet.rotate, obj_list, axis, angle, center)
+#
+# def scale(obj_list, scale_vector, center=[0, 0, 0]):
+#     return _transfrom_op(ObjectSet.scale, obj_list, scale_vector, center)
+
+
 
 def rotation_matrix(axis, theta):
     """
@@ -1242,4 +1539,4 @@ def rotation_matrix(axis, theta):
     bc, ad, ac, ab, bd, cd = b * c, a * d, a * c, a * b, b * d, c * d
     return np.array([[aa + bb - cc - dd, 2 * (bc + ad), 2 * (bd - ac)],
                      [2 * (bc - ad), aa + cc - bb - dd, 2 * (cd + ab)],
-                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
\ No newline at end of file
+                     [2 * (bd + ac), 2 * (cd - ab), aa + dd - bb - cc]])
diff --git a/src/bgem/gmsh/heal_mesh.py b/src/bgem/gmsh/heal_mesh.py
index 20134b7..9ed23e9 100644
--- a/src/bgem/gmsh/heal_mesh.py
+++ b/src/bgem/gmsh/heal_mesh.py
@@ -101,9 +101,14 @@ def gmsh_gamma(self):
             # circum radius
             a,b,c,A,B,C = self.edge_lens
             x_area = (a*A + b*B + c*C)*(a*A + b*B - c*C)*(a*A - b*B + c*C)*(-a*A + b*B + c*C)
-            assert x_area > 0, x_area
-            R = np.sqrt(x_area) / 24 / max(1e-300, V)
-            self._gamma = 3 * r/max(1e-300, R)
+            # assert x_area > 0, x_area
+            double_lim = 1e-300
+            if x_area < double_lim or V < double_lim:
+                print(f"gamma=0: z_area={x_area}, V={V}, double_lim={double_lim} (test 'x_area < double_lim or V < double_lim')")
+                self._gamma = 0
+            else:
+                R = np.sqrt(x_area) / 24 / max(double_lim, V)
+                self._gamma = 3 * r/max(double_lim, R)
         return self._gamma
 
 
diff --git a/src/bgem/gmsh/options.py b/src/bgem/gmsh/options.py
index 1c55052..c835daf 100644
--- a/src/bgem/gmsh/options.py
+++ b/src/bgem/gmsh/options.py
@@ -169,6 +169,8 @@ def __init__(self):
         self.RandomFactor = 1e-9
         # Random factor used in the 2D meshing algorithm (should be increased if RandomFactor * size(triangle)/size(model) approaches machine accuracy)
         self.RandomFactor3D = 1e-12
+        # Seed of pseudo-random number generator
+        self.RandomSeed = 1
         # Random factor used in the 3D meshing algorithm
         self.finish_init()
 
diff --git a/src/bgem/stochastic/fracture.py b/src/bgem/stochastic/fracture.py
index a915514..e78bee8 100644
--- a/src/bgem/stochastic/fracture.py
+++ b/src/bgem/stochastic/fracture.py
@@ -528,7 +528,7 @@ def copy_full_range(self):
     @classmethod
     def from_mean_area(cls, power, diam_range, p32, p32_power=None):
         """
-        Construct the distribution using the mean arrea (P32) instead of intensity.
+        Construct the distribution using the mean area (P32) instead of intensity.
         :param p32: mean area of the fractures in given `diam_range`.
         :param p32_power: if the mean area is given for different power parameter.
         :return: PowerLawSize instance.
@@ -842,21 +842,30 @@ class FrFamily:
     size: PowerLawSize
     shape_angle: VonMisesOrientation
 
+    @staticmethod
+    def _create_power_law_size(family):
+        size_range = (family["r_min"], family["r_max"])
+        if "p_32" in family:
+            power_law_size = PowerLawSize.from_mean_area(family["power"], size_range, family["p_32"])
+            assert np.isclose(family["p_32"], power_law_size.mean_area())
+            return power_law_size
+        elif "p_30" in family:
+            power_law_size = PowerLawSize(family["power"], size_range, family["p_30"])
+            return power_law_size
+        else:
+            raise Exception("Found no existing fracture power intensity (p_30/p_32)")
+
     @staticmethod
     def from_cfg_3d(family):
         fisher_orientation = FisherOrientation(family["trend"], family["plunge"], family["concentration"])
-        size_range = (family["r_min"], family["r_max"])
-        power_law_size = PowerLawSize.from_mean_area(family["power"], size_range, family["p_32"])
-        assert np.isclose(family["p_32"], power_law_size.mean_area())
+        power_law_size = FrFamily._create_power_law_size(family)
         shape_angle = VonMisesOrientation(trend=0, concentration=0)
         return FrFamily(family["name"], fisher_orientation, power_law_size, shape_angle)
 
     @staticmethod
     def from_cfg_2d(family):
         orientation = FisherOrientation(0, 90, np.inf)
-        size_range = (family["r_min"], family["r_max"])
-        power_law_size = PowerLawSize.from_mean_area(family["power"], size_range, family["p_32"])
-        assert np.isclose(family["p_32"], power_law_size.mean_area())
+        power_law_size = FrFamily._create_power_law_size(family)
         shape_angle = VonMisesOrientation(family["trend"], family["concentration"])
         return FrFamily(family["name"], orientation, power_law_size, shape_angle)
 
@@ -877,7 +886,7 @@ class Population:
 
     @property
     def volume(self):
-        return np.product([l if l>0 else 1.0 for l in self.domain])
+        return np.prod([l if l>0 else 1.0 for l in self.domain])
 
     @classmethod
     def initialize_3d(cls, families: List[Dict[str, Any]], box):
diff --git a/tests/gmsh/test_gmsh_dimtags.py b/tests/gmsh/test_gmsh_dimtags.py
new file mode 100644
index 0000000..238b082
--- /dev/null
+++ b/tests/gmsh/test_gmsh_dimtags.py
@@ -0,0 +1,96 @@
+import pytest
+
+from bgem.gmsh.gmsh import ObjectSet, Region
+
+
+def make_region(dim, tag):
+    return Region(dim=dim, id=tag, name="region_{}".format(tag))
+
+
+def make_object_set(factory, dim_tags):
+    regions = [make_region(dim, tag) for dim, tag in dim_tags]
+    obj = ObjectSet(factory=factory, dim_tags=list(dim_tags), regions=regions)
+    obj.mesh_step_size = [10 + i for i in range(len(dim_tags))]
+    return obj
+
+
+def assert_object_set_data(obj, dim_tags, regions, mesh_step_size):
+    assert obj.dim_tags == dim_tags
+    assert obj.regions == regions
+    assert obj.mesh_step_size == mesh_step_size
+
+
+def test_dt_equal_compares_unsorted_dimtags():
+    factory = object()
+    obj = make_object_set(factory, [(2, 20), (1, 10), (3, 30)])
+    same_dimtags = make_object_set(factory, [(3, 30), (2, 20), (1, 10)])
+    different_dimtags = make_object_set(factory, [(3, 30), (2, 20), (0, 1)])
+
+    assert obj.dt_equal(same_dimtags)
+    assert not obj.dt_equal(different_dimtags)
+
+
+def test_dt_copy_preserves_data_without_sharing_lists():
+    factory = object()
+    obj = make_object_set(factory, [(2, 20), (1, 10), (3, 30)])
+
+    copied = obj.dt_copy()
+
+    assert copied is not obj
+    assert_object_set_data(copied, obj.dim_tags, obj.regions, obj.mesh_step_size)
+
+    copied.dim_tags.pop()
+    copied.regions.pop()
+    copied.mesh_step_size.pop()
+
+    assert obj.dim_tags == [(2, 20), (1, 10), (3, 30)]
+    assert len(obj.regions) == 3
+    assert obj.mesh_step_size == [10, 11, 12]
+
+
+def test_dt_intersection_filters_self_by_other_dimtags_preserving_self_data():
+    factory = object()
+    obj = make_object_set(factory, [(2, 20), (1, 10), (3, 30), (2, 21)])
+    other_a = make_object_set(factory, [(3, 30), (0, 5)])
+    other_b = make_object_set(factory, [(2, 20), (3, 30)])
+
+    result = obj.dt_intersection(other_a, other_b)
+
+    assert result.factory is factory
+    assert_object_set_data(
+        result,
+        [(2, 20), (3, 30)],
+        [obj.regions[0], obj.regions[2]],
+        [10, 12],
+    )
+
+
+def test_dt_intersection_rejects_non_object_set_arguments():
+    obj = make_object_set(object(), [(2, 20)])
+
+    with pytest.raises(Exception, match="expecting ObjectSet"):
+        obj.dt_intersection([(2, 20)])
+
+
+def test_dt_drop_removes_all_matching_dimtags_in_place():
+    factory = object()
+    obj = make_object_set(factory, [(2, 20), (1, 10), (3, 30), (2, 21)])
+    drop_a = make_object_set(factory, [(1, 10), (0, 5)])
+    drop_b = make_object_set(factory, [(2, 21)])
+
+    result = obj.dt_drop(drop_a, drop_b)
+
+    assert result is obj
+    assert_object_set_data(
+        obj,
+        [(2, 20), (3, 30)],
+        [make_region(2, 20), make_region(3, 30)],
+        [10, 12],
+    )
+
+
+def test_dt_drop_rejects_non_object_set_arguments():
+    obj = make_object_set(object(), [(2, 20)])
+
+    with pytest.raises(Exception, match="expecting ObjectSet"):
+        obj.dt_drop([(2, 20)])