[cig-commits] [commit] master: py: Simplify Jacobian calculation. (98097d3)

[cig_noreply at geodynamics.org]

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

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

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

commit 98097d3f8bc7514c5f10c20d119c60c7faec1954
Author: Elliott Sales de Andrade <esalesde at physics.utoronto.ca>
Date:   Sun Jan 11 04:41:40 2015 -0500

    py: Simplify Jacobian calculation.
    
    Using some NumPy functions, we can create the 2D array rather quickly
    without double loops.


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

98097d3f8bc7514c5f10c20d119c60c7faec1954
 Python_version/gll.py | 16 ++++++----------
 1 file changed, 6 insertions(+), 10 deletions(-)

diff --git a/Python_version/gll.py b/Python_version/gll.py
index fe81357..20fe8cc 100644
--- a/Python_version/gll.py
+++ b/Python_version/gll.py
@@ -122,11 +122,9 @@ def jacobian(ticks, param):
        (and GLJ points for the first element in axisymmetric).
        Returns a matrix nSpec*(N+1) containing its value for each element and
        each points"""
-    dE=ticks[1:]-ticks[:len(ticks)-1] # dE[i] : length of element number i
-    dXdKsi=np.zeros((len(dE),len(param.ksiGLL)),dtype='d')
-    for e in np.arange(len(dE)):
-        for k in np.arange(len(param.ksiGLL)):
-            dXdKsi[e,k] = dE[e]/2 # Here it does not depend on ksi
+    Np1 = len(param.ksiGLL)
+    dE = np.diff(ticks) / 2
+    dXdKsi = np.repeat(dE, Np1).reshape((-1, Np1))
     return dXdKsi
 
 
@@ -135,9 +133,7 @@ def jacobian_inverse(ticks, param):
        (and GLJ points for the first element in axisymmetric).
        Returns a matrix nSpec*(N+1) containing its value for each element and
        each points"""
-    dE=ticks[1:]-ticks[:len(ticks)-1] # dE[i] : length of element number i
-    dKsiDx=np.zeros((len(dE),len(param.ksiGLL)),dtype='d')
-    for e in np.arange(len(dE)):
-        for k in np.arange(len(param.ksiGLL)):
-            dKsiDx[e,k] = 2/dE[e]   # Here it does not depend on ksi
+    Np1 = len(param.ksiGLL)
+    dE = 2 / np.diff(ticks)
+    dKsiDx = np.repeat(dE, Np1).reshape((-1, Np1))
     return dKsiDx