Meep geom (#56)

* trying mkdocs stuff on master branch as readthedocs seems to have trouble with other branches

* synced with stevengj/master

* add missing line to configure.ac

* add Makefile.am and other missing files

* added libmeep_geom with three working tests

* renamed meep_geom -> meepgeom; removed -std=c++11 compiler flag

* change .travis.yml to pull from libctl master instead of specific commit

* further updates to remove c++11 syntax

* removed bend-flux-ll.cpp test from libmeepgeom

* static initialization for vacuum material_type structure

.travis.yml

@@ -12,7 +12,7 @@ addons:
     - swig
 script:
   - git clone https://github.com/stevengj/libctl libctl-src
-  - (cd libctl-src && git checkout 7d51710e63952d343ce8f241cf92b5cbfd5f1ada && sh autogen.sh --prefix=$HOME/local --enable-shared && make && make install)
+  - (cd libctl-src && git checkout master && sh autogen.sh --prefix=$HOME/local --enable-shared && make && make install)
   - git clone https://github.com/stevengj/harminv
   - (cd harminv && git checkout c221b2bcbaaa761f683aa5e2c6fa7efbbecdca1f && sh autogen.sh --prefix=$HOME/local --enable-shared && make && make install)
   - (CPPFLAGS=-I${HOME}/local/include LDFLAGS=-L${HOME}/local/lib GEN_CTL_IO=${HOME}/local/bin/gen-ctl-io sh autogen.sh --verbose --enable-maintainer-mode --prefix=$HOME/local --with-libctl=$HOME/local/share/libctl)

Makefile.am

@@ -6,7 +6,7 @@ endif
 
 ACLOCAL_AMFLAGS=-I m4
 
-SUBDIRS = src $(LIBCTL) tests examples
+SUBDIRS = src $(LIBCTL) tests examples libmeepgeom 
 
 if WITH_PYTHON
   SUBDIRS += python

configure.ac

@@ -513,6 +513,7 @@ AC_CONFIG_FILES([
 	examples/Makefile
 	libctl/Makefile
 	libctl/meep.scm
+	libmeepgeom/Makefile
 	python/Makefile
 ])
 

index.md

@@ -0,0 +1 @@
+Hello, this is index.

libmeepgeom/Makefile.am

@@ -0,0 +1,28 @@
+LIBCTLGEOM   = -lctlgeom
+LIBMEEP      = $(top_builddir)/src/libmeep.la
+MEEPLIBS     = $(LIBCTLGEOM) $(LIBMEEP)
+
+lib_LTLIBRARIES = libmeepgeom.la
+pkginclude_HEADERS = meepgeom.hpp material_data.hpp
+
+#AM_CPPFLAGS = -I$(top_srcdir)/src -std=c++11
+AM_CPPFLAGS = -I$(top_srcdir)/src
+
+libmeepgeom_la_SOURCES = \
+ meepgeom.cpp meepgeom.hpp material_data.hpp
+
+libmeepgeom_la_LDFLAGS = -version-info @SHARED_VERSION_INFO@
+
+check_PROGRAMS      = pw-source-ll ring-ll
+
+pw_source_ll_SOURCES = pw-source-ll.cpp
+pw_source_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
+
+ring_ll_SOURCES = ring-ll.cpp
+ring_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)
+
+TESTS               = pw-source-ll ring-ll
+
+noinst_PROGRAMS = bend-flux-ll
+bend_flux_ll_SOURCES = bend-flux-ll.cpp
+bend_flux_ll_LDADD   = libmeepgeom.la $(MEEPLIBS)

libmeepgeom/bend-flux-ll.cpp

@@ -0,0 +1,218 @@
+/***************************************************************/
+/* C++ port of meep/examples/bend-flux.ctl, using the          */
+/* "low-level" meep C++ interface stack, which consists of     */
+/* libmeep_geom + libctlgeom + libmeep                         */
+/***************************************************************/
+
+/*
+; From the Meep tutorial: transmission around a 90-degree waveguide
+; bend in 2d.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <complex>
+
+#include "meep.hpp"
+
+#include "ctl-math.h"
+#include "ctlgeom.h"
+
+#include "meepgeom.hpp"
+
+using namespace meep;
+
+typedef std::complex<double> cdouble;
+
+vector3 v3(double x, double y=0.0, double z=0.0)
+{
+  vector3 v;
+  v.x=x; v.y=y; v.z=z;
+  return v;
+}
+
+/***************************************************************/
+/* dummy material function needed to pass to structure( )      */
+/* constructor as a placeholder before we can call             */
+/* set_materials_from_geometry                                 */
+/***************************************************************/
+double dummy_eps(const vec &) { return 1.0; }
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+void bend_flux(bool no_bend)
+{
+  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 resolution=10; // (set-param! resolution 10)
+
+  // (set! geometry-lattice (make lattice (size sx sy no-size)))
+  // (set! pml-layers (list (make pml (thickness 1.0))))
+  grid_volume gv = voltwo(sx, sy, resolution);
+  gv.center_origin();
+  structure the_structure(gv, dummy_eps, pml(1.0));
+
+  double wvg_ycen = -0.5*(sy - w - 2.0*pad); //y center of horiz. wvg
+  double wvg_xcen =  0.5*(sx - w - 2.0*pad); //x center of vert. wvg
+
+   //  (set! geometry
+   //      (if no-bend?
+   //	  (list
+   //	  (list
+   //	   (make block 
+   //	     (center wvg-xcen (* 0.5 pad))
+   //	     (size w (- sy pad) infinity)
+   //	     (material (make dielectric (epsilon 12)))))))
+   //	   (make block 
+   //	     (center 0 wvg-ycen)
+   //	     (size infinity w infinity)
+   //	     (material (make dielectric (epsilon 12)))))
+   //	   (make block 
+   //	     (center (* -0.5 pad) wvg-ycen)
+   //	     (size (- sx pad) w infinity)
+   //	     (material (make dielectric (epsilon 12))))
+   vector3 e1=v3(1.0, 0.0, 0.0);
+   vector3 e2=v3(0.0, 1.0, 0.0);
+   vector3 e3=v3(0.0, 0.0, 1.0);
+
+   material_type dielectric = meep_geom::make_dielectric(12.0);
+   if (no_bend)
+    {
+      geometric_object objects[1];
+      objects[0] = make_block(dielectric,
+                              v3(0.0, wvg_ycen),
+                              e1, e2, e3,
+                              v3(HUGE_VAL, w, HUGE_VAL)
+                             );
+      geometric_object_list g={ 1, objects };
+      meep_geom::set_materials_from_geometry(&the_structure, g);
+    }
+   else
+    {
+      geometric_object objects[2];
+      objects[0] = make_block(dielectric,
+                              v3(-0.5*pad, wvg_ycen),
+                              e1, e2, e3,
+                              v3(sx-pad, w, HUGE_VAL)
+                             );
+
+      objects[1] = make_block(dielectric,
+                              v3(wvg_xcen, 0.5*pad),
+                              e1, e2, e3,
+                              v3(w, sy-pad, HUGE_VAL)
+                             );
+
+      geometric_object_list g={ 2, objects };
+      meep_geom::set_materials_from_geometry(&the_structure, g);
+    };
+
+  fields f(&the_structure);
+
+  //  (set! sources (list
+  //	       (make source 
+  //		 (src (make gaussian-src (frequency fcen) (fwidth df)))
+  //		 (component Ez)
+  //		 (center (+ 1 (* -0.5 sx)) wvg-ycen)
+  //		 (size 0 w))))
+  //
+  double fcen = 0.15;  // ; pulse center frequency
+  double df   = 0.1;   // ; df
+  gaussian_src_time src(fcen, df);
+  volume v( vec(1.0-0.5*sx, wvg_ycen), vec(0.0,w) );
+  f.add_volume_source(Ez, src, v);
+
+  //(define trans ; transmitted flux
+  //      (add-flux fcen df nfreq
+  //		(if no-bend?
+  //		    (make flux-region
+  //		      (center (- (/ sx 2) 1.5) wvg-ycen) (size 0 (* w 2)))
+  //		    (make flux-region
+  //		      (center wvg-xcen (- (/ sy 2) 1.5)) (size (* w 2) 0)))))
+  double f_start = fcen-0.5*df;
+  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), vec(0.0, 2.0*w))
+           : new volume(vec(wvg_xcen, 0.5*sy-1.5), vec(2.0*w, 0.0));
+  volume_list trans_vl = volume_list(*trans_volume, Sz);
+  dft_flux trans=f.add_dft_flux(&trans_vl, 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), vec(0.0,2.0*w));
+  volume_list refl_vl= volume_list(refl_volume, Sz);
+  dft_flux refl=f.add_dft_flux(&refl_vl, f_start, f_end, nfreq);
+
+  //; for normal run, load negated fields to subtract incident from refl. fields
+  //(if (not no-bend?) (load-minus-flux "refl-flux" refl))
+  const char *dataname="refl-flux";
+  if (!no_bend)
+   { refl.load_hdf5(f, dataname);
+     refl.scale_dfts(-1.0);
+   };
+
+  //(run-sources+ 
+  // (stop-when-fields-decayed 50 Ez
+  //			   (if no-bend?
+  //			       (vector3 (- (/ sx 2) 1.5) wvg-ycen)
+  //			       (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);
+  double DeltaT=50.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();
+
+     // manually check fields-decayed condition
+     double absEz = abs(f.get_field(Ez, eval_point));
+     cur_max = fmax(cur_max, absEz);
+     if ( f.round_time() >= NextCheckTime )
+      { NextCheckTime += DeltaT;
+        max_abs = fmax(max_abs, cur_max);
+        if ( (max_abs>0.0) && cur_max < Tol*max_abs)
+         Done=true;
+        cur_max=0.0;
+      };
+     //printf("%.2e %.2e %.2e %.2e\n",f.round_time(),absEz,max_abs,cur_max);
+   } while(!Done);
+
+  //; for normalization run, save flux fields for refl. plane
+  //(if no-bend? (save-flux "refl-flux" refl))
+  if (no_bend)
+   refl.save_hdf5(f, dataname);
+
+  //(display-fluxes trans refl)
+  printf("%11s | %12s | %12s\n", "   Time    ", " trans flux", "  refl flux");
+  double f0=fcen-0.5*df, fstep=df/(nfreq-1);
+  double *trans_flux=trans.flux();
+  double *refl_flux=refl.flux();
+  for(int nf=0; nf<nfreq; nf++)
+   printf("%.4e | %+.4e | %+.4e\n",f0+nf*fstep,trans_flux[nf],refl_flux[nf]);
+
+}
+
+/***************************************************************/
+/***************************************************************/
+/***************************************************************/
+int main(int argc, char *argv[])
+{
+  initialize mpi(argc, argv);
+
+  bend_flux(true);
+  bend_flux(false);
+
+  // success if we made it here
+  exit(0);
+}

