稳态stokes方程的标准算例

stokes.py

# 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)