DirichletBC with subdomains causes an error

12 months ago by
Hello all,
I am trying to impose Dirichlet B.C. on a 2D space and fenics throws an error when I try to ensemble the PETSc matrix. I am posting the minimal code and the error below. If I have only a single domain in the simulation, everything works well
from dolfin import *
from mshr import *
from pylab import show,triplot
from math import pi

G1 = Constant(1.0)
G2 = Constant(2.0)
w_sim = 10
h_sim = 5
w_wg = 4
h_wg = 3

class Waveguide(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (-w_wg/2, w_wg/2)) and
                    between(x[1], (-h_wg/2,h_wg/2 )))

# define simulation space with domains    
domain = Rectangle(Point(-w_sim/2,-h_sim/2), Point(w_sim/2,h_sim/2))
domain.set_subdomain(1, Rectangle(Point(-w_sim/2,-h_sim/2), Point(w_sim/2,h_sim/2)))
domain.set_subdomain(2, Rectangle(Point(-w_wg/2,-h_wg/2), Point(w_wg/2,h_wg/2)))
mesh = generate_mesh(domain,20)

# Initialize mesh function for interior domains
waveguide = Waveguide()
domains = CellFunction("size_t", mesh)
waveguide.mark(domains, 1)

# define function space
V = FunctionSpace(mesh, "Lagrange", 1)
# Define new measures associated with the interior domains 
dx = Measure("dx")(subdomain_data=domains)

# Define basis and bilinear form
u = TrialFunction(V)
v = TestFunction(V)
a = G1*inner(grad(u), grad(v))*dx(0) 
+ G2*inner(grad(u), grad(v))*dx(1) 

# Define Dirichlet boundary conditions

def boundary(x, on_boundary):
    return on_boundary
bcs = DirichletBC(V, 0.0,boundary)

# Assemble stiffness form
A = PETScMatrix()
assemble(a, tensor=A)

And the corresponding error message:
*** 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: 2016.2.0
Community: FEniCS Project

1 Answer

12 months ago by
Just change the definition of your bilinear form by

a = G1*inner(grad(u), grad(v))*dx(0) + G2*inner(grad(u), grad(v))*dx(1)


a = G1*inner(grad(u), grad(v))*dx(0) \
  + G2*inner(grad(u), grad(v))*dx(1)​
Thank you for spotting the error! I got stuck last night.
written 12 months ago by Marcin  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »