added principal stress computation

Author: akaszynski

.gitignore

@@ -11,6 +11,7 @@
 ######################
 .fuse_hidden*
 *~
+*swp
 
 # Pip generated folders #
 #########################
@@ -24,4 +25,3 @@ pyansys/Interface.py
 
 # Testing
 Testing/
-

doc/index.rst

@@ -30,7 +30,7 @@ Dependencies:
 - ``vtkInterface``
 - ``vtk`` (optional)
 
-Minimum requirements are numpy to extract results from a results file. To convert the raw data to a VTK unstructured grid, VTK 5.0 or greater must be installed with Python bindings.
+Minimum requirements are numpy to extract results from a results file. To convert the raw data to a VTK unstructured grid, VTK 7.0 or greater must be installed with Python bindings.
 
 See `Install VTK <http://vtkinterface.readthedocs.io/en/latest/installation.html>`_ for details instructions for installing VTK.
 

doc/loading_results.rst

@@ -5,6 +5,7 @@ The ANSYS result file is a FORTRAN formatted binary file containing the results
     - Nodal DOF results from a static analysis or modal analysis.
     - Nodal DOF results from a cyclic static or modal analysis.
     - Nodal averaged component stresses (i.e. x, y, z, xy, xz, yz)
+    - Nodal principal stresses (i.e. S1, S2, S3, SEQV, SINT)
 
 We're working on adding additional plotting and retrieval functions to the code If you would like us to add an additional result type to be loaded, please open an issue in `GitHub <https://github.com/akaszynski/pyansys>`_  and include result file for the result type you wish to load.
 
@@ -84,7 +85,7 @@ The DOF solution for an analysis for each node in the analysis can be obtained u
     # normalized displacement can be plotted by excluding the direction string
     result.PlotNodalResult(0, label='Normalized')
 
-Stress can be obtained as well using the below code.  The nodal stress is computed in the same manner as ANSYS by averaging the stress evaluated at that node for all attached elements.  For now, only component stresses can be displayed.
+Stress can be obtained as well using the below code.  The nodal stress is computed in the same manner as ANSYS by averaging the stress evaluated at that node for all attached elements.
 
 .. code:: python
     
@@ -95,6 +96,10 @@ Stress can be obtained as well using the below code.  The nodal stress is comput
     # Display node averaged stress in x direction for result 6
     result.PlotNodalStress(5, 'Sx')
 
+    # Compute principal nodal stresses and plot SEQV for result 1
+    pstress = result.PrincipalNodalStress(0)
+    result.PlotPrincipalNodalStress(0, 'SEQV')
+
 
 Results from a Cyclic Analysis
 ------------------------------

pyansys/_version.py

@@ -1,5 +1,5 @@
 # major, minor, patch
-version_info = 0, 20, 0
+version_info = 0, 21, 0
 
 
 # Nice string for the version

pyansys/binary_reader.py

@@ -846,18 +846,18 @@ def StoreGeometry(self):
         self.edge_node_num_idx = np.unique(self.edge_idx)
 
         # catch the disassociated node bug
-#        try:
-#            self.edge_node_num = self.geometry['nnum'][self.edge_node_num_idx]
-#        except:
-#            logging.warning('unable to generate edge_node_num')
+        try:
+            self.edge_node_num = self.geometry['nnum'][self.edge_node_num_idx]
+        except:
+            logging.warning('unable to generate edge_node_num')
 
     def NodalStress(self, rnum):
         """
         Retrives the component stresses for each node in the solution.
 
         The order of the results corresponds to the sorted node numbering.
 
-        This algorthim, like ANSYS, computes the nodal stress by averaging the
+        This algorithm, like ANSYS, computes the nodal stress by averaging the
         stress for each element at each node.  Due to the discontinunities
         across elements, stresses will vary based on the element they are
         evaluated from.
@@ -869,12 +869,12 @@ def NodalStress(self, rnum):
 
         Returns
         -------
-        stress : array
+        stress : numpy.ndarray
             Stresses at Sx Sy Sz Sxy Syz Sxz averaged at each corner node.
-            For the corresponding node numbers, see "edge_node_num"
+            For the corresponding node numbers, see "result.edge_node_num"
+            where result is the result object.
 
         """
