Add potential interactions in 4C

[davidrudlstorfer]

In this PR I added the following things:

• A new top level directory FourC (I learned that Python modules are not allowed to start with a number, e.g., 4C) which now includes the beam_potential.py file for potential interactions.

• Therein, I've added a new class BeamPotential which helps with the header and runtime output definition as well as the boundary conditions for potential interactions

• I also added your proposed code snippet to the utilities of MeshPy, as this is a very useful tool.

The final structure is probably upon discussion, therefore, I've waited to add a test (as this only tests for simple header changes this is not yet important)

Finally I now use this code:

Click me

import numpy as np

from meshpy import (
  InputFile,
  GeometrySet,
  Function,
  MaterialReissner,
  Beam3rHerm2Line3,
  BoundaryCondition,
)

from meshpy.conf import mpy
from meshpy.FourC.beam_potential import BeamPotential
from meshpy.mesh_creation_functions.beam_basic_geometry import create_beam_mesh_helix
from meshpy.utility import find_node_on_plane


def main():
  """Minimal example of multiple helical beams with potential interactions."""

  ### define input file
  input_file = InputFile()

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

  ### determine quantities of helice
  # 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))
      )

  ### Define functions for line charge density
  fun_1 = Function("TESTFUNCTION 123")
  fun_1.n_global = 1
  fun_2 = Function("TESTFUNCTION 234")
  fun_2.n_global = 2

  ### Define beam potential
  beampotential = BeamPotential(
      input_file,
      pot_law_prefactor=[-1.0e-3, 12.45e-8],
      pot_law_exponent=[6.0, 12.0],
      pot_law_line_charge_density=[1.0, 2.0],
      pot_law_line_charge_density_funcs=[fun_1, fun_2],
  )

  ### Add header and runtime output to .dat file

  beampotential.add_header(
      cutoff_radius=10.0,
      evaluation_strategy="SingleLengthSpecific_SmallSepApprox_Simple",
      potential_reduction_length=15.0,
  )
  beampotential.add_runtime_output()

  ### Create helices and add potential charge condition to each
  for radius, twist_angle in zip(radii, twist_angles):

      helix_set = create_beam_mesh_helix(
          input_file,
          Beam3rHerm2Line3,
          mat,
          [0, 0, 1],  # axis
          [0, 0, 0],  # axis_point
          [radius, 0, radius],  # start_point
          twist_angle,  # twist_angle
          height_helix=80,
          n_el=20,
          warning_straight_line=False,
      )

      beampotential.add_potential_charge_condition(geometry_set=helix_set["line"])

  ### Add boundary condition to bottom and top node
  input_file.add(
      BoundaryCondition(
          GeometrySet(
              input_file.get_nodes_by_function(
                  find_node_on_plane, normal=[0, 0, 1], origin_distance=0.0, tol=0.1
              )
          ),
          "NUMDOF 9 ONOFF 1 1 1 1 1 1 0 0 0",
          bc_type=mpy.bc.dirichlet,
      )
  )

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

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


if __name__ == "__main__":

  main()

to create the following perfectly fine input file (.txt instead of .dat due to Github not supporting .dat files)
test_potential_interaction.txt

I would be more than happy to get your feedback regarding these changes 😊

---

Closes #72

[davidrudlstorfer]

One thing I am unsure about is how to add a global (or get) a global ID (number) for a defined function?

I specifically set it with

fun_1.n_global = 1

[isteinbrecher]

Thanks @davidrudlstorfer for the contribution. I left some comments. Also can you please add a test case to cover your changes? One (or how many you need) for the beam potential interactions and one where we check the nodes_on_plane function?

[isteinbrecher]

@davidrudlstorfer in general the script looks good to me, one comment:

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

In testing we always turn the header off, but in practice I would keep it on, as it adds some information to your input file that can later be used to help to identify the state of the script when the file was created.

[davidrudlstorfer]

@isteinbrecher thanks for your review! I've added all of your recommendations to latest changes. Also thanks for the help with the function numbering within the .dat file.

Regarding the tests:

I took the liberty and renamed testing_utility.py into utilities.py and created a testing_utility.py for the newly added node testing function. This was done to adhere with the current file naming scheme of the test files.

I've also added a test to check an exemplary beam potential interaction input file for 4C.

[isteinbrecher]

@davidrudlstorfer everything looks good to me, if you think this is ready, we can merge this.

[davidrudlstorfer]

Thanks for the review! From my side it would be ready to merge 🚀

[isteinbrecher]

Can you click on the merge button or is this for member of the group imcs-compsim only?

[davidrudlstorfer]

Unfortunately I am not authorized to push the big button

[isteinbrecher]

I invited you to the project with write acces, does it work now?

[davidrudlstorfer]

@isteinbrecher thanks for the invitation, unfortunately it still does not work

Is the rule active within this repository that contributors cannot merge if they created the PR?

[isteinbrecher]

How about now?

[davidrudlstorfer]

Yes - worked like a charm. Thanks for your help and review regarding this PR! This concludes the first steps for me regarding MeshPy. I am now going to include it within my own repo for the creation of input files. I am going to contribute further to MeshPy depending on different upcoming needs 😄

[isteinbrecher]

Thanks for your contributions and looking forward to the next ones! I might come up to you for some small improvements cleanup in the future 😉