From 4a1868878f92ad97700bb82151bd309abdf550e1 Mon Sep 17 00:00:00 2001
From: weisse02 <andreas.weissenbrunner@ptb.de>
Date: Wed, 22 Feb 2023 17:42:43 +0100
Subject: [PATCH] Changed the regGridinterpolators get_profile works

---
 Flow_class.py | 109 +++++++++++++++-----------------------------------
 1 file changed, 33 insertions(+), 76 deletions(-)

diff --git a/Flow_class.py b/Flow_class.py
index 091ddbd..cbe9045 100644
--- a/Flow_class.py
+++ b/Flow_class.py
@@ -117,10 +117,9 @@ def loadPODres(case):
     #
     case_ind = [0, 2121, 26321, 44592]
     # distinguish between the cases 0, 1, 2 for SE,DE,DSE 
-    data["coeffs"] = data["coeffs"][case_ind[ci-1]:case_ind[ci]]
-    
-    
-    
+    coeffs = data["coeffs"][case_ind[ci-1]:case_ind[ci]]
+    modes  = data["modes"]
+    umean  = data["umean"]
     # if case == "SingleElbow":
     #     # the POD res files are in npz-format
     #     ux = np.load(data_path+"/Ux_POD_SE.npz")
@@ -137,13 +136,12 @@ def loadPODres(case):
     #     uy = np.load(data_path+"/Uy_POD_DSE.npz")
     #     uz = np.load(data_path+"/Udiff_POD_DSE.npz")
     #
-    return data
-
-
-            
-
-
-
+    return coeffs,modes, umean
+#
+#
+#  ------------------------------  Start Class Definition  -----------------
+#
+#
 class Elbow_profile():
     
     def loadCoordinates(self,case,delete_last=False):
@@ -201,15 +199,10 @@ class Elbow_profile():
         # Nrk   = self.Rk.size
         # Nt    = self.dist.size * self.Rk.size
         #
-        pod = loadPODres(case)
-        self.Nmodes    = pod["modes"].shape[-1]
-        # die moden aufteilen oder nicht??
-        self.modes =  pod["modes"]
-        # the mean profiles
-        self.ux_mean  =  pod_x["umean"]
-        self.uy_mean  =  pod_y["umean"]
-        self.uz_mean  =  pod_z["umean"]
-        
+        coeffs, self.modes, self.umean  = loadPODres(case)
+        self.Nmodes    = self.modes.shape[-1]
+        # split the modes in x, y, z ? better not because of the reconstruction?
+        #
         if case=="SingleElbow":
             # save the modes in local variables and add zeros at the wall coordinates
             # self.modes_ux =  np.concatenate((pod_x["modes"].reshape(self.Nphi,self.Nr-1,self.Nmodes),np.zeros((self.Nphi,1,self.Nmodes))),axis=1).reshape(-1,self.Nmodes)
@@ -220,10 +213,7 @@ class Elbow_profile():
             # self.uy_mean  =  np.concatenate((pod_y["umean"].reshape(self.Nphi,self.Nr-1),np.zeros((self.Nphi,1))),axis=1).ravel()
             # self.uz_mean  =  np.concatenate((pod_z["umean"].reshape(self.Nphi,self.Nr-1),np.zeros((self.Nphi,1))),axis=1).ravel()
             #
-            #
-            coeffs = pod_x["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
-            coeffs_y = pod_y["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
-            coeffs_z = pod_z["coeffs"].reshape((len(self.Rk),len(self.dist),self.Nmodes))
+            coeffs = coeffs.reshape((len(self.Rk),len(self.dist),self.Nmodes))
         else: 
             # self.modes_ux =  pod_x["modes"]
             # self.modes_uy =  pod_y["modes"]
@@ -233,9 +223,7 @@ class Elbow_profile():
             # self.uy_mean  =  pod_y["umean"]
             # self.uz_mean  =  pod_z["umean"]
             #
-            coeffs_x = pod_x["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
-            coeffs_y = pod_y["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
-            coeffs_z = pod_z["coeffs"].reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
+            coeffs = coeffs.reshape((len(self.Rk),len(self.dl),len(self.dist),self.Nmodes))
         #
         # plot for test 
         # indrk    = 7
@@ -249,43 +237,23 @@ class Elbow_profile():
         # nicecontour(self.x,self.y,utest)
         
         #polydeg  = 3 
-        self.regint_ux   = []
-        self.regint_uy   = []
-        self.regint_uz   = []
+        # self.regint_ux   = []
+        # self.regint_uy   = []
+        # self.regint_uz   = []
         methodstr = 'linear' # 'linear' 'cubic' or 'quintic' are possible
         # the spline interpolation attention it can only evaluate ordered input arrays 
         # RegularGridInterpolator can also handle multiple output values --> all modes a loop is not neccessary!!!!!
-        # to do: Remove the loops, regint is not a list but a interpolater object 
+        # but only for linear interpolation not for higher order polynomials! 
         if case=="SingleElbow":
-            for i in range(0,self.Nmodes):
-                # self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                # self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                # self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                self.regint_ux.append(interpol.RegularGridInterpolator((self.Rk,self.dist), coeffs_x[:,:,i],method = methodstr, bounds_error = False,fill_value=None))#,bounds_error = False,fill_value=None ))
-                self.regint_uy.append(interpol.RegularGridInterpolator((self.Rk,self.dist), coeffs_y[:,:,i],method = methodstr, bounds_error = False,fill_value=None))#,bounds_error = True,fill_value=None ))
-                self.regint_uz.append(interpol.RegularGridInterpolator((self.Rk,self.dist), coeffs_z[:,:,i],method = methodstr, bounds_error = False,fill_value=None))#,bounds_error = True,fill_value=None ))
+            self.regint_A = interpol.RegularGridInterpolator((self.Rk,self.dist), coeffs,method = methodstr, bounds_error = False,fill_value=None)
         else: 
-            for i in range(0,self.Nmodes):
-                # self.splineint_ux.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_x["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                # self.splineint_uy.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_y["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                # self.splineint_uz.append(interpol.RectBivariateSpline(self.dist,self.Rk,self.dl,pod_z["coeffs"][:,:,i].T,kx= polydeg,ky=polydeg ))
-                self.regint_ux.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_x[:,:,:,i],method = methodstr ,bounds_error = False,fill_value=None ))
-                self.regint_uy.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_y[:,:,:,i],method = methodstr ,bounds_error = False,fill_value=None ))
-                self.regint_uz.append(interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs_z[:,:,:,i],method = methodstr ,bounds_error = False,fill_value=None ))
-            #self.coeffs_z = coeffs_z
-        del coeffs_x
-        del coeffs_y
-        del coeffs_z 
+            self.regint_A = interpol.RegularGridInterpolator((self.Rk,self.dl,self.dist), coeffs,method = methodstr ,bounds_error = False,fill_value=None )
         #
