Hi all, a problem I am working on requires me to define inhomogeneous dirichlet boundary conditions in the tangential/perpendicular directions for Nedelec/Raviart-Thomas elements. The equations I am actually trying to solve are the Time-Dependent Ginzburg-Landau equations if you happen to be familiar with them, however for a MWE I will solve the Maxwell Equations using the formulation layed out here to demonstrate the issue. First, I use the following code to generate an irregularly spaced cube mesh with Gmsh (which I will compare behaviors to FEniCS’s BoxMesh function):
import numpy as np
import meshio
import dolfin as d
import os
char_length = 0.1
xmin = 0.0
ymin = 0.0
zmin = 0.0
xmax = 1.0
ymax = 1.0
zmax = 1.0
transfinite = False
verbose = True
if transfinite:
transf_str = ""
else:
transf_str = "//"
filename = "./MeshFiles/cubemesh_test"
gmsh_code = f"""// This code was generated by direct writing of a string to file
lc = {char_length};
Point(1) = {{{xmin}, {ymin}, {zmin}, lc}};
Point(2) = {{{xmax}, {ymin}, {zmin}, lc}};
Point(3) = {{{xmax}, {ymax}, {zmin}, lc}};
Point(4) = {{{xmin}, {ymax}, {zmin}, lc}};
Point(5) = {{{xmin}, {ymin}, {zmax}, lc}};
Point(6) = {{{xmax}, {ymin}, {zmax}, lc}};
Point(7) = {{{xmax}, {ymax}, {zmax}, lc}};
Point(8) = {{{xmin}, {ymax}, {zmax}, lc}};
Line(1) = {{1, 2}}; //bottom front
Line(2) = {{2, 3}}; //bottom right
Line(3) = {{3, 4}}; //bottom back
Line(4) = {{4, 1}}; //bottom left
Line(5) = {{1, 5}}; //vertical front left
Line(6) = {{2, 6}}; //vertical front right
Line(7) = {{3, 7}}; //vertical back right
Line(8) = {{4, 8}}; //vertical back left
Line(9) = {{5, 6}}; //top front
Line(10) = {{6, 7}}; //top right
Line(11) = {{7, 8}}; //top back
Line(12) = {{8, 5}}; //top left
Curve Loop(1) = {{-1,-2,-3,-4}}; Plane Surface(1) = {{1}}; //bottom
Curve Loop(2) = {{4,5,-12,-8}}; Plane Surface(2) = {{2}}; //left
Curve Loop(3) = {{2,7,-10,-6}}; Plane Surface(3) = {{3}}; //right
Curve Loop(4) = {{1,6,-9,-5}}; Plane Surface(4) = {{4}}; //front
Curve Loop(5) = {{3,8,-11,-7}}; Plane Surface(5) = {{5}}; //back
Curve Loop(6) = {{9,10,11,12}}; Plane Surface(6) = {{6}}; //top
Surface Loop(1) = {{1,2,4,3,5,6}};
Volume(1) = {{1}};
{transf_str}Transfinite Surface{{1}};
{transf_str}Transfinite Surface{{2}};
{transf_str}Transfinite Surface{{3}};
{transf_str}Transfinite Surface{{4}};
{transf_str}Transfinite Surface{{5}};
{transf_str}Transfinite Surface{{6}};
{transf_str}Transfinite Volume{{1}};
Physical Surface("boundary", 1) = {{1,2,3,4,5,6}};
Physical Volume("domain", 1) = {{1}};
"""
geo_file = open(f"{filename}.geo", "w")
geo_file.write(gmsh_code)
geo_file.close()
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
out_mesh = meshio.Mesh(points=mesh.points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
if prune_z:
out_mesh.prune_z_0()
return out_mesh
if verbose:
os.system(f'gmsh -3 {filename}.geo -o {filename}.msh')
else:
os.system(f'gmsh -3 {filename}.geo -o {filename}.msh -v 0')
meshdivot_full = meshio.read(f"{filename}.msh")
tetra_mesh = create_mesh(meshdivot_full, "tetra")
triangle_mesh = create_mesh(meshdivot_full, "triangle")
meshio.write(f"{filename}.xdmf", tetra_mesh)
meshio.write(f"{filename}_mf.xdmf", triangle_mesh)
This will create the mesh that I will be solving here. Now here is the code I use to solve Maxwell’s equations, I assume perfect conducting boundary conditions on the faces with normal directions in the x and y directions, and just for this toy problem I am applying a function of H_a\sin(\pi x)\sin(\pi y)\hat{z}(in the z direction because it’s perpendicular to the boundary and the magnetic field uses Raviart-Thomas elements) to the magnetic field on the z=0 and z=1 faces, where I have set H_a = 0.5(1-e^{-10t}) (an exponential ramp up to 0.5). The code I solve is the following:
from dolfin import *
import numpy as np
class boundaryXY(SubDomain):
def inside(self, x, on_boundary):
"""Returns true if the the point located at x is not near the top or bottom of the mesh"""
return on_boundary and not (near(x[2], 0.0,eps=1E-12) or near(x[2], 1.0,eps=1E-12))
class boundaryZ(SubDomain):
def inside(self, x, on_boundary):
"""Returns true if the the point located at x is near the top or bottom of the mesh"""
return on_boundary and (near(x[2], 0.0,eps=1E-12) or near(x[2], 1.0,eps=1E-12))
class HA(UserExpression):
def __init__(self, ha, degree, *args, **kwargs):
super(HA, self).__init__(*args, **kwargs)
self.haCurrent = 0 # Currently applied field for saving
self.ha = ha
def updateAppliedField(self, time = 0):
"""Updates the applied field to a given time."""
self.haCurrent = self.ha*(1 - np.exp(-time*10.0))
def eval(self, values, x):
"""Return the current applied field"""
values[0]= 0.0
values[1]= 0.0
values[2]= self.haCurrent*np.sin(np.pi*x[0])*np.sin(np.pi*x[1])
def boundary(x, on_boundary):
return on_boundary
def solveMaxwell(mesh, file, T, num_steps, Deg, Ha, save_every):
"""
Solve's Maxwell's Equations for a waveguide oriented in the z-direction
Based on a formulation from a paper by Asad Anees and Lutz Angermann (https://ieeexplore.ieee.org/document/8713469)
mesh: mesh to solve on
file: pvd file to save data to
T: total time to solve (starting from 0)
num_steps: steps to take, will be linearly spaced
Deg: degree of function spaces, at least 2 is reccomended for better solutions
Ha: Applied field as a DirichletBC on the open end of the waveguide
save_every: number of steps between solution saves, set to 1 to save every step
"""
dt = T/num_steps
dbcXY = boundaryXY()
dbcZ = boundaryZ()
V = FunctionSpace(mesh, "N1curl", Deg) #Function space for the electric field
Q = FunctionSpace(mesh, "RT", Deg) #Function space for the magnetic field
W = FunctionSpace(mesh, V.ufl_element() * Q.ufl_element()) #Mixed function space
Happlied = HA(Ha, degree=Deg, domain=Q.ufl_domain(), element = Q.ufl_element())
dU = TrialFunction(W)
E, H = split(dU)
psi, phi = TestFunctions(W)
u = Function(W)
u_n = Function(W)
E0, H0 = split(u_n)
F = (inner(E, psi) - inner(E0, psi) - dt*inner(H, curl(psi)) + inner(H, phi) - inner(H0, phi) + dt*inner(curl(E), phi))*dx
a = lhs(F)
L = rhs(F)
Ebc = DirichletBC(W.sub(0), Constant((0.0,0.0,0.0)), dbcXY)
Hbc1 = DirichletBC(W.sub(1), Constant((0.0,0.0,0.0)), dbcXY)
t = 0
for n in range(num_steps):
t += dt
Happlied.updateAppliedField(time=t)
Hbc2 = DirichletBC(W.sub(1), Happlied, dbcZ)
bcs = [Ebc,Hbc1,Hbc2]
solve(a == L, u, bcs, solver_parameters={'linear_solver':'gmres'})
u_n.assign(u)
if n % save_every == 0:
(E, H) = u.split(True)
file << (E, t)
file << (H, t)
#Now we set parameters and run the solves
T = np.pi/4
num_steps = 100
dim = 10
Deg = 2
Ha = 0.5
save_every = 10
loc = "./TestingFiles/" #folder to save data into
meshfile = "./MeshFiles/cubemesh_test.xdmf"
FEniCSmesh = BoxMesh(Point(0,0,0), Point(1,1,1), dim, dim, dim)
FEniCSfile = File(loc + f"Maxwell_fenics_T-{T}_dim-{dim}_Deg-{Deg}.pvd")
solveMaxwell(FEniCSmesh, FEniCSfile, T, num_steps, Deg, Ha, save_every)
Gmshmesh = Mesh()
with XDMFFile(meshfile) as infile:
infile.read(Gmshmesh)
Gmshfile = File(loc + f"Maxwell_gmsh_T-{T}_dim-{dim}_Deg-{Deg}.pvd")
solveMaxwell(Gmshmesh, Gmshfile, T, num_steps, Deg, Ha, save_every)
Solving on the mesh generated with BoxMesh gives pretty much what I would expect:
However, the mesh constructed with gmsh does not give a remotely comparable result.
I believe that the issue is somehow with how applying inhomogeneous boundary conditions work for irregularly spaced mesh elements, I ultimately want to solve on different shaped meshes constructed by gmsh so I can’t just settle for using BoxMesh but if this issue is not a bug but is instead fundamental to how Raviart-Thomas/Nedelec elements work, then I may just need to find another formulation which does not use them.
Some things I have considered which I do not think are the issue here (though I could always be wrong):
It does not seem to be an issue with Gmsh, as if I take the mesh and make it transfinite, the solution looks like the first picture, so the issue seems to be with the irregularity in this case.
Marking the boundaries with Gmsh and creating a meshfunction to apply the boundary conditions does not change anything here, so I have not taken that approach just to keep the code more concise
I am aware that the problem setup here for Maxwell’s equations is not necessarily physical, however I am more concerned with the fact that these two otherwise identical simulations are giving wildly different results, so I have chosen these boundary conditions just to demonstrate the issue
I will admit that my understanding of finite element methods beyond Lagrange elements is a little bit black box, so it may be the case that RT/Nedelec elements should not be used on tetrahedral elements that are not evenly spaced in a grid like they are for BoxMesh in which case I need to find another formulation, but I want to make sure that this isn’t just a bug or I have messed up my implementation here first. If anyone else has encountered similar issues or has an idea what is going on here, let me know.
Thanks in advance for your time!