How to find the rigid transformation data for External Data?

Pierre Thieffry

External Data allows to enter a rigid transformation in the global coordinate system in order to orient data that could be misaligned with the geometry. Finding the proper angles and translation data is not obvious when the data is not oriented along the axis of the global CS and could require a long "trial and error" process.

How can we automate this within Mechanical? The assumption is that we have the "geometry" of the data to be imported either in geometric form or as STL or as a mesh.

Pierre Thieffry

Here's a proposal. For comprehensiveness, it is assumed that the user has been able to get a representation of the data to be imported as geometry in Mechanical (either from a CAD/STL file or a mesh). The geometry is used to visually check the alignment. Yet it is not mandatory, as the script will only compute the transformation if no body is selected.

The principle is the following: the user will define two coordinate systems to be aligned. One is created on the target geometry (the one on which data will be mapped) and another one on the data to be imported. The user then selects in the tree both coordinate systems and the body to be moved (the one correspionding the data to be imported). A Part Transform will be created to show how the alignment has been achieved. If no body is selected, a message box with the rigid transformation between the two coordinate systems is displayed.

# Compute Rigid Transformation between two coordinate systems
# Optionally create a Part transform on a body
# Created by P. Thieffry
# Last update: 2021/04/07


import math
import units

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


# Utility functions



def dotProduct(v1,v2):
    return v2[0]*v1[0]+v2[1]*v1[1]+v2[2]*v1[2]
    
def matProduct(m1,m2):
    nr1=len(m1)
    nc1=len(m1[0])
    nc2=len(m2[0])
    resM = [([0]*nc2) for i in range(nr1)]
    for i in range(0,nr1):
        for j in range(0,nc2):
            val=0.
            for k in range(0,nc1):
                val+=resM[i][j]+m1[i][k]*m2[k][j]
            resM[i][j]=val
    return resM
    
def InverseTransformationMatrix(m):
    mt=[[0.,0.,0.,0.],[0.,0.,0.,0.],[0.,0.,0.,0.],[0.,0.,0.,0.]]
    for i in range(0,4):
        for j in range(0,4):
            mt[i][j]=m[j][i]
            
    tVec=matProduct(mt,[[-m[0][3]],[-m[1][3]],[-m[2][3]],[-m[3][3]]])
    
    for i in range(0,4):
        mt[i][3]=tVec[i][0]
        mt[3][i]=0.
    mt[3][3]=1.
    return mt


def CSBodyDistance(csys,body):
    op_unit=csys.OriginX.Unit
    orX=units.ConvertUnit(csys.OriginX.Value, fromUnit=op_unit, toUnit='m')
    orY=units.ConvertUnit(csys.OriginY.Value, fromUnit=op_unit, toUnit='m')
    orZ=units.ConvertUnit(csys.OriginZ.Value, fromUnit=op_unit, toUnit='m')
    dist=(body.Centroid[0]-orX)**2+(body.Centroid[1]-orY)**2+(body.Centroid[2]-orZ)**2
    return dist

def CreateSelection(body):
    # Create selection from body
    ExtAPI.SelectionManager.ClearSelection()        
    temp_sel= ExtAPI.SelectionManager.CreateSelectionInfo(SelectionTypeEnum.GeometryEntities)                
    temp_sel.Ids=[body.Id]
    ExtAPI.SelectionManager.NewSelection(temp_sel)
    
def CreateCSOnBody(body):
    # Create coordinate system on body so axes are along the axis of inertia of the body
    CreateSelection(body)
    csGroup=ExtAPI.DataModel.GetObjectsByType(DataModelObjectCategory.CoordinateSystems)[0]
    newCS=csGroup.AddCoordinateSystem()
    newCS.Name='CS_body_'+body.Name
#    newCS.OriginDefineBy=Ansys.Mechanical.DataModel.Enums.CoordinateSystemAlignmentType.Fixed
    return newCS
    


