I wrote some examples that show how extending the domain causes the error (even if they do not describe my complete problem). Please note that the script can be run directly, the only requirement is to have the `gmsh`

python API installed, nothing else (therefore, even if they might not be the shortest working examples possible, it is really easy to run them).

In the script below, I’m using only the domain 1, and everything works fine:

```
# DOMAIN 1 ---> WORKING FINE
from dolfin import *
import numpy as np
import gmsh
import sys
import meshio
um = 10**-6
def msh_to_mf():
dim = 2
cell_str = "triangle"
facet_str = "line"
msh = meshio.read("mesh.msh")
cells = np.vstack(np.array([cells.data for cells in msh.cells
if cells.type == cell_str]))
total_mesh = meshio.Mesh(points=msh.points,
cells=[(cell_str, cells)])
total_mesh.prune_z_0()
cells_data = msh.cell_data_dict["gmsh:physical"][cell_str]
cells_mesh = meshio.Mesh(points=msh.points,
cells=[(cell_str, cells)],
cell_data={"name_to_read": [cells_data]})
cells_mesh.prune_z_0()
facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
if cells.type == facet_str]))
facet_data = msh.cell_data_dict["gmsh:physical"][facet_str]
facet_mesh = meshio.Mesh(points=msh.points,
cells=[(facet_str, facet_cells)],
cell_data={"name_to_read": [facet_data]})
facet_mesh.prune_z_0()
meshio.write("total.xdmf", total_mesh)
meshio.write("facets.xdmf", facet_mesh)
meshio.write("cells.xdmf", cells_mesh)
parameters["ghost_mode"] = "shared_vertex"
mesh = Mesh()
with XDMFFile("cells.xdmf") as infile:
infile.read(mesh)
mvc_facets = MeshValueCollection("size_t", mesh, dim-1)
with XDMFFile("facets.xdmf") as infile:
infile.read(mvc_facets, "name_to_read")
mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)
mvc_cells = MeshValueCollection("size_t", mesh, dim)
with XDMFFile("cells.xdmf") as infile:
infile.read(mvc_cells, "name_to_read")
mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)
return mesh, mf_facets, mf_cells
gmsh.initialize(sys.argv)
gmsh.model.add("geometry")
gmsh.model.occ.addCircle(0, 0, 0, 0.005*um, tag=1)
gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.addPhysicalGroup(1, [1], tag=2)
gmsh.model.mesh.setSize([(0, 1)], size=0.0001*um)
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
mesh, mf_facets, mf_cells = msh_to_mf()
mesh.init_cell_orientations(Expression(("x[0]", "x[1]", "x[2]"), degree = 2))
wl0 = 0.5*um
dx = Measure("dx", domain=mesh, subdomain_data=mf_cells)
div_el = FiniteElement('N1div', mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, div_el*div_el)
(rP, iP) = TrialFunctions(V)
(rK, iK) = TestFunctions(V)
rPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
iPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
zero_normal_real = DirichletBC(V.sub(0), rPn, mf_facets, 2, method="geometric")
zero_normal_imag = DirichletBC(V.sub(1), iPn, mf_facets, 2, method="geometric")
bcs = [zero_normal_real, zero_normal_imag]
rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
P_term = -(dot(rP, rK))*dx(1) \
-(dot(iP, iK))*dx(1) \
-(dot(rP, iK))*dx(1) \
+(dot(iP, rK))*dx(1)
Eb_term = -(dot(rEb, rK))*dx(1) \
-(dot(iEb, iK))*dx(1)
equation = P_term + Eb_term
F = equation
a, L = lhs(F), rhs(F)
U = Function(V)
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
(rPm, iPm) = U.split()
```

As soon as I add the domain 2, I’m getting the error, as running the following script show. I think this is because the weak form does not have any information on my trial functions in the domain 2.

