How to create beam bending moment (BMD) and shear force diagram (SFD) using PyMAPDL and PyDPF

Inputs and import statements

#Import statements
from ansys.mapdl.core import launch_mapdl
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.cm import ScalarMappable
from matplotlib.colors import Normalize
from matplotlib import ticker
from ansys.dpf import core as dpf
import os
run_location=r"D:\apdl\test\4"
nproc=4
override=True
mapdl = launch_mapdl(run_location=run_location, nproc=nproc, override=override)

Creation of model using pyMAPDL

def create_model_setup(SECDATA, 
                       BEAM_LENGTH, 
                       LINE_DIV,
                       LENGTH_DIV, 
                       YM, 
                       MU, 
                       SECTYPE, 
                       SECSUBTYPE):
    mapdl.clear("All")
    mapdl.prep7()
    mapdl.et('', 'BEAM188')
    mapdl.mp(lab='EX', mat='1', c0='%s'%(YM))
    mapdl.mp(lab='PRXY', mat='1', c0='%s'%(MU))
    mapdl.sectype(secid='1', type_='BEAM', subtype='%s'%(SECTYPE))
    mapdl.secdata(tuple(SECDATA))
    mapdl.seccontrol(1E15,1E15)
    mapdl.k(npt='1', x='0', y='0', z='0')
    mapdl.k(npt='2', x='%s'%(BEAM_LENGTH), y='0', z='0')
    mapdl.l(p1='1', p2='2', ndiv='%s'%(LINE_DIV))
    mapdl.lesize(nl1='ALL', size='', angsiz='', ndiv='%s'%(LENGTH_DIV))
    mapdl.lmesh(nl1='ALL')

Adding Loads and Bc and Solving them

def run_simulation(fix_loc,
                   pin_loc, 
                   concentrated_loc, 
                   concentrated_force,
                   distributed_loc,
                   distributed_force,
                   elem_plot_flag):
    mapdl.solution()
    mapdl.ddele(node='ALL', lab='ALL')
    mapdl.fdele(node='ALL', lab='ALL')
    mapdl.sfedele(elem='ALL', lkey='ALL', lab='ALL')
    for fix_lo in fix_loc:
        fix_node = mapdl.queries.node(fix_lo[0], fix_lo[1], fix_lo[2])
        mapdl.d(node='%s'%(fix_node), lab='ALL', value='0')
    for pin_lo in pin_loc:
        pin_node = mapdl.queries.node(pin_lo[0], pin_lo[1], pin_lo[2])
        mapdl.d(node='%s'%(pin_node), lab='UX', value='0')
        mapdl.d(node='%s'%(pin_node), lab='UY', value='0')
        mapdl.d(node='%s'%(pin_node), lab='UZ', value='0')
        mapdl.d(node='%s'%(pin_node), lab='ROTY', value='0')
        mapdl.d(node='%s'%(pin_node), lab='ROTZ', value='0')
    for count, force_lo in enumerate(concentrated_loc):
        force_node = mapdl.queries.node(force_lo[0], force_lo[1], force_lo[2])
        mapdl.f(node='%s'%(force_node), lab='FY', value=concentrated_force[count])
    for count, force_lo in enumerate(distributed_loc):
        mapdl.esel(type_='S', item='CENT', comp='X', vmin=force_lo[0], vmax=force_lo[1])
        mapdl.sfbeam(elem='ALL', lkey='2', lab='PRES', vali=distributed_force[count])   
        mapdl.allsel(labt='ALL', entity='ALL')
    mapdl.allsel(labt='ALL', entity='ALL')
    mapdl.outres(item='ALL', freq='ALL')
    if elem_plot_flag == True:
        mapdl.eplot(plot_bc=True, bc_labels=["UX", "UY", "UZ", "FX", "FY", "FZ"], plot_bc_labels=True)
    mapdl.solve()
    mapdl.save(fname='test', ext='db')

Extracting shear force and bending moments using DPF and plot using matplotlib

