From d58e2a41f378e0b4bb21754876160548f660d116 Mon Sep 17 00:00:00 2001 From: Alex Kaszynski Date: Sun, 25 Feb 2018 19:55:35 +0100 Subject: [PATCH] allows missing mass or stiffness matricies --- .gitignore | 2 + pyansys/_version.py | 2 +- pyansys/archive_reader.py | 26 +++++- pyansys/binary_reader.py | 152 +++++++++++++++++-------------- pyansys/cython/_parser.pyx | 12 ++- pyansys/cython/_relaxmidside.pyx | 10 +- pyansys/cython/_rstHelper.pyx | 13 +-- 7 files changed, 123 insertions(+), 94 deletions(-) diff --git a/.gitignore b/.gitignore index 44f7421244..a9b33db233 100644 --- a/.gitignore +++ b/.gitignore @@ -25,3 +25,5 @@ pyansys/Interface.py # Testing Testing/ +UnitTesting/ +TODO \ No newline at end of file diff --git a/pyansys/_version.py b/pyansys/_version.py index dac76457dc..d53231e73d 100644 --- a/pyansys/_version.py +++ b/pyansys/_version.py @@ -1,5 +1,5 @@ # major, minor, patch -version_info = 0, 22, 5 +version_info = 0, 22, 6 # Nice string for the version diff --git a/pyansys/archive_reader.py b/pyansys/archive_reader.py index b7e9739cf9..20aee1acef 100755 --- a/pyansys/archive_reader.py +++ b/pyansys/archive_reader.py @@ -94,7 +94,6 @@ def ParseVTK(self, force_linear=False, allowable_types=None): cells[cells > numref.max()] = -1 # Check for missing midside nodes - possible_merged = False if force_linear or np.all(cells != -1): nodes = self.raw['nodes'][:, :3].copy() nnum = self.raw['nnum'] @@ -116,9 +115,10 @@ def ParseVTK(self, force_linear=False, allowable_types=None): nodes[nnodes:] = temp_nodes[nnodes:] # merge nodes - from FEMORPH import utilities new_nodes = temp_nodes[nnodes:] - unique_nodes, idxA, idxB = utilities.UniqueRows(new_nodes, return_index=True, return_inverse=True) + unique_nodes, idxA, idxB = UniqueRows(new_nodes, + return_index=True, + return_inverse=True) # rewrite node numbers cells[mask] = idxB + maxnum @@ -217,3 +217,23 @@ def SaveAsVTK(self, filename): """ # place holder print('Run grid = archive.ParseVTK() and then\ngrid.WriteGrid(filename, binary)') + + +def UniqueRows(a, return_index=True, return_inverse=False): + """ Returns unique rows of a and indices of those rows """ + if not a.flags.c_contiguous: + a = np.ascontiguousarray(a) + + b = a.view(np.dtype((np.void, a.dtype.itemsize * a.shape[1]))) + _, idx, idx2 = np.unique(b, True, True) + + if return_index: + if return_inverse: + return a[idx], idx, idx2 + else: + return a[idx], idx + else: + if return_inverse: + return a[idx], idx2 + else: + return a[idx] diff --git a/pyansys/binary_reader.py b/pyansys/binary_reader.py index 040e684c11..49061d065f 100755 --- a/pyansys/binary_reader.py +++ b/pyansys/binary_reader.py @@ -165,54 +165,63 @@ def LoadKM(self, as_sparse=True, utri=True): dof_ref = np.vstack((nref, dref)).T # stack these references # Read k and m blocks (see help(ReadArray) for block description) - k_block = _rstHelper.ReadArray(self.filename, ptrSTF, ntermK, neqn, - index_arr) - - m_block = _rstHelper.ReadArray(self.filename, ptrMAS, ntermM, neqn, - index_arr) - k_diag = k_block[3] - k_data_diag = k_block[4] - - m_diag = m_block[3] - m_data_diag = m_block[4] + if ntermK: + k_block = _rstHelper.ReadArray(self.filename, ptrSTF, ntermK, neqn, + index_arr) + k_diag = k_block[3] + k_data_diag = k_block[4] + else: + warnings.warn('Missing stiffness matrix') + k_block = 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] + else: + warnings.warn('Missing mass matrix') + m_block = None self.m_block = m_block self.k_block = k_block # assemble data if utri: - # 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 - - # 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 + 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 else: - # stiffness matrix - krow = np.hstack( - (k_block[0], k_block[1], k_diag)) # row, col and diag - kcol = np.hstack((k_block[1], k_block[0], k_diag)) # col and diag - kdata = np.hstack( - (k_block[2], k_block[2], k_data_diag)) # data and diag - - # mass matrix - mrow = np.hstack( - (m_block[0], m_block[1], m_diag)) # row, col and diag - mcol = np.hstack((m_block[1], m_block[0], m_diag)) # col and diag - mdata = np.hstack( - (m_block[2], m_block[2], m_data_diag)) # data and diag + 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)) # store data for later reference - self.krow = krow - self.kcol = kcol - self.kdata = kdata - self.mrow = mrow - self.mcol = mcol - self.mdata = mdata + if k_block: + self.krow = krow + self.kcol = kcol + self.kdata = kdata + if m_block: + self.mrow = mrow + self.mcol = mcol + self.mdata = mdata # number of dimentions ndim = nref.size @@ -220,28 +229,40 @@ def LoadKM(self, as_sparse=True, utri=True): # output as a sparse matrix if as_sparse: - k = coo_matrix((ndim,) * 2) - k.data = kdata # data has to be set first - k.row = krow - k.col = kcol + if k_block: + k = coo_matrix((ndim,) * 2) + k.data = kdata # data has to be set first + k.row = krow + k.col = kcol - # convert to csc matrix (generally faster for sparse solvers) - k = csc_matrix(k) + # convert to csc matrix (generally faster for sparse solvers) + k = csc_matrix(k) + else: + k = None - m = coo_matrix((ndim,) * 2) - m.data = mdata - m.row = mrow - m.col = mcol + if m_block: + m = coo_matrix((ndim,) * 2) + m.data = mdata + m.row = mrow + m.col = mcol - # convert to csc matrix (generally faster for sparse solvers) - m = csc_matrix(m) + # convert to csc matrix (generally faster for sparse solvers) + m = csc_matrix(m) + else: + m = None else: - k = np.zeros((ndim,) * 2) - k[krow, kcol] = kdata - - m = np.zeros((ndim,) * 2) - m[mrow, mcol] = mdata + if k_block: + k = np.zeros((ndim,) * 2) + k[krow, kcol] = kdata + else: + k = None + + if m_block: + m = np.zeros((ndim,) * 2) + m[mrow, mcol] = mdata + else: + m = None # store if constrained and number of degrees of freedom per node self.const = const < 0 @@ -268,22 +289,12 @@ class ResultReader(object): """ - def __init__(self, filename, logger=False, load_geometry=True): + def __init__(self, filename, load_geometry=True): """ Loads basic result information from result file and initializes result object. """ - - # set logging level depending on settings - # if logger: - # logging.basicConfig(level=logging.DEBUG) - # else: - # logging.basicConfig(level=logging.CRITICAL) - - # store logger pointer - # log = logging.getLogger(__name__) - # Store filename result header items self.filename = filename self.resultheader = GetResultInfo(filename) @@ -307,13 +318,17 @@ def __init__(self, filename, logger=False, load_geometry=True): # store geometry for later retrival if load_geometry: self.StoreGeometry() - # import pdb; pdb.set_trace() + if self.resultheader['nSector'] > 1 and load_geometry: self.iscyclic = True # Add cyclic properties self.AddCyclicProperties() + def Plot(self): + """ plots result geometry """ + self.grid.Plot() + def AddCyclicProperties(self): """ Adds cyclic properties to result object """ # ansys's duplicate sector contains nodes from the second half of the @@ -831,8 +846,7 @@ def StoreGeometry(self): # store the reference array cell_type = ['45', '95', '185', '186', '92', '187', '154'] - # result = _parser.Parse(self.geometry, True, cell_type) # force_linear - result = _parser.Parse(self.geometry, False, cell_type) # force_linear + result = _parser.Parse(self.geometry, False, cell_type) cells, offset, cell_type, self.numref, _, _, _ = result # catch -1 diff --git a/pyansys/cython/_parser.pyx b/pyansys/cython/_parser.pyx index 14af4614f1..160a1aef91 100755 --- a/pyansys/cython/_parser.pyx +++ b/pyansys/cython/_parser.pyx @@ -41,10 +41,10 @@ typeB[1] = 187 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): + int64_t [::1] cells, uint8 [::1] cell_type, + int64_t [::1] numref, int [:, ::1] elem, int i, int lin): """ - Stores surface triagle vtk cell. Element may be quadradic or linear + Stores surface triangle vtk cell. Element may be quadradic or linear """ # Populate offset array offset[ecount[0]] = ccount[0] @@ -79,10 +79,12 @@ cdef inline void StoreSurfTri(int64_t [::1] offset, int64_t *ecount, int64_t *cc cdef inline void StoreSurfQuad(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): + int64_t [::1] cells, uint8 [::1] cell_type, + int64_t [::1] numref, int [:, ::1] elem, int i, int lin): """ Stores surface quad in vtk cell array. Element may be quadradic or linear + + """ # Populate offset array offset[ecount[0]] = ccount[0] diff --git a/pyansys/cython/_relaxmidside.pyx b/pyansys/cython/_relaxmidside.pyx index 3e38c6ffa2..340ff2ac12 100755 --- a/pyansys/cython/_relaxmidside.pyx +++ b/pyansys/cython/_relaxmidside.pyx @@ -189,24 +189,18 @@ cdef inline void RelaxMid_Hex(int64_t [::1] cellarr, int c, double [:, ::1] pts, def ResetMidside(int64_t [::1] cellarr, double [:, ::1] pts): """ - SIGNATURE (int64_t [::1] cellarr, double [:, ::1] pts) - - DESCRIPTION Resets positions of midside nodes to directly between edge nodes. - INPUTS + Parameters + ---------- cellarr (int64_t [::1]) VTK formatted cell array pts (double [:, ::1]) 3D double point array. - - - OUTPUTS - NONE """ diff --git a/pyansys/cython/_rstHelper.pyx b/pyansys/cython/_rstHelper.pyx index 26501b1428..f528b69678 100755 --- a/pyansys/cython/_rstHelper.pyx +++ b/pyansys/cython/_rstHelper.pyx @@ -282,9 +282,8 @@ def LoadStressDouble(filename, py_table_index, int64_t [::1] ele_ind_table, def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): """ - Reads stiffness or mass matrices from ANSYS fortran files - + Parameters ---------- filename : string @@ -301,7 +300,7 @@ def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): index_arr : numpy int array Indexing array - + Returns ------- rows : numpy int32 array @@ -395,21 +394,19 @@ def ReadArray(filename, int ptr, int nterm, int neqn, int [::1] index_arr): np.asarray(kdiag), np.asarray(kdata_diag) -# consider adding a sorting flag to disable sorting def FullNodeInfo(filename, int ptrDOF, int nNodes, int neqn): """ - Reads in full file details required for the assembly of the mass and stiffness matrices. - + The reference arrays are sorted by default, though this increases the bandwidth of the mass and stiffness matrices. - + Parameters ---------- filename : string Full file filename - + ptrDOF : int Location of the DOF block in the full file