FFC warning: evaluate_dof(s) for enriched element not implemented

9 weeks ago by
Good morning:

I am trying to use bubble element in poroelastic problem as shown below:

However, there is a warning 'FFC warning: evaluate_dof(s) for enriched element not implemented'. Do you have any recommendations?

from dolfin import *
import matplotlib.pyplot as plt

mesh = UnitSquareMesh(64,64)

n = FacetNormal(mesh)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and near(x[1], 1.0))

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and near(x[0], 1.0))

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and near(x[1], 0.0))

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (on_boundary and near(x[0], 0.0))

top = Top()
right = Right()
bottom = Bottom()
left = Left()

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
top.mark(boundaries, 1)
right.mark(boundaries, 2)
bottom.mark(boundaries, 3)
left.mark(boundaries, 4)

P1 = FiniteElement("CG", triangle, 1)
B3 = FiniteElement("Bubble", triangle, 3)
U = VectorElement(P1 + B3)
PS = FiniteElement("CG", triangle, 1)
PF = FiniteElement("CG", triangle, 1)
W = FunctionSpace(mesh, MixedElement([U, PS, PF]))

def E_nu_to_mu_lmbda(E, nu):
    mu = E/(2*(1.0+nu))
    lmbda = (nu*E)/((1-2*nu)*(1+nu))
    return (mu, lmbda)

E = 14.e9
nu = 0.35
(mu_l, lmbda_l) = E_nu_to_mu_lmbda(E, nu)
mu = Constant(1.e-3)
alpha = Constant(0.6)
f = Constant((0.0, 0.0))
g = Constant(0.0)
k = Constant(1.0e-16)
f_stress_x = Constant(-24.e6)
f_stress_y = Constant(-32.e6)

(u, ps, pf) = TrialFunctions(W)
(y, q1, q2) = TestFunctions(W)

w = Function(W)
w0 = Function(W)
(u0, _, pf0) = split(w0)

pressure_dir_init = Constant(19.e6)

bc0_0 = DirichletBC(W.sub(0).sub(0), 0.0, boundaries, 1)
bc1_0 = DirichletBC(W.sub(0).sub(1), 0.0, boundaries, 3)
bcp1_0 = DirichletBC(W.sub(2), pressure_dir_init, boundaries, 1)
bcp2_0 = DirichletBC(W.sub(2), pressure_dir_init, boundaries, 2)
bcp3_0 = DirichletBC(W.sub(2), pressure_dir_init, boundaries, 3)
bcp4_0 = DirichletBC(W.sub(2), pressure_dir_init, boundaries, 4)
bcs_0 = [bc0_0, bc1_0, bcp1_0, bcp2_0, bcp3_0, bcp4_0]

I = Identity(mesh.topology().dim())
def sigma(u):
    return 2*mu_l*sym(grad(u)) + lmbda_l*div(u)*I

def strain(u):
    return sym(grad(u))

a0 =  inner(2*mu_l*strain(u)+ps*I-alpha*pf*I,sym( grad(y)))*dx\
    + inner(k/mu*grad(pf),grad(q2))*dx\
    + (1/lmbda_l)*inner(ps,q1)*dx\
    - inner(div(u),q1)*dx

L0 = inner(f,y)*dx\
    + g*q2*dx\
    + dot(f_stress_x*n,y)*ds(1) \
    + dot(f_stress_x*n,y)*ds(2) \
    + dot(f_stress_y*n,y)*ds(3) \
    + dot(f_stress_y*n,y)*ds(4)

solve(a0 == L0, w, bcs_0)
Community: FEniCS Project

1 Answer

9 weeks ago by
Hi, this is discussed here: https://bitbucket.org/fenics-project/ffc/issues/69/

You can use the following instead of your Dirichlet-BCs on W.sub(0).sub(i):

zero0 = project(Constant(0), W.sub(0).sub(0).collapse())
zero1 = project(Constant(0), W.sub(0).sub(1).collapse())
bc0_0 = DirichletBC(W.sub(0).sub(0), zero0, boundaries, 1)
bc1_0 = DirichletBC(W.sub(0).sub(1), zero1, boundaries, 3)​
It works!

May I ask if it is possible to explain the reasons behind this issue so I fix it myself next time?

Thank you so much!

written 8 weeks ago by Teeratorn Kadeethum  
Hi, curiously your original code seems to run?
I don't know how FEniCS handles this. In older versions you would have gotten an error instead of a warning.

Just make sure that when you apply a DirichletBC to an enriched (sub) space (like U = VectorElement(P1 + B3)) you project the value you want to impose to the correct space beforehand.
written 8 weeks ago by David Nolte  
Sorry for a late reply. Yes, it ran with lots of warnings.

Thank you!
written 8 weeks ago by Teeratorn Kadeethum  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »