Hi all,
I’m working on a 2D electromagnetic problem described by a system of two PDEs, the unknowns are the electric field \mathbf{E} (with in-plane polarized) and the polarization density field \mathbf{P}. For the electric field, I’m using N1curl elements. For the polarization density, I am using the N1div elements, since I need to impose this condition on one of the boundaries:
\mathbf{P}\cdot\mathbf{n} = 0
\mathbf{n} is the normal with respect to the boundary.
In order to do so, I’ve written the code below, where the first code represents the .geo file for the mesh generation, the second one is the part of the code for loading the mesh in python, and the third one defines and solves the problem. The second and third scripts are part of a single file, here they’re split only for visualization purpose. Please notice that the problem is not complete, in the sense that the weak form below only has few terms (therefore it is for sure ill-posed) but it is anyway a good MWE to show what is my problem.
My questions are the following ones:
-
Is it correct to use
N1divelements when you have boundary conditions of the form \mathbf{P}\cdot\mathbf{n} = 0? In other words, with this kind of elements, does theDirichletBCfunction directly acts on the normal component of the field, as it correspondingly happens for theN1curlelements and the tangential component of the fields? -
Is the workflow for solving a system of equations defined below correct? This is the first time I do it.
-
My main problem: if I try to run this code, I get the error that is shown below. If I insert the line
mesh.init_cell_orientations(Expression(('x[0]', 'x[1]', 'x[2]'), degree = 2))after having loaded the mesh, this error goes away, but I applied it blindly, without knowing if it is the correct to do it. Therefore I’d like to ask: why is this error appearing? What is the correct way of using theinit_cell_orientations()function in order to avoid it?
I’m using Ubuntu 20.04.2 LTS and fenics 2019.1.0.
# Error I get
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)
# .geo file
//+
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;
//+
# First part of .py file
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)
# Second part of .py file
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()