这可能是因为网格的原因，但是对于所有 b 而不是 1 的值，以下代码都会崩溃。有人在 Fenics 中使用椭圆网格有经验吗？

import matplotlib
matplotlib.use('Agg')
from matplotlib import ticker
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *

def epsilon(u):
  return grad(u)+nabla_grad(u)

b=1
embryo = Ellipse(Point(0.0, 0.0), 1, b)
mesh = generate_mesh(embryo, 32)

# Define function spaces
P2 = VectorElement('CG', triangle, 5)
P1 = FiniteElement('CG', triangle, 3)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)
g = Constant(0.0)
mu = Constant(1.0)
force = Constant((0.0, 0.0))

# Specify Boundary Conditions
boundary = 'on_boundary'
flow_profile = ('-sin(atan2(x[1]*b,x[0]))*(a0+a1*sin(2*atan2(x[1]*b,x[0])))','cos(atan2(x[1]*b,x[0]))*(a0+a1*sin(2*atan2(x[1]*b,x[0])))')
bcu = DirichletBC(W.sub(0), Expression(flow_profile, degree=5, a0= -0.4202, a1 = 0.7653, b=b), boundary)
bc = [bcu]

# Define trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

a1 = inner( mu*epsilon(u)+p*Identity(2), nabla_grad(v))*dx + div(u)*q*dx
L1 = dot(force,v)*dx + g*q*dx
print('Preliminaries Done')

# Solve system
U = Function(W)
solve(a1==L1, U, bc)
print('Solving Done')

# Get sub-functions
u, p = U.split()

folder='./'
fig=plot(u, title='Velocity')
plt.colorbar(fig)
plt.savefig(folder+'vel.png', dpi=1000)
plt.close()
print('Plotting Done')