This might be because of the mesh, but the following code blows up for all values of b not 1. Does anybody have any experience working with the ellipse mesh in Fenics?

import matplotlib
from matplotlib import ticker
import matplotlib.pyplot as plt
from dolfin import *
from mshr import *

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

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

fig=plot(u, title='Velocity')
plt.savefig(folder+'vel.png', dpi=1000)
print('Plotting Done')
  • $\begingroup$ You should be asking this question on the FEniCS forum. That's where all of the people who know FEniCS well hang out. $\endgroup$ Commented May 29, 2019 at 3:11
  • $\begingroup$ 1) Does this still happen when the velocity/pressure spaces are degree 2 and 1 respectively? I'm not sure that P5/P3 is stable for the Stokes equations, see pp. 5-6 of these notes. 2) Do you have the option of choosing a direct factorization solver? $\endgroup$ Commented May 29, 2019 at 18:16
  • $\begingroup$ @WolfgangBangerth - Thanks! I couldn't find the FEniCS forum. Could you drop a link? $\endgroup$
    – user242538
    Commented May 29, 2019 at 20:39
  • $\begingroup$ @DanielShapero - 1) The same issue emerges for P2/P1. Further, I am looking at 4th derivatives for the velocity field, so I definitely need at least P5. Further, due to non-linearities, I find that P5 is more accurate. I suspect that the lack of solution is actually because of how FEniCS handles boundary nodes because this issue only emerges for certain geometries. Thanks for the notes! I will have to read them! 2) Do you have a recommendation for a direct factorization solver? I find that most of those solvers do not converge. $\endgroup$
    – user242538
    Commented May 29, 2019 at 20:43
  • 1
    $\begingroup$ @DanielShapero: $P_k \times P_{k-1}$ elements are stable for Stokes for $k\ge 2$. If you make the velocity space larger, then the inf-sup constant can only improve, so $P_{k+1} \times P_{k-1}$ must also be stable. $\endgroup$ Commented May 29, 2019 at 21:46


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.