# Mass matrix diagonalization (lumping)



In [1]:

from sympy import *
init_session()




IPython console for SymPy 1.0 (Python 3.6.0-64-bit) (ground types: python)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at http://docs.sympy.org/1.0/



## Elemental mass matrices



In [2]:

def mass_tet4():
"""Mass matrix for a 4 node tetrahedron"""
r, s, t = symbols("r s t")
N = Matrix([1 - r - s - t, r, s, t])
return (N * N.T).integrate((t, 0, 1 - r - s), (s, 0, 1 - r), (r, 0, 1))

"""Mass matrix for a 8 node quadrilateral"""
r, s = symbols("r s")
Haux = Matrix([
(1 - r**2)*(1 + s),
(1 - s**2)*(1 - r),
(1 - r**2)*(1 - s),
(1 - s**2)*(1 + r)])
N = S(1)/4*Matrix([
(1 + r)*(1 + s) - Haux[0] - Haux[3],
(1 - r)*(1 + s) - Haux[0] - Haux[1],
(1 - r)*(1 - s) - Haux[1] - Haux[2],
(1 + r)*(1 - s) - Haux[2] - Haux[3],
2*Haux[0], 2*Haux[1], 2*Haux[2], 2*Haux[3]])
return (N * N.T).integrate((s, -1, 1), (r, -1, 1))



The elemental mass matrices look like



In [3]:

mass_tet4()




Out[3]:

$$\left[\begin{matrix}\frac{1}{60} & \frac{1}{120} & \frac{1}{120} & \frac{1}{120}\\\frac{1}{120} & \frac{1}{60} & \frac{1}{120} & \frac{1}{120}\\\frac{1}{120} & \frac{1}{120} & \frac{1}{60} & \frac{1}{120}\\\frac{1}{120} & \frac{1}{120} & \frac{1}{120} & \frac{1}{60}\end{matrix}\right]$$




In [4]:




Out[4]:

$$\left[\begin{matrix}\frac{2}{15} & \frac{2}{45} & \frac{1}{15} & \frac{2}{45} & - \frac{2}{15} & - \frac{8}{45} & - \frac{8}{45} & - \frac{2}{15}\\\frac{2}{45} & \frac{2}{15} & \frac{2}{45} & \frac{1}{15} & - \frac{2}{15} & - \frac{2}{15} & - \frac{8}{45} & - \frac{8}{45}\\\frac{1}{15} & \frac{2}{45} & \frac{2}{15} & \frac{2}{45} & - \frac{8}{45} & - \frac{2}{15} & - \frac{2}{15} & - \frac{8}{45}\\\frac{2}{45} & \frac{1}{15} & \frac{2}{45} & \frac{2}{15} & - \frac{8}{45} & - \frac{8}{45} & - \frac{2}{15} & - \frac{2}{15}\\- \frac{2}{15} & - \frac{2}{15} & - \frac{8}{45} & - \frac{8}{45} & \frac{32}{45} & \frac{4}{9} & \frac{16}{45} & \frac{4}{9}\\- \frac{8}{45} & - \frac{2}{15} & - \frac{2}{15} & - \frac{8}{45} & \frac{4}{9} & \frac{32}{45} & \frac{4}{9} & \frac{16}{45}\\- \frac{8}{45} & - \frac{8}{45} & - \frac{2}{15} & - \frac{2}{15} & \frac{16}{45} & \frac{4}{9} & \frac{32}{45} & \frac{4}{9}\\- \frac{2}{15} & - \frac{8}{45} & - \frac{8}{45} & - \frac{2}{15} & \frac{4}{9} & \frac{16}{45} & \frac{4}{9} & \frac{32}{45}\end{matrix}\right]$$



## Lumping

One method for lumping is to sum the matrix per rows, i.e.

$$M^\text{(lumped)}_{ii}= \sum_{j} M_{ij}$$


In [5]:

def row_lump(mass_mat):
"""Matrix lumping by row summing"""
return diag(*[sum(mass_mat[i, :]) for i in range(mass_mat.shape[0])])



One method for lumping is to sum the matrix per rows, i.e.

$$M^\text{(lumped)}_{ii} = c M_{ii}$$

with $c$ adjusted to satisfy $\sum_j M^\text{(lumped)}_{jj} = \int_\Omega \rho d\Omega$. Particularly, we can choose $c = Tr(M)/M_\text{total}$.



In [6]:

def diag_scaling_lump(mass_mat):
"""Matrix lumping by diagonal scaling method"""
mass = sum(mass_mat)
trace = mass_mat.trace()
c = mass/trace
return diag(*[c*mass_mat[i, i] for i in range(mass_mat.shape[0])])




In [7]:

