Skip to main content

Aligning external imported data with your Mechanical Model via Python scripting

jimmy.he@ansys.com | 11.20.2024

When defining our Mechanical model, sometimes we need to apply loads (like surface pressure or forces) or other input data (such as material distribution or orientation) that comes from previous simulation results which are imported via the External Data function in Workbench. Since the previous simulation may not be conducted in the same coordinate system or even geometry as our current Mechanical analysis, an alignment process is often needed to align the externally imported data to the current Mechanical model.

Mechanical provides two methods for aligning the data under the Rigid Transformation detail:

  1. Use Origin and Euler Angles: Define the translations (in X,Y and Z directions) and rotations (in terms of Euler angles) for the transformation
  2. Use Coordinate Systems: Define the source (attached to imported data) and destination (attached to current model) coordinate systems

Regardless of which option to choose, it might not be immediately obvious to the user what numbers to put down or how to manually define the coordiante systems. This can easily become an lengthy trial-and-error process to try to eye-ball the motion needed to align the two models. In this blog, we demonstrate how we can automate this process with Python scripting within Mechanical.

General idea of the workflow

The overall idea of this python script is simple: create a source coordinate system attached to the imported data, and another one attached to the current model, then use the second option as mentioned above for alignment. We will be creating a coordinate system using three points that are not co-linear.

Coordinate system creation

First, let's define some helper functions for vector math calculations. These functions will be used to compute the information necessary to create a coordinate system from three points.

def norm(v):
    return (v[0]**2 + v[1]**2 + v[2]**2)**0.5

def unit_vec(v):
    n = norm(v)
    return [ i/n for i in v ]

def cross(a, b):
    c = [a[1]*b[2] - a[2]*b[1],
         a[2]*b[0] - a[0]*b[2],
         a[0]*b[1] - a[1]*b[0]]
    return c

def dotProduct(v1,v2):
    return v2[0]*v1[0]+v2[1]*v1[1]+v2[2]*v1[2]
    
def subtract(a , b):
    return [a[i] - b[i] for i in range(3)]

def add(a , b):
    return [a[i] + b[i] for i in range(3)]

Next, we define the function that forms a set of orthogonal coordinate basis from three random points (that are not co-linear) via cross product:

def pt2CS( pts ):
    X = unit_vec(subtract(pts[1] , pts[0]))  # X
    Y_ = subtract(pts[2] , pts[0])  # Y
    Z = unit_vec(cross(X , Y_ )) #Z
    Y = cross(Z , X)
    return [X , Y , Z]

We will need a function that calculates the 3x3 transformation matrix and convert that the Euler angles in the Z->Y->X rotation order:

def TransformationMatrix(origCS,targCS):
    # compute transformation matrix to align origCS to targCS
    # All coordinates are expressed in the global CS
    #
    tMat=[[dotProduct(origCS[0],targCS[0]),dotProduct(origCS[0],targCS[1]),dotProduct(origCS[0],targCS[2])],
    [dotProduct(origCS[1],targCS[0]),dotProduct(origCS[1],targCS[1]),dotProduct(origCS[1],targCS[2])],
    [dotProduct(origCS[2],targCS[0]),dotProduct(origCS[2],targCS[1]),dotProduct(origCS[2],targCS[2])]]
    return tMat

def rotationMatrixToEulerAngles_ZYX(R) :
    sy = math.sqrt(R[0][0] * R[0][0] +  R[1][0] * R[1][0])
    singular = sy < 1e-6
    factor = 180./math.pi
    if  not singular :
        x = math.atan2(R[2][1] , R[2][2])*factor
        y = math.atan2(-R[2][0], sy)*factor
        z = math.atan2(R[1][0], R[0][0])*factor
    else :
        x = math.atan2(-R[1][2], R[1][1])*factor
        y = math.atan2(-R[2][0], sy)*factor
        z = 0
    return [z, y, x]

Then we can define the function that uses these helpers to create a coordinate system in Mechanical from three points:

