Adaptive solver gives error: Order "1" invalid for "Real" finite element

[cansu evcin]

I'm dealing with the solution of the forward facing step flow and I'm using mixed element and function spaces for the velocity, pressure and the lagrange multiplier(for the pressure). When I use the adaptive solver it fails with the error:

Order "1" invalid for "Real" finite element.

The elements and spaces that I generate are given below:

V_element = VectorElement('CG', triangle, 2)
Q_element = FiniteElement('CG', triangle, 1)
R_element = FiniteElement('R', triangle, 0)
W_element = MixedElement([V_element,Q_element,R_element])

W = FunctionSpace(mesh, W_element)
v, q, d = TestFunctions(W)
w       = Function(W)

u, p, c = split(w)

F = ...

J   = derivative(F,w)

pde = NonlinearVariationalProblem(F,w,bcs,J)
M   = inner(u,u)*dx
solver = AdaptiveNonlinearVariationalSolver(pde,M)
solver.solve(tol)

[Chris Richardson]

Can you ask first on the Q&A forum - fenicproject.org/qa

[cansu evcin]

I did it already before writing here, but it was asked and discussed before my question and it was not solved properly. Therefore, FEniCS Expert (MiroK) suggested me to reopen this issue in bitbucket. Thanks.

[Chris Richardson]

Sorry - I will try to format your message a bit

[cansu evcin]

Thanks for your interest.

[Jan Blechta]

This looks like a bug but we can't do anything without MWE. Could you, @cevcin, post a complete example to reproduce the problem?

[cansu evcin]

"""
steady_channel_adaptive.py 

solves the forward facing step channel flow problem with 
* inflow parabolic Dirichlet BC at inflow boundary
* homogenous Dirichlet BC at upper and bottom walls
  -nu* div(grad(u)) + (u* div)*u + grad(p) = 0, in (0,1)x(0,5)/(0.5,1)x(0,0.25)
                                    div(u) = 0, in (0,1)x(0,5)/(0.5,1)x(0,0.25)

*** standart boundary implementation for u in FEM
*** Lagrange multiplier approach for pressure

"""
from dolfin import *
from mshr import *
################################ SETUP ####################################
x01   = 0.0
x02   = 0.5
xend1 = 1.0
xend2 = 1.0
y01   = 0.0
y02   = 0.0
yend1 = 0.5
yend2 = 0.25
nx    = 20

Rectangle1      = Rectangle(Point(x01,y01), Point(xend1,yend1))
Rectangle2      = Rectangle(Point(x02,y02), Point(xend2,yend2))
channel         = Rectangle1 - Rectangle2
mesh            = generate_mesh(channel,nx)

V_element = VectorElement('CG', triangle, 2)
Q_element = FiniteElement('CG', triangle, 1)
R_element = FiniteElement('R', triangle, 0)
W_element = MixedElement([V_element,Q_element,R_element])

W = FunctionSpace(mesh, W_element)
V = FunctionSpace(mesh, V_element)
Q = FunctionSpace(mesh, Q_element)
R = FunctionSpace(mesh, R_element)
v, q, d = TestFunctions(W)
w       = Function(W)
u, p, c = split(w)

##########  Dirichlet boundary ###########
def inflow_boundary(x, on_boundary):
    return on_boundary and (near(x[0], 0.0))
def outflow_boundary(x, on_boundary):
    return on_boundary and (near(x[0], 1.0))
def top_boundary(x, on_boundary):
    return on_boundary and (near(x[1], 0.5))
def bottom1_boundary(x, on_boundary):
    return on_boundary and (near(x[1], 0.0))
def bottom2_boundary(x, on_boundary):
    return on_boundary and (near(x[1], 0.25)) and (x[0]>=(0.5+tol))  
def bottom3_boundary(x, on_boundary):
    return on_boundary and (near(x[0], 0.5)) and (x[1]<=(0.25+tol))

