Computing the closest distance between two faces in Mechanical

Pierre Thieffry

This is a request I got from a customer. This could be useful in estimating distance between non parallel faces for contact. While this is available in SpaceClaim with the measurement tool, Mechanical only computes the distance based on centroids of faces.

How can this be achieved in Mechanical?

James Derrick

Hi @Pierre Thieffry This doesn't appear to be a question? Would you be able to edit it into the standard Question/Answer format?

Pierre Thieffry

The below script requires the user to pick two faces. When the script runs, it will compute the distance between the two faces using various projections methods.

The main one computes the distance between vertices of a face and the other face. To do this, I use the PointAtParam and ParamAtPoint functions. Then I check whether the projection belongs to the face (this is not always the case) by using the bounding box of the face.

The script implements a BoundingBoxData class that will retrieve some information like the lengths of the edges of the box, its centroid, and unit vectors defining its orientation.

# Compute distance between two faces
# Created by P. Thieffry
# Last update: 03/17/21

import math
import units
clr.AddReference("Ans.UI.Toolkit.Base")
clr.AddReference("Ans.UI.Toolkit")
from Ansys.UI.Toolkit import *

# Utility functions
def distance(p1,p2):
    return math.sqrt((p2[0]-p1[0])**2+(p2[1]-p1[1])**2+(p2[2]-p1[2])**2)
    
def dotProduct(v1,v2):
    return v2[0]*v1[0]+v2[1]*v1[1]+v2[2]*v1[2]
    
def vecLength(vec):
    return math.sqrt(vec[0]**2+vec[1]**2+vec[2]**2)

class BoundingBoxData:
# Class to compute geometric elements of a bounding box such as elementary vectors, edge length and midpoint
# Assumes ABCDEFGH points where A is origin, G is diagonal point, orientation triad is AB / AC /AE
    def __init__(self,face):
        bbox=face.GetBoundingBox()
        self.Face=face
        self.Origin=[bbox[0],bbox[1],bbox[2]]
        self.Diagonal=[bbox[3],bbox[4],bbox[5]]
        self.Midpoint=[0.5*(self.Diagonal[0]+self.Origin[0]),0.5*(self.Diagonal[1]+self.Origin[1]),0.5*(self.Diagonal[2]+self.Origin[2])]
        pointB=[bbox[3],bbox[1],bbox[2]]
        pointC=[bbox[0],bbox[4],bbox[2]]
        pointE=[bbox[0],bbox[1],bbox[5]]
        
        self.vecX=[pointB[0]-self.Origin[0],pointB[1]-self.Origin[1],pointB[2]-self.Origin[2]]
        self.vecY=[pointC[0]-self.Origin[0],pointC[1]-self.Origin[1],pointC[2]-self.Origin[2]]
        self.vecZ=[pointE[0]-self.Origin[0],pointE[1]-self.Origin[1],pointE[2]-self.Origin[2]]    
        self.xlength=distance(pointB,self.Origin)
        self.ylength=distance(pointC,self.Origin)
        self.zlength=distance(pointE,self.Origin)
        
        for i in range(0,3):
            if self.xlength != 0.: self.vecX[i]=self.vecX[i]/self.xlength
            if self.ylength != 0.: self.vecY[i]=self.vecY[i]/self.ylength
            if self.zlength != 0.: self.vecZ[i]=self.vecZ[i]/self.zlength
        
    def CheckPointInBox(self,point):
    # Check if a point [x,y,z] lies within the bounding box
        isInside=True
        vector=[point[0]-self.Midpoint[0],point[1]-self.Midpoint[1],point[2]-self.Midpoint[2]]
        if self.xlength != 0:
            checkX=abs(dotProduct(vector,self.vecX))-self.xlength/2<=1e-8
        else:
            checkX=abs(point[0]-self.Origin[0])<1e-8
        
        if self.ylength != 0:
            checkY=abs(dotProduct(vector,self.vecY))-self.ylength/2<=1e-8
        else:
            checkY=abs(point[1]-self.Origin[1])<1e-8
            
        if self.zlength != 0:
            checkZ=abs(dotProduct(vector,self.vecZ))-self.zlength/2<=1e-8
        else:
            checkZ=abs(point[2]-self.Origin[2])<1e-8
            
        if not(checkX and checkY and checkZ):
            isInside=False
           
        return isInside
        
def checkEdgeDistance(edge,bbox):
# Compute distance between an edge and a face
# edge and face from bbox are the ones to compare
# bbox is the boundingboxdata object related to face
# start_vertex holds the coordinates of the end vertex of the edge to start from when scanning the edge (in [x,y,z] form)
#
    mindist=1e31
    
    deltaP=0.05 # increment in parametric value on edge
    p0=0    
    
    face=bbox.Face
    
    while p0<=1:
        vert=edge.PointAtParam(p0) # get point at p0 
        pa=face.ParamAtPoint((vert[0],vert[1],vert[2])) # project on face
        pt=face.PointAtParam(pa[0],pa[1]) # from (u,v), compute projection coordinates
        if bbox.CheckPointInBox(pt):   # verify if projection is on the face or not using bounding box            
            dist=distance(vert,pt)
            if dist<mindist:
                mindist=dist
        p0=p0+deltaP
    return mindist
    

minDistSource={'f1vf2': 'between first face vertices and second face',
    'f2vf1': 'between second face vertices and first face',
    'f1vf2e': 'between first face vertices and second face edges',
    'f2vf1e': 'between second face vertices and first face edges',
    'f1vf2mod': 'between first face edges and second face',
    'f2vf1mod': 'between second face edges and first face',
}

