assembly_interfaces_xfunc returns only one face per block pair

[yck011522]

Describe the bug
I have a test data of three blocks. See attached image.

I used the following code to compute the interface:

assembly_interfaces_xfunc = XFunc('compas_assembly.datastructures.assembly_interfaces_xfunc', python='C:\ProgramData\Anaconda3\envs\compas_assembly\python.exe')
data = {'assembly': assembly.to_data(),
        'blocks': {str(key): assembly.blocks[key].to_data() for key in assembly.vertices()}}

result = assembly_interfaces_xfunc(data)
assembly.data = result['assembly']

I expected Interface 1-3 (marked in red) to return all the intersecting faces between the two blocks.
Or is this function designed to only return one face pair?

I digged ito the code of compas_assembly/datastructures/interfaces_numpy.py
But I couldn't understand the code logic, so maybe I'll ask here.
I also looked into paper Data management and modelling of complex interfaces in imperfect discrete-element assemblies but there wasn't any mention of the code behaviour

To Reproduce
Steps to reproduce the behavior:

1. Win 10 64 bit, Rhino 6
2. Sample script attached. Run file 02 and then run 03 inside Rhino Python env.
   InterfaceRhinoPython.zip

Expected behavior
I expect all of the intersecting faces to be return.

[tomvanmele]

i will have a closer look tomorrow, but, if i remember correctly, the algorithm indeed only allows for one interface between two blocks (there can be only one edge between two vertices of the assembly network). this, in part, has to do with how RBE works and therefore should not be a constraint for the general assembly data structure. however, it is also not trivial to encode this properly with the network data structure because there can be at most two edges between any two vertices (if they are in opposite directions). to have more edges than two between any two vertices we need a multigraph implementation...

[tomvanmele]

ha and furthermore...

i only just noticed that block 3 is non-convex. that is in any case something you have to handle differently.

1. split up block 3 in convex parts
2. find the interfaces as usual
3. mark the interface(s) between the parts of block 3 as "internal" to the block

[yck011522]

I think as a generic data structure to represent assembly contact points. The (bidirectional) edges structure can stay the same even for multiple interfaces. The computation of the intersection creates an edge to represent the interface, and a list of contact faces are stored within that edge.

On a separate note, If the algorithm is intended to only work for convex objects. Then maybe it is not right for what I do. I have very often 3 or more contact faces per block pair.

For what I do, I only need the normal from those faces. I don't actually need the curve intersections between the faces. I devised an algorithm to do something similar, it is based on RhinoCommons. I haven't thoroughly tested it yet. Maybe if there is enough interest, I can try plot it to COMPAS only implmentation:

"""
This function computes all the interfaces between two mesh objects.
Interfaces are defined as co-planar mesh faces that overlaps with each other.
(Edges touching is also considered as overlaps in this function)
The meshes can be non-convex.

    Inputs:
        mesh_a: One of the mesh to be detected
        mesh_b: Another one of the mesh to be detected
    Output:
        face_id_a: List of the indices of interface faces on mesh a 
        face_id_b: List of the indices of interface faces on mesh b
        block_vector_a: List of the blocking vector (face normal) 
            corrisponding to the interface faces on mesh a
        block_vector_b: List of the blocking vector (face normal) 
            corrisponding to the interface faces on mesh b
            (block_vector_a and block_vector_b have the same number
            of vectors that points to opposite directions)
            
This function is expected to run within Rhino v6 Grasshopper GHPython Object 
"""

__author__ = "yck011522"
__version__ = "2019.07.17"
threshold = tolerance

import Rhino
import rhinoscriptsyntax as rs

# Extract Face normals
mesh_a.FaceNormals.ComputeFaceNormals()
mesh_a.FaceNormals.UnitizeFaceNormals()
mesh_b.FaceNormals.ComputeFaceNormals()
mesh_b.FaceNormals.UnitizeFaceNormals()

output_face_id_a, output_face_id_b = [],[]
output_block_vector_a, output_block_vector_b = [],[]

#Helper function to check if this block vector is already inserted
def add_vector_to_list(vectorlist, vector_to_add):
    for vector in vectorlist:
        if (vector * vector_to_add > tolerance/100):
            return
    vectorlist.append(vector_to_add)
            
for a in range (mesh_a.Faces.Count): # Loop all faces in  mesh a
    for b in range (mesh_b.Faces.Count): # Loop all faces in  mesh b
        
        #Check normal to be opposite
        normal_a = mesh_a.FaceNormals[a]
        normal_b =  mesh_b.FaceNormals[b]
        if((normal_a + normal_b).Length > threshold):
            continue
            
        #Check if the faces are coplanar
        center_a = (mesh_a.Faces.GetFaceCenter(a))
        center_b = (mesh_b.Faces.GetFaceCenter(b))
        vector_a_to_b = Rhino.Geometry.Point3d.Subtract(center_b, center_a)
        coplanar_distance = Rhino.Geometry.Vector3d(normal_a) * vector_a_to_b # Dot Product
        if (abs(coplanar_distance) > threshold):
            continue
            
        #Get the vertices coordinate of the face corners
        pt_a =[0]*4
        pt_b =[0]*4
        success, pt_a[0], pt_a[1], pt_a[2], pt_a[3] = mesh_a.Faces.GetFaceVertices(a)  
        success, pt_b[0], pt_b[1], pt_b[2], pt_b[3] = mesh_b.Faces.GetFaceVertices(b)
        
        #Cast data type from point3f to point3d
        for i in range(4):
            pt_a[i] = Rhino.Geometry.Point3d(pt_a[i])
            pt_b[i] = Rhino.Geometry.Point3d(pt_b[i])
        
        #Draw polyline to represent the mesh face
        curve_a = Rhino.Geometry.Polyline(pt_a).ToPolylineCurve()
        curve_b = Rhino.Geometry.Polyline(pt_b).ToPolylineCurve()
        
        #Check actual intersection between two mesh faces
        _ , testPlane = Rhino.Geometry.Plane.FitPlaneToPoints(pt_a)
        region_containment_result = Rhino.Geometry.Curve.PlanarClosedCurveRelationship(curve_a,curve_b,testPlane,threshold)
        if (region_containment_result == Rhino.Geometry.RegionContainment.Disjoint):
            continue 
            
        #print(region_containment_result)
        #print (str(a) + "-" + str(b))
        
        #If the two mesh faces matches all the tests, all its property to output list
        output_face_id_a.append(a)
        output_face_id_b.append(b)
                
        add_vector_to_list(output_block_vector_a,normal_a)
        add_vector_to_list(output_block_vector_b,normal_b)
        
#Routing Output Data 
face_id_a = output_face_id_a
face_id_b = output_face_id_b

block_vector_a = output_block_vector_a
block_vector_b = output_block_vector_b