boundary_markers = FacetFunction("size_t", mesh) 
ds = Measure("ds", domain=mesh, subdomain_data=boundary_markers)   
n       = FacetNormal(mesh)
class Outflow_boundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0], 1.0))
outflow_boundary = Outflow_boundary()
outflow_boundary.mark(boundary_markers, 0)           
########## Dirichlet Functions ###########
Dirichlet1 = Expression(("1.0-((x[1]-0.25)/0.25)*((x[1]-0.25)/0.25)","0.0","0.0","0.0"),degree=2)
Dirichlet2 = Expression(("0.0","0.0","0.0","0.0"),degree=2)
dirichlet1 = project(Dirichlet1,W)
dirichlet2 = project(Dirichlet2,W)
v1,q1,r1 = dirichlet1.split()
v2,q2,r2 = dirichlet2.split()
########################### Dirichlet Boundary condition ##############################              
bc1 = DirichletBC(W.sub(0), v1, inflow_boundary)
bc2 = DirichletBC(W.sub(0), v2, top_boundary)
bc3 = DirichletBC(W.sub(0), v2, bottom1_boundary)
bc4 = DirichletBC(W.sub(0), v2, bottom2_boundary)
bc5 = DirichletBC(W.sub(0), v2, bottom3_boundary)

bcs = [bc1,bc2,bc3,bc4,bc5]
############################# Variational Form #############################
nu = Constant(1.0/500.0)
tol = 0.05
F = nu*inner(grad(u), grad(v))*dx\
    + inner(grad(u)*u, v)*dx \
    - p*div(v)*dx \
    + q*div(u)*dx \
    + c*q*dx\
    + p*d*dx\
    +inner(p*n,v)*ds(0)
############################### solver  ######################################
J   = derivative(F,w)
pde = NonlinearVariationalProblem(F,w,bcs,J)
M   = inner(w,w)*dx
solver = AdaptiveNonlinearVariationalSolver(pde,M)
solver.solve(tol)
solver.summary()
(u,p,c) = w.split()

[Jan Blechta]

Could you reduce it to absolutely necessary minimum to reproduce the problem? That's a point of MWE.

[cansu evcin]

from dolfin import *
mesh = UnitSquareMesh(20,20)
V_element = VectorElement('CG', triangle, 2)
Q_element = FiniteElement('CG', triangle, 1)
R_element = FiniteElement('R', triangle, 0)
W_element = MixedElement([V_element,Q_element,R_element])

W = FunctionSpace(mesh, W_element)
V = FunctionSpace(mesh, V_element)
Q = FunctionSpace(mesh, Q_element)
R = FunctionSpace(mesh, R_element)
v, q, d = TestFunctions(W)
w       = Function(W)
u, p, c = split(w)

bc = DirichletBC(W.sub(0), ("0.0","0.0"), "on_boundary")
############################# Variational Form #############################
F = inner(grad(u), grad(v))*dx\
    + inner(grad(u)*u, v)*dx \
    - p*div(v)*dx \
    + q*div(u)*dx \
    + c*q*dx\
    + p*d*dx

solve(F==0,w,bc)

[Jan Blechta]

Cannot reproduce the problem with any of latest dev version, 2017.1.0, 2016.x.0 using the last snippet.

[adequidt NA]

Hello.

Here is a MWE showing the problem with version 2017.1.0 (docker)

from dolfin import *

mesh = UnitSquareMesh(10, 10)
qelt = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
celt = FiniteElement("Real", mesh.ufl_cell(), 0)
X = FunctionSpace(mesh, qelt * celt)

w = Function(X)
u, l = split(w)
F = dot(grad(u), grad(u))/2*dx + l*(u-1)*dx
h = derivative(F, w, TestFunction(X))
bc = DirichletBC(X.sub(0), Constant(0.), "on_boundary")

M = u*dx()
tol = 1.e-6

solve(h == 0, w, bc, M=M, tol=tol)

[Jan Blechta]

Fix in https://bitbucket.org/fenics-project/ufl/pull-requests/80.

[Jan Blechta]

Fixed by https://bitbucket.org/fenics-project/ufl/commits/8ff2909ed888a5e4d8a4344850a3498f6d7d5be3.