Merge branch 'modelgeometry'

Author: tomvanmele

CHANGELOG.md

@@ -9,10 +9,25 @@ 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()`.
+
 ### 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 +50,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

docs/_static/PLACEHOLDER

@@ -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
 =======

docs/api/compas_model.interactions.rst

@@ -110,6 +110,12 @@ select = ["E", "F", "I"]
 
 [tool.ruff.lint.isort]
 force-single-line = true
+known-first-party = [
+    "compas",
+    "compas_assembly",
+    "compas_model",
+    "compas_viewer",
+]
 
 [tool.ruff.lint.pydocstyle]
 convention = "numpy"

pyproject.toml

@@ -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

scripts/test_blockmodel_arch.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/__init__.py

@@ -1,39 +1,40 @@
 from math import fabs
-from typing import List
 
 from compas.datastructures import Mesh
 from compas.geometry import Frame
 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 model_interfaces(
     model: Model,
     nmax: int = 10,
     tmax: float = 1e-6,
     amin: float = 1e-2,
     nnbrs_dims: int = 3,
-):
+) -> None:
     """Identify the interfaces between the blocks of an assembly.
 
     Parameters
     ----------
-    assembly : compas_assembly.datastructures.Assembly
-        An assembly of discrete blocks.
+    model : :class:`Model`
+        A model containing blocks with mesh geometry.
     nmax : int, optional
         Maximum number of neighbours per block.
     tmax : float, optional
@@ -43,17 +44,21 @@ def blockmodel_interfaces(
 
     Returns
     -------
-    :class:`Assembly`
+    None
+
+    Notes
+    -----
+    Interface planes are computed from the bestfit frame of the combined points of two faces.
 
     """
     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()}
@@ -68,71 +73,96 @@ def blockmodel_interfaces(
             n = index_node[j]
 
             if n == node:
-                # a block has no interfaces with itself
                 continue
 
             if model.graph.has_edge((n, node), directed=False):
-                # the interfaces between these two blocks have already been identified
                 continue
 
             nbr = blocks[j]
 
             interfaces = mesh_mesh_interfaces(block, nbr, tmax, amin)
 
             if interfaces:
+                # this can't be stored under interactions
+                # it should be in an atteibute called "interfaces"
+                # there can also be "collisions", "overlaps", "..."
                 model.graph.add_edge(node, n, interactions=interfaces)
 
-    return model
-
 
 def mesh_mesh_interfaces(
-    a: BlockGeometry,
-    b: BlockGeometry,
+    a: Mesh,
+    b: Mesh,
     tmax: float = 1e-6,
     amin: float = 1e-1,
 ) -> list[ContactInterface]:
     """Compute all face-face contact interfaces between two meshes.
 
     Parameters
     ----------
-    a : :class:`Block`
-    b : :class:`Block`
+    a : :class:`Mesh`
+        The source mesh.
+    b : :class:`Mesh`
+        The target mesh.
     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`]
+    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 bestfit frame does not align with the normal of the base source frame,
+    it will be inverted, such that it corresponds to whatever edge is created from this source to a target.
 
     """
     world = Frame.worldXY()
-    interfaces = []
-    frames = a.frames()
+    interfaces: list[ContactInterface] = []
 
     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 = 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):
                 continue
 
-            intersection = p0.intersection(p1)
+            intersection: ShapelyPolygon = p0.intersection(p1)
             area = intersection.area
 
             if area < amin:
@@ -141,6 +171,10 @@ def mesh_mesh_interfaces(
             coords = [[x, y, 0.0] for x, y, _ in intersection.exterior.coords]
             coords = transform_points(coords, matrix.inverted())[:-1]
 
+            # this is not always an accurate representation of the interface
+            # if the polygon has holes
+            # the interface is incorrect
+
             interface = ContactInterface(
                 size=area,
                 points=coords,
@@ -172,7 +206,7 @@ def merge_coplanar_interfaces(model: Model, tol: float = 1e-6) -> None:
 
     """
     for edge in model.graph.edges():
-        interfaces: List[ContactInterface] = model.graph.edge_attribute(edge, "interfaces")
+        interfaces: list[ContactInterface] = model.graph.edge_attribute(edge, "interfaces")
 
         if interfaces:
             polygons = []

src/compas_model/algorithms/interfaces.py

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

src/compas_model/algorithms/intersections.py

@@ -1,8 +1,27 @@
+import numpy.typing as npt
 from numpy import asarray
 from scipy.spatial import cKDTree
 
 
-def find_nearest_neighbours(cloud, nmax, dims=3):
+def find_nearest_neighbours(cloud: npt.ArrayLike, nmax: int, dims: int = 3) -> list[tuple[list[float], list[int]]]:
+    """Find the nearest neighbours of each point in a cloud among all other points in the cloud.
+
+    Parameters
+    ----------
+    cloud : array-like
+        The cloud of points.
+    nmax : int
+        The maximum number of neighbours per point.
+    dims : int, optional
+        The number of dimensions to include in the search.
+
+    Results
+    -------
+    list[tuple[list[float], list[int]]]
+        For each point, a tuple with the distances to the nearest neighbours,
+        and the indices to the corresponding points.
+
+    """
     cloud = asarray(cloud)[:, :dims]
     tree = cKDTree(cloud)
     nnbrs = [tree.query(root, nmax) for root in cloud]