+        del coeffs
         #
-        #self.pod_z = pod_z
-        del pod_x 
-        del pod_y 
-        del pod_z
-        
         #
         self.int_weights = int2d.get_int_weights(x=self.x,y=self.y,method="tri")
-        
+        #
         self.make_fully()
         # The parametrized function to be plotted
         #self.get_profile = interpol.RegularGridInterpolator((self.Rk,self.dist), self.uall + self.u_fully)
@@ -297,38 +265,27 @@ class Elbow_profile():
         # to evaluate the spline interpolation and get the coefficient matrix at the desired 
         # values of dist and Rk
         if self.case=="SingleElbow":
-            A_ux = np.array([int_i((Rk,dist)).T for int_i in self.regint_ux] ).T
-            A_uy = np.array([int_i((Rk,dist)).T for int_i in self.regint_uy] ).T
-            A_uz = np.array([int_i((Rk,dist)).T for int_i in self.regint_uz] ).T
+            A    = self.regint_A((Rk,dist))
+            # A_uz = np.array([int_i((Rk,dist)).T for int_i in self.regint_uz] ).T
         else:
-            A_ux = np.array([int_i((Rk,dl,dist)).T for int_i in self.regint_ux] ).T
-            A_uy = np.array([int_i((Rk,dl,dist)).T for int_i in self.regint_uy] ).T
-            A_uz = np.array([int_i((Rk,dl,dist)).T for int_i in self.regint_uz] ).T
+            A    = self.regint_A((Rk,dl,dist))
+            # A_ux = np.array([int_i((Rk,dl,dist)).T for int_i in self.regint_ux] ).T
         #
         # print("this is A for dist = " + str(dist) + " and Rk = " + str(Rk) + ":")
         # print(A_uz)
         #
-        A_ux = A_ux.reshape(-1,A_ux.shape[-1])
-        A_uy = A_uy.reshape(-1,A_uy.shape[-1])
-        A_uz = A_uz.reshape(-1,A_uz.shape[-1])
-        # now reconstruct the velocities
-        # print(A_ux.shape)
+        A    = A.reshape(-1,A.shape[-1])
         # print(A_uy.shape)
         # print(A_uz.shape)
         if asmat:
-            ux = do_reconstruct(A_ux, self.modes_ux, umean = self.ux_mean).reshape(self.xm.shape)
-            uy = do_reconstruct(A_uy, self.modes_uy, umean = self.uy_mean).reshape(self.xm.shape)
-            uz = do_reconstruct(A_uz, self.modes_uz, umean = self.uz_mean).reshape(self.xm.shape)
-            # if the dimensions of a is larger:
-            # b1 = bla.reshape(*bla.shape[0:-1],self.xm.shape)
+            # reshape to u.shape = 3 * x.shape
+            u  = do_reconstruct(A,self.modes,umean = self.umean ).reshape(3,*self.xm.shape)
         else:
-            ux = do_reconstruct(A_ux, self.modes_ux, umean = self.ux_mean)
-            uy = do_reconstruct(A_uy, self.modes_uy, umean = self.uy_mean)
-            uz = do_reconstruct(A_uz, self.modes_uz, umean = self.uz_mean)
+            u  = do_reconstruct(A, self.modes, umean = self.umean)
         if addfully:
-            uz += self.u_fully
+            u[2] += self.u_fully
         #
-        return ux,uy,uz
+        return u
     
     def get_path_profile_diam_se(self,Rk,dist,reverse =False, addfully = False,R = 1):
         # returns the velocity profile at a certain distance and coordinates x,y, (if given)
-- 
GitLab