Issue321 rectangular grid metis

tofu/__init__.py

@@ -91,7 +91,7 @@
     lsubout = ['tofu.{0}'.format(ss) for ss in lsubout]
     msg = "\nThe following subpackages are not available:"
     msg += "\n    - " + "\n    - ".join(lsubout)
-    msg += "\n  => see tofu.dsub[<subpackage>] for details."
+    msg += "\n  => see print(tofu.dsub[<subpackage>]) for details."
     warnings.warn(msg)
 
 # -------------------------------------

tofu/data/_comp.py

@@ -432,6 +432,8 @@ def func(pts, vect=None, t=None, ntall=ntall,
                             vquant[indtq[indtu[ii]], :],
                             trifinder=trifind
                         )(r, z).filled(fill_value)
+                if np.any(np.isnan(val)) and not np.isnan(fill_value):
+                    val[np.isnan(val)] = fill_value
                 return val, t
 
         # --------------------
@@ -463,6 +465,8 @@ def func(pts, vect=None, t=None, ntall=ntall,
                         ind = indt == indtu[ii]
                         val[ind, indok] = vquant[indtq[indtu[ii]],
                                                  indpts[indok]]
+                if np.any(np.isnan(val)) and not np.isnan(fill_value):
+                    val[np.isnan(val)] = fill_value
                 return val, t
 
 
@@ -518,14 +522,13 @@ def func(pts, vect=None, t=None, ntall=ntall,
                         # interpolate 1d
                         ind = indt == indtu[ii]
                         val[ind, ...] = scpinterp.interp1d(
-                            vr1[indtr1[indtu[ii]], :],
+                            vr1[indtr1[ii], :],
                             vquant[indtq[indtu[ii]], :],
                             kind='linear',
                             bounds_error=False,
                             fill_value=fill_value
                         )(np.asarray(vii))
                 val[np.isnan(val)] = fill_value
-
                 return val, t
 
         else:
@@ -563,12 +566,16 @@ def func(pts, vect=None, t=None, ntall=ntall,
                         # interpolate 1d
                         ind = indt == indtu[ii]
                         val[ind, indok] = scpinterp.interp1d(
-                            vr1[indtr1[indtu[ii]], :],
+                            vr1[indtr1[ii], :],
                             vquant[indtq[indtu[ii]], :],
                             kind='linear',
                             bounds_error=False,
                             fill_value=fill_value
                         )(vr2[indtr2[ii], indpts[indok]])
+
+                # Double check nan in case vr2 is itself nan on parts of mesh
+                if np.any(np.isnan(val)) and not np.isnan(fill_value):
+                    val[np.isnan(val)] = fill_value
                 return val, t
     return func
 
@@ -672,13 +679,13 @@ def func(pts, vect=None, t=None, ntall=ntall,
                  tall=tall, tbinall=tbinall):
 
             # Get pts in (r,z,phi)
-            r, z = np.hypot(pts[0,:],pts[1,:]), pts[2,:]
-            phi = np.arctan2(pts[1,:],pts[0,:])
+            r, z = np.hypot(pts[0, :], pts[1, :]), pts[2, :]
+            phi = np.arctan2(pts[1, :], pts[0, :])
 
             # Deduce vect in (r,z,phi)
-            vR = np.cos(phi)*vect[0,:] + np.sin(phi)*vect[1,:]
-            vphi = -np.sin(phi)*vect[0,:] + np.cos(phi)*vect[1,:]
-            vZ = vect[2,:]
+            vR = np.cos(phi)*vect[0, :] + np.sin(phi)*vect[1, :]
+            vPhi = -np.sin(phi)*vect[0, :] + np.cos(phi)*vect[1, :]
+            vZ = vect[2, :]
 
             # Prepare output
             shapeval = list(pts.shape)
@@ -689,6 +696,7 @@ def func(pts, vect=None, t=None, ntall=ntall,
 
             # Interpolate
             indpts = trifind(r, z)
+            indok = indpts > -1
             if t is None:
                 for ii in range(0,ntall):
                     valR[ii, ...] = vq2dR[indtq[ii], indpts]
@@ -697,9 +705,8 @@ def func(pts, vect=None, t=None, ntall=ntall,
                 t = tall
             else:
                 ntall, indt, indtu = plasma._get_indtu(t=t, tall=tall,
-                                                       tbinall=tbinall,
-                                                       idref1d=idref1d,
-                                                       idref2d=idref2d)[1:]
+                                                       tbinall=tbinall)[1:-2]
+
                 for ii in range(0, ntall):
                     ind = indt == indtu[ii]
                     valR[ind, ...] = vq2dR[indtq[indtu[ii]], indpts]

tofu/data/_core.py

@@ -2026,7 +2026,11 @@ def _extract_common_params(obj0, obj1=None):
 
     @staticmethod
     def _recreatefromoperator(d0, other, opfunc):
-        if type(other) in [int,float,np.int64,np.float64]:
+        if type(other) in [int, float, np.int64, np.float64]:
+            data = opfunc(d0.data, other)
+            dcom = d0._extract_common_params(d0)
+
+        elif isinstance(other, np.ndarray):
             data = opfunc(d0.data, other)
             dcom = d0._extract_common_params(d0)
 
@@ -2515,6 +2519,9 @@ def trifind(r, z,
                                 indpts = indR*nZ + indZ
                             else:
                                 indpts = indZ*nR + indR
+                            indout = ((r < R[0]) | (r > R[-1])
+                                      | (z < Z[0]) | (z > Z[-1]))
+                            indpts[indout] = -1
                             return indpts
 
                         dd[dk][k0]['R'] = R
@@ -3735,7 +3742,7 @@ def _get_finterp(self,
                      idquant=None, idref1d=None, idref2d=None,
                      idq2dR=None, idq2dPhi=None, idq2dZ=None,
                      interp_t=None, interp_space=None,
-                     fill_value=np.nan, ani=False, Type=None):
+                     fill_value=None, ani=False, Type=None):
 
         if interp_t is None:
             interp_t = 'nearest'
@@ -3853,7 +3860,7 @@ def _checkformat_qr12RPZ(self, quant=None, ref1d=None, ref2d=None,
     def get_finterp2d(self, quant=None, ref1d=None, ref2d=None,
                       q2dR=None, q2dPhi=None, q2dZ=None,
                       interp_t=None, interp_space=None,
-                      fill_value=np.nan, Type=None):
+                      fill_value=None, Type=None):
         """ Return the function interpolating (X,Y,Z) pts on a 1d/2d profile
 
         Can be used as input for tf.geom.CamLOS1D/2D.calc_signal()
@@ -3880,7 +3887,7 @@ def interp_pts2profile(self, pts=None, vect=None, t=None,
                            quant=None, ref1d=None, ref2d=None,
                            q2dR=None, q2dPhi=None, q2dZ=None,
                            interp_t=None, interp_space=None,
-                           fill_value=np.nan, Type=None):
+                           fill_value=None, Type=None):
         """ Return the value of the desired profiles_1d quantity
 
         For the desired inputs points (pts):
@@ -3911,8 +3918,23 @@ def interp_pts2profile(self, pts=None, vect=None, t=None,
                 else:
                     idmesh = [id_ for id_ in self._ddata[idref2d]['depend']
                               if self._dindref[id_]['group'] == 'mesh'][0]
-            pts = self.dmesh[idmesh]['data']['nodes']
-            pts = np.array([pts[:,0], np.zeros((pts.shape[0],)), pts[:,1]])
+            if self.dmesh[idmesh]['data']['type'] == 'rect':
+                if self.dmesh[idmesh]['data']['shapeRZ'] == ('R', 'Z'):
+                    R = np.repeat(self.dmesh[idmesh]['data']['R'],
+                                  self.dmesh[idmesh]['data']['nZ'])
+                    Z = np.tile(self.dmesh[idmesh]['data']['Z'],
+                                self.dmesh[idmesh]['data']['nR'])
+                else:
+                    R = np.tile(self.dmesh[idmesh]['data']['R'],
+                                self.dmesh[idmesh]['data']['nZ'])
+                    Z = np.repeat(self.dmesh[idmesh]['data']['Z'],
+                                  self.dmesh[idmesh]['data']['nR'])
+                pts = np.array(
+                    [R, np.zeros((self.dmesh[idmesh]['data']['size'],)), Z])
+            else:
+                pts = self.dmesh[idmesh]['data']['nodes']
+                pts = np.array(
+                    [pts[:, 0], np.zeros((pts.shape[0],)), pts[:, 1]])
 
         pts = np.atleast_2d(pts)
         if pts.shape[0] != 3:
@@ -3945,7 +3967,7 @@ def calc_signal_from_Cam(self, cam, t=None,
                              quant=None, ref1d=None, ref2d=None,
                              q2dR=None, q2dPhi=None, q2dZ=None,
                              Brightness=True, interp_t=None,
-                             interp_space=None, fill_value=np.nan,
+                             interp_space=None, fill_value=None,
                              res=0.005, DL=None, resMode='abs', method='sum',
                              ind=None, out=object, plot=True, dataname=None,
                              fs=None, dmargin=None, wintit=None, invert=True,

tofu/data/_plot.py

@@ -364,9 +364,9 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax,
     c0 = [all([dd.dlabels[kk] == lData[0].dlabels[kk] for dd in lData[1:]])
           for kk in ['t','X','data']]
     if not all(c0):
-        msg = "All Data objects must have the same:\n"
+        msg = "All Data objects do not have the same:\n"
         msg += "    dlabels['t'], dlabels['X'] and dlabels['data'] !"
-        raise Exception(msg)
+        warnings.warn(msg)
 
 
     # ---------
@@ -387,8 +387,9 @@ def _DataCam12D_plot(lData, key=None, nchMax=_nchMax, ntMax=_ntMax,
     # Check nch and X
     c0 = [dd.nch == lData[0].nch for dd in lData[1:]]
     if not all(c0):
-        msg = "All Data objects must have the same number of channels (self.nch)"
-        msg += "\nYou can set the indices of the channels with self.set_indch()"
+        msg = ("All Data objects must have the same nb. of channels\n"
+               + "\t- self.nch: {}\n".format([dd.nch for dd in lData])
+               + "\n  => use self.set_indch()")
         raise Exception(msg)
     nch = lData[0].nch
 

tofu/geom/_core.py

@@ -5968,6 +5968,9 @@ def _calc_signal_postformat(
         if Brightness is False:
             if dataname is None:
                 dataname = r"LOS-integral x Etendue"
+            if E is None or np.all(np.isnan(E)):
+                msg = "Cannot use etendue, it was not set properly !"
+                raise Exception(msg)
             if t is None or len(t) == 1 or E.size == 1:
                 sig = sig * E
             else:
@@ -6207,7 +6210,7 @@ def calc_signal_from_Plasma2D(
         self,
         plasma2d,
         t=None,
-        newcalc=False,
+        newcalc=True,
         quant=None,
         ref1d=None,
         ref2d=None,
@@ -6252,7 +6255,9 @@ def calc_signal_from_Plasma2D(
 
         if newcalc:
             # Get time vector
-            if t is None:
+            lc = [t is None, type(t) is str, type(t) is np.ndarray]
+            assert any(lc)
+            if lc[0]:
                 out = plasma2d._checkformat_qr12RPZ(
                      quant=quant,
                      ref1d=ref1d,
@@ -6262,6 +6267,8 @@ def calc_signal_from_Plasma2D(
                      q2dZ=q2dZ,
                 )
                 t = plasma2d._get_tcom(*out[:4])[0]
+            elif lc[1]:
+                t = plasma2d._ddata[t]['data']
             else:
                 t = np.atleast_1d(t).ravel()
 
@@ -6359,7 +6366,7 @@ def funcbis(*args, **kwdargs):
                 nbrep = np.r_[
                     indpts[0], np.diff(indpts), pts.shape[1] - indpts[-1]
                 ]
-                vect = np.repeat(self.u, nbrep, axis=1)
+                vect = -np.repeat(self.u, nbrep, axis=1)
             if fill_value is None:
                 fill_value = 0.