curSel=ExtAPI.SelectionManager.CurrentSelection

# Need to have two faces selected, otherwise pass
if (curSel.Ids.Count!=2):
    mess_line1 = 'Please select 2 faces'
    MessageBox.Show(mess_line1)
    pass

# Check if selection contains faces, otherwise don't do anything
f1=ExtAPI.DataModel.GeoData.GeoEntityById(curSel.Ids[0])
if (f1.GetType() == Ansys.ACT.Common.Geometry.GeoFaceWrapper):
    f2=ExtAPI.DataModel.GeoData.GeoEntityById(curSel.Ids[1])
    
    # Compute bounding box data once to be used in later checks
    bbox1=BoundingBoxData(f1)
    bbox2=BoundingBoxData(f2)
    
    edgef1=f1.Edges
    edgef2=f2.Edges
    
    edgef1Id=[]
    for edge in edgef1:
        edgef1Id.append(edge.Id)
        
    edgef2Id=[]
    for edge in edgef2:
        edgef2Id.append(edge.Id)
    
    # Collect vertices from the two faces - if no vertices, retrieve points from edges
    vertf1=[]
    if len(f1.Vertices)!=0:
        for vert in f1.Vertices:
            vertf1.append([vert.X,vert.Y,vert.Z])
    else:
        for edge in edgef1:
            for i in range(0,len(edge.Points)/3):
                vertf1.append([edge.Points[3*i],edge.Points[3*i+1],edge.Points[3*i+2]])
    
    vertf2=[]
    if len(f2.Vertices)!=0:
        for vert in f2.Vertices:
            vertf2.append([vert.X,vert.Y,vert.Z])
    else:
        for edge in edgef2:
            for i in range(0,len(edge.Points)/3):
                vertf2.append([edge.Points[3*i],edge.Points[3*i+1],edge.Points[3*i+2]])
    
    # Compute shortest distance
    mindist=1e31
    minSource=''
    
    # Run 4 trials:
    # - f1 vertices projected on f2
    # - f2 vertices projected on f1
    # - f1 vertices projected on edges of f2
    # - f2 vertices projected on edges of f1
    
    for vert in vertf1: # Compare vertices in f1 vs projection on f2        
        pa=f2.ParamAtPoint((vert[0],vert[1],vert[2])) # project vert on f2 
        pt=f2.PointAtParam(pa[0],pa[1]) # from (u,v), compute projection coordinates
        if bbox2.CheckPointInBox(pt):   # verify if projection is on the face or not using bounding box     
            dist=distance(vert,pt)
            if dist < mindist:
                mindist=dist
                minVertex=vert
                minSource='f1vf2'
                                
    if mindist!=1e31: # a minimum has been found, scan edges related to closest vertex to see if a closer location exist
        for vert in f1.Vertices:
            if vert.X==minVertex[0] and vert.Y==minVertex[1] and vert.Z==minVertex[2]:
                for edge in vert.Edges: # scan edges related to closest vertex 
                    if edge.Id in edgef1Id: # scan only edges related to face 1
                        mindist2=checkEdgeDistance(edge,bbox2)
                        if mindist2 < mindist:
                            mindist=mindist2
                            minSource='f1vf2mod'    
    
    for vert in vertf2: # Compare vertices in f2 vs projection on f1 - same as above
        pa=f1.ParamAtPoint((vert[0],vert[1],vert[2])) # from (u,v), compute projection coordinates
        pt=f1.PointAtParam(pa[0],pa[1])
        if bbox1.CheckPointInBox(pt):
            dist=distance(vert,pt)
            if dist < mindist:
                mindist=dist
                minVertex=vert
                minSource='f2vf1'
               
    if minSource=='f2vf1': # a minimum has been found, scan edges related to closest vertex to see if a closer location exist
        for vert in f2.Vertices:
            if vert.X==minVertex[0] and vert.Y==minVertex[1] and vert.Z==minVertex[2]:
                for edge in vert.Edges:
                    if edge.Id in edgef2Id:
                        mindist2=checkEdgeDistance(edge,bbox1)
                        if mindist2 < mindist:
                            mindist=mindist2
                            minSource='f2vf1mod'
                
    for vert in vertf1: # Compare vertices in f1 vs projection on edges of f2
        for edge in edgef2:
            pa=edge.ParamAtPoint((vert[0],vert[1],vert[2])) 
            if pa>=0 and pa<=1: # verify if projection is on the edge using reduced coordinate
                pt=edge.PointAtParam(pa)
                dist=distance(vert,pt)
                if dist < mindist:
                    mindist=dist
                    minSource='f1vf2e'
                        
    for vert in vertf2: # Compare vertices in f2 vs projection on edges of f1
        for edge in edgef1:
            pa=edge.ParamAtPoint((vert[0],vert[1],vert[2]))
            if pa>=0 and pa<=1:
                pt=edge.PointAtParam(pa)
                dist=distance(vert,pt)
                if dist < mindist:
                    mindist=dist                
                    minSource='f2vf1e'
                        
    length_unit=ExtAPI.DataModel.GeoData.Unit
    curUnit=ExtAPI.DataModel.CurrentUnitFromQuantityName("Length")
    
    # Convert to Mechanical units
    mindist_curUnit=units.ConvertUnit(mindist, fromUnit=length_unit, toUnit=curUnit)