How do I find the time varying center of gravity location of my body using DPF ?

Vishnu

The COG of a body might vary with time as the body deforms and you might be interested in knowing how to compute the time varying COG.

the example also highlights on how you can create your own mesh with DPF.

Example uses SOLID elements in this case to demonstrate and used DPF inside mechanical

Vishnu

NOTE : You will need 2024R1 version of Ansys for this to work

#Inputs
time_point = 3
named_selection_string = "MYBODY"

#Import statements
import mech_dpf
import Ans.DataProcessing as dpf

# Result Data
analysis1 = ExtAPI.DataModel.Project.Model.Analyses[0]
dataSource = dpf.DataSources(analysis1.Solution.ResultFilePath)

# Get mesh
my_model = dpf.Model(analysis1.Solution.ResultFilePath)
my_mesh = my_model.Mesh

# Time list
timelist = dpf.operators.metadata.time_freq_provider(data_sources=dataSource).outputs.gettime_freq_support().TimeFreqs.Data
time_index_list = [int(timelist.IndexOf(time_point)+1)]

# Named Selection Mesh scoping
ns_op = dpf.operators.scoping.on_named_selection(requested_location='Elemental', named_selection_name=named_selection_string, data_sources=dataSource)
mymeshscoping = ns_op.outputs.getmesh_scoping()

# Get Deformation
u = dpf.operators.result.displacement(time_scoping=time_index_list, data_sources=dataSource)
disp_field = u.outputs.fields_container.GetData()[0]

# Create my own mesh
created_mesh = dpf.MeshedRegion(my_mesh.Nodes.Count, 0)
for node in my_mesh.Nodes:
    displacement = disp_field.GetEntityDataById(node.Id)
    created_mesh.AddNode(node.Id, [node.X + displacement[0], node.Y + displacement[1], node.Z + displacement[2]])

createdmesh_nodeids = created_mesh.NodeIds
for elem in my_mesh.Elements:
    created_mesh.AddSolidElement(elem.Id, list(map(lambda x: createdmesh_nodeids.IndexOf(x), elem.CornerNodeIds)))

# Centroids
centroids = dpf.Operator("compute_element_centroids")
centroids.Connect(1, mymeshscoping)
centroids.Connect(7, created_mesh)
centroid = centroids.GetOutputAsField(0)


def compute_elemental_mass():
    # List to store elemental densities
    all_mats = []

    mats = my_mesh.GetPropertyField("mat")
    mat_prop = my_model.CreateOperator("mapdl_material_properties")
    mat_prop.inputs.materials.Connect(mats)

    mat_prop.inputs.properties_name.Connect("DENS")
    mat_field = mat_prop.outputs.properties_value.GetData()[0]
    for eid in my_mesh.ElementIds:
        mat_id = mats.GetEntityDataById(eid)
        density = mat_field.GetEntityDataById(mat_id[0])
        all_mats.append(density[0])

    # Create Density Field
    density = dpf.FieldsFactory.CreateScalarField(numEntities=my_mesh.ElementCount,location=dpf.locations.elemental)
    density.ScopingIds = my_mesh.ElementIds
    density.Data = all_mats

    # Element volume
    el_vol = dpf.operators.geo.elements_volume(mesh=my_mesh)

    # Calculate elemental mass (density * volume)
    elemental_mass = dpf.operators.math.generalized_inner_product(fieldA=el_vol, fieldB=density).outputs.getfield()
    return elemental_mass

elemental_mass = compute_elemental_mass()


#Function to calculate mass-weighted centroid
def calculate_mass_weighted_centroid(centroid, mass):
    num_points = centroid.ElementaryDataCount

    # Extract x, y, z components
    x_comp = dpf.operators.logic.component_selector(field=centroid, component_number=0).outputs.getfield().Data
    y_comp = dpf.operators.logic.component_selector(field=centroid, component_number=1).outputs.getfield().Data
    z_comp = dpf.operators.logic.component_selector(field=centroid, component_number=2).outputs.getfield().Data

    # Calculate mass-weighted sums
    sum_mass = sum(mass.Data)
    sum_x = sum(x * m for x, m in zip(x_comp, elemental_mass.Data))
    sum_y = sum(y * m for y, m in zip(y_comp, elemental_mass.Data))
    sum_z = sum(z * m for z, m in zip(z_comp, elemental_mass.Data))

    # Calculate mass-weighted centroid
    if sum_mass > 0:
        centroid_x = sum_x / sum_mass
        centroid_y = sum_y / sum_mass
        centroid_z = sum_z / sum_mass
    else:
        raise ValueError("Total mass is zero or negative")

    return (centroid_x, centroid_y, centroid_z)

#Print mass-weighted centroid
print(calculate_mass_weighted_centroid(centroid,elemental_mass))

Paul Profizi

Hi @Vishnu, thank you for this example!

One remark though:

This is tagged as "PyDPF", yet as it is using the API for use within Mechanical, I am afraid this will create confusion for people trying to use this with PyDPF-Core/Post. Is there a better suited tag available? Or maybe we should create one specific for use of DPF from within Mechanical, be it in Python.

Vishnu

I think so too we may need to create one for DPF from within Mechanical

James Derrick

Hey! I added the DPF tag to your post! Also can you confirm that COG = "Centre of Gravity"? I don't think you state it anywhere.