# Python Postprocessing of CFS++ Results

As an example we treat a single 2D peam problem.

## HDF5 Data Access

We use a simple beam example to demonstrate the postprocessing of data in python. We can read the maximum deflection from the hdf-database (*.cfs file). We use the h5py library and some tools for accessing CFS++ files.

### h5py basics

The library provides a dictionary-like interface to access hdf5 files. For more details consult the documentation: http://docs.h5py.org/en/latest/

# import the h5py library
import h5py

# load a hdf5-File
h5 = h5py.File('results_hdf5/myBeam.cfs','r')

# list groups of hdf5 file
list(h5.keys())
print([ var for var in h5 ])

# list sub-Groups in group 'Mesh'
for name in h5['Mesh'] :
print(name)

# show the complete strucure of a hdf5 file
def print_h5(h5f,lvl=0) :
for name in h5f:
print('  '*lvl+name)
if type(h5f[name])==h5py.Dataset :
pass
#print(h5f[name].shape) # print shape
elif type(h5f[name])==h5py.Group : # print all sub-groups
print_h5(h5f[name],lvl=lvl+1)
else :
print('It is a ',type(h5f[name]))

print()
print("Element 'Results/Mesh/MultiStep_1'")
print_h5(h5['Results/Mesh/MultiStep_1'])

h5.close()


### Reading from CFS++ result files

There are some python functions available for CFS++ result files in the library hdf5_tools.py. Use them to conveniantly access data from CFS++. For more information on the data structure see https://cfs-doc.mdmt.tuwien.ac.at/mediawiki/index.php/HDF5_Data_Structure.

# if you use your own functions in a library and change them while running the notebook,
# make sure to auto-reloead them ...

from sys import path
#path.append('/home/ftoth/cfs/CFS/share/python')
path.append('../../share/') # adapt this path so that it poits to the location of hdf5_tools.py

from hdf5_tools import *

try:
h5f = h5py.File('results_hdf5/myBeam.cfs','r')
X = get_coordinates(h5f) # all node coorinates
U = get_result(h5f,'mechDisplacement',region='S_beam',multistep=1) # all displacements
h5f.close()
except:
# the h5 file should be closed if the data access does not work
h5f.close()
raise

# find the maximum in y-direction (=minimum value)
U[:,1].min()


### Reload of changed Data Files

In order to re-load a changed data file, e.g. to plot the new result after some input change, you must * close the data file * open it again

In the iPython notebook this only works properly if you close the file object in the same cell as you opened it before. If if does not work one can restart the kernel. A few possiblities which make sure the file is closed are given in cells above and below.

# use a string denoting the filename (will open & close) -> possible IO-overhead
X = get_coordinates('results_hdf5/myBeam.cfs')
X[:3,:]

# using with to make sure the file is closed ...
with h5py.File('results_hdf5/myBeam.cfs','r') as h5f:
X = get_coordinates(h5f) # all node coorinates
U = get_result(h5f,'mechDisplacement',region='S_beam',multistep=1) # all displacements
U[:5,:]


## Python scipting Examples

### Compute the analytical solution

We compute the static deflection of a clmaned beam subject to a point load at the end.

# define the beam dimensions
import numpy as np
l = 0.5 # length in m
h = 0.05 # thickness in m
b = 1.0 # 2D-dimension
E = 70e+9 # Youngs Modulus in Pa
tau = 1.0e+6 # traction in Pa

# analytical solution
J = b*h**3/12 # areal moment of inertia
P = tau*b*h # end load
w_max = P*l**3/(3*E*J) # maximum deflection
print('analytical: %.3f mm'%(w_max*1e+3))


We read the deflection at the end of the beam and compare them to the analytical result computed above.

# load the FE result
from glob import glob
fname = glob('history/myBeam-mechDisplacement-node-*-end.hist') # node number might depend on the mesh
print('FE        : %.3f mm'%(w_fe*1e+3))
print("Error     : %.2f %%"%(((abs(w_fe)-w_max)/w_max)*100))


### Graphs

We plot the deflection of the beam along its axis.

# show plots inline
%matplotlib inline

# load matplotlib
import matplotlib as mpl
from matplotlib import pyplot as plt
# set sensible defaults
mpl.rc('figure',figsize=[2.7,1.6]) # figure size in inch
mpl.rc('figure',dpi=200) # inline dpi (=display size in browser)
mpl.rc('font',size=8.0)
#mpl.rcParams['font.sans-serif'] = 'Helvetica'
mpl.rc('lines',linewidth=0.7)
mpl.rc('axes',grid=True,linewidth=1.0,axisbelow=True,unicode_minus=False)
mpl.rc('grid',linewidth=0.3,linestyle=':')
mpl.rc('legend',fontsize='medium',framealpha=1.0,numpoints=1)
mpl.rc('svg',fonttype='none')
mpl.rc('savefig',dpi=600)

It = get_subregion_idx('results_hdf5/myBeam.cfs','S_beam','L_top') # indices of nodes in 'L_top' in 'S_beam'

fig,ax=plt.subplots()
ax.plot(X[It,0],U[It,1],label='FE (top)')
ax.set_xlabel('x coordinate in m')
ax.set_ylabel('y-dispacement in m')

# compute and plot the analytic solution
x = np.linspace(0,l,100)
w = P/(6*E*J)*(3*l*x**2 - x**3)
ax.plot(x-l/2,-w,'k--',lw=0.3,label='analytic') # shift in x and correct sign

ax.legend()
fig.tight_layout()


We read the stress tensor saved for the centroid of each element in the *.cfs file. The determine the maximum value, and plot the component $\sigma_{xx}$ over the beam thickness ($y$) at the leftmost row of elements.

sig_max = P*l/J*h/2. # maximum stress
print('analytical: %.2f MPa'%(sig_max*1e-6))

# read the stress tensor (2D) for all element centroids
h5f = 'results_hdf5/myBeam.cfs'
S = get_result(h5f,'mechStress')
C = get_centroids(h5f,region='S_beam') # get centroid coordinates

S.shape,C.shape

# get the maximum s_xx
S[:,0].max()

# compute the von-Mises stress, result columns are: s_xx, s_yy, s_xy
s_v = np.sqrt( S[:,0]**2 + S[:,1]**2 - S[:,0]*S[:,1] + 3*S[:,2]**2 )
s_v.max()

# plot s_xx for the leftmost elements

# search for elements on the left (close to minimal x-coordinate)
Ie_left = np.argwhere(np.abs(C[:,0] - C[:,0].min())<1e-6).ravel()

# get their y-coorninate
y = C[Ie_left,1]
s_xx = S[Ie_left,0]

fig,ax=plt.subplots()
ax.plot(s_xx*1e-6,y*1e+3,'.-')
ax.set_xlabel('$\sigma_{xx}$ in MPa')
ax.set_ylabel('y-coordinate in mm')
ax.set_ylim(-h/2*1e+3,h/2*1e+3)
ax.set_xlim(-sig_max*1e-6,sig_max*1e-6)


### Step values

Step values are the time steps of a transient analysis, the frequency points of an harmonic analysis, or the natural frequencies of an eigenfrequency step.

stepvals = get_step_values('results_hdf5/myBeam.cfs')
stepvals[1]


### Time Signals

Extract the results from all time steps of a transient analysis. For single points it might be beneficial to write the time history into *.hist files or as history data to the *.cfs file.

Ut = get_result('results_hdf5/myBeam.cfs','mechDisplacement',region='S_beam',step='all',multistep=3) # all displacements
Ut.shape

# find node at the right-top corner
Irt = np.argmin(np.linalg.norm(X - X.max(axis=0),axis=1))
X[Irt,:]

# plot time signal
t = stepvals[2]
fig,ax=plt.subplots()
ax.plot(t,Ut[:,Irt,1])
ax.set_xlabel('Time in s')
ax.set_ylabel('y-displacement in m')