def CreateCS(origin , xyz , name , unit='mm'):
    # Create coordinate
    testCS = ExtAPI.DataModel.Project.Model.CoordinateSystems.AddCoordinateSystem()
    testCS.OriginDefineBy= Ansys.Mechanical.DataModel.Enums.CoordinateSystemAlignmentType.Fixed
    testCS.OriginX = Quantity(origin[0],unit)
    testCS.OriginY = Quantity(origin[1],unit)
    testCS.OriginZ = Quantity(origin[2],unit)
    
    # Get rotations
    origCS = [[1,0,0],[0,1,0],[0,0,1]]
    tMat = TransformationMatrix(origCS,xyz)
    angs = rotationMatrixToEulerAngles_ZYX(tMat)
    
    # Perform transformation
    testCS.PrimaryAxisDefineBy = Ansys.Mechanical.DataModel.Enums.CoordinateSystemAlignmentType.GlobalX
    # Rotation Z
    testCS.RotateZ()
    testCS.SetTransformationValue(1,angs[0])
    # Rotation Y
    testCS.RotateY()
    testCS.SetTransformationValue(2,angs[1])
    # Rotation X
    testCS.RotateX()
    testCS.SetTransformationValue(3,angs[2])

    # Set name
    testCS.Name = name

Locating the points for coordinate system creation

Now that we have a way to define a coordinate system from three points, we need to look into how to find the coordinates of these three points. For the geometry in the Mechanical model, this is simple, since we can select vertices on the body through the GUI window. Once a vertex is selected, we can use the following script to get its object ID:

print(ExtAPI.SelectionManager.CurrentSelection.Ids[0])

Once we have the object ID, we can get its coordinates in the subsequent code, for example:

vertex_id = 1
curr_vertex = DataModel.GeoData.GeoEntityById(vertex_id)
vertex_coordinate = [curr_vertex.X , curr_vertex.Y , curr_vertex.Z]

The harder part is getting the coordinates for the imported data. Unfortunately, currently there is no way to extract those coordinates using Mechanical APIs for the imported object. This means that we will need to parse the source file ourselves to obtain the coordinates. To achieve this, we will need information on the file path and how to interpret each column in the file. This can be done by the following helper functions:

def GetPath(curr_obj):
    import wbjn

    # Get coordinate of the selected object in the tree
    coord = curr_obj.Parent.Name.split('(')[-1].split(')')[0]
    # Get the source path from Workbench
    WBcmds = '''coord = "<$coord$>"\nimport clr\nclr.AddReference("Ans.UI")\nimport Ansys.UI\nclr.AddReference("Ans.ProjectSchematic")\nimport Ansys.ProjectSchematic\n'''.replace('<$coord$>',coord)
    WBcmds += '''view1 = Ansys.UI.UIManager.Instance.GetActiveWorkspace().GetView(Ansys.ProjectSchematic.View.ProjectSchematicView.ViewName)\ncoord_map = dict(view1.CoordinateMap)\n'''
    WBcmds += '''coord_map2 = {Ansys.UI.IDManager.GetAlphabeticLabelFromCoordinate(coord_map[key]):key for key in coord_map}\nsetup1 = coord_map2[coord]\n'''
    WBcmds += '''try: setup1 = setup1.GetContainer()\nexcept: pass\nexternalLoadData1 = setup1.GetExternalLoadData()\n'''
    WBcmds += '''externalLoadFileData1 = externalLoadData1.GetExternalLoadFileData(Name="ExternalLoadFileData")\nfile = externalLoadFileData1.File.Location\nreturnValue(file)\n'''
    source_path = wbjn.ExecuteCommand(ExtAPI,WBcmds)

    # Get format of each column
    WBcmds2 = '''externalLoadFileDataProperty1=externalLoadFileData1.GetDataProperty()\nxyz = [0,0,0,externalLoadFileDataProperty1.DelimiterType,0]\nfor idx , cd in enumerate( externalLoadFileDataProperty1.ColumnsData ):\n    if 'X Coordinate' in cd.DataType: xyz[0] = idx\n    elif 'Y Coordinate' in cd.DataType: xyz[1] = idx\n    elif 'Z Coordinate' in cd.DataType:\n        xyz[2] = idx\n        xyz[4] = cd.Unit\nreturnValue(str(xyz))\n'''
    tmp = wbjn.ExecuteCommand(ExtAPI,WBcmds2) # This can only be returned as a string

    return source_path , tmp

