Creating a field and plotting over geometry in DPF Mechanical

Hello,

I am trying to create a custom field to map some calculated values over geometry. I want to do this using Python Result. The Python Result has scoping defined by user using Property Provider. Up to operator op_div everything works fine. Then I try to do extra calculations, create a new field and use it for further calculations with operators. I tried a few approaches, I think rbf mapper should work, but I am not sure.

For simple geometry, the code seems to work, please see the picture below:

For main model, it does something weird (I cannot share geometry for this):

The main models consists of multiple bodies connected using bonded contacts together with some frictionless contacts. I also have bolt pretension defined.

The code is below:

def post_started(sender, analysis):# Do not edit this line
    define_dpf_workflow(analysis)

def define_dpf_workflow(analysis):
    import mech_dpf
    import Ans.DataProcessing as dpf
    mech_dpf.setExtAPI(ExtAPI)
    dataSource = dpf.DataSources(analysis.ResultFileName)
    meshData = DataModel.MeshDataByName("Global")
    timeScop = dpf.Scoping()
    timeScop.Ids = [3]
    timeScop.Location='Elemental'

    scoping_refs = this.GetCustomPropertyByPath("Geometry/Scoping Property/Geometry Selection").Value
    nodes_loc = mech_dpf.GetNodeScopingByRefId(scoping_refs)



    #s1= dpf.operators.result.stress_principal_1() 
    #s2= dpf.operators.result.stress_principal_3() 
    #s3= dpf.operators.result.stress_principal_3()

    s1= dpf.operators.result.stress_X() 
    s2= dpf.operators.result.stress_Y() 
    s3= dpf.operators.result.stress_Z() 

    seqv = dpf.operators.result.stress_von_mises()

    Add_1_2 = dpf.operators.math.add()
    Add_1_2_3 = dpf.operators.math.add()
    op_1_3 = dpf.operators.math.scale()
    op_div = dpf.operators.math.component_wise_divide()
    op_mult_scale = dpf.operators.math.scale()
    op_div_2 = dpf.operators.math.component_wise_divide()

    s1.inputs.time_scoping.Connect(timeScop)
    s2.inputs.time_scoping.Connect(timeScop)
    s3.inputs.time_scoping.Connect(timeScop)
    seqv.inputs.time_scoping.Connect(timeScop)

    s1.inputs.data_sources.Connect(dataSource)
    s2.inputs.data_sources.Connect(dataSource)
    s3.inputs.data_sources.Connect(dataSource)
    seqv.inputs.data_sources.Connect(dataSource)

    s1.inputs.mesh_scoping.Connect(nodes_loc)
    s2.inputs.mesh_scoping.Connect(nodes_loc)
    s3.inputs.mesh_scoping.Connect(nodes_loc)
    seqv.inputs.mesh_scoping.Connect(nodes_loc)

    Add_1_2.inputs.fieldA.Connect(s1)
    Add_1_2.inputs.fieldB.Connect(s2)
    s_1_2 = Add_1_2.outputs.field.GetData()

    Add_1_2_3.inputs.fieldA.Connect(s_1_2)
    Add_1_2_3.inputs.fieldB.Connect(s3)
    s_1_2_3 = Add_1_2_3.outputs.field.GetData()

    op_1_3.inputs.field.Connect(s_1_2_3)
    op_1_3.inputs.ponderation.Connect(0.3333333333)
    #op_1_3.inputs.boolean.Connect(my_boolean)# optional
    sigmaH = op_1_3.outputs.field.GetData()

    op_div.inputs.fieldA.Connect(sigmaH)
    op_div.inputs.fieldB.Connect(seqv)
    H = op_div.outputs.field.GetData()

    Min_tol_strain = 0.05
    Rm=470 #Mpa
    Rp = 345 #MPa
    hard_factor = 0.5*(1+Rm/Rp)
    A = 0.22
    E = 210000 #MPa
    ref_strain = A+Rp/E
    tol_strain = 0

    npl_field = dpf.FieldsFactory.CreateScalarField(H.Data.Count, "Nodal")

    input_support = dpf.FieldsFactory.Create3DVectorField(H.Data.Count, "Nodal")
    input_field = dpf.FieldsFactory.CreateScalarField(H.Data.Count, "Nodal")
    lineIndex=1

    msg=Ansys.Mechanical.Application.Message(str(H.Data.Count),MessageSeverityType.Info)
    ExtAPI.Application.Messages.Add(msg)

    for node_Id in H.Scoping.Ids:
        if H.GetEntityDataById(node_Id) <=1/3:
            tol_strain = ref_strain
        else: tol_strain = Min_tol_strain + ((ref_strain-Min_tol_strain)/3)^(3*H.GetEntityDataById(node_Id))
        node = meshData.NodeById(node_Id) 
        input_support.Add(lineIndex, [node.X, node.Y, node.Z])
        # npl_field.UpdateEntityDataByEntityIndex(i,MIN(sqrt(E*tol_strain/Rp),hard_factor))
        input_field.Add(lineIndex,[min(sqrt(E*tol_strain/Rp),hard_factor)])
        lineIndex += 1
        #npl_field.UpdateEntityDataByEntityIndex(i,10.00)

    #DPF Mapper Operator
    rbf_prep = dpf.Operator("make_rbf_mapper")
    rbf_prep.Connect(0, input_support)
    rbf_apply = dpf.Operator("apply_mapping")
    rbf_apply.Connect(0, rbf_prep)
    rbf_apply.Connect(1, input_field)
    rbf_apply.Connect(2, input_support)
    #Radius to read results for mapping
    r = 0.01
    rbf_prep.Connect(1, r)
    npl_field = rbf_apply.GetOutputAsField(0)

    print(npl_field.Data)

    msg=Ansys.Mechanical.Application.Message(str(H.GetEntityDataByIndex(10))+" "+ str(npl_field.GetEntityDataByIndex(10)),MessageSeverityType.Info)
    ExtAPI.Application.Messages.Add(msg)

    op_mult_scale.inputs.field.Connect(npl_field)
    op_mult_scale.inputs.ponderation.Connect(Rp)
    sigmaSK = op_mult_scale.outputs.field.GetData()

    op_div_2.inputs.fieldA.Connect(seqv)
    op_div_2.inputs.fieldB.Connect(sigmaSK)
    UC = op_div_2.outputs.field.GetData()

    # timeScop = dpf.Scoping()
    # timeScop.Ids = [1]
    # u.inputs.time_scoping.Connect(timeScop)
    dpf_workflow = dpf.Workflow()
    dpf_workflow.Add(s1)
    dpf_workflow.Add(s2)
    dpf_workflow.Add(s3)
    dpf_workflow.Add(seqv)
    dpf_workflow.Add(Add_1_2)
    dpf_workflow.Add(Add_1_2_3)
    dpf_workflow.Add(op_1_3)
    dpf_workflow.Add(op_div)
    dpf_workflow.Add(op_mult_scale)
    dpf_workflow.Add(op_div_2)

    #dpf_workflow.Add(nrm)
    # dpf_workflow.SetInputName(u, 0, 'time')
    # dpf_workflow.Connect('time', timeScop)
    dpf_workflow.SetOutputContour(op_div_2)
    dpf_workflow.Record('wf_id', False)
    this.WorkflowId = dpf_workflow.GetRecordedId()

The main idea is that I get node Ids for field H. Then I loop over all node Ids in H, create input support field with coordinates of the nodes, populate input_field with values. Then outside of the loop use dpf mapper operator to create npl_field with output_support equal to input_support (I am mapping to the same locations as input) with data calculated in the for loop. Then I use field npl_field with new operators and at the end I want to plot the result.

Could you please advise if the code above makes sense and suggest what I should change to make it work? The most confusing thing for me is that for larger model this doesn't work but seems to work for simpler geometry.

PS MsgBox were used to check a few things, I don't need them in the final code.

Answers