Goal: Given a named selection of N faces and another of N*2 faces, add a joint between each face in the first named selection and the nearest two faces in the second named selection.

Code:

import math

#distance formula where c1 and c2 are each arrays of 3 real numbers
def dist(c1, c2):
    x = c1[0]-c2[0]
    x = x*x
    y = c1[1]-c2[1]
    y = y*y
    z = c1[2]-c2[2]
    z = z*z
    return math.sqrt(x+y+z)

jh1 = DataModel.GetObjectsByName("jh1")[0]       #named selection object of faces to become the single joint reference scope
jh2 = DataModel.GetObjectsByName("jh2")[0]       #named selection object of faces from which the two closest will become the mobile scope
group = DataModel.GetObjectsByName("Connection Group")[0]       #connection group to hold all of the joints that are added

geo = DataModel.GeoData
face_ids1 = jh1.Location.Ids       #face ids in jh1
face_ids2 = jh2.Location.Ids       #face ids in jh2

#get the centroids of all faces in face_ids2
face_centroids2 = [geo.GeoEntityById(face_id2).Centroid for face_id2 in face_ids2]

def get_min_indeces(face_centroid1):
    #array to hold the index and distance to face1_center for two closest faces in face_centroids2
    min_indeces = [None, None]

    for index in range(len(face_centroids2)):
        face_centroid2 = face_centroids2[index]
        distance = dist(face_centroid1, face_centroid2)

        #initilize with the first two iterations of this loop
        if index == 0:
            min_indeces[0] = (index, distance)
            continue
        if index == 1:
            min_indeces[1] = (index, distance)
            continue

        #get the index of the item with the larger distance
        if min_indeces[0][1] < min_indeces[1][1]:
            larger_dist_index = 1
        else:
            larger_dist_index = 0

        #replace that item with the current face if the distance is smaller
        if distance < min_indeces[larger_dist_index][1]:
            min_indeces[larger_dist_index] = (index, distance)
    return min_indeces

#use a transaction to speed up adding joints to the tree in a loop
with Transaction():
    for face_id1 in face_ids1:
        face1 = geo.GeoEntityById(face_id1)
        face_centroid1 = face1.Centroid       #get the centroid of each face in face_ids1
        min_indeces = get_min_indeces(face_centroid1)

        #get the face ids using the indices computed above
        face_ids = [face_ids2[min_indeces[0][0]], face_ids2[min_indeces[1][0]]]

        #add a joint and assign its geometry selection
        joint = group.AddJoint()

        reference_selection = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.GeometryEntities)
        reference_selection.Ids = [face_id1]

        mobile_selection = ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.GeometryEntities)
        mobile_selection.Ids = face_ids

        joint.ReferenceLocation = reference_selection
        joint.MobileLocation = mobile_selection