Block matrix reduce system for finite element

test_systemsolution.py

import numpy as np
import sympy as sp
from numpy import linalg as la

"""
https://mathoverflow.net/questions/425067/block-matrix-reduce-system-for-finite-element-method
"""

def random_sym_matrix(side):
	F = np.random.rand(side, side)
	return (F+F.T)/2

tolerance = 1e-9
ntests = 1000
testallall = True
for i in range(ntests):

	# a, b, c, d, e = 4, 4, 4, 4, 4
	a = int(1+100*np.random.rand())
	b = int(1+100*np.random.rand())
	c = int(1+100*np.random.rand())
	d = int(1+100*np.random.rand())
	e = int(1+100*np.random.rand())
	# a = b = c = d = e = 1
	n = (a+b+c+d+e)
	A = random_sym_matrix(a+b)
	B = random_sym_matrix(b+c)
	C = random_sym_matrix(c+d)
	D = random_sym_matrix(d+e)
	Aexp = np.zeros((n, n), dtype="float64")
	Bexp = np.zeros((n, n), dtype="float64")
	Cexp = np.zeros((n, n), dtype="float64")
	Dexp = np.zeros((n, n), dtype="float64")
	Eexp = np.zeros((n, n), dtype="float64")
	Aexp[:a+b,:a+b] = A
	Bexp[a:a+b+c,a:a+b+c] = B
	Cexp[a+b:a+b+c+d,a+b:a+b+c+d] = C
	Dexp[a+b+c:,a+b+c:] = D

	F1 = np.random.rand(a)
	F2 = np.random.rand(b)
	F3 = np.random.rand(c)
	F3 = np.zeros(c)
	F4 = np.random.rand(d)
	F5 = np.random.rand(e)
	Kexp = Aexp + Bexp + Cexp + Dexp
	Fexp = np.zeros(n)
	Fexp[:a] += F1
	Fexp[a:a+b] += F2
	Fexp[a+b:a+b+c]+=F3
	Fexp[a+b+c:a+b+c+d]+=F4
	Fexp[a+b+c+d:]+=F5

	Uexp = la.solve(Kexp, Fexp)
	U1e = Uexp[:a]
	U2e = Uexp[a:a+b]
	U3e = Uexp[a+b:a+b+c]
	U4e = Uexp[a+b+c:a+b+c+d]
	U5e = Uexp[a+b+c+d:]



	M11 = B[:b,:b]
	M12 = B[:b,b:]
	M22 = B[b:,b:] + C[:c,:c]
	M23 = C[:c,c:]
	M33 = C[c:,c:]
	M22inv = la.inv(M22)

	W = np.zeros((b+d, b+d), dtype="float64")
	W[:b, :b] += M11 - M12 @ M22inv @ M12.T
	W[:b, b:] -= M12 @ M22inv @ M23
	W[b:, :b] -= M23.T @ M22inv @ M12.T
	W[b:, b:] += M33 - M23.T @ M22inv @ M23

	diff = np.abs(W-W.T)
	isSymm = np.all(diff < tolerance)
	# print("Is W symmetric ? ", isSymm)

	m = a+b+d+e
	Ared = np.zeros((m, m))
	Dred = np.zeros((m, m))
	Wred = np.zeros((m, m))
	Fred = np.zeros(m)
	Ared[:a+b, :a+b] += A
	Wred[a:a+b+d, a:a+b+d] += W
	Dred[a+b:, a+b:] += D
	Kred = Ared + Wred + Dred
	Fred[:a] += F1
	Fred[a:a+b] += F2
	Fred[a+b:a+b+d] += F4
	Fred[a+b+d:] += F5

	Ured = la.solve(Kred, Fred)
	U1r = Ured[:a]
	U2r = Ured[a:a+b]
	U4r = Ured[a+b:a+b+d]
	U5r = Ured[a+b+d:]

	diffU1 = U1e - U1r
	diffU2 = U2e - U2r
	diffU4 = U4e - U4r
	diffU5 = U5e - U5r
	test1 = np.all(np.abs(diffU1) < tolerance)
	test2 = np.all(np.abs(diffU2) < tolerance)
	test4 = np.all(np.abs(diffU4) < tolerance)
	test5 = np.all(np.abs(diffU5) < tolerance)
	testall = test1 * test2 * test4 * test5
	testallall *= testall
	if not testall:
		print("For i = %d: %s" % (i, testall))
		if not test1:
			print("Max diff1= ", np.max(np.abs(diffU1)))
		if not test2:
			print("Max diff2= ", np.max(np.abs(diffU2)))
		if not test4:
			print("Max diff4= ", np.max(np.abs(diffU4)))
		if not test5:
			print("Max diff5= ", np.max(np.abs(diffU5)))

print("Final result = ", testallall)