```
# DOMAIN 1 + DOMAIN 2 ---> Solution failed to converge in 0 iterations
from dolfin import *
import numpy as np
import gmsh
import sys
import meshio
um = 10**-6
def msh_to_mf():
dim = 2
cell_str = "triangle"
facet_str = "line"
msh = meshio.read("mesh.msh")
cells = np.vstack(np.array([cells.data for cells in msh.cells
if cells.type == cell_str]))
total_mesh = meshio.Mesh(points=msh.points,
cells=[(cell_str, cells)])
total_mesh.prune_z_0()
cells_data = msh.cell_data_dict["gmsh:physical"][cell_str]
cells_mesh = meshio.Mesh(points=msh.points,
cells=[(cell_str, cells)],
cell_data={"name_to_read": [cells_data]})
cells_mesh.prune_z_0()
facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
if cells.type == facet_str]))
facet_data = msh.cell_data_dict["gmsh:physical"][facet_str]
facet_mesh = meshio.Mesh(points=msh.points,
cells=[(facet_str, facet_cells)],
cell_data={"name_to_read": [facet_data]})
facet_mesh.prune_z_0()
meshio.write("total.xdmf", total_mesh)
meshio.write("facets.xdmf", facet_mesh)
meshio.write("cells.xdmf", cells_mesh)
parameters["ghost_mode"] = "shared_vertex"
mesh = Mesh()
with XDMFFile("cells.xdmf") as infile:
infile.read(mesh)
mvc_facets = MeshValueCollection("size_t", mesh, dim-1)
with XDMFFile("facets.xdmf") as infile:
infile.read(mvc_facets, "name_to_read")
mf_facets = cpp.mesh.MeshFunctionSizet(mesh, mvc_facets)
mvc_cells = MeshValueCollection("size_t", mesh, dim)
with XDMFFile("cells.xdmf") as infile:
infile.read(mvc_cells, "name_to_read")
mf_cells = cpp.mesh.MeshFunctionSizet(mesh, mvc_cells)
return mesh, mf_facets, mf_cells
gmsh.initialize(sys.argv)
gmsh.model.add("geometry")
gmsh.model.occ.addCircle(0, 0, 0, 0.005*um, tag=1)
gmsh.model.occ.addCircle(0, 0, 0, 0.05*um, tag=2) #added
gmsh.model.occ.addCurveLoop([1], tag=1)
gmsh.model.occ.addPlaneSurface([1], tag=1)
gmsh.model.occ.addCurveLoop([2], tag=2) #added
gmsh.model.occ.addCurveLoop([1], tag=3) #added
gmsh.model.occ.addPlaneSurface([2, 3], tag=2) #added
gmsh.model.occ.synchronize()
gmsh.model.addPhysicalGroup(2, [1], tag=1)
gmsh.model.addPhysicalGroup(1, [1], tag=2)
gmsh.model.addPhysicalGroup(2, [2], tag=3) #added
gmsh.model.mesh.setSize([(0, 1)], size=0.0001*um)
gmsh.model.mesh.setSize([(0, 2)], size=0.02*um) #added
gmsh.model.mesh.generate(2)
gmsh.write("mesh.msh")
mesh, mf_facets, mf_cells = msh_to_mf()
mesh.init_cell_orientations(Expression(("x[0]", "x[1]", "x[2]"), degree = 2))
wl0 = 0.5*um
dx = Measure("dx", domain=mesh, subdomain_data=mf_cells)
div_el = FiniteElement('N1div', mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, div_el*div_el)
(rP, iP) = TrialFunctions(V)
(rK, iK) = TestFunctions(V)
rPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
iPn = Expression(("0.0", "0.0"), domain=mesh, degree=2)
zero_normal_real = DirichletBC(V.sub(0), rPn, mf_facets, 2, method="geometric")
zero_normal_imag = DirichletBC(V.sub(1), iPn, mf_facets, 2, method="geometric")
bcs = [zero_normal_real, zero_normal_imag]
rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])"),pi=np.pi, wl=wl0, domain=mesh, degree=2)
P_term = -(dot(rP, rK))*dx(1) \
-(dot(iP, iK))*dx(1) \
-(dot(rP, iK))*dx(1) \
+(dot(iP, rK))*dx(1)
Eb_term = -(dot(rEb, rK))*dx(1) \
-(dot(iEb, iK))*dx(1)
equation = P_term + Eb_term
F = equation
a, L = lhs(F), rhs(F)
U = Function(V)
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
(rPm, iPm) = U.split()
```

The error is this one:

```
Traceback (most recent call last):
File "n1div_single_equation_not_working.py", line 125, in <module>
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 233, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solving.py", line 273, in _solve_varproblem
solver.solve()
RuntimeError:
*** -------------------------------------------------------------------------
*** 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 solve linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_PC_FAILED, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0
***
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
```

Therefore I’m asking again: how can I specify to `dolfin`

that my unknown is `= 0`

in the domain 2? Or is there any other way to make the simulation run? Thank you.