Removing experimental example codes

# Conflicts:
#	examples/ex21p.cpp
#	fem/complex_fem.cpp
#	miniapps/electromagnetics/makefile

fem/coefficient.cpp

@@ -285,6 +285,22 @@ void VectorRestrictedCoefficient::Eval(
    }
 }
 
+void UnitNormalCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                 const IntegrationPoint &ip)
+{
+   V.SetSize(vdim);
+   V = 0.0;
+
+   const DenseMatrix & J = T.Jacobian();
+   if (J.Width() == J.Height() - 1)
+   {
+      CalcOrtho(J, V);
+      double norm = V.Norml2();
+      MFEM_ASSERT(norm > 0.0, "Length of normal vector is non-positive!");
+      V /= norm;
+   }
+}
+
 void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {

fem/coefficient.hpp

@@ -685,6 +685,16 @@ public:
                      const IntegrationRule &ir);
 };
 
+/// VectorCoefficient which computes unit normal vector on the mesh boundary
+class UnitNormalCoefficient : public VectorCoefficient
+{
+public:
+   UnitNormalCoefficient(int dim) : VectorCoefficient(dim) {}
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
 
 /// Base class for Matrix Coefficients that optionally depend on time and space.
 class MatrixCoefficient

miniapps/common/fem_extras.cpp

@@ -56,9 +56,45 @@ RT_FESpace::~RT_FESpace()
    delete FEC_;
 }
 
+L2_FESpace::L2_FESpace(Mesh *m, const int p, const int space_dim,
+                       int vdim, int order)
+   : FiniteElementSpace(m, new L2_FECollection(p,space_dim),vdim,order)
+{
+   FEC_ = this->FiniteElementSpace::fec;
+}
+
+L2_FESpace::~L2_FESpace()
+{
+   delete FEC_;
+}
+
+DiscreteInterpolationOperator::~DiscreteInterpolationOperator()
+{}
+
+DiscreteGradOperator::DiscreteGradOperator(FiniteElementSpace *dfes,
+                                           FiniteElementSpace *rfes)
+   : DiscreteInterpolationOperator(dfes, rfes)
+{
+   this->AddDomainInterpolator(new GradientInterpolator);
+}
+
+DiscreteCurlOperator::DiscreteCurlOperator(FiniteElementSpace *dfes,
+                                           FiniteElementSpace *rfes)
+   : DiscreteInterpolationOperator(dfes, rfes)
+{
+   this->AddDomainInterpolator(new CurlInterpolator);
+}
+
+DiscreteDivOperator::DiscreteDivOperator(FiniteElementSpace *dfes,
+                                         FiniteElementSpace *rfes)
+   : DiscreteInterpolationOperator(dfes, rfes)
+{
+   this->AddDomainInterpolator(new DivergenceInterpolator);
+}
+
 void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
                    Mesh &mesh, const char *title,
-                   int x, int y, int w, int h, const char * keys, bool vec)
+                   int x, int y, int w, int h, const char *keys, bool vec)
 {
    bool newly_opened = false;
    int connection_failed;
@@ -93,7 +129,7 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
 
 void VisualizeField(socketstream &sock, const char *vishost, int visport,
                     GridFunction &gf, const char *title,
-                    int x, int y, int w, int h, const char * keys, bool vec)
+                    int x, int y, int w, int h, const char *keys, bool vec)
 {
    Mesh &mesh = *gf.FESpace()->GetMesh();
 

miniapps/common/fem_extras.hpp

@@ -66,6 +66,50 @@ private:
 };
 
 
+/** The L2_FESpace class is a FiniteElementSpace which automatically
+    allocates and destroys its own FiniteElementCollection, in this
+    case an L2_FECollection object.
+*/
+class L2_FESpace : public FiniteElementSpace
+{
+public:
+   L2_FESpace(Mesh *m, const int p, const int space_dim,
+              int vdim = 1, int order = Ordering::byNODES);
+   ~L2_FESpace();
+private:
+   const FiniteElementCollection *FEC_;
+};
+
+class DiscreteInterpolationOperator : public DiscreteLinearOperator
+{
+public:
+   DiscreteInterpolationOperator(FiniteElementSpace *dfes,
+                                 FiniteElementSpace *rfes)
+      : DiscreteLinearOperator(dfes, rfes) {}
+   virtual ~DiscreteInterpolationOperator();
+};
+
+class DiscreteGradOperator : public DiscreteInterpolationOperator
+{
+public:
+   DiscreteGradOperator(FiniteElementSpace *dfes,
+                        FiniteElementSpace *rfes);
+};
+
+class DiscreteCurlOperator : public DiscreteInterpolationOperator
+{
+public:
+   DiscreteCurlOperator(FiniteElementSpace *dfes,
+                        FiniteElementSpace *rfes);
+};
+
+class DiscreteDivOperator : public DiscreteInterpolationOperator
+{
+public:
+   DiscreteDivOperator(FiniteElementSpace *dfes,
+                       FiniteElementSpace *rfes);
+};
+
 /// Visualize the given mesh object, using a GLVis server on the
 /// specified host and port. Set the visualization window title, and optionally,
 /// its geometry.

miniapps/electromagnetics/CMakeLists.txt

@@ -9,6 +9,12 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
+add_mfem_miniapp(hertz_ser
+  MAIN hertz_ser.cpp
+  EXTRA_SOURCES hertz_ser_solver.cpp
+  EXTRA_HEADERS hertz_ser_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+  LIBRARIES mfem mfem_miniapps_common)
+
 if (MFEM_USE_MPI)
   add_mfem_miniapp(tesla
     MAIN tesla.cpp
@@ -22,6 +28,12 @@ if (MFEM_USE_MPI)
     volta_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
     LIBRARIES mfem mfem-common)
 
+ add_mfem_miniapp(hertz
+   MAIN hertz.cpp
+   EXTRA_SOURCES hertz_solver.cpp
+   EXTRA_HEADERS hertz_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+   LIBRARIES mfem mfem_miniapps_common)
+
  add_mfem_miniapp(joule
    MAIN joule.cpp
    joule_solver.cpp

miniapps/electromagnetics/hertz.cpp

@@ -0,0 +1,672 @@
+// 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.
+//
+//   -----------------------------------------------------------------------
+//   Hertz Miniapp:  Simple Frequency-Domain Electromagnetic Simulation Code
+//   -----------------------------------------------------------------------
+//
+//   Assumes that all sources and boundary conditions oscillate with the same
+//   frequency although not necessarily in phase with one another.  This
+//   assumptions implies that we can factor out the time dependence which we
+//   take to be of the form exp(i omega t).  With these assumptions we can
+//   write the Maxwell equations in the form:
+//
+//   i omega epsilon E = Curl mu^{-1} B - J - sigma E
+//   i omega B         = - Curl E
+//
+//   Which combine to yield:
+//
+//   Curl mu^{-1} Curl E - omega^2 epsilon E + i omega sigma E = - i omega J
+//
+//   We discretize this equation with H(Curl) a.k.a Nedelec basis
+//   functions.  The curl curl operator must be handled with
+//   integration by parts which yields a surface integral:
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               + (W, n x (mu^{-1} Curl E))_{\Gamma}
+//
+//   or
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               - i omega (W, n x H)_{\Gamma}
+//
+//   For plane waves
+//     omega B = - k x E
+//     omega D = k x H, assuming n x k = 0 => n x H = omega epsilon E / |k|
+//
+//   c = omega/|k|
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               - i omega sqrt{epsilon/mu} (W, E)_{\Gamma}
+//
+//
+// Compile with: make hertz
+//
+// Sample runs:
+//
+//   By default the sources and fields are all zero
+//     mpirun -np 4 hertz
+//
+//   Current source in a metal sphere
+//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
+//                        -dbcs '-1' -f 3e8 -herm
+//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//
+//   Current source in a sphere with absorbing boundary conditions
+//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
+//                        -abcs '-1' -f 3e8
+//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//
+//   Current source in a metal sphere with dielectric and conducting materials
+//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
+//                        -dbcs '-1' -f 3e8
+//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//                        -cs '0.0 0.0 -0.5 .2 10'
+//                        -ds '0.0 0.0 0.5 .2 10'
+//
+//   Current source in a metal box
+//     mpirun -np 4 hertz -m ../../data/fichera.mesh -rs 3
+//                        -dbcs '-1' -f 3e8
+//                        -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//
+//   Current source with a mixture of absorbing and reflecting boundaries
+//     mpirun -np 4 hertz -m ../../data/fichera.mesh -rs 3
+//                        -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//                        -dbcs '4 8 19 21' -abcs '5 18' -f 3e8
+//
+
+#include "hertz_solver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::electromagnetics;
+
+// Permittivity Functions
+Coefficient * SetupPermittivityCoefficient();
+
+static Vector pw_eps_(0);     // Piecewise permittivity values
+static Vector ds_params_(0);  // Center, Radius, and Permittivity
+//                               of dielectric sphere
+double dielectric_sphere(const Vector &);
+
+// Permeability Function
+Coefficient * SetupInvPermeabilityCoefficient();
+
+static Vector pw_mu_(0);      // Piecewise permeability values
+static Vector pw_mu_inv_(0);  // Piecewise inverse permeability values
+static Vector ms_params_(0);  // Center, Inner and Outer Radii, and
+//                               Permeability of magnetic shell
+double magnetic_shell(const Vector &);
+double magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
+
+// Conductivity Functions
+Coefficient * SetupConductivityCoefficient();
+
+static Vector pw_sigma_(0);   // Piecewise conductivity values
+static Vector cs_params_(0);  // Center, Radius, and Conductivity
+//                               of conductive sphere
+double conductive_sphere(const Vector &);
+
+// Impedance
+Coefficient * SetupAdmittanceCoefficient(const Mesh & mesh,
+                                         const Array<int> & abcs);
+
+static Vector pw_eta_(0);      // Piecewise impedance values
+static Vector pw_eta_inv_(0);  // Piecewise inverse impedance values
+
+// Current Density Function
+static Vector do_params_(0);  // Axis Start, Axis End, Rod Radius,
+//                               Total Current of Rod
+void dipole_oscillator(const Vector &x, Vector &j);
+void j_src(const Vector &x, Vector &j) { dipole_oscillator(x, j); }
+
+// Electric Field Boundary Condition: The following function returns zero but
+// any function could be used.
+void e_bc_r(const Vector &x, Vector &E);
+void e_bc_i(const Vector &x, Vector &E);
+
+static double freq_ = 1.0;
+
+// Prints the program's logo to the given output stream
+void display_banner(ostream & os);
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   if ( mpi.Root() ) { display_banner(cout); }
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/ball-nurbs.mesh";
+   int order = 1;
+   int maxit = 100;
+   int serial_ref_levels = 0;
+   int parallel_ref_levels = 0;
+   int sol = 3;
+   int prec = 4;
+   bool herm_conv = false;
+   bool visualization = true;
+   bool visit = true;
+
+   Array<int> abcs;
+   Array<int> dbcs;
+
+   SolverOptions solOpts;
+   solOpts.maxIter = 1000;
+   solOpts.kDim = 50;
+   solOpts.printLvl = 1;
+   solOpts.relTol = 1e-4;
+   solOpts.euLvl = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
+                  "Number of serial refinement levels.");
+   args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
+                  "Number of parallel refinement levels.");
+   args.AddOption(&freq_, "-f", "--frequency",
+                  "Frequency in Hertz (of course...)");
+   args.AddOption(&prec, "-pc", "--precond",
+                  "Preconditioner: 1 - Diagonal Scaling, 2 - ParaSails, "
+                  "3 - Euclid, 4 - AMS");
+   args.AddOption(&sol, "-s", "--solver",
+                  "Solver: 1 - GMRES, 2 - FGMRES, 3 - MINRES"
+#ifdef MFEM_USE_SUPERLU
+                  ", 4 - SuperLU"
+#endif
+#ifdef MFEM_USE_STRUMPACK
+                  ", 5 - STRUMPACK"
+#endif
+                 );
+   args.AddOption(&solOpts.maxIter, "-sol-it", "--solver-iterations",
+                  "Maximum number of solver iterations.");
+   args.AddOption(&solOpts.kDim, "-sol-k-dim", "--solver-krylov-dimension",
+                  "Krylov space dimension for GMRES and FGMRES.");
+   args.AddOption(&solOpts.relTol, "-sol-tol", "--solver-tolerance",
+                  "Relative tolerance for GMRES or FGMRES.");
+   args.AddOption(&solOpts.printLvl, "-sol-prnt-lvl", "--solver-print-level",
+                  "Logging level for solvers.");
+   args.AddOption(&solOpts.euLvl, "-eu-lvl", "--euclid-level",
+                  "Euclid factorization level for ILU(k).");
+   args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
+                  "Piecewise values of Permittivity");
+   args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
+                  "Center, Radius, and Permittivity of Dielectric Sphere");
+   args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
+                  "Piecewise values of Permeability");
+   args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
+                  "Center, Inner Radius, Outer Radius, "
+                  "and Permeability of Magnetic Shell");
+   args.AddOption(&pw_sigma_, "-pws", "--piecewise-sigma",
+                  "Piecewise values of Conductivity");
+   args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
+                  "Center, Radius, and Conductivity of Conductive Sphere");
+   args.AddOption(&pw_eta_, "-pwz", "--piecewise-eta",
+                  "Piecewise values of Impedance (one value per abc surface)");
+   args.AddOption(&do_params_, "-do", "--dipole-oscillator-params",
+                  "Axis End Points, Radius, and Amplitude");
+   args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
+                  "Absorbing Boundary Condition Surfaces");
+   args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
+                  "Dirichlet Boundary Condition Surfaces");
+   args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
+                  "Max number of iterations in the main AMR loop.");
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
+                  "Enable or disable VisIt visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(cout);
+   }
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 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 = new Mesh(mesh_file, 1, 1);
+
+   if (mpi.Root())
+   {
+      cout << "Starting initialization." << endl;
+   }
+
+   // Project a NURBS mesh to a piecewise-quadratic curved mesh
+   if (mesh->NURBSext)
+   {
+      mesh->UniformRefinement();
+      if (serial_ref_levels > 0) { serial_ref_levels--; }
+
+      mesh->SetCurvature(2);
+   }
+
+   // Ensure that quad and hex meshes are treated as non-conforming.
+   mesh->EnsureNCMesh();
+
+   // Refine the serial mesh on all processors to increase the resolution. In
+   // this example we do 'ref_levels' of uniform refinement.
+   for (int l = 0; l < serial_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 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(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // Refine this mesh in parallel to increase the resolution.
+   int par_ref_levels = parallel_ref_levels;
+   for (int l = 0; l < par_ref_levels; l++)
+   {
+      pmesh.UniformRefinement();
+   }
+
+   // Create a coefficient describing the dielectric permittivity
+   Coefficient * epsCoef = SetupPermittivityCoefficient();
+
+   // Create a coefficient describing the magnetic permeability
+   Coefficient * muInvCoef = SetupInvPermeabilityCoefficient();
+
+   // Create a coefficient describing the electrical conductivity
+   Coefficient * sigmaCoef = SetupConductivityCoefficient();
+
+   // Create a coefficient describing the surface admittance
+   Coefficient * etaInvCoef = SetupAdmittanceCoefficient(pmesh, abcs);
+
+   // Create the Magnetostatic solver
+   HertzSolver Hertz(pmesh, order, freq_,
+                     (HertzSolver::SolverType)sol, solOpts,
+                     (HertzSolver::PrecondType)prec,
+                     conv, *epsCoef, *muInvCoef, sigmaCoef, etaInvCoef,
+                     abcs, dbcs,
+                     e_bc_r, e_bc_i,
+                     (do_params_.Size() > 0 ) ? j_src : NULL, NULL
+                    );
+
+   // Initialize GLVis visualization
+   if (visualization)
+   {
+      Hertz.InitializeGLVis();
+   }
+
+   // Initialize VisIt visualization
+   VisItDataCollection visit_dc("Hertz-AMR-Parallel", &pmesh);
+
+   if ( visit )
+   {
+      Hertz.RegisterVisItFields(visit_dc);
+   }
+   if (mpi.Root()) { cout << "Initialization done." << endl; }
+
+   // The main AMR loop. In each iteration we solve the problem on the current
+   // mesh, visualize the solution, estimate the error on all elements, refine
+   // the worst elements and update all objects to work with the new mesh. We
+   // refine until the maximum number of dofs in the Nedelec finite element
+   // space reaches 10 million.
+   const int max_dofs = 10000000;
+   for (int it = 1; it <= maxit; it++)
+   {
+      if (mpi.Root())
+      {
+         cout << "\nAMR Iteration " << it << endl;
+      }
+
+      // Display the current number of DoFs in each finite element space
+      Hertz.PrintSizes();
+
+      // Assemble all forms
+      Hertz.Assemble();
+
+      // Solve the system and compute any auxiliary fields
+      Hertz.Solve();
+
+      // Determine the current size of the linear system
+      int prob_size = Hertz.GetProblemSize();
+
+      // Write fields to disk for VisIt
+      if ( visit )
+      {
+         Hertz.WriteVisItFields(it);
+      }
+
+      // Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         Hertz.DisplayToGLVis();
+      }
+
+      if (mpi.Root())
+      {
+         cout << "AMR iteration " << it << " complete." << endl;
+      }
+
+      // Check stopping criteria
+      if (prob_size > max_dofs)
+      {
+         if (mpi.Root())
+         {
+            cout << "Reached maximum number of dofs, exiting..." << endl;
+         }
+         break;
+      }
+      if ( it == maxit )
+      {
+         break;
+      }
+
+      // Wait for user input. Ask every 10th iteration.
+      char c = 'c';
+      if (mpi.Root() && (it % 10 == 0))
+      {
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+      }
+      MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+      if (c != 'c')
+      {
+         break;
+      }
+
+      // Estimate element errors using the Zienkiewicz-Zhu error estimator.
+      Vector errors(pmesh.GetNE());
+      Hertz.GetErrorEstimates(errors);
+
+      double local_max_err = errors.Max();
+      double global_max_err;
+      MPI_Allreduce(&local_max_err, &global_max_err, 1,
+                    MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+
+      // Refine the elements whose error is larger than a fraction of the
+      // maximum element error.
+      const double frac = 0.5;
+      double threshold = frac * global_max_err;
+      if (mpi.Root()) { cout << "Refining ..." << endl; }
+      pmesh.RefineByError(errors, threshold);
+
+      // Update the magnetostatic solver to reflect the new state of the mesh.
+      Hertz.Update();
+
+      if (pmesh.Nonconforming() && mpi.WorldSize() > 1 && false)
+      {
+         if (mpi.Root()) { cout << "Rebalancing ..." << endl; }
+         pmesh.Rebalance();
+
+         // Update again after rebalancing
+         Hertz.Update();
+      }
+   }
+
+   // Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      Hertz.DisplayAnimationToGLVis();
+   }
+
+   delete epsCoef;
+   delete muInvCoef;
+   delete sigmaCoef;
+   delete etaInvCoef;
+
+   return 0;
+}
+
+// Print the Hertz ascii logo to the given ostream
+void display_banner(ostream & os)
+{
+   os << "     ____  ____              __           " << endl
+      << "    /   / /   / ____________/  |_________ " << endl
+      << "   /   /_/   /_/ __ \\_  __ \\   __\\___   / " << endl
+      << "  /   __    / \\  ___/|  | \\/|  |  /   _/  " << endl
+      << " /___/ /_  /   \\___  >__|   |__| /_____ \\ " << endl
+      << "         \\/        \\/                  \\/ " << endl << flush;
+}
+
+// The Permittivity is a required coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupPermittivityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( ds_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(dielectric_sphere);
+   }
+   else if ( pw_eps_.Size() > 0 )
+   {
+      coef = new PWConstCoefficient(pw_eps_);
+   }
+   else
+   {
+      coef = new ConstantCoefficient(epsilon0_);
+   }
+
+   return coef;
+}
+
+// The Permeability is a required coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupInvPermeabilityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( ms_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(magnetic_shell_inv);
+   }
+   else if ( pw_mu_.Size() > 0 )
+   {
+      pw_mu_inv_.SetSize(pw_mu_.Size());
+      for (int i = 0; i < pw_mu_.Size(); i++)
+      {
+         MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
+         pw_mu_inv_[i] = 1.0/pw_mu_[i];
+      }
+      coef = new PWConstCoefficient(pw_mu_inv_);
+   }
+   else
+   {
+      coef = new ConstantCoefficient(1.0/mu0_);
+   }
+
+   return coef;
+}
+
+// The Conductivity is an optional coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupConductivityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( cs_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(conductive_sphere);
+   }
+   else if ( pw_sigma_.Size() > 0 )
+   {
+      coef = new PWConstCoefficient(pw_sigma_);
+   }
+
+   return coef;
+}
+
+// The Admittance is an optional coefficient defined on boundary surfaces which
+// can be used in conjunction with absorbing boundary conditions.
+Coefficient *
+SetupAdmittanceCoefficient(const Mesh & mesh, const Array<int> & abcs)
+{
+   Coefficient * coef = NULL;
+
+   if ( pw_eta_.Size() > 0 )
+   {
+      MFEM_VERIFY(pw_eta_.Size() == abcs.Size(),
+                  "Each impedance value must be associated with exactly one "
+                  "absorbing boundary surface.");
+
+      pw_eta_inv_.SetSize(mesh.bdr_attributes.Size());
+
+      if ( abcs[0] == -1 )
+      {
+         pw_eta_inv_ = 1.0 / pw_eta_[0];
+      }
+      else
+      {
+         pw_eta_inv_ = 0.0;
+
+         for (int i=0; i<pw_eta_.Size(); i++)
+         {
+            pw_eta_inv_[abcs[i]-1] = 1.0 / pw_eta_[i];
+         }
+      }
+      coef = new PWConstCoefficient(pw_eta_inv_);
+   }
+
+   return coef;
+}
+
+// A sphere with constant permittivity.  The sphere has a radius,
+// center, and permittivity specified on the command line and stored
+// in ds_params_.
+double dielectric_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
+   }
+
+   if ( sqrt(r2) <= ds_params_(x.Size()) )
+   {
+      return ds_params_(x.Size()+1) * epsilon0_;
+   }
+   return epsilon0_;
+}
+
+// A spherical shell with constant permeability.  The sphere has inner
+// and outer radii, center, and relative permeability specified on the
+// command line and stored in ms_params_.
+double magnetic_shell(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i = 0; i < x.Size(); i++)
+   {
+      r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
+   }
+
+   if ( sqrt(r2) >= ms_params_(x.Size()) &&
+        sqrt(r2) <= ms_params_(x.Size()+1) )
+   {
+      return mu0_*ms_params_(x.Size()+2);
+   }
+   return mu0_;
+}
+
+// A sphere with constant conductivity.  The sphere has a radius,
+// center, and conductivity specified on the command line and stored
+// in ls_params_.
+double conductive_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
+   }
+
+   if ( sqrt(r2) <= cs_params_(x.Size()) )
+   {
+      return cs_params_(x.Size()+1);
+   }
+   return 0.0;
+}
+
+// A cylindrical rod of current density.  The rod has two axis end
+// points, a radus, a current amplitude in Amperes.  All of these
+// parameters are stored in do_params_.
+void dipole_oscillator(const Vector &x, Vector &j)
+{
+   MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
+
+   j.SetSize(x.Size());
+   j = 0.0;
+
+   Vector  v(x.Size());  // Normalized Axis vector
+   Vector xu(x.Size());  // x vector relative to the axis end-point
+
+   xu = x;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      xu[i] -= do_params_[i];
+      v[i]   = do_params_[x.Size()+i] - do_params_[i];
+   }
+
+   double h = v.Norml2();
+
+   if ( h == 0.0 )
+   {
+      return;
+   }
+   v /= h;
+
+   double r = do_params_[2*x.Size()+0];
+   double a = do_params_[2*x.Size()+1];
+
+   double xv = xu * v;
+
+   // Compute perpendicular vector from axis to x
+   xu.Add(-xv, v);
+
+   double xp = xu.Norml2();
+
+   if ( xv >= 0.0 && xv <= h && xp <= r )
+   {
+      j.Add(a, v);
+   }
+}
+
+void e_bc_r(const Vector &x, Vector &E)
+{
+   E.SetSize(3);
+   E = 0.0;
+}
+
+void e_bc_i(const Vector &x, Vector &E)
+{
+   E.SetSize(3);
+   E = 0.0;
+}

