diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2208a9a3..d6f4d643 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -20,6 +20,15 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Added `compas_model.elements.Element.parent` as alias for `compas_model.elements.Element.tree_node.parent`.
 * Added missing graph node reference to elements during deserialisation process.
+* Added a base `BlockModel`.
+* Added reference to model `Element.model` to `Element`.
+* Added `Element.modelgeometry` as the cached geometry of an element in model coordinates, taking into account the modifying effect of interactions with other elements.
+* Added `Element.modeltransformation` as the cached transformation from element to model coordinates.
+* Added `Element.compute_elementgeometry()`.
+* Added `Element.compute_modelgeometry()` to replace `Element.compute_geometry()`.
+* Added `Element.compute_modeltransformation()` to replace `Element.compute_worldtransformation()`.
+* Added `Element.is_dirty` that is changed together with neighbor elements.  
+* Added tests for element attributes.
 
 ### Changed
 
@@ -29,11 +38,16 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Changed `compas_model.elements.Element.compute_worldtransformation` to include frame of model.
 * Changed `compas_model.models.elementnode.ElementNode` to include children (previous functionality of `GroupNode`).
 * Changed root of element tree to element node instead of group node.
-* Changed deserialisation process of model according to removal of group node.
+* Changed deserialisation process of model according to removal of group node.\
+* Changed `Element.graph_node` to `Element.graphnode`.
+* Changed `Element.tree_node` to `Element.treenode`.
+* Changed `blockmodel_interfaces` to use the bestfit frame shared by two aligned interfaces instead of the frame of first face of the pair.
 
 ### Removed
 
 * Removed `compas_model.models.groupnode.GroupNode`.
+* Removed model reference `ElementTree.model` from `ElementTree`.
+* Removed `InterfaceElement` from elements.
 
 ## [0.4.4] 2024-06-13
 
diff --git a/docs/examples/grid_model/000_one_unit.py b/docs/examples/grid_model/000_one_unit.py
new file mode 100644
index 00000000..a459bf6e
--- /dev/null
+++ b/docs/examples/grid_model/000_one_unit.py
@@ -0,0 +1,5 @@
+# Example file to do the following:
+# Create a model with 4 columns, 4 column-heads and 1 plate.
+# Compute interactions between the elements. TODO: Add interactions and the base element level when geometry is computed.
+# Change one column head and check if the elements are dirty.
+# Recompute the modelgeometry.
\ No newline at end of file
diff --git a/scripts/test_blockmodel_arch.py b/scripts/test_blockmodel_arch.py
index af2c9d70..f8ab1c1f 100644
--- a/scripts/test_blockmodel_arch.py
+++ b/scripts/test_blockmodel_arch.py
@@ -3,18 +3,19 @@
 import compas
 from compas.colors import Color
 from compas_assembly.geometry import Arch
+from compas_viewer import Viewer
+
 from compas_model.algorithms import blockmodel_interfaces
 from compas_model.analysis import cra_penalty_solve
 from compas_model.elements import BlockElement
 from compas_model.interactions import ContactInterface
 from compas_model.models import Model
-from compas_viewer import Viewer
 
 # =============================================================================
 # Block model
 # =============================================================================
 
-template = Arch(rise=3, span=10, thickness=0.2, depth=0.5, n=30)
+template = Arch(rise=3, span=10, thickness=0.2, depth=0.5, n=200)
 
 model = Model()
 
@@ -31,7 +32,7 @@
 # Equilibrium
 # =============================================================================
 
-elements: list[BlockElement] = sorted(model.elements(), key=lambda e: e.geometry.centroid().z)[:2]
+elements: list[BlockElement] = sorted(model.elements(), key=lambda e: e.modelgeometry.centroid().z)[:2]
 
 for element in elements:
     element.is_support = True
@@ -61,7 +62,8 @@
         color = Color(0.8, 0.8, 0.8)
         show_faces = False
 
-    viewer.scene.add(element.geometry, show_points=False, show_faces=show_faces, facecolor=color)
+    viewer.scene.add(element.modelgeometry, show_points=False, show_faces=show_faces, facecolor=color)
+
 
 for interaction in model.interactions():
     interaction: ContactInterface
diff --git a/src/compas_model/algorithms/interfaces.py b/src/compas_model/algorithms/interfaces.py
index a2998090..3dc40ac1 100644
--- a/src/compas_model/algorithms/interfaces.py
+++ b/src/compas_model/algorithms/interfaces.py
@@ -6,14 +6,16 @@
 from compas.geometry import Plane
 from compas.geometry import Polygon
 from compas.geometry import Transformation
+from compas.geometry import Vector
+from compas.geometry import bestfit_frame_numpy
 from compas.geometry import centroid_polygon
 from compas.geometry import is_colinear
 from compas.geometry import is_coplanar
+from compas.geometry import is_parallel_vector_vector
 from compas.geometry import transform_points
 from compas.itertools import window
 from shapely.geometry import Polygon as ShapelyPolygon
 
-# from compas_model.elements import BlockElement
 from compas_model.elements import BlockGeometry
 from compas_model.interactions import ContactInterface
 from compas_model.models import Model
@@ -21,6 +23,14 @@
 from .nnbrs import find_nearest_neighbours
 
 
