[CIG-SHORT] Principle stress

[Brad Aagaard]

On 12/04/2013 10:11 PM, Ma, Xiao wrote:
> Hi,
> I am still trying to output the priciple stress in pylith, now when I am using the pylith_genxdmf --file=stress.h5
> I get this error:
> _____________________________________
> assertion "0" failed: file "topology/FieldBase.cc", line 87, function: static py
> lith::topology::FieldBase::VectorFieldEnum pylith::topology::FieldBase::parseVec
> torFieldString(const char*)
> Aborted (core dumped)

You are getting this error because the pylith_genxdmf script isn't able 
to parse the vector_field_attribute of a field (probably pstress) into a 
known value. Use the h5dump command on the HDF5 file to verify that the 
vector_field_type attribute for the pstress field is set to 'other' and 
that the field data values are correct.

> _______________________________________
> What I am doing is like the following:
> ________________________________________
> import numpy
> import h5py
> import matplotlib.pyplot as plt
> # Load in change in tractions from coseismic simulation
> h5 = h5py.File("output/stress.h5", 'r', driver="sec2")
>
> vertices = h5['geometry/vertices'][:]
> #tractions_change = h5['vertex_fields/traction_change'][0,:,:]
> stress=h5['cell_fields/stress'][:]
> h5.close()
> pstress=(stress[:,:,0]+stress[:,:,1])/2+numpy.sqrt(((stress[:,:,0]-stress[:,:,1])/2)**2+(stress[:,:,2])**2)
> h5 = h5py.File("output/stress.h5", 'a', driver="sec2")
> field=h5.create_dataset('cell_fields/pstress',data=pstress[:,:])
> field.attrs['vector_field_type']='other'
> h5.close()
> __________________________________________
> I can see that in the .h5 file the pstress is already add into the cell fields.
> Also for example the principle stress pstress  dimension is 500,10000, which is good , because the 500 time steps, and 10000 elments, however the vertices number is 10302, because there is a fault interface in the model, there are extra nodes.
> What I want to do is that I want plot contours of pricple stress all over the mesh, using python, but the vertrices have repeated coordinates , because the fault interface.
> How can I implement that ?

The principal stresses are a field over the cells, so you should be able 
to visualize it in ParaView once the Xdmf/HDF5 files are setup properly. 
The changes to the topology as a result of the insertion of the fault 
don't matter matter as long as any post-processing you do maintains 
fields consistent with the topology in the output files.

To make contours, you need to convert the cell data to point data in 
ParaView (see the Filters) and then generate the contours using isosurfaces.

If you have further questions about ParaView, please see the ParaView 
documentation, http://www.paraview.org/paraview/help/documentation.html. 
If you are using another visualization package, see its documentation.

Regards,
Brad