Plot contact pressure along deformed edge

Pernelle Marone-Hitz

I have a model as follows:

I'd like to plot the pressure variation for nodes in contact with the upper plate, along with the deformation profile.

Is this possible?

Pernelle Marone-Hitz

Credits should go to @Pierre Thieffry who created the example.
This can be achieved thanks to a Python Code object in Mechanical.

The logic behind is as follows:

• Create a named selection to get the nodes of the profile: - Note - when passed to the solver, the name of the named selection gets converted to capital letters.
• Use DPF to retrieve:
• *the contact pressure on the contact elements. This is done by selecting contact elements (type 172 in MAPDL) and extracting SMISC 5 to get the contact pressure.
• • the node coordinates in deformed position for the profile

• The abscissa and pressure values are then exported to a file in the user_files directory of the project
  
  

• To plot the graph in Mechanical, several options are possible. One way is to use matplotlib to do this, but matplotlib is a CPython library whereas Python Code operates in IronPython. The trick is to write a python script in the user_files directory that will be executed by a CPython engine thanks to subprocess.
  
  This script generates an image in the user_files directory
  
  
  and Mechanical can then grab this image and add it to the Mechanical tree.
  
  The code is a bit long so I will have to post this in several comments because of restrictions in the number of characters.

Pernelle Marone-Hitz

Code up to helper functions:

def after_post(this, solution):# Do not edit this line
    """
    Called after post processing.
    Keyword Arguments : 
        this -- the datamodel object instance of the python code object you are currently editing in the tree
        solution -- Solution
    """

    import math
    import os
    import mech_dpf
    import Ans.DataProcessing as dpf

    def get_user_files_directory():
        user_files = ''
        if ExtAPI.Context == 'Project':
            user_files = GetUserFilesDirectory()
        elif ExtAPI.Context == 'Mechanical':
            import wbjn
            user_files = wbjn.ExecuteCommand(ExtAPI, "returnValue(GetUserFilesDirectory())")
        return user_files

    def findClosestNode(ref_coord,node_coords):
        mindist=1e31
        minnode=-1
        for nd,coord in node_coords.items():
            dist=(ref_coord[0]-coord[0])**2+(ref_coord[1]-coord[1])**2+(ref_coord[2]-coord[2])**2
            if dist<mindist:
                mindist=dist
                minnode=nd
        return minnode,mindist

    def getCurvilinearAbscissa(node_coords):
        local_node_coords={nid:coord for nid,coord in node_coords.items()}
        curv_absc=[]
        distances={}
        # Find maximum node location in Y direction
        ymax=-1e31
        for nd,coord in local_node_coords.items():
            if ymax<coord[1]:
                ymax=coord[1]
        # Filter nodes to those at Ymax - this will focus the results to the target metal sheet 
        # where maximum pressure should occur
        to_remove=[]
        for ni,coord in local_node_coords.items():
            if abs(coord[1]-ymax)>5e-3:
                to_remove.append(ni)
        for ni in to_remove:                
            del local_node_coords[ni]
        # Now find chain of nodes from lowest X to largest X and
        # compute distance between nodes
        # find node with lowest X
        xmin=1e31
        for nd,coord in local_node_coords.items():
            if xmin>coord[0]:
                xmin=coord[0]
                ndmin=nd
        dist_idx=0
        ref_coord=local_node_coords[ndmin]
        distances[dist_idx]=[ndmin,0.]
        dist_idx+=1
        del local_node_coords[ndmin]
        # Now compute distance between nodes
        while len(local_node_coords.keys())>0:
            next_node,dist=findClosestNode(ref_coord,local_node_coords)
            distances[dist_idx]=[next_node,math.sqrt(dist)]
            ref_coord=node_coords[next_node]
            del local_node_coords[next_node]
            dist_idx+=1
        curv_absc=[[distances[0][0],xmin]]
        # Create curvilinear abscissa by adding distances between nodes
        for i in range(1,max(distances.keys())+1):
            curv_absc.append([distances[i][0],curv_absc[-1][1]+distances[i][1]])
        return curv_absc

Pernelle Marone-Hitz

