debug meshdesign.getgz function

fluidfoam/meshdesign.py

@@ -16,46 +16,52 @@ def getgz(h, dz1, N):
     in blockMesk and the z and dz vectors.
     Usage: z,dz,gz=getgz(h,dz1,N)
     """
-    # polynomial p of the sequence terms summation between 0 and N-1
-    p = np.zeros(N)
-    p[0] = h / dz1 - 1
-    p[1] = -h / dz1
-    p[N - 1] = 1
-    # Find the roots of the polynomial
-    sol = np.roots(p)
-    # Keep only the real non trivial solution
-    toto = np.where(np.imag(sol) == 0)
-    solreal = sol[toto]
-    titi = np.where(np.real(solreal) > 0)
-    titi = np.where(
-        np.logical_and(np.real(solreal) > 0, np.real(solreal) < 0.9999)
-    )
-    solgood = solreal[titi]
-    # Compute the common ratio of the sequence
-    try:
-        r = 1.0 / np.real(np.real(solgood[0]))
-    except IndexError:
-        print("dz1 so high compared to the number of cell and the domain size")
-        print("gz would be superior to 1 : NOT IMPLEMENTED")
-        r = 1.0
-        gz = 0
-    # Compute the grading factor for blockMesh
-    gz = r ** (N - 2)
-    # Compute the grid points position and the associated spacing
-    z = np.zeros(N)
-    dz = np.zeros(N - 1)
-    for i in range(N - 1):
-        z[i + 1] = z[i] + r ** i * dz1
-        dz[i] = z[i + 1] - z[i]
+    if dz1 == h/float(N):
+        print("Uniform mesh")
+        r = 1
+        gz = 1
+    else:
+        # polynomial p of the sequence terms summation between 0 and N-1
+        p = np.zeros(N+1)
+        p[N] =  h / dz1 - 1
+        p[N-1] = -h / dz1
+        p[0] = 1
+        # Find the roots of the polynomial
+        sol = np.roots(p)
+        # Keep only the real non trivial solution
+        toto = np.where(np.imag(sol) == 0)
+        solreal = sol[toto]
+        titi = np.where(np.real(solreal) > 0)
+        #print("Real solutions:",solreal)
+        titi = np.where(
+                np.logical_and(np.real(solreal) > 0, np.abs(np.real(solreal)-1) > 1e-6)
+                )
+        solgood = solreal[titi]
+        #print("Good solution:",solgood)
 
+        # Compute the common ratio of the sequence
+        try:
+            r = np.real(np.real(solgood[0]))
+        except IndexError:
+            print("There is a problem with input data or a bug")
+        # Compute the grading factor for blockMesh
+        gz = r ** (N - 1)
+    # Compute the grid points position and the associated spacing
+    z = np.zeros(N+1)
+    dz = np.zeros(N)
+    dz[0] = dz1
+    for i in range(N-1):
+        dz[i + 1] = r **(i+1) * dz1
+        z[i + 1] = z[i] + dz[i]
+    z[N] = z[N - 1] + dz[N - 1]
     # print some output
     print("grid sizes, common ratio and grading factor")
-    print("z[N-1]=", z[N - 1])
-    print("dz[0]=", dz[0])
-    print("dz[N-2]=", dz[N - 2])
-    print("common ratio: r=", r)
-    print("gz=", gz)
-    print("1/gz=", 1.0 / gz)
+    print("z[N]   =", round(z[N],4))
+    print("dz[0]  =", round(dz[0],4))
+    print("dz[N-1]=", round(dz[N - 1],4))
+    print("common ratio: r=", round(r,4))
+    print("gz=", round(gz,4))
+    print("1/gz=", round(1.0 / gz,4))
     return z, dz, gz