[cig-commits] [commit] master: py: Simplify mass matrix creation slightly. (420d361)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/specfem1d

On branch  : master
Link       : https://github.com/geodynamics/specfem1d/compare/2ccfba13481e26296eebf034518cd31d2c727906...33b512de01491f2c70de206939855d95b22febf4

>---------------------------------------------------------------

commit 420d3619ff4d8c4884330ddcfb6f902a1610f837
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date:   Sun Jan 11 03:30:25 2015 -0500

    py: Simplify mass matrix creation slightly.
    
    Drop the extra cases which are redundant at the moment. If the commented
    stuff needs to be re-enabled, it can be added later. We can drop the
    inner loop as well just to save a tiny bit.


>---------------------------------------------------------------

420d3619ff4d8c4884330ddcfb6f902a1610f837
 Python_version/functions.py | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/Python_version/functions.py b/Python_version/functions.py
index 049ec09..f49b7de 100644
--- a/Python_version/functions.py
+++ b/Python_version/functions.py
@@ -55,16 +55,18 @@ def make_stiffness_matrix(grid, param):
     return Ke
 
 
-def make_mass_matrix(grid,param):
+def make_mass_matrix(grid, param):
     """Computation of global mass matrix"""
-    M=np.zeros(param.nGlob,dtype='d')
-    for e in np.arange(param.nSpec):
-        for i in np.arange(param.nGLL):
-            if param.axisym and e != 0:
-                Me=param.wGLL[i]*grid.rho[e,i]*grid.dXdKsi[e,i]#*project_inverse(param.ksiGLL[i],e,grid.ticks) # Axisym
-            if param.axisym and e == 0:
-                Me=param.wGLJ[i]*grid.rho[e,i]*grid.dXdKsi[e,i]
-            else:
-                Me=param.wGLL[i]*grid.rho[e,i]*grid.dXdKsi[e,i]#*project_inverse(param.ksiGLL[i],e,grid.ticks)
-            M[param.ibool[e,i]]+=Me
+    M = np.zeros(param.nGlob)
+    # NOTE: We cannot simply drop this outer loop because param.ibool contains
+    # repeated indices into M and the addition would not work correctly if we
+    # tried to slice it in.
+    for e in range(param.nSpec):
+        if param.axisym and e == 0:
+            Me = param.wGLJ * grid.rho[e,:] * grid.dXdKsi[e,:]
+        else:
+            Me = param.wGLL * grid.rho[e,:] * grid.dXdKsi[e,:]
+        # We can still slice along the second index because the individual GLL
+        # points are assumed to not repeat a point within an spectral element.
+        M[param.ibool[e,:]] += Me
     return M