def anglesFromTransMatrix(tmat):
    pi=math.pi
    sY=-tmat[2][0]
    
    if sY < 1.:
        if (sY > -1.):
            rotY=math.asin(sY)
            rotX = math.atan2(tmat[2][1],tmat[2][2])
            rotZ = math.atan2(tmat[1][0],tmat[0][0])
        else: # sY = -1
            rotY= -pi/2
            rotX = -math.atan2(tmat[1][0],tmat[1][1])
            rotZ = 0.
    else: # sY=1
        rotY= pi/2
        rotX = -math.atan2(tmat[0][1],tmat[1][1])
        rotZ = 0.
    
    return rotX*180./pi,rotY*180./pi,rotZ*180./pi
    
def TransformationMatrix(origCS,targCS):
    # compute transformation matrix to align origCS to targCS
    # All coordinates are expressed in the global CS
    #
    # Retrieve CS origins in m
    op_unit=origCS.OriginX.Unit
    orX=units.ConvertUnit(origCS.OriginX.Value, fromUnit=op_unit, toUnit='m')
    orY=units.ConvertUnit(origCS.OriginY.Value, fromUnit=op_unit, toUnit='m')
    orZ=units.ConvertUnit(origCS.OriginZ.Value, fromUnit=op_unit, toUnit='m')
    taX=units.ConvertUnit(targCS.OriginX.Value, fromUnit=op_unit, toUnit='m')
    taY=units.ConvertUnit(targCS.OriginY.Value, fromUnit=op_unit, toUnit='m')
    taZ=units.ConvertUnit(targCS.OriginZ.Value, fromUnit=op_unit, toUnit='m')
    
    tMat=[[dotProduct(origCS.XAxis,targCS.XAxis),dotProduct(origCS.XAxis,targCS.YAxis),dotProduct(origCS.XAxis,targCS.ZAxis),-orX+taX],
    [dotProduct(origCS.YAxis,targCS.XAxis),dotProduct(origCS.YAxis,targCS.YAxis),dotProduct(origCS.YAxis,targCS.ZAxis),-orY+taY],
    [dotProduct(origCS.ZAxis,targCS.XAxis),dotProduct(origCS.ZAxis,targCS.YAxis),dotProduct(origCS.ZAxis,targCS.ZAxis),-orZ+taZ],
    [0.,0.,0.,1.]]
   
    
    return tMat
    
def TransformationData(origCS,targCS):
    # compute angles and translations to be used in part transforms
    # based on transformation matrix between two CS
    #
    
    csGroup=ExtAPI.DataModel.GetObjectsByType(DataModelObjectCategory.CoordinateSystems)[0]
    globalCS=csGroup.Children[0] # get global CS as this is the one we want angles and translations expressed in
    
    # First transform from globalCS to origCS
    mGlobalToOrig=TransformationMatrix(globalCS,origCS)

    # Compute inverse transformation
    mOrigToGlobal=InverseTransformationMatrix(mGlobalToOrig)
    # Transfom from global to targ
    mGlobalToTarg=TransformationMatrix(globalCS,targCS)
    # Now combine them
    mOrigToTarg=matProduct(mGlobalToTarg,mOrigToGlobal)
    
    # Now retrieve angle for transformation in global CS
    rotX,rotY,rotZ=anglesFromTransMatrix(mOrigToTarg)
    # Get translations in global CS
    trX=mOrigToTarg[0][3]
    trY=mOrigToTarg[1][3]
    trZ=mOrigToTarg[2][3]
    

    return trX,trY,trZ,rotX,rotY,rotZ

    
    
