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$ May 29 '19 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$ May 29 '19 at 18:16
  • $\begingroup$ @WolfgangBangerth - Thanks! I couldn't find the FEniCS forum. Could you drop a link? $\endgroup$
    – user242538
    May 29 '19 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
    May 29 '19 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$ May 29 '19 at 21:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.