def plot_bending_moment_shear_force(rstpath, 
                                    smiscI, 
                                    smiscJ, 
                                    plotname):
    # data source
    my_data_sources = dpf.DataSources(rstpath)
    # model
    my_model = dpf.Model(rstpath)
    # mesh
    my_mesh = my_model.metadata.meshed_region
    # Get result operator
    result_operator = dpf.operators.result.__getattribute__("smisc")
    result_field_container = result_operator(
        time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh, item_index=smiscI
    ).outputs.fields_container.get_data()
    # First Field usually field corresponding to first time scoping
    resultfield = result_field_container[0]
    element_ids = resultfield.scoping.ids
    i_nodes_coordinates = [my_mesh.elements.element_by_id(
        x).nodes[1].coordinates[0] for x in element_ids]
    BMY_T_I = resultfield.data_as_list
    result_field_container = result_operator(
        time_scoping=[1], data_sources=my_data_sources, mesh=my_mesh, item_index=smiscJ
    ).outputs.fields_container.get_data()
    # First Field usually field corresponding to first time scoping
    resultfield = result_field_container[0]
    element_ids = resultfield.scoping.ids
    j_nodes_coordinates = [my_mesh.elements.element_by_id(
        x).nodes[1].coordinates[0] for x in element_ids]
    BMY_T_J = resultfield.data_as_list
    # Assuming i_nodes_coordinates, j_nodes_coordinates, BMY_T_I, BMY_T_J are defined
    x = np.concatenate((i_nodes_coordinates, j_nodes_coordinates))
    y = np.concatenate((BMY_T_I, BMY_T_J))
    # Create a colormap
    cmap = plt.get_cmap('viridis')
    norm = Normalize(vmin=np.min(y), vmax=np.max(y))
    # Create a ScalarMappable
    sm = ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])
    # Create the bar chart
    fig, ax = plt.subplots(figsize=(15, 5))  # Adjust the figure size
    ax.bar(x, y, color=cmap(norm(y)))
    ax.set_xlabel('X_distance')
    ax.set_title(plotname+' Diagram')
    # Set y-axis limits
    ax.set_ylim([np.min(y), np.max(y)])
    # Create an array of 10 evenly spaced values between min(y) and max(y)
    ticks = np.linspace(np.min(y), np.max(y), 10)
    # Add a colorbar with specific ticks
    cbar = fig.colorbar(sm, ax=ax, ticks=ticks)
    cbar.set_label(plotname)
    # Format the colorbar tick labels to have two decimal places
    cbar.formatter = ticker.FuncFormatter(lambda x, pos: '{:.4f}'.format(x))
    cbar.update_ticks()
    plt.show()

Function calls Below

### Call the defined functions####
# create model
create_model_setup(
    SECDATA=[5, 8, 2],  # SECDATA = [BEAM_RADIUS,CIRCUM_DIV,RADIAL_DIV]
    BEAM_LENGTH=100,
    LINE_DIV=50,
    LENGTH_DIV=1,
    YM=2E8,
    MU=0.3,
    SECTYPE="CSOLID",
    SECSUBTYPE="CIRBEAM",
)
#Run simulation
run_simulation(
    fix_loc=[[10, 0, 0], [90, 0, 0]],
    pin_loc=[[50, 0, 0]],
    concentrated_loc = [[30, 0, 0], [80, 0, 0],[100, 0, 0]],
    concentrated_force = [1000, 500,500],
    distributed_loc=[[0,100]],
    distributed_force = [1000],
    elem_plot_flag = True
)
# Bending Moment diagram
plot_bending_moment_shear_force(
    rstpath=os.path.join(run_location, "file.rst"),
    smiscI=32,
    smiscJ=37,
    plotname="Bending Moment"
)
# Shear Force diagram
plot_bending_moment_shear_force(
    rstpath=r"D:\apdl\test\4\file.rst",
    smiscI=6,
    smiscJ=19,
    plotname="Shear Force"
)
mapdl.exit()

Also for further reference please refer BEAM188 documentation below,
https://ansyshelp.ansys.com/account/secured?returnurl=/Views/Secured/corp/v231/en/ans_elem/Hlp_E_BEAM188.html

@deepak
In recent releases its
dpf.operators.result.smisc
in earlier releases its
dpf.operators.result.mapdl.smisc

Please check these