diff --git a/.gitignore b/.gitignore index abaf0a7bbd..44f7421244 100644 --- a/.gitignore +++ b/.gitignore @@ -11,6 +11,7 @@ ###################### .fuse_hidden* *~ +*swp # Pip generated folders # ######################### @@ -24,4 +25,3 @@ pyansys/Interface.py # Testing Testing/ - diff --git a/doc/index.rst b/doc/index.rst index 153a35d90f..0ed2423256 100644 --- a/doc/index.rst +++ b/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 `_ for details instructions for installing VTK. diff --git a/doc/loading_results.rst b/doc/loading_results.rst index aff579346a..1c1a373382 100644 --- a/doc/loading_results.rst +++ b/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 `_ 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 ------------------------------ diff --git a/pyansys/_version.py b/pyansys/_version.py index 5c2b1d7198..0b7ac6a041 100644 --- a/pyansys/_version.py +++ b/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 diff --git a/pyansys/binary_reader.py b/pyansys/binary_reader.py index e0d4d75114..1801af8544 100755 --- a/pyansys/binary_reader.py +++ b/pyansys/binary_reader.py @@ -846,10 +846,10 @@ 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): """ @@ -857,7 +857,7 @@ def NodalStress(self, rnum): 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 diff --git a/pyansys/cython/_rstHelper.pyx b/pyansys/cython/_rstHelper.pyx index 35c821d018..f0d956fc25 100755 --- a/pyansys/cython/_rstHelper.pyx +++ b/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) \ No newline at end of file + 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) diff --git a/setup.py b/setup.py index dfdc179bfd..3fbcd8d7c2 100644 --- a/setup.py +++ b/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():