1- import numpy as np
21import itertools
2+
3+ import numpy as np
34import scipy .sparse
4- from .weights import weights
5- from . import operators as ops
6- from ..enums import bc , direction
5+
6+ from ..enums import bc
77from ..evaluation import evaluate
8+ from . import operators as ops
9+ from .weights import weights
810
911
1012def shift (grid , x , bcs ):
@@ -76,24 +78,25 @@ def interpolation(xs, grid, bcs):
7678 else :
7779 raise Exception ("interpolation not implemented for ndim > 3" )
7880 for it in tuples :
79- I = 0 # the index to push back
80- V = 1 # the value to push back
81+ II = 0 # the index to push back
82+ VV = 1 # the value to push back
8183 for kk in range (0 , grid .ndim ):
82- I = (I * grid .N [kk ] + nn [kk ]) * grid .n [kk ] + it [kk ]
83- V = V * px [kk ][it [kk ]]
84+ II = (II * grid .N [kk ] + nn [kk ]) * grid .n [kk ] + it [kk ]
85+ VV = VV * px [kk ][it [kk ]]
8486 # print( i, I, V)
8587 rows .append (i )
86- cols .append (round (I ))
88+ cols .append (round (II ))
8789 if not negative :
88- vals .append (V )
90+ vals .append (VV )
8991 else :
90- vals .append (- V )
92+ vals .append (- VV )
9193 # sort
92- rows , cols , vals = zip (* sorted (zip (rows , cols , vals )))
93- return scipy .sparse .coo_matrix ((vals , (rows , cols )), shape = (len (xs [0 ]),grid .size ()))
94+ rows , cols , vals = zip (* sorted (zip (rows , cols , vals , strict = False )), strict = False )
95+ return scipy .sparse .coo_matrix ((vals , (rows , cols )), shape = (len (xs [0 ]),
96+ grid .size ()))
9497
9598
96- def projection ( grid_new , grid_old ):
99+ def projection (grid_new , grid_old ):
97100 """ Create a projection between two grids
98101
99102 This matrix can be applied to vectors defined on the old (fine) grid to obtain
@@ -106,26 +109,28 @@ def projection( grid_new, grid_old):
106109 """
107110 ndim = grid_old .ndim
108111 if grid_new .ndim != ndim :
109- raise Exception ( "Cannot project between grids with different dimensions" )
110- for i in range ( 0 , ndim ):
111- if grid_old .N [i ] % grid_new .N [i ] != 0 :
112- print ( f"WARNING you project between incompatible grids!! old N: { grid_old .N [i ]} new N { grid_new .N [i ]} " )
113- if grid_old .n [i ] < grid_new .n [i ] :
114- print ( f"WARNING you project between incompatible grids!! old n: { grid_old .n [i ]} new n { grid_new .n [i ]} " )
115- wf = scipy .sparse .diags (weights ( grid_old ))
116- points = list ()
117- bcs = [bc .PER for i in range (0 ,ndim )]
112+ raise Exception ("Cannot project between grids with different dimensions" )
113+ for i in range (0 , ndim ):
114+ if grid_old .N [i ] % grid_new .N [i ] != 0 :
115+ print (f"WARNING you project between incompatible grids!! old N: \
116+ { grid_old .N [i ]} new N { grid_new .N [i ]} " )
117+ if grid_old .n [i ] < grid_new .n [i ]:
118+ print (f"WARNING you project between incompatible grids!! old n: \
119+ { grid_old .n [i ]} new n { grid_new .n [i ]} " )
120+ wf = scipy .sparse .diags (weights (grid_old ))
121+ points = []
122+ bcs = [bc .PER for i in range (0 , ndim )]
118123 if ndim == 1 :
119- points .append ( evaluate ( lambda x : x , grid_old ))
124+ points .append (evaluate (lambda x : x , grid_old ))
120125 elif ndim == 2 :
121- points .append ( evaluate ( lambda y ,x : y , grid_old ))
122- points .append ( evaluate ( lambda y ,x : x , grid_old ))
126+ points .append (evaluate (lambda y , x : y , grid_old ))
127+ points .append (evaluate (lambda y , x : x , grid_old ))
123128 elif ndim == 3 :
124- points .append ( evaluate ( lambda z ,y , x : z , grid_old ))
125- points .append ( evaluate ( lambda z ,y , x : y , grid_old ))
126- points .append ( evaluate ( lambda z ,y , x : x , grid_old ))
129+ points .append (evaluate (lambda z , y , x : z , grid_old ))
130+ points .append (evaluate (lambda z , y , x : y , grid_old ))
131+ points .append (evaluate (lambda z , y , x : x , grid_old ))
127132 else :
128133 raise Exception ("Projection not implemented for ndim > 3" )
129- A = interpolation ( points , grid_new , bcs )
130- vc = scipy .sparse .diags (1. / weights ( grid_new ))
134+ A = interpolation (points , grid_new , bcs )
135+ vc = scipy .sparse .diags (1. / weights (grid_new ))
131136 return vc @ A .transpose () @ wf
0 commit comments