def min_dist_lump(mass_mat):
"""
Matrix lumping by minimizing the Frobenius norm subject
to a constraint of conservation of mass.
"""
num = mass_mat.shape[0]
mass = sum(mass_mat)
lamda = symbols("lambda")
Ms = symbols('M0:%d'%num)
var = list(Ms)
mass_diag = diag(*var)
C = mass_mat - mass_diag
fun = (C.T*C).trace() + lamda*(mass - sum(mass_diag))
var.append(lamda)
grad = [diff(fun, x) for x in var]
return diag(*list(sol.values())[:-1])



We can compare the methods for the tetrahedron



In [8]:

row_lump(mass_tet4())




Out[8]:

$$\left[\begin{matrix}\frac{1}{24} & 0 & 0 & 0\\0 & \frac{1}{24} & 0 & 0\\0 & 0 & \frac{1}{24} & 0\\0 & 0 & 0 & \frac{1}{24}\end{matrix}\right]$$




In [9]:

diag_scaling_lump(mass_tet4())




Out[9]:

$$\left[\begin{matrix}\frac{1}{24} & 0 & 0 & 0\\0 & \frac{1}{24} & 0 & 0\\0 & 0 & \frac{1}{24} & 0\\0 & 0 & 0 & \frac{1}{24}\end{matrix}\right]$$




In [10]:

min_dist_lump(mass_tet4())




Out[10]:

$$\left[\begin{matrix}\frac{1}{24} & 0 & 0 & 0\\0 & \frac{1}{24} & 0 & 0\\0 & 0 & \frac{1}{24} & 0\\0 & 0 & 0 & \frac{1}{24}\end{matrix}\right]$$



We can compare the methods for the serendipity quadrilaterals. For this type of element we can't use the row lumping method since it leads to negative masses.



In [11]:




Out[11]:

$$\left[\begin{matrix}- \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & - \frac{1}{3} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & - \frac{1}{3} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & - \frac{1}{3} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \frac{4}{3} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{4}{3} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & \frac{4}{3} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{4}{3}\end{matrix}\right]$$




In [12]:




Out[12]:

$$\left[\begin{matrix}\frac{3}{19} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{3}{19} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{3}{19} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & \frac{3}{19} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \frac{16}{19} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{16}{19} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & \frac{16}{19} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{16}{19}\end{matrix}\right]$$




In [13]:




Out[13]:

$$\left[\begin{matrix}\frac{19}{90} & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & \frac{19}{90} & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & \frac{19}{90} & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & \frac{19}{90} & 0 & 0 & 0 & 0\\0 & 0 & 0 & 0 & \frac{71}{90} & 0 & 0 & 0\\0 & 0 & 0 & 0 & 0 & \frac{71}{90} & 0 & 0\\0 & 0 & 0 & 0 & 0 & 0 & \frac{71}{90} & 0\\0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{71}{90}\end{matrix}\right]$$




In [ ]:




In [14]:

from IPython.core.display import HTML
def css_styling():
return HTML(styles)
css_styling()




Out[14]:

/* Based on Lorena Barba template available at: https://github.com/barbagroup/AeroPython/blob/master/styles/custom.css*/
@font-face {
font-family: "Computer Modern";
src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');
}
div.cell{
width:800px;
margin-left:16% !important;
margin-right:auto;
}
h1 {
font-family: 'Alegreya Sans', sans-serif;
}
h2 {
font-family: 'Fenix', serif;
}
h3{
font-family: 'Fenix', serif;
margin-top:12px;
margin-bottom: 3px;
}
h4{
font-family: 'Fenix', serif;
}
h5 {
font-family: 'Alegreya Sans', sans-serif;
}
div.text_cell_render{
font-family: 'Alegreya Sans',Computer Modern, "Helvetica Neue", Arial, Helvetica, Geneva, sans-serif;
line-height: 135%;
font-size: 120%;
width:600px;
margin-left:auto;
margin-right:auto;
}
.CodeMirror{
font-family: "Source Code Pro";
font-size: 90%;
}
/* .prompt{
display: None;
}*/
.text_cell_render h1 {
font-weight: 200;
font-size: 50pt;
line-height: 100%;
color:#CD2305;
margin-bottom: 0.5em;
margin-top: 0.5em;
display: block;
}
.text_cell_render h5 {
font-weight: 300;
font-size: 16pt;
color: #CD2305;
font-style: italic;
margin-bottom: .5em;
margin-top: 0.5em;
display: block;
}
.warning{
color: rgb( 240, 20, 20 )
}

MathJax.Hub.Config({
TeX: {
extensions: ["AMSmath.js"]
},
tex2jax: {
inlineMath: [ ['$','$'], ["\$","\$"] ],
displayMath: [ ['$$','$$'], ["\$","\$"] ]
},
displayAlign: 'center', // Change this to 'center' to center equations.
"HTML-CSS": {
styles: {'.MathJax_Display': {"margin": 4}}
}
});




In [ ]: