Skip to content

Commit

Permalink
Instantiation of geometries from GDSII files (NanoComp#357)
Browse files Browse the repository at this point in the history
  • Loading branch information
Homer Reid authored and stevengj committed Jun 7, 2018
1 parent f34293e commit 61f32a7
Show file tree
Hide file tree
Showing 5 changed files with 283 additions and 22 deletions.
7 changes: 7 additions & 0 deletions configure.ac
Expand Up @@ -428,6 +428,13 @@ AM_CONDITIONAL(WITH_LIBCTLGEOM, test x"$have_libctlgeom" = "xyes")
LIBCTLGEOM_LIBS="-lctlgeom"
AC_SUBST(LIBCTLGEOM_LIBS)

##############################################################################
# check for libGDSII
AC_CHECK_HEADER(libGDSII.h, [have_gdsii=maybe], [have_gdsii=no])
if test $have_gdsii = maybe; then
AC_CHECK_LIB(GDSII, libGDSIIExists)
fi

##############################################################################
# The following function is used only for debugging. Note that
# we must test for it *after* setting the compiler flags (which
Expand Down
190 changes: 190 additions & 0 deletions libmeepgeom/GDSIIgeom.cpp
@@ -0,0 +1,190 @@
/* Copyright (C) 2005-2015 Massachusetts Institute of Technology
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2, or (at your option)
% any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details. %
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software Foundation,
% Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
*/

/***************************************************************/
/* GDSII.cpp -- libmeepgeom code to support meep geometry */
/* -- definitions from GDSII files */
/* homer reid -- 5/2018 */
/***************************************************************/

#include <vector>
#include <string>
#include "meepgeom.hpp"

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#ifdef HAVE_LIBGDSII
#include <libGDSII.h>
#endif

namespace meep_geom {

#ifdef HAVE_LIBGDSII

void get_polygon_bounding_box(dVec vertex_coordinates, meep::vec &max_corner, meep::vec &min_corner)
{
double xmax=vertex_coordinates[2*0+0], xmin=xmax;
double ymax=vertex_coordinates[2*0+1], ymin=ymax;
for(size_t nv=0; nv<vertex_coordinates.size()/2; nv++)
{ double x = vertex_coordinates[2*nv+0], y=vertex_coordinates[2*nv+1];
xmax=fmax(xmax,x); ymax=fmax(ymax,y);
xmin=fmin(xmin,x); ymin=fmin(ymin,y);
}
max_corner.set_direction(meep::X, xmax);
max_corner.set_direction(meep::Y, ymax);
min_corner.set_direction(meep::X, xmin);
min_corner.set_direction(meep::Y, ymin);
}

void get_polygon_center_size(dVec vertex_coordinates, meep::vec &center, meep::vec &size)
{
meep::vec max_corner, min_corner;
get_polygon_bounding_box(vertex_coordinates, max_corner, min_corner);

center.set_direction(meep::X, 0.5*(max_corner.in_direction(meep::X) + min_corner.in_direction(meep::X)));
center.set_direction(meep::Y, 0.5*(max_corner.in_direction(meep::Y) + min_corner.in_direction(meep::Y)));
center.set_direction(meep::Z, 0.0);

size.set_direction(meep::X, max_corner.in_direction(meep::X) - min_corner.in_direction(meep::X));
size.set_direction(meep::Y, max_corner.in_direction(meep::Y) - min_corner.in_direction(meep::Y));
size.set_direction(meep::Z, 0.0);
}

// Search the geometry for a polygon on a given layer containing
// (the reference point of) a given text label.
// If Text==NULL, find any polygon on the given layer.
// If Layer==-1, search all layers.
// If multiple matching polygons are found, choose one arbitrarily.
dVec get_polygon(const char *GDSIIFile, const char *Text, int Layer=-1)
{
PolygonList polygons = libGDSII::GetPolygons(GDSIIFile, Text, Layer);

char Description[100];
if (Text)
snprintf(Description,100,"with label %s",Text);
else
snprintf(Description,100,"on layer %i",Layer);

if ( polygons.size()==0 )
meep::abort("%s: found no polygons %s",GDSIIFile, Description);
if ( polygons.size()>1 )
fprintf(stderr,"warning: %s: found multiple polygons %s (choosing arbitrarily)\n",GDSIIFile,Description);

return polygons[0];
}

meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, const char *Text, int Layer, double zsize)
{
dVec polygon = get_polygon(GDSIIFile, Text, Layer);

meep::vec center, size;
get_polygon_center_size(polygon, center, size);

geometry_lattice.size.x=size.in_direction(meep::X);
geometry_lattice.size.y=size.in_direction(meep::Y);
geometry_lattice.size.z=zsize;
meep::grid_volume gv =
zsize==0.0 ? meep::vol2d(geometry_lattice.size.x, geometry_lattice.size.y, resolution)
: meep::vol3d(geometry_lattice.size.x, geometry_lattice.size.y, zsize, resolution);
gv.center_origin();
return gv;
}

meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, int Layer, double zsize)
{ return set_geometry_from_GDSII(resolution, GDSIIFile, 0, Layer, zsize); }

geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, const char *Text, int Layer, double zmin, double zmax)
{
dVec polygon = get_polygon(GDSIIFile, Text, Layer);

int num_vertices = polygon.size()/2;
vector3 *vertices = new vector3[num_vertices];
for(int nv=0; nv<num_vertices; nv++)
{ vertices[nv].x=polygon[2*nv+0];
vertices[nv].y=polygon[2*nv+1];
vertices[nv].z=zmin;
}

double height = zmax-zmin;
vector3 zaxis={0,0,1};
return make_prism(material, vertices, num_vertices, height, zaxis);
}

geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, int Layer, double zmin, double zmax)
{ return get_GDSII_prism(material, GDSIIFile, 0, Layer, zmin, zmax); }

meep::volume get_GDSII_volume(const char *GDSIIFile, const char *Text, int Layer, double zmin, double zmax)
{
dVec polygon = get_polygon(GDSIIFile, Text, Layer);
meep::ndim di = ((zmin==0.0 && zmax==0.0) ? meep::D2 : meep::D3);
di = meep::D2;
meep::vec max_corner=meep::zero_vec(di), min_corner=meep::zero_vec(di);
get_polygon_bounding_box(polygon, max_corner, min_corner);
max_corner.set_direction(meep::Z, zmax);
min_corner.set_direction(meep::Z, zmin);
return meep::volume(max_corner, min_corner);
}

meep::volume get_GDSII_volume(const char *GDSIIFile, int Layer, double zmin, double zmax)
{ return get_GDSII_volume(GDSIIFile, 0, Layer, zmin, zmax); }

/***************************************************************/
/* stubs for compilation without libGDSII **********************/
/***************************************************************/
#else // HAVE_LIBGDSII

void GDSIIError(const char *Routine)
{ meep::abort("Meep must be configured/compiled with libGDSII for %s",Routine); }

meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, const char *Text, int Layer, double zsize)
{ (void) resolution; (void) GDSIIFile; (void) Text; (void) Layer; (void) zsize;
GDSIIError("set_geometry_from_GDSII");
return meep::grid_volume();

}
meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, int Layer, double zsize)
{ (void) resolution; (void) GDSIIFile; (void) Layer; (void) zsize;
GDSIIError("set_geometry_from_GDSII");
return meep::grid_volume();
}
geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, const char *Text, int Layer, double zmin, double zmax)
{ (void) material; (void) GDSIIFile; (void) Text; (void) Layer; (void) zmin; (void) zmax;
GDSIIError("get_GDSII_prism");
vector3 center={0.0, 0.0, 0.0};
return make_sphere(0,center,0.0);
}
geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, int Layer, double zmin, double zmax)
{ (void) material; (void) GDSIIFile; (void) Layer; (void) zmin; (void) zmax;
GDSIIError("get_GDSII_prism");
vector3 center={0.0, 0.0, 0.0};
return make_sphere(0,center,0.0);
}
meep::volume get_GDSII_volume(const char *GDSIIFile, const char *Text, int Layer, double zmin, double zmax)
{ (void) GDSIIFile; (void) Text; (void) Layer; (void) zmin; (void) zmax;
GDSIIError("get_GDSII_volume");
return meep::volume(meep::vec());
}
meep::volume get_GDSII_volume(const char *GDSIIFile, int Layer, double zmin, double zmax)
{ (void) GDSIIFile; (void) Layer; (void) zmin; (void) zmax;
GDSIIError("get_GDSII_volume");
return meep::volume(meep::vec());
}

#endif // HAVE_LIBGDSII

} // namespace meep_geom
2 changes: 1 addition & 1 deletion libmeepgeom/Makefile.am
Expand Up @@ -9,7 +9,7 @@ pkginclude_HEADERS = meepgeom.hpp material_data.hpp
AM_CPPFLAGS = -I$(top_srcdir)/src -DDATADIR=\"$(srcdir)/\"

libmeepgeom_la_SOURCES = \
meepgeom.cpp meepgeom.hpp material_data.hpp
meepgeom.cpp GDSIIgeom.cpp meepgeom.hpp material_data.hpp

libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@
libmeepgeom_la_LIBADD = $(LIBMEEP)
Expand Down
96 changes: 75 additions & 21 deletions libmeepgeom/bend-flux-ll.cpp
Expand Up @@ -39,16 +39,54 @@ vector3 v3(double x, double y=0.0, double z=0.0)
double dummy_eps(const vec &) { return 1.0; }

/***************************************************************/
/* define geometry from GDSII file *****************************/
/***************************************************************/
#define GEOM_LAYER 0 // hard-coded indices of GDSII layers
#define STRAIGHT_WVG_LAYER 1 // on which various geometric entities are defined
#define BENT_WVG_LAYER 2
#define SOURCE_LAYER 3
#define RFLUX_LAYER 4
#define STRAIGHT_TFLUX_LAYER 5
#define BENT_TFLUX_LAYER 6

structure create_structure_from_GDSII(char *GDSIIFile, bool no_bend,
volume &vsrc, volume &vrefl, volume &vtrans)
{
// set computational cell
double dpml = 1.0;
double resolution = 10.0;
grid_volume gv=meep_geom::set_geometry_from_GDSII(resolution, GDSIIFile, GEOM_LAYER);
structure the_structure(gv, dummy_eps, pml(dpml));

// define waveguide
geometric_object objects[1];
meep_geom::material_type dielectric = meep_geom::make_dielectric(12.0);
int GDSIILayer = (no_bend ? STRAIGHT_WVG_LAYER : BENT_WVG_LAYER);
objects[0] = meep_geom::get_GDSII_prism(dielectric, GDSIIFile, GDSIILayer);
geometric_object_list g={ 1, objects };
meep_geom::set_materials_from_geometry(&the_structure, g);

// define volumes for source and flux-monitor regions
vsrc = meep_geom::get_GDSII_volume(GDSIIFile, SOURCE_LAYER);
vrefl = meep_geom::get_GDSII_volume(GDSIIFile, RFLUX_LAYER);
vtrans = meep_geom::get_GDSII_volume(GDSIIFile, (no_bend ? STRAIGHT_TFLUX_LAYER : BENT_TFLUX_LAYER));

return the_structure;

}

/***************************************************************/
/***************************************************************/
void bend_flux(bool no_bend, bool use_prisms)
/***************************************************************/
structure create_structure_by_hand(bool no_bend, bool use_prisms,
volume &vsrc, volume &vrefl, volume &vtrans)
{
double sx=16.0; // size of cell in X direction
double sy=32.0; // size of cell in Y direction
double pad=4.0; // padding distance between waveguide and cell edge
double w=1.0; // width of waveguide
double dpml=1.0; // PML thickness
double resolution=5;
double resolution=10.0;

// (set! geometry-lattice (make lattice (size sx sy no-size)))
// (set! pml-layers (list (make pml (thickness 1.0))))
Expand Down Expand Up @@ -137,6 +175,28 @@ void bend_flux(bool no_bend, bool use_prisms)
}
}

vsrc = volume( vec(-0.5*sx + dpml, wvg_ycen-0.5*w), vec(-0.5*sx + dpml, wvg_ycen+0.5*w) );
vrefl = volume( vec(-0.5*sx + 1.5, wvg_ycen-w ), vec(-0.5*sx+1.5, wvg_ycen+w ) );
vtrans = (
no_bend ? volume( vec(0.5*sx-1.5, wvg_ycen-w), vec(0.5*sx-1.5, wvg_ycen+w ) )
: volume( vec(wvg_xcen-w, 0.5*sy-1.5), vec(wvg_xcen+w, 0.5*sy-1.5 ) )
);

return the_structure;

}

/***************************************************************/
/***************************************************************/
/***************************************************************/
void bend_flux(bool no_bend, char *GDSIIFile, bool use_prisms)
{
vec v0;
volume vsrc(v0), vrefl(v0), vtrans(v0);
structure the_structure =
GDSIIFile ? create_structure_from_GDSII(GDSIIFile, no_bend, vsrc, vrefl, vtrans)
: create_structure_by_hand(no_bend, use_prisms, vsrc, vrefl, vtrans);

fields f(&the_structure);

// (set! sources (list
Expand All @@ -149,8 +209,7 @@ void bend_flux(bool no_bend, bool use_prisms)
double fcen = 0.15; // ; pulse center frequency
double df = 0.1; // ; df
gaussian_src_time src(fcen, df);
volume v( vec(-0.5*sx + dpml, wvg_ycen-0.5*w), vec(-0.5*sx + dpml,wvg_ycen+0.5*w) );
f.add_volume_source(Ez, src, v);
f.add_volume_source(Ez, src, vsrc);

//(define trans ; transmitted flux
// (add-flux fcen df nfreq
Expand All @@ -163,19 +222,15 @@ void bend_flux(bool no_bend, bool use_prisms)
double f_end = fcen+0.5*df;
int nfreq = 100; // number of frequencies at which to compute flux

volume *trans_volume=
no_bend ? new volume(vec(0.5*sx-1.5, wvg_ycen-w), vec(0.5*sx-1.5, wvg_ycen+w))
: new volume(vec(wvg_xcen-w, 0.5*sy-1.5), vec(wvg_xcen+w, 0.5*sy-1.5));
direction trans_dir = no_bend ? X : Y;
dft_flux trans=f.add_dft_flux(trans_dir, *trans_volume, f_start, f_end, nfreq);
dft_flux trans=f.add_dft_flux(trans_dir, vtrans, f_start, f_end, nfreq);

//(define refl ; reflected flux
// (add-flux fcen df nfreq
// (make flux-region
// (center (+ (* -0.5 sx) 1.5) wvg-ycen) (size 0 (* w 2)))))
//
volume refl_volume( vec(-0.5*sx+1.5, wvg_ycen-w), vec(-0.5*sx+1.5,wvg_ycen+w));
dft_flux refl=f.add_dft_flux(NO_DIRECTION, refl_volume, f_start, f_end, nfreq);
dft_flux refl=f.add_dft_flux(NO_DIRECTION, vrefl, f_start, f_end, nfreq);

char *prefix = const_cast<char *>(no_bend ? "straight" : "bent");
f.output_hdf5(Dielectric, f.total_volume(), 0, false, true, prefix);
Expand All @@ -195,12 +250,11 @@ void bend_flux(bool no_bend, bool use_prisms)
// (vector3 wvg-xcen (- (/ sy 2) 1.5)))
// 1e-3)
// (at-beginning output-epsilon))
vec eval_point = no_bend ? vec(0.5*sx-1.5, wvg_ycen) : vec(wvg_xcen, 0.5*sy - 1.5);
vec eval_point = vtrans.center();
double DeltaT=10.0, NextCheckTime = f.round_time() + DeltaT;
double Tol=1.0e-3;
double max_abs=0.0, cur_max=0.0;
bool Done=false;
/*
do
{
f.step();
Expand All @@ -218,13 +272,6 @@ void bend_flux(bool no_bend, bool use_prisms)

//printf("%.2e %.2e %.2e %.2e\n",f.round_time(),absEz,max_abs,cur_max);
} while(!Done);
*/

f.step();
master_printf("howdage! refl: \n");
f.output_dft(refl,"refl");
master_printf("howdage! trans: \n");
f.output_dft(trans,"trans");

//; for normalization run, save flux fields for refl. plane
//(if no-bend? (save-flux "refl-flux" refl))
Expand All @@ -249,12 +296,19 @@ int main(int argc, char *argv[])
initialize mpi(argc, argv);

bool use_prisms=false;
char *GDSIIFile=0;
for(int narg=1; narg<argc; narg++)
if (!strcasecmp(argv[narg],"--use-prisms"))
use_prisms=true;
for(int narg=1; narg<argc-1; narg++)
if (!strcasecmp(argv[narg],"--GDSIIFile"))
GDSIIFile=argv[narg+1];

if ( GDSIIFile!=0 && use_prisms==true )
fprintf(stderr,"warning: --use-prisms is ignored if --GDSIIFile is specified\n");

bend_flux(true, use_prisms);
bend_flux(false, use_prisms);
bend_flux(true, GDSIIFile, use_prisms);
bend_flux(false, GDSIIFile, use_prisms);

// success if we made it here
return 0;
Expand Down
10 changes: 10 additions & 0 deletions libmeepgeom/meepgeom.hpp
Expand Up @@ -122,6 +122,16 @@ bool is_medium(void* md, medium_struct **m);
bool is_metal(meep::field_type ft, const material_type *material);
void check_offdiag(medium_struct *m);

/***************************************************************/
/* routines in GDSIIgeom.cc ************************************/
/***************************************************************/
meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, const char *Text, int Layer=-1, double zsize=0.0);
meep::grid_volume set_geometry_from_GDSII(double resolution, const char *GDSIIFile, int Layer, double zsize=0.0);
geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, const char *Text, int Layer=-1, double zmin=0.0, double zmax=0.0);
geometric_object get_GDSII_prism(material_type material, const char *GDSIIFile, int Layer, double zmin=0.0, double zmax=0.0);
meep::volume get_GDSII_volume(const char *GDSIIFile, const char *Text, int Layer=-1, double zmin=0.0, double zmax=0.0);
meep::volume get_GDSII_volume(const char *GDSIIFile, int Layer, double zmin=0.0, double zmax=0.0);

}; // namespace meep_geom

#endif // #ifndef MEEP_GEOM_H

0 comments on commit 61f32a7

Please sign in to comment.