RuntimeError: cell orientation must be defined (not -1)

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