miniapps/electromagnetics/hertz_ser.cpp

@@ -0,0 +1,638 @@
+// 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.
+//
+//   -----------------------------------------------------------------------
+//   Hertz Miniapp:  Simple Frequency-Domain Electromagnetic Simulation Code
+//   -----------------------------------------------------------------------
+//
+//   Assumes that all sources and boundary conditions oscillate with the same
+//   frequency although not necessarily in phase with one another.  This
+//   assumptions implies that we can factor out the time dependence which we
+//   take to be of the form exp(i omega t).  With these assumptions we can
+//   write the Maxwell equations in the form:
+//
+//   i omega epsilon E = Curl mu^{-1} B - J - sigma E
+//   i omega B         = - Curl E
+//
+//   Which combine to yield:
+//
+//   Curl mu^{-1} Curl E - omega^2 epsilon E + i omega sigma E = - i omega J
+//
+//   We discretize this equation with H(Curl) a.k.a Nedelec basis
+//   functions.  The curl curl operator must be handled with
+//   integration by parts which yields a surface integral:
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               + (W, n x (mu^{-1} Curl E))_{\Gamma}
+//
+//   or
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               - i omega (W, n x H)_{\Gamma}
+//
+//   For plane waves
+//     omega B = - k x E
+//     omega D = k x H, assuming n x k = 0 => n x H = omega epsilon E / |k|
+//
+//   c = omega/|k|
+//
+//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
+//               - i omega sqrt{epsilon/mu} (W, E)_{\Gamma}
+//
+//
+// Compile with: make hertz
+//
+// Sample runs:
+//
+//   By default the sources and fields are all zero
+//     hertz_ser
+//
+//   Current source in a sphere with absorbing boundary conditions
+//     hertz_ser -m ../../data/ball-nurbs.mesh -rs 2
+//               -abcs '-1'
+//               -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//
+//   Current source in a metal sphere with dielectric and conducting materials
+//     hertz_ser -m ../../data/ball-nurbs.mesh -rs 2
+//               -dbcs '-1'
+//               -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
+//               -cs '0.0 0.0 -0.5 .2 10'
+//               -ds '0.0 0.0 0.5 .2 10'
+//
+//   Current source in a metal box
+//     hertz_ser -m ../../data/fichera.mesh -rs 3
+//               -dbcs '-1'
+//               -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//
+//   Current source with a mixture of absorbing and reflecting boundaries
+//     hertz_ser -m ../../data/fichera.mesh -rs 3
+//               -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
+//               -dbcs '4 8 19 21' -abcs '5 18'
+//
+
+#include "hertz_ser_solver.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::electromagnetics;
+
+// Permittivity Functions
+Coefficient * SetupPermittivityCoefficient();
+
+static Vector pw_eps_(0);     // Piecewise permittivity values
+static Vector ds_params_(0);  // Center, Radius, and Permittivity
+//                               of dielectric sphere
+double dielectric_sphere(const Vector &);
+
+// Permeability Function
+Coefficient * SetupInvPermeabilityCoefficient();
+
+static Vector pw_mu_(0);      // Piecewise permeability values
+static Vector pw_mu_inv_(0);  // Piecewise inverse permeability values
+static Vector ms_params_(0);  // Center, Inner and Outer Radii, and
+//                               Permeability of magnetic shell
+double magnetic_shell(const Vector &);
+double magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
+
+// Conductivity Functions
+Coefficient * SetupConductivityCoefficient();
+
+static Vector pw_sigma_(0);   // Piecewise conductivity values
+static Vector cs_params_(0);  // Center, Radius, and Conductivity
+//                               of conductive sphere
+double conductive_sphere(const Vector &);
+
+// Impedance
+Coefficient * SetupAdmittanceCoefficient(const Mesh & mesh,
+                                         const Array<int> & abcs);
+
+static Vector pw_eta_(0);      // Piecewise impedance values
+static Vector pw_eta_inv_(0);  // Piecewise inverse impedance values
+
+// Current Density Function
+static Vector do_params_(0);  // Axis Start, Axis End, Rod Radius,
+//                               Total Current of Rod
+void dipole_oscillator(const Vector &x, Vector &j);
+void j_src(const Vector &x, Vector &j) { dipole_oscillator(x, j); }
+
+// Electric Field Boundary Condition: The following function returns zero but
+// any function could be used.
+void e_bc_r(const Vector &x, Vector &E);
+void e_bc_i(const Vector &x, Vector &E);
+
+static double freq_ = 1.0;
+
+// Prints the program's logo to the given output stream
+void display_banner(ostream & os);
+
+int main(int argc, char *argv[])
+{
+   display_banner(cout);
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/ball-nurbs.mesh";
+   int order = 1;
+   int maxit = 100;
+   int serial_ref_levels = 0;
+   int sol = 1;
+   bool herm_conv = false;
+   bool visualization = true;
+   bool visit = true;
+
+   Array<int> abcs;
+   Array<int> dbcs;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
+                  "Number of serial refinement levels.");
+   args.AddOption(&freq_, "-f", "--frequency",
+                  "Frequency in Hertz (of course...)");
+   args.AddOption(&sol, "-s", "--solver",
+                  "Solver: 1 - GMRES");
+   args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
+                  "Piecewise values of Permittivity");
+   args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
+                  "Center, Radius, and Permittivity of Dielectric Sphere");
+   args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
+                  "Piecewise values of Permeability");
+   args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
+                  "Center, Inner Radius, Outer Radius, "
+                  "and Permeability of Magnetic Shell");
+   args.AddOption(&pw_sigma_, "-pws", "--piecewise-sigma",
+                  "Piecewise values of Conductivity");
+   args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
+                  "Center, Radius, and Conductivity of Conductive Sphere");
+   args.AddOption(&pw_eta_, "-pwz", "--piecewise-eta",
+                  "Piecewise values of Impedance (one value per abc surface)");
+   args.AddOption(&do_params_, "-do", "--dipole-oscillator-params",
+                  "Axis End Points, Radius, and Amplitude");
+   args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
+                  "Absorbing Boundary Condition Surfaces");
+   args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
+                  "Dirichlet Boundary Condition Surfaces");
+   /*
+   args.AddOption(&dbcv, "-dbcv", "--dirichlet-bc-vals",
+                  "Dirichlet Boundary Condition Values");
+   args.AddOption(&dbcg, "-dbcg", "--dirichlet-bc-gradient",
+                  "-no-dbcg", "--no-dirichlet-bc-gradient",
+                  "Dirichlet Boundary Condition Gradient (phi = -z)");
+   args.AddOption(&nbcs, "-nbcs", "--neumann-bc-surf",
+                  "Neumann Boundary Condition Surfaces");
+   args.AddOption(&nbcv, "-nbcv", "--neumann-bc-vals",
+                  "Neumann Boundary Condition Values");
+   args.AddOption(&kbcs, "-kbcs", "--surface-current-bc",
+                  "Surfaces for the Surface Current (K) Boundary Condition");
+   args.AddOption(&vbcs, "-vbcs", "--voltage-bc-surf",
+                  "Voltage Boundary Condition Surfaces (to drive K)");
+   args.AddOption(&vbcv, "-vbcv", "--voltage-bc-vals",
+                  "Voltage Boundary Condition Values (to drive K)");
+   */
+   args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
+                  "Max number of iterations in the main AMR loop.");
+   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
+                  "--no-hermitian", "Use convention for Hermitian operators.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
+                  "Enable or disable VisIt visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   ComplexOperator::Convention conv =
+      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
+
+   // 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 = new Mesh(mesh_file, 1, 1);
+
+   cout << "Starting initialization." << endl;
+
+   // Project a NURBS mesh to a piecewise-quadratic curved mesh
+   if (mesh->NURBSext)
+   {
+      mesh->UniformRefinement();
+      if (serial_ref_levels > 0) { serial_ref_levels--; }
+
+      mesh->SetCurvature(2);
+   }
+
+   // Ensure that quad and hex meshes are treated as non-conforming.
+   mesh->EnsureNCMesh();
+
+   // Refine the serial mesh on all processors to increase the resolution. In
+   // this example we do 'ref_levels' of uniform refinement.
+   for (int l = 0; l < serial_ref_levels; l++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // If values for Voltage BCs were not set issue a warning and exit
+   /*
+   if ( ( vbcs.Size() > 0 && kbcs.Size() == 0 ) ||
+        ( kbcs.Size() > 0 && vbcs.Size() == 0 ) ||
+        ( vbcv.Size() < vbcs.Size() ) )
+   {
+      if ( mpi.Root() )
+      {
+         cout << "The surface current (K) boundary condition requires "
+              << "surface current boundary condition surfaces (with -kbcs), "
+              << "voltage boundary condition surface (with -vbcs), "
+              << "and voltage boundary condition values (with -vbcv)."
+              << endl;
+      }
+      return 3;
+   }
+   */
+   // Create a coefficient describing the dielectric permittivity
+   Coefficient * epsCoef = SetupPermittivityCoefficient();
+
+   // Create a coefficient describing the magnetic permeability
+   Coefficient * muInvCoef = SetupInvPermeabilityCoefficient();
+
+   // Create a coefficient describing the electrical conductivity
+   Coefficient * sigmaCoef = SetupConductivityCoefficient();
+
+   // Create a coefficient describing the surface admittance
+   Coefficient * etaInvCoef = SetupAdmittanceCoefficient(*mesh, abcs);
+
+   // Create the Magnetostatic solver
+   HertzSolver Hertz(*mesh, order, freq_, (HertzSolver::SolverType)sol,
+                     conv, *epsCoef, *muInvCoef, sigmaCoef, etaInvCoef,
+                     abcs, dbcs,
+                     e_bc_r, e_bc_i,
+                     (do_params_.Size() > 0 ) ? j_src : NULL, NULL
+                    );
+
+   //(b_uniform_.Size() > 0 ) ? a_bc_uniform  : NULL,
+   //(cr_params_.Size() > 0 ) ? current_ring  : NULL,
+   //(bm_params_.Size() > 0 ) ? bar_magnet    :
+   //(ha_params_.Size() > 0 ) ? halbach_array : NULL);
+
+   // Initialize GLVis visualization
+   if (visualization)
+   {
+      Hertz.InitializeGLVis();
+   }
+
+   // Initialize VisIt visualization
+   VisItDataCollection visit_dc("Hertz-AMR-Serial", mesh);
+
+   if ( visit )
+   {
+      Hertz.RegisterVisItFields(visit_dc);
+   }
+
+   // ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
+   ZienkiewiczZhuEstimator estimator(*Hertz.GetErrorInteg(),
+                                     *Hertz.GetErrorField(),
+                                     Hertz.GetErrorFluxFES());
+   ThresholdRefiner refiner(estimator);
+   refiner.SetTotalErrorFraction(0.5);
+
+   cout << "Initialization done." << endl;
+
+   // The main AMR loop. In each iteration we solve the problem on the current
+   // mesh, visualize the solution, estimate the error on all elements, refine
+   // the worst elements and update all objects to work with the new mesh. We
+   // refine until the maximum number of dofs in the Nedelec finite element
+   // space reaches 10 million.
+   const int max_dofs = 10000000;
+   for (int it = 1; it <= maxit; it++)
+   {
+      cout << "\nAMR Iteration " << it << endl;
+
+      // Display the current number of DoFs in each finite element space
+      Hertz.PrintSizes();
+
+      // Assemble all forms
+      Hertz.Assemble();
+
+      // Solve the system and compute any auxiliary fields
+      Hertz.Solve();
+
+      // Determine the current size of the linear system
+      int prob_size = Hertz.GetProblemSize();
+
+      // Write fields to disk for VisIt
+      if ( visit )
+      {
+         Hertz.WriteVisItFields(it);
+      }
+
+      // Send the solution by socket to a GLVis server.
+      if (visualization)
+      {
+         Hertz.DisplayToGLVis();
+      }
+
+      cout << "AMR iteration " << it << " complete." << endl;
+
+      // Check stopping criteria
+      if (prob_size > max_dofs)
+      {
+         cout << "Reached maximum number of dofs, exiting..." << endl;
+         break;
+      }
+      if ( it == maxit )
+      {
+         break;
+      }
+
+      // Wait for user input. Ask every 10th iteration.
+      char c = 'c';
+      if (it % 10 == 0)
+      {
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+      }
+
+      if (c != 'c')
+      {
+         break;
+      }
+
+      // Estimate element errors using the Zienkiewicz-Zhu error estimator.
+      refiner.Apply(*mesh);
+      if (refiner.Stop())
+      {
+         cout << "Stopping criterion satisfied. Stop." << endl;
+         break;
+      }
+      /*
+      Vector errors(mesh.GetNE());
+      Hertz.GetErrorEstimates(errors);
+
+      double global_max_err = errors.Max();
+
+      // Refine the elements whose error is larger than a fraction of the
+      // maximum element error.
+      const double frac = 0.5;
+      double threshold = frac * global_max_err;
+      cout << "Refining ..." << endl;
+      mesh.RefineByError(errors, threshold);
+      */
+
+      // Update the magnetostatic solver to reflect the new state of the mesh.
+      Hertz.Update();
+
+   }
+
+   delete epsCoef;
+   delete muInvCoef;
+   delete sigmaCoef;
+
+   return 0;
+}
+
+// Print the Hertz ascii logo to the given ostream
+void display_banner(ostream & os)
+{
+   os << "     ____  ____              __           " << endl
+      << "    /   / /   / ____________/  |_________ " << endl
+      << "   /   /_/   /_/ __ \\_  __ \\   __\\___   / " << endl
+      << "  /   __    / \\  ___/|  | \\/|  |  /   _/  " << endl
+      << " /___/ /_  /   \\___  >__|   |__| /_____ \\ " << endl
+      << "         \\/        \\/                  \\/ " << endl << flush;
+}
+
+// The Permittivity is a required coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupPermittivityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( ds_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(dielectric_sphere);
+   }
+   else if ( pw_eps_.Size() > 0 )
+   {
+      coef = new PWConstCoefficient(pw_eps_);
+   }
+   else
+   {
+      coef = new ConstantCoefficient(epsilon0_);
+   }
+
+   return coef;
+}
+
+// The Permeability is a required coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupInvPermeabilityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( ms_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(magnetic_shell_inv);
+   }
+   else if ( pw_mu_.Size() > 0 )
+   {
+      pw_mu_inv_.SetSize(pw_mu_.Size());
+      for (int i = 0; i < pw_mu_.Size(); i++)
+      {
+         MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
+         pw_mu_inv_[i] = 1.0/pw_mu_[i];
+      }
+      coef = new PWConstCoefficient(pw_mu_inv_);
+   }
+   else
+   {
+      coef = new ConstantCoefficient(1.0/mu0_);
+   }
+
+   return coef;
+}
+
+// The Conductivity is an optional coefficient which may be defined in
+// various ways so we'll determine the appropriate coefficient type here.
+Coefficient *
+SetupConductivityCoefficient()
+{
+   Coefficient * coef = NULL;
+
+   if ( cs_params_.Size() > 0 )
+   {
+      coef = new FunctionCoefficient(conductive_sphere);
+   }
+   else if ( pw_sigma_.Size() > 0 )
+   {
+      coef = new PWConstCoefficient(pw_sigma_);
+   }
+
+   return coef;
+}
+
+// The Admittance is an optional coefficient defined on boundary surfaces which
+// can be used in conjunction with absorbing boundary conditions.
+Coefficient *
+SetupAdmittanceCoefficient(const Mesh & mesh, const Array<int> & abcs)
+{
+   Coefficient * coef = NULL;
+
+   if ( pw_eta_.Size() > 0 )
+   {
+      MFEM_VERIFY(pw_eta_.Size() == abcs.Size(),
+                  "Each impedance value must be associated with exactly one "
+                  "absorbing boundary surface.");
+
+      pw_eta_inv_.SetSize(mesh.bdr_attributes.Size());
+
+      if ( abcs[0] == -1 )
+      {
+         pw_eta_inv_ = 1.0 / pw_eta_[0];
+      }
+      else
+      {
+         pw_eta_inv_ = 0.0;
+
+         for (int i=0; i<pw_eta_.Size(); i++)
+         {
+            pw_eta_inv_[abcs[i]-1] = 1.0 / pw_eta_[i];
+         }
+      }
+      coef = new PWConstCoefficient(pw_eta_inv_);
+   }
+
+   return coef;
+}
+
+// A sphere with constant permittivity.  The sphere has a radius,
+// center, and permittivity specified on the command line and stored
+// in ds_params_.
+double dielectric_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
+   }
+
+   if ( sqrt(r2) <= ds_params_(x.Size()) )
+   {
+      return ds_params_(x.Size()+1) * epsilon0_;
+   }
+   return epsilon0_;
+}
+
+// A spherical shell with constant permeability.  The sphere has inner
+// and outer radii, center, and relative permeability specified on the
+// command line and stored in ms_params_.
+double magnetic_shell(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i = 0; i < x.Size(); i++)
+   {
+      r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
+   }
+
+   if ( sqrt(r2) >= ms_params_(x.Size()) &&
+        sqrt(r2) <= ms_params_(x.Size()+1) )
+   {
+      return mu0_*ms_params_(x.Size()+2);
+   }
+   return mu0_;
+}
+
+// A sphere with constant conductivity.  The sphere has a radius,
+// center, and conductivity specified on the command line and stored
+// in ls_params_.
+double conductive_sphere(const Vector &x)
+{
+   double r2 = 0.0;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
+   }
+
+   if ( sqrt(r2) <= cs_params_(x.Size()) )
+   {
+      return cs_params_(x.Size()+1);
+   }
+   return 0.0;
+}
+
+// A cylindrical rod of current density.  The rod has two axis end
+// points, a radus, a current amplitude in Amperes.  All of these
+// parameters are stored in do_params_.
+void dipole_oscillator(const Vector &x, Vector &j)
+{
+   MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
+
+   j.SetSize(x.Size());
+   j = 0.0;
+
+   Vector  v(x.Size());  // Normalized Axis vector
+   Vector xu(x.Size());  // x vector relative to the axis end-point
+
+   xu = x;
+
+   for (int i=0; i<x.Size(); i++)
+   {
+      xu[i] -= do_params_[i];
+      v[i]   = do_params_[x.Size()+i] - do_params_[i];
+   }
+
+   double h = v.Norml2();
+
+   if ( h == 0.0 )
+   {
+      return;
+   }
+   v /= h;
+
+   double r = do_params_[2*x.Size()+0];
+   double a = do_params_[2*x.Size()+1];
+
+   double xv = xu * v;
+
+   // Compute perpendicular vector from axis to x
+   xu.Add(-xv, v);
+
+   double xp = xu.Norml2();
+
+   if ( xv >= 0.0 && xv <= h && xp <= r )
+   {
+      j.Add(a, v);
+   }
+}
+
+void e_bc_r(const Vector &x, Vector &E)
+{
+   E.SetSize(3);
+   E = 0.0;
+
+}
+
+void e_bc_i(const Vector &x, Vector &E)
+{
+   E.SetSize(3);
+   E = 0.0;
+}