Main part of the code (up to exporting the values in text format):

    analysis=this.Parent.Parent
    mech_dpf.setExtAPI(ExtAPI)

    dataSource = dpf.DataSources(analysis.ResultFileName)

    # Retrieve scoping for all contact elements (type 172 in MAPDL)
    scoping_contact_elem=dpf.operators.scoping.on_property(requested_location='Elemental',
        property_name='mapdl_element_type',property_id=172,data_sources=dataSource)

    # Nodes from user named selection
    scoping_from_ns=dpf.operators.scoping.on_named_selection(requested_location='Nodal',
        named_selection_name='EDGEFORPLOT',data_sources=dataSource,int_inclusive=0)

    # Actual scopings
    elem_scoping=scoping_contact_elem.outputs.mesh_scoping.GetData()
    nodal_scoping=scoping_from_ns.outputs.mesh_scoping.GetData()

    # Operator to retrieve contact pressure (SMISC 5 for conta172 elements)
    cont_pres_elem= dpf.operators.result.mapdl.smisc(data_sources=dataSource,mesh_scoping=elem_scoping,item_index=5)
    cont_pres_nodal=dpf.operators.averaging.elemental_nodal_to_nodal_fc(cont_pres_elem)

    # Retrieve location of user_files directory for this project
    userFilesDir = get_user_files_directory()

    # Get pressure for all contact elements
    pressure_all_contacts=cont_pres_nodal.outputs.getfields_container()[0]

    # Rescope to nodes from selections
    pressure_selection_op=dpf.operators.scoping.rescope(fields=pressure_all_contacts,
        mesh_scoping=nodal_scoping)

    # Now this is the pressure field we want to plot over curvilinear abscissa
    pl=pressure_selection_op.outputs.getfields_as_field()

    # Get node ids and pressure values
    nodes=pl.ScopingIds
    values=pl.Data

    # Find node coordinates (in deformed position1)
    value_by_node={}
    node_coords={}
    rst_mesh=dpf.Model(dataSource).Mesh
    # Displacement field to update node coordinates with actual displacements
    u=dpf.operators.result.displacement(data_sources=dataSource,mesh_scoping=nodal_scoping)
    u_field=u.outputs.getfields_container()[0]

    for i in range(values.Count):
        val=values[i]
        value_by_node[nodes[i]]=val
        u_node=u_field.GetEntityDataById(nodes[i])
        nd=rst_mesh.NodeById(nodes[i])
        node_coords[nodes[i]]=[nd.X+u_node[0],nd.Y+u_node[1],nd.Z+u_node[2]]

    # Compute curvilinear abscissa
    curv_abs=getCurvilinearAbscissa(node_coords)

    # Reorder node coordinates on X values
    nids=[]
    xcoords=[]
    ycoords=[]
    for nid, coord in node_coords.items():
        nids.append(nid)
        xcoords.append(coord[0])
        ycoords.append(coord[1])

    sorted_coords = zip(xcoords,ycoords)
    sorted_coords.sort()


    # Export pressure vs abscissa
    myPath = os.path.join(userFilesDir,'pressure_on_path_{0}.txt'.format(dp))
    with open(myPath,'w') as f:
        f.write('Abscissa\tPressure}\n')
        for ca in curv_abs:
            f.write('{0:12.6e}\t{1:12.6e}\n'.format(ca[1],value_by_node[ca[0]]))
        # In same file, export node profile after deformation
        f.write('$\n')
        for row in sorted_coords:
            f.write('{0:12.6e}\t{1:12.6e}\n'.format(row[0],row[1]))

Pernelle Marone-Hitz

Code to create the plot with Matplotlib and add the image to Mechanical:

    # Create plot using CPython code
    pngPath=os.path.join(userFilesDir,'pressure_on_path_{0}.png'.format(dp))

    c_python_script='''    
import os
import matplotlib.pyplot as plt

with open(r"'''+myPath+'''","r") as f:
    x_for_plot=[]
    y_for_plot=[]

    x_profile=[]
    y_profile=[]

    line=f.readline()
    line=f.readline()
    while not('$' in line):
        sp_line=line.strip("\\n").split("\\t")
        x_for_plot.append(float(sp_line[0]))
        y_for_plot.append(float(sp_line[1]))    
        line=f.readline()

    line=f.readline()
    while line:
        sp_line=line.strip("\\n").split("\\t")
        x_profile.append(float(sp_line[0]))
        y_profile.append(float(sp_line[1]))    
        line=f.readline()

fig, ax = plt.subplots(figsize = (10, 5))    

# using the twinx() for creating another
# axes object for secondary y-Axis
ax2 = ax.twinx()
ax.plot(x_for_plot,y_for_plot,'b')
ax2.plot(x_profile,y_profile,'r-')

# giving labels to the axises
ax.set_xlabel('abscissa')
ax.set_ylabel('pressure',color='b')
ax2.set_ylabel('y location of node',color='r')

plt.tight_layout()

myPath = r"'''+pngPath+'''"
fig.savefig(myPath)
'''    

    scriptPath = os.path.join(userFilesDir,'generate_plot.py'.format(dp))
    with open(scriptPath,'w') as f:
        f.write(c_python_script)

    import sys,subprocess
    from System.Diagnostics import Process

    currentVersion=DataModel.Project.ProductVersion
    versionSplit=currentVersion.split(' ')    
    root_folder='AWP_ROOT'+versionSplit[0][2:]+versionSplit[-1][-1]    
    path1 = os.path.join( os.getenv(root_folder),r'commonfiles\CPython\3_7\winx64\Release\python', 'python.exe') # path to Mechanical's Python
    path2 = scriptPath
    # Argument for Python.exe call
    CMDARGS = '"{0}"'.format(path2)
    # Command for python exe including path to exe
    command = '"{0}"'.format(path1)

    ExtAPI.Log.WriteMessage(command)
    ExtAPI.Log.WriteMessage(CMDARGS)
    # Launch command with arguments
    proc = Process()
    proc.StartInfo.FileName = command
    proc.StartInfo.Arguments = CMDARGS
    proc.Start()
    proc.WaitForExit()

    # Add plot to Mechanical Tree
    with Transaction():
        images=ExtAPI.DataModel.GetObjectsByType(Ansys.ACT.Automation.Mechanical.Image)
        try:
            for im in images:
                if im.Name=="Pressure on path":
                    im.Delete()
        except:
            pass
        image = solution.AddImage(pngPath)
        image.Name="Pressure on path"