Merge remote-tracking branch 'origin/master' into bravais-dev

# Conflicts:
#	fem/fe_coll.cpp
#	fem/geom.hpp
#	fem/intrules.cpp
#	mesh/mesh.cpp
#	mesh/mesh.hpp
#	mesh/pyramid.cpp
#	mesh/pyramid.hpp

miniapps/common/makefile

@@ -41,7 +41,7 @@ ifneq (clean,$(MAKECMDGOALS))
    override shared = $(if $(MFEM_SHARED:YES=),,YES)
 endif
 
-SEQ_MINIOBJS = mesh_extras.o fem_extras.o
+SEQ_MINIOBJS = mesh_extras.o fem_extras.o bravais.o
 ifeq ($(MFEM_USE_MPI),NO)
    MINIOBJS = $(SEQ_MINIOBJS)
 else

miniapps/meta-material/makefile

@@ -0,0 +1,81 @@
+# Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at the
+# Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights reserved.
+# See file COPYRIGHT for details.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability see http://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the GNU Lesser General Public License (as published by the Free
+# Software Foundation) version 2.1 dated February 1999.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/meta-material/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_MINIAPPS = 
+PAR_MINIAPPS = meta_material
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+COMMON_O=../common/bravais.o ../common/pfem_extras.o
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) $(COMMON_O)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ \
+	$(COMMON_O) $(MFEM_LIBS)
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+$(MINIAPPS): \
+%: $(SRC)%.cpp %_solver.o $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $@_solver.o $(COMMON_O) $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+MFEM_TESTS = MINIAPPS
+include $(MFEM_TEST_MK)
+
+# Testing: Parallel vs. serial runs
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) 4
+volta-test-par: meta_material
+	@$(call mfem-test,$<, $(RUN_MPI), Meta-Material miniapp,\
+	-bl 1 -rs 0 -rp 0)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ test_meta_mat meta_material
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -f mobius-strip.mesh klein-bottle.mesh mesh-explorer.mesh
+	@rm -f partitioning.txt shaper.mesh

miniapps/meta-material/meta_material.cpp

@@ -0,0 +1,692 @@
+#include "mfem.hpp"
+#include "../common/bravais.hpp"
+#include <fstream>
+#include <iostream>
+#include <cerrno>      // errno
+
+#include "meta_material_solver.hpp"
+
+#ifndef _WIN32
+#include <sys/stat.h>  // mkdir
+#else
+#include <direct.h>    // _mkdir
+#define mkdir(dir, mode) _mkdir(dir)
+#endif
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::miniapps;
+using namespace mfem::bravais;
+
+// Volume Fraction Coefficient
+static int prob_ = -1;
+double vol_frac_coef(const Vector &);
+
+int CreateDirectory(const string &dir_name, MPI_Comm & comm, int myid);
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Comm comm = MPI_COMM_WORLD;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(comm, &num_procs);
+   MPI_Comm_rank(comm, &myid);
+
+   // 2. Parse command-line options.
+   int bl_type = 1;
+   string lattice_label = "";
+   int order = 1;
+   // int sr = 0, pr = 2;
+   int logging = 0;
+
+   bool visualization = 1;
+   bool visit = true;
+   bool densityCalc = true;
+   bool stiffnessCalc = false;
+   bool bandGapCalc = false;
+   bool band_gap_mid_pts = false;
+
+   double a = -1.0, b = -1.0, c = -1.0;
+   double alpha = -1.0, beta = -1.0, gamma = -1.0;
+   double alpha_deg = -1.0, beta_deg = -1.0, gamma_deg = -1.0;
+   double lcf = 0.3;
+   double density_tol = 0.05;
+   double stiffness_tol = 0.05;
+   double band_gap_tol = 0.05;
+   int    band_gap_max_ref = 2;
+   int    band_gap_samp_pow = 2;
+   // double lambda = 2.07748e+9;
+   // double mu = 0.729927e+9;
+   // Gallium Arsenide at T=300K
+   // double lambda = 5.34e+11;
+   // double mu = 3.285e+11;
+   // double rho0 = 0.0;
+
+   // Acrylonitrile Butadiene Styrene (ABS)
+   // double    rho = 1110.0; // Mass density 1.075 kg/m^3
+   double    rho = 1000.0; // Mass density 1.0 kg/m^3
+   double      E = 2.0e9;  // Young's Modulus 2GPa
+   double     nu = 0.4064; // Poisson's Ratio
+   double epsRel = 10.0;   // Relative Dielectric Permittivity
+   double  muRel = 1.0;    // Relative Magnetic Permeability
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&bl_type, "-bl", "--bravais-lattice",
+                  "Bravais Lattice Type: \n"
+                  "  1 - Primitive Cubic (a),\n"
+                  "  2 - Face-Centered Cubic (a),\n"
+                  "  3 - Body-Centered Cubic (a),\n"
+                  "  4 - Tetragonal (a, c),\n"
+                  "  5 - Body-Centered Tetragonal (a, c),\n"
+                  "  6 - Orthorhombic (a < b < c),\n"
+                  "  7 - Face-Centered Orthorhombic (a < b < c),\n"
+                  "  8 - Body-Centered Orthorhombic (a < b < c),\n"
+                  "  9 - C-Centered Orthorhombic (a < b, c),\n"
+                  " 10 - Hexagonal Prism (a, c),\n"
+                  " 11 - Rhombohedral (a, 0 < alpha < pi),\n"
+                  " 12 - Monoclinic (a, b <= c, 0 < alpha < pi/2),\n"
+                  " 13 - C-Centered Monoclinic (a, b <= c, 0 < alpha < pi/2),\n"
+                  " 14 - Triclinic (0 < alpha, beta, gamma < pi)\n"
+                 );
+   args.AddOption(&a, "-a", "--lattice-a",
+                  "Lattice spacing a");
+   args.AddOption(&b, "-b", "--lattice-b",
+                  "Lattice spacing b");
+   args.AddOption(&c, "-c", "--lattice-c",
+                  "Lattice spacing c");
+   args.AddOption(&alpha, "-alpha", "--lattice-alpha",
+                  "Lattice angle alpha");
+   args.AddOption(&beta, "-beta", "--lattice-beta",
+                  "Lattice angle beta");
+   args.AddOption(&gamma, "-gamma", "--lattice-gamma",
+                  "Lattice angle gamma");
+   args.AddOption(&alpha_deg, "-alpha-deg", "--lattice-alpha-degrees",
+                  "Lattice angle alpha in degrees");
+   args.AddOption(&beta_deg, "-beta-deg", "--lattice-beta-degrees",
+                  "Lattice angle beta in degrees");
+   args.AddOption(&gamma_deg, "-gamma-deg", "--lattice-gamma-degrees",
+                  "Lattice angle gamma in degrees");
+   args.AddOption(&lcf, "-lcf", "--lattice-coef-frac",
+                  "Fraction of inscribed circle radius for rods");
+   args.AddOption(&density_tol, "-rtol", "--density-tolerance",
+                  "Stopping tolerance specified as a relative difference "
+                  "in computed density");
+   args.AddOption(&stiffness_tol, "-ctol", "--stiffness-tolerance",
+                  "Stopping tolerance specified as a relative difference "
+                  "in the 2-norm of the computed stiffness tensor");
+   args.AddOption(&band_gap_tol, "-bgtol", "--band-gap-tolerance",
+                  "Stopping tolerance specified as a relative difference "
+                  "in the computed band gap");
+   args.AddOption(&band_gap_max_ref, "-bgmr", "--band-gap-max-ref",
+                  "Maximum number of uniform mesh refinements to perform "
+                  "when computing the band gap");
+   args.AddOption(&band_gap_samp_pow, "-bgsp", "--band-gap-samp-pow",
+                  "Number of sampling points per segment is 2^samp_pow "
+                  "when computing the band gap");
+   args.AddOption(&band_gap_mid_pts, "-bgmp", "--band-gap-mid-pts",
+                  "-no-bgmp", "--no-band-gap-mid-pts",
+                  "Whether to include mid points in the band gap computation");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   // args.AddOption(&sr, "-sr", "--serial-refinement",
+   //               "Number of serial refinement levels.");
+   // args.AddOption(&pr, "-pr", "--parallel-refinement",
+   //               "Number of parallel refinement levels.");
+   // args.AddOption(&prob_, "-p", "--problem-type",
+   //                "Problem Geometry.");
+   // args.AddOption(&lambda, "-l", "--lambda",
+   //               "Lambda");
+   // args.AddOption(&mu, "-m", "--mu",
+   //               "Mu");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit",
+                  "--no-visualization",
+                  "Enable or disable VisIt visualization.");
+   args.AddOption(&densityCalc, "-rho", "--density", "-no-rho", "--no-density",
+                  "Enable or disable density calculation.");
+   args.AddOption(&stiffnessCalc, "-C", "--stiffness",
+                  "-no-C", "--no-stiffness",
+                  "Enable or disable stiffness tensor calculation.");
+   // args.AddOption(&dispersionPlot, "-disp", "--dispersion",
+   //             "-no-disp", "--no-dispersion",
+   //             "Enable or disable dispersion plot calculation.");
+   args.AddOption(&bandGapCalc, "-bg", "--band-gap",
+                  "-no-bg", "--no-band-gap",
+                  "Enable or disable band gap calculation.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   if ( alpha_deg > 0.0 ) { alpha = alpha_deg * M_PI / 180.0; }
+   if (  beta_deg > 0.0 ) { beta =  beta_deg * M_PI / 180.0; }
+   if ( gamma_deg > 0.0 ) { gamma = gamma_deg * M_PI / 180.0; }
+
+   BRAVAIS_LATTICE_TYPE lattice_type = (BRAVAIS_LATTICE_TYPE)(bl_type + 5);
+   BravaisLattice * bravais = BravaisLatticeFactory(lattice_type,
+                                                    a, b, c,
+                                                    alpha, beta, gamma,
+                                                    logging);
+   BravaisLattice3D * bravais3d = dynamic_cast<BravaisLattice3D*>(bravais);
+
+   lattice_label = bravais->GetLatticeTypeLabel();
+
+   if ( myid == 0 )
+   {
+      cout << "Lattice Type: " << lattice_label << endl;
+   }
+
+   bravais3d->GetAxialLengths(a, b, c);
+   bravais3d->GetInteraxialAngles(alpha, beta, gamma);
+
+   ostringstream oss_prefix;
+   oss_prefix << "Meta-Material-" << lattice_label
+              << "-" << (int)round(100.0*a)
+              << "-" << (int)round(100.0*b)
+              << "-" << (int)round(100.0*c)
+              << "-" << (int)round(180.0*alpha / M_PI)
+              << "-" << (int)round(180.0*beta  / M_PI)
+              << "-" << (int)round(180.0*gamma / M_PI)
+              << "-" << (int)round(1000.0*lcf);
+   // << "-r" << sr + pr;
+
+   CreateDirectory(oss_prefix.str(),comm,myid);
+
+   // 3. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   /*
+   Mesh * mesh = bravais->GetPeriodicWignerSeitzMesh();
+
+   int euler = mesh->EulerNumber();
+   if ( myid == 0 ) { cout << "Initial Euler Number:  " << euler << endl; }
+   mesh->CheckElementOrientation(false);
+   mesh->CheckBdrElementOrientation(false);
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement.
+   {
+      int ref_levels = sr;
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+         int euler = mesh->EulerNumber();
+         if ( myid == 0 )
+         {
+            cout << l+1 << ", Refined Euler Number:  " << euler << endl;
+         }
+         mesh->CheckElementOrientation(false);
+         mesh->CheckBdrElementOrientation(false);
+      }
+   }
+   */
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   // ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   /*
+   delete mesh;
+   {
+      int par_ref_levels = pr;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   */
+   // L2_ParFESpace * L2FESpace    = new L2_ParFESpace(pmesh, 0,
+   //                                               pmesh->Dimension());
+   /*
+   // int nElems = L2FESpace->GetVSize();
+   int nElems = mesh->GetNE();
+   cout << myid << ": nElems = " << nElems << endl;
+   delete mesh;
+   */
+   // LatticeCoefficient latCoef(*bravais, lcf);
+
+   // ParGridFunction * vf0 = new ParGridFunction(L2FESpace);
+   // ParGridFunction * vf1 = new ParGridFunction(L2FESpace);
+
+   // FunctionCoefficient vfFunc(vol_frac_coef);
+   // vf0->ProjectCoefficient(latCoef);
+   /*
+   vf1->ProjectCoefficient(vfFunc);
+
+   double vf13 = (*vf0)[13];
+   double dvf = 0.01 * (0.5 - vf13);
+   (*vf1)[13] += dvf;
+   */
+   VisData vd("localhost", 19916, 1440, 900, 238, 238, 10, 45);
+   /*
+   if (visualization)
+   {
+      socketstream vf_sock;
+      VisualizeField(vf_sock, *vf0, "Volume Fraction 0", vd);
+      vd.IncrementWindow();
+      // socketstream vf1_sock;
+      // VisualizeField(vf1_sock, *vf1, "Volume Fraction 1", vd);
+      // vd.IncrementWindow();
+   }
+   */
+   /*
+   meta_material::Density density(*pmesh, bravais->GetUnitCellVolume(),
+                                  0.0, rho);
+
+   density.SetVolumeFraction(*vf0);
+   */
+   if ( myid == 0 )
+   {
+      cout << "Volume of the unit cell: "
+           << bravais->GetUnitCellVolume() << endl;
+   }
+
+   if (densityCalc)
+   {
+      LatticeCoefficient rhoCoef(*bravais, lcf, 0.0, rho);
+
+      // The density computation does not require a periodic mesh
+      Mesh * mesh_rho = bravais->GetWignerSeitzMesh();
+      mesh_rho->EnsureNCMesh();
+      mesh_rho->UniformRefinement();
+      mesh_rho->UniformRefinement();
+      ParMesh *pmesh_rho = new ParMesh(MPI_COMM_WORLD, *mesh_rho);
+      delete mesh_rho;
+
+      meta_material::Density density(*pmesh_rho, rho,
+                                     bravais->GetUnitCellVolume(),
+                                     rhoCoef, density_tol);
+
+      vector<double> effective_rho;
+      density.GetHomogenizedProperties(effective_rho);
+      if ( myid == 0 )
+      {
+         ostringstream oss;
+         oss << oss_prefix.str() << "/density.dat";
+         ofstream ofs(oss.str().c_str());
+
+         cout << "Effective Density:  ";
+         for (unsigned int i=0; i<effective_rho.size(); i++)
+         {
+            cout << effective_rho[i]; ofs << effective_rho[i];
+            if ( i < effective_rho.size()-1 ) { cout << ", "; ofs << "\t"; }
+         }
+         cout << endl; ofs << endl;
+         ofs.close();
+      }
+
+      if (visualization)
+      {
+         density.InitializeGLVis(vd);
+         density.DisplayToGLVis();
+      }
+      if ( visit )
+      {
+         density.WriteVisItFields(oss_prefix.str(), "Density");
+      }
+      delete pmesh_rho;
+   }
+
+   if (stiffnessCalc)
+   {
+      double lambda = E * nu / ( (1.0 + nu) * (1.0 - 2.0 * nu) );
+      double mu = 0.5 * E / (1.0 + nu);
+      double mat_scale = 1.0e-6;
+
+      LatticeCoefficient lambdaCoef(*bravais, lcf, lambda * mat_scale, lambda);
+      LatticeCoefficient muCoef(*bravais, lcf, mu * mat_scale, mu);
+      /*
+      meta_material::StiffnessTensor elasticity(*pmesh,
+                                                bravais->GetUnitCellVolume(),
+                                                lambda * mat_scale, mu * mat_scale,
+                                                lambda, mu);
+
+      elasticity.SetVolumeFraction(*vf0);
+      */
+      Mesh * mesh_C = bravais->GetPeriodicWignerSeitzMesh();
+
+      if ( mesh_C->EulerNumber() != 0 )
+      {
+         MFEM_ABORT("Euler number equal to " << mesh_C->EulerNumber()
+                    << ". Periodic Bravais Lattice meshes "
+                    "should have Euler number 0!");
+      }
+
+      // mesh_C->UniformRefinement();
+      mesh_C->UniformRefinement();
+      mesh_C->EnsureNCMesh();
+
+      ParMesh *pmesh_C = new ParMesh(MPI_COMM_WORLD, *mesh_C);
+      delete mesh_C;
+
+      meta_material::StiffnessTensor elasticity(*pmesh_C,
+                                                bravais->GetUnitCellVolume(),
+                                                lambdaCoef, muCoef,
+                                                stiffness_tol);
+
+      vector<double> elas;
+      elasticity.GetHomogenizedProperties(elas);
+      if ( myid == 0 )
+      {
+         ostringstream oss;
+         oss << oss_prefix.str() << "/stiffness_tensor.dat";
+         ofstream ofs(oss.str().c_str());
+
+         cout << "Effective Elasticity Tensor:  " << endl;
+         int k = 0;
+         for (unsigned int i=0; i<6; i++)
+         {
+            for (unsigned int j=0; j<i; j++)
+            {
+               cout << " -----------";
+               ofs << elas[(11 - j) * j / 2 + i] << "\t";
+            }
+            for (unsigned int j=i; j<6; j++)
+            {
+               cout << " " << elas[k]; ofs << elas[k];
+               if ( k < 20 ) { ofs << "\t"; }
+               k++;
+            }
+            cout << endl; ofs << endl;
+         }
+         cout << endl;
+         ofs.close();
+      }
+      if (visualization)
+      {
+         elasticity.InitializeGLVis(vd);
+         elasticity.DisplayToGLVis();
+      }
+      if ( visit )
+      {
+         elasticity.WriteVisItFields(oss_prefix.str(), "StiffnessTensor");
+      }
+      delete pmesh_C;
+   }
+
+   if (bandGapCalc)
+   {
+      LatticeCoefficient epsCoef(*bravais, lcf, 1.0, epsRel);
+      LatticeCoefficient  muCoef(*bravais, lcf, 1.0,  muRel);
+
+      Mesh * mesh_bg = bravais->GetPeriodicWignerSeitzMesh();
+
+      if ( mesh_bg->EulerNumber() != 0 )
+      {
+         MFEM_ABORT("Euler number equal to " << mesh_bg->EulerNumber()
+                    << ". Periodic Bravais Lattice meshes "
+                    "should have Euler number 0!");
+      }
+      // mesh_bg->UniformRefinement();
+
+      ParMesh *pmesh_bg = new ParMesh(MPI_COMM_WORLD, *mesh_bg);
+      delete mesh_bg;
+
+      meta_material::MaxwellBandGap maxwell_bg(*pmesh_bg, *bravais,
+                                               band_gap_samp_pow,
+                                               epsCoef, muCoef,
+                                               band_gap_mid_pts,
+                                               band_gap_max_ref,
+                                               band_gap_tol);
+
+      vector<double> bg;
+      maxwell_bg.GetHomogenizedProperties(bg);
+
+      if ( myid == 0 )
+      {
+         ostringstream oss;
+         oss << oss_prefix.str() << "/disp.dat";
+         ofstream ofs(oss.str().c_str());
+
+         maxwell_bg.PrintDispersionPlot(ofs);
+         ofs.close();
+      }
+      if (visualization)
+      {
+         maxwell_bg.InitializeGLVis(vd);
+         maxwell_bg.DisplayToGLVis();
+      }
+      if ( visit )
+      {
+         maxwell_bg.WriteVisItFields(oss_prefix.str(), "MaxwellBandGap");
+      }
+      delete pmesh_bg;
+   }
+
+   // delete vf0;
+   // delete vf1;
+   // delete L2FESpace;
+   // delete pmesh;
+   delete bravais;
+
+   MPI_Finalize();
+
+   if ( myid == 0 )
+   {
+      cout << "Exiting Main" << endl;
+   }
+
+   return 0;
+}
+
+int CreateDirectory(const string &dir_name, MPI_Comm & comm, int myid)
+{
+   int err;
+#ifndef MFEM_USE_MPI
+   err = mkdir(dir_name.c_str(), 0775);
+   err = (err && (errno != EEXIST)) ? 1 : 0;
+#else
+   if (myid == 0)
+   {
+      err = mkdir(dir_name.c_str(), 0775);
+      err = (err && (errno != EEXIST)) ? 1 : 0;
+      MPI_Bcast(&err, 1, MPI_INT, 0, comm);
+   }
+   else
+   {
+      // Wait for rank 0 to create the directory
+      MPI_Bcast(&err, 1, MPI_INT, 0, comm);
+   }
+#endif
+   return err;
+}
+
+double
+distToLine(double ox, double oy, double oz,
+           double tx, double ty, double tz, const Vector & x)
+{
+   double xo_data[3];
+   double xt_data[3];
+   Vector xo(xo_data, 3);
+   Vector xt(xt_data, 3);
+
+   // xo = x - {ox,oy,oz}
+   xo[0] = x[0] - ox;
+   xo[1] = x[1] - oy;
+   xo[2] = x[2] - oz;
+
+   // xt = cross_product({tx,ty,tz}, xo)
+   xt[0] = ty * xo[2] - tz * xo[1];
+   xt[1] = tz * xo[0] - tx * xo[2];
+   xt[2] = tx * xo[1] - ty * xo[0];
+
+   return xt.Norml2();
+}
+double
+vol_frac_coef(const Vector & x)
+{
+   switch ( prob_ )
+   {
+      case -1:
+         // Uniform
+         return 1.0;
+         break;
+      case 0:
+         // Slab
+         if ( fabs(x(0)) <= 0.25 ) { return 1.0; }
+         break;
+      case 1:
+         // Cylinder
+         if ( sqrt(x(0)*x(0)+x(1)*x(1)) <= 0.5 ) { return 1.0; }
+         break;
+      case 2:
+         // Sphere
+         if ( x.Norml2() <= 0.5 ) { return 1.0; }
+         break;
+      case 3:
+         // Sphere and 3 Rods
+      {
+         double r1 = 0.14, r2 = 0.36, r3 = 0.105;
+         if ( x.Norml2() <= r1 ) { return 0.0; }
+         if ( x.Norml2() <= r2 ) { return 1.0; }
+         if ( sqrt(x(1)*x(1)+x(2)*x(2)) <= r3 ) { return 1.0; }
+         if ( sqrt(x(2)*x(2)+x(0)*x(0)) <= r3 ) { return 1.0; }
+         if ( sqrt(x(0)*x(0)+x(1)*x(1)) <= r3 ) { return 1.0; }
+      }
+      break;
+      case 4:
+         // Sphere and 4 Rods
+      {
+         double r1 = 0.14, r2 = 0.28, r3 = 0.1;
+         if ( x.Norml2() <= r1 ) { return 0.0; }
+         if ( x.Norml2() <= r2 ) { return 1.0; }
+
+         Vector y = x;
+         y[0] -= 0.5; y[1] -= 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] -= 0.5; y[1] -= 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] -= 0.5; y[1] += 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] -= 0.5; y[1] += 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] -= 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] -= 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] += 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] += 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+
+         double a = r3;
+         double b = 1.0/sqrt(3.0);
+         if ( distToLine(0.0, 0.0, 0.0, b, b, b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0,-b, b, b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0,-b,-b, b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, b,-b, b, x) <= a ) { return 1.0; }
+      }
+      break;
+      case 5:
+         // Two spheres in a BCC configuration
+         if ( x.Norml2() <= 0.3 )
+         {
+            return 1.0;
+         }
+         else
+         {
+            for (int i=0; i<8; i++)
+            {
+               int i1 = i%2;
+               int i2 = (i/2)%2;
+               int i4 = i/4;
+
+               Vector u = x;
+               u(0) -= i1?-0.5:0.5;
+               u(1) -= i2?-0.5:0.5;
+               u(2) -= i4?-0.5:0.5;
+
+               if ( u.Norml2() <= 0.2 ) { return 1.0; }
+            }
+         }
+         break;
+      case 6:
+         // Sphere and 6 Rods
+      {
+         double r1 = 0.12, r2 = 0.19, r3 = 0.08;
+         if ( x.Norml2() <= r1 ) { return 0.0; }
+         if ( x.Norml2() <= r2 ) { return 1.0; }
+
+         Vector y = x;
+         y[0] -= 0.5; y[1] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] -= 0.5; y[1] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[0] += 0.5; y[1] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+
+         y = x; y[1] -= 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[1] -= 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[1] += 0.5; y[2] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[1] += 0.5; y[2] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+
+         y = x; y[2] -= 0.5; y[0] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[2] -= 0.5; y[0] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[2] += 0.5; y[0] -= 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+         y = x; y[2] += 0.5; y[0] += 0.5;
+         if ( y.Norml2() <= r1 ) { return 0.0; }
+         if ( y.Norml2() <= r2 ) { return 1.0; }
+
+         double a = r3;
+         double b = 1.0/sqrt(2.0);
+         if ( distToLine(0.0, 0.0, 0.0, b, b, 0, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, b,-b, 0, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, b, 0, b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, b, 0,-b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, 0, b, b, x) <= a ) { return 1.0; }
+         if ( distToLine(0.0, 0.0, 0.0, 0, b,-b, x) <= a ) { return 1.0; }
+
+      }
+      break;
+      case 7:
+         if ( fabs(x(1)) + fabs(x(2)) < 0.25 ||
+              fabs(x(0)) + fabs(x(2)) < 0.25 ||
+              fabs(x(0)) + fabs(x(1) - 0.5) < 0.25 )
+         {
+            return 1.0;
+         }
+         break;
+   }
+   return 0.0;
+}