def GetExtData(curr_obj):
    # Get file path and column formatting
    source_path , tmp = GetPath(curr_obj)

    # Decode the column formatting
    source_format = tmp[1:-1].split(',')
    for i in range(5):
        if i < 3:
            source_format[i] = int(source_format[i])
        else:
            source_format[i] = source_format[i].replace('[','').replace(']','').replace("'",'')
    if 'Tab' in source_format[3]:
        delimit = '\t'
    elif 'Comma' in source_format[3]:
        delimit = ','
    else:
        delimit = source_format[3]
    unit = source_format[-1]

    # Read the source file
    f = open(source_path , 'r')
    L = f.readlines()
    f.close()
    sXYZ = [  ]
    for l in L:
        try:
            tmp = l.split(delimit)
            cc = []
            for i in range(2):
                value = float(tmp[source_format[i]])
                cc.append(value)
            cc.append(0.)
            sXYZ.append(cc)
        except:
            pass
    print('Read ' + str(len(sXYZ)) + ' points from source file')

    return sXYZ , unit

Creating two coordinate systems and align external data

Now that we have all the helper functions defined, we can layout our workflow. We first start by locating the object in the tree that uses an external imported data, and we does so via its name (although other means are possible). Note that we also display the source points and their IDs in the GUI window to allow the user to identify the best-match source point.

curr_obj = DataModel.GetObjectsByName('Imported Pressure')[0] # Here we use an imported pressure as example

# Turn on graphics display and show the IDs of the source points
curr_obj.Activate()
curr_obj.DisplaySourcePoints=True
curr_obj.DisplaySourcePointIds=True

# Read source points and unit
sXYZ , source_unit = GetExtData(curr_obj)

Then, the user chooses three pairs of best-matching points. For each pair, first select a vertex in the geometry (and get its object ID), then locate its best-matching counterpart from the graphics display of the source points, and enter its point ID as displayed on screen:

# Pair 1
v1_id = 1 # Enter the object ID of the selected vertex
s1_id = 2 # Enter the point ID of the best matching source point as showed on screen

# Pair 2
v2_id = 3
s2_id = 4

# Pair 3
v3_id = 5
s3_id = 6

Finally, we can create the coordinate systems from the selected points and perform alignment:

obj_name = curr_obj.Name 

# Get selected pts
source_pts , target_pts = [] , []
for vid , sid in zip( [v1_id,v2_id,v3_id] , [s1_id,s2_id,s3_id] ):
    curr_v = DataModel.GeoData.GeoEntityById(vid)
    target_pts.append([curr_v.X , curr_v.Y , curr_v.Z])
    
    pt = sid - 1
    source_pts.append([sXYZ[pt][0] , sXYZ[pt][1] , sXYZ[pt][2]])

# Build and set coordinate systems
target_cs = pt2CS(target_pts)
CreateCS(target_pts[0] , target_cs , obj_name+'_target' , CADunit)

source_cs = pt2CS(source_pts)
CreateCS(ource_pts[0] , source_cs , obj_name+'_source' , source_unit)

# Align data
curr_obj.RigidBodyTransformationType = Ansys.Mechanical.DataModel.Enums.RigidBodyTransformationType.CoordinateSystem
curr_obj.RigidTransformSourceCoordinateSystem = DataModel.GetObjectsByName(obj_name+'_source')[0]
curr_obj.RigidTransformTargetCoordinateSystem = DataModel.GetObjectsByName(obj_name+'_target')[0]

Conclusion

In this blog, we demonstrated how to create coordinate systems from the coordinates of three non-colinear points, and how we can use two coordinate systems to align the external imported data with the current Mechanical model. With this tool, we can eliminate the guessing work from the alignment process and improve the accuracy of data alignment. If you are looking for further resources for PyMechanical here are some recommendations: