- Goal
The following script enables you to calculate the shortest distance between two faces.
- Code
# Compute distance between two faces 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 mess_line1 = 'Please select 2 faces' if (curSel.Ids.Count!=2): MessageBox.Show(mess_line1) pass face=ExtAPI.DataModel.GeoData.GeoEntityById(curSel.Ids[0]) if face.Type != GeoCellTypeEnum.GeoFace: MessageBox.Show(mess_line1) pass # Check if selection contains faces, otherwise don't do anything f1=ExtAPI.DataModel.GeoData.GeoEntityById(curSel.Ids[0]) 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) # Display result mess_line = 'Computed distance (approx) = ' + str(round(mindist_curUnit,3)) + ' ' + str(curUnit)+'\n' #mess_line += '(found '+minDistSource[minSource]+')\n' MessageBox.Show(mess_line)