I am trying to solve a test problem related to the concept of electromagnetics. In this problem, I aimed at solving the following expression,
\nabla \times (\frac{1}{\mu} \nabla \times E) + \epsilon_r \dfrac{\partial^2 E}{\partial t^2}= - \dfrac{\partial J}{\partial t}
To obtain the electric vector field. For this purpose, I introduce a “Nedelec 1st kind H(curl)”. In the case of the structured mesh, the results, which ParaView views, have a homogenous and uniform distribution. However, in the case of the unstructured mesh, oscillations are observed in the post-process. Does anybody know what to do to obtain a uniform distribution of the electric field when the domain has an unstructured mesh?
You can find the code attached below:
from dolfin import *
from mshr import *
import matplotlib.pyplot as plt
T_max = 2.
nSteps = 100
dTime = T_max/nSteps
Geometry = Rectangle(Point(-2.,-2.), Point(2.,2.))
mesh = generate_mesh(Geometry,50)
# mesh = RectangleMesh(Point(-2.,-2.), Point(2.,2.), 50, 50)
subdomain = CompiledSubDomain("x[1] >=-0.2 && x[1] <= 0.2 && x[0] >= -1.2 && x[0] <= 1.2")
markers = MeshFunction("size_t", mesh, mesh.topology().dim(), mesh.domains())
dx = Measure('dx', domain=mesh, subdomain_data=markers)
markers.set_all(0)
subdomain.mark(markers, 1)
plot(markers)
plt.savefig("mesh.png")
###
mesh_data = open('Structured_Mesh.data', 'w')
numel = mesh.num_cells()
mesh_hmax = mesh.hmax()
mesh_hmin = mesh.hmin()
print("Number_of_Elements=", mesh.num_cells())
print("Max_Element_Size=",mesh.hmax())
print("Min_Element_Size=",mesh.hmin())
mesh_data.write( \
"Numer of Elements=%i\n Max. Element Size=%e\n Min. Element Size=%e" \
%(numel, mesh_hmax, mesh_hmin))
mesh_data.close()
###
# Define function spaces
V = FunctionSpace(mesh, "N1curl", 2)
##V = VectorFunctionSpace(mesh, "P", 1)
##W = VectorFunctionSpace(mesh, "DG", 1)
# Define test and trial functions
v = TestFunction(V)
u = TrialFunction(V)
zero = Constant((0.0,0.0))
u_n =Function(V)
u_n_1 =Function(V)
u_n = interpolate(zero, V)
u_n_1 = interpolate(zero, V)
TOL = 1E-5
# Define functions for boundary condiitons
def LeftBoundary(x, on_boundary):
return near(x[0],-2.,TOL) and on_boundary
def RightBoundary(x, on_boundary):
return near(x[0],+2.,TOL) and on_boundary
def TopBoundary(x, on_boundary):
return near(x[1],-2.,TOL) and on_boundary
def BottomBoundary(x, on_boundary):
return near(x[1],+2.,TOL) and on_boundary
BC_Left = DirichletBC(V, zero, LeftBoundary)
BC_Right = DirichletBC(V, zero, RightBoundary)
BC_Bottom = DirichletBC(V, zero, TopBoundary)
BC_Top = DirichletBC(V, zero, BottomBoundary)
# class DirichletBoundary(SubDomain):
# def inside(self, x, on_boundary):
# return on_boundary
BCs = [BC_Left, BC_Right, BC_Top, BC_Bottom]
## BCs = DirichletBC(V, zero, DirichletBoundary())
Jt = Expression(("x[0]*t","x[1]*t"), degree=2, t=0)
# Define magnetic permeability
class Permittivity(UserExpression):
def __init__(self, markers, **kwargs):
super(Permittivity, self).__init__(**kwargs)
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 0:
values[0] = 1.0 # vacuum
elif self.markers[cell.index] == 1:
values[0] = 3.5 # iron (should really be 6.3e-3)
eps_r = Permittivity(markers, degree=1)
mu = 1.
#
a = dTime*dTime*inner(curl(v), curl(u)/mu)*dx + dot(v,eps_r*u)*dx
L = 2.*dot(eps_r*u_n,v)*dx - dot(eps_r*u_n_1,v)*dx + dTime*dTime*dot(v,Jt)*dx(1)
E_Field = Function(V)
VTKFile = File("./Results/Electric_Field.pvd")
Time = 0.
for n in range(nSteps):
Time+=dTime
Jt.t = Time
solve(a == L, E_Field, BCs)
u_n_1.assign(u_n)
u_n.assign(E_Field)
E_Field.rename("Electric_Field","1")
VTKFile << E_Field
