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 ...
%load_ext autoreload
%autoreload 2

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

Read history results¶

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
w_fe = np.abs(np.loadtxt(fname[0],usecols=[2]))
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()

Reading element results¶

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

Additional Tutorials¶

Please have a look at this scripting tutorial with python
Please have a look at this scripting tutorial with python