Hi all,

I wrote a post some days ago, receiving no answers. Probably my post was not clearly understandable by other users, and I apologize for that. I try to divide it in little bits in order to be more understandable to others.

I’m trying to solve a system of two PDEs where the unknowns are the imaginary and real part of two different vector fields. For the first one, I’m using `N1curl`

elements, while for the other one I am using `N1div`

elements.

When I try to run the code, I get this error:

```
Traceback (most recent call last):
File "n1div.py", line 96, 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: cell orientation must be defined (not -1)
```

This error goes away if I add this line after having loaded the `.msh`

file:

`mesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]'), degree = 2))`

I do not understand the role of the `init_cell_orientations()`

function. Why should I use that? And mainly, how should I use that? Does the `Expression`

inside this function depend on the geometry and/or the mesh in any way?

Thank you very much. If you need further details, I’m posting the code below.

```
//+
SetFactory("OpenCASCADE");
//+
Circle(1) = {0, 0, 0, 0.040, 0, 2*Pi};
//+
Circle(2) = {0, 0, 0, 0.05, 0, 2*Pi};
//+
Circle(3) = {0, 0, 0, 1.2, 0, 2*Pi};
//+
Circle(4) = {0, 0, 0, 1.8, 0, 2*Pi};
Curve Loop(1) = {1};
//+
Plane Surface(1) = {1};
//+
Curve Loop(2) = {2};
//+
Curve Loop(3) = {1};
//+
Plane Surface(2) = {2, 3};
//+
Curve Loop(4) = {3};
//+
Curve Loop(5) = {2};
//+
Plane Surface(3) = {4, 5};
//+
Curve Loop(6) = {4};
//+
Curve Loop(7) = {3};
//+
Plane Surface(4) = {6, 7};
//+
Physical Surface(1) = {2, 1};
//+
Physical Surface(2) = {3};
//+
Physical Surface(3) = {4};
//+
Physical Curve(4) = {2};
//+
Physical Curve(5) = {3};
//+
MeshSize {1} = 0.02;
//+
MeshSize {2} = 0.02;
//+
MeshSize {3} = 0.08;
//+
MeshSize {4} = 0.08;
//+
```

```
from dolfin import *
import numpy as np
import meshio
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],
dtype='object'))
total_mesh = meshio.Mesh(points=msh.points,
cells=[(cell_str, cells)])
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]})
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]})
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)
dx = Measure("dx", domain=mesh, subdomain_data=mf_cells)
curl_el = FiniteElement('N1curl', mesh.ufl_cell(), 2)
div_el = FiniteElement('N1div', mesh.ufl_cell(), 2)
# 4 elements since there are real and imaginary parts
element = MixedElement([curl_el, curl_el, div_el, div_el])
V = FunctionSpace(mesh, element)
W = TrialFunction(V)
# real and imaginary parts for the fields and test functions
rE, iE, rP, iP = split(W)
rV, iV, rK, iK = TestFunctions(V)
# real and imaginary DirichletBC
rPn = Expression(("0.0", "0.0", "0.0"), domain=mesh, degree=2)
iPn = Expression(("0.0", "0.0", "0.0"), domain=mesh, degree=2)
# DirichletBC
zero_normal_real = DirichletBC(V.sub(2), rPn, mf_facets, 4, method="geometric")
zero_normal_imag = DirichletBC(V.sub(3), iPn, mf_facets, 4, method="geometric")
bcs = [zero_normal_real, zero_normal_imag]
first_equation = inner(curl(rE), curl(rV))*(dx(1)+dx(2)) \
+ inner(curl(iE), curl(iV))*(dx(1)+dx(2))
second_equation = -inner(rP, rK)*dx(1) + \
-inner(iP, iK)*dx(1)
F = first_equation + second_equation
a, L = lhs(F), rhs(F)
U = Function(V)
solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
rEs, iEs, rPm, iPm = U.split()
```