Skip to content

Instantly share code, notes, and snippets.

@shaoyaoqian
Last active November 14, 2023 08:09
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save shaoyaoqian/664fb1ba512f802baa8a8643d8325a79 to your computer and use it in GitHub Desktop.
稳态stokes方程的标准算例
# import matplotlib.pyplot as plt
from dolfin import *
mesh = UnitSquareMesh(160, 160)
# Define function spaces
P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
TH = P2 * P1
W = FunctionSpace(mesh, TH)
# No-slip boundary condition for velocity
# Boundaries
def right(x, on_boundary): return x[0] > (1.0 - DOLFIN_EPS)
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def bottom(x, on_boundary): return x[1] < DOLFIN_EPS
def top(x, on_boundary): return x[1] > (1.0 - DOLFIN_EPS)
# No-slip boundary condition for velocity
noslip = Constant((0.0, 0.0))
bc3 = DirichletBC(W.sub(0), noslip, right)
bc1 = DirichletBC(W.sub(0), noslip, left)
bc2 = DirichletBC(W.sub(0), noslip, bottom)
# Inflow boundary condition for velocity
topflow = Expression(("1", "0.0"), degree=1)
bc0 = DirichletBC(W.sub(0), topflow, top)
# Collect boundary conditions
bcs = [bc0, bc1, bc2, bc3]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
a = (inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx
L = inner(f, v)*dx
# Compute solution
w = Function(W)
solve(a == L, w, bcs)
# Split the mixed solution using deepcopy
# (needed for further computation on coefficient vector)
(u, p) = w.split(True)
File("velocity.pvd") << u
File("pressure.pvd") << p
Nx = 100
Ny = 100
i = 5
# for i in range(Nx):
for j in range(Ny):
# print(i/Nx+0.5/Nx, j/Ny)
result = "{:.16e}".format(u(i/Nx+0.5/Nx, j/Ny)[1])
print(result)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment