MeshPy with Potential Interactions

[davidrudlstorfer]

Currently I am further investigating the usage of MeshPy for the creation of multiple helices (based on #60) within a single input file including beam potential interactions.

In a simplified setting I have multiple helices with different start points, different axes, different twist angles and differing heights which are combined into a single 4C input file. Further, I want to enable potential interactions between individual fibers which results in the necessity to define each helix as an individual line.

As a simple proof of concept I use the following Python Code.

Code

import numpy as np

from meshpy import (
  Mesh,
  InputFile,
  MaterialReissner,
  Beam3rHerm2Line3,
)

from meshpy.mesh_creation_functions.beam_basic_geometry import create_beam_mesh_helix


def main():
  """Minimal example of multiple helical beams"""

  # define input file
  input_file = InputFile()

  # define material
  mat = MaterialReissner(youngs_modulus=1000, radius=0.5, shear_correction=1.0)

  # define radii of differing helical segments
  radii = np.array([0, 3, 6, 9, 12, 15, 18])

  # determine twist angles if helical segments
  outermost_twist_angle = np.deg2rad(45)
  twist_angles = []
  for radius in radii:
      if np.isclose(radius, 0.0):
          twist_angles.append(np.pi / 2)
          continue
      twist_angles.append(
          np.arctan((radii[-1] / radius) * np.tan(outermost_twist_angle))
      )

  # create helices
  for radius, twist_angle in zip(radii, twist_angles):

      tmp_mesh = Mesh()

      helix_set = create_beam_mesh_helix(
          tmp_mesh,
          Beam3rHerm2Line3,
          mat,
          [0, 0, 1],  # axis
          [0, 0, 0],  # axis_point
          [radius, 0, radius],  # start_point
          twist_angle,  # twist_angle
          height_helix=2 * np.pi * 10,
          n_el=20,
          warning_straight_line=False,
      )

      input_file.add_mesh(tmp_mesh)
      input_file.add(helix_set["line"])

  input_file.write_input_file("_test_output/test.dat", header=False)
  input_file.write_vtk("test", "_test_output/")

  input_file.display_pyvista(
      beam_nodes=True,
      beam_tube=True,
      beam_cross_section_directors=True,
      resolution=50,
      parallel_projection=True,
  )


if __name__ == "__main__":

  main()

This provides the necessary geometries and also adds each individual helix as a DLINE-NODE TOPOLOGY. To enable potential interactions within 4C I additionally need to add the header settings for the potential interactions:

--------------------------------------------------------------BEAM POTENTIAL
POT_LAW_PREFACTOR               -1e-03 1e-06
POT_LAW_EXPONENT                6.0 12.0
BEAMPOTENTIAL_TYPE              Volume
CUTOFF_RADIUS                   30
STRATEGY                        SingleLengthSpecific_SmallSepApprox_Simple
NUM_INTEGRATION_SEGMENTS        1
NUM_GAUSSPOINTS                 50
POTENTIAL_REDUCTION_LENGTH      1.5
AUTOMATIC_DIFFERENTIATION       no
CHOICE_MASTER_SLAVE             higher_eleGID_is_slave
----------------------------------BEAM POTENTIAL/RUNTIME VTK OUTPUT
VTK_OUTPUT_BEAM_POTENTIAL           yes
INTERVAL_STEPS                      1
EVERY_ITERATION                     no
FORCES                              yes
MOMENTS                             yes
WRITE_FORCE_MOMENT_PER_ELEMENTPAIR  no // maybe change to yes

Finally, I need to add the interaction potential to each fiber as follows (each potential law needs to be added to each line, e.g., law 1 to line 1 and 2 and law 2 to line 1 and 2):

--------------------------------DESIGN LINE BEAM POTENTIAL CHARGE CONDITIONS
DLINE                           2
E 1 - POTLAW 1 VAL 13.123 FUNCT none
E 1 - POTLAW 2 VAL 13.123 FUNCT none
E 2 - POTLAW 1 VAL 13.123 FUNCT none
E 2 - POTLAW 2 VAL 13.123 FUNCT none

Regarding boundary conditions I want to fix all nodes at certain planes in all directions, e.g. fix all nodes at z=0.0 and final geometry length z=100.0.

This leads to the following questions where I am unsure on how to proceed:

1. Would it be the best option to simply add the overall BEAM POTENTIAL settings via input_file.add("""BEAM POTENTIAL ...""")? Or is it useful to somehow implement this as a fixed function into MeshPy?

2. Would it be ok to add a new runtime output for potential interactions to header_functions.py ->set_runtime_output() ?

3. Regarding the DESIGN LINE BEAM POTENTIAL CHARGE CONDITIONS; as a preferred way I would just add a helix_set.add_potential_charge_condition(potlaw=1, value=13.123, function=none). Is it a good/doable thing within PyMesh? It would be nice if this could directly happen within the geometry creation as each geometry set already has all necessary information on its corresponding DLINE. Currently I am unsure on how to proceed within MeshPy, would it make sense to create a new boundary condition?

4. Regarding the boundary conditions; Is it possible to select specific nodes at a plane in space (e.g., z=100.0 or z=100.0+-0.01) and add a boundary condition to them? If not possible, would it make sense to include this option in MeshPy?

Thanks in advance for your help!

[isteinbrecher]

1. Would it be the best option to simply add the overall BEAM POTENTIAL settings via input_file.add("""BEAM POTENTIAL ...""")? Or is it useful to somehow implement this as a fixed function into MeshPy?
2. Would it be ok to add a new runtime output for potential interactions to header_functions.py ->set_runtime_output()

The utility functions for creating header sections / input file parameters in 4C have a historic background. They were introduced to ease the creation of input files by wrapping common parameters. On one hand, this has the advantages of the MeshPy scripts being relative robust regarding changes in 4C. On the other hand, this requires a constant update of MeshPy to keep the parameter naming in 4C up to date. I personally only use these functions on rare occasions any more. For me it turned out to be simpler to just add all input file parameters directly as strings to the input file. However, I am not creating input scripts on a regular basis as I used to, so maybe I am biased there.

Long story short, do as you seem best fit. In the near future I can envision a better encapsulation of the 4C specific functionality in MeshPy, similar to the ABAQUS functionaliy stored in meshpy/abaqus. So feel free to add a meshpy/four_c/beam_potential.py file that contains the beam potential specific header functions (would also be a good place to store beam potential boundary condition functions).

3. Regarding the DESIGN LINE BEAM POTENTIAL CHARGE CONDITIONS

You have two options here:

1. You can add your boundary condition to the boundary condition definitions of MeshPy.
2. You can simply state the section name as a string when defining the boundary conditions

   for i_pot_law in range(2):
       input_file.add(
           BoundaryCondition(
               helix_set["line"],
               f"POTLAW {i_pot_law+1} VAL 13.123 FUNCT none",
               bc_type="DESIGN LINE BEAM POTENTIAL CHARGE CONDITIONS",
           )
       )

4. Regarding the boundary conditions; Is it possible to select specific nodes at a plane in space (e.g., z=100.0 or z=100.0+-0.01) and add a boundary condition to them? If not possible, would it make sense to include this option in MeshPy?

Yes, that is easily possible. Try the following code snippet to fix all points at z=100:

    def find_node_on_plane(node, *, normal=None, origin_distance=None, tol=mpy.eps_pos):
        """Return the nodes lie on the plane defined by normal and origin_distance"""
        projection = np.dot(node.coordinates, normal) / np.linalg.norm(normal)
        return np.abs(projection - origin_distance) < tol

    nodes_at_100 = input_file.get_nodes_by_function(
        find_node_on_plane, normal=[0, 0, 1], origin_distance=100.0, tol=0.1
    )
    input_file.add(
        BoundaryCondition(
            GeometrySet(nodes_at_100),
            "NUMDOF 9 ONOFF 1 1 1 1 1 1 0 0 0",
            bc_type=mpy.bc.dirichlet,
        )
    )

If this turns out to be used in many cases, we can add a file containing "node finding functions". To store commonly used functions like this one.

[isteinbrecher]

Regarding the script you posted, looks good overall, two comments:

1. You almost never have to add geometry sets explicitly to the mesh like you do in input_file.add(helix_set["line"]). The sets in 4C only make sense when they are used within a boundary condition, and in this case MeshPy adds the set automatically. This makes however sense for testing, when you want to check the created sets without the need to define a boundary condition.
2. I missed one important detail when we added the helix function. The "add geometry" functions should not alter the existing geometry in the given mesh. Basically this means that you don't have to create the tmp_mesh, this should be done internally. Can you please add a fixup for this? (If you are unsure, you can have a look how this in done for example in create_beam_mesh_honeycomb)

[davidrudlstorfer]

Thanks a lot @isteinbrecher for your elaborate and helpful input! Without that it would not have been possible to get such a good solution! I have pushed a first implementation to include potential interactions within MeshPy for 4C in #74