Hi all,
I am trying to write a basic module to allow for simple 3D shapes in
MEEP (for now, polyhedra). I am most interested in rectangular pyramids.
I have two main questions:
* first, is this an appropriate way to go about writing pyramid
structures? All I did was make a class for the Pyramid (in my actual
code, it inherits a Shape class), along with a material. Then I
convert a list of these objects with MEEP's GeometricObject to a
material_func for my Simulation object. I'm running into performance
problems, especially for subpixel averaging. Is there a better way
to go about this? It would be nice to just make a more general
GeometricObject (say, a polyhedron defined by a set of planes) but I
don't know how to do this.
* I have subpixel averaging disabled in my code by using
eps_average=False. I thought this would fix my problem, but it does
not. Even with this extra line if I run the code it still gives
"subpixel-averaging is 97.7767% done, 0.0910229 s remaining" etc.
Why is it still doing this, and how should I get around this?
I scraped the rest of the code and left a bare-bones example in the .py
file attached. The only output is some 2D cross-sectional plots of the
pyramid in question. Output seems to be right but it takes a long time
(and this is just a minimal example -- more realistically I'll have more
shapes and a much higher resolution).
Best regards,
Philip
"""
:author JPHaupt
minimal example
the actual code in use is in a module
"""
# %% imports
import meep as mp
from meep import Vector3 as V # lazy shorthand
from meep.materials import cSi
from matplotlib import pyplot as plt
# %% function/class definitions
# these are actually in a separate module
class Plane(object):
"""
a plane parametrised by a normal vector and a point in the plane, useful for
polyhedron shapes.
"""
def __init__(self, normal, point):
self.n = normal
self.p = point
def sameSide(self, point):
v = point - self.p
proj = v.cdot(self.n)
return (proj > 0)
class Pyramid(object):
def __init__(self, apex, b1, b2, b3, b4, material=mp.Medium()):
"""
initialises a rectangular pyramid object with points on the pyramid
given by vectors b1,b2,b3,b4 and the apex vector
geometric object for a pyramid with some medium and consisting of points
b1,b2,b3,b4 (on the rectangular base) and apex (the apex). Supports
oblique rectangular pyramids as well
:param Shape: [description]
:type Shape: [type]
:param apex: the apex of the pyramid
:type apex: mp.Vector3
:param b1: one of the vertices for the rectangular base of the pyramid
:type b1: mp.Vector3
:param b2: one of the vertices for the rectangular base of the pyramid
:type b2: mp.Vector3
:param b3: one of the vertices for the rectangular base of the pyramid
:type b3: mp.Vector3
:param b4: one of the vertices for the rectangular base of the pyramid
:type b4: mp.Vector3
.. note:: currently, we require that the points on the base of the
pyramid be given counterclockwise relative to outside (i.e. side
**without** the apex).
"""
self.material = material
# TODO implement such that we don't need to add them counterclockwise
self.planes = []
# make base normal (pointing outside), partner with arbitrary point
base_normal = (b2-b1).cross(b3-b1)
self.planes.append(Plane(base_normal, b1))
# make triangular sides' normals, partner with (apex) point
n1 = (apex-b1).cross(b2-b1)
# n1 /= n1.norm()
self.planes.append(Plane(n1, apex))
n2 = (apex-b2).cross(b3-b2)
# n2 /= n2.norm()
self.planes.append(Plane(n2, apex))
n3 = (apex-b3).cross(b4-b3)
# n3 /= n3.norm()
self.planes.append(Plane(n3, apex))
n4 = (apex-b4).cross(b1-b4)
# n4 /= n4.norm()
self.planes.append(Plane(n4, apex))
# do the rest, i.e. material
# super(Pyramid, self).__init__(**kwargs)
# override the contains function
def contains(self, point):
"""
checks to see if the point is in the pyramid or not
this is achieved by taking the projection of the vector connecting a
point on a given plane to the point in question onto the normal vector
of that same plane. If all of these projections are nonpositive, then
the point is inside the pyramid. If any one of them is positive, then it
is outside the pyramid.
:param point: a point to be checked if inside or outside the pyramid
:type point: mp.Vector3
"""
for plane in self.planes:
if plane.sameSide(point):
return False
return True
def geom2med(geometries, default_material=mp.Medium(epsilon=1.0)):
"""
takes a mixed list of mp.GeometricObject and Shape (defined in this module)
and returns a function of mp.Vector3 that returns the medium at the point
:param geometries: list of geometric objects/shapes
:type geometries: [mp.GeometricObject or Shape]
:default_material: the material to be chosen if point not in any of the
objects. Default is air (epsilon=1.0)
:type default_material: mp.Medium
.. note:: objects at the end of the list are checked first (i.e. if two
objects contain the point, then it is considered to be in the last
one).
"""
return lambda v: _geom2med_helper(geometries, v, default_material)
def _geom2med_helper(geometries, point, default_material):
"""
Function that geom2med wraps around
"""
for geom in reversed(geometries):
if isinstance(geom, mp.GeometricObject):
if mp.is_point_in_object(point, geom):
return geom.material
# elif isinstance(geom, shapes.Shape):
elif isinstance(geom, Pyramid):
if geom.contains(point):
return geom.material
return default_material
# %% actual program
dpml = 0.5 # PML thickness
k = mp.Vector3(0,0,0) # makes PBC
pml_layers = [mp.PML(dpml, direction=mp.Z, side=mp.ALL)]
# domain size
sx = 5
sy = 5
sz = 5
L = 2.5
resolution = 30
# recall: **must** put in anticlockwise order!
b1 = V(L/2, L/2, 0)
b2 = V(-L/2, L/2, 0)
b3 = V(-L/2, -L/2, 0)
b4 = V(L/2, -L/2, 0)
apex = V(0,0,-L/3)
tmp_pyr = Pyramid(apex, b1, b2, b3, b4, material=cSi)
# wavelengths
# 400-1100 nm => 0.4-1.1 um
wvl_min = 0.3 # minimum wavelength
wvl_max = 1.2 # maximum wavelength
freq_min = 1/wvl_max
freq_max = 1/wvl_min
freq_centre = (freq_max+freq_min)/2
freq_width = freq_max-freq_min
tmp_source = [mp.Source(mp.GaussianSource(frequency=freq_centre, width=freq_width),
component=mp.Ex,
center=mp.Vector3(0, 0, 0.5*sz-dpml-0.25),
# size=mp.Vector3(sx-2*pml_width,sy-2*pml_width,0))]
size=mp.Vector3(sx,sy,0))]
sim = mp.Simulation(cell_size=V(sx,sy,sz),
boundary_layers=pml_layers,
resolution=resolution,
k_point=k,
sources=tmp_source,
material_function=geom2med([tmp_pyr]),
eps_averaging=False,
extra_materials=[cSi])
# visualise
num = 4
print('num=',num)
plt.figure()
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(sx,sy,0)))
plt.savefig('tmp_xy_' + str(num) + '.png')
plt.figure()
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(sx,0,sz)))
plt.savefig('tmp_xz_' + str(num) + '.png')
plt.figure()
sim.plot2D(output_plane=mp.Volume(center=mp.Vector3(0,0,0), size=mp.Vector3(0,sy,sz)))
plt.savefig('tmp_yz_' + str(num) + '.png')_______________________________________________
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
