em_simulation.py

from dolfin import *
import meshio
import numpy as np
import gmsh
import sys

wl = 0.4

r_cyl = 0.025
r_dom = 0.5

#--------- Start Mesh ---------#

gmsh.initialize(sys.argv)
if MPI.comm_world.rank == 0:
    gmsh.model.add("geometry")
    gmsh.model.occ.addCircle(0, 0, 0, r_cyl, angle1 = 0, angle2 = pi*2, tag=1)
    gmsh.model.occ.addCircle(0, 0, 0, r_dom, angle1 = 0, angle2 = pi*2, tag=2)
    gmsh.model.occ.addCurveLoop([1], tag=1)
    gmsh.model.occ.addCurveLoop([2], tag=2)
    gmsh.model.occ.addPlaneSurface([1], tag=1)
    gmsh.model.occ.addPlaneSurface([2, 1], tag=2)
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(2, [1], tag=1)
    gmsh.model.addPhysicalGroup(2, [2], tag=2)
    gmsh.model.addPhysicalGroup(1, [2], tag=3)
    gmsh.model.mesh.setSize([(0, 1)], size = 0.0001)
    gmsh.model.mesh.setSize([(0, 2)], size = 0.01)
    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")
gmsh.finalize()
if MPI.comm_world.rank == 0:
    msh = meshio.read("mesh.msh")

    cells = np.vstack(np.array([cells.data for cells in msh.cells
                                if cells.type=="triangle"]))

    cells_mesh = meshio.Mesh(points=msh.points,
                                cells=[("triangle", cells)])

    triangle_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                    if cells.type=="triangle"])) 

    triangle_data = msh.cell_data_dict["gmsh:physical"]["triangle"]

    triangle_mesh = meshio.Mesh(points=msh.points,
                                cells=[("triangle", triangle_cells)],
                                cell_data={"name_to_read": [triangle_data]})

    facet_cells = np.vstack(np.array([cells.data for cells in msh.cells
                                    if cells.type=="line"]))

    facet_data = msh.cell_data_dict["gmsh:physical"]["line"]

    facet_mesh = meshio.Mesh(points=msh.points,
                            cells=[("line", facet_cells)],
                            cell_data={"name_to_read": [facet_data]})

    meshio.write("total.xdmf", cells_mesh)
    meshio.write("boundaries.xdmf", facet_mesh)
    meshio.write("volume.xdmf", triangle_mesh)
MPI.comm_world.barrier()

parameters["ghost_mode"] = "shared_vertex"
comm = MPI.comm_world
mesh = Mesh(MPI.comm_world)

with XDMFFile(MPI.comm_world, "volume.xdmf") as infile:
    infile.read(mesh)

mvc_2d = MeshValueCollection("size_t", mesh, 1)

with XDMFFile(MPI.comm_world, "boundaries.xdmf") as infile:
    infile.read(mvc_2d, "name_to_read")

mf_2d = cpp.mesh.MeshFunctionSizet(mesh, mvc_2d)

mvc_3d = MeshValueCollection("size_t", mesh, 2)

with XDMFFile(MPI.comm_world, "volume.xdmf") as infile:
    infile.read(mvc_3d, "name_to_read")

mf_3d = cpp.mesh.MeshFunctionSizet(mesh, mvc_3d)

#--------- End Mesh ---------#

ds = Measure("ds", domain=mesh, subdomain_data=mf_2d)
dx = Measure("dx", domain=mesh, subdomain_data=mf_3d)

n = FacetNormal(mesh)

curl_el = FiniteElement('N1curl', mesh.ufl_cell(), 2)
V = FunctionSpace(mesh, curl_el * curl_el) 

rEb = Expression(("0.0", "cos(-2*pi/wl*x[0])", "0.0"),pi=np.pi, wl=wl, domain=mesh, degree=2)
iEb = Expression(("0.0", "sin(-2*pi/wl*x[0])", "0.0"),pi=np.pi, wl=wl, domain=mesh, degree=2)

k0 = Constant(2*np.pi/wl)

reps = Constant(-1.078245)
ieps = Constant(5.808852)

bcs = []

(rEs, iEs) = TrialFunctions(V)
(rV, iV) = TestFunctions(V)

a =  inner(curl(rEs), curl(rV))*dx - k0**2*inner(rEs,rV)*dx(2) - reps*k0**2*inner(rEs,rV)*dx(1) + ieps*k0**2*inner(iEs,rV)*dx(1) \
    +inner(curl(iEs), curl(iV))*dx - k0**2*inner(iEs,iV)*dx(2) - reps*k0**2*inner(iEs,iV)*dx(1) - ieps*k0**2*inner(rEs,iV)*dx(1) \
    -inner(k0*cross(iEs,n),cross(rV,n))*ds(3) \
    +inner(k0*cross(rEs,n),cross(iV,n))*ds(3)

L = -(-(reps-1)*k0**2*dot(rEb,rV) + ieps*k0**2*dot(iEb,rV))*dx(1) \
    -(-(reps-1)*k0**2*dot(iEb,iV) - ieps*k0**2*dot(rEb,iV))*dx(1)

U = Function(V)

solve(a == L, U, bcs, solver_parameters={'linear_solver':'mumps'})
(real_Es, imag_Es) = U.split()

list_timings(TimingClear.clear, [TimingType.wall])