def CreateTransform(body,bodyCS,targCS):    
    # Transform body from bodyCS to targCS 
    # bodyCS and targCS should be oriented along axes of inertia
    # Angles and translations to be computed based on global Coordinate System so they can be reused in External Data
    #
    csGroup=ExtAPI.DataModel.GetObjectsByType(DataModelObjectCategory.CoordinateSystems)[0]
    globalCS=csGroup.Children[0] # get global CS as this is the one we want angles and translations expressed in
    
    # get translations and angles
    tX,tY,tZ,rotX,rotY,rotZ=TransformationData(bodyCS,targCS)
    
    # create part transform
    # first check if Transforms group exist, if not create it with a dummy part transform
    trGroup=ExtAPI.DataModel.GetObjectsByType(DataModelObjectCategory.PartTransformGroup)
    if len(trGroup)==0:
        dummyT=ExtAPI.DataModel.Project.Model.AddPartTransform()
        dummyT.Delete()
    trGroup=ExtAPI.DataModel.GetObjectsByType(DataModelObjectCategory.PartTransformGroup)[0]
    # create part transform
    newTR=trGroup.AddPartTransform()
    # Set geometry scoping on body
    CreateSelection(body)
    newTR.Location=ExtAPI.SelectionManager.CurrentSelection
    
    # All data expressed in globalCS for reuse in External Data
    newTR.CoordinateSystem=globalCS
    
    # Set transformation data
    newTR.TranslationX=Quantity(tX,'m')
    newTR.TranslationY=Quantity(tY,'m')
    newTR.TranslationZ=Quantity(tZ,'m')

    newTR.RotationX=Quantity(rotX,'deg')
    newTR.RotationY=Quantity(rotY,'deg')
    newTR.RotationZ=Quantity(rotZ,'deg')
        
    # transform gemetry
    newTR.TransformGeometry()
    
    return newTR

curSel=ExtAPI.DataModel.Tree.ActiveObjects

# Need to have two CS selected, otherwise stop
mess_line1 = 'Please select at least 2 coordinate systems in the tree and a maximum of one body if you want to visually check the transformation'
if (curSel.Count<2):
    MessageBox.Show(mess_line1)
    pass
else:        
    bodies=[]
    cSyst=[]
    for sel in curSel:
        if sel.GetType()==Ansys.ACT.Automation.Mechanical.Body:
            bodies.append(sel.GetGeoBody())
        if sel.GetType()==Ansys.ACT.Automation.Mechanical.CoordinateSystem:
            cSyst.append(sel)
            
    if cSyst.Count!=2:
        MessageBox.Show(mess_line1)
        pass
    
    if cSyst.Count==2 and bodies.Count>1:
        MessageBox.Show(mess_line1)
        pass

    if bodies.Count == 1: # a body has been selected, we'll create a part transform
        #find CS closest to selected body
        
        d1=CSBodyDistance(cSyst[0],bodies[0])
        d2=CSBodyDistance(cSyst[1],bodies[0])
        
        if d1<d2:
            origCS=cSyst[0]
            targCS=cSyst[1]
        
        else:
            origCS=cSyst[1]
            targCS=cSyst[0]
                
        CreateTransform(bodies[0],origCS,targCS)
    else: # no body has been selected, just compute transformation data and display in a message box
        origCS=cSyst[0] # first selected CS is the one corresponding to external data
        targCS=cSyst[1] # second CS is the one corresponding to actual geometry to be mapped on
        tX,tY,tZ,rotX,rotY,rotZ=TransformationData(origCS,targCS)
        op_Unit=origCS.OriginX.Unit
        
        mess_line = 'Computed data for rigid transformation:\n\n'
        
        mess_line+='  Orign X : '+str(Quantity(units.ConvertUnit(tX, fromUnit='m', toUnit=op_Unit),op_Unit))+'\n'
        mess_line+='  Orign Y : '+str(Quantity(units.ConvertUnit(tY, fromUnit='m', toUnit=op_Unit),op_Unit))+'\n'
        mess_line+='  Orign Z : '+str(Quantity(units.ConvertUnit(tZ, fromUnit='m', toUnit=op_Unit),op_Unit))+'\n'
        mess_line+='\n'
        mess_line+='  Rotation X : '+ str(Quantity(rotX,'deg'))+'\n'
        mess_line+='  Rotation Y : '+ str(Quantity(rotY,'deg'))+'\n'
        mess_line+='  Rotation Z : '+ str(Quantity(rotZ,'deg'))+'\n'        
        
        MessageBox.Show(mess_line)