### DirichletBC with subdomains causes an error

227

views

0

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

And the corresponding error message:

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)
domains.set_all(0)
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)
bcs.apply(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

1

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)`

or

```
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
9 months ago by
Marcin

Please login to add an answer/comment or follow this question.