Skip to content

Commit

Permalink
allows missing mass or stiffness matricies
Browse files Browse the repository at this point in the history
  • Loading branch information
@akaszynski
akaszynski committed Feb 25, 2018
1 parent 8817d49 commit d58e2a4
Show file tree
Hide file tree
Showing 7 changed files with 123 additions and 94 deletions.
2 changes: 2 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,3 +25,5 @@ pyansys/Interface.py

# Testing
Testing/
UnitTesting/
TODO
2 changes: 1 addition & 1 deletion pyansys/_version.py
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# major, minor, patch
version_info = 0, 22, 5
version_info = 0, 22, 6


# Nice string for the version
Expand Down
26 changes: 23 additions & 3 deletions pyansys/archive_reader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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']
Expand All @@ -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
Expand Down Expand Up @@ -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]
152 changes: 83 additions & 69 deletions pyansys/binary_reader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -165,83 +165,104 @@ 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

# 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
Expand All @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down
12 changes: 7 additions & 5 deletions pyansys/cython/_parser.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down Expand Up @@ -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]
Expand Down
10 changes: 2 additions & 8 deletions pyansys/cython/_relaxmidside.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
"""

Expand Down
13 changes: 5 additions & 8 deletions pyansys/cython/_rstHelper.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit d58e2a4

Please sign in to comment.