
###############################################################################
# Post-processing of the results using PyDPF
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# PyDPF can be used in this context to visualise the crack opening throughout
# the simulation as an animation. The steps for producing that are provided in
# below.

# Generating DPF model
data_src = dpf.DataSources(os.path.join(mapdl.directory, "file.rst"))
model = Model(data_src)

# Getting mesh of the model
meshed_region = model.metadata.meshed_region
Mesh_field = meshed_region.field_of_properties(dpf.common.nodal_properties.coordinates)

# Index of NMISC results
nmisc_index = 70

# Generating damage object
dam_op = dpf.operators.result.nmisc()
dam_op.inputs.data_sources(data_src)

# Generating displacement object
disp_op = dpf.Operator("U")
disp_op.inputs.data_sources.connect(model.metadata.data_sources)

# Summing displacements to obtain total at step n
op = dpf.operators.math.add()

my_fieldA = dpf.Field()
my_fieldA.location = locations.nodal
op.inputs.fieldA.connect(Mesh_field)

# Generating the plotter object
plotter = pv.Plotter()
plotter.open_gif("dcb.gif")

# Generating grid of the whole beam
mesh_beam = meshed_region.grid
plotter.add_mesh(
    mesh_beam,
    lighting=False,
    show_edges=True,
    scalar_bar_args={"title": "Damage"},
    clim=[0, 1],
    opacity=0.3,
)

# Generating mesh of the contact part only 
# mesh_contact = ? # We need a way to get the mesh of the contact results only named selection "_ELEMCM"
# plotter.add_mesh(
#         mesh_contact,
#         opacity=0.3,
#         scalar_bar_args={"title": "Damage"},
#         clim=[0, 1]
#     )


for i in range(1, 100):

    # Getting displacements
    disp = model.results.displacement(time_scoping=i).eval()

    my_fieldB = dpf.Field()
    my_fieldB.location = locations.nodal
    op.inputs.fieldB.connect(disp[0])

    # Getting sum
    result_field = op.outputs.field()

    plotter.update_coordinates(result_field.data, mesh=mesh_beam, render=False)
    # plotter.update_scalars(disp_values, mesh=mesh_beam, render=False)

    # pLotting contact
    # Setting time step
    dam_op.inputs.time_scoping([i])
    dam_op.inputs.item_index.connect(nmisc_index)

    # Getting the values
    # dam_op ...
    # ...
    # result_contact =

    # Updating
    # plotter.update_scalars(result_contact, mesh=mesh_contact, render=False)

    plotter.write_frame()

plotter.show()

