From d7f711b7fa05c6bb956c377b5def7e7d60305cb3 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Wed, 7 Mar 2018 12:56:54 +0100 Subject: [PATCH] fixed full reader --- .gitignore | 1 + doc/loading_km.rst | 36 ++++---- pyansys/_version.py | 2 +- pyansys/binary_reader.py | 156 +++++++++++++++++++--------------- pyansys/cython/_rstHelper.pyx | 138 +++++++++++------------------- pyansys/examples/examples.py | 36 +++----- 6 files changed, 172 insertions(+), 197 deletions(-) diff --git a/.gitignore b/.gitignore index a9b33db233..9c2d3ececa 100644 --- a/.gitignore +++ b/.gitignore @@ -6,6 +6,7 @@ *.cpp *.so *.o +*.cache # OS generated files # ###################### diff --git a/doc/loading_km.rst b/doc/loading_km.rst index 08e9793dbb..9a65c852f5 100644 --- a/doc/loading_km.rst +++ b/doc/loading_km.rst @@ -7,30 +7,39 @@ Reading a Full File ------------------- This example reads in the mass and stiffness matrices associated with the above example. ``LoadKM`` sorts degrees of freedom such that the nodes are ordered from minimum to maximum, and each degree of freedom (i.e. X, Y, Z), are sorted within each node. The matrices ``k`` and ``m`` are sparse by default, but if ``scipy`` is not installed, or if the optional parameter ``as_sparse=False`` then they will be full numpy arrays. -By default ``LoadKM`` outputs the upper triangle of both matrices, to output the full matrix, set ``utri=False``. Additionally, the constrained nodes of the analysis can be identified by accessing ``fobj.const`` where the constrained degrees of freedom are True and all others are False. This corresponds to the degrees of reference in ``dof_ref``. +By default ``LoadKM`` outputs the upper triangle of both matrices. The constrained nodes of the analysis can be identified by accessing ``fobj.const`` where the constrained degrees of freedom are True and all others are False. This corresponds to the degrees of reference in ``dof_ref``. + +By default dof_ref is unsorted. To sort these values, set ``sort==True``. It is enabled for this example to allow for plotting of the values later on. .. code:: python # Load pyansys import pyansys + from pyansys import examples # Create result reader object and read in full file - full = pyansys.FullReader(pyansys.examples.fullfile) - dof_ref, k, m = full.LoadKM(utri=False) # return the full matrix + full = pyansys.FullReader(examples.fullfile) + dof_ref, k, m = full.LoadKM(sort=True) + -If you have ``scipy`` installed, you can solve solve for the natural frequencies and mode shapes of a system. Realize that constrained degrees of freedom must be removed from the ``k`` and ``m`` matrices for the correct solution. +ANSYS only stores the upper triangular matrix in the full file. To make the full matrix: .. code:: python - import numpy as np - # remove the constrained degrees of freedom - # NOTE: There are more efficient way to remove these indices - free = np.logical_not(full.const).nonzero()[0] - k = k[free][:, free] - m = m[free][:, free] + k += sparse.triu(k, 1).T + m += sparse.triu(m, 1).T + +If you have ``scipy`` installed, you can solve solve for the natural frequencies and mode shapes of a system. + +.. code:: python + import numpy as np from scipy.sparse import linalg + # condition the k matrix + # to avoid getting the "Factor is exactly singular" error + k += sparse.diags(np.random.random(k.shape[0])/1E20, shape=k.shape) + # Solve w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) @@ -46,6 +55,7 @@ If you have ``scipy`` installed, you can solve solve for the natural frequencies First four natural frequencies 1283.200 Hz 1283.200 Hz + 5781.975 Hz 6919.399 Hz @@ -59,11 +69,7 @@ You can also plot the mode shape of this finite element model. Since the constr 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(dof_ref.shape[0]) - full_mode_shape[np.logical_not(full.const)] = mode_shape + full_mode_shape = v[:, 3] # x, y, z displacement for each node # reshape and compute the normalized displacement disp = full_mode_shape.reshape((-1, 3)) diff --git a/pyansys/_version.py b/pyansys/_version.py index 33489c8d95..70132a3489 100644 --- a/pyansys/_version.py +++ b/pyansys/_version.py @@ -1,5 +1,5 @@ # major, minor, patch -version_info = 0, 23, 0 +version_info = 0, 24, 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 a0ed63a554..7a796a292e 100755 --- a/pyansys/binary_reader.py +++ b/pyansys/binary_reader.py @@ -108,7 +108,7 @@ def __init__(self, filename): raise Exception( "Unable to read an unsymmetric mass/stiffness matrix.") - def LoadKM(self, as_sparse=True, utri=True): + def LoadKM(self, as_sparse=True, sort=False): """ Load and construct mass and stiffness matrices from an ANSYS full file. @@ -118,21 +118,21 @@ def LoadKM(self, as_sparse=True, utri=True): Outputs the mass and stiffness matrices as scipy csc sparse arrays when True by default. - utri : bool, optional - Outputs only the upper triangle of both the mass and stiffness - arrays. + sort : bool, optional + Rearranges the k and m matrices such that the rows correspond to + to the sorted rows and columns in dor_ref. Also sorts dor_ref. Returns ------- dof_ref : (n x 2) np.int32 array This array contains the node and degree corresponding to each row - and column in the mass and stiffness matrices. When the sort - parameter is set to True this array will be sorted by node number - first and then by the degree of freedom. In a 3 DOF analysis the - intergers will correspond to: + and column in the mass and stiffness matrices. In a 3 DOF + analysis the dof intergers will correspond to: 0 - x 1 - y 2 - z + Sort these values by node number and DOF by enabling the sort + parameter k : (n x n) np.float or scipy.csc array Stiffness array @@ -140,7 +140,6 @@ def LoadKM(self, as_sparse=True, utri=True): m : (n x n) np.float or scipy.csc array Mass array """ - # check file exists if not os.path.isfile(self.filename): raise Exception('%s not found' % self.filename) @@ -157,84 +156,89 @@ def LoadKM(self, as_sparse=True, utri=True): ntermK = self.header[9] # number of terms in stiffness matrix ptrSTF = self.header[19] # Location of stiffness matrix ptrMAS = self.header[27] # Location in file to mass matrix - nNodes = self.header[33] # Number of nodes considered by assembly + # nNodes = self.header[33] # Number of nodes considered by assembly ntermM = self.header[34] # number of terms in mass matrix ptrDOF = self.header[36] # pointer to DOF info - # get details for reading the mass and stiffness arrays - node_info = _rstHelper.FullNodeInfo(self.filename, ptrDOF, nNodes, - neqn) + # DOF information + ptrDOF = self.header[36] # pointer to DOF info + with open(self.filename, 'rb') as f: + ReadTable(f, skip=True) # standard header + ReadTable(f, skip=True) # full header + ReadTable(f, skip=True) # number of degrees of freedom + neqv = ReadTable(f) # Nodal equivalence table + + f.seek(ptrDOF*4) + ndof = ReadTable(f) + const = ReadTable(f) - nref, dref, index_arr, const, ndof = node_info - dof_ref = np.vstack((nref, dref)).T # stack these references + # dof_ref = np.vstack((ndof, neqv)).T # stack these references + dof_ref = [ndof, neqv] # Read k and m blocks (see help(ReadArray) for block description) if ntermK: - k_block = _rstHelper.ReadArray(self.filename, ptrSTF, ntermK, neqn, - index_arr) - k_diag = k_block[3] - k_data_diag = k_block[4] + krow, kcol, kdata = _rstHelper.ReadArray(self.filename, + ptrSTF, + ntermK, + neqn, + const) else: warnings.warn('Missing stiffness matrix') - k_block = None + kdata = None if ntermM: - m_block = _rstHelper.ReadArray(self.filename, ptrMAS, ntermM, neqn, - index_arr) - m_diag = m_block[3] - m_data_diag = m_block[4] + mrow, mcol, mdata = _rstHelper.ReadArray(self.filename, + ptrMAS, + ntermM, + neqn, + const) else: warnings.warn('Missing mass matrix') - m_block = None - - self.m_block = m_block - self.k_block = k_block - - # assemble data - if utri: - if k_block: - # stiffness matrix - krow = np.hstack((k_block[1], k_diag)) # row and diag - kcol = np.hstack((k_block[0], k_diag)) # col and diag - kdata = np.hstack((k_block[2], k_data_diag)) # data and diag - - if m_block: - # mass matrix - mrow = np.hstack((m_block[1], m_diag)) # row and diag - mcol = np.hstack((m_block[0], m_diag)) # col and diag - mdata = np.hstack((m_block[2], m_data_diag)) # data and diag - + mdata = None + + # remove constrained entries + if np.any(const < 0): + if kdata is not None: + remove = np.nonzero(const < 0)[0] + mask = ~np.logical_or(np.in1d(krow, remove), np.in1d(kcol, remove)) + krow = krow[mask] + kcol = kcol[mask] + kdata = kdata[mask] + + if mdata is not None: + mask = ~np.logical_or(np.in1d(mrow, remove), np.in1d(mcol, remove)) + mrow = mrow[mask] + mcol = mcol[mask] + mdata = mdata[mask] + + dof_ref, index, nref, dref = _rstHelper.SortNodalEqlv(neqn, neqv, ndof) + if sort: # make sorting the same as ANSYS rdfull would output + # resort to make in upper triangle + krow = index[krow] + kcol = index[kcol] + krow, kcol = np.sort(np.vstack((krow, kcol)), 0) + + mrow = index[mrow] + mcol = index[mcol] + mrow, mcol = np.sort(np.vstack((mrow, mcol)), 0) else: - if k_block: - # stiffness matrix - krow = np.hstack((k_block[0], k_block[1], k_diag)) - kcol = np.hstack((k_block[1], k_block[0], k_diag)) - kdata = np.hstack((k_block[2], k_block[2], k_data_diag)) - - if m_block: - # mass matrix - mrow = np.hstack((m_block[0], m_block[1], m_diag)) - mcol = np.hstack((m_block[1], m_block[0], m_diag)) - mdata = np.hstack((m_block[2], m_block[2], m_data_diag)) + dof_ref = np.vstack((nref, dref)).T # store data for later reference - if k_block: + if kdata is not None: self.krow = krow self.kcol = kcol self.kdata = kdata - if m_block: + if mdata is not None: self.mrow = mrow self.mcol = mcol self.mdata = mdata - # number of dimentions - ndim = nref.size - # output as a sparse matrix if as_sparse: - if k_block: - k = coo_matrix((ndim,) * 2) + if kdata is not None: + k = coo_matrix((neqn,) * 2) k.data = kdata # data has to be set first k.row = krow k.col = kcol @@ -244,8 +248,8 @@ def LoadKM(self, as_sparse=True, utri=True): else: k = None - if m_block: - m = coo_matrix((ndim,) * 2) + if mdata is not None: + m = coo_matrix((neqn,) * 2) m.data = mdata m.row = mrow m.col = mcol @@ -256,14 +260,14 @@ def LoadKM(self, as_sparse=True, utri=True): m = None else: - if k_block: - k = np.zeros((ndim,) * 2) + if kdata is not None: + k = np.zeros((neqn,) * 2) k[krow, kcol] = kdata else: k = None - if m_block: - m = np.zeros((ndim,) * 2) + if mdata is not None: + m = np.zeros((neqn,) * 2) m[mrow, mcol] = mdata else: m = None @@ -376,8 +380,7 @@ def AddCyclicProperties(self): vtkappend.AddInputData(sector) - # Combine meshes and add VTK_Utilities functions - # vtkappend.MergePointsOn() + # Combine meshes vtkappend.Update() self.rotor = vtkInterface.UnstructuredGrid(vtkappend.GetOutput()) @@ -1454,3 +1457,18 @@ def delete_row_csc(mat, i): mat.indptr[i:] -= n mat.indptr = mat.indptr[:-1] mat._shape = (mat._shape[0] - 1, mat._shape[1]) + + +def ReadTable(f, dtype='i', skip=False): + """ read fortran style table """ + tablesize = np.fromfile(f, 'i', 1)[0] + f.seek(4, 1) # skip padding + if skip: + f.seek((tablesize + 1)*4, 1) + return + else: + if dtype == 'double': + tablesize //= 2 + table = np.fromfile(f, dtype, tablesize) + f.seek(4, 1) # skip padding + return table diff --git a/pyansys/cython/_rstHelper.pyx b/pyansys/cython/_rstHelper.pyx index 531ab6e42c..111a378301 100755 --- a/pyansys/cython/_rstHelper.pyx +++ b/pyansys/cython/_rstHelper.pyx @@ -347,7 +347,7 @@ def LoadStressDouble(filename, py_table_index, int64_t [::1] ele_ind_table, fclose(cfile) -def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): +def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] const): """ Reads stiffness or mass matrices from ANSYS fortran files @@ -365,8 +365,8 @@ def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): neqn : int Number of equations - index_arr : numpy int array - Indexing array + const : numpy int array + If DOF is fixed Returns ------- @@ -390,12 +390,10 @@ def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): ----- Function signature ReadArray(filename, int ptrSTF, int nread, int nterm, int neqn, - int [::1] index_arr) - + int [::1] const) """ cdef int i, j, k, ind - cdef bytes buf with open(filename, "rb") as f: f.seek(ptr*4) @@ -404,64 +402,55 @@ def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): cdef char * p = buf # python to c character array pointer cdef int loc = 0 # location long buffer - # Half of the data - cdef int ntermadj = nterm - neqn # upper triangle (sorta) - cdef int [::1] krow = np.empty(ntermadj, np.int32) - cdef int [::1] kcol = np.empty(ntermadj, np.int32) - cdef double [::1] kdata = np.empty(ntermadj) - + # Number of terms is the number of terms stored in the upper triangle + cdef int [::1] krow = np.empty(nterm, np.int32) + cdef int [::1] kcol = np.empty(nterm, np.int32) + cdef double [::1] kdata = np.empty(nterm) cdef int [::1] kdiag = np.empty(neqn, np.int32) cdef double [::1] kdata_diag = np.empty(neqn) cdef int c = 0 # index counter cdef int d = 0 # data counter + cdef int c_diag = 0 # diag counter cdef int row, col, nitems, intval cdef double val + for i in range(neqn): - col = index_arr[i] + col = i # number of items to read nitems = GetInt(&p[loc]); loc += 4 loc += 4 - - # Indices: read in all but diagional term - for j in range(nitems - 1): + + # Read indices + for j in range(nitems): # get row number row = GetInt(&p[loc]) - 1; loc += 4 # convert to c indexing - - # store upper triangle - if index_arr[row] > col: - krow[c] = index_arr[row] + + if row < col: + krow[c] = row kcol[c] = col else: krow[c] = col - kcol[c] = index_arr[row] - + kcol[c] = row c += 1 - - # store diagional term - row = GetInt(&p[loc]) - 1; loc += 4 - kdiag[i] = index_arr[row] + loc += 12 - - # Data: read in all but diagional term - for j in range(nitems - 1): + + # Read data + for j in range(nitems): # Store data kdata[d] = GetDouble(&p[loc]); loc += 8 d += 1 - - # store diagional data term - kdata_diag[i] = GetDouble(&p[loc]); loc += 8 - + # seek past end of data loc += 4 - return np.asarray(krow), np.asarray(kcol), np.asarray(kdata), \ - np.asarray(kdiag), np.asarray(kdata_diag) - - -def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn): + return np.asarray(krow)[:c], np.asarray(kcol)[:c], np.asarray(kdata)[:c] + + +def SortNodalEqlv(int neqn, int [::1] neqv, int [::1] ndof): """ Reads in full file details required for the assembly of the mass and stiffness matrices. @@ -471,57 +460,28 @@ def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn): Parameters ---------- - filename : string - Full file filename + neqn : int + Number of equations in full file. - ptrDOF : int - Location of the DOF block in the full file - - nNodes : - Number of nodes in full file + ndof : int [::1] + Degrees of freedom for each node. + + neqv : int [::1] + Nodal equivalance array. - neqn : - Number of equations in full file - - Returns ------- - nref : numpy np.int32 array - Sorted nodal reference array - - dref: numpy np.int32 array + dof_ref: numpy np.int32 array Sorted degree of freedom reference array. index_arr : numpy np.int32 array Index array to sort rows and columns. - - const : numpy np.int32 array - Negative if a node's DOF is constrained - - ndof : numpy np.int32 array - Number of degrees of freedom for each node in nref - """ + """ cdef int i, j, ind - cdef int [::1] neqv = np.empty(nNodes, np.int32) - cdef int [::1] ndof = np.empty(nNodes, np.int32) - cdef int [::1] const = np.empty(neqn, np.int32) - - with open(filename, "rb") as f: - # nodal equivalency - f.seek((212 + 2)*4) - neqv = np.fromfile(f, 'i', nNodes) - - # Total DOFs - f.seek((ptrDOF + 2)*4) - ndof = np.fromfile(f, 'i', nNodes) - - # Constrained DOFs - f.seek((ptrDOF + 5 + nNodes)*4) - const = np.fromfile(f, 'i', neqn) - # create sorting array + cdef int nNodes = ndof.size cdef int [::1] cumdof = np.empty(nNodes, np.int32) cdef int csum = 0 for i in range(nNodes): @@ -539,15 +499,15 @@ def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn): nref[c] = val dref[c] = j c += 1 - + # sort nodal equivalance array cdef int [::1] sidx = np.argsort(neqv).astype(np.int32) cdef int [::1] ndof_sort = np.empty(nNodes, np.int32) for i in range(nNodes): ndof_sort[i] = ndof[sidx[i]] - + cdef int d = 0 - # create an index array. this tells the array readers down the line where + # create an index array. this tells the array readers where # to place each row and col when it's sorted cdef int [::1] index_arr = np.empty(neqn, np.int32) for i in range(nNodes): @@ -557,18 +517,16 @@ def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn): s_neqv_dof[d] = c + j index_arr[c + j] = d d += 1 - + # sort node and dof references - cdef int [::1] nref_sort = np.empty(neqn, np.int32) - cdef int [::1] dref_sort = np.empty(neqn, np.int32) - cdef int [::1] const_sort = np.empty(neqn, np.int32) + cdef int [:, ::1] dof_ref = np.empty((neqn, 2), np.int32) for i in range(neqn): - nref_sort[i] = nref[s_neqv_dof[i]] - dref_sort[i] = dref[s_neqv_dof[i]] - 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) + ind = s_neqv_dof[i] + dof_ref[i, 0] = nref[ind] + dof_ref[i, 1] = dref[ind] + + return np.asarray(dof_ref), np.asarray(index_arr), np.asarray(nref), \ + np.asarray(dref) def ComputePrincipalStress(float [:, ::1] stress): diff --git a/pyansys/examples/examples.py b/pyansys/examples/examples.py index eb41ca4c0e..d7056bc069 100755 --- a/pyansys/examples/examples.py +++ b/pyansys/examples/examples.py @@ -91,7 +91,7 @@ def LoadKM(): # Create file reader object fobj = pyansys.FullReader(fullfile) - dofref, k, m = fobj.LoadKM(utri=False) + dofref, k, m = fobj.LoadKM() # print results ndim = k.shape[0] @@ -102,13 +102,13 @@ def LoadKM(): # compute natural frequencies if installed try: from scipy.sparse import linalg + from scipy import sparse 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] + k += sparse.triu(k, 1).T + m += sparse.triu(m, 1).T + k += sparse.diags(np.random.random(k.shape[0])/1E20, shape=k.shape) # Solve w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) @@ -126,18 +126,19 @@ def SolveKM(): Loads and solves a mass and stiffness matrix from an ansys full file """ from scipy.sparse import linalg + from scipy import sparse # load the mass and stiffness matrices full = pyansys.FullReader(pyansys.examples.fullfile) - dofref, k, m = full.LoadKM(utri=False) + dofref, k, m = full.LoadKM(sort=True) - # removing constrained nodes (this is not efficient) - free = np.logical_not(full.const).nonzero()[0] - k = k[free][:, free] - m = m[free][:, free] + # make symmetric + k += sparse.triu(k, 1).T + m += sparse.triu(m, 1).T + k += sparse.diags(np.random.random(k.shape[0])/1E20, shape=k.shape) # Solve - w, v = linalg.eigsh(k, k=20, M=m, sigma=10000) + w, v = linalg.eigsh(k, k=20, M=m, sigma=1000) # System natural frequencies f = (np.real(w))**0.5 / (2 * np.pi) @@ -145,12 +146,8 @@ def SolveKM(): # %% 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 - + full_mode_shape = v[:, 3] # x, y, z displacement for each node + # reshape and compute the normalized displacement disp = full_mode_shape.reshape((-1, 3)) n = (disp * disp).sum(1)**0.5 @@ -160,9 +157,6 @@ def SolveKM(): archive = pyansys.ReadArchive(pyansys.examples.hexarchivefile) grid = archive.ParseVTK() - # plot the normalized displacement - # grid.Plot(scalars=n) - # Fancy plot the displacement plobj = vtkInterface.PlotClass() @@ -170,13 +164,11 @@ def SolveKM(): 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'):