+def invert(self):
+    self._yaxis = self._yaxis * -1
+    self._zaxis = self._zaxis * -1
+
+
+Frame.invert = invert
+
+
 def blockmodel_interfaces(
     model: Model,
     nmax: int = 10,
@@ -49,7 +59,7 @@ def blockmodel_interfaces(
     node_index = {node: index for index, node in enumerate(model.graph.nodes())}
     index_node = {index: node for index, node in enumerate(model.graph.nodes())}
 
-    blocks: List[BlockGeometry] = [model.graph.node_element(node).geometry for node in model.graph.nodes()]
+    blocks: List[BlockGeometry] = [model.graph.node_element(node).modelgeometry for node in model.graph.nodes()]
 
     nmax = min(nmax, len(blocks))
 
@@ -106,27 +116,51 @@ def mesh_mesh_interfaces(
     -------
     List[:class:`ContactInterface`]
 
+    Notes
+    -----
+    For equilibrium calculations with CRA, it is important that interface frames are aligned
+    with the direction of the (interaction) edges on which they are stored.
+
+    This means that if the
+
     """
     world = Frame.worldXY()
     interfaces = []
-    frames = a.frames()
+    # frames = a.frames()
 
     for face in a.faces():
-        points = a.face_coordinates(face)
-        frame = frames[face]
-        matrix = Transformation.from_change_of_basis(world, frame)
-        projected = transform_points(points, matrix)
-        p0 = ShapelyPolygon(projected)
+        a_points = a.face_coordinates(face)
+        a_normal = a.face_normal(face)
 
         for test in b.faces():
-            points = b.face_coordinates(test)
-            projected = transform_points(points, matrix)
-            p1 = ShapelyPolygon(projected)
+            b_points = b.face_coordinates(test)
+            b_normal: Vector = b.face_normal(test)
+
+            if not is_parallel_vector_vector(a_normal, b_normal):
+                continue
+
+            # this ensures that a shared frame is used to do the interface calculations
+            # the frame should be oriented along the normal of the "a" face
+            # this will align the interface frame with the resulting interaction edge
+            # whgich is important for calculations with solvers such as CRA
+            frame = Frame(*bestfit_frame_numpy(a_points + b_points))
+            if frame.zaxis.dot(a_normal) < 0:
+                frame.invert()
+
+            matrix = Transformation.from_change_of_basis(world, frame)
+
+            a_projected = transform_points(a_points, matrix)
+            p0 = ShapelyPolygon(a_projected)
+
+            b_projected = transform_points(b_points, matrix)
+            p1 = ShapelyPolygon(b_projected)
+
+            projected = a_projected + b_projected
 
             if not all(fabs(point[2]) < tmax for point in projected):
                 continue
 
-            if p1.area < amin:
+            if p0.area < amin or p1.area < amin:
                 continue
 
             if not p0.intersects(p1):
diff --git a/src/compas_model/analysis/cra.py b/src/compas_model/analysis/cra.py
index b26434e7..c5bc2b43 100644
--- a/src/compas_model/analysis/cra.py
+++ b/src/compas_model/analysis/cra.py
@@ -20,10 +20,10 @@ def cra_penalty_solve(
     element_block = {}
 
     for element in model.elements():
-        block: Block = element.geometry.copy(cls=Block)
+        block: Block = element.modelgeometry.copy(cls=Block)
         x, y, z = block.centroid()
         node = assembly.add_block(block, x=x, y=y, z=z, is_support=element.is_support)
-        element_block[element.graph_node] = node
+        element_block[element.graphnode] = node
 
     for edge in model.graph.edges():
         interactions: list[ContactInterface] = model.graph.edge_interactions(edge)
diff --git a/src/compas_model/elements/__init__.py b/src/compas_model/elements/__init__.py
index 03b18ef2..d150ec1a 100644
--- a/src/compas_model/elements/__init__.py
+++ b/src/compas_model/elements/__init__.py
@@ -4,8 +4,6 @@
 from .block import BlockElement
 from .block import BlockFeature
 from .block import BlockGeometry
-from .interface import InterfaceElement
-from .interface import InterfaceFeature
 from .plate import PlateElement
 from .plate import PlateFeature
 
@@ -17,8 +15,6 @@
     "BlockElement",
     "BlockFeature",
     "BlockGeometry",
-    "InterfaceElement",
-    "InterfaceFeature",
     "PlateElement",
     "PlateFeature",
 ]
diff --git a/src/compas_model/elements/block.py b/src/compas_model/elements/block.py
index a35b51a7..cb1e7189 100644
--- a/src/compas_model/elements/block.py
+++ b/src/compas_model/elements/block.py
@@ -3,10 +3,10 @@
 from compas.geometry import Box
 from compas.geometry import Frame
 from compas.geometry import Point
+from compas.geometry import bestfit_frame_numpy
 from compas.geometry import bounding_box
 from compas.geometry import centroid_points
 from compas.geometry import centroid_polyhedron
-from compas.geometry import cross_vectors
 from compas.geometry import dot_vectors
 from compas.geometry import oriented_bounding_box
 from compas.geometry import volume_polyhedron
@@ -15,30 +15,16 @@
 from compas_model.elements import Feature
 
 
-class BlockGeometry(Mesh):
-    def centroid(self):
-        """Compute the centroid of the block.
-
-        Returns
-        -------
-        :class:`compas.geometry.Point`
-
-        """
-        x, y, z = centroid_points([self.vertex_coordinates(key) for key in self.vertices()])
-        return Point(x, y, z)
+def invert(self):
+    self._yaxis = self._yaxis * -1
+    self._zaxis = self._zaxis * -1
 
-    def frames(self):
-        """Compute the local frame of each face of the block.
 
-        Returns
-        -------
-        dict
-            A dictionary mapping face identifiers to face frames.
+Frame.invert = invert
 
-        """
-        return {face: self.frame(face) for face in self.faces()}
 
-    def frame(self, face):
+class BlockGeometry(Mesh):
+    def face_frame(self, face):
         """Compute the frame of a specific face.
 
         Parameters
@@ -52,11 +38,12 @@ def frame(self, face):
 
         """
         xyz = self.face_coordinates(face)
-        o = self.face_center(face)
-        w = self.face_normal(face)
-        u = [xyz[1][i] - xyz[0][i] for i in range(3)]  # align with longest edge instead?
-        v = cross_vectors(w, u)
-        return Frame(o, u, v)
+        normal = self.face_normal(face)
+        o, u, v = bestfit_frame_numpy(xyz)
+        frame = Frame(o, u, v)
+        if frame.zaxis.dot(normal) < 0:
+            frame.invert()
+        return frame
 
     def top(self):
         """Identify the *top* face of the block.
@@ -72,6 +59,17 @@ def top(self):
         normals = [self.face_normal(face) for face in faces]
         return sorted(zip(faces, normals), key=lambda x: dot_vectors(x[1], z))[-1][0]
 
+    def centroid(self):
+        """Compute the centroid of the block.
+
+        Returns
+        -------
+        :class:`compas.geometry.Point`
+
+        """
+        x, y, z = centroid_points([self.vertex_coordinates(key) for key in self.vertices()])
+        return Point(x, y, z)
+
     def center(self):
         """Compute the center of mass of the block.
 
@@ -138,10 +136,12 @@ class BlockElement(Element):
 
     """
 
+    elementgeometry: BlockGeometry
+    modelgeometry: BlockGeometry
+
     @property
-    def __data__(self):
-        # type: () -> dict
-        data = super(BlockElement, self).__data__
+    def __data__(self) -> dict:
+        data = super().__data__
         data["shape"] = self.shape
         data["features"] = self.features
         data["is_support"] = self.is_support
@@ -149,34 +149,46 @@ def __data__(self):
 
     def __init__(self, shape, features=None, is_support=False, frame=None, transformation=None, name=None):
         # type: (Mesh | BlockGeometry, list[BlockFeature] | None, bool, compas.geometry.Frame | None, compas.geometry.Transformation | None, str | None) -> None
+        super().__init__(frame=frame, transformation=transformation, name=name)
 
-        super(BlockElement, self).__init__(frame=frame, transformation=transformation, name=name)
         self.shape = shape if isinstance(shape, BlockGeometry) else shape.copy(cls=BlockGeometry)
         self.features = features or []  # type: list[BlockFeature]
         self.is_support = is_support
 
-    # don't like this
-    # but want to test the collider
-    @property
-    def face_polygons(self):
-        # type: () -> list[compas.geometry.Polygon]
-        return [self.geometry.face_polygon(face) for face in self.geometry.faces()]  # type: ignore
-
     # =============================================================================
     # Implementations of abstract methods
     # =============================================================================
 
-    def compute_geometry(self, include_features=False):
+    def compute_elementgeometry(self):
         geometry = self.shape
-        if include_features:
-            if self.features:
-                for feature in self.features:
-                    geometry = feature.apply(geometry)
-        geometry.transform(self.worldtransformation)
+        # apply features?
+        return geometry
+
+    def compute_modelgeometry(self):
+        if not self.model:
+            raise Exception
+
+        geometry = self.elementgeometry.transformed(self.modeltransformation)
+
+        # apply effect of interactions?
+        node = self.graphnode
+        nbrs = self.model.graph.neighbors_in(node)
+        for nbr in nbrs:
+            element = self.model.graph.node_element(nbr)
+            if element:
+                for interaction in self.model.graph.edge_interactions((nbr, node)):
+                    # example interactions are
+                    # cutters, boolean operations, slicers, ...
+                    if hasattr(interaction, "apply"):
+                        try:
+                            interaction.apply(geometry)
+                        except Exception:
+                            pass
+
         return geometry
 
     def compute_aabb(self, inflate=0.0):
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
+        points = self.modelgeometry.vertices_attributes("xyz")  # type: ignore
         box = Box.from_bounding_box(bounding_box(points))
         box.xsize += inflate
         box.ysize += inflate
@@ -184,7 +196,7 @@ def compute_aabb(self, inflate=0.0):
         return box
 
     def compute_obb(self, inflate=0.0):
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
+        points = self.modelgeometry.vertices_attributes("xyz")  # type: ignore
         box = Box.from_bounding_box(oriented_bounding_box(points))
         box.xsize += inflate
         box.ysize += inflate
@@ -195,7 +207,7 @@ def compute_collision_mesh(self):
         # TODO: (TvM) make this a pluggable with default implementation in core and move import to top
         from compas.geometry import convex_hull_numpy
 
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
+        points = self.modelgeometry.vertices_attributes("xyz")  # type: ignore
         vertices, faces = convex_hull_numpy(points)
         vertices = [points[index] for index in vertices]  # type: ignore
         return Mesh.from_vertices_and_faces(vertices, faces)
diff --git a/src/compas_model/elements/element.py b/src/compas_model/elements/element.py
index 053fa779..0adf8649 100644
--- a/src/compas_model/elements/element.py
+++ b/src/compas_model/elements/element.py
@@ -1,30 +1,30 @@
-import compas
-
-if not compas.IPY:
-    from typing import TYPE_CHECKING
-    from typing import Union  # noqa: F401
-
-    if TYPE_CHECKING:
-        from compas_model.models import ElementNode  # noqa: F401
-
 from functools import reduce
 from functools import wraps
 from operator import mul
+from typing import TYPE_CHECKING
+from typing import Optional
+from typing import Union
 
-import compas.datastructures  # noqa: F401
+import compas
+import compas.datastructures
 import compas.geometry
 from compas.data import Data
 from compas.geometry import Transformation
 
+if TYPE_CHECKING:
+    pass
+
 
 def reset_computed(f):
     @wraps(f)
     def wrapper(*args, **kwargs):
-        self = args[0]
+        self: Element = args[0]
         self._aabb = None
         self._obb = None
         self._collision_mesh = None
         self._geometry = None
+        self._modelgeometry = None
+        self._modeltransformation = None
         return f(*args, **kwargs)
 
     return wrapper
@@ -41,16 +41,14 @@ class Feature(Data):
     """
 
     @property
-    def __data__(self):
-        # type: () -> dict
+    def __data__(self) -> dict:
         return {}
 
-    def __init__(self, name=None):
-        # type: (str | None) -> None
-        super(Feature, self).__init__(name=name)
+    def __init__(self, name: Optional[str] = None) -> None:
+        super().__init__(name=name)
 
-    def apply(self, shape):
-        # type: (compas.datastructures.Mesh | compas.geometry.Brep) -> compas.datastructures.Mesh | compas.geometry.Brep
+    def apply(self, shape: Union[compas.datastructures.Mesh, compas.geometry.Brep]) -> Union[compas.datastructures.Mesh, compas.geometry.Brep]:
+        """Apply the feature to the given shape, which represents the base geometry of the host element of the feature."""
         raise NotImplementedError
 
 
@@ -68,9 +66,9 @@ class Element(Data):
 
     Attributes
     ----------
-    graph_node : int
+    graphnode : int
         The identifier of the corresponding node in the interaction graph of the parent model.
-    tree_node : :class:`compas.datastructures.TreeNode`
+    treenode : :class:`compas.datastructures.TreeNode`
         The node in the hierarchical element tree of the parent model.
     frame : :class:`compas.geometry.Frame`
         The local coordinate frame of the element.
@@ -94,8 +92,7 @@ class Element(Data):
     """
 
     @property
-    def __data__(self):
-        # type: () -> dict
+    def __data__(self) -> dict:
         # note that the material can/should not be added here,
         # because materials should be added by/in the context of a model
         # and becaue this would also require a custom "from_data" classmethod.
@@ -110,24 +107,33 @@ def __init__(
         geometry=None,  # type: compas.geometry.Shape | compas.geometry.Brep | compas.datastructures.Mesh | None
         frame=None,  # type: compas.geometry.Frame | None
         transformation=None,  # type: compas.geometry.Transformation | None
+        features=None,  # type: list[Feature]
         name=None,  # type: str | None
     ):  # type: (...) -> None
         super(Element, self).__init__(name=name)
-        self.graph_node = None  # type: int | None
-        self.tree_node = None  # type: ElementNode | None
-        self._aabb = None
-        self._obb = None
-        self._collision_mesh = None
-        self._geometry = geometry
+        self.model = None  # type: Model | None
+
+        self.graphnode = None  # type: int | None
+        self.treenode = None  # type: ElementNode | None
+
         self._frame = frame
         self._transformation = transformation
-        self._worldtransformation = None
+        self._geometry = geometry  # this is same as elementgeometry
+        self._features = features or []  # type: list[Feature]
         self._material = None
-        self.features = []  # type: list[Feature]
+
+        self._aabb = None
+        self._obb = None
+        self._collision_mesh = None
+        self._modelgeometry = None
+        self._modeltransformation = None
+
         self.include_features = False
         self.inflate_aabb = 0.0
         self.inflate_obb = 0.0
 
+        self._is_dirty = True
+
     # this is not entirely correct
     def __repr__(self):
         # type: () -> str
@@ -143,30 +149,18 @@ def frame(self):
         return self._frame
 
     @frame.setter
+    @reset_computed
     def frame(self, frame):
-        self._worldtransformation = None
-        self._aabb = None
-        self._obb = None
-        self._collision_mesh = None
-        self._geometry = None
         self._frame = frame
 
     @property
     def transformation(self):
         # type: () -> compas.geometry.Transformation | None
-        """
-        Transformation of the element wrt its own frame.
-        When this property is changed, all computed properties have to be recomputed.
-        """
         return self._transformation
 
     @transformation.setter
+    @reset_computed
     def transformation(self, transformation):
-        self._worldtransformation = None
-        self._aabb = None
-        self._obb = None
-        self._collision_mesh = None
-        self._geometry = None
         self._transformation = transformation
 
     @property
@@ -174,28 +168,42 @@ def material(self):
         return self._material
 
     @property
-    def parent(self):
-        return self.tree_node.parent
+    def is_dirty(self):
+        return self._is_dirty
+
+    @is_dirty.setter
+    def is_dirty(self, value):
+        self._is_dirty = value
+
+        if value:
+            elements = list(self.model.elements())
+            for neighbor in self.model.graph.neighbors_out(self.graphnode):
+                elements[neighbor].is_dirty = value
 
     # ==========================================================================
     # Computed attributes
     # ==========================================================================
 
-    # it might be easier to just always compute these
+    @property
+    def elementgeometry(self):
+        # type: () -> ...
+        if self._geometry is None:
+            self._geometry = self.compute_elementgeometry()
+        return self._geometry
 
     @property
-    def worldtransformation(self):
+    def modeltransformation(self):
         # type: () -> compas.geometry.Transformation
-        if self._worldtransformation is None:
-            self._worldtransformation = self.compute_worldtransformation()
-        return self._worldtransformation
+        if self._modeltransformation is None:
+            self._modeltransformation = self.compute_modeltransformation()
+        return self._modeltransformation
 
     @property
-    def geometry(self):
+    def modelgeometry(self):
         # type: () -> ...
-        if self._geometry is None:
-            self._geometry = self.compute_geometry()
-        return self._geometry
+        if self._modelgeometry is None:
+            self._modelgeometry = self.compute_modelgeometry()
+        return self._modelgeometry
 
     @property
     def aabb(self):
@@ -227,9 +235,34 @@ def collision_mesh(self):
     # Abstract methods
     # ==========================================================================
 
-    def compute_worldtransformation(self):
+    def compute_elementgeometry(self):
+        # type: () -> compas.datastructures.Mesh | compas.geometry.Brep
+        """Compute the geometry of the element in local coordinates.
+
+        This is the parametric representation of the element,
+        without considering its location in the model or its interaction(s) with connected elements.
+
+        Returns
+        -------
+        :class:`compas.datastructures.Mesh` | :class:`compas.geometry.Brep`
+
+        """
+        raise NotImplementedError
+
+    def apply_features(self):
+        # type: () -> None
+        """Apply the features to the (base) geometry.
+
+        Returns
+        -------
+        None
+
+        """
+        raise NotImplementedError
+
+    def compute_modeltransformation(self):
         # type: () -> compas.geometry.Transformation
-        """Compute the transformation to world coordinates of this element
+        """Compute the transformation to model coordinates of this element
         based on its position in the spatial hierarchy of the model.
 
         Returns
@@ -250,28 +283,21 @@ def compute_worldtransformation(self):
                     frame_stack.append(parent.element.frame)
             parent = parent.parent
 
-        frame = self.tree_node.tree.model.frame
-        if frame:
-            frame_stack.append(frame)
-
         matrices = [Transformation.from_frame(f) for f in frame_stack]
 
         if matrices:
-            worldtransformation = reduce(mul, matrices[::-1])
+            modeltransformation = reduce(mul, matrices[::-1])
         else:
-            worldtransformation = Transformation()
+            modeltransformation = Transformation()
 
         if self.transformation:
-            worldtransformation = worldtransformation * self.transformation
+            modeltransformation = modeltransformation * self.transformation
 
-        return worldtransformation
+        return modeltransformation
 
-    def compute_geometry(self):
+    def compute_modelgeometry(self):
         # type: () -> compas.datastructures.Mesh | compas.geometry.Brep
-        """Compute the geometry of the element.
-
-        Implementations of this method should transform the geometry to world coordinates,
-        using `self.worldtransformation`.
+        """Compute the geometry of the element in model coordinates and taking into account the effect of interactions with connected elements.
 
         Returns
         -------
@@ -282,9 +308,7 @@ def compute_geometry(self):
 
     def compute_aabb(self):
         # type: (float | None) -> compas.geometry.Box
-        """Computes the Axis Aligned Bounding Box (AABB) of the element.
-
-        Implementations of this method should base the computation on the geometry in world coordinates.
+        """Computes the Axis Aligned Bounding Box (AABB) of the geometry of the element.
 
         Returns
         -------
@@ -296,9 +320,7 @@ def compute_aabb(self):
 
     def compute_obb(self):
         # type: (float | None) -> compas.geometry.Box
-        """Computes the Oriented Bounding Box (OBB) of the element.
-
-        Implementations of this method should base the computation on the geometry in world coordinates.
+        """Computes the Oriented Bounding Box (OBB) of the geometry of the element.
 
         Returns
         -------
@@ -310,9 +332,7 @@ def compute_obb(self):
 
     def compute_collision_mesh(self):
         # type: () -> compas.datastructures.Mesh
-        """Computes the collision geometry of the element.
-
-        Implementations of this method should base the computation on the geometry in world coordinates.
+        """Computes the collision geometry of the geometry of the element.
 
         Returns
         -------
@@ -378,4 +398,4 @@ def add_feature(self, feature):
         None
 
         """
-        self.features.append(feature)
+        self._features.append(feature)
diff --git a/src/compas_model/elements/interface.py b/src/compas_model/elements/interface.py
deleted file mode 100644
index 88869579..00000000
--- a/src/compas_model/elements/interface.py
+++ /dev/null
@@ -1,101 +0,0 @@
-import compas.datastructures  # noqa: F401
-from compas.datastructures import Mesh
-from compas.geometry import Box
-from compas.geometry import bounding_box
-from compas.geometry import oriented_bounding_box
-
-from compas_model.elements import Element
-from compas_model.elements import Feature
-
-
-class InterfaceFeature(Feature):
-    pass
-
-
-class InterfaceElement(Element):
-    """Class representing a phyisical interface between two other elements.
-
-    Parameters
-    ----------
-    polygon : :class:`compas.geometry.Polygon`
-        A polygon that represents the outline of the interface.
-    thickness : float
-        The thickness of the interface.
-    features : list[:class:`InterfaceFeature`], optional
-        Additional interface features.
-    name : str, optional
-        The name of the element.
-
-    Attributes
-    ----------
-    shape : :class:`compas.datastructure.Mesh`
-        The base shape of the interface.
-    features : list[:class:`BlockFeature`]
-        A list of additional interface features.
-
-    Notes
-    -----
-    The shape of the interface is calculated automatically from the input polygon and thickness.
-    The frame of the element is the frame of the polygon.
-
-    """
-
-    @property
-    def __data__(self):
-        # type: () -> dict
-        data = super(InterfaceElement, self).__data__
-        return data
-
-    def __init__(self, polygon, thickness, features=None, frame=None, name=None):
-        # type: (compas.geometry.Polygon, float, list[InterfaceFeature] | None, compas.geometry.Frame | None, str | None) -> None
-
-        frame = frame or polygon.frame
-
-        super(InterfaceElement, self).__init__(frame=frame, name=name)
-        self._polygon = polygon
-        self._thickness = thickness
-        self.shape = self.compute_shape()
-        self.features = features or []  # type: list[InterfaceFeature]
-
-    def compute_shape(self):
-        # type: () -> compas.datastructures.Mesh
-        mesh = Mesh.from_polygons([self._polygon])
-        return mesh.thickened(self._thickness)
-
-    # =============================================================================
-    # Implementations of abstract methods
-    # =============================================================================
-
-    def compute_geometry(self, include_features=False):
-        geometry = self.shape
-        if include_features:
-            if self.features:
-                for feature in self.features:
-                    geometry = feature.apply(geometry)
-        geometry.transform(self.worldtransformation)
-        return geometry
-
-    def compute_aabb(self, inflate=0.0):
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
-        box = Box.from_bounding_box(bounding_box(points))
-        box.xsize += inflate
-        box.ysize += inflate
-        box.zsize += inflate
-        return box
-
-    def compute_obb(self, inflate=0.0):
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
-        box = Box.from_bounding_box(oriented_bounding_box(points))
-        box.xsize += inflate
-        box.ysize += inflate
-        box.zsize += inflate
-        return box
-
-    def compute_collision_mesh(self):
-        # TODO: (TvM) make this a pluggable with default implementation in core and move import to top
-        from compas.geometry import convex_hull_numpy
-
-        points = self.geometry.vertices_attributes("xyz")  # type: ignore
-        vertices, faces = convex_hull_numpy(points)
-        vertices = [points[index] for index in vertices]  # type: ignore
-        return Mesh.from_vertices_and_faces(vertices, faces)
diff --git a/src/compas_model/models/__init__.py b/src/compas_model/models/__init__.py
index f652fff4..ca79a1c2 100644
--- a/src/compas_model/models/__init__.py
+++ b/src/compas_model/models/__init__.py
@@ -2,6 +2,7 @@
 from .elementtree import ElementTree
 from .interactiongraph import InteractionGraph
 from .model import Model
+from .blockmodel import BlockModel
 
 
 __all__ = [
@@ -9,4 +10,5 @@
     "ElementTree",
     "InteractionGraph",
     "Model",
+    "BlockModel",
 ]
diff --git a/src/compas_model/models/blockmodel.py b/src/compas_model/models/blockmodel.py
index e69de29b..383e1464 100644
--- a/src/compas_model/models/blockmodel.py
+++ b/src/compas_model/models/blockmodel.py
@@ -0,0 +1,78 @@
+from compas.geometry import Brep
+from compas.tolerance import Tolerance
+from compas_occ.brep import OCCBrepFace as BrepFace
+
+from compas_model.algorithms.nnbrs import find_nearest_neighbours
+from compas_model.elements import Element
+from compas_model.interactions import ContactInterface
+
+from .model import Model
+
+
+class BlockModel(Model):
+    def __init__(self, name=None):
+        super().__init__(name)
+
+    def find_intersections(self):
+        pass
+
+    def find_overlaps(self, deflection=None, tolerance=1, max_distance=50, min_area=0):
+        pass
+
+    def identify_interactions(self, deflection=None, tolerance=1, max_distance=50, min_area=0, nmax=10):
+        deflection = deflection or Tolerance().lineardeflection
+
+        node_index = {node: index for index, node in enumerate(self.graph.nodes())}
+        index_node = {index: node for index, node in enumerate(self.graph.nodes())}
+
+        geometries: list[Brep] = [self.graph.node_element(node).geometry for node in self.graph.nodes()]
+
+        cloud = [geometry.centroid for geometry in geometries]
+        nnbrs = find_nearest_neighbours(cloud, nmax, dims=3)
+
+        for u in self.graph.nodes():
+            A: Element = self.graph.node_element(u)  # type: ignore
+
+            i = node_index[u]
+            nbrs = nnbrs[i][1]
+
+            for j in nbrs:
+                v = index_node[j]
+
+                if u == v:
+                    continue
+
+                if self.graph.has_edge((u, v), directed=False):
+                    continue
+
+                B: Element = self.graph.node_element(v)  # type: ignore
+
+                faces_A, faces_B = A.geometry.overlap(B.geometry, deflection=deflection, tolerance=tolerance)  # type: ignore
+                faces_A: list[BrepFace]
+                faces_B: list[BrepFace]
+
+                if faces_A and faces_B:
+                    for face_A in faces_A:
+                        brep_A = Brep.from_brepfaces([face_A])
+
+                        if brep_A.area < min_area:
+                            continue
+
+                        for face_B in faces_B:
+                            brep_B = Brep.from_brepfaces([face_B])
+
+                            if brep_B.area < min_area:
+                                continue
+
+                            brep_C: Brep = Brep.from_boolean_intersection(brep_A, brep_B)
+
+                            if brep_C.area < min_area:
+                                continue
+
+                            poly_C = brep_C.to_polygons()[0]
+                            mesh_C = brep_C.to_tesselation()[0]
+
+                            interaction = ContactInterface(points=poly_C.points, mesh=mesh_C)
+
+                            # do something with the interactions
+                            self.add_interaction(A, B, interaction=interaction)
diff --git a/src/compas_model/models/elementnode.py b/src/compas_model/models/elementnode.py
index bdcec7ef..7f863332 100644
--- a/src/compas_model/models/elementnode.py
+++ b/src/compas_model/models/elementnode.py
@@ -1,7 +1,12 @@
+from typing import TYPE_CHECKING
+
 from compas.datastructures import TreeNode
 
 from compas_model.elements import Element  # noqa: F401
 
+if TYPE_CHECKING:
+    from .elementtree import ElementTree
+
 
 class ElementNode(TreeNode):
     """Class representing nodes containing elements in an element tree.
@@ -23,6 +28,8 @@ class ElementNode(TreeNode):
 
     """
 
+    tree: "ElementTree"
+
     @property
     def __data__(self):
         # type: () -> dict
@@ -39,7 +46,7 @@ def __init__(self, element=None, **kwargs):
         # type: (Element | None, str | None) -> None
         super(ElementNode, self).__init__(**kwargs)
         if element:
-            element.tree_node = self
+            element.treenode = self
         self.element = element
 
     def __getitem__(self, index):
diff --git a/src/compas_model/models/elementtree.py b/src/compas_model/models/elementtree.py
index 0d8d3484..4e97fad5 100644
--- a/src/compas_model/models/elementtree.py
+++ b/src/compas_model/models/elementtree.py
@@ -38,10 +38,9 @@ def __from_data__(cls, data):
         # type: (dict) -> ElementTree
         raise Exception("Serialisation outside model context not allowed.")
 
-    def __init__(self, model, name=None):
-        # type: (compas_model.models.Model, str | None) -> None
+    def __init__(self, name=None):
+        # type: (str | None) -> None
         super(ElementTree, self).__init__(name=name)
-        self.model = model
         root = ElementNode(name="root")
         self.add(root)
 
diff --git a/src/compas_model/models/interactiongraph.py b/src/compas_model/models/interactiongraph.py
index 48982a3b..1019973f 100644
--- a/src/compas_model/models/interactiongraph.py
+++ b/src/compas_model/models/interactiongraph.py
@@ -51,7 +51,7 @@ def __from_data__(cls, data, guid_element):
         for node, attr in graph.nodes(data=True):
             element = guid_element[attr["element"]]
             attr["element"] = element  # type: ignore
-            element.graph_node = node
+            element.graphnode = node
         return graph
 
     def copy(self):
diff --git a/src/compas_model/models/model.py b/src/compas_model/models/model.py
index 3f72d76e..1fffb7be 100644
--- a/src/compas_model/models/model.py
+++ b/src/compas_model/models/model.py
@@ -1,20 +1,18 @@
 from collections import OrderedDict
 from collections import deque
+from typing import Generator  # noqa: F401
+from typing import Type  # noqa: F401
 
 import compas
-
-if not compas.IPY:
-    from typing import Generator  # noqa: F401
-    from typing import Type  # noqa: F401
-
 import compas.datastructures  # noqa: F401
 import compas.geometry  # noqa: F401
 from compas.datastructures import Datastructure
 from compas.geometry import Frame
+from compas.geometry import Transformation
 
-from compas_model.elements import Element  # noqa: F401
+from compas_model.elements import Element
 from compas_model.interactions import Interaction  # noqa: F401
-from compas_model.materials import Material  # noqa: F401
+from compas_model.materials import Material
 
 from .elementnode import ElementNode
 from .elementtree import ElementTree
@@ -51,12 +49,13 @@ class Model(Datastructure):
     """
 
     @property
-    def __data__(self):
+    def __data__(self) -> dict:
         # in their data representation,
         # the element tree and the interaction graph
         # refer to model elements by their GUID, to avoid storing duplicate data representations of those elements
         # the elements are stored in a global list
         data = {
+            "frame": self.frame,
             "tree": self._tree.__data__,
             "graph": self._graph.__data__,
             "elements": list(self.elements()),
@@ -66,19 +65,17 @@ def __data__(self):
         return data
 
     @classmethod
-    def __from_data__(cls, data):
+    def __from_data__(cls, data: dict) -> "Model":
         model = cls()
         model._guid_material = {str(material.guid): material for material in data["materials"]}
         model._guid_element = {str(element.guid): element for element in data["elements"]}
 
         for e, m in data["element_material"].items():
-            element = model._guid_element[e]
-            material = model._guid_material[m]
+            element: Element = model._guid_element[e]
+            material: Material = model._guid_material[m]
             element._material = material
 
-        def add(nodedata, parentnode):
-            # type: (dict, ElementNode) -> None
-
+        def add(nodedata: dict, parentnode: ElementNode) -> None:
             if "children" in nodedata:
                 for childdata in nodedata["children"]:
                     guid = childdata["element"]
@@ -102,11 +99,12 @@ def add(nodedata, parentnode):
         return model
 
     def __init__(self, name=None):
-        super(Model, self).__init__(name=name)
+        super().__init__(name=name)
+
         self._frame = None
         self._guid_material = {}
         self._guid_element = OrderedDict()
-        self._tree = ElementTree(model=self)
+        self._tree = ElementTree()
         self._graph = InteractionGraph()
         self._graph.update_default_node_attributes(element=None)
         self._graph.update_default_edge_attributes(interactions=None)
@@ -140,13 +138,11 @@ def __str__(self):
     # such that the user can't make unintended and potentially breaking changes
 
     @property
-    def tree(self):
-        # type: () -> ElementTree
+    def tree(self) -> ElementTree:
         return self._tree
 
     @property
-    def graph(self):
-        # type: () -> InteractionGraph
+    def graph(self) -> InteractionGraph:
         return self._graph
 
     # A model should have a coordinate system.
@@ -158,22 +154,20 @@ def graph(self):
     # the model can compute the transformations of all of the elements in the tree.
 
     @property
-    def frame(self):
-        # type: () -> compas.geometry.Frame
+    def frame(self) -> Frame:
         if not self._frame:
             self._frame = Frame.worldXY()
         return self._frame
 
     @frame.setter
-    def frame(self, frame):
+    def frame(self, frame: Frame):
         self._frame = frame
 
     # =============================================================================
     # Datastructure "abstract" methods
     # =============================================================================
 
-    def transform(self, transformation):
-        # type: (compas.geometry.Transformation) -> None
+    def transform(self, transformation: Transformation) -> None:
         """Transform the model and all that it contains.
 
         Parameters
@@ -194,8 +188,7 @@ def transform(self, transformation):
     # Methods
     # =============================================================================
 
-    def has_element(self, element):
-        # type: (Element) -> bool
+    def has_element(self, element: Element) -> bool:
         """Returns True if the model contains the given element.
 
         Parameters
@@ -211,8 +204,7 @@ def has_element(self, element):
         guid = str(element.guid)
         return guid in self._guid_element
 
-    def has_interaction(self, a, b):
-        # type: (Element, Element) -> bool
+    def has_interaction(self, a: Element, b: Element) -> bool:
         """Returns True if two elements have an interaction set between them.
 
         Parameters
@@ -227,15 +219,14 @@ def has_interaction(self, a, b):
         bool
 
         """
-        edge = a.graph_node, b.graph_node
+        edge = a.graphnode, b.graphnode
         result = self.graph.has_edge(edge)
         if not result:
-            edge = b.graph_node, a.graph_node
+            edge = b.graphnode, a.graphnode
             result = self.graph.has_edge(edge)
         return result
 
-    def has_material(self, material):
-        # type: (Material) -> bool
+    def has_material(self, material: Material) -> bool:
         """Verify that the model contains a specific material.
 
         Parameters
@@ -284,15 +275,15 @@ def add_element(self, element, parent=None, material=None):
             raise Exception("Element already in the model.")
         self._guid_element[guid] = element
 
-        element.graph_node = self.graph.add_node(element=element)
+        element.graphnode = self.graph.add_node(element=element)
 
         if not parent:
             parent = self._tree.root
 
         if isinstance(parent, Element):
-            if parent.tree_node is None:
+            if parent.treenode is None:
                 raise ValueError("The parent element is not part of this model.")
-            parent = parent.tree_node
+            parent = parent.treenode
 
         if not isinstance(parent, ElementNode):
             raise ValueError("Parent should be an Element or ElementNode of the current model.")
@@ -306,6 +297,8 @@ def add_element(self, element, parent=None, material=None):
         if material:
             self.assign_material(material=material, element=element)
 
+        element.model = self
+
         return element_node
 
     def add_elements(self, elements, parent=None):
@@ -377,8 +370,8 @@ def add_interaction(self, a, b, interaction=None):
         if not self.has_element(a) or not self.has_element(b):
             raise Exception("Please add both elements to the model first.")
 
-        node_a = a.graph_node
-        node_b = b.graph_node
+        node_a = a.graphnode
+        node_b = b.graphnode
 
         if not self.graph.has_node(node_a) or not self.graph.has_node(node_b):
             raise Exception("Something went wrong: the elements are not in the interaction graph.")
@@ -411,8 +404,8 @@ def remove_element(self, element):
             raise Exception("Element not in the model.")
         del self._guid_element[guid]
 
-        self.graph.delete_node(element.graph_node)
-        self.tree.remove(element.tree_node)
+        self.graph.delete_node(element.graphnode)
+        self.tree.remove(element.treenode)
 
     def remove_interaction(self, a, b, interaction=None):
         # type: (Element, Element, Interaction) -> None
@@ -432,12 +425,12 @@ def remove_interaction(self, a, b, interaction=None):
         if interaction:
             raise NotImplementedError
 
-        edge = a.graph_node, b.graph_node
+        edge = a.graphnode, b.graphnode
         if self.graph.has_edge(edge):
             self.graph.delete_edge(edge)
             return
 
-        edge = b.graph_node, a.graph_node
+        edge = b.graphnode, a.graphnode
         if self.graph.has_edge(edge):
             self.graph.delete_edge(edge)
             return
diff --git a/src/compas_model/notebook/scene/__init__.py b/src/compas_model/notebook/scene/__init__.py
index 3893b126..e219c18d 100644
--- a/src/compas_model/notebook/scene/__init__.py
+++ b/src/compas_model/notebook/scene/__init__.py
@@ -18,7 +18,8 @@
 def register_scene_objects():
     register(BlockElement, ThreeBlockObject, context="Notebook")
     register(Model, ThreeModelObject, context="Notebook")
-    print("PyThreeJS Model elements registered.")
+
+    # print("PyThreeJS Model elements registered.")
 
 
 __all__ = [
diff --git a/tests/test_element.py b/tests/test_element.py
new file mode 100644
index 00000000..acbf651a
--- /dev/null
+++ b/tests/test_element.py
@@ -0,0 +1,52 @@
+from pytest import fixture
+
+from compas_model.models import Model
+from compas_model.elements import BlockElement
+from compas_model.interactions import Interaction
+from compas.datastructures import Mesh
+
+
+@fixture
+def my_model():
+    model = Model()
+    a = BlockElement(name="a", shape=Mesh.from_polyhedron(6))
+    b = BlockElement(name="b", shape=Mesh.from_polyhedron(6))
+    c = BlockElement(name="c", shape=Mesh.from_polyhedron(6))
+    d = BlockElement(name="d", shape=Mesh.from_polyhedron(6))
+    model.add_element(a)
+    model.add_element(b)
+    model.add_element(c)
+    model.add_element(d)
+    a.is_dirty = False
+    b.is_dirty = False
+    c.is_dirty = False
+    d.is_dirty = False
+    model.add_interaction(a, c, interaction=Interaction(name="i_a_c"))
+    model.add_interaction(a, b, interaction=Interaction(name="i_b_c"))
+    return model
+
+
+def test_is_dirty(my_model):
+    elements = list(my_model.elements())
+
+    print(elements[0].is_dirty)
+    print(elements[1].is_dirty)
+    print(elements[2].is_dirty)
+    print(elements[3].is_dirty)
+
+    assert not elements[0].is_dirty
+    assert not elements[1].is_dirty
+    assert not elements[2].is_dirty
+    assert not elements[3].is_dirty
+
+    elements[0].is_dirty = True
+
+    print(elements[0].is_dirty)
+    print(elements[1].is_dirty)
+    print(elements[2].is_dirty)
+    print(elements[3].is_dirty)
+
+    assert elements[0].is_dirty
+    assert elements[1].is_dirty
+    assert elements[2].is_dirty
+    assert not elements[3].is_dirty