miniapps/meta-material/meta_material_solver.hpp

@@ -0,0 +1,731 @@
+// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
+// reserved. See file COPYRIGHT for details.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability see http://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the GNU Lesser General Public License (as published by the Free
+// Software Foundation) version 2.1 dated February 1999.
+
+#ifndef MFEM_META_MATERIAL_SOLVER
+#define MFEM_META_MATERIAL_SOLVER
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "../common/pfem_extras.hpp"
+#include "../common/bravais.hpp"
+
+extern "C"
+{
+   void dsygv_(int *ITYPE, char *JOBZ, char *UPLO, int *N,
+               double *A, int * LDA, double *B, int *LDB, double *W,
+               double *WORK, int *LWORK, int *INFO);
+}
+
+namespace mfem
+{
+
+using common::H1_ParFESpace;
+using common::ND_ParFESpace;
+using common::RT_ParFESpace;
+using common::L2_ParFESpace;
+using common::ParDiscreteGradOperator;
+using common::ParDiscreteCurlOperator;
+using common::ParDiscreteInterpolationOperator;
+using common::VisData;
+using bravais::BravaisLattice;
+
+namespace meta_material
+{
+
+class LinearCoefficient : public GridFunctionCoefficient
+{
+public:
+   LinearCoefficient(GridFunction * gf, double a0, double a1);
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+
+   double GetSensitivity(ElementTransformation &T,
+                         const IntegrationPoint &ip)
+   { return (c1_ - c0_); }
+
+private:
+   double c0_;
+   double c1_;
+};
+/*
+class PenaltyCoefficient : public GridFunctionCoefficient
+{
+public:
+   PenaltyCoefficient(GridFunction * gf, int penalty, double a0, double a1);
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+
+   double GetSensitivity(ElementTransformation &T,
+                         const IntegrationPoint &ip);
+
+private:
+   int penalty_;
+   double c0_;
+   double c1_;
+};
+*/
+class Homogenization
+{
+public:
+   Homogenization(MPI_Comm comm);
+   virtual ~Homogenization() {}
+
+   void SetVolumeFraction(ParGridFunction & vf) { vf_ = &vf; newVF_ = true; }
+
+   virtual void GetHomogenizedProperties(std::vector<double> & p) = 0;
+   virtual void GetPropertySensitivities(std::vector<ParGridFunction> & dp)
+      = 0;
+
+   void InitializeGLVis(VisData & vd) {}
+
+   void DisplayToGLVis() {}
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label) {}
+
+protected:
+   MPI_Comm comm_;
+   int      myid_;
+   int      numProcs_;
+   bool     newVF_;
+
+   ParGridFunction * vf_;
+};
+
+class Density : public Homogenization
+{
+public:
+   // Density(ParMesh & pmesh, double vol,
+   //       double rho0, double rho1);
+   Density(ParMesh & pmesh, double refDensity, double vol,
+           Coefficient &rhoCoef, double tol = 0.05);
+   ~Density();
+
+   // void SetVolumeFraction(ParGridFunction & vf);
+
+   void GetHomogenizedProperties(std::vector<double> & p);
+   void GetPropertySensitivities(std::vector<ParGridFunction> & dp) {}
+
+   void InitializeGLVis(VisData & vd);
+
+   void DisplayToGLVis();
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label);
+
+private:
+   void updateRho();
+
+   ParMesh         * pmesh_;
+   L2_ParFESpace   * L2FESpace_;
+   ParGridFunction * rho_;
+   ParGridFunction * divGradRho_;
+   ParLinearForm   * cellVol_;
+
+   // LinearCoefficient   rhoCoef_;
+   Coefficient * rhoCoef_;
+   ConstantCoefficient one_;
+
+   double refDensity_;
+   double vol_;
+   double tol_;
+
+   VisData      * vd_;
+   socketstream * sock_;
+   socketstream * sock2_;
+};
+
+class StiffnessTensor : public Homogenization
+{
+public:
+   /*
+   StiffnessTensor(ParMesh & pmesh, double vol,
+                    double lambda0, double mu0,
+                    double lambda1, double mu1);
+   */
+   StiffnessTensor(ParMesh & pmesh, double vol,
+                   Coefficient &lambdaCoef, Coefficient &muCoef,
+                   double tol = 0.05);
+   ~StiffnessTensor();
+
+   // void SetVolumeFraction(ParGridFunction & vf);
+
+   void GetHomogenizedProperties(std::vector<double> & p);
+   void GetPropertySensitivities(std::vector<ParGridFunction> & dp) {}
+
+   void InitializeGLVis(VisData & vd);
+
+   void DisplayToGLVis();
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label);
+
+private:
+
+   void TensorGradient(const Vector & x, Vector & y);
+   void TensorGradientTranspose(const Vector & x, Vector & y);
+   void TensorMassMatrix(const Vector & x, Vector & y);
+   void RestrictedTensorMassMatrix(int r, const Vector & x, Vector & y);
+   void RestrictedVectorAdd(int r, const Vector & x, Vector & y);
+
+   // Produces a 2-tensor coefficient for the diagonal portion of an
+   // elasticity operator.  The coefficient is defined by:
+   //     / lambda + 2 * mu, i == j == axis
+   // K = |              mu, i == j != axis
+   //     \               0, i != j
+   class DiagElasticityCoef : public MatrixCoefficient
+   {
+   public:
+      DiagElasticityCoef(Coefficient & lambda,
+                         Coefficient & mu,
+                         int axis);
+
+      virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                        const IntegrationPoint &ip);
+
+   private:
+      int axis_;
+      Coefficient * lambda_;
+      Coefficient * mu_;
+   };
+
+   // Produces a 2-tensor coefficient for the off-diagonal portion of an
+   // elasticity operator.  The coefficient is defined by:
+   //     / lambda, i == axis0, j == axis1
+   // K = |     mu, i == axis1, j == axis0
+   //     \      0, otherwise
+   class OffDiagElasticityCoef : public MatrixCoefficient
+   {
+   public:
+      OffDiagElasticityCoef(Coefficient & lambda,
+                            Coefficient & mu,
+                            int axis0, int axis1);
+
+      virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                        const IntegrationPoint &ip);
+
+   private:
+      int axis0_;
+      int axis1_;
+      Coefficient * lambda_;
+      Coefficient * mu_;
+   };
+
+   void solve(const Vector & E, Vector & Chi);
+
+   int dim_;
+   int irOrder_;
+   int geom_;
+   bool amg_elast_;
+
+   ParMesh * pmesh_;
+
+   L2_ParFESpace * L2FESpace_;
+   L2_ParFESpace * L2VFESpace_;
+   H1_ParFESpace * H1FESpace_;
+   H1_ParFESpace * H1VFESpace_;
+   ND_ParFESpace * HCurlFESpace_;
+   ND_ParFESpace * HCurlVFESpace_;
+   RT_ParFESpace * HDivFESpace_;
+
+   // LinearCoefficient lambdaCoef_;
+   // LinearCoefficient muCoef_;
+   Coefficient * lambdaCoef_;
+   Coefficient * muCoef_;
+
+   DiagElasticityCoef    xxCoef_;
+   DiagElasticityCoef    yyCoef_;
+   DiagElasticityCoef    zzCoef_;
+   OffDiagElasticityCoef yzCoef_;
+   OffDiagElasticityCoef xzCoef_;
+   OffDiagElasticityCoef xyCoef_;
+
+   ParGridFunction * lambda_;
+   ParGridFunction * mu_;
+
+   ParBilinearForm * a_;
+   ParBilinearForm * m_[6];
+   ParGridFunction * Chi_[6];
+   ParGridFunction * E_[3];
+   ParGridFunction * F_[6];
+   ParLinearForm   * MF_[6];
+
+   ParDiscreteGradOperator * grad_;
+
+   BilinearFormIntegrator * diffInteg_[3];
+   ParGridFunction * errSol_;
+   ErrorEstimator  * errEst_[3];
+   ParGridFunction * errors_[4];
+
+   ParGridFunction * b_;
+   // HypreParVector  * tmp1_;
+
+   double vol_;
+   double tol_;
+
+   VisData      * vd_;
+   socketstream   socks_[8];
+   socketstream   err_socks_[4];
+
+   int seqVF_;
+};
+
+class ParDiscreteVectorProductOperator
+   : public ParDiscreteInterpolationOperator
+{
+public:
+   ParDiscreteVectorProductOperator(ParFiniteElementSpace *dfes,
+                                    ParFiniteElementSpace *rfes,
+                                    const Vector & v);
+
+private:
+   VectorConstantCoefficient vCoef_;
+};
+
+class ParDiscreteVectorCrossProductOperator
+   : public ParDiscreteInterpolationOperator
+{
+public:
+   ParDiscreteVectorCrossProductOperator(ParFiniteElementSpace *dfes,
+                                         ParFiniteElementSpace *rfes,
+                                         const Vector & v);
+private:
+   VectorConstantCoefficient vCoef_;
+};
+
+//class MaxwellBlochWaveEquation;
+//class MaxwellBlochWaveEquation::MaxwellBlochWaveProjector;
+
+class MaxwellBlochWaveEquation
+{
+public:
+   MaxwellBlochWaveEquation(ParMesh & pmesh,
+                            int order);
+   ~MaxwellBlochWaveEquation();
+
+   // Where kappa is the phase shift vector
+   void SetKappa(const Vector & kappa);
+
+   // Where beta*zeta = kappa
+   void SetBeta(double beta);
+   void SetZeta(const Vector & zeta);
+
+   // void SetAzimuth(double alpha_a);
+   // void SetInclination(double alpha_i);
+
+   // void SetOmega(double omega);
+   void SetAbsoluteTolerance(double atol);
+   void SetNumEigs(int nev);
+   void SetMassCoef(Coefficient & m);
+   void SetStiffnessCoef(Coefficient & k);
+
+   void Setup();
+
+   void SetInitialVectors(int num_vecs, HypreParVector ** vecs);
+
+   // void SetBravaisLattice(BravaisLattice & bravais) { bravais_ = &bravais; }
+
+   void Update();
+
+   /// Solve the eigenproblem
+   void Solve();
+
+   /// Collect the converged eigenvalues
+   void GetEigenvalues(std::vector<double> & eigenvalues);
+
+   /// A convenience method which combines six methods into one
+   void GetEigenvalues(int nev, const Vector & kappa,
+                       std::vector<HypreParVector*> & init_vecs,
+                       std::vector<double> & eigenvalues);
+
+   /// Extract a single eigenvector
+   HypreParVector * ReturnEigenvector(unsigned int i);
+
+   void CopyEigenvector(unsigned int i,
+                        HypreParVector & V);
+
+   void GetEigenvector(unsigned int i,
+                       HypreParVector & Er,
+                       HypreParVector & Ei,
+                       HypreParVector & Br,
+                       HypreParVector & Bi);
+   void GetEigenvectorE(unsigned int i,
+                        HypreParVector & Er,
+                        HypreParVector & Ei);
+   void GetEigenvectorB(unsigned int i,
+                        HypreParVector & Br,
+                        HypreParVector & Bi);
+
+   BlockOperator * GetAOperator() { return A_; }
+   BlockOperator * GetMOperator() { return M_; }
+
+   Solver   * GetPreconditioner() { return Precond_; }
+   Operator * GetSubSpaceProjector() { return SubSpaceProj_; }
+
+   ParFiniteElementSpace * GetH1FESpace() { return H1FESpace_; }
+   ParFiniteElementSpace * GetHCurlFESpace() { return HCurlFESpace_; }
+   ParFiniteElementSpace * GetHDivFESpace()  { return HDivFESpace_; }
+
+   // void TestVector(const HypreParVector & v);
+
+   ParGridFunction * GetEigenvectorEnergy(unsigned int i) { return energy_[i]; }
+
+   void GetFourierCoefficients(HypreParVector & Vr,
+                               HypreParVector & Vi,
+                               Array2D<double> &f);
+
+   void IdentifyDegeneracies(double zero_tol, double rel_tol,
+                             std::vector<std::set<int> > & degen);
+
+   void GetFieldAverages(unsigned int i,
+                         Vector & Er, Vector & Ei,
+                         Vector & Br, Vector & Bi,
+                         Vector & Dr, Vector & Di,
+                         Vector & Hr, Vector & Hi);
+
+   void ComputeHomogenizedCoefs();
+
+   void DetermineBasis(const Vector & v1, std::vector<Vector> & e);
+
+   void WriteVisitFields(const std::string & prefix,
+                         const std::string & label);
+
+   void GetSolverStats(double &meanTime, double &stdDevTime,
+                       double &meanIter, double &stdDevIter,
+                       int &nSolves);
+
+   void TestProjector() const;
+
+private:
+
+   MPI_Comm comm_;
+   int myid_;
+   int hcurl_loc_size_;
+   int hdiv_loc_size_;
+   int nev_;
+
+   bool newBeta_;
+   bool newZeta_;
+   bool newOmega_;
+   bool newMCoef_;
+   bool newKCoef_;
+
+   ParMesh        * pmesh_;
+   H1_ParFESpace  * H1FESpace_;
+   ND_ParFESpace  * HCurlFESpace_;
+   RT_ParFESpace  * HDivFESpace_;
+   // L2_ParFESpace  * L2FESpace_;
+
+   // BravaisLattice     * bravais_;
+   // HCurlFourierSeries * fourierHCurl_;
+
+   double           atol_;
+   double           beta_;
+
+   Vector           zeta_;
+   Vector           kappa_;
+
+   Coefficient    * mCoef_;
+   Coefficient    * kCoef_;
+
+   Array<int>       block_offsets_;
+   Array<int>       block_trueOffsets_;
+   Array<int>       block_trueOffsets2_;
+   Array<HYPRE_Int> tdof_offsets_;
+
+   BlockOperator  * A_;
+   BlockOperator  * M_;
+   BlockOperator  * C_;
+
+   BlockVector    * blkHCurl_;
+   BlockVector    * blkHDiv_;
+
+   HypreParMatrix * M1_;
+   HypreParMatrix * M2_;
+   HypreParMatrix * S1_;
+   HypreParMatrix * T1_;
+   HypreParMatrix * T12_;
+   HypreParMatrix * Z12_;
+
+   HypreParMatrix * DKZ_;
+   // HypreParMatrix * DKZT_;
+
+   HypreAMS       * T1Inv_;
+
+   ParDiscreteCurlOperator * Curl_;
+   ParDiscreteVectorCrossProductOperator * Zeta_;
+
+   BlockDiagonalPreconditioner * BDP_;
+
+   Solver   * Precond_;
+   //MaxwellBlochWaveProjector * SubSpaceProj_;
+   Operator * SubSpaceProj_;
+
+   HypreParVector ** vecs_;
+   HypreParVector * vec0_;
+
+   HypreLOBPCG * lobpcg_;
+   HypreAME    * ame_;
+
+   ParGridFunction ** energy_;
+   /*
+    HypreParVector * AvgHCurl_coskx_[3];
+    HypreParVector * AvgHCurl_sinkx_[3];
+    HypreParVector * AvgHDiv_coskx_[3];
+    HypreParVector * AvgHDiv_sinkx_[3];
+
+    HypreParVector * AvgHCurl_eps_coskx_[3];
+    HypreParVector * AvgHCurl_eps_sinkx_[3];
+    HypreParVector * AvgHDiv_muInv_coskx_[3];
+    HypreParVector * AvgHDiv_muInv_sinkx_[3];
+   */
+   std::vector<double> solve_times_;
+   std::vector<int>    solve_iters_;
+
+   class MaxwellBlochWavePrecond : public Solver
+   {
+   public:
+      MaxwellBlochWavePrecond(ParFiniteElementSpace & HCurlFESpace,
+                              BlockDiagonalPreconditioner & BDP,
+                              Operator & subSpaceProj,
+                              double w);
+
+      ~MaxwellBlochWavePrecond();
+
+      void Mult(const Vector & x, Vector & y) const;
+
+      void SetOperator(const Operator & A);
+
+   private:
+      int myid_;
+
+      // ParFiniteElementSpace * HCurlFESpace_;
+      BlockDiagonalPreconditioner * BDP_;
+      const Operator * A_;
+      Operator * subSpaceProj_;
+      // mutable HypreParVector *r_, *u_, *v_;
+      mutable HypreParVector *u_;
+      // double w_;
+   };
+
+   class MaxwellBlochWaveProjector : public Operator
+   {
+   public:
+      MaxwellBlochWaveProjector(ParFiniteElementSpace & HCurlFESpace,
+                                ParFiniteElementSpace & H1FESpace,
+                                BlockOperator & M,
+                                double beta, const Vector & zeta);
+      ~MaxwellBlochWaveProjector();
+
+      void SetBeta(double beta);
+      void SetZeta(const Vector & zeta);
+
+      void Setup();
+      virtual void Mult(const Vector &x, Vector &y) const;
+
+   private:
+      int myid_;
+      int locSize_;
+
+      bool newBeta_;
+      bool newZeta_;
+
+      ParFiniteElementSpace * HCurlFESpace_;
+      ParFiniteElementSpace * H1FESpace_;
+
+      double beta_;
+      Vector zeta_;
+
+      HypreParMatrix * T01_;
+      HypreParMatrix * Z01_;
+      HypreParMatrix * A0_;
+      HypreParMatrix * DKZ_;
+
+      MINRESSolver   * minres_;
+
+      Array<int>       block_offsets0_;
+      Array<int>       block_offsets1_;
+      Array<int>       block_trueOffsets0_;
+      Array<int>       block_trueOffsets1_;
+
+      BlockOperator * S0_;
+      BlockOperator * M_;
+      BlockOperator * G_;
+
+      mutable HypreParVector * urDummy_;
+      mutable HypreParVector * uiDummy_;
+      mutable HypreParVector * vrDummy_;
+      mutable HypreParVector * viDummy_;
+
+      mutable BlockVector * u0_;
+      mutable BlockVector * v0_;
+      mutable BlockVector * u1_;
+      mutable BlockVector * v1_;
+   };
+};
+
+class InverseCoefficient : public TransformedCoefficient
+{
+public:
+   InverseCoefficient(Coefficient * q) : TransformedCoefficient(q, inv_) {}
+private:
+   static double inv_(double v) { return 1.0/v; }
+};
+
+class MaxwellBlochWaveSolver
+{
+public:
+   MaxwellBlochWaveSolver(ParMesh & pmesh, BravaisLattice & bravais,
+                          Coefficient & epsCoef, Coefficient & muCoef,
+                          int max_ref = 2, int nev = 24, double tol = 0.05);
+   ~MaxwellBlochWaveSolver();
+
+   // Where kappa is the phase shift vector
+   void SetKappa(const Vector & kappa);
+
+   // Where beta*zeta = kappa and |zeta| = 1
+   void SetBeta(double beta);
+   void SetZeta(const Vector & zeta);
+
+   void GetEigenfrequencies(std::vector<double> & omega);
+
+   MaxwellBlochWaveEquation * GetFineSolver()
+   { return mbwe_[mbwe_.size()-1]; }
+
+   HypreParVector * ReturnFineEigenvector(int i);
+
+   void InitializeGLVis(VisData & vd);
+
+   void DisplayToGLVis();
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label);
+
+private:
+
+   void createPartitioning(ParFiniteElementSpace & pfes, HYPRE_Int *& part);
+
+   int max_lvl_;
+   int nev_;
+   double tol_;
+
+   std::vector<ParMesh*> pmesh_;
+   std::vector<MaxwellBlochWaveEquation*> mbwe_;
+   std::vector<const Operator*> refineOp_;
+   std::vector<std::pair<HypreParVector*, HypreParVector*> > EField_;
+   //std::vector<std::vector<std::pair<ParGridFunction,
+   //                ParGridFunction> > > efield_;
+   std::vector<std::pair<ParGridFunction*, ParGridFunction*> > efield_;
+   std::vector<std::vector<HypreParVector*> > initialVecs_;
+   std::vector<int> locSize_;
+   std::vector<HYPRE_Int*> part_;
+
+   Vector kappa_;
+   Coefficient * epsCoef_;
+   InverseCoefficient muInvCoef_;
+
+   // Coefficient * muCoef_;
+
+};
+
+class MaxwellDispersion
+{
+public:
+   MaxwellDispersion(ParMesh & pmesh, BravaisLattice & bravais,
+                     int sample_power,
+                     Coefficient & epsCoef, Coefficient & muCoef,
+                     bool midPts = false, int max_ref = 2,
+                     int nev = 24, double tol = 0.05);
+   ~MaxwellDispersion();
+
+   const std::vector<std::vector<std::map<int,std::vector<double> > > > &
+   GetDispersionData();
+
+   void PrintDispersionPlot(std::ostream & os);
+
+   void InitializeGLVis(VisData & vd);
+
+   void DisplayToGLVis();
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label);
+
+private:
+
+   void buildRawBasis();
+
+   void approxEigenfrequencies(std::vector<double> & omega);
+
+   void traverseBrillouinZone();
+
+   std::string modLabel(const std::string & label) const;
+
+   void findAndReplace(const std::string & f, const std::string & r,
+                       std::string & str) const;
+
+   BravaisLattice         * bravais_;
+   MaxwellBlochWaveSolver * mbws_;
+
+   HypreParVector * Ax_;
+   HypreParVector * Mx_;
+
+   DenseMatrix A_;
+   DenseMatrix M_;
+
+   std::vector<HypreParVector*> rawBasis_;
+   std::vector<HypreParVector*> projBasis_;
+
+   std::map<std::string,std::vector<double> > sp_eigs_;
+   std::vector<std::vector<std::map<int,std::vector<double> > > > seg_eigs_;
+
+   int n_pow_;
+   int n_div_;
+   int samp_pow_;
+   int nev_;
+
+   bool midPts_;
+};
+
+class MaxwellBandGap : public Homogenization
+{
+public:
+   MaxwellBandGap(ParMesh & pmesh, BravaisLattice & bravais,
+                  int samp_pow,
+                  Coefficient & epsCoef, Coefficient & muCoef,
+                  bool midPts = false, int max_ref = 2, double tol = 0.05);
+   ~MaxwellBandGap();
+
+   void GetHomogenizedProperties(std::vector<double> & p);
+   void GetPropertySensitivities(std::vector<ParGridFunction> & dp) {}
+
+   void PrintDispersionPlot(std::ostream & os);
+
+   void InitializeGLVis(VisData & vd);
+
+   void DisplayToGLVis();
+
+   void WriteVisItFields(const std::string & prefix,
+                         const std::string & label);
+private:
+
+   MaxwellDispersion * disp_;
+};
+
+} // namespace meta_material
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_META_MATERIAL_SOLVER