Computing the closest distance between two faces in Mechanical

Options
Pierre Thieffry
Pierre Thieffry Member, Moderator, Employee Posts: 101
First Anniversary Ansys Employee Solution Developer Community of Practice Member Photogenic
edited June 2023 in Structures

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:

Answers

  • James Derrick
    James Derrick Administrator, Employee Posts: 248
    First Anniversary Solution Developer Community of Practice Member Ansys Employee 5 Up Votes
    admin
    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?

  • Pierre Thieffry
    Pierre Thieffry Member, Moderator, Employee Posts: 101
    First Anniversary Ansys Employee Solution Developer Community of Practice Member Photogenic
    Answer ✓
    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
    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)