auto adaptive solver for mixed vector function spaces

Question

I've been playing around with the auto adaptive solver option, as explained in the demo Auto adaptive Poisson equation. I've been able to expand it to cases where there is a mixed function space, such as in the demo Poisson equation with pure Neumann boundary conditions, or when a vector function space is used. However, when I tried to use it for a problem that involved a mixed function space consisting of vector function spaces I get the error "Order "1" invalid for "Real" finite element." Here is the code:

# Create mesh
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
R = VectorFunctionSpace(mesh, "R", 0)
P = MixedFunctionSpace([V, R])

# define part of boundary
class BDY1(SubDomain):
def inside(self, x, on_boundary):
    return near(x[1], 0)

bdy1 = BDY1()

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
bdy1.mark(boundaries, 1)

ds = Measure("ds")[boundaries]

# Define functions
(u, c) = TrialFunctions(P)
(v, d) = TestFunctions(P)

t = Constant((1.0, 0.0))  # Traction force on the boundary

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

# Set up linear problem
a = (inner(grad(u),grad(v)) + inner(c, v) + inner(u,d))*dx

L = dot(t, v)*ds(1)

# Solve variational problem
p = Function(P)
(u, c) = split(p)

M = u[i]*u[i]*dx

solve(a == L, p, tol=0.001, M=M)

# Plot and hold solution
plot(p.root_node().sub(0), mode = "displacement")
plot(p.leaf_node().sub(0), mode = "displacement")

interactive()

If anyone knows why this doesn't work and how to fix it, I would be grateful if they let me know.

Answer 1

Hi, this is a bug that I'll report on bitbucket. Meanwhile, here's a workaround

V = VectorFunctionSpace(mesh, "Lagrange", 1)
R = FunctionSpace(mesh, "R", 0)
RR = MixedFunctionSpace([R, R])
P = MixedFunctionSpace([V, RR])

Increase the 'tol' to e.g 1e-6 to observe the adaptivity in action.

Comments

Thanks, that definitely does the trick. However, I don't understand why it works.

---

Issue reported here. The code there should clarify things for you.

---

I ran the code in the issue and I can't even run:

mesh = UnitTriangleMesh()

V = FunctionSpace(mesh, 'R', 0)
Vi = increase_order(V)

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

It can't increase the order for a Real element, has the implementation changed in these two years?

---

Hi, previously increase_order was a 'do nothing' operation on spaces of reals. Apparently, this has changed and the function is trying to construct real space with degree 1.