-
         # Get the header information from the header dictionary
         endian = self.resultheader['endian']
         rpointers = self.resultheader['rpointers']
@@ -960,6 +960,90 @@ def NodalStress(self, rnum):
 
         return s_node
 
+    def PrincipalNodalStress(self, rnum):
+        """
+        Computes the principal component stresses for each node in the
+        solution.
+
+        The order of the results corresponds to the sorted node numbering.
+        See result.edge_node_num
+
+        Parameters
+        ----------
+        rnum : interger
+            Result set to load using zero based indexing.
+
+        Returns
+        -------
+        pstress : numpy.ndarray
+            Principal stresses, stress intensity, and equivalant stress.
+            [sigma1, sigma2, sigma3, sint, seqv]
+
+        Notes
+        -----
+        ANSYS equivalant of:
+        PRNSOL, S, PRIN
+
+        which returns:
+        S1, S2, S3 principal stresses, SINT stress intensity, and SEQV
+        equivalent stress.
+
+        """
+        # get component stress
+        stress = self.NodalStress(rnum)
+
+        # compute principle stress
+        if stress.dtype != np.float32:
+            stress = stress.astype(np.float32)
+        return _rstHelper.ComputePrincipalStress(stress)
+
+    def PlotPrincipalNodalStress(self, rnum, stype):
+        """
+        Plot the principal stress at each node in the solution.
+
+        Parameters
+        ----------
+        rnum : interger
+            Result set using zero based indexing.
+
+        stype : string
+            Stress type to plot.  S1, S2, S3 principal stresses, SINT stress
+            intensity, and SEQV equivalent stress.
+
+            Stress type must be a string from the following list:
+
+            ['S1', 'S2', 'S3', 'SINT', 'SEQV']
+
+        Returns
+        -------
+        cpos : list
+            VTK camera position.
+
+        """
+        stress_types = ['S1', 'S2', 'S3', 'SINT', 'SEQV']
+        if stype not in stress_types:
+            raise Exception('Stress type not in \n' + str(stress_types))
+
+        sidx = stress_types.index(stype)
+
+        # create a zero array for all nodes
+        stress = np.zeros(self.resultheader['nnod'])
+
+        # Populate with nodal stress at edge nodes
+        pstress = self.PrincipalNodalStress(rnum)
+        stress[self.edge_node_num_idx] = pstress[:, sidx]
+        stitle = 'Nodal Stress\n{:s}'.format(stype)
+
+        # Generate plot
+        plobj = vtkInterface.PlotClass()
+        plobj.AddMesh(self.grid, scalars=stress, stitle=stitle,
+                      flipscalars=True, interpolatebeforemap=True)
+
+        text = 'Result {:d} at {:f}'.format(rnum + 1, self.tvalues[rnum])
+        plobj.AddText(text)
+
+        return plobj.Plot()  # store camera position
+
     def PlotNodalStress(self, rnum, stype):
         """
         Plots the stresses at each node in the solution.
@@ -1002,7 +1086,6 @@ def PlotNodalStress(self, rnum, stype):
         plobj.AddText(text)
 
         cpos = plobj.Plot()  # store camera position
-        del plobj
 
         return cpos
 

pyansys/cython/_rstHelper.pyx

@@ -6,6 +6,7 @@ import ctypes
 import numpy as np
 cimport numpy as np
 
+from libc.math cimport sqrt, fabs
 from libc.stdio cimport fopen, FILE, fclose, fread, fseek
 from libc.stdio cimport SEEK_CUR, ftell, SEEK_SET
 from libc.string cimport memcpy
@@ -425,4 +426,88 @@ def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn):
         const_sort[i] = const[s_neqv_dof[i]]
 
     return np.asarray(nref_sort), np.asarray(dref_sort), np.asarray(index_arr), \
-           np.asarray(const_sort), np.asarray(ndof)
+           np.asarray(const_sort), np.asarray(ndof)
+
+
+def ComputePrincipalStress(float [:, ::1] stress):
+    """
+    Returns the principal stresses based on component stresses.
+
+    Parameters
+    ----------
+    stress : numpy.ndarray
+        Stresses at Sx Sy Sz Sxy Syz Sxz averaged at each corner node.
+
+    Returns
+    -------
+    pstress : numpy.ndarray
+        Principal stresses, stress intensity, and equivalant stress.
+        [sigma1, sigma2, sigma3, sint, seqv]
+
+    Notes
+    -----
+    ANSYS equivalant of:
+    PRNSOL, S, PRIN
+
+    Which returns:
+    S1, S2, S3 principal stresses, SINT stress intensity, and SEQV
+    equivalent stress.
+    """
+    # reshape the stress array into 3x3 stress tensor arrays
+    cdef int nnode = stress.shape[0]
+    cdef float [:, :, ::1] stress_tensor = np.empty((nnode, 3, 3), np.float32)
+    cdef float s_xx, x_yy, s_zz, s_xy, s_yz, s_xz
+
+    for i in range(nnode):
+        s_xx = stress[i, 0]
+        s_yy = stress[i, 1]
+        s_zz = stress[i, 2]
+        s_xy = stress[i, 3]
+        s_yz = stress[i, 4]
+        s_xz = stress[i, 5]
+
+        # populate stress tensor
+        stress_tensor[i, 0, 0] = s_xx
+        stress_tensor[i, 0, 1] = s_xy
+        stress_tensor[i, 0, 2] = s_xz
+        stress_tensor[i, 1, 0] = s_xy
+        stress_tensor[i, 1, 1] = s_yy
+        stress_tensor[i, 1, 2] = s_yz
+        stress_tensor[i, 2, 0] = s_xz
+        stress_tensor[i, 2, 1] = s_yz
+        stress_tensor[i, 2, 2] = s_zz
+
+    # compute principle stresses
+    w, v = np.linalg.eig(np.asarray(stress_tensor))
+    w[:, ::-1].sort(1)
+
+    temp = np.empty((nnode, 5), np.float32)
+    temp[:, :3] = w
+
+    cdef float [:, ::1] pstress = temp
+    cdef float p1, p2, p3, c1, c2, c3
+
+    # compute stress intensity and von mises (equivalent) stress
+    for i in range(nnode):
+        p1 = pstress[i, 0]
+        p2 = pstress[i, 1]
+        p3 = pstress[i, 2]
+
+        c1 = fabs(p1 - p2)
+        c2 = fabs(p2 - p3)
+        c3 = fabs(p3 - p1)
+
+        if c1 > c2:
+            if c1 > c3:
+                pstress[i, 3] = c1
+            else:
+                pstress[i, 3] = c3
+        else:
+            if c2 > c3:
+                pstress[i, 3] = c2
+            else:
+                pstress[i, 3] = c3
+
+        pstress[i, 4] = sqrt(0.5*(c1**2 + c2**2 + c3**2))
+
+    return np.asarray(pstress)

setup.py

@@ -6,17 +6,17 @@
 from io import open as io_open
 
 from setuptools import setup, Extension
-from setuptools.command.build_ext import build_ext
-
-# from setuptools.command.build_ext import build_ext as _build_ext
-# # Create a build class that includes numpy directory
-# class build_ext(_build_ext):
-#     def finalize_options(self):
-#         _build_ext.finalize_options(self)
-#         # Prevent numpy from thinking it is still in its setup process:
-#         __builtins__.__NUMPY_SETUP__ = False
-#         import numpy
-#         self.include_dirs.append(numpy.get_include())
+#from setuptools.command.build_ext import build_ext
+
+from setuptools.command.build_ext import build_ext as _build_ext
+# Create a build class that includes numpy directory
+class build_ext(_build_ext):
+    def finalize_options(self):
+        _build_ext.finalize_options(self)
+        # Prevent numpy from thinking it is still in its setup process:
+        __builtins__.__NUMPY_SETUP__ = False
+        import numpy
+        self.include_dirs.append(numpy.get_include())
 
 
 def compilerName():