diff --git a/doc/ansys_examples.rst b/doc/ansys_examples.rst new file mode 100644 index 0000000000..cbcc108ed4 --- /dev/null +++ b/doc/ansys_examples.rst @@ -0,0 +1,298 @@ +ANSYS Examples +============== +These examples are used to demonstrate ANSYS examples and reading them in using ``pyansys``. + + +Torsional Load on a Bar using SURF154 Elements +---------------------------------------------- +This ANSYS APDL script builds a bar and applies torque to it using SURF154 elements. This is a static analysis example. + +.. code:: + + !---------------------------------------- + ! Input torque applied (moment) + ! Input radius, height, element size... + !---------------------------------------- + TORQUE = 100 + RADIUS = 2 + H_TIP = 2 + HEIGHT = 20 + ELEMSIZE = 1 + PI = acos(-1) + FORCE = 100/RADIUS + PRESSURE = FORCE/(H_TIP*2*PI*RADIUS) + + !---------------------------------------- + ! Define higher-order SOLID186 + ! Define surface effect elements SURF154 + ! which is used to apply torque + ! as a tangential pressure + !---------------------------------------- + /prep7 + et, 1, 186 + et, 2, 154 + r,1, + r,2, + + !---------------------------------------- + ! Aluminum properties (or something) + !---------------------------------------- + mp,ex,1,10e6 + mp,nuxy,1,.3 + mp,dens,1,.1/386.1 + mp,dens,2,0 + + !---------------------------------------- + ! Simple cylinder + !---------------------------------------- + *do, ICOUNT, 1, 4 + cylind,RADIUS,,HEIGHTH_TIP,HEIGHT,90*(ICOUNT-1),90*ICOUNT + *enddo + + + nummrg,kp + lsel,s,loc,x,0 + + lsel,r,loc,y,0 + lsel,r,loc,z,0,HEIGHT-H_TIP + lesize,all,ELEMSIZE*2 + mshape,0 + mshkey,1 + esize,ELEMSIZE + allsel,all + VSWEEP, ALL + csys,1 + asel,s,loc,z,HEIGHT-H_TIP+0.0001,HEIGHT0.0001 + asel,r,loc,x,RADIUS + local,11,1 + csys,0 + aatt,2,2,2,11 + amesh,all + finish + /solu + antype,static,new + eqslv,pcg,1e-8 + + !---------------------------------------- + ! Apply tangential pressure + !---------------------------------------- + esel,s,type,,2 + sfe,all,2,pres,,PRESSURE + + !---------------------------------------- + ! Constrain bottom of cylinder/rod + !---------------------------------------- + asel,s,loc,z,0 + nsla,s,1 + d,all,all + allsel,all + /psf,pres,,2 + /pbc,u,1 + /title, Simple torsional example + solve + finish + /post1 + set,last + fsum + esel,u,type,,2 + SAVE + +Read and plot the results within python using: +.. code:: python + + import pyansys + result = pyansys.ResultReader('file.rst') + + # node numbers and nodal stress + nodennum, stress = result.NodalStress(0) + + # stress at each element + element_stress, elemnum, enode = result.ElementStress(0) + + # plot result + result.PlotNodalResult(0) + + + +Spotweld SHELL181 Example +------------------------- +This ANSYS APDL example demonstrates how to model spot welding on three thin sheets of metal. + +.. code:: + + !---------------------------------------- + ! Example problem for demonstrating + ! Spotweld technology + !---------------------------------------- + ! + !---------------------------------------- + ! Originated in 9.0 JJDoyle 2004/09/01 + !---------------------------------------- + /prep7 + /num,0 + /pnum,area,1 + + k,1,2,10, + k,2,10,10 + k,3,10,0.15 + k,4,14,0.15 + ! + l,1,2 + l,2,3 + l,3,4 + lfillt,1,2,3 + lfillt,2,3,2 + ! + k,9,, + k,10,11, + k,11,15, + l,9,10 + l,10,11 + + k,12,,10 + lsel,s,,,6,7 + AROTAT,all,,,,,,9,12,12,1, + + lsel,s,,,1,5 + AROTAT,all,,,,,,9,12,12,1, + areverse,1 + areverse,2 + + asel,s,,,3,7 + ARSYM,Y,all, , , ,0,0 + allsel + + !******** + !define weld location with hardpoint + !******** + HPTCREATE,AREA,7,0,COORD,12.9,0.15,-1.36, + + /view,1,1,1,1 + gplo + ! + et,1,181 + r,1,0.15 + r,2,0.1 + ! + mp,ex,1,30e6 + mp,prxy,1,0.3 + ! + esize,0.25 + real,1 + amesh,1 + amesh,2 + real,2 + asel,s,,,3,12 + amesh,all + ! + lsel,s,,,1,9 + lsel,a,,,12,17 + lsel,a,,,26,38,3 + lsel,a,,,24,36,3 + nsll,s,1 + wpstyle,0.05,0.1,-1,1,0.003,0,0,,5 + WPSTYLE,,,,,,,,1 + wpro,,-90.000000, + CSWPLA,11,1,1,1, + csys,11 + nrotat,all + d,all,uy + d,all,rotx + + csys,0 + + lsel,s,,,23 + nsll,s,1 + d,all,uz + + lsel,s,,,17 + nsll,s,1 + d,all,uz,4 + + ALLSEL + /view,1,1,1,1 + /eshape,1 + ksel,s,,,33 + nslk,s,1 + *get,sw_node,node,,num,max + + /solu + allsel + nlgeom,on + time,4 + nsubst,10,25,5 + outres,all,all + fini + + !------------------------------------ + !build flex spotweld with BEAM188, run the solution, + !and post process results + !------------------------------------ + fini + allsel + /prep7 + mp,ex,2,28e6 + mp,prxy,2,0.3 + ! + SECTYPE,2,beam,csolid + SECDATA,0.25 + ! + et,2,188 + type,2 + mat,2 + secnum,2 + + SWGEN,sweld1,0.50,7,2,sw_node,, + SWADD,sweld1,,12 + + /solu + allsel + nlgeom,on + time,4 + nsubst,10,25,5 + outres,all,all + solve + FINISH + +Here's the Python script using ``pyansys`` to access the results after running the ANSYS analysis. + +.. code:: python + + import pyansys + + # Open the result file and plot the displacement of time step 3 + resultfile = os.path.join('file.rst') + result = pyansys.ResultReader(resultfile) + result.PlotNodalResult(2) + +.. figure:: ./images/spot_disp.png + :width: 300pt + + Spot Weld: Displacement + +Get the nodal and element component stress at time step 0. Plot the stress in the Z direction. + +.. code:: python + + nodenum, stress = result.NodalStress(0) + element_stress, elemnum, enode = result.ElementStress(0) + + # plot the Z direction stress (the stress at the contact element simulating + # the spot weld) + result.PlotNodalStress(0, 'Sz') + +.. figure:: ./images/spot_sz.png + :width: 300pt + + Spot Weld: Z Stress + +.. code:: python + + # Get the principal nodal stress and plot the von Mises Stress + nnum, pstress = result.PrincipalNodalStress(0) + result.PlotPrincipalNodalStress(0, 'SEQV') + +.. figure:: ./images/spot_seqv.png + :width: 300pt + + Spot Weld: von Mises Stress diff --git a/doc/conf.py b/doc/conf.py index 302e4dc1c8..c37c86bad9 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -51,10 +51,10 @@ # built documents. # # The short X.Y version. -version = u'0.18' +version = u'0.23' # The full version, including alpha/beta/rc tags. -release = u'0.18.0' +release = u'0.23.0' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/doc/images/spot_disp.png b/doc/images/spot_disp.png new file mode 100644 index 0000000000..b4dd959363 Binary files /dev/null and b/doc/images/spot_disp.png differ diff --git a/doc/images/spot_seqv.png b/doc/images/spot_seqv.png new file mode 100644 index 0000000000..6112280452 Binary files /dev/null and b/doc/images/spot_seqv.png differ diff --git a/doc/images/spot_sz.png b/doc/images/spot_sz.png new file mode 100644 index 0000000000..a093e3301d Binary files /dev/null and b/doc/images/spot_sz.png differ diff --git a/doc/index.rst b/doc/index.rst index 0ed2423256..8ce867b4f9 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -12,6 +12,7 @@ Contents :maxdepth: 2 examples + ansys_examples loading_results loading_km @@ -34,6 +35,7 @@ Minimum requirements are numpy to extract results from a results file. To conver See `Install VTK `_ for details instructions for installing VTK. + License ------- This module is licensed under the MIT license. See the license file for more details. diff --git a/doc/loading_results.rst b/doc/loading_results.rst index 63e8f80aab..d9e829b59c 100644 --- a/doc/loading_results.rst +++ b/doc/loading_results.rst @@ -91,16 +91,16 @@ Stress can be obtained as well using the below code. The nodal stress is comput # obtain the component node averaged stress for the first result # organized with one [Sx, Sy Sz, Sxy, Syz, Sxz] entry for each node - stress = result.NodalStress(0) # results in a np array (nnod x 6) + nodenum, stress = result.NodalStress(0) # results in a np array (nnod x 6) # 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) + nodenum, pstress = result.PrincipalNodalStress(0) result.PlotPrincipalNodalStress(0, 'SEQV') -Element stress can be obtained using the following segment of code. Ensure that the element results are expanded within ANSYS with:: +Element stress can be obtained using the following segment of code. Ensure that the element results are expanded for a modal analysis within ANSYS with:: /SOLU MXPAND, ALL, , , YES diff --git a/pyansys/_version.py b/pyansys/_version.py index d53231e73d..33489c8d95 100644 --- a/pyansys/_version.py +++ b/pyansys/_version.py @@ -1,6 +1,5 @@ # major, minor, patch -version_info = 0, 22, 6 - +version_info = 0, 23, 0 # Nice string for the version __version__ = '.'.join(map(str, version_info)) diff --git a/pyansys/binary_reader.py b/pyansys/binary_reader.py index 49061d065f..a0ed63a554 100755 --- a/pyansys/binary_reader.py +++ b/pyansys/binary_reader.py @@ -55,7 +55,7 @@ ptrEPT - pointer to element temperatures ptrESF - pointer to element surface stresses ptrEDI - pointer to diffusion strains -ptrETB - pointer to ETABLE items(post1 only +ptrETB - pointer to ETABLE items(post1 only) ptrECT - pointer to contact data ptrEXY - pointer to integration point locations ptrEBA - pointer to back stresses @@ -63,9 +63,13 @@ ptrMNL - pointer to material nonlinear record """ +# element types with stress outputs +validENS = [45, 92, 95, 181, 185, 186, 187] + class FullReader(object): - """Object to store the results of an ANSYS full file. + """ + Stores the results of an ANSYS full file. Parameters ---------- @@ -469,8 +473,8 @@ def PlotCyclicNodalResult( sec_sol = np.real(r) * np.cos(i * hindex * alpha + phase) -\ np.imag(r) * np.sin(i * hindex * alpha + phase) - # adjust the "direction" of the x and y vectors as they're being - # rotated + # adjust the "direction" of the x and y vectors as they're + # being rotated s_x = sec_sol[:, 0] * np.cos(alpha * i + phase) -\ sec_sol[:, 1] * np.sin(alpha * i + phase) s_y = sec_sol[:, 0] * np.sin(alpha * i + phase) +\ @@ -614,17 +618,13 @@ def PlotNodalResult(self, rnum, comp='norm', as_abs=False, label=''): if as_abs: d = np.abs(d) - # Generate plot -# text = 'Result {:d} at {:f}'.format(rnum + 1, self.tvalues[rnum]) - # setup text ls_table = self.resultheader['ls_table'] - freq = self.GetTimeValues()[rnum] + timevalue = self.GetTimeValues()[rnum] text = 'Cumulative Index: {:3d}\n'.format(ls_table[rnum, 2]) - text += 'Loadstep: {:3d}\n'.format(ls_table[rnum, 0]) - text += 'Substep: {:3d}\n'.format(ls_table[rnum, 1]) -# text += 'Harmonic Index: {:3d}\n'.format(hindex) - text += 'Frequency: {:10.4f} Hz'.format(freq) + text += 'Loadstep: {:3d}\n'.format(ls_table[rnum, 0]) + text += 'Substep: {:3d}\n'.format(ls_table[rnum, 1]) + text += 'Time Value: {:10.4f}'.format(timevalue) plobj = vtkInterface.PlotClass() plobj.AddMesh(self.grid, scalars=d, stitle=stitle, @@ -845,8 +845,8 @@ def StoreGeometry(self): 'e_rcon': np.ones_like(enum)} # store the reference array - cell_type = ['45', '95', '185', '186', '92', '187', '154'] - result = _parser.Parse(self.geometry, False, cell_type) + valid_types = ['45', '95', '185', '186', '92', '187', '154', '181'] + result = _parser.Parse(self.geometry, False, valid_types) cells, offset, cell_type, self.numref, _, _, _ = result # catch -1 @@ -854,11 +854,20 @@ def StoreGeometry(self): # Create vtk object if vtk installed if vtkloaded: + + element_type = np.zeros_like(etype) + for key, typekey in ekey: + # if str(typekey) in valid_types: + element_type[etype == key] = typekey + + # validmask = element_type != 0 + # element_type = element_type[validmask] nodes = nloc[:, :3] self.quadgrid = vtkInterface.UnstructuredGrid(offset, cells, cell_type, nodes) self.quadgrid.AddCellScalars(enum, 'ANSYS_elem_num') self.quadgrid.AddPointScalars(nnum, 'ANSYSnodenum') + self.quadgrid.AddCellScalars(element_type, 'Element Type') self.grid = self.quadgrid.LinearGridCopy() # get edge nodes @@ -893,60 +902,50 @@ def ElementSolutionHeader(self, rnum): self.nsets)) # Read a result - f = open(self.filename, 'rb') - - f.seek((rpointers[rnum] + 1) * 4) # item 20 - # solheader = np.fromfile(f, endian + 'i', 200) - - # Seek to result table and to get pointer to DOF results of result - # f.seek((rpointers[rnum] + 13) * 4) # item 12 32-bit pointer to element solution - # ptrESL = np.fromfile(f, endian + 'i', 1)[0] - - # key to extrapolate integration - f.seek((rpointers[rnum] + 17) * 4) # item 16 - rxtrap = np.fromfile(f, endian + 'i', 1)[0] - # point results to nodes - # = 0 - move - # = 1 - extrapolate unless active - # non-linear - # = 2 - extrapolate always - # print(rxtrap) - if rxtrap == 0: - warnings.warn('Strains and stresses are moved and not ' + - 'extrapolated nodal stress calculations will ' + - 'be incorrect') - - f.seek((rpointers[rnum] + 120) * 4) # item 122 64-bit pointer to element solutio - ptrESL = np.fromfile(f, endian + 'i8', 1)[0] - - if not ptrESL: - f.close() - raise Exception('No element solution in result set %d' % (rnum + 1)) - - # Seek to element result header - element_rst_ptr = rpointers[rnum] + ptrESL + 2 - f.seek(element_rst_ptr * 4) - - # element index table - ele_ind_table = np.fromfile(f, endian + 'i8', nelm) - ele_ind_table += element_rst_ptr - # Each element header contains 25 records for the individual results - # get the location of the nodal stress - table_index = e_table.index('ptrENS') - - # check number of records to read (differs with each version) - f.seek((ele_ind_table[0] + table_index) * 4) - ptrENS = np.fromfile(f, endian + 'i', 1)[0] - - nnode_elem = nodstr[etype[0]] - f.seek((ele_ind_table[0] + ptrENS - 2) * 4) - nitem = np.fromfile(f, endian + 'i', 1)[0] / nnode_elem - f.close() - - return table_index, ele_ind_table, nodstr, etype, nitem + with open(self.filename, 'rb') as f: - def NodalStress(self, rnum, debug=False): + f.seek((rpointers[rnum] + 1) * 4) # item 20 + # solheader = np.fromfile(f, endian + 'i', 200) + + # key to extrapolate integration + f.seek((rpointers[rnum] + 17) * 4) # item 16 + rxtrap = np.fromfile(f, endian + 'i', 1)[0] + # point results to nodes + # = 0 - move + # = 1 - extrapolate unless active + # non-linear + # = 2 - extrapolate always + # print(rxtrap) + if rxtrap == 0: + warnings.warn('Strains and stresses are moved and not ' + + 'extrapolated nodal stress calculations will ' + + 'be incorrect') + + # item 122 64-bit pointer to element solution + f.seek((rpointers[rnum] + 120) * 4) + ptrESL = np.fromfile(f, endian + 'i8', 1)[0] + + if not ptrESL: + f.close() + raise Exception('No element solution in result set %d' % (rnum + 1)) + + # Seek to element result header + element_rst_ptr = rpointers[rnum] + ptrESL + 2 + f.seek(element_rst_ptr * 4) + + # element index table + ele_ind_table = np.fromfile(f, endian + 'i8', nelm) + ele_ind_table += element_rst_ptr + # Each element header contains 25 records for the individual + # results. Get the location of the nodal stress + table_index = e_table.index('ptrENS') + + return table_index, ele_ind_table, nodstr, etype + + def NodalStress(self, rnum): """ + Equivalent ANSYS command: PRNSOL, S + Retrives the component stresses for each node in the solution. The order of the results corresponds to the sorted node numbering. @@ -963,57 +962,47 @@ def NodalStress(self, rnum, debug=False): Returns ------- + nodenum : numpy.ndarray + Node numbers of the result. + stress : numpy.ndarray Stresses at Sx Sy Sz Sxy Syz Sxz averaged at each corner node. For the corresponding node numbers, see "result.edge_node_num" where result is the result object. """ - table_index, ele_ind_table, nodstr, etype, nitem = self.ElementSolutionHeader(rnum) + element_stress, elemnum, enode = self.ElementStress(0) - # number of nodes - nnod = self.resultheader['nnod'] - if nitem == 22: # double precision < v14.5 - nitem = 11 - ele_data_arr = np.zeros((nnod, 6), np.float64) - _rstHelper.LoadStressDouble(self.filename, table_index, - ele_ind_table, - nodstr.astype(c_int64), - etype.astype(c_int64), nitem, - ele_data_arr, - self.edge_idx.astype(c_int64)) - elif nitem == 11 or nitem == 6: # single precision < v14.5 - ele_data_arr = np.zeros((nnod, 6), np.float32) - _rstHelper.LoadStress(self.filename, table_index, ele_ind_table, - nodstr.astype(c_int64), - etype.astype(c_int64), nitem, - ele_data_arr, - self.edge_idx.astype(c_int64)) + # nodes for each element + nnum = np.hstack(enode) + nodenum = np.unique(nnum) - else: - raise Exception('Invalid nitem. Unable to load nodal stresses') + # stack the element stresses + arr = np.vstack(element_stress) - # Average based on the edges of elements - enode = self.edge_node_num_idx - ntimes = np.bincount(self.edge_idx)[enode] - s_node = ele_data_arr[enode] - s_node /= ntimes.reshape((-1, 1)) + # determine the number of times each node occurs in the results + arr_ones = np.ones_like(arr) + ncount = np.bincount(nnum, weights=arr_ones[:, 0]) + mask = ncount != 0 - if debug: - return s_node, ele_data_arr - else: - return s_node + # sum and weight the stress at each node + stress = np.empty((nodenum.size, 6)) + for i in range(6): + stress[:, i] = np.bincount(nnum, weights=arr[:, i])[mask] + + stress /= ncount[mask].reshape(-1, 1) + return nodenum, stress def ElementStress(self, rnum): """ - Retrives the component stresses for each node in the solution. + Equivalent ANSYS command: PRESOL, S - The order of the results corresponds to the sorted node numbering. + Retrives the component stresses for each node in the solution. 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. + stress for each element at each node. Due to the discontinuities + across elements, stresses at nodes will vary based on the element + they are evaluated from. Parameters ---------- @@ -1032,55 +1021,74 @@ def ElementStress(self, rnum): Node numbers corresponding to each element's stress results. One list entry for each element + Notes + ----- + Only elements that output stress will have their stresses reported. + See the ANSYS element guide. + """ - table_index, ele_ind_table, nodstr, etype, nitem = self.ElementSolutionHeader(rnum) - nnod = self.resultheader['nnod'] + result = self.ElementSolutionHeader(rnum) + table_index, ele_ind_table, nodstr, etype = result - n = nodstr[etype].sum() # nodes per element * elements - overflow = 100 # if nitem > 6 - if nitem == 22: # double precision < v14.5 - raise Exception - # nitem = 11 - ele_data_arr = np.zeros((n + overflow, 6), np.float64) # Sx Sy Sz Sxy Syz Sxz - _rstHelper.LoadElementStressDouble(self.filename, table_index, - ele_ind_table, - nodstr.astype(c_int64), - etype.astype(c_int64), nitem, - ele_data_arr, - self.edge_idx.astype(c_int64)) - elif nitem == 11 or nitem == 6: # single precision < v14.5 - ele_data_arr = np.zeros((n + overflow, 6), np.float32) # Sx Sy Sz Sxy Syz Sxz - _rstHelper.LoadElementStress(self.filename, table_index, + # certain element types do not output stress + elemtype = self.grid.GetCellScalars('Element Type') + validmask = np.in1d(elemtype, validENS) + validmask[:] = True + + ele_ind_table = ele_ind_table[validmask] + etype = etype[validmask].astype(c_int64) + + # load in raw results + nnode = nodstr[etype] + nelemnode = nnode.sum() + ver = float(self.resultheader['verstring']) + if ver >= 14.5: + if self.resultheader['rstsprs'] != 0: + nitem = 6 + else: + nitem = 11 + ele_data_arr = np.empty((nelemnode, nitem), np.float32) + _rstHelper.ReadElementStress(self.filename, table_index, ele_ind_table, nodstr.astype(c_int64), - etype.astype(c_int64), nitem, + etype, ele_data_arr, - self.edge_idx.astype(c_int64)) + self.edge_idx.astype(c_int64), + nitem) + if nitem != 6: + ele_data_arr = ele_data_arr[:, :6] else: - raise Exception('Invalid nitem. Unable to load nodal stresses') + ele_data_arr = np.empty((nelemnode, 6), np.float64) + _rstHelper.ReadElementStressDouble(self.filename, table_index, + ele_ind_table, + nodstr.astype(c_int64), + etype, + ele_data_arr, + self.edge_idx.astype(c_int64)) - # return result without overflow - element_stress = ele_data_arr[:n] - nnode = nodstr[etype] - element_stress = [] + splitind = np.cumsum(nnode) + element_stress = np.split(ele_data_arr, splitind[:-1]) + + # reorder list using sorted indices + enum = self.grid.GetCellScalars('ANSYS_elem_num')[validmask] + sidx = np.argsort(enum) + element_stress = [element_stress[i] for i in sidx] + + elem = self.geometry['elem'][validmask] enode = [] - count = 0 - for i in range(etype.size): - element_stress.append(ele_data_arr[count:count + nnode[i]]) - count += nnode[i] - enode.append(self.geometry['elem'][i, :nnode[i]]) + for i in sidx: + enode.append(elem[i, :nnode[i]]) - return element_stress, self.geometry['enum'], enode + # Get element numbers + elemnum = self.geometry['enum'][self.sidx_elem] + return element_stress, elemnum, enode 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 @@ -1088,6 +1096,9 @@ def PrincipalNodalStress(self, rnum): Returns ------- + nodenum : numpy.ndarray + Node numbers of the result. + pstress : numpy.ndarray Principal stresses, stress intensity, and equivalant stress. [sigma1, sigma2, sigma3, sint, seqv] @@ -1103,12 +1114,12 @@ def PrincipalNodalStress(self, rnum): """ # get component stress - stress = self.NodalStress(rnum) + nodenum, stress = self.NodalStress(rnum) # compute principle stress if stress.dtype != np.float32: stress = stress.astype(np.float32) - return _rstHelper.ComputePrincipalStress(stress) + return nodenum, _rstHelper.ComputePrincipalStress(stress) def PlotPrincipalNodalStress(self, rnum, stype): """ @@ -1132,6 +1143,26 @@ def PlotPrincipalNodalStress(self, rnum, stype): cpos : list VTK camera position. + """ + stress = self.PrincipleStressForPlotting(rnum, stype) + + # Generate plot + stitle = 'Nodal Stress\n{:s}'.format(stype) + plobj = vtkInterface.PlotClass() + plobj.AddMesh(self.grid, scalars=stress, stitle=stitle, + flipscalars=True, interpolatebeforemap=True) + + text = 'Result %d at %f' % (rnum + 1, self.tvalues[rnum]) + plobj.AddText(text) + + return plobj.Plot(), stress + + def PrincipleStressForPlotting(self, rnum, stype): + """ + returns stress used to plot + + further documentation forthcoming + """ stress_types = ['S1', 'S2', 'S3', 'SINT', 'SEQV'] if stype not in stress_types: @@ -1139,23 +1170,20 @@ def PlotPrincipalNodalStress(self, rnum, stype): sidx = stress_types.index(stype) - # create a zero array for all nodes - stress = np.zeros(self.resultheader['nnod']) + # don't display elements that can't store stress (!) + # etype = self.grid.GetCellScalars('Element Type') + # valid = (np.in1d(etype, validENS)).nonzero()[0] + # grid = self.grid.ExtractSelectionCells(valid) + grid = self.grid # bypassed (for now) # 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) + nodenum = grid.GetPointScalars('ANSYSnodenum') + stress_nnum, edge_stress = self.PrincipalNodalStress(rnum) + temp_arr = np.zeros((nodenum.max() + 1, 5)) + temp_arr[stress_nnum] = edge_stress - # 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 + # find matching edge nodes + return temp_arr[nodenum, sidx] def PlotNodalStress(self, rnum, stype): """ @@ -1171,41 +1199,49 @@ def PlotNodalStress(self, rnum, stype): ---------- rnum : interger Result set using zero based indexing. + stype : string Stress type from the following list: [Sx Sy Sz Sxy Syz Sxz] + Returns + ------- + cpos : list + 3 x 3 vtk camera position. """ - stress_types = ['Sx', 'Sy', 'Sz', 'Sxy', 'Syz', 'Sxz', 'Seqv'] + stress_types = ['sx', 'sy', 'sz', 'sxy', 'syz', 'sxz', 'seqv'] + stype = stype.lower() if stype not in stress_types: - raise Exception('Stress type not in \n' + - "['Sx', 'Sy', 'Sz', 'Sxy', 'Syz', 'Sxz']") - - sidx = ['Sx', 'Sy', 'Sz', 'Sxy', 'Syz', 'Sxz'].index(stype) + 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']) + # don't display elements that can't store stress + # etype = self.grid.GetCellScalars('Element Type') + # valid = (np.in1d(etype, validENS)).nonzero()[0] + # grid = self.grid.ExtractSelectionCells(valid) + grid = self.grid # bypassed for now # Populate with nodal stress at edge nodes - stress[self.edge_node_num_idx] = self.NodalStress(rnum)[:, sidx] - stitle = 'Nodal Stress\n{:s}'.format(stype) + nodenum = grid.GetPointScalars('ANSYSnodenum') + stress_nnum, edge_stress = self.NodalStress(rnum) + temp_arr = np.zeros((nodenum.max() + 1, 6)) + temp_arr[stress_nnum] = edge_stress + stress = temp_arr[nodenum, sidx] + + # stress[mask] = edge_stress[:, sidx] + stitle = 'Nodal Stress\n{:s}'.format(stype.capitalize()) # Generate plot plobj = vtkInterface.PlotClass() - plobj.AddMesh(self.grid, scalars=stress, stitle=stitle, + plobj.AddMesh(grid, scalars=stress, stitle=stitle, flipscalars=True, interpolatebeforemap=True) - - text = 'Result {:d} at {:f}'.format(rnum + 1, self.tvalues[rnum]) + text = 'Result %d at %f' % (rnum + 1, self.tvalues[rnum]) plobj.AddText(text) - - cpos = plobj.Plot() # store camera position - - return cpos + return plobj.Plot() def SaveAsVTK(self, filename, binary=True): """ - Writes all appends all results to an unstructured grid and writes it to - disk. + Appends all results to an unstructured grid and writes it to disk. The file extension will select the type of writer to use. '.vtk' will use the legacy writer, while '.vtu' will select the VTK XML writer. @@ -1216,6 +1252,7 @@ def SaveAsVTK(self, filename, binary=True): Filename of grid to be written. The file extension will select the type of writer to use. '.vtk' will use the legacy writer, while '.vtu' will select the VTK XML writer. + binary : bool, optional Writes as a binary file by default. Set to False to write ASCII @@ -1229,15 +1266,18 @@ def SaveAsVTK(self, filename, binary=True): grid = self.grid.Copy() for i in range(self.nsets): - # Nodal Results + # Nodal results val = self.GetNodalResult(i) grid.AddPointScalars(val, 'NodalResult{:03d}'.format(i)) - # Nodal Stress values are only valid - stress = self.NodalStress(i) - val = np.zeros((grid.GetNumberOfPoints(), stress.shape[1])) - val[self.edge_node_num_idx] = stress - grid.AddPointScalars(val, 'NodalStress{:03d}'.format(i)) + # Populate with nodal stress at edge nodes + nodenum = self.grid.GetPointScalars('ANSYSnodenum') + stress_nnum, edge_stress = self.NodalStress(i) + temp_arr = np.zeros((nodenum.max() + 1, 6)) + temp_arr[stress_nnum] = edge_stress + stress = temp_arr[nodenum] + + grid.AddPointScalars(stress, 'NodalStress{:03d}'.format(i)) grid.Write(filename) @@ -1254,9 +1294,9 @@ def GetResultInfo(filename): Returns ------- resultheader : dictionary + Result header """ - f = open(filename, 'rb') # initialize result header dictionary @@ -1340,6 +1380,9 @@ def GetResultInfo(filename): # cyclic symmetry coordinate system resultheader['csCord'] = rheader[21] + # bitmask for suppressed items (item 40) + resultheader['rstsprs'] = rheader[40 - 1] + # Read nodal equivalence table f.seek((ptrNODl + 2) * 4) # Start of pointer, then empty, then data resultheader['neqv'] = np.fromfile( @@ -1411,38 +1454,3 @@ def delete_row_csc(mat, i): mat.indptr[i:] -= n mat.indptr = mat.indptr[:-1] mat._shape = (mat._shape[0] - 1, mat._shape[1]) - - -# ============================================================================= -# load stress debug using numpy -# ============================================================================= -# #%% numpy debug -# f = open(self.filename) -# nstresses = self.edge_idx.size -# stresses = np.empty((nstresses, 6), np.float32) -# c = 0 -# for i in range(len(ele_ind_table)): -# # Element nodal stresses, ptrENS, is the third item in the table -# f.seek((ele_ind_table[i] + table_index)*4) -# ptrENS = np.fromfile(f, endian + 'i', 1)[0] -# -# # read the stresses evaluated at the intergration points or nodes -# nnode_elem = nodstr[etype[i]] -# -# f.seek((ele_ind_table[i] + ptrENS)*4) -# stress = np.fromfile(f, endian + 'f', nnode_elem*nitem).reshape((-1, nitem))#[:, sidx] -# -# # store stresses -# stresses[c:c + nnode_elem] = stress[:, :6] -# c += nnode_elem -# -# # close file -# f.close() -# -# # Average the stresses for each element at each node -# enode = self.edge_node_num_idx -# s_node = np.empty((enode.size, 6), np.float32) -# for i in range(6): -# s_node[:, i] = np.bincount(self.edge_idx, weights=stresses[:, i])[enode] -# ntimes = np.bincount(self.edge_idx)[enode] -# s_node /= ntimes.reshape((-1, 1)) diff --git a/pyansys/cython/_parser.pyx b/pyansys/cython/_parser.pyx index 160a1aef91..acaf29dae8 100755 --- a/pyansys/cython/_parser.pyx +++ b/pyansys/cython/_parser.pyx @@ -10,6 +10,9 @@ from libc.stdint cimport int32_t, int64_t ctypedef unsigned char uint8 # VTK numbering for vtk cells +cdef uint8 VTK_EMPTY_CELL = 0 +cdef uint8 VTK_VERTEX = 1 +cdef uint8 VTK_LINE = 3 cdef uint8 VTK_TRIANGLE = 5 cdef uint8 VTK_QUAD = 9 cdef uint8 VTK_QUADRATIC_TRIANGLE = 22 @@ -40,6 +43,27 @@ typeB[0] = 92 typeB[1] = 187 +cdef inline void StoreLine(int64_t [::1] offset, int64_t *ecount, int64_t *ccount, + int64_t [::1] cells, uint8 [::1] cell_type, + int64_t [::1] numref, int [:, ::1] elem, int i, int lin): + """ + Stores surface triangle vtk cell. Element may be quadradic or linear + """ + # Populate offset array + offset[ecount[0]] = ccount[0] + + if lin: + cells[ccount[0]] = 2; ccount[0] += 1 + + # Populate cell array while renumbering nodes + for j in range(2): + cells[ccount[0]] = numref[elem[i, j]]; ccount[0] += 1 + + # Populate cell type array + cell_type[ecount[0]] = VTK_LINE + + ecount[0] += 1 + cdef inline void StoreSurfTri(int64_t [::1] offset, int64_t *ecount, int64_t *ccount, int64_t [::1] cells, uint8 [::1] cell_type, int64_t [::1] numref, int [:, ::1] elem, int i, int lin): @@ -387,6 +411,7 @@ def Parse(raw, pyforce_linear, allowable_types): Parses raw cdb data from downstream conversion to a vtk unstructured grid """ cdef int force_linear = pyforce_linear + cdef int i, j, k, lin # ANSYS element type definitions cdef int [4] typeA @@ -425,11 +450,19 @@ def Parse(raw, pyforce_linear, allowable_types): else: typeB[1] = -1 - cdef int [1] typeC - if '154' in allowable_types: - typeC[0] = 154 - else: - typeC[0] = -1 + # shell types + cdef int [2] typeC + for i, atype in enumerate(['154', '181']): + if atype in allowable_types: + typeC[i] = int(atype) + else: + typeC[i] = -1 + + cdef int [1] linetype + for i, atype in enumerate(['188']): + if atype in allowable_types: + linetype[i] = int(atype) + cdef int [:, ::1] ekey = raw['ekey'].astype(ctypes.c_int) cdef int [:, ::1] elem = raw['elem'].astype(ctypes.c_int) @@ -438,7 +471,6 @@ def Parse(raw, pyforce_linear, allowable_types): cdef int [::1] raw_enum = raw['enum'].astype(ctypes.c_int) cdef int [::1] raw_rcon = raw['e_rcon'].astype(ctypes.c_int) - cdef int i, j, k, lin cdef int nelem = elem.shape[0] cdef int nnode = nnum.shape[0] @@ -492,9 +524,12 @@ def Parse(raw, pyforce_linear, allowable_types): # can read cdef int64_t ccount = 0 # cell/offset counter cdef int64_t ecount = 0 # element number counter + cdef int64_t cstart cdef int elem_etype for i in range(nelem): elem_etype = elem_type[etype[i]] + cstart = ccount + for j in range(4): if elem_etype == typeA[j]: enum[ecount] = raw_enum[i] @@ -526,7 +561,7 @@ def Parse(raw, pyforce_linear, allowable_types): elem, i, lin) break # Continue to next element - + # Test for element type B for j in range(2): if elem_etype == typeB[j]: @@ -542,26 +577,60 @@ def Parse(raw, pyforce_linear, allowable_types): StoreTet_TypeB(offset, &ecount, &ccount, cells, cell_type, numref, elem, i, lin) - break # Continue to next element + break # Continue to next element + + # test if surface element + for j in range(2): + if elem_etype == typeC[j]: + enum[ecount] = raw_enum[i] + etype_out[ecount] = elem_etype + rcon[ecount] = raw_rcon[i] - # test if surface element SURF154 - if elem_etype == typeC[0]: + # if force_linear: + lin = 1 + # else: + # lin = elem[i, 4] == -1 + + # check if this is a triangle + if elem[i, 2] == elem[i, 3]: + StoreSurfTri(offset, &ecount, &ccount, cells, cell_type, + numref, elem, i, lin) + else: + StoreSurfQuad(offset, &ecount, &ccount, cells, cell_type, + numref, elem, i, lin) + break # Continue to next element + # continue + + # test if line element + for j in range(1): + if elem_etype == linetype[j]: + enum[ecount] = raw_enum[i] + etype_out[ecount] = elem_etype + rcon[ecount] = raw_rcon[i] + + lin = 1 + StoreLine(offset, &ecount, &ccount, cells, cell_type, + numref, elem, i, lin) + + # add null element if element isn't in the allowable types + if cstart == ccount: enum[ecount] = raw_enum[i] etype_out[ecount] = elem_etype rcon[ecount] = raw_rcon[i] + offset[ecount] = ccount - if force_linear: - lin = 1 - else: - lin = elem[i, 4] == -1 + # add null cell + cells[ccount] = 0; ccount += 1 + + cell_type[ecount] = VTK_EMPTY_CELL + ecount += 1 + + # Populate cell array while renumbering nodes + # for j in range(2): + # cells[ccount[0]] = numref[elem[i, j]]; ccount[0] += 1 + + # Populate cell type array - # check if this is a triangle - if elem[i, 2] == elem[i, 3]: - StoreSurfTri(offset, &ecount, &ccount, cells, cell_type, - numref, elem, i, lin) - else: - StoreSurfQuad(offset, &ecount, &ccount, cells, cell_type, - numref, elem, i, lin) return np.asarray(cells[:ccount]), np.asarray(offset[:ecount]), \ diff --git a/pyansys/cython/_rstHelper.pyx b/pyansys/cython/_rstHelper.pyx index f528b69678..531ab6e42c 100755 --- a/pyansys/cython/_rstHelper.pyx +++ b/pyansys/cython/_rstHelper.pyx @@ -194,6 +194,73 @@ def LoadElementStressDouble(filename, py_table_index, int64_t [::1] ele_ind_tabl fclose(cfile) +def ReadElementStress(filename, py_table_index, int64_t [::1] ele_ind_table, + int64_t [::1] nodstr,int64_t [::1] etype, + float [:, ::1] ele_data_arr, int64_t [::1] edge_idx, + int nitem): + """ Read element results from ANSYS directly into a numpy array """ + cdef int64_t i, j, k, ind, nread + cdef int64_t table_index = py_table_index + + cdef FILE* cfile + cdef bytes py_bytes = filename.encode() + cdef char* c_filename = py_bytes + cfile = fopen(c_filename, 'rb') + + cdef int64_t ele_table, ptrENS, nnode_elem + cdef int64_t c = 0 + for i in range(len(ele_ind_table)): + # get location of pointers to element data + ele_table = ele_ind_table[i] + fseek(cfile, (ele_table + table_index)*4, SEEK_SET) + fread(&ptrENS, sizeof(int), 1, cfile) + + # Get the nodes in the element + nnode_elem = nodstr[etype[i]] + + # read the stresses evaluated at the intergration points or nodes + fseek(cfile, (ele_table + ptrENS)*4, SEEK_SET) + + nread = nnode_elem*nitem + fread(&ele_data_arr[c, 0], sizeof(float), nread, cfile) + c += nnode_elem + + fclose(cfile) + + +def ReadElementStressDouble(filename, py_table_index, int64_t [::1] ele_ind_table, + int64_t [::1] nodstr,int64_t [::1] etype, + double [:, ::1] ele_data_arr, int64_t [::1] edge_idx): + """ Read element results from ANSYS directly into a numpy array """ + cdef int64_t i, j, k, ind, nread + cdef int64_t table_index = py_table_index + + cdef FILE* cfile + cdef bytes py_bytes = filename.encode() + cdef char* c_filename = py_bytes + cfile = fopen(c_filename, 'rb') + + cdef int64_t ele_table, ptrENS, nnode_elem + cdef int64_t c = 0 + for i in range(len(ele_ind_table)): + # get location of pointers to element data + ele_table = ele_ind_table[i] + fseek(cfile, (ele_table + table_index)*4, SEEK_SET) + fread(&ptrENS, sizeof(int), 1, cfile) + + # Get the nodes in the element + nnode_elem = nodstr[etype[i]] + + # read the stresses evaluated at the intergration points or nodes + fseek(cfile, (ele_table + ptrENS)*4, SEEK_SET) + + nread = nnode_elem*6 + fread(&ele_data_arr[c, 0], sizeof(double), nread, cfile) + c += nread + + fclose(cfile) + + def LoadStress(filename, py_table_index, int64_t [::1] ele_ind_table, int64_t [::1] nodstr,int64_t [::1] etype, int64_t nitem, float [:, ::1] ele_data_arr, int64_t [::1] edge_idx): diff --git a/pyansys/examples/examples.py b/pyansys/examples/examples.py index e3cf64e3cc..eb41ca4c0e 100755 --- a/pyansys/examples/examples.py +++ b/pyansys/examples/examples.py @@ -1,216 +1,210 @@ -""" -pyansys examples - -""" - -from __future__ import print_function -import os -import inspect -import sys - -import pyansys - - -# get location of this folder and the example files -dir_path = os.path.dirname(os.path.realpath(__file__)) -rstfile = os.path.join(dir_path, 'file.rst') -hexarchivefile = os.path.join(dir_path, 'HexBeam.cdb') -tetarchivefile = os.path.join(dir_path, 'TetBeam.cdb') -fullfile = os.path.join(dir_path, 'file.full') - - -def RunAll(): - """ Runs all the functions within this module """ - testfunctions = [] - for name, obj in inspect.getmembers(sys.modules[__name__]): - if inspect.isfunction(obj) and name != 'RunAll': - testfunctions.append(obj) - - # run all the functions - any(f() for f in testfunctions) - - -def DisplayHexBeam(): - """ Displays a hex beam mesh """ - # Load an archive file - archive = pyansys.ReadArchive(hexarchivefile) - grid = archive.ParseVTK() - grid.Plot() - - -def LoadResult(): - """ - Loads a result file and prints out the displacement of all the nodes from - a modal analysis. - - """ - - # Load result file - result = pyansys.ResultReader(rstfile) - print('Loaded result file with {:d} result sets'.format(result.nsets)) - print('Contains {:d} nodes'.format(len(result.nnum))) - - # display result - disp = result.GetNodalResult(0) - - print('Nodal displacement for nodes 30 to 40 is:') - - for i in range(29, 40): - node = result.nnum[i] - x = disp[i, 0] - y = disp[i, 1] - z = disp[i, 2] - print('{:2d} {:10.6f} {:10.6f} {:10.6f}'.format(node, x, y, z)) - - -def DisplayDisplacement(): - """ Load and plot 1st bend of a hexahedral beam """ - - # get location of this file - fobj = pyansys.ResultReader(rstfile) - - print('Displaying ANSYS Mode 1') - fobj.PlotNodalResult(0, label='Displacement') - - -def DisplayStress(): - """ Load and plot 1st bend of a hexahedral beam """ - - # get location of this file - result = pyansys.ResultReader(rstfile) - - print('Displaying node averaged stress in x direction for Mode 6') - result.PlotNodalStress(5, 'Sx') - - -def LoadKM(): - """ Loads m and k matrices from a full file """ - import numpy as np - - # Create file reader object - fobj = pyansys.FullReader(fullfile) - dofref, k, m = fobj.LoadKM(utri=False) - - # print results - ndim = k.shape[0] - print('Loaded {:d} x {:d} mass and stiffness matrices'.format(ndim, ndim)) - print('\t k has {:d} entries'.format(k.indices.size)) - print('\t m has {:d} entries'.format(m.indices.size)) - - # compute natural frequencies if installed - try: - from scipy.sparse import linalg - except BaseException: - return - - # removing constrained nodes (this is not efficient) - free = np.logical_not(fobj.const).nonzero()[0] - k = k[free][:, free] - m = m[free][:, free] - - import numpy as np - # Solve - w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) - - # System natural frequencies - f = np.real(w)**0.5 / (2 * np.pi) - - print('First four natural frequencies:') - for i in range(4): - print('{:.3f} Hz'.format(f[i])) - - -def SolveKM(): - """ - Loads and solves a mass and stiffness matrix from an ansys full file - """ - import numpy as np - from scipy.sparse import linalg - import pyansys - - # load the mass and stiffness matrices - full = pyansys.FullReader(pyansys.examples.fullfile) - dofref, k, m = full.LoadKM(utri=False) - - # removing constrained nodes (this is not efficient) - free = np.logical_not(full.const).nonzero()[0] - k = k[free][:, free] - m = m[free][:, free] - - # Solve - w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) - - # System natural frequencies - f = (np.real(w))**0.5 / (2 * np.pi) - - #%% Plot result - import vtkInterface - - # Get the 4th mode shape - mode_shape = v[:, 3] # x, y, z displacement for each node - - # create the full mode shape including the constrained nodes - full_mode_shape = np.zeros(dofref.shape[0]) - full_mode_shape[np.logical_not(full.const)] = mode_shape - - # reshape and compute the normalized displacement - disp = full_mode_shape.reshape((-1, 3)) - n = (disp * disp).sum(1)**0.5 - n /= n.max() # normalize - - # load an archive file and create a vtk unstructured grid - archive = pyansys.ReadArchive(pyansys.examples.hexarchivefile) - grid = archive.ParseVTK() - - # plot the normalized displacement -# grid.Plot(scalars=n) - - # Fancy plot the displacement - plobj = vtkInterface.PlotClass() - - # add two meshes to the plotting class. Meshes are copied on load - plobj.AddMesh(grid, style='wireframe') - plobj.AddMesh(grid, scalars=n, stitle='Normalized\nDisplacement', - flipscalars=True) - - # Update the coordinates by adding the mode shape to the grid - plobj.UpdateCoordinates(grid.GetNumpyPoints() + disp / 80, render=False) - plobj.AddText('Cantliver Beam 4th Mode Shape at {:.4f}'.format(f[3]), - fontsize=30) - plobj.Plot() - del plobj - - -def DisplayCellQual(meshtype='tet'): - """ - Displays minimum scaled jacobian of a sample mesh - - For an ANSYS analysis to run properly, the minimum scaled jacobian of each - cell must be greater than 0. - - Parameters - ---------- - meshtype string, optional - Set to 'hex' to display cell quality of a hexahedral meshed beam. - - Returns - ------- - None - - """ - - # load archive file and parse for subsequent FEM queries - if meshtype == 'hex': - archive = pyansys.ReadArchive(hexarchivefile) - else: - archive = pyansys.ReadArchive(tetarchivefile) - - # create vtk object - grid = archive.ParseVTK() - - # get cell quality - qual = pyansys.CellQuality(grid) - - # plot cell quality - grid.Plot(scalars=qual, stitle='Cell Minimum Scaled\nJacobian', - rng=[0, 1]) +""" +pyansys examples + +""" + +from __future__ import print_function +import os +import inspect +import sys + +import numpy as np + +import pyansys +import vtkInterface + + +# get location of this folder and the example files +dir_path = os.path.dirname(os.path.realpath(__file__)) +rstfile = os.path.join(dir_path, 'file.rst') +hexarchivefile = os.path.join(dir_path, 'HexBeam.cdb') +tetarchivefile = os.path.join(dir_path, 'TetBeam.cdb') +fullfile = os.path.join(dir_path, 'file.full') + + +def RunAll(): + """ Runs all the functions within this module """ + testfunctions = [] + for name, obj in inspect.getmembers(sys.modules[__name__]): + if inspect.isfunction(obj) and name != 'RunAll': + testfunctions.append(obj) + + # run all the functions + any(f() for f in testfunctions) + + +def DisplayHexBeam(): + """ Displays a hex beam mesh """ + # Load an archive file + archive = pyansys.ReadArchive(hexarchivefile) + grid = archive.ParseVTK() + grid.Plot() + + +def LoadResult(): + """ + Loads a result file and prints out the displacement of all the nodes from + a modal analysis. + + """ + + # Load result file + result = pyansys.ResultReader(rstfile) + print('Loaded result file with {:d} result sets'.format(result.nsets)) + print('Contains {:d} nodes'.format(len(result.nnum))) + + # display result + disp = result.GetNodalResult(0) + + print('Nodal displacement for nodes 30 to 40 is:') + + for i in range(29, 40): + node = result.nnum[i] + x = disp[i, 0] + y = disp[i, 1] + z = disp[i, 2] + print('{:2d} {:10.6f} {:10.6f} {:10.6f}'.format(node, x, y, z)) + + +def DisplayDisplacement(): + """ Load and plot 1st bend of a hexahedral beam """ + + # get location of this file + fobj = pyansys.ResultReader(rstfile) + + print('Displaying ANSYS Mode 1') + fobj.PlotNodalResult(0, label='Displacement') + + +def DisplayStress(): + """ Load and plot 1st bend of a hexahedral beam """ + + # get location of this file + result = pyansys.ResultReader(rstfile) + + print('Displaying node averaged stress in x direction for Mode 6') + result.PlotNodalStress(5, 'Sx') + + +def LoadKM(): + """ Loads m and k matrices from a full file """ + + # Create file reader object + fobj = pyansys.FullReader(fullfile) + dofref, k, m = fobj.LoadKM(utri=False) + + # print results + ndim = k.shape[0] + print('Loaded {:d} x {:d} mass and stiffness matrices'.format(ndim, ndim)) + print('\t k has {:d} entries'.format(k.indices.size)) + print('\t m has {:d} entries'.format(m.indices.size)) + + # compute natural frequencies if installed + try: + from scipy.sparse import linalg + except BaseException: + return + + # removing constrained nodes (this is not efficient) + free = np.logical_not(fobj.const).nonzero()[0] + k = k[free][:, free] + m = m[free][:, free] + + # Solve + w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) + + # System natural frequencies + f = np.real(w)**0.5 / (2 * np.pi) + + print('First four natural frequencies:') + for i in range(4): + print('{:.3f} Hz'.format(f[i])) + + +def SolveKM(): + """ + Loads and solves a mass and stiffness matrix from an ansys full file + """ + from scipy.sparse import linalg + + # load the mass and stiffness matrices + full = pyansys.FullReader(pyansys.examples.fullfile) + dofref, k, m = full.LoadKM(utri=False) + + # removing constrained nodes (this is not efficient) + free = np.logical_not(full.const).nonzero()[0] + k = k[free][:, free] + m = m[free][:, free] + + # Solve + w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) + + # System natural frequencies + f = (np.real(w))**0.5 / (2 * np.pi) + + # %% Plot result + + # Get the 4th mode shape + mode_shape = v[:, 3] # x, y, z displacement for each node + + # create the full mode shape including the constrained nodes + full_mode_shape = np.zeros(dofref.shape[0]) + full_mode_shape[np.logical_not(full.const)] = mode_shape + + # reshape and compute the normalized displacement + disp = full_mode_shape.reshape((-1, 3)) + n = (disp * disp).sum(1)**0.5 + n /= n.max() # normalize + + # load an archive file and create a vtk unstructured grid + archive = pyansys.ReadArchive(pyansys.examples.hexarchivefile) + grid = archive.ParseVTK() + + # plot the normalized displacement + # grid.Plot(scalars=n) + + # Fancy plot the displacement + plobj = vtkInterface.PlotClass() + + # add two meshes to the plotting class + plobj.AddMesh(grid.Copy(), style='wireframe') + plobj.AddMesh(grid, scalars=n, stitle='Normalized\nDisplacement', + flipscalars=True) + + # Update the coordinates by adding the mode shape to the grid + plobj.UpdateCoordinates(grid.GetNumpyPoints() + disp / 80, render=False) + plobj.AddText('Cantliver Beam 4th Mode Shape at {:.4f}'.format(f[3]), + fontsize=30) + plobj.Plot() + del plobj + + +def DisplayCellQual(meshtype='tet'): + """ + Displays minimum scaled jacobian of a sample mesh + + For an ANSYS analysis to run properly, the minimum scaled jacobian of each + cell must be greater than 0. + + Parameters + ---------- + meshtype string, optional + Set to 'hex' to display cell quality of a hexahedral meshed beam. + + """ + + # load archive file and parse for subsequent FEM queries + if meshtype == 'hex': + archive = pyansys.ReadArchive(hexarchivefile) + else: + archive = pyansys.ReadArchive(tetarchivefile) + + # create vtk object + grid = archive.ParseVTK() + + # get cell quality + qual = pyansys.CellQuality(grid) + + # plot cell quality + grid.Plot(scalars=qual, stitle='Cell Minimum Scaled\nJacobian', + rng=[0, 1]) diff --git a/setup.py b/setup.py index 0b449ec8be..3566d8a01d 100644 --- a/setup.py +++ b/setup.py @@ -1,146 +1,146 @@ -""" -Installation file for pyansys -""" -import os -import sys -from io import open as io_open - -from setuptools import setup, Extension -from setuptools.command.build_ext import build_ext as _build_ext - - -class build_ext(_build_ext): - """ build class that includes numpy directory """ - 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(): - """ Check compiler and assign compile arguments accordingly """ - import re - import distutils.ccompiler - comp = distutils.ccompiler.get_default_compiler() - getnext = False - - for a in sys.argv[2:]: - if getnext: - comp = a - getnext = False - continue - # separated by space - if a == '--compiler' or re.search('^-[a-z]*c$', a): - getnext = True - continue - # without space - m = re.search('^--compiler=(.+)', a) - if m is None: - m = re.search('^-[a-z]*c(.+)', a) - if m: - comp = m.group(1) - - return comp - - -# Assign arguments based on compiler -compiler = compilerName() -if compiler == 'unix': - cmp_arg = ['-O3'] -else: - cmp_arg = ['/Ox'] - - -# Get version from version info -__version__ = None -version_file = os.path.join( - os.path.dirname(__file__), - 'pyansys', - '_version.py') -with io_open(version_file, mode='r') as fd: - # execute file from raw string - exec(fd.read()) - - -# Actual setup -setup( - name='pyansys', - packages=['pyansys', 'pyansys.examples'], - - # Version - version=__version__, - - description='Pythonic interface to ANSYS binary files', - long_description=open('README.rst').read(), - - # Author details - author='Alex Kaszynski', - author_email='akascap@gmail.com', - - license='MIT', - classifiers=[ - 'Development Status :: 4 - Beta', - 'Intended Audience :: Science/Research', - 'Topic :: Scientific/Engineering :: Information Analysis', - 'License :: OSI Approved :: MIT License', - 'Programming Language :: Python :: 2.7', - 'Programming Language :: Python :: 3.5', - 'Programming Language :: Python :: 3.6', - ], - - # Website - url='https://github.com/akaszynski/pyansys', - - # Build cython modules - cmdclass={'build_ext': build_ext}, - ext_modules=[Extension("pyansys._parsefull", - ['pyansys/cython/_parsefull.pyx', - 'pyansys/cython/parsefull.c'], - extra_compile_args=cmp_arg, - language='c'), - - Extension("pyansys._parser", - ["pyansys/cython/_parser.pyx"], - extra_compile_args=cmp_arg, - language='c'), - - Extension('pyansys._reader', - ['pyansys/cython/_reader.pyx', - 'pyansys/cython/reader.c'], - extra_compile_args=cmp_arg, - language='c',), - - Extension("pyansys._relaxmidside", - ["pyansys/cython/_relaxmidside.pyx"], - extra_compile_args=cmp_arg, - language='c'), - - # cell quality module - Extension("pyansys._cellqual", - ["pyansys/cython/_cellqual.pyx"], - extra_compile_args=cmp_arg, - language='c'), - - # cell quality module - Extension("pyansys._cellqualfloat", - ["pyansys/cython/_cellqualfloat.pyx"], - extra_compile_args=cmp_arg, - language='c'), - - Extension("pyansys._rstHelper", - ["pyansys/cython/_rstHelper.pyx"], - extra_compile_args=cmp_arg, - language='c'), - - ], - - keywords='vtk ANSYS cdb full rst', - package_data={'pyansys.examples': ['TetBeam.cdb', 'HexBeam.cdb', - 'file.rst', 'file.full']}, - - # Might work with earlier versions - install_requires=['numpy>1.9.3', 'cython>0.23.1', 'vtkInterface>=0.4.0'] - -) +""" +Installation file for pyansys +""" +import os +import sys +from io import open as io_open + +from setuptools import setup, Extension +from setuptools.command.build_ext import build_ext as _build_ext + + +class build_ext(_build_ext): + """ build class that includes numpy directory """ + 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(): + """ Check compiler and assign compile arguments accordingly """ + import re + import distutils.ccompiler + comp = distutils.ccompiler.get_default_compiler() + getnext = False + + for a in sys.argv[2:]: + if getnext: + comp = a + getnext = False + continue + # separated by space + if a == '--compiler' or re.search('^-[a-z]*c$', a): + getnext = True + continue + # without space + m = re.search('^--compiler=(.+)', a) + if m is None: + m = re.search('^-[a-z]*c(.+)', a) + if m: + comp = m.group(1) + + return comp + + +# Assign arguments based on compiler +compiler = compilerName() +if compiler == 'unix': + cmp_arg = ['-O3'] +else: + cmp_arg = ['/Ox'] + + +# Get version from version info +__version__ = None +version_file = os.path.join( + os.path.dirname(__file__), + 'pyansys', + '_version.py') +with io_open(version_file, mode='r') as fd: + # execute file from raw string + exec(fd.read()) + + +# Actual setup +setup( + name='pyansys', + packages=['pyansys', 'pyansys.examples'], + + # Version + version=__version__, + + description='Pythonic interface to ANSYS binary files', + long_description=open('README.rst').read(), + + # Author details + author='Alex Kaszynski', + author_email='akascap@gmail.com', + + license='MIT', + classifiers=[ + 'Development Status :: 4 - Beta', + 'Intended Audience :: Science/Research', + 'Topic :: Scientific/Engineering :: Information Analysis', + 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python :: 2.7', + 'Programming Language :: Python :: 3.5', + 'Programming Language :: Python :: 3.6', + ], + + # Website + url='https://github.com/akaszynski/pyansys', + + # Build cython modules + cmdclass={'build_ext': build_ext}, + ext_modules=[Extension("pyansys._parsefull", + ['pyansys/cython/_parsefull.pyx', + 'pyansys/cython/parsefull.c'], + extra_compile_args=cmp_arg, + language='c'), + + Extension("pyansys._parser", + ["pyansys/cython/_parser.pyx"], + extra_compile_args=cmp_arg, + language='c'), + + Extension('pyansys._reader', + ['pyansys/cython/_reader.pyx', + 'pyansys/cython/reader.c'], + extra_compile_args=cmp_arg, + language='c',), + + Extension("pyansys._relaxmidside", + ["pyansys/cython/_relaxmidside.pyx"], + extra_compile_args=cmp_arg, + language='c'), + + # cell quality module + Extension("pyansys._cellqual", + ["pyansys/cython/_cellqual.pyx"], + extra_compile_args=cmp_arg, + language='c'), + + # cell quality module + Extension("pyansys._cellqualfloat", + ["pyansys/cython/_cellqualfloat.pyx"], + extra_compile_args=cmp_arg, + language='c'), + + Extension("pyansys._rstHelper", + ["pyansys/cython/_rstHelper.pyx"], + extra_compile_args=cmp_arg, + language='c'), + + ], + + keywords='vtk ANSYS cdb full rst', + package_data={'pyansys.examples': ['TetBeam.cdb', 'HexBeam.cdb', + 'file.rst', 'file.full']}, + + # Might work with earlier versions + install_requires=['numpy>1.9.3', 'cython>0.23.1', 'vtkInterface>=0.4.0'] + +) diff --git a/uploadpypi.sh b/uploadpypi.sh new file mode 100755 index 0000000000..bf6699ebd2 --- /dev/null +++ b/uploadpypi.sh @@ -0,0 +1 @@ +twine upload dist/*