# Scripting with python

All files in this tutorial can be downloaded here: pythonScript.zip and the jupyter-notebook

This tutorial will cover how scripting with python/jupyter-notebook is done. Therefore, we automatically generate a mesh with variable mesh discretization, as well as running a simulation with variable input parameters. Our Testexample will be a simple cube, where we apply a sinoidal deformation on top.

After the automatic simulation we will do a short postprocessing and visualize the results.

## Lets go through the script

First, lets import all the necessary libaries. Dont forget to specify the path where those libaries are located. In our case they are in the same folder.

Click to see code-snippet
%load_ext autoreload

import sys
sys.path.append("./")

import matplotlib.pyplot as plt

import hdf5_tools
import numpy as np
import cfsResultClass as c


Now in this section, we first define our parameters for the simulation and mesh and generate a mesh-journal file and simulation-input xml based on templates. We read in the templates and replacing specific strings with our parameters.

Click to see code-snippet
steps = 10
deltaT = 1
amplitude = 0.001
interval = 3

#Write the simulation input out of a template file

#creating a unique name for this parameter combination
simulationName = f"sim_Steps{steps}_dT{deltaT}_Amp{amplitude}_Int{interval}"

sim_tmp = open("sim_tmp.xml", "r")
sim = open(simulationName + ".xml", "w")
for line in sim_tmp.readlines():
if "STEPS" in line:
line = line.replace("STEPS", str(steps))
if "DELTAT" in line:
line = line.replace("DELTAT",str(deltaT))
if "AMPLITUDE" in line:
line = line.replace("AMPLITUDE", str(amplitude))

sim.writelines(line)
sim_tmp.close()
sim.close()

#Write the mesh-journal file
jou_tmp = open("UnitCube_tmp.jou", "r")
journalName = f"UnitCube_f{interval}.jou"
jou = open(journalName, "w")
for line in jou_tmp.readlines():
if "INTERVAL" in line:
line = line.replace("INTERVAL", str(interval))
jou.writelines(line)
jou_tmp.close()
jou.close()



One could aswell use a xml-libary to utilize the full potential of the xml-input. This could be realized with the following code snippet:

Click to see code-snippet
import xml.etree.ElementTree as ET
ns = {'cfs': "http://www.cfs++.org/simulation"}
#ET.register_namespace('', "http://www.cfs++.org/simulation")

tree = ET.parse("sim_tmp.xml")
root = tree.getroot()

#Using XPath to reduce code
for numSteps in root.findall(".//cfs:numSteps",ns):
numSteps.text = str(steps)

for dT in root.findall(".//cfs:deltaT",ns):
dT.text = str(deltaT)

for displacement in root.findall(".//*[@name='S_T']/cfs:comp",ns):
displacement.set("value",f"{amplitude}*sin(t)")

# Without XPATH the replacement of the amplitude would look like this:
# for load in root.iter('{http://www.cfs++.org/simulation}displacement'):
#     if load.get("name") == "S_T":
#         for child in list(load):
#             if child.get("dof")=="z":
#                 child.set("value",f"{amplitude}*sin(t)")

tree.write("newXML.xml")


After generating the journal file and the simulation-xml-input we have to create the mesh and run the simulation. This is done in the next section. Here we use jupyter-notebook cell magic with %%capture output to supress the output of the cell and the ipython command ! to run commando commands. If you want to see the terminal command, comment the first line out.

If you are unsure about the local path of your program, just execute the command type Programname in the terminal. This command will return the path were the executable is stored.

After simulation, we delete all created files to not pollute this folder. This is especially important if storage-place is critical.

Click to see code-snippet
%%capture output
# To supress all output of this cell, MUST be ontop of the cell

#Create mesh, insert your local path
!/opt/programs/Coreform-Cubit-2021.11/bin/coreform_cubit -batch -nojournal -nographics $journalName #Simulate, insert your local path to cfs !/home/alex/Devel/CFS_BIN/EclipseDebugBuild/bin/cfs$simulationName

#Delete simulation input files, to not pullute folder
!rm *.info.xml
!rm $simulationName* !rm$journalName*


After running the simulation we read in the cfs-result with as a cfsResult. This is a class whichs automatically reads in all results. However its currently limited to only read 1 sequence step and 1 region at a time, otherwise an error occurs. But this functionallity could be easly added (by you?).

Click to see code-snippet
# Read in the result with the help of cfsResultClass. (Currently this is only limited for 1 sequence step and one region!)
result = c.cfsResult(f"./results_hdf5/{simulationName}.cfs",multistep =1, region=None)


After reading in the result we do a bit of postprocessing and plot the displacement of one node on the bottom and top of the unitcube. We also display the vonMises-stress of an element which is placed on top of the cube.

Click to see code-snippet
#Read in Values
#Displacement in Z
#u_z[Nt, Node] = value
u_z = result.mechDisplacement.Values.z
#coord_u[Node,xyz]
coord_u = result.mechDisplacement.Values.coord
#Get all indices from nodes which are placed on top and bottom
idx_u_zTop = np.where(coord_u[:,2] == coord_u[:,2].max())[0]
idx_u_zBot = np.where(coord_u[:,2] == coord_u[:,2].min())[0]

#Get the timesteps
t = hdf5_tools.get_step_values(f"./results_hdf5/{simulationName}.cfs")[0]

#VonMises Stress
#vonMises[Nt, Element] = value
vonMises = result.vonMisesStress.Values.val
coord_vonMises = result.vonMisesStress.Values.coord
#Get all indices from elements which are placed on top
idx_vonMisesTop = np.where(coord_vonMises[:,2] == coord_vonMises[:,2].max())[0]

# Plot results
fig, ax = plt.subplots(2, figsize = (5,2*3))

#plot mechDisplacement

ax[0].plot(t, u_z[:,idx_u_zTop[0]], label = "One node placed on top")
ax[0].plot(t, u_z[:,idx_u_zBot[0]], label = "One node on bottom")
ax[0].set_ylabel("Displacement in m")
ax[0].set_xlabel("Timesteps")
ax[0].set_title("Displacement over timesteps")
ax[0].legend()

ax[1].plot(t, vonMises[:,idx_vonMisesTop[0]], label = "One element placed on top")
ax[1].set_ylabel("vonMises-Stress in N/m²")
ax[1].set_xlabel("Time in s")
ax[1].set_title("VonMises-Stress over time")
ax[1].legend()

fig.suptitle("Results", fontsize = 16)

#Makes plots prettier (eg. labels dont stand in other plots)
plt.tight_layout()
fig.savefig("resultPlot.png", dpi = 100)


## Further suggestions

• Increase the time discretiation by changing deltaT, to get a smoother result
• Create a for-loop for different timesteps
• Save those results in a compact dictinary
• Do a convergence study:
• Compare the difference of the solution at different mesh-discretization
• for different time step sizes