How to fix the "Unable to set given (local) rows to identity matrix" error

1

Hi everyone,

I am trying to simulate steady flow on an upper half of a square domain. I am using the MeshView to create submesh and I use FEniCS on docker mixed-dimensional branch on Docker.
I got the following error,

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to set given (local) rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  f0a90b8fffc958ed612484cd3465e59994b695f9
*** -------------------------------------------------------------------------

My MWE looks like the below snippet.

from dolfin import *

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)
    
class Left1(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.0) and between(x[1], (0.5, 1))

class Right1(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 1.0) and between(x[1], (0.5, 1))
 
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5)
        
class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)
    
# Create meshes
n = 20
mesh = UnitSquareMesh(n, n)

# Domain markers
domain = MeshFunction("size_t", mesh, mesh.geometry().dim(), 0)
for c in cells(mesh):
  domain[c] = c.midpoint().y() > 0.5

# Interface markers
facets = MeshFunction("size_t", mesh, mesh.geometry().dim()-1, 0)
Interface().mark(facets,1)
  
# Meshes
mesh_U = MeshView.create(domain, 1)
mesh_L = MeshView.create(domain, 0)
mesh_I = MeshView.create(facets, 1)

# Function spaces
VE = VectorElement("CG", mesh.ufl_cell(), 1)   # velocity
PE = FiniteElement("CG", mesh.ufl_cell(), 1)   # Pressure
U_U  = FunctionSpace(mesh_U,VE)    # Velocity on the upper domain
P_U  = FunctionSpace(mesh_U,PE)    # Pressure on the upper domain
W = MixedFunctionSpace(U_U, P_U)#, T_U, T_L, L_I_temp)

# Markers for Dirichlet and Neuman bcs
ff_U = MeshFunction("size_t", mesh_U, mesh_U.geometry().dim()-1, 0)
Top().mark(ff_U, 1)
Left1().mark(ff_U, 2)
Right1().mark(ff_U, 3)
Interface().mark(ff_U, 4)

# Boundary conditions
bc_U1 = DirichletBC(W.sub_space(0), Constant((0,0)), ff_U, 1)
bc_U2 = DirichletBC(W.sub_space(0), Constant((0,0)), ff_U, 4)
bc_U3 = DirichletBC(W.sub_space(1), Constant(80), ff_U, 2)
bc_U4 = DirichletBC(W.sub_space(1), Constant(0), ff_U, 3)
bcs = [bc_U1, bc_U2, bc_U3, bc_U4]

# Measures
dx_U  = Measure("dx", domain=W.sub_space(0).mesh())
ds_U  = Measure("ds", domain=W.sub_space(0).mesh(), subdomain_data=ff_U)

# Function and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)

# Parameters for the variational problem
mu  = 0.89

# create function space
w = Function(W)
u = w.sub(0)
p = w.sub(1)

# Variational problem
F = rho*inner(grad(u)*u, v)*dx_U \
    + mu*inner(grad(u), grad(v))*dx_U \
    + inner(grad(p),v)*dx_U \
    - div(u)*q*dx_U \

# solving the problem     
solve(F == 0, w, bcs)

Can anyone point me out what is the cause of this error?

Thank you