libmeepgeom/material_data.hpp

@@ -0,0 +1,100 @@
+/* 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.
+*/
+#ifndef MATERIAL_DATA_H
+#define MATERIAL_DATA_H
+
+/***************************************************************/
+/* the void *data field in geometric_object_struct points to   */
+/* a material_data structure, defined below.                   */
+/***************************************************************/
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <meep.hpp>
+#include <ctlgeom.h>
+
+#include "meep.hpp"
+
+namespace meep_geom {
+
+typedef struct susceptibility_struct {
+  vector3 sigma_offdiag;
+  vector3 sigma_diag;
+  double frequency;
+  double gamma;
+  double noise_amp;
+  bool drude;
+  bool is_self;
+} susceptibility;
+
+typedef struct {
+  int num_items;
+  susceptibility *items;
+} susceptibility_list;
+
+// prototype for user-supplied material function
+// TODO: this prototype only allows user-defined scalar dielectric permeabilities
+//       (which seems to be all that is possible via the libctl interface as well).
+//       extend to allow more complicated user-specified materials?
+//       
+typedef std::complex<double> (*user_material_func)(vector3 x, void *user_data);
+
+typedef struct medium_struct {
+  vector3 epsilon_diag;
+  vector3 epsilon_offdiag;
+  vector3 mu_diag;
+  vector3 mu_offdiag;
+  susceptibility_list E_susceptibilities;
+  susceptibility_list H_susceptibilities;
+  vector3 E_chi2_diag;
+  vector3 E_chi3_diag;
+  vector3 H_chi2_diag;
+  vector3 H_chi3_diag;
+  vector3 D_conductivity_diag;
+  vector3 B_conductivity_diag;
+} medium_struct;
+
+typedef struct material_type_list
+ { material_type *items;
+   int num_items;
+ } material_type_list;
+
+typedef struct material_data_struct
+ {
+   enum { MATERIAL_TYPE_SELF, MATERIAL_FUNCTION, PERFECT_METAL, MEDIUM } which_subclass;
+
+   // these fields used only if which_subclass==MATERIAL_FUNCTION
+   user_material_func user_func;
+   void *             user_data;
+
+   // used only if which_subclass==MEDIUM
+   medium_struct *medium;
+
+ } material_data;
+
+// global variables and exported functions
+extern material_type vacuum;
+material_type make_dielectric(double epsilon);
+
+void read_epsilon_file(const char *eps_input_file);
+
+
+}; // namespace meep_geom
+
+#endif // #ifndef MATERIAL_DATA_H