from IPython.display import display, Image from IPython.display import display_pretty, display_html, display_jpeg, display_png, display_json, display_latex, display_svg from IPython.display import Latex from IPython.display import HTML Image(filename='logo3.jpg', width=100) uN = np.array(x*0) vN = np.array(avg(x)*0) uS = np.array(x*0) vS = np.array(avg(x)*0) uW = np.array(avg(y)*0+1.0) vW = np.array(y*0) uE = np.array(avg(y)*0+1.0) vE = np.array(y*0) Ue = op.cncat((uW, U, uE)) Ue = op.cncat((2 * op.vec2col(uS) - Ue[:, 0], Ue, 2 * op.vec2col(uN)-Ue[:, -1]), dim='h') Ve = op.cncat((vS.T, V, vN.T), dim='h') Ve = op.cncat((2 * op.vec2row(vW) - Ve[1, :], Ve, op.vec2row(vE) - Ve[:, -1])) Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond='Neumann')) + \ sparse.kron(Lap(ny, hy, boudary_cond='Neumann'), sparse.eye(nx)) from scipy.sparse.linalg import spsolve def generate_laplacians(ubx='Neumann', uby="Neumann", vbx="Neumann", vby="Neumann", pbx="Neumann", pby="Neumann", qbx="Dirichlet_between", qby="Dirichlet_between"): """Returns laplacians for vector field (U, V) and scalar fields P and stream function Q Arguments are just names of boundary conditions """ Lu = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny), Lap(nx-1, hx, boudary_cond=ubx)) + sparse.kron(Lap(ny, hy, boudary_cond=uby), sparse.eye(nx-1))) Lv = sparse.eye(nx*(ny-1))+dt / Re*(sparse.kron(sparse.eye(ny-1), Lap(nx, hx, boudary_cond=vbx)) + sparse.kron(Lap(ny-1, hy, boudary_cond=vby), sparse.eye(nx))) Lp = sparse.kron(sparse.eye(ny), Lap(nx, hx, boudary_cond=pbx)) + \ sparse.kron(Lap(ny, hy, boudary_cond=pby), sparse.eye(nx)) Lq = sparse.kron(sparse.eye(ny-1), Lap(nx-1, hx, boudary_cond=qbx)) + \ sparse.kron(Lap(ny-1, hy, boudary_cond=qby), sparse.eye(nx-1)) return Lu, Lv, Lp, Lq #velocity: rhs = np.ravel(U) U = spsolve(Lu, rhs, permc_spec='COLAMD') #colamd is one of preparation method used in spsolve U = np.reshape(U, (nx-1, ny)) rhs = np.ravel(V) V = spsolve(Lv, rhs, permc_spec='COLAMD') V = np.reshape(V, (nx, ny-1)) #pressure: rhs = np.ravel((np.diff(op.cncat((uW, U, uE)), axis=0)/hx + np.diff(op.cncat((vS.T, V, vN.T), dim='h'))/hy)) #gradient P = spsolve(Lp, -rhs, permc_spec='COLAMD') P = np.reshape(P[:], (nx, ny)) # 2.2a U = U - np.diff(P, axis=0)/hx # 2.2b V = V - np.diff(P)/hy rhs = np.diff(U)/hy - np.diff(V, axis=0)/hx rhs = np.ravel(rhs) q = spsolve(Lq, rhs, permc_spec='COLAMD') Q = np.zeros((nx+1, ny+1)) Q[1:-1, 1:-1] = np.reshape(q, (nx-1, ny-1)) rhs = np.diff(U)/hy - np.diff(V, axis=0)/hx rhs = np.ravel(rhs) q = spsolve(Lq, rhs, permc_spec='COLAMD') Q = np.zeros((nx+1, ny+1)) Q[1:-1, 1:-1] = np.reshape(q, (nx-1, ny-1)) [x_avg, y_avg] = np.meshgrid(avg(x), avg(y)) cxf = plt.contourf(x_avg, y_avg, P.T) cx = plt.contour(X, Y, Q.T) Ue = op.cncat((op.vec2col(uS), avg(op.cncat((uW, U, uE)).T).T, op.vec2col(uN)), dim='h') Ve = np.concatenate((op.vec2row(vW), avg(op.cncat((op.vec2col(vS), V, op.vec2col(vN)), dim='h')), op.vec2row(vE))) #Len = (Ue**2 + Ve**2 + eps)**0.5 #normalizing factor #ax = plt.quiver(x, y, (Ue/Len).T, (Ve/Len).T, scale=50, headwidth=4, headlength=3, headaxislength=3)#with normalisation #st = plt.streamplot(x, y, (Ue/Len).T, (Ve/Len).T) #stream function ax = plt.quiver(x, y, Ue.T, Ve.T, scale=50, headwidth=4, headlength=3, headaxislength=3) tx = fig.add_subplot(111) for obs in obstacles: #obstacles tx.add_patch(obs.get_patch()[0]) position=fig.add_axes([0.88,0.15,0.02,0.7]) cbar = plt.colorbar(cxf, cax = position) plt.draw() plt.clf() class Obstacle(): def __init__(self, lx, ly, sx, sy, boundary_cond='Dirichlet'): self.sx = sx self.sy = sy self.lx = lx self.ly = ly self.x0 = self.sx - self.lx/2.0 self.y0 = self.sy - self.ly/2.0 self.mesh = [] self.bc = boundary_cond self.polygon = plt.Rectangle((self.x0, self.y0), self.lx, self.ly, fill=True) self.patches = [self.polygon] def set_points(self, hx, hy, nx, ny): map = np.zeros((nx, ny)) map[round(self.x0/hx):round((self.x0+self.lx)/hx), round(self.y0/hy):round((self.y0+self.ly)/hy)] = 1 self.mesh = map def get_points(self): return self.mesh def insert_obj(self, U, V): #set velocity to zero within an object U[self.mesh[:-1, :]==1] = 0 V[self.mesh[:, :-1]==1] = 0 def get_patch(self): return self.patches obstacles = [] obstacles.append(Obstacle(0.1, 0.3, 0.2, 0.4)) for obs in obstacles: obs.set_points(hx, hy, nx, ny) #inside a loop: for obs in obstacles: obs.insert_obj(U, V) from IPython.display import YouTubeVideo YouTubeVideo("jV8QKn6Xcuc") from IPython.display import HTML HTML('')