diff --git a/fluidfoam/meshdesign.py b/fluidfoam/meshdesign.py index caab35d..648f02c 100644 --- a/fluidfoam/meshdesign.py +++ b/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