Modelgeometry

CHANGELOG.md

@@ -9,10 +9,26 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ### Added
 
+* 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`.
+
 ### Changed
 
+* 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 model reference `ElementTree.model` from `ElementTree`.
+* Removed `InterfaceElement` from elements.
+
 
 ## [0.4.5] 2024-12-11
 
@@ -35,6 +51,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 * Removed `compas_model.models.groupnode.GroupNode`.
 
+
 ## [0.4.4] 2024-06-13
 
 ### Added

README.md

@@ -21,7 +21,7 @@ pip install -e ".[dev]"
 ## Documentation
 
 For further "getting started" instructions, a tutorial, examples, and an API reference,
-please check out the online documentation here: [COMPAS DR docs](https://blockresearchgroup.github.io/compas_model)
+please check out the online documentation here: [COMPAS MODEL docs](https://blockresearchgroup.github.io/compas_model)
 
 ## Issue Tracker
 

docs/api/compas_model.interactions.rst

@@ -4,6 +4,16 @@ compas_model.interactions
 
 .. currentmodule:: compas_model.interactions
 
+This module provides classes for defining the type of interaction that exists between two elements.
+The interaction type could determine, for example, how forces are transferred from one element to the other.
+The interaction type could also determine whether an interaction is permanent or temporary;
+for example, for designing construction sequences.
+The different types of interactions will have to be interpreted by the context in which the model is used.
+
+Interactions do not define the geometry of a joint or interface, but rather how the elements are connected.
+In the case of a wood joint, for example, an interaction could define whether the joinery is dry, glued, or mechanical,
+and what the properties of this connection are.
+
 
 Classes
 =======

scripts/test_blockmodel_arch.py

@@ -3,18 +3,19 @@
 import compas
 from compas.colors import Color
 from compas_assembly.geometry import Arch
-from compas_model.algorithms import blockmodel_interfaces
+from compas_viewer import Viewer
+
+from compas_model.algorithms import model_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()
 
@@ -25,13 +26,13 @@
 # Interfaces
 # =============================================================================
 
-blockmodel_interfaces(model, amin=0.01)
+model_interfaces(model, amin=0.01)
 
 # =============================================================================
 # 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

src/compas_model/algorithms/__init__.py

@@ -1,15 +1,19 @@
 from .collisions import is_aabb_aabb_collision
 from .collisions import is_box_box_collision
 from .collisions import is_face_to_face_collision
-from .collisions import get_collision_pairs
 
-from .interfaces import blockmodel_interfaces
+from .collisions import get_collision_pairs  # rename to model_collisions
+from .interfaces import model_interfaces
+from .intersections import model_intersections
+from .overlaps import model_overlaps
 
 
 __all__ = [
     "is_aabb_aabb_collision",
     "is_box_box_collision",
     "is_face_to_face_collision",
     "get_collision_pairs",
-    "blockmodel_interfaces",
+    "model_interfaces",
+    "model_intersections",
+    "model_overlaps",
 ]

src/compas_model/algorithms/interfaces.py

@@ -6,22 +6,32 @@
 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
 
 from .nnbrs import find_nearest_neighbours
 
 
-def blockmodel_interfaces(
+def invert(self):
+    self._yaxis = self._yaxis * -1
+    self._zaxis = self._zaxis * -1
+
+
+Frame.invert = invert
+
+
+def model_interfaces(
     model: Model,
     nmax: int = 10,
     tmax: float = 1e-6,
@@ -49,11 +59,11 @@ 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))
 
-    block_cloud = [block.centroid() for block in blocks]
+    block_cloud = [block.centroid for block in blocks]
     block_nnbrs = find_nearest_neighbours(block_cloud, nmax, dims=nnbrs_dims)
 
     model.graph.edge = {node: {} for node in model.graph.nodes()}
@@ -86,8 +96,8 @@ def blockmodel_interfaces(
 
 
 def mesh_mesh_interfaces(
-    a: BlockGeometry,
-    b: BlockGeometry,
+    a: Mesh,
+    b: Mesh,
     tmax: float = 1e-6,
     amin: float = 1e-1,
 ) -> list[ContactInterface]:
@@ -106,27 +116,50 @@ 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()
 
     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):

src/compas_model/algorithms/intersections.py

@@ -0,0 +1,2 @@
+def model_intersections():
+    pass

src/compas_model/algorithms/overlaps.py

@@ -0,0 +1,137 @@
+from compas.geometry import Brep
+from compas.tolerance import TOL
+from compas_occ.brep import OCCBrepFace as BrepFace
+
+from compas_model.interactions import ContactInterface
+from compas_model.models import Model
+
+from .nnbrs import find_nearest_neighbours
+
+
+def model_overlaps(
+    model: Model,
+    deflection=None,
+    tolerance=1,
+    nmax: int = 10,
+    tmax: float = 1e-6,
+    amin: float = 1e-2,
+    nnbrs_dims: int = 3,
+):
+    """Identify the interfaces between the blocks of an assembly.
+
+    Parameters
+    ----------
+    assembly : compas_assembly.datastructures.Assembly
+        An assembly of discrete blocks.
+    nmax : int, optional
+        Maximum number of neighbours per block.
+    tmax : float, optional
+        Maximum deviation from the perfectly flat interface plane.
+    amin : float, optional
+        Minimum area of a "face-face" interface.
+
+    Returns
+    -------
+    :class:`Assembly`
+
+    """
+    deflection = deflection or TOL.lineardeflection
+
+    node_index = {node: index for index, node in enumerate(model.graph.nodes())}
+    index_node = {index: node for index, node in enumerate(model.graph.nodes())}
+
+    geometries: list[Brep] = [model.graph.node_element(node).modelgeometry for node in model.graph.nodes()]
+
+    nmax = min(nmax, len(geometries))
+    cloud = [geometry.centroid for geometry in geometries]
+    nnbrs = find_nearest_neighbours(cloud, nmax, dims=nnbrs_dims)
+
+    model.graph.edge = {node: {} for node in model.graph.nodes()}
+
+    for u in model.graph.nodes():
+        i = node_index[u]
+        A = geometries[i]
+
+        nbrs = nnbrs[i][1]
+
+        for j in nbrs:
+            v = index_node[j]
+
+            if u == v:
+                continue
+
+            if model.graph.has_edge((u, v), directed=False):
+                continue
+
+            B = geometries[j]
+
+            overlaps = brep_brep_overlaps(A, B, deflection=deflection, tolerance=tolerance, tmax=tmax, amin=amin)
+
+            if overlaps:
+                model.graph.add_edge(u, v, interactions=overlaps)
+
+    return model
+
+
+def brep_brep_overlaps(
+    A: Brep,
+    B: Brep,
+    deflection=None,
+    tolerance=1,
+    tmax: float = 1e-6,
+    amin: float = 1e-2,
+) -> list[ContactInterface]:
+    """Compute all face-face contact interfaces between two meshes.
+
+    Parameters
+    ----------
+    a : :class:`Block`
+    b : :class:`Block`
+    tmax : float, optional
+        Maximum deviation from the perfectly flat interface plane.
+    amin : float, optional
+        Minimum area of a "face-face" interface.
+
+    Returns
+    -------
+    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
+
+    """
+    faces_A, faces_B = A.overlap(B, deflection=deflection, tolerance=tolerance)
+    faces_A: list[BrepFace]
+    faces_B: list[BrepFace]
+
+    overlaps: list[ContactInterface] = []
+
+    if faces_A and faces_B:
+        for face_A in faces_A:
+            brep_A = Brep.from_brepfaces([face_A])
+
+            if brep_A.area < amin:
+                continue
+
+            for face_B in faces_B:
+                brep_B = Brep.from_brepfaces([face_B])
+
+                if brep_B.area < amin:
+                    continue
+
+                brep_C: Brep = Brep.from_boolean_intersection(brep_A, brep_B)
+
+                if brep_C.area < amin:
+                    continue
+
+                poly_C = brep_C.to_polygons()[0]
+                mesh_C = brep_C.to_tesselation()[0]
+
+                overlap = ContactInterface(points=poly_C.points, mesh=mesh_C)
+                overlaps.append(overlap)
+
+    return overlaps