Mixed-dimensional branch - Error: Unable to set given (local) rows to identity matrix

Legacy DOLFIN
1

Hello! I am using the mixed-dimensional branch from @cdaversin in docker to solve a nonlinear problem on a 2D rectangle that is divided into 2 subdomains. I tried my best to create a MWE and I apologize in advance if it’s still too long:

from dolfin import *

# dimensions of rectangle
Length = 3.15
Height = 1.0

class TopDomain(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= Height/2

class BottomDomain(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= Height/2
    
class UpperBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], Height, eps=DOLFIN_EPS)
    
class LowerBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0.0, eps=DOLFIN_EPS)

# create mesh
mesh = RectangleMesh(Point(0.0, 0.0), Point(Length, Height), 10, 10, 'crossed')

# set domains
domains = MeshFunction('size_t', mesh, mesh.topology().dim())
domains.set_all(0)
top_domain = TopDomain()
top_domain.mark(domains, 1)
bottom_domain = BottomDomain()
bottom_domain.mark(domains, 2)

# create submeshes
top_mesh = MeshView.create(domains, 1)
bottom_mesh = MeshView.create(domains, 2)

# mark boundaries on submeshes
top_boundary_markers = MeshFunction('size_t', top_mesh, top_mesh.topology().dim()-1)
top_boundary_markers.set_all(0)
top_boundary = UpperBoundary()
top_boundary.mark(top_boundary_markers, 1)

bottom_boundary_markers = MeshFunction('size_t', bottom_mesh, bottom_mesh.topology().dim()-1)
bottom_boundary_markers.set_all(0)
bottom_boundary = LowerBoundary()
bottom_boundary.mark(bottom_boundary_markers, 2)

# function spaces
V1 = VectorFunctionSpace(top_mesh, 'Lagrange', 1)
V2 = VectorFunctionSpace(bottom_mesh, 'Lagrange', 1)
V = MixedFunctionSpace(V1, V2)

# Dirichlet BCs
b = Expression(("0.0", "0.0"), degree=1)
t = Expression(("0.0", "1.5"), degree=1)
bct = DirichletBC(V.sub_space(0), t, top_boundary)
bcb = DirichletBC(V.sub_space(1), b, bottom_boundary)
bcs = [bct, bcb]

u = Function(V)
u1, u2 = u.split()

dx1 = Measure('dx', domain=V.sub_space(0).mesh())
dx2 = Measure('dx', domain=V.sub_space(1).mesh())

# some simplified math
strain1 = 0.5 * (nabla_grad(u1) + nabla_grad(u1).T)
strain2 = 0.5 * (nabla_grad(u2) + nabla_grad(u2).T)
e_vgt1 = as_vector([strain1[0,0],strain1[1,1],2*strain1[0,1]])
e_vgt2 = as_vector([strain2[0,0],strain2[1,1],2*strain2[0,1]])
Pi = dot(e_vgt1,e_vgt1)*dx1 + dot(e_vgt2, e_vgt2)*dx2

Pi_blocks = extract_blocks(Pi)
F = []
for P in Pi_blocks:
    for ui in u.split():
        d = derivative(P, ui, TestFunctions(V))
        F.append(d)

J = []
for Fi in F:
    for ui in u.split():
        d = derivative(Fi, ui)
        J.append(d)

problem = MixedNonlinearVariationalProblem(F, u.split(), bcs, J)
solver = MixedNonlinearVariationalSolver(problem)    
solver.solve()

and the error I’m getting is:

*** 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:  68d9ed51bee6723163e7c4fbfffe31e6e3079bca

Edit: I have tried using the newer scientific computing from the docker ghcr.io/scientificcomputing/fenics:2024-05-13 and am still running into the same error

2

The mixed dimensional branch you are using is known to be buggy and to not work. Even the authors of the associated paper have stopped using it. They now rely on multifenics/multifenicsx for any multiphysics simulations. The unofficial recommendation is to follow their example :shushing_face:

3

Hi @Pelexi, welcome to the community.

As the author of multiphenics and multiphenicsx I appreciate you mentioning them, and I am very happy that you found them useful (I imagine, since you mentioned them).

At the same time, as an administrator of this forum and as a member of the FEniCS steering council, I would urge you to pay attention on how you phrase your posts. As it is currently stated, it may sound disrespectful towards the developers who worked on the mixed dimensional branch. You could have easily conveyed the same message with a simple “I have had similar issues in the past with the mixed dimensional branch. I can’t help you with solving them, since then I have switched to multiphenics/multiphenicsx: would you be willing to give it a try?”

2 Likes