# Computing the closest distance between two faces in Mechanical

Options
Member, Moderator, Employee Posts: 101
edited June 2023

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?

Tagged:

Options

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

• Member, Moderator, Employee Posts: 101
Options

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
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)
```