Reputation: 271
build Scipy matrix with diagonal overlays
I would like to overlay several arrays of equal 4x4 size at offset positions creating a large array where the overlapping elements are added together.
[a] 0 0 0
0 [a+b] 0 0
0 0 [b+c] 0
0 0 0 [c+...
This is for the creation of a stiffness matrix for continuous beam analysis. This is easy enough to create with a loop, but is there a better way using Scipy sparse matrix routines? block_diag is close but does not provide a method for overlays.
Members = [[100.0, 1000.0],[100.0, 2000.0],[200.0, 2000.0]]
SSM = np.zeros((2*M+2,2*M+2), dtype = float)
i=0
for mem in Members:
a = SMX(mem)
print a
print "*"*40
SSM[i:i+4,i:i+4] = SSM[i:i+4,i:i+4] + a
print SSM
i +=2
def SMX(member):
""" Prismatic beam member stiffness function
Purpose: Create a member stiffness matrix for a given set of member parameters.
Input: Modulus of elasticity - E
Moment of inertia with respect to the z axis - Iz
Length of the span - L
Returns: 4 x 4 matrix representing the member stiffness with respect to end conditions
"""
L = member[0]
Iz = member[1]
pbms = np.zeros((4,4), dtype = float)
pbms[0,0] = (12.0 * E * Iz)/L**2
pbms[0,1] = (6.0 * E * Iz)/L**2
pbms[0,2] = -pbms[0,0]
pbms[0,3] = pbms[0,1]
pbms[1,0] = pbms[0,1]
pbms[1,1] = (4.0 * E * Iz)/L
pbms[1,2] = -pbms[0,1]
pbms[1,3] = (2.0 * E * Iz)/L
pbms[2,0] = pbms[0,2]
pbms[2,1] = -pbms[0,1]
pbms[2,2] = pbms[0,0]
pbms[2,3] = -pbms[0,1]
pbms[3,0] = pbms[0,1]
pbms[3,1] = pbms[1,3]
pbms[3,2] = -pbms[0,1]
pbms[3,3] = pbms[1,1]
return pbms
The resulting output:
[[ 12000. 6000. -12000. 6000.]
[ 6000. 400000. -6000. 200000.]
[ -12000. -6000. 12000. -6000.]
[ 6000. 200000. -6000. 400000.]]
****************************************
[[ 12000. 6000. -12000. 6000. 0. 0. 0. 0.]
[ 6000. 400000. -6000. 200000. 0. 0. 0. 0.]
[ -12000. -6000. 12000. -6000. 0. 0. 0. 0.]
[ 6000. 200000. -6000. 400000. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0.]
[ 0. 0. 0. 0. 0. 0. 0. 0.]]
[[ 24000. 12000. -24000. 12000.]
[ 12000. 800000. -12000. 400000.]
[ -24000. -12000. 24000. -12000.]
[ 12000. 400000. -12000. 800000.]]
****************************************
[[ 12000. 6000. -12000. 6000. 0. 0. 0.
0.]
[ 6000. 400000. -6000. 200000. 0. 0. 0.
0.]
[ -12000. -6000. 36000. 6000. -24000. 12000. 0.
0.]
[ 6000. 200000. 6000. 1200000. -12000. 400000. 0.
0.]
[ 0. 0. -24000. -12000. 24000. -12000. 0.
0.]
[ 0. 0. 12000. 400000. -12000. 800000. 0.
0.]
[ 0. 0. 0. 0. 0. 0. 0.
0.]
[ 0. 0. 0. 0. 0. 0. 0.
0.]]
[[ 6000. 3000. -6000. 3000.]
[ 3000. 400000. -3000. 200000.]
[ -6000. -3000. 6000. -3000.]
[ 3000. 200000. -3000. 400000.]]
****************************************
[[ 12000. 6000. -12000. 6000. 0. 0. 0.
0.]
[ 6000. 400000. -6000. 200000. 0. 0. 0.
0.]
[ -12000. -6000. 36000. 6000. -24000. 12000. 0.
0.]
[ 6000. 200000. 6000. 1200000. -12000. 400000. 0.
0.]
[ 0. 0. -24000. -12000. 30000. -9000. -6000.
3000.]
[ 0. 0. 12000. 400000. -9000. 1200000. -3000.
200000.]
[ 0. 0. 0. 0. -6000. -3000. 6000.
-3000.]
[ 0. 0. 0. 0. 3000. 200000. -3000.
400000.]]
Just noticed that I left out M = 3 above.
Upvotes: 3
Views: 491
Answers (2)
Reputation: 114841
I don't see the problem with scipy.sparse.block_diag, but I probably don't understand exactly what you mean by "overlay". For example, the output you illustrated in the question can be formed like this:
In [41]: from scipy import sparse
In [42]: a = np.array([[1, 2], [3, 4]])
In [43]: b = np.array([[5, 6], [7, 8]])
In [44]: c = np.array([[9, 10], [11, 12]])
In [45]: m = sparse.block_diag([a, a+b, b+c])
In [46]: m.A
Out[46]:
array([[ 1, 2, 0, 0, 0, 0],
[ 3, 4, 0, 0, 0, 0],
[ 0, 0, 6, 8, 0, 0],
[ 0, 0, 10, 12, 0, 0],
[ 0, 0, 0, 0, 14, 16],
[ 0, 0, 0, 0, 18, 20]])
Edit: After browsing through some material on beam analysis, I think I understand better how you want to overlay the matrices.
For example, suppose there are four 2x2 arrays a, b, c, and d. An example of an "overlay" is the following:
a[0,0] a[0,1] 0 0 0
a[1,0] a[1,1]+b[0,0] b[0,1] 0 0
0 b[1,0] b[1,1]+c[0,0] c[0,1] 0
0 0 c[1,0] c[1,1]+d[0,0] d[0,1]
0 0 0 d[1,0] d[1,1]
One way to form this with block_diag is demonstrated here. It involves an expression with a term for each 2x2 array, which could be rewritten as a loop, so it doesn't really answer your question.
In [92]: a,b,c,d
Out[92]:
(array([[1, 2],
[3, 4]]),
array([[5, 6],
[7, 8]]),
array([[ 9, 10],
[11, 12]]),
array([[13, 14],
[15, 16]]))
In [93]: bd = sparse.block_diag # A convenient alias
In [94]: m = bd([a,0,0,0])+bd([0,b,0,0])+bd([0,0,c,0])+bd([0,0,0,d])
In [95]: m.A
Out[95]:
array([[ 1, 2, 0, 0, 0],
[ 3, 9, 6, 0, 0],
[ 0, 7, 17, 10, 0],
[ 0, 0, 11, 25, 14],
[ 0, 0, 0, 15, 16]])
Upvotes: 1
Reputation: 32521
You could use a coo_matrix. You need to supply three arrays: the data, the row indices and the column indices as per the second example in the link. The key point being:
By default when converting to CSR or CSC format, duplicate (i,j) entries will be summed together. This facilitates efficient construction of finite element matrices and the like. (see example)
Upvotes: 0
Related Questions
- diagonal matrix of a matrix with numpy and scipy
- Diagonal stacking in numpy?
- 2D array with 2D arrays on the diagonal
- Make special diagonal matrix in Numpy
- Efficient way of creating a matrix in numpy by duplicating a given array along a diagonal
- Stacking matrices to make a matrix with sites of parent matrices map as block diagonals
- Stack numpy array onto diagonal
- Construct (N+1)-dimensional diagonal matrix from values in N-dimensional array
- Python code to produce diagonal matrix with 1 above the main diagonal
- create a matrix from array of elements under diagonal in numpy