Map body force density on Mechanical Mesh

A Kumar

Map body force density on Mechanical Mesh

A Kumar

import os

import numpy as np

import ansys.dpf.core as dpf


server = dpf.start_local_server(ansys_path=r"C:\Program Files\ANSYS Inc\v252")

print(server)


def read_lorentz_forces(fld_path: str):


    data = np.loadtxt(fld_path, skiprows=2)   # shape (N, 6)


    coords = data[:, :3]   # x, y, z

    forces = data[:, 3:]   # fx, fy, fz


    n_nodes = data.shape[0]

    node_ids = list(range(1, n_nodes + 1))


    # --- coordinates field ---

    coords_field = dpf.fields_factory.create_3d_vector_field(

        num_entities=n_nodes,

        location=dpf.locations.nodal,

    )

    coords_field.scoping.ids = node_ids

    coords_field.data = coords

    coords_field.unit = "m"


    # source mesh

    source_mesh = dpf.MeshedRegion(num_nodes=len(coords_field.scoping.ids))

    source_mesh.unit = coords_field.unit

    for nid in coords_field.scoping.ids:

        source_mesh.nodes.add_node(nid, coords_field.get_entity_data_by_id(nid))


    # --- force density field ---

    force_density_field = dpf.fields_factory.create_3d_vector_field(

        num_entities=n_nodes,

        location=dpf.locations.nodal,

    )

    force_density_field.scoping.ids = node_ids

    force_density_field.data = forces

    force_density_field.unit = "N/m^3"


    return coords_field, source_mesh, force_density_field


wkdir = r"…"

fld_file = r"…\lorentz_forces.fld"

coords_field, source_mesh, force_density_field = read_lorentz_forces(fld_file)


rst = r"…\file.rst"

model = dpf.Model(rst)


sm = source_mesh

tm = model.metadata.meshed_region


report_path = os.path.join(wkdir, "mapping_report.html")


wf = dpf.Operator("mechanical_native_mapping::prepare::point_cloud")

wf.connect(0, sm)

wf.connect(1, tm)

wf.connect(5, dpf.locations.nodal)

wf.connect(6, 3)

wf.connect(8, report_path)

wf.connect(10, "triangulation")  # optional, default is "triangulation"


s_fc = dpf.FieldsContainer()

s_fc.labels = ["time"]

s_fc.add_field({'time' : 0}, force_density_field)


native_mech_map = dpf.Operator("mechanical_native_mapping::apply")

for index in range(5):

    native_mech_map.connect(index, wf, index)

native_mech_map.connect(5, s_fc)

native_mech_map_f = native_mech_map.eval()[0]


# Get Lorentz forces

lf = dpf.operators.geo.integrate_over_elements(field=native_mech_map_f, mesh=tm)


# Total Lorentz force

lf_sum = dpf.operators.math.accumulate(lf).eval()


print(lf_sum)


# Plot bfd

tm.plot(native_mech_map_f)