miniapps/electromagnetics/hertz_ser_solver.hpp

@@ -0,0 +1,248 @@
+// 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_HERTZ_SOLVER
+#define MFEM_HERTZ_SOLVER
+
+#include "../common/fem_extras.hpp"
+
+#include <string>
+#include <map>
+
+namespace mfem
+{
+
+using common::H1_FESpace;
+using common::ND_FESpace;
+using common::RT_FESpace;
+using common::DiscreteGradOperator;
+using common::DiscreteCurlOperator;
+//using common::DivergenceFreeProjector;
+
+namespace electromagnetics
+{
+
+// Physical Constants
+// Permittivity of Free Space (units F/m)
+static const double epsilon0_ = 8.8541878176e-12;
+// Permeability of Free Space (units H/m)
+static const double mu0_ = 4.0e-7*M_PI;
+
+//class SurfaceCurrent;
+class HertzSolver
+{
+public:
+
+   enum SolverType
+   {
+      INVALID   = -1,
+      GMRES     =  1,
+      FGMRES    =  2,
+      MINRES    =  3,
+      SUPERLU   =  4,
+      STRUMPACK =  5
+   };
+
+   HertzSolver(Mesh & mesh, int order, double freq,
+               HertzSolver::SolverType s,
+               ComplexOperator::Convention conv,
+               Coefficient & epsCoef,
+               Coefficient & muInvCoef,
+               Coefficient * sigmaCoef,
+               Coefficient * etaInvCoef,
+               Array<int> & abcs,
+               Array<int> & dbcs,
+               void   (*e_r_bc )(const Vector&, Vector&),
+               void   (*e_i_bc )(const Vector&, Vector&),
+               void   (*j_r_src)(const Vector&, Vector&),
+               void   (*j_i_src)(const Vector&, Vector&));
+   ~HertzSolver();
+
+   int GetProblemSize();
+
+   void PrintSizes();
+
+   void Assemble();
+
+   void Update();
+
+   void Solve();
+
+   //void GetErrorEstimates(Vector & errors);
+
+   void RegisterVisItFields(VisItDataCollection & visit_dc);
+
+   void WriteVisItFields(int it = 0);
+
+   void InitializeGLVis();
+
+   void DisplayToGLVis();
+
+   BilinearFormIntegrator * GetErrorInteg() { return err_integ_; }
+
+   GridFunction * GetErrorField() { return &e_->real(); }
+
+   FiniteElementSpace * GetErrorFluxFES()
+   {
+      return new FiniteElementSpace(mesh_, err_flux_fec_);
+   }
+   // const ParGridFunction & GetVectorPotential() { return *a_; }
+
+private:
+
+   int order_;
+   int logging_;
+
+   SolverType sol_;
+
+   ComplexOperator::Convention conv_;
+
+   bool ownsEtaInv_;
+
+   double freq_;
+
+   Mesh * mesh_;
+
+   // H1_ParFESpace * H1FESpace_;
+   ND_FESpace * HCurlFESpace_;
+   // RT_ParFESpace * HDivFESpace_;
+
+   Array<int> blockTrueOffsets_;
+
+   // ParSesquilinearForm * a0_;
+   SesquilinearForm * a1_;
+   BilinearForm * b1_;
+
+   // ParGridFunction * e_r_;  // Real part of electric field (HCurl)
+   // ParGridFunction * e_i_;  // Imaginary part of electric field (HCurl)
+   ComplexGridFunction * e_;  // Complex electric field (HCurl)
+   ComplexGridFunction * j_;  // Complex current density (HCurl)
+   // ParGridFunction * j_i_;  // Imaginary part of current density (HCurl)
+
+   ComplexLinearForm   * jd_; // Dual of complex current density (HCurl)
+   // ParLinearForm   * jd_r_; // Dual of real part of current density (HCurl)
+   // ParLinearForm   * jd_i_; // Dual of imaginary part of current density (HCurl)
+
+   /*
+     ParBilinearForm * curlMuInvCurl_;
+     ParBilinearForm * hCurlMass_;
+     ParMixedBilinearForm * hDivHCurlMuInv_;
+     ParMixedBilinearForm * weakCurlMuInv_;
+    */
+   // ParDiscreteGradOperator * grad_;
+   // ParDiscreteCurlOperator * curl_;
+   /*
+    ParGridFunction * a_;  // Vector Potential (HCurl)
+    ParGridFunction * b_;  // Magnetic Flux (HDiv)
+    ParGridFunction * h_;  // Magnetic Field (HCurl)
+    ParGridFunction * jr_; // Raw Volumetric Current Density (HCurl)
+    ParGridFunction * j_;  // Volumetric Current Density (HCurl)
+    ParGridFunction * k_;  // Surface Current Density (HCurl)
+    ParGridFunction * m_;  // Magnetization (HDiv)
+    ParGridFunction * bd_; // Dual of B (HCurl)
+    ParGridFunction * jd_; // Dual of J, the rhs vector (HCurl)
+   */
+   // DivergenceFreeProjector * DivFreeProj_;
+   // SurfaceCurrent          * SurfCur_;
+
+   Coefficient       * epsCoef_;   // Dielectric Material Coefficient
+   Coefficient       * muInvCoef_; // Dia/Paramagnetic Material Coefficient
+   Coefficient       * sigmaCoef_; // Electrical Conductivity Coefficient
+   Coefficient       * etaInvCoef_; // Admittance Coefficient
+
+   Coefficient * omegaCoef_;     // omega expressed as a Coefficient
+   Coefficient * negOmegaCoef_;  // -omega expressed as a Coefficient
+   Coefficient * omega2Coef_;    // omega^2 expressed as a Coefficient
+   Coefficient * negOmega2Coef_; // -omega^2 expressed as a Coefficient
+   Coefficient * massCoef_;      // -omega^2 epsilon
+   Coefficient * posMassCoef_;   // omega^2 epsilon
+   Coefficient * lossCoef_;      // -omega sigma
+   // Coefficient * gainCoef_;    // omega sigma
+   Coefficient * abcCoef_;       // -omega eta^{-1}
+
+   // VectorCoefficient * aBCCoef_;   // Vector Potential BC Function
+   VectorCoefficient * jrCoef_;     // Volume Current Density Function
+   VectorCoefficient * jiCoef_;     // Volume Current Density Function
+   VectorCoefficient * erCoef_;     // Electric Field Boundary Condition
+   VectorCoefficient * eiCoef_;     // Electric Field Boundary Condition
+   // VectorCoefficient * mCoef_;     // Magnetization Vector Function
+
+   BilinearFormIntegrator * err_integ_;
+   FiniteElementCollection * err_flux_fec_;
+
+   // void   (*a_bc_ )(const Vector&, Vector&);
+   void   (*j_r_src_)(const Vector&, Vector&);
+   void   (*j_i_src_)(const Vector&, Vector&);
+   // void   (*m_src_)(const Vector&, Vector&);
+
+   // Array of 0's and 1's marking the location of absorbing surfaces
+   Array<int> abc_marker_;
+
+   // Array of 0's and 1's marking the location of Dirichlet boundaries
+   Array<int> dbc_marker_;
+   void   (*e_r_bc_)(const Vector&, Vector&);
+   void   (*e_i_bc_)(const Vector&, Vector&);
+
+   Array<int> * dbcs_;
+   Array<int> ess_bdr_;
+   Array<int> ess_bdr_tdofs_;
+   Array<int> non_k_bdr_;
+
+   VisItDataCollection * visit_dc_;
+
+   std::map<std::string,socketstream*> socks_;
+};
+/*
+class SurfaceCurrent
+{
+public:
+ SurfaceCurrent(ParFiniteElementSpace & H1FESpace,
+                ParDiscreteGradOperator & Grad,
+                Array<int> & kbcs, Array<int> & vbcs, Vector & vbcv);
+ ~SurfaceCurrent();
+
+ void InitSolver() const;
+
+ void ComputeSurfaceCurrent(ParGridFunction & k);
+
+ void Update();
+
+ ParGridFunction * GetPsi() { return psi_; }
+
+private:
+ int myid_;
+
+ ParFiniteElementSpace   * H1FESpace_;
+ ParDiscreteGradOperator * grad_;
+ Array<int>              * kbcs_;
+ Array<int>              * vbcs_;
+ Vector                  * vbcv_;
+
+ ParBilinearForm * s0_;
+ ParGridFunction * psi_;
+ ParGridFunction * rhs_;
+
+ HypreParMatrix  * S0_;
+ mutable Vector Psi_;
+ mutable Vector RHS_;
+
+ mutable HypreBoomerAMG  * amg_;
+ mutable HyprePCG        * pcg_;
+
+ Array<int> ess_bdr_, ess_bdr_tdofs_;
+ Array<int> non_k_bdr_;
+};
+*/
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_HERTZ_SOLVER

miniapps/electromagnetics/hertz_solver.cpp

@@ -0,0 +1,711 @@
+// 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.
+
+#include "hertz_solver.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace std;
+namespace mfem
+{
+using namespace common;
+
+namespace electromagnetics
+{
+
+// Used for combining scalar coefficients
+double prodFunc(double a, double b) { return a * b; }
+
+HertzSolver::HertzSolver(ParMesh & pmesh, int order, double freq,
+                         HertzSolver::SolverType sol, SolverOptions & sOpts,
+                         HertzSolver::PrecondType prec,
+                         ComplexOperator::Convention conv,
+                         Coefficient & epsCoef,
+                         Coefficient & muInvCoef,
+                         Coefficient * sigmaCoef,
+                         Coefficient * etaInvCoef,
+                         Array<int> & abcs,
+                         Array<int> & dbcs,
+                         void (*e_r_bc )(const Vector&, Vector&),
+                         void (*e_i_bc )(const Vector&, Vector&),
+                         void (*j_r_src)(const Vector&, Vector&),
+                         void (*j_i_src)(const Vector&, Vector&))
+   : myid_(0),
+     num_procs_(1),
+     order_(order),
+     logging_(1),
+     sol_(sol),
+     solOpts_(sOpts),
+     prec_(prec),
+     conv_(conv),
+     ownsEtaInv_(etaInvCoef == NULL),
+     freq_(freq),
+     pmesh_(&pmesh),
+     HCurlFESpace_(NULL),
+     a1_(NULL),
+     b1_(NULL),
+     e_(NULL),
+     e_t_(NULL),
+     j_(NULL),
+     jd_(NULL),
+     epsCoef_(&epsCoef),
+     muInvCoef_(&muInvCoef),
+     sigmaCoef_(sigmaCoef),
+     etaInvCoef_(etaInvCoef),
+     omegaCoef_(new ConstantCoefficient(2.0 * M_PI * freq_)),
+     negOmegaCoef_(new ConstantCoefficient(-2.0 * M_PI * freq_)),
+     omega2Coef_(new ConstantCoefficient(pow(2.0 * M_PI * freq_, 2))),
+     negOmega2Coef_(new ConstantCoefficient(-pow(2.0 * M_PI * freq_, 2))),
+     massCoef_(NULL),
+     posMassCoef_(NULL),
+     lossCoef_(NULL),
+     abcCoef_(NULL),
+     posAbcCoef_(NULL),
+     jrCoef_(NULL),
+     jiCoef_(NULL),
+     erCoef_(NULL),
+     eiCoef_(NULL),
+     j_r_src_(j_r_src),
+     j_i_src_(j_i_src),
+     e_r_bc_(e_r_bc),
+     e_i_bc_(e_i_bc),
+     dbcs_(&dbcs),
+     visit_dc_(NULL)
+{
+   // Initialize MPI variables
+   MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
+   MPI_Comm_rank(pmesh_->GetComm(), &myid_);
+
+   // Define compatible parallel finite element spaces on the parallel
+   // mesh. Here we use arbitrary order Nedelec finite elements.
+   HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
+
+   // Set the size of the 2x2 block representation of the complex linear system
+   blockTrueOffsets_.SetSize(3);
+   blockTrueOffsets_[0] = 0;
+   blockTrueOffsets_[1] = HCurlFESpace_->TrueVSize();
+   blockTrueOffsets_[2] = HCurlFESpace_->TrueVSize();
+   blockTrueOffsets_.PartialSum();
+
+   // Setup Dirichlet BC
+   ess_bdr_.SetSize(pmesh.bdr_attributes.Max());
+   if ( dbcs_ != NULL )
+   {
+      if ( dbcs_->Size() == 1 && (*dbcs_)[0] == -1 )
+      {
+         ess_bdr_ = 1;
+      }
+      else
+      {
+         ess_bdr_ = 0;
+         for (int i=0; i<dbcs_->Size(); i++)
+         {
+            ess_bdr_[(*dbcs_)[i]-1] = 1;
+         }
+      }
+      HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
+
+      if (e_r_bc_)
+      {
+         erCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
+                                                 e_r_bc_);
+         if (e_i_bc_ == NULL)
+         {
+            Vector e(3); e = 0.0;
+            eiCoef_ = new VectorConstantCoefficient(e);
+         }
+      }
+      if (e_i_bc_)
+      {
+         eiCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
+                                                 e_i_bc_);
+         if (e_r_bc_ == NULL)
+         {
+            Vector e(3); e = 0.0;
+            erCoef_ = new VectorConstantCoefficient(e);
+         }
+      }
+   }
+
+   // Setup various coefficients
+   massCoef_ = new TransformedCoefficient(negOmega2Coef_, epsCoef_, prodFunc);
+   posMassCoef_ = new TransformedCoefficient(omega2Coef_, epsCoef_, prodFunc);
+   if ( sigmaCoef_ )
+   {
+      lossCoef_ = new TransformedCoefficient(omegaCoef_, sigmaCoef_, prodFunc);
+   }
+
+   // Impedance of free space for the Absorbing boundary condition
+   if ( abcs.Size() > 0 )
+   {
+      if ( myid_ == 0 && logging_ > 0 )
+      {
+         cout << "Creating Admittance Coefficient" << endl;
+      }
+
+      abc_marker_.SetSize(pmesh.bdr_attributes.Max());
+      if ( abcs.Size() == 1 && abcs[0] < 0 )
+      {
+         // Mark all boundaries as absorbing
+         abc_marker_ = 1;
+      }
+      else
+      {
+         // Mark select boundaries as absorbing
+         abc_marker_ = 0;
+         for (int i=0; i<abcs.Size(); i++)
+         {
+            abc_marker_[abcs[i]-1] = 1;
+         }
+      }
+      if ( etaInvCoef_ == NULL )
+      {
+         etaInvCoef_ = new ConstantCoefficient(sqrt(epsilon0_/mu0_));
+      }
+      abcCoef_ = new TransformedCoefficient(negOmegaCoef_, etaInvCoef_,
+                                            prodFunc);
+      posAbcCoef_ = new TransformedCoefficient(omegaCoef_, etaInvCoef_,
+                                               prodFunc);
+   }
+
+   // Volume Current Density
+   if ( j_r_src_ != NULL )
+   {
+      jrCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
+                                              j_r_src_);
+   }
+   else
+   {
+      Vector j(3); j = 0.0;
+      jrCoef_ = new VectorConstantCoefficient(j);
+   }
+   if ( j_i_src_ != NULL )
+   {
+      jiCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
+                                              j_i_src_);
+   }
+   else
+   {
+      Vector j(3); j = 0.0;
+      jiCoef_ = new VectorConstantCoefficient(j);
+   }
+
+   // Bilinear Forms
+   // Primary system operator
+   a1_ = new ParSesquilinearForm(HCurlFESpace_, conv_);
+   a1_->AddDomainIntegrator(new CurlCurlIntegrator(*muInvCoef_), NULL);
+   a1_->AddDomainIntegrator(new VectorFEMassIntegrator(*massCoef_), NULL);
+   if ( lossCoef_ )
+   {
+      a1_->AddDomainIntegrator(NULL, new VectorFEMassIntegrator(*lossCoef_));
+   }
+   if ( abcCoef_ )
+   {
+      a1_->AddBoundaryIntegrator(NULL, new VectorFEMassIntegrator(*abcCoef_),
+                                 abc_marker_);
+   }
+
+   // Operator used with the perconditioner
+   b1_ = new ParBilinearForm(HCurlFESpace_);
+   b1_->AddDomainIntegrator(new CurlCurlIntegrator(*muInvCoef_));
+   b1_->AddDomainIntegrator(new VectorFEMassIntegrator(*posMassCoef_));
+   if ( lossCoef_ )
+   {
+      b1_->AddDomainIntegrator(new VectorFEMassIntegrator(*lossCoef_));
+   }
+   if ( abcCoef_ )
+   {
+      b1_->AddBoundaryIntegrator(new VectorFEMassIntegrator(*posAbcCoef_),
+                                 abc_marker_);
+   }
+
+   // Build grid functions
+   // The solution vector is the Electric field
+   e_   = new ParComplexGridFunction(HCurlFESpace_);
+   e_t_ = new ParGridFunction(HCurlFESpace_);
+   if (erCoef_ && eiCoef_)
+   {
+      e_->ProjectCoefficient(*erCoef_, *eiCoef_);
+   }
+   else
+   {
+      *e_ = 0.0;
+   }
+
+   // A GridFunction to visualize the volumetric current density
+   j_ = new ParComplexGridFunction(HCurlFESpace_);
+   j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
+
+   // A LineatForm representation of the current denisty for the RHS
+   jd_ = new ParComplexLinearForm(HCurlFESpace_, conv_);
+   jd_->AddDomainIntegrator(new VectorFEDomainLFIntegrator(*jrCoef_),
+                            new VectorFEDomainLFIntegrator(*jiCoef_));
+   jd_->real().Vector::operator=(0.0);
+   jd_->imag().Vector::operator=(0.0);
+}
+
+HertzSolver::~HertzSolver()
+{
+   delete jrCoef_;
+   delete jiCoef_;
+   delete erCoef_;
+   delete eiCoef_;
+   delete massCoef_;
+   delete posMassCoef_;
+   delete lossCoef_;
+   delete abcCoef_;
+   delete posAbcCoef_;
+   if ( ownsEtaInv_ ) { delete etaInvCoef_; }
+   delete omegaCoef_;
+   delete negOmegaCoef_;
+   delete omega2Coef_;
+   delete negOmega2Coef_;
+
+   delete e_;
+   delete e_t_;
+   delete j_;
+   delete jd_;
+
+   delete a1_;
+   delete b1_;
+
+   delete HCurlFESpace_;
+
+   map<string,socketstream*>::iterator mit;
+   for (mit=socks_.begin(); mit!=socks_.end(); mit++)
+   {
+      delete mit->second;
+   }
+}
+
+HYPRE_Int
+HertzSolver::GetProblemSize()
+{
+   return 2 * HCurlFESpace_->GlobalTrueVSize();
+}
+
+void
+HertzSolver::PrintSizes()
+{
+   HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
+   if (myid_ == 0)
+   {
+      cout << "Number of H(Curl) unknowns: " << size_nd << endl;
+   }
+}
+
+void
+HertzSolver::Assemble()
+{
+   if ( myid_ == 0 && logging_ > 0 ) { cout << "Assembling ..." << flush; }
+
+   a1_->Assemble();
+   a1_->Finalize();
+
+   b1_->Assemble();
+   b1_->Finalize();
+
+   jd_->Assemble();
+
+   if ( myid_ == 0 && logging_ > 0 ) { cout << " done." << endl; }
+}
+
+void
+HertzSolver::Update()
+{
+   if ( myid_ == 0 && logging_ > 0 ) { cout << "Updating ..." << endl; }
+
+   // Inform the spaces that the mesh has changed
+   HCurlFESpace_->Update();
+
+   if ( ess_bdr_.Size() > 0 )
+   {
+      HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
+   }
+
+   blockTrueOffsets_[0] = 0;
+   blockTrueOffsets_[1] = HCurlFESpace_->TrueVSize();
+   blockTrueOffsets_[2] = HCurlFESpace_->TrueVSize();
+   blockTrueOffsets_.PartialSum();
+
+   // Inform the grid functions that the space has changed.
+   e_->Update();
+   if (erCoef_ && eiCoef_)
+   {
+      e_->ProjectCoefficient(*erCoef_, *eiCoef_);
+   }
+
+   j_->Update();
+   j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
+
+   jd_->Update();
+
+   // Inform the bilinear forms that the space has changed.
+   a1_->Update();
+   b1_->Update();
+}
+
+void
+HertzSolver::Solve()
+{
+   if ( myid_ == 0 && logging_ > 0 ) { cout << "Running solver ... " << endl; }
+
+   OperatorHandle A1;
+   Vector E, RHS;
+
+   a1_->FormLinearSystem(ess_bdr_tdofs_, *e_, *jd_, A1, E, RHS);
+
+   OperatorHandle PCOp;
+   b1_->FormSystemMatrix(ess_bdr_tdofs_, PCOp);
+
+   tic_toc.Clear();
+   tic_toc.Start();
+
+   Operator * pcr = NULL;
+   Operator * pci = NULL;
+   BlockDiagonalPreconditioner * BDP = NULL;
+
+   if (sol_ == FGMRES || sol_ == MINRES)
+   {
+      switch (prec_)
+      {
+         case INVALID_PC:
+            if ( myid_ == 0 && logging_ > 0 )
+            {
+               cout << "No Preconditioner Requested" << endl;
+            }
+            break;
+         case DIAG_SCALE:
+            if ( myid_ == 0 && logging_ > 0 )
+            {
+               cout << "Diagonal Scaling Preconditioner Requested" << endl;
+            }
+            pcr = new HypreDiagScale(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
+            break;
+         case PARASAILS:
+            if ( myid_ == 0 && logging_ > 0 )
+            {
+               cout << "ParaSails Preconditioner Requested" << endl;
+            }
+            pcr = new HypreParaSails(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
+            dynamic_cast<HypreParaSails*>(pcr)->SetSymmetry(1);
+            break;
+         case EUCLID:
+            if ( myid_ == 0 && logging_ > 0 )
+            {
+               cout << "Euclid Preconditioner Requested" << endl;
+            }
+            pcr = new HypreEuclid(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
+            if (solOpts_.euLvl != 1)
+            {
+               HypreSolver * pc = dynamic_cast<HypreSolver*>(pcr);
+               HYPRE_EuclidSetLevel(*pc, solOpts_.euLvl);
+            }
+            break;
+         case AMS:
+            if ( myid_ == 0 && logging_ > 0 )
+            {
+               cout << "AMS Preconditioner Requested" << endl;
+            }
+            pcr = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
+                               HCurlFESpace_);
+            break;
+         default:
+            MFEM_ABORT("Requested preconditioner is not available.");
+            break;
+      }
+      pci = pcr;
+
+      if (pcr)
+      {
+         BDP = new BlockDiagonalPreconditioner(blockTrueOffsets_);
+         BDP->SetDiagonalBlock(0, pcr);
+         BDP->SetDiagonalBlock(1, pci);
+         BDP->owns_blocks = 0;
+      }
+   }
+
+   switch (sol_)
+   {
+      case GMRES:
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "GMRES Solver Requested" << endl;
+         }
+         GMRESSolver gmres(HCurlFESpace_->GetComm());
+         gmres.SetOperator(*A1.Ptr());
+         gmres.SetRelTol(solOpts_.relTol);
+         gmres.SetMaxIter(solOpts_.maxIter);
+         gmres.SetKDim(solOpts_.kDim);
+         gmres.SetPrintLevel(solOpts_.printLvl);
+
+         gmres.Mult(RHS, E);
+      }
+      break;
+      case FGMRES:
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "FGMRES Solver Requested" << endl;
+         }
+         FGMRESSolver fgmres(HCurlFESpace_->GetComm());
+         if (BDP) { fgmres.SetPreconditioner(*BDP); }
+         fgmres.SetOperator(*A1.Ptr());
+         fgmres.SetRelTol(solOpts_.relTol);
+         fgmres.SetMaxIter(solOpts_.maxIter);
+         fgmres.SetKDim(solOpts_.kDim);
+         fgmres.SetPrintLevel(solOpts_.printLvl);
+
+         fgmres.Mult(RHS, E);
+      }
+      break;
+      case MINRES:
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "MINRES Solver Requested" << endl;
+         }
+         MINRESSolver minres(HCurlFESpace_->GetComm());
+         if (BDP) { minres.SetPreconditioner(*BDP); }
+         minres.SetOperator(*A1.Ptr());
+         minres.SetRelTol(solOpts_.relTol);
+         minres.SetMaxIter(solOpts_.maxIter);
+         minres.SetPrintLevel(solOpts_.printLvl);
+
+         minres.Mult(RHS, E);
+      }
+      break;
+#ifdef MFEM_USE_SUPERLU
+      case SUPERLU:
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "SuperLU Solver Requested" << endl;
+         }
+         ComplexHypreParMatrix * A1Z = A1.As<ComplexHypreParMatrix>();
+         HypreParMatrix * A1C = A1Z->GetSystemMatrix();
+         SuperLURowLocMatrix A_SuperLU(*A1C);
+         SuperLUSolver solver(MPI_COMM_WORLD);
+         solver.SetOperator(A_SuperLU);
+         solver.Mult(RHS, E);
+         delete A1C;
+      }
+      break;
+#endif
+#ifdef MFEM_USE_STRUMPACK
+      case STRUMPACK:
+      {
+         if ( myid_ == 0 && logging_ > 0 )
+         {
+            cout << "STRUMPACK Solver Requested" << endl;
+         }
+         ComplexHypreParMatrix * A1Z = A1.As<ComplexHypreParMatrix>();
+         HypreParMatrix * A1C = A1Z->GetSystemMatrix();
+         STRUMPACKRowLocMatrix A_STRUMPACK(*A1C);
+         STRUMPACKSolver solver(0, NULL, MPI_COMM_WORLD);
+         solver.SetPrintFactorStatistics(true);
+         solver.SetPrintSolveStatistics(false);
+         solver.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         solver.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+         solver.DisableMatching();
+         solver.SetOperator(A_STRUMPACK);
+         solver.SetFromCommandLine();
+         solver.Mult(RHS, E);
+         delete A1C;
+      }
+      break;
+#endif
+      default:
+         break;
+   };
+
+   tic_toc.Stop();
+
+   e_->Distribute(E);
+
+   delete BDP;
+   if (pci != pcr) { delete pci; }
+   delete pcr;
+
+   if ( myid_ == 0 && logging_ > 0 )
+   {
+      cout << " Solver done in " << tic_toc.RealTime() << " seconds." << endl;
+   }
+}
+
+void
+HertzSolver::GetErrorEstimates(Vector & errors)
+{
+   if ( myid_ == 0 && logging_ > 0 )
+   { cout << "Estimating Error ... " << flush; }
+
+   // Space for the discontinuous (original) flux
+   CurlCurlIntegrator flux_integrator(*muInvCoef_);
+   RT_FECollection flux_fec(order_-1, pmesh_->SpaceDimension());
+   ParFiniteElementSpace flux_fes(pmesh_, &flux_fec);
+
+   // Space for the smoothed (conforming) flux
+   double norm_p = 1;
+   ND_FECollection smooth_flux_fec(order_, pmesh_->Dimension());
+   ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
+
+   L2ZZErrorEstimator(flux_integrator, e_->real(),
+                      smooth_flux_fes, flux_fes, errors, norm_p);
+
+   if ( myid_ == 0 && logging_ > 0 ) { cout << "done." << endl; }
+}
+
+void
+HertzSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
+{
+   visit_dc_ = &visit_dc;
+
+   visit_dc.RegisterField("Re(E)", &e_->real());
+   visit_dc.RegisterField("Im(E)", &e_->imag());
+
+   if ( j_ )
+   {
+      visit_dc.RegisterField("Re(J)", &j_->real());
+      visit_dc.RegisterField("Im(J)", &j_->imag());
+   }
+}
+
+void
+HertzSolver::WriteVisItFields(int it)
+{
+   if ( visit_dc_ )
+   {
+      if (myid_ == 0) { cout << "Writing VisIt files ..." << flush; }
+
+      if ( j_ )
+      {
+         j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
+      }
+
+      HYPRE_Int prob_size = this->GetProblemSize();
+      visit_dc_->SetCycle(it);
+      visit_dc_->SetTime(prob_size);
+      visit_dc_->Save();
+
+      if (myid_ == 0) { cout << " done." << endl; }
+   }
+}
+
+void
+HertzSolver::InitializeGLVis()
+{
+   if ( myid_ == 0 ) { cout << "Opening GLVis sockets." << endl; }
+
+   socks_["Er"] = new socketstream;
+   socks_["Er"]->precision(8);
+
+   socks_["Ei"] = new socketstream;
+   socks_["Ei"]->precision(8);
+
+   if ( j_ )
+   {
+      socks_["Jr"] = new socketstream;
+      socks_["Jr"]->precision(8);
+
+      socks_["Ji"] = new socketstream;
+      socks_["Ji"]->precision(8);
+   }
+
+   if ( myid_ == 0 ) { cout << "GLVis sockets open." << endl; }
+}
+
+void
+HertzSolver::DisplayToGLVis()
+{
+   if (myid_ == 0) { cout << "Sending data to GLVis ..." << flush; }
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+
+   int Wx = 0, Wy = 0; // window position
+   int Ww = 350, Wh = 350; // window size
+   int offx = Ww+10, offy = Wh+45; // window offsets
+
+   VisualizeField(*socks_["Er"], vishost, visport,
+                  e_->real(), "Electric Field, Re(E)", Wx, Wy, Ww, Wh);
+   Wx += offx;
+
+   VisualizeField(*socks_["Ei"], vishost, visport,
+                  e_->imag(), "Electric Field, Im(E)", Wx, Wy, Ww, Wh);
+
+   Wx = 0; Wy += offy; // next line
+
+   if ( j_ )
+   {
+      j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
+
+      VisualizeField(*socks_["Jr"], vishost, visport,
+                     j_->real(), "Current Density, Re(J)", Wx, Wy, Ww, Wh);
+      Wx += offx;
+      VisualizeField(*socks_["Ji"], vishost, visport,
+                     j_->imag(), "Current Density, Im(J)", Wx, Wy, Ww, Wh);
+   }
+
+   Wx = 0; Wy += offy; // next line
+
+   if (myid_ == 0) { cout << " done." << endl; }
+}
+
+void
+HertzSolver::DisplayAnimationToGLVis()
+{
+   if (myid_ == 0) { cout << "Sending animation data to GLVis ..." << flush; }
+
+   Vector zeroVec(3); zeroVec = 0.0;
+   VectorConstantCoefficient zeroCoef(zeroVec);
+
+   double norm_r = e_->real().ComputeMaxError(zeroCoef);
+   double norm_i = e_->imag().ComputeMaxError(zeroCoef);
+
+   *e_t_ = e_->real();
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock << "parallel " << num_procs_ << " " << myid_ << "\n";
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << *pmesh_ << *e_t_
+            << "window_title 'Harmonic Solution (t = 0.0 T)'"
+            << "valuerange 0.0 " << max(norm_r, norm_i) << "\n"
+            << "autoscale off\n"
+            << "keys cvvv\n"
+            << "pause\n" << flush;
+   if (myid_ == 0)
+      cout << "GLVis visualization paused."
+           << " Press space (in the GLVis window) to resume it.\n";
+   int num_frames = 24;
+   int i = 0;
+   while (sol_sock)
+   {
+      double t = (double)(i % num_frames) / num_frames;
+      ostringstream oss;
+      oss << "Harmonic Solution (t = " << t << " T)";
+
+      add( cos( 2.0 * M_PI * t), e_->real(),
+           sin( 2.0 * M_PI * t), e_->imag(), *e_t_);
+      sol_sock << "parallel " << num_procs_ << " " << myid_ << "\n";
+      sol_sock << "solution\n" << *pmesh_ << *e_t_
+               << "window_title '" << oss.str() << "'" << flush;
+      i++;
+   }
+}
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI

miniapps/electromagnetics/hertz_solver.hpp

@@ -0,0 +1,194 @@
+// 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_HERTZ_SOLVER
+#define MFEM_HERTZ_SOLVER
+
+#include "../common/pfem_extras.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include <string>
+#include <map>
+
+namespace mfem
+{
+
+using common::H1_ParFESpace;
+using common::ND_ParFESpace;
+using common::RT_ParFESpace;
+using common::ParDiscreteGradOperator;
+using common::ParDiscreteCurlOperator;
+using common::DivergenceFreeProjector;
+
+namespace electromagnetics
+{
+
+// Physical Constants
+// Permittivity of Free Space (units F/m)
+static const double epsilon0_ = 8.8541878176e-12;
+// Permeability of Free Space (units H/m)
+static const double mu0_ = 4.0e-7*M_PI;
+
+// Solver options
+struct SolverOptions
+{
+   int maxIter;
+   int kDim;
+   int printLvl;
+   double relTol;
+
+   // Euclid Options
+   int euLvl;
+};
+
+//class SurfaceCurrent;
+class HertzSolver
+{
+public:
+
+   enum PrecondType
+   {
+      INVALID_PC  = -1,
+      DIAG_SCALE  =  1,
+      PARASAILS   =  2,
+      EUCLID      =  3,
+      AMS         =  4
+   };
+
+   enum SolverType
+   {
+      INVALID   = -1,
+      GMRES     =  1,
+      FGMRES    =  2,
+      MINRES    =  3,
+      SUPERLU   =  4,
+      STRUMPACK =  5
+   };
+
+   HertzSolver(ParMesh & pmesh, int order, double freq,
+               HertzSolver::SolverType s, SolverOptions & sOpts,
+               HertzSolver::PrecondType p,
+               ComplexOperator::Convention conv,
+               Coefficient & epsCoef,
+               Coefficient & muInvCoef,
+               Coefficient * sigmaCoef,
+               Coefficient * etaInvCoef,
+               Array<int> & abcs,
+               Array<int> & dbcs,
+               void (*e_r_bc )(const Vector&, Vector&),
+               void (*e_i_bc )(const Vector&, Vector&),
+               void (*j_r_src)(const Vector&, Vector&),
+               void (*j_i_src)(const Vector&, Vector&));
+   ~HertzSolver();
+
+   HYPRE_Int GetProblemSize();
+
+   void PrintSizes();
+
+   void Assemble();
+
+   void Update();
+
+   void Solve();
+
+   void GetErrorEstimates(Vector & errors);
+
+   void RegisterVisItFields(VisItDataCollection & visit_dc);
+
+   void WriteVisItFields(int it = 0);
+
+   void InitializeGLVis();
+
+   void DisplayToGLVis();
+
+   void DisplayAnimationToGLVis();
+
+private:
+
+   int myid_;
+   int num_procs_;
+   int order_;
+   int logging_;
+
+   SolverType sol_;
+   SolverOptions & solOpts_;
+   PrecondType prec_;
+
+   ComplexOperator::Convention conv_;
+
+   bool ownsEtaInv_;
+
+   double freq_;
+
+   ParMesh * pmesh_;
+
+   ND_ParFESpace * HCurlFESpace_;
+
+   Array<HYPRE_Int> blockTrueOffsets_;
+
+   ParSesquilinearForm * a1_;
+   ParBilinearForm * b1_;
+
+   ParComplexGridFunction * e_;   // Complex electric field (HCurl)
+   ParGridFunction        * e_t_; // Real electric field (HCurl)
+   ParComplexGridFunction * j_;   // Complex current density (HCurl)
+
+   ParComplexLinearForm   * jd_; // Dual of complex current density (HCurl)
+
+   Coefficient       * epsCoef_;    // Dielectric Material Coefficient
+   Coefficient       * muInvCoef_;  // Dia/Paramagnetic Material Coefficient
+   Coefficient       * sigmaCoef_;  // Electrical Conductivity Coefficient
+   Coefficient       * etaInvCoef_; // Admittance Coefficient
+
+   Coefficient * omegaCoef_;     // omega expressed as a Coefficient
+   Coefficient * negOmegaCoef_;  // -omega expressed as a Coefficient
+   Coefficient * omega2Coef_;    // omega^2 expressed as a Coefficient
+   Coefficient * negOmega2Coef_; // -omega^2 expressed as a Coefficient
+   Coefficient * massCoef_;      // -omega^2 epsilon
+   Coefficient * posMassCoef_;   // omega^2 epsilon
+   Coefficient * lossCoef_;      // -omega sigma
+   Coefficient * abcCoef_;       // -omega eta^{-1}
+   Coefficient * posAbcCoef_;    // omega eta^{-1}
+
+   VectorCoefficient * jrCoef_;     // Volume Current Density Function
+   VectorCoefficient * jiCoef_;     // Volume Current Density Function
+   VectorCoefficient * erCoef_;     // Electric Field Boundary Condition
+   VectorCoefficient * eiCoef_;     // Electric Field Boundary Condition
+
+   void (*j_r_src_)(const Vector&, Vector&);
+   void (*j_i_src_)(const Vector&, Vector&);
+
+   // Array of 0's and 1's marking the location of absorbing surfaces
+   Array<int> abc_marker_;
+
+   // Array of 0's and 1's marking the location of Dirichlet boundaries
+   Array<int> dbc_marker_;
+   void (*e_r_bc_)(const Vector&, Vector&);
+   void (*e_i_bc_)(const Vector&, Vector&);
+
+   Array<int> * dbcs_;
+   Array<int> ess_bdr_;
+   Array<int> ess_bdr_tdofs_;
+   Array<int> non_k_bdr_;
+
+   VisItDataCollection * visit_dc_;
+
+   std::map<std::string,socketstream*> socks_;
+};
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_HERTZ_SOLVER

miniapps/electromagnetics/makefile

@@ -25,8 +25,8 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_MINIAPPS =
-PAR_MINIAPPS = volta tesla maxwell joule
+SEQ_MINIAPPS = hertz_ser
+PAR_MINIAPPS = volta tesla maxwell hertz joule
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -118,4 +118,6 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_*
+	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_* \
+	Hertz-AMR-Parallel* Hertz-AMR-Serial*
+

miniapps/electromagnetics/volta_solver.cpp

@@ -448,7 +448,7 @@ VoltaSolver::Solve()
    divEpsGrad_->RecoverFEMSolution(Phi, *rhod_, *phi_);
 
    // Compute the negative Gradient of the solution vector.  This is
-   // the magnetic field corresponding to the scalar potential
+   // the electric field corresponding to the scalar potential
    // represented by phi.
    grad_->Mult(*phi_, *e_); *e_ *= -1.0;