Create custom field in Mechanical DPF

Hello,

I have a script which I would like to use to plot certain parameter I calculate. To do that, I figured I could create a field and populate it in the loop. The script seems to work for simple geometry, please see picture below:

However, when I run it on actual geometry I get this

The script uses Property Provided to provide geometry selected by the user, in this case a single body. For some reason, when script is run on the more complicated geometry I get more bodies than 1 (originally the bottom plate is selected) and no results are displayed (they will not be 0).

Please see the script below:

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

# Uncomment this function to enable retrieving results from the table/chart
# def table_retrieve_result(value):# Do not edit this line
    # import mech_dpf
    # import Ans.DataProcessing as dpf
    # wf = dpf.Workflow(this.WorkflowId)
    # wf.Connect('contour_selector', value)
    # this.Evaluate()

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



    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)


    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()

Operator op_div seems to display correctly so there is something going on with npl_field. To be honest, I am not sure if rbf mapper is the way to go but it makes sense in my head:
1. Create a field of nodes positions input_support (will be the same nodes as other H field's scoping
2. Loop over all the nodes and populate input_field
3. Outside of the loop, use rbf mapper to create npl_field. Output support is equal to input support (I am mapping over identical locations).
4. Feed npl_field into some operators to finish of the calculations.
5. Plot op_div_2 operator results.

Can you please help me figure out if this is the correct approach and indicate why this script doesn't work for larger geometry? Larger geometry also has contacts and preloaded bolts so I wonder if this has some effect? Unfortunately I cannot share geometry model.

Best Answer

  • Rohith Patchigolla
    Rohith Patchigolla Member, Moderator, Employee Posts: 178
    100 Comments 25 Answers 25 Likes First Anniversary
    ✭✭✭✭
    Answer ✓

    Hello @ArnoldH

    Please find a modified script here for your reference.
    Some comments on changes:
    • When defining integer numbers which may be involved in a division later on, always add .0 after the number (Eg: Rm = 470.0), else the division will lead to rounded integer
    • Time scoping corresponds to result set number and not the time.
    • rbf map operator is not needed
    • To plot a field created in a python result, please use forward_field operator
    • For scoping of the new field, I used the scoping IDs coming in directly from H to maintain the correct order of nodes.

    def post_started(sender, analysis):# Do not edit this line
        define_dpf_workflow(analysis)
    
    # Uncomment this function to enable retrieving results from the table/chart
    # def table_retrieve_result(value):# Do not edit this line
        # import mech_dpf
        # import Ans.DataProcessing as dpf
        # wf = dpf.Workflow(this.WorkflowId)
        # wf.Connect('contour_selector', value)
        # this.Evaluate()
    
    def define_dpf_workflow(analysis):
        import mech_dpf
        import Ans.DataProcessing as dpf
        import math
    
        mech_dpf.setExtAPI(ExtAPI)
        dataSource = dpf.DataSources(analysis.ResultFileName)
    
        timeScop = dpf.Scoping()
        timeScop.Ids = [3]
    
        scoping_refs = this.GetCustomPropertyByPath("Geometry/Scoping Property/Geometry Selection").Value
        nodes_loc = mech_dpf.GetNodeScopingByRefId(scoping_refs)
    
        #ns_op = dpf.operators.scoping.on_named_selection(data_sources = dataSource, requested_location = "Nodal", named_selection_name = "SELECTION")
    
        s1= dpf.operators.result.stress_X(data_sources = dataSource, time_scoping = timeScop, mesh_scoping = nodes_loc) 
        s2= dpf.operators.result.stress_Y(data_sources = dataSource, time_scoping = timeScop, mesh_scoping = nodes_loc) 
        s3= dpf.operators.result.stress_Z(data_sources = dataSource, time_scoping = timeScop, mesh_scoping = nodes_loc) 
        seqv = dpf.operators.result.stress_von_mises(data_sources = dataSource, time_scoping = timeScop, mesh_scoping = nodes_loc)
    
        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()
        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)
        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.0 #Mpa
        Rp = 345.0 #MPa
        hard_factor = 0.5*(1+Rm/Rp)
        A = 0.22
        E = 210000 #MPa
        ref_strain = A+Rp/E
        tol_strain = 1
    
        npl_field = dpf.FieldsFactory.CreateScalarField(H.Data.Count, "Nodal")
        npl_field_Data = []
        for node_Id in H.Scoping.Ids:
            if H.GetEntityDataById(node_Id)[0] <=1/3:
                tol_strain = ref_strain
            else:
                tol_strain = Min_tol_strain + math.pow((ref_strain-Min_tol_strain)/3,(3*H.GetEntityDataById(node_Id)[0]))
            npl_field_Data.append(min(sqrt(E*tol_strain/Rp),hard_factor))
    
        npl_field.Data = npl_field_Data
        #npl_field.Data = H.Data
    
        my_scoping = dpf.Scoping()
        my_scoping.Location = "Nodal"
        my_scoping.Ids = H.Scoping.Ids
        npl_field.Scoping = my_scoping
    
        # forward field to an operator to plot it
        forward = dpf.operators.utility.forward_field()
        forward.inputs.field.Connect(npl_field)
    
        dpf_workflow = dpf.Workflow()
        dpf_workflow.Add(forward)
        # dpf_workflow.SetInputName(u, 0, 'time')
        # dpf_workflow.Connect('time', timeScop)
        dpf_workflow.SetOutputContour(forward)
        dpf_workflow.Record('wf_id', False)
        this.WorkflowId = dpf_workflow.GetRecordedId()
    

Answers