# Source code for assembly

```"""
Assemble a series of sparse matrices using two gathering matrices. For
`[F. Cuvelier, C. Japhet and G. Scarella, An efficient way to assemble
finite element matrices in vector languages, 2016]
<http://doi.org/10.1007/s10543-015-0587-4>`_.

⭕ To access the source code, click on the [source] button at the right
side or click on
Dependence may exist. In case of error, check import and install required

"""

import numpy as np
from scipy import sparse as spspa

[docs]def assemble(Ms, Gr, Gc=None):
"""A (not very fast) routine for assembling matrices and
vectors.

:param Ms: The matrices/vectors to be assembled. ``Ms[i]`` refers
to the ith (for example, in element No. i) matrix/vector. The
matrices/vectors need to be sparse matrices/vectors. A sparse
vector is a csc_matrix of shape :math:`(n, 1)`.
:param Gr: The row gathering matrix.
:type Gr: np.array
:param Gc: (default: ``None``) The column gathering matrix. When it
is ``None``, it means we are assembling vectors. Therefore we
will only need the row gathering matrix, ``Gr``.
:type Gc: None, np.array
:return: A csc matrix representing the assembled matrix/vector.

:example:

>>> from crazy_mesh import CrazyMeshGlobalNumbering
>>> from scipy import sparse as spspa
>>> K = 2
>>> N = 3
>>> GM = CrazyMeshGlobalNumbering(K, N)
>>> GM_F = GM.FP
>>> GM_V = GM.VP
>>> Ms = [spspa.random(N**3, 3*(N+1)*N**2, 0.1, 'csc')
...     for _ in range(K**3)] # generate a series of sparse matrices
>>> M = assemble(Ms, GM_V, GM_F) # assemble matrices
>>> M.shape
(216, 756)
>>> Vs = [spspa.random(N**3, 1, 0.5, 'csc')
...     for _ in range(K**3)] # generate a series of sparse vectors
>>> V = assemble(Vs, GM_V)
>>> V.shape
(216, 1)

"""
if Gc is None:  # we are assembling vectors
DEP = int(np.max(Gr) + 1)
ROW = list()
DAT = list()
for i, Vi in enumerate(Ms):
assert Vi.__class__.__name__ == 'csc_matrix' and \
np.shape(Vi)[1] == 1, f"V[{i}] is not a sparse vector."
indices = Vi.indices
data = Vi.data
ROW.extend(Gr[i][indices])
DAT.extend(data)
return spspa.csc_matrix((DAT, ROW, [0, len(ROW)]), shape=(DEP, 1))

else:
DEP = int(np.max(Gr) + 1)
WID = int(np.max(Gc) + 1)
ROW = list()
COL = list()
DAT = list()
for i, Mi in enumerate(Ms):
assert Mi.__class__.__name__ in ('csc_matrix', 'csr_matrix')
indices = Mi.indices
indptr = Mi.indptr
data = Mi.data
nums = np.diff(indptr)
if Mi.__class__.__name__ == 'csc_matrix':
for j, num in enumerate(nums):
idx = indices[indptr[j]:indptr[j + 1]]
ROW.extend(Gr[i][idx])
COL.extend([Gc[i][j], ] * num)
elif Mi.__class__.__name__ == 'csr_matrix':
for j, num in enumerate(nums):
idx = indices[indptr[j]:indptr[j + 1]]
ROW.extend([Gr[i][j], ] * num)
COL.extend(Gc[i][idx])
else:
raise Exception("I can not handle %r." % Mi)
DAT.extend(data)

return spspa.csc_matrix((DAT, (ROW, COL)),shape=(DEP, WID))

if __name__ == '__main__':
import doctest
doctest.testmod()
```