Cleans up thermal dielectric components.

# Conflicts:
#	doc/CodeDocumentation.dox
#	makefile

CHANGELOG

@@ -312,6 +312,102 @@ Miscellaneous
 
 - Various other simplifications, extensions, and bugfixes in the code.
 
+- Added PA support for MixedScalarCurlIntegrator in 2D and
+  MixedVectorGradientIntegrator in 2D and 3D, as well as their transposes.
+
+- Added hipSPARSE support for sparse mat-vec multiplications.
+
+- Added support for using the HYPRE library built with HIP support. Similar to
+  the HYPRE + CUDA support added earlier, most of the MFEM examples and miniapps
+  work transparently with HYPRE + HIP builds. This includes the BoomerAMG, AMS,
+  and ADS solvers.
+
+- More explicit and consistent formating of the output of iterative solvers
+  with the new IterativeSolver::PrintLevel options. See linalg/solvers.hpp.
+
+- Added a miniapp for PDE-based extrapolation of finite element functions. See
+  miniapps/shifted/extrapolate.cpp.
+
+- Added support for automatic differentiation. Users can select between native
+  implementation and external library implementation during configuration. One
+  parallel and two serial examples are implemented in the miniapps/autodiff/
+  directory.
+
+- GridFunctionCoefficient (and the related vector, gradient, divergence, and
+  curl classes) now work properly with LORDiscretization and LORSolver.
+
+- Added support for mesh preprocessing to resolve fine scale problem data
+  before simulation. This feature uses adaptive mesh refinement to control the
+  associated data oscillation error. See the new Example 30/30p.
+
+- Switched from Artistic Style (astyle) version 2.05.1 to version 3.1 for code
+  formatting. See the "make style" target.
+
+- Split the fem/fe.?pp files into separate files in the new fem/fe/ directory
+  to simplify and clarify the organization of FiniteElement classes.
+
+- Added support for hr-adaptivity using TMOP-based error estimator.
+
+- Coefficient::SetTime now propagates the new time into internally stored
+  Coefficient objects.
+
+- Added initial support for google-benchmarks in the tests/benchmarks directory.
+  It can be enabled with MFEM_USE_BENCHMARK=YES.
+
+- Added Binder (mybinder.org) configuration files for C++ MFEM Jupyter Notebooks
+  with inline GLVis visualization as well as a new examples/jupyter/ directory
+  with a sample notebook based on Example 1. Implementation based on xeus-cling,
+  github.com/jupyter-xeus/xeus-cling + xeus-glvis, github.com/GLVis/xeus-glvis.
+
+- Added 'double' atomicAdd implementation for previous versions of CUDA.
+
+- Adding lowest order Nedelec and Raviart-Thomas basis functions on wedge
+  shaped elements.
+
+- Added initial support for meshes with pyramidal elements, including several
+  pyramidal meshes in the data/ directory and support for the lowest order H1,
+  Nedelec, Raviart-Thomas, and L2 basis functions on pyramids.
+
+- Updated the hypre interface according to changes in hypre-2.22.1. The ADS
+  solver is now fully working on GPUs.
+
+- Tetrahedral meshes no longer need to be reordered to support high order
+  Nedelec basis functions. This will allow future support for Nedelec basis
+  functions on wedges and pyramids which are not amenable to reordering. The
+  ReorientTetMesh method of the Mesh and ParMesh classes has been deprecated.
+
+- Gmsh meshes where all elements have zero physical tag (the default Gmsh
+  output format if no physical groups are defined) are now successfully loaded,
+  and elements are reassigned attribute number 1.
+
+- Added new miniapps that use the ParELAG library, its hybrid smoothers, and the
+  hierarchy of spaces created by the element-based AMG (AMGe) methodology in
+  ParELAG to build multigrid solvers for H(curl) and H(div) forms. See the
+  miniapps/parelag directory for more details.
+
+- Fixed several MinGW build issues on Windows.
+
+- Remove the 'u' flag in the ar command, to update all files in the archive,
+  avoiding file name collisions from different subdirectories.
+
+- Added initial TMOP-based capabilities for surface fitting and tangential
+  relaxation in the mesh-optimizer and pmesh-optimizer miniapps.
+
+- Added ParMesh Adjaceny Set (adjset) creation support to the Conduit Mesh
+  Blueprint MFEM wrapper functions in ConduitDataCollection.
+
+- `HypreParVector` and `Vector` now support move semantics, and the copy
+  constructor for `HypreParVector` now copies the local vector data.
+
+- The HPC versions of ex1 and ex1p (in miniapps/performance) now support
+  runtime selection of either 2D or 3D meshes.
+
+- Added arbitrary order Nedelec and Raviart-Thomas basis functions for
+  wedge-shaped elements.
+
+- Added ParaView visualization of `QuadratureFunction` fields, through both
+  `QuadratureFunction::SaveVTU` and `ParaViewDataCollection::RegisterQField`.
+
 
 Version 4.4, released on March 21, 2022
 =======================================
@@ -493,6 +589,12 @@ Discretization improvements
 
 - Added support for nonscalar coefficient with VectorDiffusionIntegrator.
 
+- Added support for Partial Assembly and Element Assembly with Discontinuous
+  Galerkin methods on nonconforming meshes.
+
+- Added a simpler interface to request face information: see
+  `Mesh::FaceInformation` and `Mesh::GetFaceInformation`.
+
 Linear and nonlinear solvers
 ----------------------------
 - Added support for AMG preconditioners on GPUs based on the hypre library

fem/coefficient.cpp

@@ -1041,6 +1041,73 @@ double DeterminantCoefficient::Eval(ElementTransformation &T,
    return ma.Det();
 }
 
+VectorComponentCoefficient::VectorComponentCoefficient(VectorCoefficient &A,
+                                                       int c)
+   : a(&A), va(A.GetVDim())
+{
+   SetComponent(c);
+}
+
+void VectorComponentCoefficient::SetComponent(int c)
+{
+   MFEM_ASSERT(c < a->GetVDim() && c >= 0,
+               "VectorComponentCoefficient:  "
+               "Index not in range.");
+
+   component = c;
+}
+
+void VectorComponentCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+double VectorComponentCoefficient::Eval(ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(va, T, ip);
+   return va[component];
+}
+
+MatrixComponentCoefficient::MatrixComponentCoefficient(MatrixCoefficient &A,
+                                                       int ri, int ci)
+   : a(&A), ma(A.GetHeight(), A.GetWidth())
+{
+   SetRowIndex(ri);
+   SetColumnIndex(ci);
+}
+
+void MatrixComponentCoefficient::SetRowIndex(int ri)
+{
+   MFEM_ASSERT(ri < a->GetHeight() && ri >= 0,
+               "MatrixComponentCoefficient:  "
+               "Row index not in range.");
+
+   row_idx = ri;
+}
+
+void MatrixComponentCoefficient::SetColumnIndex(int ci)
+{
+   MFEM_ASSERT(ci < a->GetWidth() && ci >= 0,
+               "MatrixComponentCoefficient:  "
+               "Column index not in range.");
+   col_idx = ci;
+}
+
+void MatrixComponentCoefficient::SetTime(double t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+double MatrixComponentCoefficient::Eval(ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   return ma(row_idx,col_idx);
+}
+
 VectorSumCoefficient::VectorSumCoefficient(int dim)
    : VectorCoefficient(dim),
      ACoef(NULL), BCoef(NULL),
@@ -1196,6 +1263,7 @@ void MatrixVectorProductCoefficient::SetTime(double t)
 void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
                                           const IntegrationPoint &ip)
 {
+   V.SetSize(ma.Height());
    a->Eval(ma, T, ip);
    b->Eval(vb, T, ip);
    V.SetSize(vdim);

fem/coefficient.hpp

@@ -1672,6 +1672,62 @@ public:
    { return pow(a->Eval(T, ip), p); }
 };
 
+/// Coefficient which returns (k*x) or func(k*x) where k is a vector
+class PhaseCoefficient : public Coefficient
+{
+private:
+   double(*func_)(double);
+   VectorCoefficient * k_;
+   mutable Vector kVec_;
+
+public:
+   PhaseCoefficient(Vector & k, double(*func)(double) = NULL)
+      : func_(func), k_(NULL), kVec_(k) {}
+
+   PhaseCoefficient(VectorCoefficient & k, double(*func)(double) = NULL)
+      : func_(func), k_(&k), kVec_(k.GetVDim()) {}
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      if (k_) { k_->Eval(kVec_, T, ip); }
+      return (func_) ? (*func_)(kVec_ * transip) : (kVec_ * transip);
+   }
+};
+
+/// Coefficient which returns func(kr*x)*exp(-ki*x) where kr and ki are vectors
+class ComplexPhaseCoefficient : public Coefficient
+{
+private:
+   double(*func_)(double);
+   VectorCoefficient * kr_;
+   VectorCoefficient * ki_;
+   mutable Vector krVec_;
+   mutable Vector kiVec_;
+
+public:
+   ComplexPhaseCoefficient(Vector & kr, Vector & ki, double(&func)(double))
+      : func_(&func), kr_(NULL), ki_(NULL), krVec_(kr), kiVec_(ki) {}
+
+   ComplexPhaseCoefficient(VectorCoefficient & kr, VectorCoefficient & ki,
+                           double(&func)(double))
+      : func_(&func), kr_(&kr), ki_(&ki),
+        krVec_(kr.GetVDim()), kiVec_(ki.GetVDim()) {}
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+      T.Transform(ip, transip);
+      if (kr_) { kr_->Eval(krVec_, T, ip); }
+      if (ki_) { ki_->Eval(kiVec_, T, ip); }
+      return (*func_)(krVec_ * transip)*exp(-(kiVec_ * transip));
+   }
+};
 
 /// Scalar coefficient defined as the inner product of two vector coefficients
 class InnerProductCoefficient : public Coefficient
@@ -1761,6 +1817,85 @@ public:
                        const IntegrationPoint &ip);
 };
 
+/// Scalar coefficient defined as component of a vector coefficient
+class VectorComponentCoefficient : public Coefficient
+{
+private:
+   VectorCoefficient *a = nullptr;
+
+   mutable Vector va;
+   int component;
+
+public:
+   /// Construct with a vector coefficient.
+   VectorComponentCoefficient(VectorCoefficient &A)
+      : a(&A), va(A.GetVDim()), component(0) {};
+
+   VectorComponentCoefficient(VectorCoefficient &A, int c);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(double t) override;
+
+   /// Reset the vector coefficient
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+
+   /// Return the vector coefficient
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Set the component
+   void SetComponent(int c);
+
+   /// Return the component
+   int GetComponent() const { return component; }
+
+   /// Evaluate the trace coefficient at @a ip.
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip) override;
+};
+
+/// Scalar coefficient defined as component of a matrix coefficient
+class MatrixComponentCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient *a = nullptr;
+
+   mutable DenseMatrix ma;
+   int row_idx,col_idx;
+
+public:
+   MatrixComponentCoefficient(MatrixCoefficient &A)
+      : a(&A), ma(A.GetHeight(), A.GetWidth()), row_idx(0), col_idx(0) {};
+
+   /// Construct with the matrix coefficient.
+   MatrixComponentCoefficient(MatrixCoefficient &A, int ri, int ci);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(double t) override;
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Reset the index
+   void SetRowIndex(int ri);
+
+   /// Return the index
+   int GetRowIndex() const { return row_idx; }
+
+   /// Reset the index
+   void SetColumnIndex(int ci);
+
+   /// Return the index
+   int GetColumnIndex() const { return col_idx; }
+
+
+   /// Evaluate the trace coefficient at @a ip.
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip) override;
+};
+
 /// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {

fem/complex_fem.cpp

@@ -1415,6 +1415,192 @@ ParSesquilinearForm::Update(FiniteElementSpace *nfes)
    if ( pblfi ) { pblfi->Update(nfes); }
 }
 
+bool ParMixedSesquilinearForm::RealInteg()
+{
+   int nint = pblfr->GetTFBFI()->Size() + pblfr->GetDBFI()->Size() +
+              pblfr->GetBBFI()->Size() + pblfr->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+bool ParMixedSesquilinearForm::ImagInteg()
+{
+   int nint = pblfi->GetTFBFI()->Size() + pblfi->GetDBFI()->Size() +
+              pblfi->GetBBFI()->Size() + pblfi->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                                                   ParFiniteElementSpace *te_pf,
+                                                   ComplexOperator::Convention
+                                                   convention)
+   : conv(convention),
+     pblfr(new ParMixedBilinearForm(tr_pf, te_pf)),
+     pblfi(new ParMixedBilinearForm(tr_pf, te_pf))
+{}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                                                   ParFiniteElementSpace *te_pf,
+                                                   ParMixedBilinearForm *pbfr,
+                                                   ParMixedBilinearForm *pbfi,
+                                                   ComplexOperator::Convention convention)
+   : conv(convention),
+     pblfr(new ParMixedBilinearForm(tr_pf,te_pf,pbfr)),
+     pblfi(new ParMixedBilinearForm(tr_pf,te_pf,pbfi))
+{}
+
+ParMixedSesquilinearForm::~ParMixedSesquilinearForm()
+{
+   delete pblfr;
+   delete pblfi;
+}
+
+void ParMixedSesquilinearForm::Mult(const ParComplexGridFunction & x,
+                                    ParComplexLinearForm & y) const
+{
+   pblfr->Mult(x.real(), y.real());
+   pblfi->AddMult(x.imag(), y.real(), -1.0);
+
+   pblfr->Mult(x.imag(), y.imag());
+   pblfi->AddMult(x.real(), y.imag(), 1.0);
+
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      y.imag() *= -1.0;
+   }
+
+}
+
+void ParMixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator
+                                                   *bfi_real,
+                                                   BilinearFormIntegrator
+                                                   *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddDomainIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddDomainIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator
+                                                *bfi_real,
+                                                BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator
+                                                *bfi_real,
+                                                BilinearFormIntegrator *bfi_imag,
+                                                Array<int> & bdr_marker)
+{
+   if (bfi_real) { pblfr->AddBoundaryIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { pblfi->AddBoundaryIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParMixedSesquilinearForm::AddTraceFaceIntegrator(BilinearFormIntegrator
+                                                 *bfi_real,
+                                                 BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                    *bfi_real,
+                                                    BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real) { pblfr->AddBdrTraceFaceIntegrator(bfi_real); }
+   if (bfi_imag) { pblfi->AddBdrTraceFaceIntegrator(bfi_imag); }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                    *bfi_real,
+                                                    BilinearFormIntegrator *bfi_imag,
+                                                    Array<int> &bdr_marker)
+{
+   if (bfi_real) { pblfr->AddBdrTraceFaceIntegrator(bfi_real, bdr_marker); }
+   if (bfi_imag) { pblfi->AddBdrTraceFaceIntegrator(bfi_imag, bdr_marker); }
+}
+
+void
+ParMixedSesquilinearForm::Assemble(int skip_zeros)
+{
+   pblfr->Assemble(skip_zeros);
+   pblfi->Assemble(skip_zeros);
+}
+
+void
+ParMixedSesquilinearForm::Finalize(int skip_zeros)
+{
+   pblfr->Finalize(skip_zeros);
+   pblfi->Finalize(skip_zeros);
+}
+
+ComplexHypreParMatrix *
+ParMixedSesquilinearForm::ParallelAssemble()
+{
+   return new ComplexHypreParMatrix(pblfr->ParallelAssemble(),
+                                    pblfi->ParallelAssemble(),
+                                    true, true, conv);
+}
+
+void ParMixedSesquilinearForm::FormRectangularSystemMatrix(
+   const Array<int> &trial_tdof_list,
+   const Array<int> &test_tdof_list,
+   OperatorHandle &A)
+{
+   OperatorHandle A_r, A_i;
+   if (RealInteg())
+   {
+      pblfr->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_r);
+   }
+   if (ImagInteg())
+   {
+      pblfi->FormRectangularSystemMatrix(trial_tdof_list, test_tdof_list, A_i);
+   }
+   if (!RealInteg() && !ImagInteg())
+   {
+      MFEM_ABORT("Both Real and Imaginary part of the MixedSesquilinear form are empty");
+   }
+
+   // A = A_r + i A_i
+   A.Clear();
+   if ( A_r.Type() == Operator::Hypre_ParCSR ||
+        A_i.Type() == Operator::Hypre_ParCSR )
+   {
+      ComplexHypreParMatrix * A_hyp =
+         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
+                                   A_i.As<HypreParMatrix>(),
+                                   A_r.OwnsOperator(),
+                                   A_i.OwnsOperator(),
+                                   conv);
+      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.As<Operator>(),
+                             A_i.As<Operator>(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
+void
+ParMixedSesquilinearForm::Update()
+{
+   if ( pblfr ) { pblfr->Update(); }
+   if ( pblfi ) { pblfi->Update(); }
+}
+
 #endif // MFEM_USE_MPI
 
 }

fem/complex_fem.hpp

@@ -682,6 +682,141 @@ public:
    virtual ~ParSesquilinearForm();
 };
 
+/** Class for a parallel mixed sesquilinear form
+
+    A sesquilinear form is a generalization of a mixed bilinear form to
+    complex-valued fields. Sesquilinear forms are linear in the second argument
+    but the first argument involves a complex conjugate in the sense that:
+
+                a(alpha u, beta v) = conj(alpha) beta a(u, v)
+
+    The @a convention argument in the class's constructor is documented in the
+    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
+
+    When supplying integrators to the ParSesquilinearForm either the real or
+    imaginary integrator can be NULL. This indicates that the corresponding
+    portion of the complex-valued material coefficient is equal to zero.
+*/
+class ParMixedSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv;
+
+   ParMixedBilinearForm *pblfr;
+   ParMixedBilinearForm *pblfi;
+
+   /* These methods check if the real/imag parts of the sesqulinear form are not
+      empty */
+   bool RealInteg();
+   bool ImagInteg();
+
+public:
+   ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                            ParFiniteElementSpace *te_pf,
+                            ComplexOperator::Convention
+                            convention = ComplexOperator::HERMITIAN);
+
+   /** @brief Create a ParMixedSesquilinearForm on the ParFiniteElementSpaces
+       @a tr_pf and @a te_pf, using the same integrators as the
+       ParMixedBilinearForms @a pbfr and @a pbfi .
+
+       The pointer @a pf is not owned by the newly constructed object.
+
+       The integrators are copied as pointers and they are not owned by the
+       newly constructed ParSesquilinearForm. */
+   ParMixedSesquilinearForm(ParFiniteElementSpace *tr_pf,
+                            ParFiniteElementSpace *te_pf,
+                            ParMixedBilinearForm *pbfr,
+                            ParMixedBilinearForm *pbfi,
+                            ComplexOperator::Convention
+                            convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv; }
+   void SetConvention(const ComplexOperator::Convention &
+                      convention) { conv = convention; }
+
+   /// Set the desired assembly level.
+   /** Valid choices are:
+
+       - AssemblyLevel::LEGACYFULL (default)
+       - AssemblyLevel::FULL
+       - AssemblyLevel::PARTIAL
+       - AssemblyLevel::ELEMENT
+       - AssemblyLevel::NONE
+
+       This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level)
+   {
+      pblfr->SetAssemblyLevel(assembly_level);
+      pblfi->SetAssemblyLevel(assembly_level);
+   }
+
+   ParMixedBilinearForm & real() { return *pblfr; }
+   ParMixedBilinearForm & imag() { return *pblfi; }
+   const ParMixedBilinearForm & real() const { return *pblfr; }
+   const ParMixedBilinearForm & imag() const { return *pblfi; }
+
+   /// Matrix multiplication: \f$ y = M x \f$
+   void Mult(const ParComplexGridFunction & x,
+             ParComplexLinearForm & y) const;
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi_real,
+                            BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag);
+
+   /** @brief Adds new boundary Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The array @a bdr_marker is stored internally as a pointer to the given
+       Array<int> object. */
+   void AddBoundaryIntegrator(BilinearFormIntegrator *bfi_real,
+                              BilinearFormIntegrator *bfi_imag,
+                              Array<int> &bdr_marker);
+
+   /// Adds new Face Integrator. Assumes ownership of @a bfi.
+   void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                               BilinearFormIntegrator *bfi_imag);
+
+   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag);
+
+   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The array @a bdr_marker is stored internally as a pointer to the given
+       Array<int> object. */
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi_real,
+                                  BilinearFormIntegrator *bfi_imag,
+                                  Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexHypreParMatrix *ParallelAssemble();
+
+   void FormRectangularSystemMatrix(const Array<int> &trial_tdof_list,
+                                    const Array<int> &test_tdof_list,
+                                    OperatorHandle &A);
+
+   virtual void Update();
+
+   virtual ~ParMixedSesquilinearForm();
+};
+
 #endif // MFEM_USE_MPI
 
 }

fem/estimators.cpp

@@ -498,4 +498,64 @@ void LpErrorEstimator::ComputeEstimates()
    current_sequence = sol->FESpace()->GetMesh()->GetSequence();
 }
 
+void ComplexLpErrorEstimator::ComputeEstimates()
+{
+   MFEM_VERIFY(real_coef != NULL || real_vcoef != NULL,
+               "ComplexLpErrorEstimator has no coefficient "
+               "for the real part!  "
+               "Call SetRealCoef first.");
+   MFEM_VERIFY(imag_coef != NULL || imag_vcoef != NULL,
+               "ComplexLpErrorEstimator has no coefficient "
+               "for the imaginary part!  "
+               "Call SetImagCoef first.");
+
+   int ne = 0;
+   if (sol) { ne = sol->FESpace()->GetMesh()->GetNE(); }
+#ifdef MFEM_USE_MPI
+   if (par_sol) { ne = par_sol->FESpace()->GetMesh()->GetNE(); }
+#endif
+
+   error_estimates.SetSize(ne);
+
+   const Vector & real_errors = real_estimator.GetLocalErrors();
+   const Vector & imag_errors = imag_estimator.GetLocalErrors();
+
+   if (local_norm_p < infinity())
+   {
+      for (int i=0; i<ne; i++)
+      {
+         const double re = pow(real_errors[i], local_norm_p);
+         const double ie = pow(imag_errors[i], local_norm_p);
+
+         error_estimates[i] = pow(re + ie, 1./local_norm_p);
+      }
+   }
+   else
+   {
+      for (int i=0; i<ne; i++)
+      {
+         error_estimates[i] = std::max(real_errors[i], imag_errors[i]);
+      }
+   }
+   /*
+   #ifdef MFEM_USE_MPI
+   total_error = error_estimates.Sum();
+   auto pfes = dynamic_cast<ParFiniteElementSpace*>(sol->FESpace());
+   if (pfes)
+   {
+      auto process_local_error = total_error;
+      MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
+                    MPI_SUM, pfes->GetComm());
+   }
+   #endif // MFEM_USE_MPI
+   total_error = pow(total_error, 1.0/local_norm_p);
+   */
+   current_sequence = -1;
+   if (sol) { current_sequence = sol->FESpace()->GetMesh()->GetSequence(); }
+#ifdef MFEM_USE_MPI
+   if (par_sol)
+   { current_sequence = par_sol->FESpace()->GetMesh()->GetSequence(); }
+#endif
+}
+
 } // namespace mfem

fem/estimators.hpp

@@ -17,6 +17,7 @@
 #include "../config/config.hpp"
 #include "../linalg/vector.hpp"
 #include "bilinearform.hpp"
+#include "complex_fem.hpp"
 #ifdef MFEM_USE_MPI
 #include "pgridfunc.hpp"
 #endif
@@ -519,6 +520,200 @@ public:
 };
 
 
+/** @brief The ComplexLpErrorEstimator class compares the solution to a known
+    coefficient.
+
+    This class can be used, for example, to adapt a mesh to a non-trivial
+    initial condition in a time-dependent simulation. It can also be used to
+    force refinement in the neighborhood of small features before switching to a
+    more traditional error estimator.
+
+    The ComplexLpErrorEstimator supports either complex-valued scalar or vector\    coefficients and works both in serial and in parallel.
+*/
+class ComplexLpErrorEstimator : public ErrorEstimator
+{
+protected:
+   long current_sequence;
+   int local_norm_p;
+   Vector error_estimates;
+
+   // double total_error = 0.0;
+
+   Coefficient * real_coef;
+   Coefficient * imag_coef;
+   VectorCoefficient * real_vcoef;
+   VectorCoefficient * imag_vcoef;
+   ComplexGridFunction * sol;
+#ifdef MFEM_USE_MPI
+   ParComplexGridFunction * par_sol;
+#endif
+
+   LpErrorEstimator real_estimator;
+   LpErrorEstimator imag_estimator;
+
+   /// Check if the mesh of the solution was modified.
+   bool MeshIsModified()
+   {
+      long mesh_sequence = 0;
+      if (sol) { mesh_sequence = sol->FESpace()->GetMesh()->GetSequence(); }
+#ifdef MFEM_USE_MPI
+      if (par_sol)
+      { mesh_sequence = par_sol->FESpace()->GetMesh()->GetSequence(); }
+#endif
+      MFEM_ASSERT(mesh_sequence >= current_sequence, "");
+      return (mesh_sequence > current_sequence);
+   }
+
+   /// Compute the element error estimates.
+   void ComputeEstimates();
+
+public:
+   /** @brief Construct a new ComplexLpErrorEstimator object for a scalar field.
+       @param p    Integer which selects which Lp norm to use.
+       @param sol  The ComplexGridFunction representation of the scalar field.
+       Note: the coefficient must be set before use with the SetCoef method.
+   */
+   ComplexLpErrorEstimator(int p, ComplexGridFunction &sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(NULL), imag_coef(NULL),
+        real_vcoef(NULL), imag_vcoef(NULL), sol(&sol),
+#ifdef MFEM_USE_MPI
+        par_sol(NULL),
+#endif
+        real_estimator(p, sol.real()), imag_estimator(p, sol.imag()) { }
+
+   /** @brief Construct a new ComplexLpErrorEstimator object for a scalar field.
+       @param p    Integer which selects which Lp norm to use.
+       @param real_coef The scalar Coefficient to compare to the real part of
+                        the solution.
+       @param imag_coef The scalar Coefficient to compare to the imaginary part
+                        of the solution.
+       @param sol  The ComplexGridFunction representation of the scalar field.
+   */
+   ComplexLpErrorEstimator(int p,
+                           Coefficient &real_coef, Coefficient &imag_coef,
+                           ComplexGridFunction &sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(&real_coef), imag_coef(&imag_coef),
+        real_vcoef(NULL), imag_vcoef(NULL), sol(&sol),
+#ifdef MFEM_USE_MPI
+        par_sol(NULL),
+#endif
+        real_estimator(p, real_coef, sol.real()),
+        imag_estimator(p, imag_coef, sol.imag()) { }
+
+   /** @brief Construct a new ComplexLpErrorEstimator object for a vector field.
+       @param p    Integer which selects which Lp norm to use.
+       @param real_coef The vector VectorCoefficient to compare to the real
+                        part of the solution.
+       @param imag_coef The vector VectorCoefficient to compare to the
+                        imaginary part of the solution.
+       @param sol  The ComplexGridFunction representation of the vector field.
+   */
+   ComplexLpErrorEstimator(int p,
+                           VectorCoefficient &real_coef,
+                           VectorCoefficient &imag_coef,
+                           ComplexGridFunction &sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(NULL), imag_coef(NULL),
+        real_vcoef(&real_coef), imag_vcoef(&imag_coef), sol(&sol),
+#ifdef MFEM_USE_MPI
+        par_sol(NULL),
+#endif
+        real_estimator(p, real_coef, sol.real()),
+        imag_estimator(p, imag_coef, sol.imag()) { }
+
+#ifdef MFEM_USE_MPI
+   /** @brief Construct a new ComplexLpErrorEstimator object for a scalar field.
+       @param p    Integer which selects which Lp norm to use.
+       @param sol  The ComplexGridFunction representation of the scalar field.
+       Note: the coefficient must be set before use with the SetCoef method.
+   */
+   ComplexLpErrorEstimator(int p, ParComplexGridFunction &par_sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(NULL), imag_coef(NULL),
+        real_vcoef(NULL), imag_vcoef(NULL),
+        sol(NULL), par_sol(&par_sol),
+        real_estimator(p, par_sol.real()), imag_estimator(p, par_sol.imag()) { }
+
+   /** @brief Construct a new ComplexLpErrorEstimator object for a scalar field.
+       @param p    Integer which selects which Lp norm to use.
+       @param real_coef The scalar Coefficient to compare to the real part of
+                        the solution.
+       @param imag_coef The scalar Coefficient to compare to the imaginary part
+                        of the solution.
+       @param sol  The ComplexGridFunction representation of the scalar field.
+   */
+   ComplexLpErrorEstimator(int p,
+                           Coefficient &real_coef, Coefficient &imag_coef,
+                           ParComplexGridFunction &par_sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(&real_coef), imag_coef(&imag_coef),
+        real_vcoef(NULL), imag_vcoef(NULL),
+        sol(NULL), par_sol(&par_sol),
+        real_estimator(p, real_coef, par_sol.real()),
+        imag_estimator(p, imag_coef, par_sol.imag()) { }
+
+   /** @brief Construct a new ComplexLpErrorEstimator object for a vector field.
+       @param p    Integer which selects which Lp norm to use.
+       @param real_coef The vector VectorCoefficient to compare to the real
+                        part of the solution.
+       @param imag_coef The vector VectorCoefficient to compare to the
+                        imaginary part of the solution.
+       @param sol  The ComplexGridFunction representation of the vector field.
+   */
+   ComplexLpErrorEstimator(int p,
+                           VectorCoefficient &real_coef,
+                           VectorCoefficient &imag_coef,
+                           ParComplexGridFunction &par_sol)
+      : current_sequence(-1), local_norm_p(p),
+        error_estimates(0),
+        real_coef(NULL), imag_coef(NULL),
+        real_vcoef(&real_coef), imag_vcoef(&imag_coef),
+        sol(NULL), par_sol(&par_sol),
+        real_estimator(p, real_coef, par_sol.real()),
+        imag_estimator(p, imag_coef, par_sol.imag()) { }
+#endif
+
+   /** @brief Set the exponent, p, of the Lp norm used for computing the local
+       element errors. */
+   void SetLocalErrorNormP(int p)
+   {
+      local_norm_p = p;
+      real_estimator.SetLocalErrorNormP(p);
+      imag_estimator.SetLocalErrorNormP(p);
+   }
+
+   void SetRealCoef(Coefficient &A)
+   { real_coef = &A; real_estimator.SetCoef(A); }
+   void SetImagCoef(Coefficient &A)
+   { imag_coef = &A; imag_estimator.SetCoef(A); }
+   void SetRealCoef(VectorCoefficient &A)
+   { real_vcoef = &A; real_estimator.SetCoef(A); }
+   void SetImagCoef(VectorCoefficient &A)
+   { imag_vcoef = &A; imag_estimator.SetCoef(A); }
+
+   /// Reset the error estimator.
+   virtual void Reset() override
+   { current_sequence = -1; real_estimator.Reset(); imag_estimator.Reset(); }
+
+   /// Get a Vector with all element errors.
+   virtual const Vector &GetLocalErrors() override
+   {
+      if (MeshIsModified()) { ComputeEstimates(); }
+      return error_estimates;
+   }
+
+   /// Destructor
+   virtual ~ComplexLpErrorEstimator() {}
+};
+
+
 /** @brief The KellyErrorEstimator class provides a fast error indication
     strategy for smooth scalar parallel problems.
 

fem/gridfunc.cpp

@@ -2738,7 +2738,6 @@ void GridFunction::ProjectBdrCoefficientNormal(
    Array<int> dofs;
    int dim = vcoeff.GetVDim();
    Vector vc(dim), nor(dim), lvec;
-
    for (int i = 0; i < fes->GetNBE(); i++)
    {
       if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)

linalg/CMakeLists.txt

@@ -27,6 +27,7 @@ list(APPEND SRCS
   sparsemat.cpp
   sparsesmoothers.cpp
   vector.cpp
+  vector_operator.cpp
   )
 
 list(APPEND HDRS
@@ -57,6 +58,7 @@ list(APPEND HDRS
   tensor.hpp
   dual.hpp
   vector.hpp
+  vector_operator.hpp
   )
 
 if (MFEM_USE_MPI)

linalg/blockoperator.cpp

@@ -361,4 +361,133 @@ BlockLowerTriangularPreconditioner::~BlockLowerTriangularPreconditioner()
    }
 }
 
+SchurComplimentOperator::SchurComplimentOperator(Solver & AInv, Operator * B,
+                                                 Operator * C, Operator & D)
+   : Operator(),
+     APtr(NULL), BPtr(B), CPtr(C), DPtr(&D), AInvPtr(&AInv), DInvPtr(NULL),
+     sizeA(AInv.Height()), sizeD(D.Height())
+{
+   height = sizeD;
+   width  = height;
+
+   rhs.SetSize(sizeD);
+
+   y2.SetSize(sizeD);
+   x1.SetSize(sizeA);
+   rhs1.SetSize(sizeA);
+}
+
+SchurComplimentOperator::SchurComplimentOperator(Operator & A, Operator * B,
+                                                 Operator * C, Solver & DInv)
+   : APtr(&A), BPtr(B), CPtr(C), DPtr(NULL), AInvPtr(NULL), DInvPtr(&DInv),
+     sizeA(A.Height()), sizeD(DInv.Height())
+{
+   height = sizeA;
+   width  = height;
+
+   rhs.SetSize(sizeA);
+
+   y1.SetSize(sizeA);
+   x2.SetSize(sizeD);
+   rhs2.SetSize(sizeD);
+}
+
+const Vector & SchurComplimentOperator::GetRHSVector(const Vector & a,
+                                                     const Vector & b)
+{
+   if (DInvPtr)
+   {
+      if (BPtr)
+      {
+         DInvPtr->Mult(b, x2);
+         BPtr->Mult(x2, rhs);
+         rhs *= -1.0;
+         rhs.Add(1.0, a);
+      }
+      else
+      {
+         rhs.Set(1.0, a);
+      }
+   }
+   else
+   {
+      if (CPtr)
+      {
+         AInvPtr->Mult(a, x1);
+         CPtr->Mult(x1, rhs);
+         rhs *= -1.0;
+         rhs.Add(1.0, b);
+      }
+      else
+      {
+         rhs.Set(1.0, b);
+      }
+   }
+
+   return rhs;
+}
+
+void SchurComplimentOperator::Mult(const Vector & x, Vector & y) const
+{
+   if (DInvPtr)
+   {
+      APtr->Mult(x, y);
+
+      if (BPtr && CPtr)
+      {
+         CPtr->Mult(x, rhs2);
+         DInvPtr->Mult(rhs2, x2);
+         BPtr->Mult(x2, y1);
+
+         y.Add(-1.0, y1);
+      }
+   }
+   else
+   {
+      DPtr->Mult(x, y);
+
+      if (BPtr && CPtr)
+      {
+         BPtr->Mult(x, rhs1);
+         AInvPtr->Mult(rhs1, x1);
+         CPtr->Mult(x1, y2);
+
+         y.Add(-1.0, y2);
+      }
+   }
+}
+
+void SchurComplimentOperator::Solve(const Vector & b, const Vector & x,
+                                    Vector & y)
+{
+   if (DInvPtr)
+   {
+      if (CPtr)
+      {
+         CPtr->Mult(x, rhs2);
+         rhs2 *= -1.0;
+         rhs2.Add(1.0, b);
+      }
+      else
+      {
+         rhs2.Set(1.0, b);
+      }
+      DInvPtr->Mult(rhs2, y);
+   }
+   else
+   {
+      if (BPtr)
+      {
+         BPtr->Mult(x, rhs1);
+         rhs1 *= -1.0;
+         rhs1.Add(1.0, b);
+      }
+      else
+      {
+         rhs1.Set(1.0, b);
+      }
+      AInvPtr->Mult(rhs1, y);
+   }
+}
+
 }

linalg/blockoperator.hpp

@@ -290,6 +290,42 @@ private:
    mutable Vector tmp2;
 };
 
+class SchurComplimentOperator : public Operator
+{
+private:
+   Operator * APtr;
+   Operator * BPtr;
+   Operator * CPtr;
+   Operator * DPtr;
+
+   Solver * AInvPtr;
+   Solver * DInvPtr;
+
+   int sizeA, sizeD;
+
+   Vector rhs;
+
+   mutable Vector x1;
+   mutable Vector x2;
+   mutable Vector y1;
+   mutable Vector y2;
+   mutable Vector rhs1;
+   mutable Vector rhs2;
+
+public:
+   SchurComplimentOperator(Solver & AInv, Operator * B,
+                           Operator * C, Operator & D);
+
+   SchurComplimentOperator(Operator & A, Operator * B,
+                           Operator * C, Solver & DInv);
+
+   const Vector & GetRHSVector(const Vector & a, const Vector & b);
+
+   void Mult(const Vector & x, Vector & y) const;
+
+   void Solve(const Vector & b, const Vector & x, Vector & y);
+};
+
 }
 
 #endif /* MFEM_BLOCKOPERATOR */

linalg/complex_operator.cpp

@@ -877,6 +877,434 @@ ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
    delete [] stat;
 }
 
+#ifdef MFEM_USE_MUMPS
+
+void ComplexMUMPSSolver::SetOperator(const Operator &op)
+{
+   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
+
+   MFEM_VERIFY(APtr, "Not compatible matrix type");
+   height = op.Height();
+   width = op.Width();
+
+   conv = APtr->GetConvention();
+   comm = APtr->real().GetComm();
+   MPI_Comm_size(comm, &numProcs);
+   MPI_Comm_rank(comm, &myid);
+
+   auto parcsr_op_r = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>
+                      (APtr->real());
+   auto parcsr_op_i = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>
+                      (APtr->imag());
+
+   hypre_CSRMatrix *csr_op_r = hypre_MergeDiagAndOffd(parcsr_op_r);
+   hypre_CSRMatrix *csr_op_i = hypre_MergeDiagAndOffd(parcsr_op_i);
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op_r);
+   hypre_CSRMatrixBigJtoJ(csr_op_i);
+#endif
+   MFEM_VERIFY(csr_op_r->num_nonzeros == csr_op_i->num_nonzeros,
+               "Incompatible sparsity partters");
+
+   int *Iptr = csr_op_r->i;
+   int *Jptr = csr_op_r->j;
+   int n_loc = csr_op_r->num_rows;
+   row_start = parcsr_op_i->first_row_index;
+   MUMPS_INT8 nnz = csr_op_r->num_nonzeros;
+
+   int * I = new int[nnz];
+   int * J = new int[nnz];
+
+   // Fill in I and J arrays for
+   // COO format in 1-based indexing
+   int k = 0;
+   double * data_r = csr_op_r->data;
+   double * data_i = csr_op_i->data;
+   mumps_double_complex *zdata = new mumps_double_complex[nnz];
+   for (int i = 0; i < n_loc; i++)
+   {
+      for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+      {
+         I[k] = row_start + i + 1;
+         J[k] = Jptr[k] + 1;
+         zdata[k].r = data_r[k];
+         zdata[k].i = data_i[k];
+         k++;
+      }
+   }
+
+   // new MUMPS object
+   if (id)
+   {
+      id->job = -2;
+      zmumps_c(id);
+      delete id;
+   }
+   id = new ZMUMPS_STRUC_C;
+   // C to Fortran communicator
+   id->comm_fortran = (MUMPS_INT) MPI_Comm_c2f(comm);
+
+   // Host is involved in computation
+   id->par = 1;
+
+   id->sym = 0;
+
+   // MUMPS init
+   id->job = -1;
+   zmumps_c(id);
+
+   // Set MUMPS default parameters
+   SetParameters();
+
+   id->n = parcsr_op_r->global_num_rows;
+
+   id->nnz_loc = nnz;
+
+   id->irn_loc = I;
+
+   id->jcn_loc = J;
+
+   id->a_loc = zdata;
+
+   // MUMPS Analysis
+   id->job = 1;
+   zmumps_c(id);
+
+   // MUMPS Factorization
+   id->job = 2;
+   zmumps_c(id);
+
+   hypre_CSRMatrixDestroy(csr_op_r);
+   hypre_CSRMatrixDestroy(csr_op_i);
+   delete [] I;
+   delete [] J;
+   delete [] zdata;
+
+#if MFEM_MUMPS_VERSION >= 530
+   delete [] irhs_loc;
+   irhs_loc = new int[n_loc];
+   for (int i = 0; i < n_loc; i++)
+   {
+      irhs_loc[i] = row_start + i + 1;
+   }
+   row_starts.SetSize(numProcs);
+   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
+#else
+   if (myid == 0)
+   {
+      delete [] rhs_glob;
+      delete [] recv_counts;
+      global_num_rows = parcsr_op_r->global_num_rows;
+      rhs_glob = new mumps_double_complex[global_num_rows];
+      recv_counts = new int[numProcs];
+   }
+   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
+   if (myid == 0)
+   {
+      delete [] displs;
+      displs = new int[numProcs];
+      displs[0] = 0;
+      int s = 0;
+      for (int k = 0; k < numProcs-1; k++)
+      {
+         s += recv_counts[k];
+         displs[k+1] = s;
+      }
+   }
+#endif
+}
+
+void ComplexMUMPSSolver::Mult(const Vector &x, Vector &y) const
+{
+   int n = x.Size()/2;
+   double * datax = x.GetData();
+   double * datay = y.GetData();
+   Vector ximag;
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      ximag.SetDataAndSize(&datax[n],n);
+      ximag *=-1.0;
+   }
+#if MFEM_MUMPS_VERSION >= 530
+   id->nloc_rhs = n;
+   id->lrhs_loc = n;
+   mumps_double_complex *zx = new mumps_double_complex[n];
+   for (int i = 0; i<n; i++)
+   {
+      zx[i].r = x[i];
+      zx[i].i = x[n+i];
+   }
+   id->rhs_loc = zx;
+   id->irhs_loc = irhs_loc;
+
+   id->lsol_loc = id->MUMPSC_INFO(23);
+   id->isol_loc = new int[id->MUMPSC_INFO(23)];
+   id->sol_loc = new mumps_double_complex[id->MUMPSC_INFO(23)];
+
+   // MUMPS solve
+   id->job = 3;
+   zmumps_c(id);
+
+   double *zy = new double[2*id->MUMPSC_INFO(23)];
+   for (int i = 0; i<id->MUMPSC_INFO(23); i++)
+   {
+      zy[i] = id->sol_loc[i].r;
+      zy[id->MUMPSC_INFO(23)+i] = id->sol_loc[i].i;
+   }
+
+   RedistributeSol(id->isol_loc, zy, y.GetData());
+
+   delete [] zy;
+   delete [] zx;
+   delete [] id->sol_loc;
+   delete [] id->isol_loc;
+#else
+   // real
+   double * rhs_glob_r = nullptr;
+   double * rhs_glob_i = nullptr;
+   if (myid == 0)
+   {
+      rhs_glob_r = new double[global_num_rows];
+      rhs_glob_i = new double[global_num_rows];
+   }
+   double * xdata = x.GetData();
+   MPI_Gatherv(xdata, n, MPI_DOUBLE,
+               rhs_glob_r, recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+   MPI_Gatherv(&xdata[n], n, MPI_DOUBLE,
+               rhs_glob_i, recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+
+   if (myid == 0)
+   {
+      for (int i = 0; i<global_num_rows; i++)
+      {
+         rhs_glob[i].r = rhs_glob_r[i];
+         rhs_glob[i].i = rhs_glob_i[i];
+      }
+      id->rhs = rhs_glob;
+   }
+   // MUMPS solve
+   id->job = 3;
+   zmumps_c(id);
+   if (myid == 0)
+   {
+      for (int i = 0; i<global_num_rows; i++)
+      {
+         rhs_glob_r[i] = rhs_glob[i].r;
+         rhs_glob_i[i] = rhs_glob[i].i;
+      }
+   }
+   double * ydata = y.GetData();
+   MPI_Scatterv(rhs_glob_r, recv_counts, displs,
+                MPI_DOUBLE, ydata, n,
+                MPI_DOUBLE, 0, comm);
+   MPI_Scatterv(rhs_glob_i, recv_counts, displs,
+                MPI_DOUBLE, &ydata[n], n,
+                MPI_DOUBLE, 0, comm);
+
+   if (myid == 0)
+   {
+      delete [] rhs_glob_r;
+      delete [] rhs_glob_i;
+   }
+
+#endif
+
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      ximag *=-1.0;
+   }
+}
+
+void ComplexMUMPSSolver::SetPrintLevel(int print_lvl)
+{
+   print_level = print_lvl;
+}
+
+ComplexMUMPSSolver::~ComplexMUMPSSolver()
+{
+   if (id)
+   {
+#if MFEM_MUMPS_VERSION >= 530
+      delete [] irhs_loc;
+#else
+      delete [] recv_counts;
+      delete [] displs;
+      delete [] rhs_glob;
+#endif
+      id->job = -2;
+      zmumps_c(id);
+      delete id;
+   }
+}
+
+void ComplexMUMPSSolver::SetParameters()
+{
+   // output stream for error messages
+   id->ICNTL(1) = 6;
+   // output stream for diagnosting printing local to each proc
+   id->ICNTL(2) = 6;
+   // output stream for global info
+   id->ICNTL(3) = 6;
+   // Level of error printing
+   id->ICNTL(4) = print_level;
+   //input matrix format (assembled)
+   id->ICNTL(5) = 0;
+   // Use A or A^T
+   id->ICNTL(9) = 1;
+   // Iterative refinement (disabled)
+   id->ICNTL(10) = 0;
+   // Error analysis-statistics (disabled)
+   id->ICNTL(11) = 0;
+   // Use of ScaLAPACK (Parallel factorization on root)
+   id->ICNTL(13) = 0;
+   // Percentage increase of estimated workspace (default = 20%)
+   id->ICNTL(14) = 20;
+   // Number of OpenMP threads (default)
+   id->ICNTL(16) = 0;
+   // Matrix input format (distributed)
+   id->ICNTL(18) = 3;
+   // Schur complement (no Schur complement matrix returned)
+   id->ICNTL(19) = 0;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS
+   id->ICNTL(20) = 10;
+   // Distributed Sol
+   id->ICNTL(21) = 1;
+#else
+   // Centralized RHS
+   id->ICNTL(20) = 0;
+   // Centralized Sol
+   id->ICNTL(21) = 0;
+#endif
+   // Out of core factorization and solve (disabled)
+   id->ICNTL(22) = 0;
+   // Max size of working memory (default = based on estimates)
+   id->ICNTL(23) = 0;
+}
+
+#if MFEM_MUMPS_VERSION >= 530
+int ComplexMUMPSSolver::GetRowRank(int i, const Array<int> &row_starts_) const
+{
+   if (row_starts_.Size() == 1)
+   {
+      return 0;
+   }
+   auto up = std::upper_bound(row_starts_.begin(), row_starts_.end(), i);
+   return std::distance(row_starts_.begin(), up) - 1;
+}
+
+void ComplexMUMPSSolver::RedistributeSol(const int * row_map,
+                                         const double * x, double * y) const
+{
+   int size = id->MUMPSC_INFO(23);
+   int n = id->nloc_rhs;
+   int * send_count = new int[numProcs]();
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank) { continue; }
+      send_count[row_rank]++;
+   }
+
+   int * recv_count = new int[numProcs];
+   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
+
+   int * send_displ = new int [numProcs]; send_displ[0] = 0;
+   int * recv_displ = new int [numProcs]; recv_displ[0] = 0;
+   int sbuff_size = send_count[numProcs-1];
+   int rbuff_size = recv_count[numProcs-1];
+   for (int k = 0; k < numProcs - 1; k++)
+   {
+      send_displ[k + 1] = send_displ[k] + send_count[k];
+      recv_displ[k + 1] = recv_displ[k] + recv_count[k];
+      sbuff_size += send_count[k];
+      rbuff_size += recv_count[k];
+   }
+
+   int * sendbuf_index = new int[sbuff_size];
+   double * sendbuf_values_r = new double[sbuff_size];
+   double * sendbuf_values_i = new double[sbuff_size];
+   int * soffs = new int[numProcs]();
+
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank)
+      {
+         int local_index = j - row_start;
+         y[local_index] = x[i];
+         y[local_index+n] = x[i+size];
+      }
+      else
+      {
+         int k = send_displ[row_rank] + soffs[row_rank];
+         sendbuf_index[k] = j;
+         sendbuf_values_r[k] = x[i];
+         sendbuf_values_i[k] = x[i+size];
+         soffs[row_rank]++;
+      }
+   }
+
+   int * recvbuf_index = new int[rbuff_size];
+   double * recvbuf_values_r = new double[rbuff_size];
+   double * recvbuf_values_i = new double[rbuff_size];
+   MPI_Alltoallv(sendbuf_index,
+                 send_count,
+                 send_displ,
+                 MPI_INT,
+                 recvbuf_index,
+                 recv_count,
+                 recv_displ,
+                 MPI_INT,
+                 comm);
+   MPI_Alltoallv(sendbuf_values_r,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_values_r,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+   MPI_Alltoallv(sendbuf_values_i,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_values_i,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+
+   // Unpack recv buffer
+   for (int i = 0; i < rbuff_size; i++)
+   {
+      int local_index = recvbuf_index[i] - row_start;
+      y[local_index] = recvbuf_values_r[i];
+      y[local_index+n] = recvbuf_values_i[i];
+   }
+
+   delete [] recvbuf_values_r;
+   delete [] recvbuf_values_i;
+   delete [] recvbuf_index;
+   delete [] soffs;
+   delete [] sendbuf_values_r;
+   delete [] sendbuf_values_i;
+   delete [] sendbuf_index;
+   delete [] recv_displ;
+   delete [] send_displ;
+   delete [] recv_count;
+   delete [] send_count;
+}
+#endif // MUMPS VERSION
+#endif // MFEM_USE_CMUMPS
+
+
 #endif // MFEM_USE_MPI
 
 }

linalg/complex_operator.hpp

@@ -22,6 +22,11 @@
 #include <umfpack.h>
 #endif
 
+#ifdef MFEM_USE_MUMPS
+#include "zmumps_c.h"
+#include <vector>
+#endif
+
 namespace mfem
 {
 
@@ -289,6 +294,49 @@ private:
    int nranks_;
 };
 
+#ifdef MFEM_USE_MUMPS
+
+class ComplexMUMPSSolver : public mfem::Solver
+{
+public:
+   ComplexMUMPSSolver() {}
+   void SetOperator(const Operator &op);
+   void Mult(const Vector &x, Vector &y) const;
+   void SetPrintLevel(int print_lvl);
+   ~ComplexMUMPSSolver();
+private:
+   MPI_Comm comm;
+   ComplexOperator::Convention conv;
+   int numProcs;
+   int myid;
+   int print_level = 0;
+   int row_start;
+#define ICNTL(I) icntl[(I) -1]
+#define MUMPSC_INFO(I) info[(I) -1]
+   ZMUMPS_STRUC_C *id=nullptr;
+   void SetParameters();
+
+#if MFEM_MUMPS_VERSION >= 530
+
+   Array<int> row_starts;
+   int * irhs_loc = nullptr;
+
+   int GetRowRank(int i, const Array<int> &row_starts_) const;
+
+   void RedistributeSol(const int * row_map,
+                        const double * x,
+                        double * y) const;
+#else
+   int global_num_rows;
+   int * recv_counts = nullptr;
+   int * displs = nullptr;
+   mumps_double_complex * rhs_glob = nullptr;
+#endif
+
+}; // mfem::ComplexMUMPSSolver class
+
+#endif // MFEM_USE_CMUMPS
+
 #endif // MFEM_USE_MPI
 
 }

linalg/mumps.cpp

@@ -446,6 +446,10 @@ void MUMPSSolver::SetParameters()
    id->MUMPS_ICNTL(4) = print_level;
    // Input matrix format (assembled)
    id->MUMPS_ICNTL(5) = 0;
+   // Overiding default reordering to PARMETIS
+   id->MUMPS_ICNTL(7) = 5;
+   id->MUMPS_ICNTL(28) = 2;
+   id->MUMPS_ICNTL(29) = 2;
    // Use A or A^T
    id->MUMPS_ICNTL(9) = 1;
    // Iterative refinement (disabled)

linalg/vector_operator.cpp

@@ -0,0 +1,80 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "vector_operator.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_MPI
+
+ParVectorOperator::ParVectorOperator(MPI_Comm comm,
+                                     int myid,
+                                     int local_vec_size,
+                                     int num_vecs)
+   : Operator((myid == 0) ? num_vecs : 0, local_vec_size),
+     comm(comm),
+     myid(myid),
+     vecs(num_vecs),
+     coefs(num_vecs),
+     owns(num_vecs)
+{
+   vecs = NULL;
+   coefs = 1.0;
+   owns = false;
+}
+
+ParVectorOperator::~ParVectorOperator()
+{
+   for (int i=0; i < vecs.Size(); i++)
+   {
+      if (owns[i]) { delete vecs[i]; }
+      vecs[i] = NULL;
+   }
+}
+
+void ParVectorOperator::SetVector(int idx, Vector *vec,
+                                  double c, bool own_vec)
+{
+   MFEM_VERIFY(idx >= 0 && idx < vecs.Size(),
+               "ParVectorOperator: Index out of range");
+
+   vecs[idx] = vec;
+   coefs[idx] = c;
+   owns[idx] = own_vec;
+}
+
+void ParVectorOperator::Mult(const Vector &x, Vector &y) const
+{
+   for (int i=0; i<vecs.Size(); i++)
+   {
+      double vo = coefs[i] * (*vecs[i] * x);
+      double vi = 0.0;
+      MPI_Reduce(&vo, &vi, 1, MPI_DOUBLE, MPI_SUM, 0, comm);
+      if (myid == 0) { y[i] = vi; }
+   }
+}
+
+/// Action of the transpose operator: `y=A^t(x)`.
+void ParVectorOperator::MultTranspose(const Vector &x, Vector &y) const
+{
+   y = 0.0;
+   for (int i=0; i<vecs.Size(); i++)
+   {
+      double xi = (myid == 0) ? x[i] : 0.0;
+      MPI_Bcast(&xi, 1, MPI_DOUBLE, 0, comm);
+      y.Add(xi * coefs[i], *vecs[i]);
+   }
+}
+
+}
+
+#endif // MFEM_USE_MPI

linalg/vector_operator.hpp

@@ -0,0 +1,65 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_VECTOR_OPERATOR
+#define MFEM_VECTOR_OPERATOR
+
+#include "operator.hpp"
+#include "densemat.hpp"
+#include "vector.hpp"
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_MPI
+
+class ParVectorOperator : public Operator
+{
+private:
+   MPI_Comm comm;
+   int myid;
+
+   Array<Vector*> vecs;
+   Array<double> coefs;
+   Array<bool> owns;
+
+public:
+   ParVectorOperator(MPI_Comm comm,
+                     int myid,
+                     int local_vec_size,
+                     int num_vecs);
+
+   ~ParVectorOperator();
+
+   void SetVector(int idx, Vector *vec,
+                  double c = 1.0, bool own_vec = false);
+
+   /// Operator application: `y=A(x)`.
+   void Mult(const Vector &x, Vector &y) const;
+
+   /// Action of the transpose operator: `y=A^t(x)`.
+   void MultTranspose(const Vector &x, Vector &y) const;
+
+   /// Compute LQ factorization of this operator
+   /** This operator, A, represents a matrix with very few rows but
+       many columns which are distributed across multiple
+       processors. The LQ factorization LQ = A is related to the QR
+       factorization of the transpose of A with L = R^T and the two Q
+       operators being transposes of eachother.
+   */
+   void GetLQFactors(DenseMatrix &L, ParVectorOperator &Q);
+};
+
+#endif // MFEM_USE_MPI
+
+} // namespace mfem
+
+#endif // MFEM_VECTOR_OPERATOR

makefile

@@ -125,7 +125,7 @@ EXAMPLE_TEST_DIRS := examples
 
 MINIAPP_SUBDIRS = common electromagnetics meshing navier performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag autodiff hooke \
- multidomain dpg hdiv-linear-solver spde
+ multidomain dpg hdiv-linear-solver spde plasma
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \

miniapps/CMakeLists.txt

@@ -22,6 +22,7 @@ add_subdirectory(electromagnetics)
 add_subdirectory(navier)
 add_subdirectory(meshing)
 add_subdirectory(performance)
+add_subdirectory(plasma)
 add_subdirectory(tools)
 add_subdirectory(toys)
 add_subdirectory(nurbs)

miniapps/common/mesh_extras.cpp

@@ -10,6 +10,8 @@
 // CONTRIBUTING.md for details.
 
 #include "mesh_extras.hpp"
+#include <set>
+#include <map>
 
 using namespace std;
 
@@ -213,6 +215,301 @@ MergeMeshNodes(Mesh * mesh, int logging)
    }
 }
 
+void
+IdentifyPeriodicMeshVertices(const Mesh & mesh,
+                             const vector<Vector> & trans_vecs,
+                             Array<int> & v2v,
+                             int logging)
+{
+   int sdim = mesh.SpaceDimension();
+
+   double tol = 1.0e-8;
+   double dia = -1.0;
+
+   // map<int,map<int,map<int,int> > > c2v;
+   set<int> v;
+   set<int>::iterator si, sj, sk;
+   map<int,int>::iterator mi;
+   map<int,set<int> >::iterator msi;
+
+   Vector coord(NULL, sdim);
+
+   // map<int,vector<double> > bnd_vtx;
+   // map<int,vector<double> > shft_bnd_vtx;
+
+   // int d = 5;
+   Vector xMax(sdim), xMin(sdim), xDiff(sdim);
+   xMax = xMin = xDiff = 0.0;
+
+   for (int be=0; be<mesh.GetNBE(); be++)
+   {
+      Array<int> dofs;
+      mesh.GetBdrElementVertices(be,dofs);
+
+      for (int i=0; i<dofs.Size(); i++)
+      {
+         v.insert(dofs[i]);
+
+         coord.SetData(const_cast<double*>(mesh.GetVertex(dofs[i])));
+         for (int j=0; j<sdim; j++)
+         {
+            xMax[j] = max(xMax[j],coord[j]);
+            xMin[j] = min(xMin[j],coord[j]);
+         }
+      }
+   }
+   add(xMax, -1.0, xMin, xDiff);
+   dia = xDiff.Norml2();
+
+   if ( logging > 0 )
+   {
+      cout << "Number of Boundary Vertices:  " << v.size() << endl;
+
+      cout << "xMin: ";
+      xMin.Print(cout,sdim);
+      cout << "xMax: ";
+      xMax.Print(cout,sdim);
+      cout << "xDiff: ";
+      xDiff.Print(cout,sdim);
+   }
+
+   if ( logging > 0 )
+   {
+      for (si=v.begin(); si!=v.end(); si++)
+      {
+         cout << *si << ": ";
+         coord.SetData(const_cast<double*>(mesh.GetVertex(*si)));
+         coord.Print(cout);
+      }
+   }
+
+   map<int,int>        slaves;
+   map<int,set<int> > masters;
+
+   for (si=v.begin(); si!=v.end(); si++) { masters[*si]; }
+
+   Vector at(sdim);
+   Vector dx(sdim);
+
+   for (unsigned int i=0; i<trans_vecs.size(); i++)
+   {
+      int c = 0;
+      if ( logging > 0 )
+      {
+         cout << "trans_vecs = ";
+         trans_vecs[i].Print(cout,sdim);
+      }
+
+      for (si=v.begin(); si!=v.end(); si++)
+      {
+         coord.SetData(const_cast<double*>(mesh.GetVertex(*si)));
+
+         add(coord, trans_vecs[i], at);
+
+         for (sj=v.begin(); sj!=v.end(); sj++)
+         {
+            coord.SetData(const_cast<double*>(mesh.GetVertex(*sj)));
+            add(at, -1.0, coord, dx);
+
+            if ( dx.Norml2() > dia * tol )
+            {
+               continue;
+            }
+
+            int master = *si;
+            int slave  = *sj;
+
+            bool mInM = masters.find(master) != masters.end();
+            bool sInM = masters.find(slave)  != masters.end();
+
+            if ( mInM && sInM )
+            {
+               // Both vertices are currently masters
+               //   Demote "slave" to be a slave of master
+               if ( logging > 0 )
+               {
+                  cout << "Both " << master << " and " << slave
+                       << " are masters." << endl;
+               }
+               masters[master].insert(slave);
+               slaves[slave] = master;
+               for (sk=masters[slave].begin();
+                    sk!=masters[slave].end(); sk++)
+               {
+                  masters[master].insert(*sk);
+                  slaves[*sk] = master;
+               }
+               masters.erase(slave);
+            }
+            else if ( mInM && !sInM )
+            {
+               // "master" is already a master and "slave" is already a slave
+               // Make "master" and its slaves slaves of "slave"'s master
+               if ( logging > 0 )
+               {
+                  cout << master << " is already a master and " << slave
+                       << " is already a slave of " << slaves[slave]
+                       << "." << endl;
+               }
+               if ( master != slaves[slave] )
+               {
+                  masters[slaves[slave]].insert(master);
+                  slaves[master] = slaves[slave];
+                  for (sk=masters[master].begin();
+                       sk!=masters[master].end(); sk++)
+                  {
+                     masters[slaves[slave]].insert(*sk);
+                     slaves[*sk] = slaves[slave];
+                  }
+                  masters.erase(master);
+               }
+            }
+            else if ( !mInM && sInM )
+            {
+               // "master" is currently a slave and
+               // "slave" is currently a master
+               // Make "slave" and its slaves slaves of "master"'s master
+               if ( logging > 0 )
+               {
+                  cout << master << " is currently a slave of "
+                       << slaves[master]<< " and " << slave
+                       << " is currently a master." << endl;
+               }
+               if ( slave != slaves[master] )
+               {
+                  masters[slaves[master]].insert(slave);
+                  slaves[slave] = slaves[master];
+                  for (sk=masters[slave].begin();
+                       sk!=masters[slave].end(); sk++)
+                  {
+                     masters[slaves[master]].insert(*sk);
+                     slaves[*sk] = slaves[master];
+                  }
+                  masters.erase(slave);
+               }
+            }
+            else
+            {
+               // Both vertices are currently slaves
+               // Make "slave" and its fellow slaves slaves
+               // of "master"'s master
+               if ( logging > 0 )
+               {
+                  cout << "Both " << master << " and " << slave
+                       << " are slaves of " << slaves[master] << " and "
+                       << slaves[slave] << " respectively." << endl;
+               }
+
+               int master_of_master = slaves[master];
+               int master_of_slave  = slaves[slave];
+
+               // Move slave and its fellow slaves to master_of_master
+               if ( slaves[master] != slaves[slave] )
+               {
+                  for (sk=masters[master_of_slave].begin();
+                       sk!=masters[master_of_slave].end(); sk++)
+                  {
+                     masters[master_of_master].insert(*sk);
+                     slaves[*sk] = master_of_master;
+                  }
+                  masters.erase(master_of_slave);
+                  slaves[master_of_slave] = master_of_master;
+               }
+            }
+            c++;
+            break;
+         }
+      }
+      if ( logging > 0 )
+      {
+         cout << "Found " << c << " possible node";
+         if ( c != 1 ) { cout << "s"; }
+         cout <<" to project." << endl;
+      }
+   }
+   if ( logging > 0 )
+   {
+      cout << "Number of Master Vertices:  " << masters.size() << endl;
+      cout << "Number of Slave Vertices:   " << slaves.size() << endl;
+      cout << "Master to slave mapping:" << endl;
+      for (msi=masters.begin(); msi!=masters.end(); msi++)
+      {
+         cout << msi->first << " ->";
+         for (si=msi->second.begin(); si!=msi->second.end(); si++)
+         {
+            cout << " " << *si;
+         }
+         cout << endl;
+      }
+      cout << "Slave to master mapping:" << endl;
+      for (mi=slaves.begin(); mi!=slaves.end(); mi++)
+      {
+         cout << mi->first << " <- " << mi->second << endl;
+      }
+   }
+
+   v2v.SetSize(mesh.GetNV());
+
+   for (int i=0; i<v2v.Size(); i++)
+   {
+      v2v[i] = i;
+   }
+
+   for (mi=slaves.begin(); mi!=slaves.end(); mi++)
+   {
+      v2v[mi->first] = mi->second;
+   }
+}
+
+Mesh *
+MakePeriodicMesh(Mesh * mesh, const Array<int> & v2v,
+                 int logging)
+{
+   int dim  = mesh->Dimension();
+
+   if ( logging > 0 )
+      cout << "Euler Number of Initial Mesh:  "
+           << ((dim==3)?mesh->EulerNumber():mesh->EulerNumber2D()) << endl;
+
+   Mesh *per_mesh = new Mesh(*mesh, true);
+
+   per_mesh->SetCurvature(1, true);
+
+   // renumber elements
+   for (int i = 0; i < per_mesh->GetNE(); i++)
+   {
+      Element *el = per_mesh->GetElement(i);
+      int *v = el->GetVertices();
+      int nv = el->GetNVertices();
+      for (int j = 0; j < nv; j++)
+      {
+         v[j] = v2v[v[j]];
+      }
+   }
+   // renumber boundary elements
+   for (int i = 0; i < per_mesh->GetNBE(); i++)
+   {
+      Element *el = per_mesh->GetBdrElement(i);
+      int *v = el->GetVertices();
+      int nv = el->GetNVertices();
+      for (int j = 0; j < nv; j++)
+      {
+         v[j] = v2v[v[j]];
+      }
+   }
+
+   per_mesh->RemoveUnusedVertices();
+   // per_mesh->RemoveInternalBoundaries();
+
+   if ( logging > 0 )
+   {
+      cout << "Euler Number of Final Mesh:    "
+           << ((dim==3)?per_mesh->EulerNumber():per_mesh->EulerNumber2D())
+           << endl;
+   }
+   return per_mesh;
+}
+
 void AttrToMarker(int max_attr, const Array<int> &attrs, Array<int> &marker)
 {
    MFEM_ASSERT(attrs.Max() <= max_attr, "Invalid attribute number present.");

miniapps/common/mesh_extras.hpp

@@ -14,6 +14,7 @@
 
 #include "mfem.hpp"
 #include <sstream>
+#include <vector>
 
 namespace mfem
 {
@@ -28,7 +29,16 @@ public:
 };
 
 /// Merges vertices which lie at the same location
-void MergeMeshNodes(Mesh * mesh, int logging);
+void MergeMeshNodes(Mesh * mesh, int logging = 0);
+
+void
+IdentifyPeriodicMeshVertices(const Mesh & mesh,
+                             const std::vector<Vector> & trans_vecs,
+                             Array<int> & v2v,
+                             int logging);
+Mesh *
+MakePeriodicMesh(Mesh * mesh, const Array<int> & v2v,
+                 int logging = 0);
 
 /// Convert a set of attribute numbers to a marker array
 /** The marker array will be of size max_attr and it will contain only zeroes

miniapps/meshing/mesh-explorer.cpp

@@ -36,6 +36,7 @@
 #include <fstream>
 #include <limits>
 #include <cstdlib>
+#include <math.h>
 
 using namespace mfem;
 using namespace std;
@@ -44,18 +45,88 @@ using namespace std;
 void transformation(const Vector &p, Vector &v)
 {
    // simple shear transformation
-   double s = 0.1;
+   //double s = 0.1;
+   double h_bump = 0.4;
 
    if (p.Size() == 3)
    {
+       /*
       v(0) = p(0) + s*p(1) + s*p(2);
       v(1) = p(1) + s*p(2) + s*p(0);
       v(2) = p(2);
+        */
+       if (p(0) < 5.4){v(0) = p(0);}
+       else
+       {
+        v(0) = 5.4 + ((p(0) - 5.4)/(0.6))*(0.6 - h_bump*exp((-1.0*pow((p(1) - 0.4), 2.0))/pow(0.1, 2.0)));
+       }
+       v(1) = p(1);
+       v(2) = p(2);
    }
    else if (p.Size() == 2)
    {
-      v(0) = p(0) + s*p(1);
-      v(1) = p(1) + s*p(0);
+       /*
+       double mag = pow(pow(0.5,2) + pow(1.5,2),0.5);
+       v(0) = p(0);
+       v(1) = (0.5/mag)*p(0)+(1.5/mag)*p(1);
+       */
+       /*
+      if (p(0) < 0.4){v(0) = p(0);}
+      else
+      {
+      double endx = 0.7;
+       v(0) = ((p(0) - 0.4)/(endx - 0.4))*((endx - 0.4) - (endx - 0.5)) + 0.4;
+      }
+       v(1) = p(1);
+       
+       
+       if (p(0) < 5.4){v(0) = p(0);}
+       else
+       {
+        v(0) = 5.4 + ((p(0) - 5.4)/(0.6))*(0.6 - h_bump*exp((-1.0*pow((p(1) - 0.4), 2.0))/pow(0.1, 2.0)));
+       }
+       v(1) = p(1);
+        */
+
+      // Kohno 2015:
+      /*
+      double r = sqrt(pow(p[0],2.0)+pow(p[1],2.0));
+      double theta_rad = atan2(p[1], p[0]);
+
+      double scale = 1.0 - fabs(r-0.3)/fabs(0.3);
+      double h = 0.02;
+      double decay = 1.2;
+      double theta_on = 160.0;
+      double theta_off = 200.0;
+
+      double theta = theta_rad*(180.0/M_PI);
+      if (theta < 0){theta += 360.0;}
+      double delta_r = h*0.5*((1+tanh((theta-theta_on)/decay))-(1+tanh((theta-theta_off)/decay)));
+      double new_r = (0.3 - delta_r)*scale;
+      
+      v(0) = new_r*cos(theta_rad);
+      v(1) = new_r*sin(theta_rad); 
+      */
+      // SPARC, Z = 0 Limiter
+   
+      if (p(0) > 1.2685 && p(0) < 1.32 && p(1) > -0.5 && p(1) < 0.5)
+      {
+      //double xp = (p(1)-6.76128)/(-5.74468);
+      //double yp = -5.74468*p(0)+6.76128;
+
+      double scale = 1.0 - fabs(p(0)-1.269)/fabs(1.32-1.269);
+      double h = 0.03;
+      double decay = 0.01;  
+
+      double bump1 = h*0.5*((1+tanh((p(1)-0.45)/decay))-(1+tanh((p(1)-0.3)/decay)));
+      double bump2 = h*0.5*((1+tanh((p(1)-0.05)/decay))-(1+tanh((p(1)+0.05)/decay)));
+      double bump3 = h*0.5*((1+tanh((p(1)+0.3)/decay))-(1+tanh((p(1)+0.45)/decay)));
+      
+      v(0) = p(0) - (bump1+bump2+bump3)*scale;
+      v(1) = p(1);  
+      }
+      else {v = p;}
+
    }
    else
    {

miniapps/plasma/CMakeLists.txt

@@ -0,0 +1,43 @@
+# 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.
+
+if (MFEM_USE_MPI)
+  add_mfem_miniapp(stix1d
+    MAIN stix1d.cpp
+    EXTRA_SOURCES cold_plasma_dielectric_solver.cpp cold_plasma_dielectric_coefs.cpp g_eqdsk_data.cpp
+    EXTRA_HEADERS cold_plasma_dielectric_solver.hpp cold_plasma_dielectric_coefs.hpp g_eqdsk_data.hpp plasma.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(stix2d
+    MAIN stix2d.cpp
+    EXTRA_SOURCES cold_plasma_dielectric_solver.cpp cold_plasma_dielectric_coefs.cpp g_eqdsk_data.cpp interp_data.cpp
+    EXTRA_HEADERS cold_plasma_dielectric_solver.hpp cold_plasma_dielectric_coefs.hpp g_eqdsk_data.hpp interp_data.hpp plasma.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(stix3d
+    MAIN stix3d.cpp
+    EXTRA_SOURCES cold_plasma_dielectric_solver.cpp cold_plasma_dielectric_coefs.cpp g_eqdsk_data.cpp interp_data.cpp
+    EXTRA_HEADERS cold_plasma_dielectric_solver.hpp cold_plasma_dielectric_coefs.hpp g_eqdsk_data.hpp interp_data.hpp plasma.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(stix1d_dh
+    MAIN stix1d_dh.cpp
+    EXTRA_SOURCES cold_plasma_dielectric_solver.cpp cold_plasma_dielectric_coefs.cpp cold_plasma_dielectric_dh_solver.cpp g_eqdsk_data.cpp
+    EXTRA_HEADERS cold_plasma_dielectric_solver.hpp cold_plasma_dielectric_coefs.hpp cold_plasma_dielectric_dh_solver.hpp g_eqdsk_data.hpp plasma.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem mfem-common)
+
+  add_mfem_miniapp(stix2d_dh
+    MAIN stix2d_dh.cpp
+    EXTRA_SOURCES cold_plasma_dielectric_solver.cpp cold_plasma_dielectric_coefs.cpp cold_plasma_dielectric_dh_solver.cpp g_eqdsk_data.cpp
+    EXTRA_HEADERS cold_plasma_dielectric_solver.hpp cold_plasma_dielectric_coefs.hpp cold_plasma_dielectric_dh_solver.hpp g_eqdsk_data.hpp plasma.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
+    LIBRARIES mfem mfem-common)
+
+endif()

miniapps/plasma/Faddeeva.hh

@@ -0,0 +1,62 @@
+/* Copyright (c) 2012 Massachusetts Institute of Technology
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
+ */
+
+/* Available at: http://ab-initio.mit.edu/Faddeeva
+
+   Header file for Faddeeva.cc; see that file for more information. */
+
+#ifndef FADDEEVA_HH
+#define FADDEEVA_HH 1
+
+#include <complex>
+
+namespace Faddeeva {
+
+// compute w(z) = exp(-z^2) erfc(-iz) [ Faddeeva / scaled complex error func ]
+extern std::complex<double> w(std::complex<double> z,double relerr=0);
+extern double w_im(double x); // special-case code for Im[w(x)] of real x
+
+// Various functions that we can compute with the help of w(z)
+
+// compute erfcx(z) = exp(z^2) erfc(z)
+extern std::complex<double> erfcx(std::complex<double> z, double relerr=0);
+extern double erfcx(double x); // special case for real x
+
+// compute erf(z), the error function of complex arguments
+extern std::complex<double> erf(std::complex<double> z, double relerr=0);
+extern double erf(double x); // special case for real x
+
+// compute erfi(z) = -i erf(iz), the imaginary error function
+extern std::complex<double> erfi(std::complex<double> z, double relerr=0);
+extern double erfi(double x); // special case for real x
+
+// compute erfc(z) = 1 - erf(z), the complementary error function
+extern std::complex<double> erfc(std::complex<double> z, double relerr=0);
+extern double erfc(double x); // special case for real x
+
+// compute Dawson(z) = sqrt(pi)/2  *  exp(-z^2) * erfi(z)
+extern std::complex<double> Dawson(std::complex<double> z, double relerr=0);
+extern double Dawson(double x); // special case for real x
+
+} // namespace Faddeeva
+
+#endif // FADDEEVA_HH

miniapps/plasma/cold_plasma_dielectric_dh_solver.hpp

@@ -0,0 +1,867 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_COLD_PLASMA_DIELECTRIC_DH_SOLVER
+#define MFEM_COLD_PLASMA_DIELECTRIC_DH_SOLVER
+
+#include "../common/pfem_extras.hpp"
+#include "cold_plasma_dielectric_solver.hpp"
+#include "cold_plasma_dielectric_coefs.hpp"
+#include "plasma.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::L2_ParFESpace;
+using common::ParDiscreteGradOperator;
+using common::ParDiscreteCurlOperator;
+
+namespace plasma
+{
+/*
+// Solver options
+struct SolverOptions
+{
+   int maxIter;
+   int kDim;
+   int printLvl;
+   double relTol;
+
+   // Euclid Options
+   int euLvl;
+};
+*/
+/*
+struct ComplexCoefficientByAttr : public AttributeArrays
+{
+   Coefficient * real;
+   Coefficient * imag;
+};
+
+struct ComplexVectorCoefficientByAttr
+{
+   Array<int> attr;
+   Array<int> attr_marker;
+   VectorCoefficient * real;
+   VectorCoefficient * imag;
+};
+
+class ElectricEnergyDensityCoef : public Coefficient
+{
+public:
+   ElectricEnergyDensityCoef(VectorCoefficient &Er, VectorCoefficient &Ei,
+                             MatrixCoefficient &epsr, MatrixCoefficient &epsi);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   MatrixCoefficient &epsrCoef_;
+   MatrixCoefficient &epsiCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Dr_;
+   mutable Vector Di_;
+   mutable DenseMatrix eps_r_;
+   mutable DenseMatrix eps_i_;
+};
+
+class MagneticEnergyDensityCoef : public Coefficient
+{
+public:
+   MagneticEnergyDensityCoef(double omega,
+                             VectorCoefficient &dEr, VectorCoefficient &dEi,
+                             Coefficient &muInv);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   Coefficient &muInvCoef_;
+
+   mutable Vector Br_;
+   mutable Vector Bi_;
+};
+
+class EnergyDensityCoef : public Coefficient
+{
+public:
+   EnergyDensityCoef(double omega,
+                     VectorCoefficient &Er, VectorCoefficient &Ei,
+                     VectorCoefficient &dEr, VectorCoefficient &dEi,
+                     MatrixCoefficient &epsr, MatrixCoefficient &epsi,
+                     Coefficient &muInv);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   MatrixCoefficient &epsrCoef_;
+   MatrixCoefficient &epsiCoef_;
+   Coefficient &muInvCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Dr_;
+   mutable Vector Di_;
+   mutable Vector Br_;
+   mutable Vector Bi_;
+   mutable DenseMatrix eps_r_;
+   mutable DenseMatrix eps_i_;
+};
+*/
+class PoyntingVectorReCoefDH : public VectorCoefficient
+{
+public:
+   PoyntingVectorReCoefDH(double omega,
+                          VectorCoefficient &Er, VectorCoefficient &Ei,
+                          VectorCoefficient &Hr, VectorCoefficient &Hi);
+
+   void Eval(Vector &S, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &HrCoef_;
+   VectorCoefficient &HiCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+};
+
+class PoyntingVectorImCoefDH : public VectorCoefficient
+{
+public:
+   PoyntingVectorImCoefDH(double omega,
+                          VectorCoefficient &Er, VectorCoefficient &Ei,
+                          VectorCoefficient &Hr, VectorCoefficient &Hi);
+
+   void Eval(Vector &S, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &HrCoef_;
+   VectorCoefficient &HiCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+};
+
+class nxGradIntegrator : public BilinearFormIntegrator
+{
+private:
+
+   Coefficient *Q;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector nor, nxj;
+   DenseMatrix test_shape;
+   DenseMatrix trial_dshape;
+#endif
+
+public:
+   nxGradIntegrator() : Q(NULL) {}
+   nxGradIntegrator(Coefficient &q) : Q(&q) {}
+
+   int GetIntegrationOrder(const FiniteElement & trial_fe,
+                           const FiniteElement & test_fe,
+                           ElementTransformation &Trans)
+   { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
+
+   void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                               const FiniteElement &test_fe,
+                               ElementTransformation &Trans,
+                               DenseMatrix &elmat);
+};
+
+class nxkIntegrator : public BilinearFormIntegrator
+{
+private:
+
+   Coefficient *Q;
+   VectorCoefficient *K;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector nor, nxj, k;
+   DenseMatrix test_shape;
+   Vector trial_shape;
+#endif
+
+public:
+   nxkIntegrator(VectorCoefficient &k) : Q(NULL), K(&k) {}
+   nxkIntegrator(VectorCoefficient &k, Coefficient &q) : Q(&q), K(&k) {}
+
+   int GetIntegrationOrder(const FiniteElement & trial_fe,
+                           const FiniteElement & test_fe,
+                           ElementTransformation &Trans)
+   { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
+
+   void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                               const FiniteElement &test_fe,
+                               ElementTransformation &Trans,
+                               DenseMatrix &elmat);
+};
+
+class zkxIntegrator : public BilinearFormIntegrator
+{
+private:
+
+   Coefficient *Z;
+   VectorCoefficient *K;
+   double a;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector nor, nxj, k;
+   Vector test_shape;
+   DenseMatrix trial_shape;
+#endif
+
+public:
+   zkxIntegrator(Coefficient &z, VectorCoefficient &k, double _a = 1.0)
+      : Z(&z), K(&k), a(_a) {}
+
+   int GetIntegrationOrder(const FiniteElement & trial_fe,
+                           const FiniteElement & test_fe,
+                           ElementTransformation &Trans)
+   { return trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW(); }
+
+   void AssembleElementMatrix2(const FiniteElement &trial_fe,
+                               const FiniteElement &test_fe,
+                               ElementTransformation &Trans,
+                               DenseMatrix &elmat);
+};
+
+/// Cold Plasma Dielectric Solver
+class CPDSolverDH
+{
+public:
+
+   enum PrecondType
+   {
+      INVALID_PC  = -1,
+      DIAG_SCALE  =  1,
+      PARASAILS   =  2,
+      EUCLID      =  3,
+      AMS         =  4
+   };
+
+   enum SolverType
+   {
+      INVALID_SOL = -1,
+      GMRES       =  1,
+      FGMRES      =  2,
+      MINRES      =  3,
+      SUPERLU     =  4,
+      STRUMPACK   =  5,
+      DMUMPS      =  6,
+      ZMUMPS      =  7
+   };
+
+   CPDSolverDH(ParMesh & pmesh, int order, double omega,
+               CPDSolverDH::SolverType s, SolverOptions & sOpts,
+               CPDSolverDH::PrecondType p,
+               ComplexOperator::Convention conv,
+               VectorCoefficient & BCoef,
+               MatrixCoefficient & epsInvReCoef,
+               MatrixCoefficient & epsInvImCoef,
+               MatrixCoefficient & epsAbsCoef,
+               MatrixCoefficient & susceptReCoef,
+               MatrixCoefficient & susceptImCoef,
+               MatrixCoefficient & susceptReCoef_e,
+               MatrixCoefficient & susceptImCoef_e,
+               MatrixCoefficient & susceptReCoef_i1,
+               MatrixCoefficient & susceptImCoef_i1,
+               MatrixCoefficient * susceptReCoef_i2,
+               MatrixCoefficient * susceptImCoef_i2,
+               MatrixCoefficient * susceptReCoef_i3,
+               MatrixCoefficient * susceptImCoef_i3,
+               MatrixCoefficient & muReCoef,
+               MatrixCoefficient & muImCoef,
+               Coefficient & muCoef,
+               Coefficient * etaCoef,
+               VectorCoefficient * kReCoef,
+               VectorCoefficient * kImCoef,
+               Array<int> & abcs,
+               Array<ComplexVectorCoefficientByAttr*> & dbcs,
+               Array<ComplexVectorCoefficientByAttr*> & nbcs,
+               Array<ComplexCoefficientByAttr*> & sbcs,
+               void (*j_r_src)(const Vector&, Vector&),
+               void (*j_i_src)(const Vector&, Vector&),
+               bool vis_u = false,
+               bool pa = false,
+               bool dim2 = false);
+   ~CPDSolverDH();
+
+   HYPRE_Int GetProblemSize();
+
+   void PrintSizes();
+
+   void Assemble();
+
+   void Update();
+
+   void Solve();
+
+   double GetEFieldError(const VectorCoefficient & EReCoef,
+                         const VectorCoefficient & EImCoef) const;
+
+   double GetHFieldError(const VectorCoefficient & HReCoef,
+                         const VectorCoefficient & HImCoef) const;
+
+   void GetErrorEstimates(Vector & errors);
+
+   double GetVolume() const;
+
+   double GetGlobalDissipation() const;
+
+   double GetCoreDissipation() const;
+
+   double GetSOLDissipation() const;
+
+   double GetSheathDissipation() const;
+
+   void RegisterVisItFields(VisItDataCollection & visit_dc);
+
+   void WriteVisItFields(int it = 0);
+
+   void InitializeGLVis();
+
+   void DisplayToGLVis();
+
+   void DisplayAnimationToGLVis();
+
+   // const ParGridFunction & GetVectorPotential() { return *a_; }
+
+private:
+
+   class kekCoefficient : public MatrixCoefficient
+   {
+   private:
+      VectorCoefficient * krCoef_;
+      VectorCoefficient * kiCoef_;
+      MatrixCoefficient * erCoef_;
+      MatrixCoefficient * eiCoef_;
+
+      bool realPart_;
+      double a_;
+
+      mutable Vector kr;
+      mutable Vector ki;
+      mutable DenseMatrix er;
+      mutable DenseMatrix ei;
+
+      void kek(double a,
+               const Vector & kl, const DenseMatrix & e, const Vector &kr,
+               DenseMatrix & M)
+      {
+         for (int i=0; i<3; i++)
+         {
+            int i1 = (i+1)%3;
+            int i2 = (i+2)%3;
+            for (int j=0; j<3; j++)
+            {
+               int j1 = (j+1)%3;
+               int j2 = (j+2)%3;
+               M(i,j) +=
+                  a * (kl(i2) * e(i1,j2) * kr(j1) -
+                       kl(i2) * e(i1,j1) * kr(j2) -
+                       kl(i1) * e(i2,j2) * kr(j1) +
+                       kl(i1) * e(i2,j1) * kr(j2)
+                      );
+            }
+         }
+      }
+
+   public:
+      kekCoefficient(VectorCoefficient *krCoef, VectorCoefficient *kiCoef,
+                     MatrixCoefficient *erCoef, MatrixCoefficient *eiCoef,
+                     bool realPart, double a = 1.0)
+         : MatrixCoefficient(3),
+           krCoef_(krCoef), kiCoef_(kiCoef),
+           erCoef_(erCoef), eiCoef_(eiCoef),
+           realPart_(realPart),
+           a_(a), kr(3), ki(3), er(3), ei(3)
+      { kr = 0.0; ki = 0.0; er = 0.0; ei = 0.0; }
+
+      void Eval(DenseMatrix &M, ElementTransformation &T,
+                const IntegrationPoint &ip)
+      {
+         M.SetSize(3);
+         M = 0.0;
+         if ((krCoef_ == NULL && kiCoef_ == NULL) ||
+             (erCoef_ == NULL && eiCoef_ == NULL))
+         {
+            return;
+         }
+         if (krCoef_) { krCoef_->Eval(kr, T, ip); }
+         if (kiCoef_) { kiCoef_->Eval(ki, T, ip); }
+         if (erCoef_) { erCoef_->Eval(er, T, ip); }
+         if (eiCoef_) { eiCoef_->Eval(ei, T, ip); }
+
+         if (realPart_)
+         {
+            if (krCoef_ && erCoef_) { kek(1.0, kr, er, kr, M); }
+            if (kiCoef_ && erCoef_) { kek(-1.0, ki, er, ki, M); }
+            if (krCoef_ && eiCoef_ && kiCoef_) { kek(-1.0, kr, ei, ki, M); }
+            if (kiCoef_ && eiCoef_ && krCoef_) { kek(-1.0, ki, ei, kr, M); }
+         }
+         else
+         {
+            if (krCoef_ && eiCoef_) { kek(1.0, kr, ei, kr, M); }
+            if (kiCoef_ && eiCoef_) { kek(-1.0, ki, ei, ki, M); }
+            if (krCoef_ && erCoef_ && kiCoef_) { kek(1.0, kr, er, ki, M); }
+            if (kiCoef_ && erCoef_ && krCoef_) { kek(1.0, ki, er, kr, M); }
+         }
+         if (a_ != 1.0) { M *= a_; }
+      }
+   };
+
+   class ekCoefficient : public MatrixCoefficient
+   {
+   private:
+      VectorCoefficient * krCoef_;
+      VectorCoefficient * kiCoef_;
+      MatrixCoefficient * erCoef_;
+      MatrixCoefficient * eiCoef_;
+
+      bool realPart_;
+      double a_;
+
+      mutable Vector kr;
+      mutable Vector ki;
+      mutable DenseMatrix er;
+      mutable DenseMatrix ei;
+
+      void ek(double a,
+              const DenseMatrix & e, const Vector &k,
+              DenseMatrix & M)
+      {
+         for (int i=0; i<3; i++)
+         {
+            for (int j=0; j<3; j++)
+            {
+               int j1 = (j+1)%3;
+               int j2 = (j+2)%3;
+               M(i,j) += a * (e(i,j1) * k(j2) - e(i,j2) * k(j1));
+            }
+         }
+      }
+
+   public:
+      ekCoefficient(VectorCoefficient *krCoef, VectorCoefficient *kiCoef,
+                    MatrixCoefficient *erCoef, MatrixCoefficient *eiCoef,
+                    bool realPart,
+                    double a = 1.0)
+         : MatrixCoefficient(3),
+           krCoef_(krCoef), kiCoef_(kiCoef),
+           erCoef_(erCoef), eiCoef_(eiCoef),
+           realPart_(realPart),
+           a_(a), kr(3), ki(3), er(3), ei(3)
+      { kr = 0.0; ki = 0.0; er = 0.0; ei = 0.0; }
+
+      void Eval(DenseMatrix &M, ElementTransformation &T,
+                const IntegrationPoint &ip)
+      {
+         M.SetSize(3);
+         M = 0.0;
+         if ((krCoef_ == NULL && kiCoef_ == NULL) ||
+             (erCoef_ == NULL && eiCoef_ == NULL))
+         {
+            return;
+         }
+         if (krCoef_) { krCoef_->Eval(kr, T, ip); }
+         if (kiCoef_) { kiCoef_->Eval(ki, T, ip); }
+         if (erCoef_) { erCoef_->Eval(er, T, ip); }
+         if (eiCoef_) { eiCoef_->Eval(ei, T, ip); }
+
+         if (realPart_)
+         {
+            if (erCoef_ && krCoef_) { ek(1.0, er, kr, M); }
+            if (eiCoef_ && kiCoef_) { ek(-1.0, ei, ki, M); }
+         }
+         else
+         {
+            if (eiCoef_ && krCoef_) { ek(1.0, ei, kr, M); }
+            if (erCoef_ && kiCoef_) { ek(1.0, er, ki, M); }
+         }
+         if (a_ != 1.0) { M *= a_; }
+      }
+   };
+
+   class keCoefficient : public MatrixCoefficient
+   {
+   private:
+      VectorCoefficient * krCoef_;
+      VectorCoefficient * kiCoef_;
+      MatrixCoefficient * erCoef_;
+      MatrixCoefficient * eiCoef_;
+
+      bool realPart_;
+      double a_;
+
+      mutable Vector kr;
+      mutable Vector ki;
+      mutable DenseMatrix er;
+      mutable DenseMatrix ei;
+
+      void ke(double a,
+              const Vector &k, const DenseMatrix & e,
+              DenseMatrix & M)
+      {
+         for (int i=0; i<3; i++)
+         {
+            int i1 = (i+1)%3;
+            int i2 = (i+2)%3;
+            for (int j=0; j<3; j++)
+            {
+               M(i,j) += a * (k(i1) * e(i2,j) - k(i2) * e(i1,j));
+            }
+         }
+      }
+
+   public:
+      keCoefficient(VectorCoefficient *krCoef, VectorCoefficient *kiCoef,
+                    MatrixCoefficient *erCoef, MatrixCoefficient *eiCoef,
+                    bool realPart,
+                    double a = 1.0)
+         : MatrixCoefficient(3),
+           krCoef_(krCoef), kiCoef_(kiCoef),
+           erCoef_(erCoef), eiCoef_(eiCoef),
+           realPart_(realPart),
+           a_(a), kr(3), ki(3), er(3), ei(3)
+      { kr = 0.0; ki = 0.0; er = 0.0; ei = 0.0; }
+
+      void Eval(DenseMatrix &M, ElementTransformation &T,
+                const IntegrationPoint &ip)
+      {
+         M.SetSize(3);
+         M = 0.0;
+         if ((krCoef_ == NULL && kiCoef_ == NULL) ||
+             (erCoef_ == NULL && eiCoef_ == NULL))
+         {
+            return;
+         }
+         if (krCoef_) { krCoef_->Eval(kr, T, ip); }
+         if (kiCoef_) { kiCoef_->Eval(ki, T, ip); }
+         if (erCoef_) { erCoef_->Eval(er, T, ip); }
+         if (eiCoef_) { eiCoef_->Eval(ei, T, ip); }
+
+         if (realPart_)
+         {
+            if (krCoef_ && erCoef_) { ke(1.0, kr, er, M); }
+            if (kiCoef_ && eiCoef_) { ke(-1.0, ki, ei, M); }
+         }
+         else
+         {
+            if (krCoef_ && eiCoef_) { ke(1.0, kr, ei, M); }
+            if (kiCoef_ && erCoef_) { ke(1.0, ki, er, M); }
+         }
+         if (a_ != 1.0) { M *= a_; }
+      }
+   };
+
+   void collectBdrAttributes(const Array<AttributeArrays*> & aa,
+                             Array<int> & attr_marker);
+
+   void locateTrueDBCDofs(const Array<int> & dbc_bdr_marker,
+                          Array<int> & dbc_nd_tdofs);
+
+   void locateTrueSBCDofs(const Array<int> & sbc_bdr_marker,
+                          Array<int> & non_sbc_h1_tdofs,
+                          Array<int> & sbc_nd_tdofs);
+
+   void computeD(const ParComplexGridFunction & h,
+                 const ParComplexGridFunction & j,
+                 ParComplexGridFunction & d);
+
+   void computeE(const ParComplexGridFunction & d,
+                 ParComplexGridFunction & e);
+
+   int myid_;
+   int num_procs_;
+   int order_;
+   int logging_;
+
+   SolverType sol_;
+   SolverOptions & solOpts_;
+   PrecondType prec_;
+
+   ComplexOperator::Convention conv_;
+
+   bool ownsEta_;
+   bool vis_u_;
+   bool pa_;
+   bool dim2_;
+
+   double omega_;
+   int H_iter_;
+
+   // double solNorm_;
+
+   ParMesh * pmesh_;
+
+   L2_ParFESpace * L2FESpace_;
+   L2_ParFESpace * L2FESpace2p_;
+   L2_ParFESpace * L2VFESpace_;
+   H1_ParFESpace * H1FESpace_;
+   ND_ParFESpace * HCurlFESpace_;
+   RT_ParFESpace * HDivFESpace_;
+   RT_ParFESpace * HDivFESpace2p_;
+
+   Array<HYPRE_Int> blockTrueOffsets_;
+
+   // ParSesquilinearForm * a0_;
+   ParSesquilinearForm * a1_;
+   // ParBilinearForm * b1_;
+   ParMixedSesquilinearForm * nxD01_;
+   ParMixedSesquilinearForm * d21EpsInv_;
+
+   ParSesquilinearForm * m1_;
+   ParMixedSesquilinearForm * m21EpsInv_;
+
+   // ParBilinearForm * m0_;
+   // ParMixedBilinearForm * n20ZRe_;
+   // ParMixedBilinearForm * n20ZIm_;
+
+   ConstantCoefficient negOneCoef_;
+   ParSesquilinearForm * m0_;
+   ParMixedSesquilinearForm * nzD12_;
+   ParBilinearForm * m3_;
+   ParBilinearForm * m4r_;
+   ParBilinearForm * m4i_;
+   ParBilinearForm * m4cr_;
+   ParBilinearForm * m4ci_;
+   ParBilinearForm * m4solr_;
+   ParBilinearForm * m4soli_;
+
+   ParDiscreteGradOperator * grad_; // For Computing E from phi
+   ParDiscreteCurlOperator * curl_; // For Computing D from H
+   ParDiscreteLinearOperator * kReCross_;
+   ParDiscreteLinearOperator * kImCross_;
+
+   ParComplexGridFunction * h_;   // Complex magnetic field (HCurl)
+   ParGridFunction * hr_;   // Real component to magnetic field (HCurl)
+   ParGridFunction * hi_;   // Imag component to magnetic field (HCurl)
+   //ParGridFunction * ht_real_;   // Tangential real component to magnetic field (HCurl)
+   //ParGridFunction * ht_imag_;   // Tangential imag component to magnetic field (HCurl)
+   ParComplexGridFunction * e_;   // Complex electric field (HCurl)
+   ParComplexGridFunction * d_;   // Complex electric flux (HDiv)
+   ParComplexGridFunction * j_;   // Complex current density (HDiv)
+   ParComplexLinearForm   * curlj_; // Curl of current density (HCurl)
+   ParComplexGridFunction * phi_; // Complex sheath potential (H1)
+   ParComplexGridFunction * prev_phi_; // Complex sheath potential (H1)
+
+   // ParComplexGridFunction * e_tmp_; // Temporary complex electric field (HCurl)
+   // ParGridFunction * temp_; // Temporary grid function (HCurl)
+   // ParComplexGridFunction * phi_tmp_; // Complex sheath potential temporary (H1)
+   ParGridFunction * rectPot_; // Real valued rectified potential (H1)
+   ParGridFunction * sheath_pow_; // Real valued sheath power dissipation (H1)
+   ParGridFunction * Bn_; // Real valued angle of B field into boundary(H1)
+   // ParComplexGridFunction * j_;   // Complex current density (HCurl)
+   ParComplexLinearForm   * rhs1_; // RHS of magnetic field eqn (HCurl)
+   ParComplexLinearForm   * rhs0_; // RHS of sheath potential eqn (H1)
+   ParGridFunction        * e_t_; // Time dependent Electric field
+   ParComplexGridFunction * e_b_; // Complex parallel magnetic field (L2)
+   ParComplexGridFunction * h_b_; // Complex parallel electric field (L2)
+   ParComplexGridFunction * e_perp_; // Complex perpendicular electric field (L2)
+   ParComplexGridFunction * e_plus_; // Complex + polarized electric field (L2)
+   ParComplexGridFunction * e_min_; // Complex - polarized electric field (L2)
+   ParGridFunction        * power_absorp_t_; // Real valued total power absorption (H1)
+   ParGridFunction        * power_absorp_e_; // Real valued electron power absorption (H1)
+   ParGridFunction        * power_absorp_i1_; // Real valued ion 1 absorption (H1)
+   ParGridFunction        * power_absorp_i2_; // Real valued ion 2 power absorption (H1)
+   ParGridFunction        * power_absorp_i3_; // Real valued ion 3 power absorption (H1)
+   ParComplexGridFunction * h_v_; // Complex magnetic field (L2^d)
+   ParComplexGridFunction * e_v_; // Complex electric field (L2^d)
+   ParComplexGridFunction * d_v_; // Complex electric flux (L2^d)
+   ParComplexGridFunction * phi_v_; // Complex sheath potential (L2)
+   ParComplexGridFunction * j_v_; // Complex current density (L2^d)
+   ParGridFunction        * b_hat_; // Unit vector along B (HDiv)
+   // ParGridFunction        * u_;   // Energy density (L2)
+   // ParGridFunction        * uE_;  // Electric Energy density (L2)
+   // ParGridFunction        * uB_;  // Magnetic Energy density (L2)
+   ParComplexGridFunction * S_;  // Poynting Vector (HCurl)
+   ParComplexGridFunction * StixS_; // Stix S Coefficient (L2)
+   ParComplexGridFunction * StixD_; // Stix D Coefficient (L2)
+   ParComplexGridFunction * StixP_; // Stix P Coefficient (L2)
+   ParComplexGridFunction * EpsPara_; // B^T eps B / |B|^2 Coefficient (L2)
+
+   HypreParMatrix * M3_;
+   HypreParMatrix * M4r_;
+   HypreParMatrix * M4i_;
+   HypreParMatrix * M4cr_;
+   HypreParMatrix * M4ci_;
+   HypreParMatrix * M4solr_;
+   HypreParMatrix * M4soli_;
+   HypreParVector * PHIr_; 
+   HypreParVector * PHIi_; 
+   mutable HypreParVector * RHSr1_;
+   mutable HypreParVector * RHSi1_;
+   mutable HypreParVector * RHSr2_;
+   mutable HypreParVector * RHSi2_;
+   mutable HypreParVector * RHSr3_;
+   mutable HypreParVector * RHSi3_;
+   mutable HypreParVector * RHSr4_;
+   mutable HypreParVector * RHSi4_;
+   mutable HypreParVector * TMPr2_;
+   mutable HypreParVector * TMPi2_;
+   mutable HypreParVector * TMPr3_;
+   mutable HypreParVector * TMPi3_;
+   mutable HypreParVector * TMPr4_;
+   mutable HypreParVector * TMPi4_;
+   HypreParVector * Er_; 
+   HypreParVector * Ei_;  
+
+   VectorCoefficient * BCoef_;        // B Field Unit Vector
+   // MatrixCoefficient * epsReCoef_;    // Dielectric Material Coefficient
+   // MatrixCoefficient * epsImCoef_;    // Dielectric Material Coefficient
+   MatrixCoefficient * epsInvReCoef_;    // Dielectric Material Coefficient
+   MatrixCoefficient * epsInvImCoef_;    // Dielectric Material Coefficient
+   MatrixCoefficient * susceptReCoef_;    // Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_;    // Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_e_;    // Electron Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_e_;    // Electron Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i1_;    // Ion 1 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i1_;    // Ion 1 Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i2_;    // Ion 2 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i2_;    // Ion 2 Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i3_;    // Ion 3 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i3_;    // Ion 3 Imag Susceptibility Coefficient
+   // MatrixCoefficient * epsAbsCoef_;   // Dielectric Material Coefficient
+   Coefficient * muCoef_;       // Real Dia/Paramagnetic Material Coefficient
+   MatrixCoefficient * muReCoef_;       // Real Dia/Paramagnetic Material Coefficient
+   MatrixCoefficient * muImCoef_;       // Imag Dia/Paramagnetic Material Coefficient
+   PowerCoefficient    muInvReCoef_;    // Dia/Paramagnetic Material Coefficient
+   Coefficient       * etaCoef_;      // Impedance Coefficient
+   SheathPower       * sheathPowCoef_; // Sheath Width * Admittance Coefficient
+   VectorCoefficient * kReCoef_;        // Wave Vector
+   VectorCoefficient * kImCoef_;        // Wave Vector
+
+   Coefficient * SReCoef_; // Stix S Coefficient
+   Coefficient * SImCoef_; // Stix S Coefficient
+   Coefficient * DReCoef_; // Stix D Coefficient
+   Coefficient * DImCoef_; // Stix D Coefficient
+   Coefficient * PReCoef_; // Stix P Coefficient
+   Coefficient * PImCoef_; // Stix P 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 * abcCoef_;       // -omega eta
+   // Coefficient * sbcReCoef_;     //  omega Im(eta^{-1})
+   // Coefficient * sbcImCoef_;     // -omega Re(eta^{-1})
+   Coefficient * sinkx_;         // sin(ky * y + kz * z)
+   Coefficient * coskx_;         // cos(ky * y + kz * z)
+   Coefficient * negsinkx_;      // -sin(ky * y + kz * z)
+   // Coefficient * negMuInvCoef_;  // -1.0 / mu
+
+   Coefficient * massCoef_;  // -omega^2 mu
+   MatrixCoefficient * massReCoef_;  // -omega^2 mu_r
+   MatrixCoefficient * massImCoef_;  // -omega^2 mu_i
+   Coefficient * posMassCoef_; // omega^2 mu
+   // MatrixCoefficient * negMuInvkxkxCoef_; // -\vec{k}\times\vec{k}\times/mu
+
+   // VectorCoefficient * negMuInvkCoef_; // -\vec{k}/mu
+
+   kekCoefficient kekReCoef_;
+   kekCoefficient kekImCoef_;
+   keCoefficient keReCoef_;
+   keCoefficient keImCoef_;
+   ekCoefficient ekReCoef_;
+   ekCoefficient ekImCoef_;
+
+   VectorCoefficient * jrCoef_;     // Volume Current Density Function
+   VectorCoefficient * jiCoef_;     // Volume Current Density Function
+   VectorCoefficient * rhsrCoef_;     // Volume Current Density Function
+   VectorCoefficient * rhsiCoef_;     // Volume Current Density Function
+
+   VectorGridFunctionCoefficient erCoef_;
+   VectorGridFunctionCoefficient eiCoef_;
+
+   VectorGridFunctionCoefficient hrCoef_;
+   VectorGridFunctionCoefficient hiCoef_;
+
+   // EnergyDensityCoef     uCoef_;
+   // ElectricEnergyDensityCoef uECoef_;
+   // MagneticEnergyDensityCoef uBCoef_;
+   PoyntingVectorReCoefDH SrCoef_;
+   PoyntingVectorImCoefDH SiCoef_;
+
+   // const VectorCoefficient & erCoef_;     // Electric Field Boundary Condition
+   // const 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_bdr_marker_;
+
+   Array<ComplexVectorCoefficientByAttr*> * dbcs_;
+   Array<int> dbc_bdr_marker_;
+   Array<int> dbc_nd_tdofs_;
+   Array<int> non_k_bdr_;
+
+   Array<ComplexVectorCoefficientByAttr*> * nbcs_; // Surface current BCs
+   Array<ComplexVectorCoefficientByAttr*> * nkbcs_; // Neumann BCs (-i*omega*K)
+
+   Array<ComplexCoefficientByAttr*> * sbcs_; // Sheath BCs
+   Array<int> sbc_bdr_marker_;
+   Array<int> non_sbc_h1_tdofs_;
+   Array<int> sbc_nd_tdofs_;
+
+   Array<int> core_attr_marker_;
+   Array<int> sol_attr_marker_;
+
+   VisItDataCollection * visit_dc_;
+
+   std::map<std::string,socketstream*> socks_;
+};
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_COLD_PLASMA_DIELECTRIC_DH_SOLVER

miniapps/plasma/cold_plasma_dielectric_solver.hpp

@@ -0,0 +1,751 @@
+// 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_COLD_PLASMA_DIELECTRIC_SOLVER
+#define MFEM_COLD_PLASMA_DIELECTRIC_SOLVER
+
+#include "../common/pfem_extras.hpp"
+#include "plasma.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::L2_ParFESpace;
+using common::ParDiscreteGradOperator;
+using common::ParDiscreteCurlOperator;
+
+namespace plasma
+{
+
+// Solver options
+struct SolverOptions
+{
+   int maxIter;
+   int kDim;
+   int printLvl;
+   double relTol;
+   double absTol;
+
+   // Euclid Options
+   int euLvl;
+};
+
+struct AttributeArrays
+{
+   Array<int> attr;
+   Array<int> attr_marker;
+
+};
+
+struct ComplexCoefficientByAttr : public AttributeArrays
+{
+   Coefficient * real;
+   Coefficient * imag;
+};
+
+struct ComplexVectorCoefficientByAttr : public AttributeArrays
+{
+   VectorCoefficient * real;
+   VectorCoefficient * imag;
+};
+
+// Used for combining scalar coefficients
+double prodFunc(double a, double b);
+
+class ElectricEnergyDensityCoef : public Coefficient
+{
+public:
+   ElectricEnergyDensityCoef(VectorCoefficient &Er, VectorCoefficient &Ei,
+                             MatrixCoefficient &epsr, MatrixCoefficient &epsi);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   MatrixCoefficient &epsrCoef_;
+   MatrixCoefficient &epsiCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Dr_;
+   mutable Vector Di_;
+   mutable DenseMatrix eps_r_;
+   mutable DenseMatrix eps_i_;
+};
+
+class MagneticEnergyDensityCoef : public Coefficient
+{
+public:
+   MagneticEnergyDensityCoef(double omega,
+                             VectorCoefficient &dEr, VectorCoefficient &dEi,
+                             MatrixCoefficient &muInvRe, MatrixCoefficient &muInvIm);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   MatrixCoefficient &muInvReCoef_;
+   MatrixCoefficient &muInvImCoef_;
+
+   mutable Vector Br_;
+   mutable Vector Bi_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+   mutable DenseMatrix muInvRe_;
+   mutable DenseMatrix muInvIm_;
+};
+
+class EnergyDensityCoef : public Coefficient
+{
+public:
+   EnergyDensityCoef(double omega,
+                     VectorCoefficient &Er, VectorCoefficient &Ei,
+                     VectorCoefficient &dEr, VectorCoefficient &dEi,
+                     MatrixCoefficient &epsr, MatrixCoefficient &epsi,
+                     MatrixCoefficient &muInvRe, MatrixCoefficient &muInvIm);
+
+   double Eval(ElementTransformation &T,
+               const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   MatrixCoefficient &epsrCoef_;
+   MatrixCoefficient &epsiCoef_;
+   MatrixCoefficient &muInvReCoef_;
+   MatrixCoefficient &muInvImCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Dr_;
+   mutable Vector Di_;
+   mutable Vector Br_;
+   mutable Vector Bi_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+   mutable DenseMatrix eps_r_;
+   mutable DenseMatrix eps_i_;
+   mutable DenseMatrix muInvRe_;
+   mutable DenseMatrix muInvIm_;
+};
+
+class PoyntingVectorReCoef : public VectorCoefficient
+{
+public:
+   PoyntingVectorReCoef(double omega,
+                        VectorCoefficient &Er, VectorCoefficient &Ei,
+                        VectorCoefficient &dEr, VectorCoefficient &dEi,
+                        MatrixCoefficient &muInvRe, MatrixCoefficient &muInvIm);
+
+   void Eval(Vector &S, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   MatrixCoefficient &muInvReCoef_;
+   MatrixCoefficient &muInvImCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Br_;
+   mutable Vector Bi_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+   mutable DenseMatrix muInvRe_;
+   mutable DenseMatrix muInvIm_;
+};
+
+class PoyntingVectorImCoef : public VectorCoefficient
+{
+public:
+   PoyntingVectorImCoef(double omega,
+                        VectorCoefficient &Er, VectorCoefficient &Ei,
+                        VectorCoefficient &dEr, VectorCoefficient &dEi,
+                        MatrixCoefficient &muInvRe, MatrixCoefficient &muInvIm);
+
+   void Eval(Vector &S, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+private:
+   double omega_;
+
+   VectorCoefficient &ErCoef_;
+   VectorCoefficient &EiCoef_;
+   VectorCoefficient &dErCoef_;
+   VectorCoefficient &dEiCoef_;
+   MatrixCoefficient &muInvReCoef_;
+   MatrixCoefficient &muInvImCoef_;
+
+   mutable Vector Er_;
+   mutable Vector Ei_;
+   mutable Vector Br_;
+   mutable Vector Bi_;
+   mutable Vector Hr_;
+   mutable Vector Hi_;
+   mutable DenseMatrix muInvRe_;
+   mutable DenseMatrix muInvIm_;
+};
+
+/// Cold Plasma Dielectric Solver
+class CPDSolver
+{
+public:
+
+   enum PrecondType
+   {
+      INVALID_PC  = -1,
+      DIAG_SCALE  =  1,
+      PARASAILS   =  2,
+      EUCLID      =  3,
+      AMS         =  4
+   };
+
+   enum SolverType
+   {
+      INVALID_SOL = -1,
+      GMRES       =  1,
+      FGMRES      =  2,
+      MINRES      =  3,
+      SUPERLU     =  4,
+      STRUMPACK   =  5,
+      DMUMPS      =  6,
+      ZMUMPS      =  7
+   };
+
+   CPDSolver(ParMesh & pmesh, int order, double omega,
+             CPDSolver::SolverType s, SolverOptions & sOpts,
+             CPDSolver::PrecondType p,
+             ComplexOperator::Convention conv,
+             VectorCoefficient & BCoef,
+             MatrixCoefficient & epsReCoef,
+             MatrixCoefficient & epsImCoef,
+             MatrixCoefficient & epsAbsCoef,
+             MatrixCoefficient & susceptReCoef,
+             MatrixCoefficient & susceptImCoef,
+             MatrixCoefficient & susceptReCoef_e,
+             MatrixCoefficient & susceptImCoef_e,
+             MatrixCoefficient & susceptReCoef_i1,
+             MatrixCoefficient & susceptImCoef_i1,
+             MatrixCoefficient * susceptReCoef_i2,
+             MatrixCoefficient * susceptImCoef_i2,
+             MatrixCoefficient * susceptReCoef_i3,
+             MatrixCoefficient * susceptImCoef_i3,
+             MatrixCoefficient & muInvReCoef,
+             MatrixCoefficient & muInvImCoef,
+             Coefficient * etaInvCoef,
+             VectorCoefficient * kReCoef,
+             VectorCoefficient * kImCoef,
+             Array<int> & abcs,
+             Array<ComplexVectorCoefficientByAttr*> & dbcs,
+             Array<ComplexVectorCoefficientByAttr*> & nbcs,
+             Array<ComplexCoefficientByAttr*> & sbcs,
+             void (*j_r_src)(const Vector&, Vector&),
+             void (*j_i_src)(const Vector&, Vector&),
+             bool vis_u = false,
+             bool pa = false,
+             bool dim2 = false);
+   ~CPDSolver();
+
+   HYPRE_Int GetProblemSize();
+
+   void PrintSizes();
+
+   void Assemble();
+
+   void Update();
+
+   void Solve();
+
+   double GetError(const VectorCoefficient & EReCoef,
+                   const VectorCoefficient & EImCoef) const;
+
+   void GetErrorEstimates(Vector & errors);
+
+   double GetVolume() const;
+
+   double GetGlobalDissipation() const;
+
+   double GetElectronDissipation() const;
+
+   double GetIon1Dissipation() const;
+
+   double GetIon2Dissipation() const;
+
+   double GetIon3Dissipation() const;
+
+   //double GetCoreDissipation() const;
+
+   //double GetSOLDissipation() const;
+
+   void RegisterVisItFields(VisItDataCollection & visit_dc);
+
+   void WriteVisItFields(int it = 0);
+
+   void InitializeGLVis();
+
+   void DisplayToGLVis();
+
+   void DisplayAnimationToGLVis();
+
+   // const ParGridFunction & GetVectorPotential() { return *a_; }
+
+private:
+
+   class kmkCoefficient : public MatrixCoefficient
+   {
+   private:
+      VectorCoefficient * krCoef_;
+      VectorCoefficient * kiCoef_;
+      MatrixCoefficient * mReCoef_;
+      MatrixCoefficient * mImCoef_;
+
+      bool realPart_;
+      double a_;
+
+      mutable Vector kr;
+      mutable Vector ki;
+      mutable Vector mukr;
+      mutable Vector muki;
+      mutable DenseMatrix muInvRe_;
+      mutable DenseMatrix muInvIm_;
+
+      void kmk(double a,
+               const Vector & kl, const Vector &kr,
+               DenseMatrix & M)
+      {
+         double kk = kl * kr;
+         for (int i=0; i<3; i++)
+         {
+            for (int j=0; j<3; j++)
+            {
+               M(i,j) += a * kl(j) * kr(i);
+            }
+            M(i,i) -= a * kk;
+         }
+      }
+
+   public:
+      kmkCoefficient(VectorCoefficient *krCoef, VectorCoefficient *kiCoef,
+                     MatrixCoefficient *mReCoef, MatrixCoefficient *mImCoef,
+                     bool realPart, double a = 1.0)
+         : MatrixCoefficient(3),
+           krCoef_(krCoef), kiCoef_(kiCoef),
+           mReCoef_(mReCoef),
+           mImCoef_(mImCoef),
+           realPart_(realPart),
+           a_(a), kr(3), ki(3), mukr(3), muki(3), muInvRe_(3), muInvIm_(3)
+      { kr = 0.0; ki = 0.0; mukr = 0.0; muki = 0.0; muInvRe_ = 0.0; muInvIm_ = 0.0;}
+
+      void Eval(DenseMatrix &M, ElementTransformation &T,
+                const IntegrationPoint &ip)
+      {
+         M.SetSize(3);
+         M = 0.0;
+         if ((krCoef_ == NULL && kiCoef_ == NULL) ||
+             mReCoef_ == NULL)
+         {
+            return;
+         }
+
+         if (krCoef_) { krCoef_->Eval(kr, T, ip); }
+         if (kiCoef_) { kiCoef_->Eval(ki, T, ip); }
+         if (mReCoef_) { mReCoef_->Eval(muInvRe_, T, ip); }
+         if (mImCoef_) { mImCoef_->Eval(muInvIm_, T, ip); }
+
+         // real
+         muInvRe_.Mult(kr, mukr);
+         muInvIm_.AddMult_a(-1.0, ki, mukr);
+
+         // imag
+         muInvIm_.Mult(kr, muki);
+         muInvRe_.AddMult(ki, muki);
+
+         if (realPart_)
+         {
+            if (krCoef_) { kmk(1.0, kr, mukr, M); }
+            if (kiCoef_) { kmk(-1.0, ki, muki, M); }
+         }
+         else
+         {
+            if (krCoef_ && kiCoef_) { kmk(1.0, kr, muki, M); }
+            if (kiCoef_ && krCoef_) { kmk(1.0, ki, mukr, M); }
+         }
+         if (a_ != 1.0) { M *= a_; }
+      }
+   };
+
+   class CrossCoefficient : public MatrixCoefficient
+   {
+   private:
+      VectorCoefficient * krCoef_;
+      VectorCoefficient * kiCoef_;
+      MatrixCoefficient * mReCoef_;
+      MatrixCoefficient * mImCoef_;
+
+      bool realPart_;
+      double a_;
+
+      mutable Vector kr;
+      mutable Vector ki;
+      mutable Vector mukr;
+      mutable Vector muki;
+      mutable DenseMatrix muInvRe_;
+      mutable DenseMatrix muInvIm_;
+
+   public:
+      CrossCoefficient(VectorCoefficient *krCoef, VectorCoefficient *kiCoef,
+                       MatrixCoefficient *mReCoef, MatrixCoefficient *mImCoef,
+                       bool realPart, double a = 1.0)
+         : MatrixCoefficient(3),
+           krCoef_(krCoef),
+           kiCoef_(kiCoef),
+           mReCoef_(mReCoef),
+           mImCoef_(mImCoef),
+           realPart_(realPart),
+           a_(a), kr(3), ki(3), mukr(3), muki(3), muInvRe_(3), muInvIm_(3)
+      { kr = 0.0; ki = 0.0; mukr = 0.0; muki = 0.0; muInvRe_ = 0.0; muInvIm_ = 0.0;}
+
+      void Eval(DenseMatrix &M, ElementTransformation &T,
+                const IntegrationPoint &ip)
+      {
+         M.SetSize(3);
+         M = 0.0;
+
+         if (krCoef_) { krCoef_->Eval(kr, T, ip); }
+         if (kiCoef_) { kiCoef_->Eval(ki, T, ip); }
+         if (mReCoef_) { mReCoef_->Eval(muInvRe_, T, ip); }
+         if (mImCoef_) { mImCoef_->Eval(muInvIm_, T, ip); }
+
+         // real
+         muInvRe_.Mult(kr, mukr);
+         muInvIm_.AddMult_a(-1.0, ki, mukr);
+
+         // imag
+         muInvIm_.Mult(kr, muki);
+         muInvRe_.AddMult(ki, muki);
+
+         if (realPart_)
+         {
+            M(2,1) = a_ * mukr(0);
+            M(0,2) = a_ * mukr(1);
+            M(1,0) = a_ * mukr(2);            
+         }
+         else
+         {
+            M(2,1) = a_ * muki(0);
+            M(0,2) = a_ * muki(1);
+            M(1,0) = a_ * muki(2);
+         }
+
+         M(1,2) = -M(2,1);
+         M(2,0) = -M(0,2);
+         M(0,1) = -M(1,0);
+      }
+   };
+
+   void computeB(const ParComplexGridFunction & e,
+                 ParComplexGridFunction & b);
+
+   void computeD(const ParComplexGridFunction & e,
+                 ParComplexGridFunction & d);
+
+   int myid_;
+   int num_procs_;
+   int order_;
+   int logging_;
+
+   SolverType sol_;
+   SolverOptions & solOpts_;
+   PrecondType prec_;
+
+   ComplexOperator::Convention conv_;
+
+   bool ownsEtaInv_;
+   bool vis_u_;
+   bool pa_;
+   bool dim2_;
+
+   double omega_;
+
+   // double solNorm_;
+
+   ParMesh * pmesh_;
+
+   L2_ParFESpace * L2FESpace_;
+   L2_ParFESpace * L2FESpace2p_;
+   L2_ParFESpace * L2VFESpace_;
+   H1_ParFESpace * H1FESpace_;
+   ND_ParFESpace * HCurlFESpace_;
+   RT_ParFESpace * HDivFESpace_;
+   RT_ParFESpace * HDivFESpace2p_;
+
+   Array<HYPRE_Int> blockTrueOffsets_;
+
+   // ParSesquilinearForm * a0_;
+   ParSesquilinearForm * a1_;
+   ParBilinearForm * b1_;
+
+   ParBilinearForm * m2_;
+   ParMixedBilinearForm * m12EpsRe_;
+   ParMixedBilinearForm * m12EpsIm_;
+
+   ParDiscreteCurlOperator * curl_; // For Computing D from H
+   ParDiscreteLinearOperator * kReCross_;
+   ParDiscreteLinearOperator * kImCross_;
+
+   ParBilinearForm * m0_;
+   ParMixedBilinearForm * n20ZRe_;
+   ParMixedBilinearForm * n20ZIm_;
+   ParBilinearForm * m4r_;
+   ParBilinearForm * m4i_;
+   ParBilinearForm * m4er_;
+   ParBilinearForm * m4ei_;
+   ParBilinearForm * m4i1r_;
+   ParBilinearForm * m4i1i_;
+   ParBilinearForm * m4i2r_;
+   ParBilinearForm * m4i2i_;
+   ParBilinearForm * m4i3r_;
+   ParBilinearForm * m4i3i_;
+   
+   /*
+   ParBilinearForm * m4cr_;
+   ParBilinearForm * m4ci_;
+   
+   ParBilinearForm * m4solr_;
+   ParBilinearForm * m4soli_;
+   */
+
+   ParComplexGridFunction * e_;   // Complex electric field (HCurl)
+   ParComplexGridFunction * e_tmp_; // Temporary complex electric field (HCurl)
+   ParComplexGridFunction * d_;   // Complex electric flux (HDiv)
+   ParComplexGridFunction * b_;   // Complex magnetic flux (HDiv)
+
+   ParGridFunction * temp_; // Temporary grid function (HCurl)
+   ParDiscreteGradOperator * grad_; // For Computing E = Grad phi
+   ParDiscreteLinearOperator * kOpr_; // E += i k phi
+   ParDiscreteLinearOperator * kOpi_; // E += i (ik) phi
+   ParComplexGridFunction * phi_; // Complex sheath potential (H1)
+   ParComplexGridFunction * prev_phi_; // Complex sheath potential temporary (H1)
+   ParComplexGridFunction * next_phi_; // Complex sheath potential temporary (H1)
+   ParComplexGridFunction * z_; // Complex sheath potential (H1)
+   ParGridFunction        * power_absorp_t_; // Real valued total power absorption (H1)
+   ParGridFunction        * power_absorp_e_; // Real valued electron power absorption (H1)
+   ParGridFunction        * power_absorp_i1_; // Real valued ion 1 absorption (H1)
+   ParGridFunction        * power_absorp_i2_; // Real valued ion 2 power absorption (H1)
+   ParGridFunction        * power_absorp_i3_; // Real valued ion 3 power absorption (H1)
+
+   ParGridFunction * rectPot_; // Real valued rectified potential (H1)
+   ParComplexGridFunction * j_;   // Complex current density (HCurl)
+   ParComplexLinearForm   * rhs_; // Dual of complex current density (HCurl)
+   ParGridFunction        * e_t_; // Time dependent Electric field
+   ParComplexGridFunction * e_b_; // Complex parallel electric field (L2)
+   //ParComplexGridFunction * e_perp_; // Complex perpendicular electric field (L2)
+   ParComplexGridFunction * e_plus_; // Complex + polarized electric field (L2)
+   ParComplexGridFunction * e_min_; // Complex - polarized electric field (L2)
+   ParComplexGridFunction * e_v_; // Complex electric field (L2^d)
+   ParComplexGridFunction * d_v_; // Complex electric flux (L2^d)
+   ParComplexGridFunction * phi_v_; // Complex sheath potential (L2)
+   ParComplexGridFunction * j_v_; // Complex current density (L2^d)
+   ParGridFunction        * b_hat_; // Unit vector along B (HDiv)
+   ParGridFunction        * u_;   // Energy density (L2)
+   ParGridFunction        * uE_;  // Electric Energy density (L2)
+   ParGridFunction        * uB_;  // Magnetic Energy density (L2)
+   ParComplexGridFunction * S_;  // Poynting Vector (HDiv)
+   ParComplexGridFunction * StixS_; // Stix S Coefficient (L2)
+   ParComplexGridFunction * StixD_; // Stix D Coefficient (L2)
+   ParComplexGridFunction * StixP_; // Stix P Coefficient (L2)
+   //ParComplexGridFunction * EpsPara_; // B^T eps B / |B|^2 Coefficient (L2)
+
+   HypreParMatrix * M4r_;
+   HypreParMatrix * M4i_;
+   HypreParMatrix * M4er_;
+   HypreParMatrix * M4ei_;
+   HypreParMatrix * M4i1r_;
+   HypreParMatrix * M4i1i_;
+   HypreParMatrix * M4i2r_;
+   HypreParMatrix * M4i2i_;
+   HypreParMatrix * M4i3r_;
+   HypreParMatrix * M4i3i_;
+   
+   /*
+   HypreParMatrix * M4cr_;
+   HypreParMatrix * M4ci_;
+   HypreParMatrix * M4solr_;
+   HypreParMatrix * M4soli_;
+   */
+   mutable HypreParVector * RHSr1_;
+   mutable HypreParVector * RHSi1_;
+   mutable HypreParVector * RHSr2_;
+   mutable HypreParVector * RHSi2_;
+   mutable HypreParVector * RHSr3_;
+   mutable HypreParVector * RHSi3_;
+   mutable HypreParVector * RHSr4_;
+   mutable HypreParVector * RHSi4_;
+
+   mutable HypreParVector * RHSre_;
+   mutable HypreParVector * RHSie_;
+   mutable HypreParVector * RHSri1_;
+   mutable HypreParVector * RHSii1_;
+   mutable HypreParVector * RHSri2_;
+   mutable HypreParVector * RHSii2_;
+   mutable HypreParVector * RHSri3_;
+   mutable HypreParVector * RHSii3_;
+
+   mutable HypreParVector * TMPr2_;
+   mutable HypreParVector * TMPi2_;
+   mutable HypreParVector * TMPr3_;
+   mutable HypreParVector * TMPi3_;
+   mutable HypreParVector * TMPr4_;
+   mutable HypreParVector * TMPi4_;
+
+   mutable HypreParVector * TMPre_;
+   mutable HypreParVector * TMPie_;
+   mutable HypreParVector * TMPri1_;
+   mutable HypreParVector * TMPii1_;
+   mutable HypreParVector * TMPri2_;
+   mutable HypreParVector * TMPii2_;
+   mutable HypreParVector * TMPri3_;
+   mutable HypreParVector * TMPii3_;
+
+   HypreParVector * Er_; 
+   HypreParVector * Ei_;  
+
+   VectorCoefficient * BCoef_;        // B Field Unit Vector
+   MatrixCoefficient * epsReCoef_;    // Dielectric Material Coefficient
+   MatrixCoefficient * epsImCoef_;    // Dielectric Material Coefficient
+   MatrixCoefficient * epsAbsCoef_;   // Dielectric Material Coefficient
+   MatrixCoefficient * susceptReCoef_;    // Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_;    // Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_e_;    // Electron Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_e_;    // Electron Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i1_;    // Ion 1 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i1_;    // Ion 1 Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i2_;    // Ion 2 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i2_;    // Ion 2 Imag Susceptibility Coefficient
+   MatrixCoefficient * susceptReCoef_i3_;    // Ion 3 Real Susceptibility Coefficient
+   MatrixCoefficient * susceptImCoef_i3_;    // Ion 3 Imag Susceptibility Coefficient
+   MatrixCoefficient * muInvReCoef_;    // Real Dia/Paramagnetic Material Coefficient
+   MatrixCoefficient * muInvImCoef_;    // Imag Dia/Paramagnetic Material Coefficient
+   Coefficient       * etaInvCoef_;   // Admittance Coefficient
+   VectorCoefficient * kReCoef_;        // Wave Vector
+   VectorCoefficient * kImCoef_;        // Wave Vector
+
+   Coefficient * SReCoef_; // Stix S Coefficient
+   Coefficient * SImCoef_; // Stix S Coefficient
+   Coefficient * DReCoef_; // Stix D Coefficient
+   Coefficient * DImCoef_; // Stix D Coefficient
+   Coefficient * PReCoef_; // Stix P Coefficient
+   Coefficient * PImCoef_; // Stix P 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 * abcCoef_;       // -omega eta^{-1}
+   // Coefficient * sbcReCoef_;     //  omega Im(eta^{-1})
+   // Coefficient * sbcImCoef_;     // -omega Re(eta^{-1})
+   Coefficient * sinkx_;         // sin(ky * y + kz * z)
+   Coefficient * coskx_;         // cos(ky * y + kz * z)
+   Coefficient * negsinkx_;      // -sin(ky * y + kz * z)
+   // Coefficient * negMuInvCoef_;  // -1.0 / mu
+
+   MatrixCoefficient * massReCoef_;  // -omega^2 Re(epsilon)
+   MatrixCoefficient * massImCoef_;  // omega^2 Im(epsilon)
+   MatrixCoefficient * posMassCoef_; // omega^2 Abs(epsilon)
+   // MatrixCoefficient * negMuInvkxkxCoef_; // -\vec{k}\times\vec{k}\times/mu
+
+   kmkCoefficient kmkReCoef_;
+   kmkCoefficient kmkImCoef_;
+   CrossCoefficient kmReCoef_;
+   CrossCoefficient kmImCoef_;
+
+   VectorCoefficient * jrCoef_;     // Volume Current Density Function
+   VectorCoefficient * jiCoef_;     // Volume Current Density Function
+   VectorCoefficient * rhsrCoef_;     // Volume Current Density Function
+   VectorCoefficient * rhsiCoef_;     // Volume Current Density Function
+
+   VectorGridFunctionCoefficient erCoef_;
+   VectorGridFunctionCoefficient eiCoef_;
+
+   CurlGridFunctionCoefficient derCoef_;
+   CurlGridFunctionCoefficient deiCoef_;
+
+   EnergyDensityCoef     uCoef_;
+   ElectricEnergyDensityCoef uECoef_;
+   MagneticEnergyDensityCoef uBCoef_;
+   PoyntingVectorReCoef SrCoef_;
+   PoyntingVectorImCoef SiCoef_;
+
+   // const VectorCoefficient & erCoef_;     // Electric Field Boundary Condition
+   // const 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_bdr_marker_;
+
+   // Array of 0's and 1's marking the location of sheath surfaces
+   // Array<int> sbc_marker_;
+
+   // Array of 0's and 1's marking the location of Dirichlet boundaries
+   Array<int> dbc_bdr_marker_;
+   // void   (*e_r_bc_)(const Vector&, Vector&);
+   // void   (*e_i_bc_)(const Vector&, Vector&);
+
+   // Array<int> * dbcs_;
+   Array<ComplexVectorCoefficientByAttr*> * dbcs_;
+   Array<int> ess_bdr_;
+   Array<int> ess_bdr_tdofs_;
+   Array<int> non_k_bdr_;
+   Array<int> core_attr_marker_;
+   Array<int> sol_attr_marker_;
+
+   Array<ComplexVectorCoefficientByAttr*> * nbcs_; // Surface current BCs
+   Array<ComplexVectorCoefficientByAttr*> * nkbcs_; // Neumann BCs (-i*omega*K)
+
+   Array<ComplexCoefficientByAttr*> * sbcs_; // Sheath BCs
+
+   VisItDataCollection * visit_dc_;
+
+   std::map<std::string,socketstream*> socks_;
+};
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif // MFEM_COLD_PLASMA_DIELECTRIC_SOLVER

miniapps/plasma/g_eqdsk_data.cpp

@@ -0,0 +1,913 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "g_eqdsk_data.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+G_EQDSK_Data::G_EQDSK_Data(istream &is)
+   : init_flag_(0)
+{
+   double XDUM = 0.0;
+
+   const int buflen = 1024;
+   char buf[buflen];
+   is.getline(buf, buflen);
+   istringstream iss(buf);
+   string word;
+   iss >> std::ws;
+   while (!iss.eof())
+   {
+      iss >> word;
+      CASE_.push_back(word);
+      iss >> std::ws;
+   }
+
+   NW_ = to_int(CASE_[CASE_.size()-2]);
+   NH_ = to_int(CASE_[CASE_.size()-1]);
+
+   is >> RDIM_ >> ZDIM_ >> RCENTR_ >> RLEFT_ >> ZMID_;
+   is >> RMAXIS_ >> ZMAXIS_ >> SIMAG_ >> SIBRY_ >> BCENTR_;
+   is >> CURRENT_ >> SIMAG_ >> XDUM >> RMAXIS_ >> XDUM;
+   is >> ZMAXIS_ >> XDUM >> SIBRY_ >> XDUM >> XDUM;
+
+   FPOL_.resize(NW_);
+   PRES_.resize(NW_);
+   FFPRIM_.resize(NW_);
+   PPRIME_.resize(NW_);
+   PSIRZ_.resize(NW_ * NH_);
+   QPSI_.resize(NW_);
+
+   for (int i=0; i<NW_; i++) { is >> FPOL_[i]; }
+   for (int i=0; i<NW_; i++) { is >> PRES_[i]; }
+   for (int i=0; i<NW_; i++) { is >> FFPRIM_[i]; }
+   for (int i=0; i<NW_; i++) { is >> PPRIME_[i]; }
+   for (int j=0; j<NH_; j++)
+   {
+      for (int i=0; i<NW_; i++)
+      {
+         is >> PSIRZ_[NH_ * i + j];
+      }
+   }
+   for (int i=0; i<NW_; i++) { is >> QPSI_[i]; }
+
+   is >> NBBBS_ >> LIMITR_;
+
+   RBBBS_.resize(NBBBS_);
+   ZBBBS_.resize(NBBBS_);
+   RLIM_.resize(LIMITR_);
+   ZLIM_.resize(LIMITR_);
+
+   for (int i=0; i<NBBBS_; i++) { is >> RBBBS_[i] >> ZBBBS_[i]; }
+   for (int i=0; i<LIMITR_; i++) { is >> RLIM_[i] >> ZLIM_[i]; }
+
+   dr_ = RDIM_ / (NW_ - 1);
+   dz_ = ZDIM_ / (NH_ - 1);
+
+   double psi_bry = checkPsiBoundary();
+   if ((SIBRY_ - SIMAG_) < 1e-2 * (psi_bry - SIMAG_)) { SIBRY_ = psi_bry; }
+
+   dpsi_ = (SIBRY_ - SIMAG_) / (NW_ - 1);
+}
+
+void G_EQDSK_Data::PrintInfo(ostream & out) const
+{
+   out << endl << "G EQDSK File Info:" << endl;
+   out << "Size of grid: " << NW_ << " x " << NH_ << endl;
+   out << "Number of boundary points: " << NBBBS_ << endl;
+   out << "Number of limiter points:  " << LIMITR_ << endl;
+   out << endl;
+   out << "Range of R: " << RLEFT_ << " -> " << RLEFT_ + RDIM_ << endl;
+   out << "Range of Z: " << ZMID_ - 0.5 * ZDIM_
+       << " -> " << ZMID_ + 0.5 * ZDIM_ << endl;
+   out << "Location of magnetic axis: "
+       << "(" << RMAXIS_ << "," << ZMAXIS_ << ")" << endl;
+   out << "Poloidal flux at magnetic axis:   " << SIMAG_ << endl;
+   out << "Poloidal flux at plasma boundary: " << SIBRY_ << endl;
+   out << "R in meter of vacuum toroidal magnetic field BCENTR: "
+       << RCENTR_ << endl;
+   out << "Vacuum toroidal magnetic field in Tesla at RCENTR: "
+       << BCENTR_ << endl;
+   out << "Plasma current in Ampere: " << CURRENT_ << endl << endl;
+}
+
+void G_EQDSK_Data::DumpGnuPlotData(const string &file) const
+{
+   ostringstream oss_dat, oss_inp;
+   oss_inp << file << ".inp";
+   oss_dat << file << ".dat";
+   ofstream ofs_inp(oss_inp.str().c_str());
+   ofstream ofs_dat(oss_dat.str().c_str());
+
+   for (int i=0; i<NW_; i++)
+   {
+      ofs_dat << RLEFT_ + RDIM_ * i / (NW_ - 1)
+              << '\t' << FPOL_[i]
+              << '\t' << PRES_[i]
+              << '\t' << FFPRIM_[i]
+              << '\t' << PPRIME_[i]
+              << '\t' << QPSI_[i]
+              << '\n';
+   }
+
+   ofs_dat << "\n\n";
+   for (int i=0; i<NW_; i++)
+   {
+      for (int j=0; j<NH_; j++)
+      {
+         ofs_dat << RLEFT_ + RDIM_ * i / (NW_ - 1)
+                 << '\t' << ZMID_ - 0.5 * ZDIM_ + ZDIM_ * j / (NH_ - 1)
+                 << '\t' << PSIRZ_[NH_ * i + j]
+                 //      << '\t' << PSIRZ_[NH_ * i + j]
+                 << '\n';
+      }
+      ofs_dat << '\n';
+   }
+   ofs_dat << "\n\n";
+   for (int i=0; i<NBBBS_; i++)
+   {
+      ofs_dat << RBBBS_[i] << '\t' << ZBBBS_[i] << '\n';
+   }
+   ofs_dat << "\n\n";
+   for (int i=0; i<LIMITR_; i++)
+   {
+      ofs_dat << RLIM_[i] << '\t' << ZLIM_[i] << '\n';
+   }
+   ofs_dat.close();
+
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 0 using 1:2 w l t 'FPOL';\n";
+   ofs_inp << "set size noratio 1,1;\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 0 using 1:3 w l t 'PRES';\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 0 using 1:4 w l t 'FFPRIME';\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 0 using 1:5 w l t 'PPRIME';\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 0 using 1:6 w l t 'QPSI';\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "set view map;\n";
+   ofs_inp << "unset surface;\n";
+   ofs_inp << "set contour base;\n";
+   ofs_inp << "set cntrparam levels 20;\n";
+   ofs_inp << "set size ratio -1;\n";
+   ofs_inp << "set nokey;\n";
+   ofs_inp << "splot '" << oss_dat.str()
+           << "' index 1 with lines pal t 'PSIRZ';\n";
+   ofs_inp << "set key;\n";
+   ofs_inp << "pause -1;\n";
+   ofs_inp << "set size ratio -1;\n";
+   ofs_inp << "plot '" << oss_dat.str()
+           << "' index 2 using 1:2 w l t 'BOUNDARY',";
+   ofs_inp << " '" << oss_dat.str()
+           << "' index 3 using 1:2 w l t 'LIMITER';\n";
+
+   ofs_inp.close();
+}
+
+double G_EQDSK_Data::checkPsiBoundary()
+{
+   double psi_mid = 0.0;
+   double psi_min = DBL_MAX;
+   double psi_max = DBL_MIN;
+
+   Vector rz(2);
+   double psi = 0.0;
+   for (int i=0; i<NBBBS_; i++)
+   {
+      rz[0] = RBBBS_[i];
+      rz[1] = ZBBBS_[i];
+      psi = this->InterpPsiRZ(rz);
+
+      psi_min = std::min(psi, psi_min);
+      psi_max = std::max(psi, psi_max);
+
+      if (NBBBS_ % 2 == 1)
+      {
+         if (i == (NBBBS_ - 1) / 2) { psi_mid = psi; }
+      }
+      else
+      {
+         if (i == NBBBS_ / 2 || i + 1 == NBBBS_ / 2) { psi_mid += 0.5 * psi; }
+      }
+   }
+   //cout << psi_min << " <= psi <= " << psi_max << endl;
+   //cout << "psi_mid = " << psi_mid << endl;
+
+   return psi_mid;
+}
+/*
+double G_EQDSK_Data::InterpFPol(double r)
+{
+   if (!checkFlag(FPOL))
+   {
+      initInterpR(FPOL_, FPOL_t_);
+      setFlag(FPOL);
+   }
+   return interpR(r, FPOL_, FPOL_t_);
+}
+double G_EQDSK_Data::InterpPres(double r)
+{
+   if (!checkFlag(PRES))
+   {
+      initInterpR(PRES_, PRES_t_);
+      setFlag(PRES);
+   }
+   return interpR(r, PRES_, PRES_t_);
+}
+double G_EQDSK_Data::InterpFFPrime(double r)
+{
+   if (!checkFlag(FFPRIM))
+   {
+      initInterpR(FFPRIM_, FFPRIM_t_);
+      setFlag(FFPRIM);
+   }
+   return interpR(r, FFPRIM_, FFPRIM_t_);
+}
+double G_EQDSK_Data::InterpPPrime(double r)
+{
+   if (!checkFlag(PPRIME))
+   {
+      initInterpR(PPRIME_, PPRIME_t_);
+      setFlag(PPRIME);
+   }
+   return interpR(r, PPRIME_, PPRIME_t_);
+}
+double G_EQDSK_Data::InterpQPsi(double r)
+{
+   if (!checkFlag(QPSI))
+   {
+      initInterpR(QPSI_, QPSI_t_);
+      setFlag(QPSI);
+   }
+   return interpR(r, QPSI_, QPSI_t_);
+}
+double G_EQDSK_Data::InterpBTor(double r)
+{
+   if (!checkFlag(BTOR))
+   {
+      initInterpR(BTOR_, BTOR_t_);
+      setFlag(BTOR);
+   }
+   return interpR(r, BTOR_, BTOR_t_);
+}
+*/
+double G_EQDSK_Data::InterpFPolRZ(const Vector &rz)
+{
+   double psi = InterpPsiRZ(rz);
+
+   if (!checkFlag(FPOL))
+   {
+      initInterpPsi(FPOL_, FPOL_t_);
+      setFlag(FPOL);
+   }
+
+   return interpPsi(psi, FPOL_, FPOL_t_);
+}
+
+double G_EQDSK_Data::InterpPresRZ(const Vector &rz)
+{
+   double psi = InterpPsiRZ(rz);
+
+   if (!checkFlag(PRES))
+   {
+      initInterpPsi(PRES_, PRES_t_);
+      setFlag(PRES);
+   }
+
+   return interpPsi(psi, PRES_, PRES_t_);
+}
+
+double G_EQDSK_Data::InterpFFPrimeRZ(const Vector &rz)
+{
+   double psi = InterpPsiRZ(rz);
+
+   if (!checkFlag(FFPRIM))
+   {
+      initInterpPsi(FFPRIM_, FFPRIM_t_);
+      setFlag(FFPRIM);
+   }
+   return interpPsi(psi, FFPRIM_, FFPRIM_t_);
+}
+
+double G_EQDSK_Data::InterpPPrimeRZ(const Vector &rz)
+{
+   double psi = InterpPsiRZ(rz);
+
+   if (!checkFlag(PPRIME))
+   {
+      initInterpPsi(PPRIME_, PPRIME_t_);
+      setFlag(PPRIME);
+   }
+   return interpPsi(psi, PPRIME_, PPRIME_t_);
+}
+
+double G_EQDSK_Data::InterpPsiRZ(const Vector &rz)
+{
+   if (!checkFlag(PSIRZ))
+   {
+      initInterpRZ(PSIRZ_, PSIRZ_c_, PSIRZ_d_, PSIRZ_e_);
+      setFlag(PSIRZ);
+   }
+   return interpRZ(rz, PSIRZ_, PSIRZ_c_, PSIRZ_d_, PSIRZ_e_);
+}
+
+double G_EQDSK_Data::InterpQRZ(const Vector &rz)
+{
+   double psi = InterpPsiRZ(rz);
+
+   if (!checkFlag(QPSI))
+   {
+      initInterpPsi(QPSI_, QPSI_t_);
+      setFlag(QPSI);
+   }
+
+   return interpPsi(psi, QPSI_, QPSI_t_);
+}
+
+void G_EQDSK_Data::InterpNxGradPsiRZ(const Vector &rz, Vector &nxdp)
+{
+   if (!checkFlag(PSIRZ))
+   {
+      initInterpRZ(PSIRZ_, PSIRZ_c_, PSIRZ_d_, PSIRZ_e_);
+      setFlag(PSIRZ);
+   }
+   interpNxGradRZ(rz, PSIRZ_, PSIRZ_c_, PSIRZ_d_, PSIRZ_e_, nxdp);
+}
+
+void G_EQDSK_Data::InterpBPolRZ(const Vector &rz, Vector &bpol)
+{
+   InterpNxGradPsiRZ(rz, bpol);
+   if (rz[0] > 1e-6 * RDIM_) { bpol /= rz[0]; }
+}
+
+double G_EQDSK_Data::InterpBTorRZ(const Vector &rz)
+{
+   if (rz[0] > 1e-6 * RDIM_)
+   {
+      return InterpFPolRZ(rz) / rz[0];
+   }
+   else
+   {
+      return 0.0;
+   }
+}
+
+double G_EQDSK_Data::InterpJTorRZ(const Vector &rz)
+{
+   if (rz[0] > 1e-6 * RDIM_)
+   {
+      return InterpPPrimeRZ(rz) * rz[0] + InterpFFPrimeRZ(rz) / rz[0];
+   }
+   else
+   {
+      return 0.0;
+   }
+}
+
+void G_EQDSK_Data::initInterpR(const std::vector<double> &v,
+                               std::vector<double> &t)
+{
+   // Initialize the divided differences
+   ShiftedVector m(NW_-1, 2); m = 0.0;
+
+   m(-2) = -2.0 * v[2] + 5.0 * v[1] - 3.0 * v[0];
+   m(-1) = -1.0 * v[2] + 3.0 * v[1] - 2.0 * v[0];
+   for (int i=0; i<NW_-1; i++)
+   {
+      m(i) = v[i+1] - v[i];
+   }
+   m(NW_-1) = 2.0 * v[NW_-1] - 3.0 * v[NW_-2] + v[NW_-3];
+   m(NW_)   = 3.0 * v[NW_-1] - 5.0 * v[NW_-2] + 2.0 * v[NW_-3];
+
+   // Initialize the Slopes
+   t.resize(NW_);
+
+   for (int i=0; i<NW_; i++)
+   {
+      if (m(i+1) == m(i) && m(i-1) == m(i-2))
+      {
+         if (m(i) == m(i-1))
+         {
+            t[i] = m(i) * dr_;
+         }
+         else
+         {
+            t[i] = 0.5 * (m(i-1) + m(i)) * dr_;
+         }
+      }
+      else
+      {
+         t[i] = (fabs(m(i+1) - m(i)) * m(i-1) +
+                 fabs(m(i-1) - m(i-2)) * m(i)) * dr_ /
+                (fabs(m(i+1) - m(i)) + fabs(m(i-1) - m(i-2)));
+      }
+   }
+}
+
+double G_EQDSK_Data::interpR(double r, const vector<double> &v,
+                             const vector<double> &t)
+{
+   double rs = (r - RLEFT_) / RDIM_;
+
+   int i = std::max(0, std::min((int)floor(double(NW_-1) * rs), NW_-2));
+
+   // Compute ends of local patch
+   double r0 = RLEFT_ + RDIM_ * i / (NW_ - 1);
+   double r1 = r0 + RDIM_ / (NW_ - 1);
+
+   // Prepare position dependent factors
+   double wra = (r1 - r) / dr_;
+   double wrb = (r - r0) / dr_;
+   double wrc = (1.0 + 2.0 * wra);
+   double wrd = (1.0 + 2.0 * wrb);
+   double wra2 = wra * wra;
+   double wrb2 = wrb * wrb;
+
+   // Extract variable values at ends of local patch
+   const double &p0 = v[i];
+   const double &p1 = v[i+1];
+
+   double var = p0 * wra2 * wrd + p1 * wrb2 * wrc;
+
+   // Extract dvar/dx at ends of local patch
+   const double &px0 = t[i];
+   const double &px1 = t[i+1];
+
+   double varx = px0 * wra2 * wrb - px1 * wrb2 * wra;
+
+   var += varx * dr_;
+
+   return var;
+}
+
+void G_EQDSK_Data::initInterpRZ(const std::vector<double> &v,
+                                ShiftedDenseMatrix &c,
+                                ShiftedDenseMatrix &d,
+                                ShiftedDenseMatrix &e)
+{
+   ExtendedDenseMatrix ve(&v[0], NW_, NH_);
+
+   c.SetSize(NW_ + 3, NH_ + 2); c.SetShifts(2, 1); c = 0.0;
+   d.SetSize(NW_ + 2, NH_ + 3); d.SetShifts(1, 2); d = 0.0;
+   e.SetSize(NW_ + 1, NH_ + 1); e.SetShifts(1, 1); e = 0.0;
+
+   // x-directed divided differences
+   for (int i=-1; i<NW_; i++)
+   {
+      c(i,-1) = (ve(i+1,-1) - ve(i,-1)) / dr_;
+   }
+   for (int j=0; j<NH_; j++)
+   {
+      for (int i=-2; i<=NW_; i++)
+      {
+         c(i,j) = (ve(i+1,j) - ve(i,j)) / dr_;
+      }
+   }
+   for (int i=-1; i<NW_; i++)
+   {
+      c(i,NH_) = (ve(i+1,NH_) - ve(i,NH_)) / dr_;
+   }
+
+   // y-directed divided differences
+   for (int j=-1; j<NH_; j++)
+   {
+      d(-1,j) = (ve(-1,j+1) - ve(-1,j)) / dz_;
+   }
+   for (int i=0; i<NW_; i++)
+   {
+      for (int j=-2; j<=NH_; j++)
+      {
+         d(i,j) = (ve(i,j+1) - ve(i,j)) / dz_;
+      }
+   }
+   for (int j=-1; j<NH_; j++)
+   {
+      d(NW_,j) = (ve(NW_,j+1) - ve(NW_,j)) / dz_;
+   }
+
+   // Second order divided differences
+   for (int i=-1; i<NW_; i++)
+   {
+      for (int j=-1; j<NH_; j++)
+      {
+         e(i,j) = (c(i,j+1) - c(i,j)) / dz_;
+      }
+   }
+}
+
+double G_EQDSK_Data::interpRZ(const Vector &rz,
+                              const std::vector<double> &v,
+                              const ShiftedDenseMatrix &c,
+                              const ShiftedDenseMatrix &d,
+                              const ShiftedDenseMatrix &e)
+{
+   double r = rz[0];
+   double z = rz[1];
+
+   double rs = (r - RLEFT_) / RDIM_;
+   double zs = (z - ZMID_ + 0.5 * ZDIM_) / ZDIM_;
+
+   int i = std::max(0, std::min((int)floor(double(NW_-1) * rs), NW_-2));
+   int j = std::max(0, std::min((int)floor(double(NH_-1) * zs), NH_-2));
+
+   // Compute corners of local patch
+   double r0 = RLEFT_ + RDIM_ * i / (NW_ - 1);
+   double r1 = r0 + RDIM_ / (NW_ - 1);
+   double z0 = ZMID_ - 0.5 * ZDIM_ + ZDIM_ * j / (NH_ - 1);
+   double z1 = z0 + ZDIM_ / (NH_ - 1);
+
+   // Prepare position dependent factors
+   double wra = (r1 - r) / dr_;
+   double wrb = (r - r0) / dr_;
+   double wrc = (1.0 + 2.0 * wra);
+   double wrd = (1.0 + 2.0 * wrb);
+   double wra2 = wra * wra;
+   double wrb2 = wrb * wrb;
+
+   double wza = (z1 - z) / dz_;
+   double wzb = (z - z0) / dz_;
+   double wzc = (1.0 + 2.0 * wza);
+   double wzd = (1.0 + 2.0 * wzb);
+   double wza2 = wza * wza;
+   double wzb2 = wzb * wzb;
+
+   // Extract variable values at corners of local patch
+   double p00 = v[NH_ * i + j];
+   double p10 = v[NH_ * (i + 1) + j];
+   double p01 = v[NH_ * i + j + 1];
+   double p11 = v[NH_ * (i + 1) + j + 1];
+
+   double var = p00 * wra2 * wrd * wza2 * wzd
+                + p10 * wrb2 * wrc * wza2 * wzd
+                + p01 * wra2 * wrd * wzb2 * wzc
+                + p11 * wrb2 * wrc * wzb2 * wzc;
+
+   // Compute dvar/dx at corners of local patch
+   double wx00a = fabs(c(i-1,j) - c(i-2,j));
+   double wx00b = fabs(c(i+1,j) - c(i,j));
+
+   double wx10a = fabs(c(i,j) - c(i-1,j));
+   double wx10b = fabs(c(i+2,j) - c(i+1,j));
+
+   double wx01a = fabs(c(i-1,j+1) - c(i-2,j+1));
+   double wx01b = fabs(c(i+1,j+1) - c(i,j+1));
+
+   double wx11a = fabs(c(i,j+1) - c(i-1,j+1));
+   double wx11b = fabs(c(i+2,j+1) - c(i+1,j+1));
+
+   if (wx00a == 0.0 && wx00b == 0.0) { wx00a = 1.0; wx00b = 1.0; }
+   if (wx10a == 0.0 && wx10b == 0.0) { wx10a = 1.0; wx10b = 1.0; }
+   if (wx01a == 0.0 && wx01b == 0.0) { wx01a = 1.0; wx01b = 1.0; }
+   if (wx11a == 0.0 && wx11b == 0.0) { wx11a = 1.0; wx11b = 1.0; }
+
+   double px00 = (wx00b * c(i-1,j) + wx00a * c(i,j)) / (wx00b + wx00a);
+   double px10 = (wx10b * c(i,j) + wx10a * c(i+1,j)) / (wx10b + wx10a);
+   double px01 = (wx01b * c(i-1,j+1) + wx01a * c(i,j+1)) / (wx01b + wx01a);
+   double px11 = (wx11b * c(i,j+1) + wx11a * c(i+1,j+1)) / (wx11b + wx11a);
+
+   double varx = px00 * wra2 * wrb * wza2 * wzd
+                 - px10 * wrb2 * wra * wza2 * wzd
+                 + px01 * wrb * wra2 * wzb2 * wzc
+                 - px11 * wra * wrb2 * wzb2 * wzc;
+   var += varx * dr_;
+
+   // Compute dvar/dy at corners of local patch
+   double wy00a = fabs(d(i,j-1) - d(i,j-2));
+   double wy00b = fabs(d(i,j+1) - d(i,j));
+
+   double wy10a = fabs(d(i+1,j-1) - d(i+1,j-2));
+   double wy10b = fabs(d(i+1,j+1) - d(i+1,j));
+
+   double wy01a = fabs(d(i,j) - d(i,j-1));
+   double wy01b = fabs(d(i,j+2) - d(i,j+1));
+
+   double wy11a = fabs(d(i+1,j) - d(i+1,j-1));
+   double wy11b = fabs(d(i+1,j+2) - d(i+1,j+1));
+
+   if (wy00a == 0.0 && wy00b == 0.0) { wy00a = 1.0; wy00b = 1.0; }
+   if (wy10a == 0.0 && wy10b == 0.0) { wy10a = 1.0; wy10b = 1.0; }
+   if (wy01a == 0.0 && wy01b == 0.0) { wy01a = 1.0; wy01b = 1.0; }
+   if (wy11a == 0.0 && wy11b == 0.0) { wy11a = 1.0; wy11b = 1.0; }
+
+   double py00 = (wy00b * d(i,j-1) + wy00a * d(i,j)) / (wy00b + wy00a);
+   double py10 = (wy10b * d(i+1,j-1) + wy10a * d(i+1,j)) / (wy10b + wy10a);
+   double py01 = (wy01b * d(i,j) + wy01a * d(i,j+1)) / (wy01b + wy01a);
+   double py11 = (wy11b * d(i+1,j) + wy11a * d(i+1,j)) / (wy11b + wy11a);
+
+   double vary = py00 * wra2 * wrd * wza2 * wzb
+                 + py10 * wrb2 * wrc * wza2 * wzb
+                 - py01 * wra2 * wrd * wza * wzb2
+                 - py11 * wrb2 * wrc * wza * wzb2;
+   var += vary * dz_;
+
+   // Compute d^2var/dxdy at corners of local patch
+   double pxy00 = (wx00b * (wy00b * e(i-1,j-1) + wy00a * e(i-1,j)) +
+                   wx00a * (wy00b * e(i,j-1) + wy00a * e(i,j))) /
+                  ((wx00b + wx00a) * (wy00b + wy00a));
+   double pxy10 = (wx10b * (wy10b * e(i,j-1) + wy10a * e(i,j)) +
+                   wx10a * (wy10b * e(i+1,j-1) + wy10a * e(i+1,j))) /
+                  ((wx10b + wx10a) * (wy10b + wy10a));
+   double pxy01 = (wx01b * (wy01b * e(i-1,j) + wy01a * e(i-1,j+1)) +
+                   wx01a * (wy01b * e(i,j) + wy01a * e(i,j+1))) /
+                  ((wx01b + wx01a) * (wy01b + wy01a));
+   double pxy11 = (wx11b * (wy11b * e(i,j) + wy11a * e(i,j+1)) +
+                   wx11a * (wy11b * e(i+1,j) + wy11a * e(i+1,j+1))) /
+                  ((wx11b + wx11a) * (wy11b + wy11a));
+
+   double varxy = pxy00 * wra2 * wrb * wza2 * wzb
+                  - pxy10 * wra * wrb2 * wza2 * wzb
+                  - pxy01 * wra2 * wrb * wza * wzb2
+                  + pxy11 * wra * wrb2 * wza * wzb2;
+
+   var += dr_ * dz_ * varxy;
+
+   return var;
+}
+
+void G_EQDSK_Data::interpNxGradRZ(const Vector &rz,
+                                  const std::vector<double> &v,
+                                  const ShiftedDenseMatrix &c,
+                                  const ShiftedDenseMatrix &d,
+                                  const ShiftedDenseMatrix &e,
+                                  Vector &b)
+{
+   b.SetSize(2);
+   b = 0.0;
+
+   double r = rz[0];
+   double z = rz[1];
+
+   double rs = (r - RLEFT_) / RDIM_;
+   double zs = (z - ZMID_ + 0.5 * ZDIM_) / ZDIM_;
+
+   int i = std::max(0, std::min((int)floor(double(NW_-1) * rs), NW_-2));
+   int j = std::max(0, std::min((int)floor(double(NH_-1) * zs), NH_-2));
+
+   // Compute corners of local patch
+   double r0 = RLEFT_ + RDIM_ * i / (NW_ - 1);
+   double r1 = r0 + RDIM_ / (NW_ - 1);
+   double z0 = ZMID_ - 0.5 * ZDIM_ + ZDIM_ * j / (NH_ - 1);
+   double z1 = z0 + ZDIM_ / (NH_ - 1);
+
+   // Prepare position dependent factors
+   double wra = (r1 - r) / dr_, dwra = -1.0 / dr_;
+   double wrb = (r - r0) / dr_, dwrb = 1.0 / dr_;
+   double wrc = (1.0 + 2.0 * wra), dwrc = 2.0 * dwra;
+   double wrd = (1.0 + 2.0 * wrb), dwrd = 2.0 * dwrb;
+   double wra2 = wra * wra, dwra2 = 2.0 * wra * dwra;
+   double wrb2 = wrb * wrb, dwrb2 = 2.0 * wrb * dwrb;
+
+   double wza = (z1 - z) / dz_, dwza = -1.0 / dz_;
+   double wzb = (z - z0) / dz_, dwzb = 1.0 / dz_;
+   double wzc = (1.0 + 2.0 * wza), dwzc = 2.0 * dwza;
+   double wzd = (1.0 + 2.0 * wzb), dwzd = 2.0 * dwzb;
+   double wza2 = wza * wza, dwza2 = 2.0 * wza * dwza;
+   double wzb2 = wzb * wzb, dwzb2 = 2.0 * wzb * dwzb;
+
+   // Extract var values at corners of local patch
+   double p00 = v[NH_ * i + j];
+   double p10 = v[NH_ * (i + 1) + j];
+   double p01 = v[NH_ * i + j + 1];
+   double p11 = v[NH_ * (i + 1) + j + 1];
+
+   b[0] -=
+      (p00 * wra2 * wrd + p10 * wrb2 * wrc ) * (dwza2 * wzd + wza2 * dwzd)
+      + (p01 * wra2 * wrd + p11 * wrb2 * wrc) * (dwzb2 * wzc + wzb2 * dwzc);
+   b[1] +=
+      (p00 * wza2 * wzd + p01 * wzb2 * wzc) * (dwra2 * wrd + wra2 * dwrd)
+      + (p10 * wza2 * wzd + p11 * wzb2 * wzc) * (dwrb2 * wrc + wrb2 * dwrc);
+
+   // Compute dvar/dx at corners of local patch
+   double wx00a = fabs(c(i-1,j) - c(i-2,j));
+   double wx00b = fabs(c(i+1,j) - c(i,j));
+
+   double wx10a = fabs(c(i,j) - c(i-1,j));
+   double wx10b = fabs(c(i+2,j) - c(i+1,j));
+
+   double wx01a = fabs(c(i-1,j+1) - c(i-2,j+1));
+   double wx01b = fabs(c(i+1,j+1) - c(i,j+1));
+
+   double wx11a = fabs(c(i,j+1) - c(i-1,j+1));
+   double wx11b = fabs(c(i+2,j+1) - c(i+1,j+1));
+
+   if (wx00a == 0.0 && wx00b == 0.0) { wx00a = 1.0; wx00b = 1.0; }
+   if (wx10a == 0.0 && wx10b == 0.0) { wx10a = 1.0; wx10b = 1.0; }
+   if (wx01a == 0.0 && wx01b == 0.0) { wx01a = 1.0; wx01b = 1.0; }
+   if (wx11a == 0.0 && wx11b == 0.0) { wx11a = 1.0; wx11b = 1.0; }
+
+   double px00 = (wx00b * c(i-1,j) + wx00a * c(i,j)) / (wx00b + wx00a);
+   double px10 = (wx10b * c(i,j) + wx10a * c(i+1,j)) / (wx10b + wx10a);
+   double px01 = (wx01b * c(i-1,j+1) + wx01a * c(i,j+1)) / (wx01b + wx01a);
+   double px11 = (wx11b * c(i,j+1) + wx11a * c(i+1,j+1)) / (wx11b + wx11a);
+
+   b[0] -= dr_ *
+           ((px00 * wra2 * wrb - px10 * wrb2 * wra) *
+            (dwza2 * wzd + wza2 * dwzd) +
+            (px01 * wrb * wra2 - px11 * wra * wrb2) *
+            (dwzb2 * wzc + wzb2 * dwzc));
+
+   b[1] += dr_ *
+           ((px00 * wza2 * wzd + px01 * wzb2 * wzc) *
+            (dwra2 * wrb + wra2 * dwrb ) -
+            (px10 * wza2 * wzd + px11 * wzb2 * wzc) *
+            (dwra * wrb2 + wra * dwrb2));
+
+   // Compute dvar/dy at corners of local patch
+   double wy00a = fabs(d(i,j-1) - d(i,j-2));
+   double wy00b = fabs(d(i,j+1) - d(i,j));
+
+   double wy10a = fabs(d(i+1,j-1) - d(i+1,j-2));
+   double wy10b = fabs(d(i+1,j+1) - d(i+1,j));
+
+   double wy01a = fabs(d(i,j) - d(i,j-1));
+   double wy01b = fabs(d(i,j+2) - d(i,j+1));
+
+   double wy11a = fabs(d(i+1,j) - d(i+1,j-1));
+   double wy11b = fabs(d(i+1,j+2) - d(i+1,j+1));
+
+   if (wy00a == 0.0 && wy00b == 0.0) { wy00a = 1.0; wy00b = 1.0; }
+   if (wy10a == 0.0 && wy10b == 0.0) { wy10a = 1.0; wy10b = 1.0; }
+   if (wy01a == 0.0 && wy01b == 0.0) { wy01a = 1.0; wy01b = 1.0; }
+   if (wy11a == 0.0 && wy11b == 0.0) { wy11a = 1.0; wy11b = 1.0; }
+
+   double py00 = (wy00b * d(i,j-1) + wy00a * d(i,j)) / (wy00b + wy00a);
+   double py10 = (wy10b * d(i+1,j-1) + wy10a * d(i+1,j)) / (wy10b + wy10a);
+   double py01 = (wy01b * d(i,j) + wy01a * d(i,j+1)) / (wy01b + wy01a);
+   double py11 = (wy11b * d(i+1,j) + wy11a * d(i+1,j)) / (wy11b + wy11a);
+
+   b[0] -= dz_ *
+           ((py00 * wra2 * wrd + py10 * wrb2 * wrc) *
+            (dwza2 * wzb + wza2 * dwzb) -
+            (py01 * wra2 * wrd + py11 * wrb2 * wrc) *
+            (dwza * wzb2 + wza * dwzb2));
+   b[1] += dz_ *
+           ((py00 * wza2 * wzb - py01 * wza * wzb2) *
+            (dwra2 * wrd + wra2 * dwrd) +
+            (py10 * wza2 * wzb - py11 * wza * wzb2) *
+            (dwrb2 * wrc + wrb2 * dwrc));
+
+   // Compute d^2var/dxdy at corners of local patch
+   double pxy00 = (wx00b * (wy00b * e(i-1,j-1) + wy00a * e(i-1,j)) +
+                   wx00a * (wy00b * e(i,j-1) + wy00a * e(i,j))) /
+                  ((wx00b + wx00a) * (wy00b + wy00a));
+   double pxy10 = (wx10b * (wy10b * e(i,j-1) + wy10a * e(i,j)) +
+                   wx10a * (wy10b * e(i+1,j-1) + wy10a * e(i+1,j))) /
+                  ((wx10b + wx10a) * (wy10b + wy10a));
+   double pxy01 = (wx01b * (wy01b * e(i-1,j) + wy01a * e(i-1,j+1)) +
+                   wx01a * (wy01b * e(i,j) + wy01a * e(i,j+1))) /
+                  ((wx01b + wx01a) * (wy01b + wy01a));
+   double pxy11 = (wx11b * (wy11b * e(i,j) + wy11a * e(i,j+1)) +
+                   wx11a * (wy11b * e(i+1,j) + wy11a * e(i+1,j+1))) /
+                  ((wx11b + wx11a) * (wy11b + wy11a));
+
+   b[0] -= dr_ * dz_ * ((pxy00 * wra2 * wrb - pxy10 * wra * wrb2)
+                        * (dwza2 * wzb + wza2 * dwzb) +
+                        (pxy11 * wra * wrb2 - pxy01 * wra2 * wrb)
+                        * (dwza * wzb2 + wza * dwzb2));
+   b[1] += dr_ * dz_ * ((pxy00 * wza2 * wzb - pxy01 * wza * wzb2)
+                        * (dwra2 * wrb + wra2 * dwrb) +
+                        (pxy11 * wza * wzb2 - pxy10 * wza2 * wzb)
+                        * (dwra * wrb2 + wra * dwrb2));
+}
+
+void G_EQDSK_Data::initInterpPsi(const std::vector<double> &v,
+                                 std::vector<double> &t)
+{
+   // Initialize the divided differences
+   ShiftedVector m(NW_-1, 2); m = 0.0;
+
+   m(-2) = -2.0 * v[2] + 5.0 * v[1] - 3.0 * v[0];
+   m(-1) = -1.0 * v[2] + 3.0 * v[1] - 2.0 * v[0];
+   for (int i=0; i<NW_-1; i++)
+   {
+      m(i) = v[i+1] - v[i];
+   }
+   m(NW_-1) = 2.0 * v[NW_-1] - 3.0 * v[NW_-2] + v[NW_-3];
+   m(NW_)   = 3.0 * v[NW_-1] - 5.0 * v[NW_-2] + 2.0 * v[NW_-3];
+
+   // Initialize the Slopes
+   t.resize(NW_);
+
+   for (int i=0; i<NW_; i++)
+   {
+      if (m(i+1) == m(i) && m(i-1) == m(i-2))
+      {
+         if (m(i) == m(i-1))
+         {
+            t[i] = m(i) * dpsi_;
+         }
+         else
+         {
+            t[i] = 0.5 * (m(i-1) + m(i)) * dpsi_;
+         }
+      }
+      else
+      {
+         t[i] = (fabs(m(i+1) - m(i)) * m(i-1) +
+                 fabs(m(i-1) - m(i-2)) * m(i)) * dpsi_ /
+                (fabs(m(i+1) - m(i)) + fabs(m(i-1) - m(i-2)));
+      }
+   }
+}
+
+double G_EQDSK_Data::interpPsi(double psi, const vector<double> &v,
+                               const vector<double> &t)
+{
+   double psic = std::max(SIMAG_, std::min(psi, SIBRY_));
+
+   double psis = (psic - SIMAG_) / (SIBRY_ - SIMAG_);
+
+   int i = std::max(0, std::min((int)floor(double(NW_-1) * psis), NW_-2));
+
+   // Compute ends of local patch
+   double psi0 = SIMAG_ + (SIBRY_ - SIMAG_) * i / (NW_ - 1);
+   double psi1 = psi0 + (SIBRY_ - SIMAG_) / (NW_ - 1);
+
+   // Prepare position dependent factors
+   double wra = (psi1 - psic) / dpsi_;
+   double wrb = (psic - psi0) / dpsi_;
+   double wrc = (1.0 + 2.0 * wra);
+   double wrd = (1.0 + 2.0 * wrb);
+   double wra2 = wra * wra;
+   double wrb2 = wrb * wrb;
+
+   // Extract variable values at ends of local patch
+   const double &p0 = v[i];
+   const double &p1 = v[i+1];
+
+   double var = p0 * wra2 * wrd + p1 * wrb2 * wrc;
+
+   // Extract dvar/dx at ends of local patch
+   const double &px0 = t[i];
+   const double &px1 = t[i+1];
+
+   double varx = px0 * wra2 * wrb - px1 * wrb2 * wra;
+
+   var += varx * dpsi_;
+
+   return var;
+}
+
+void G_EQDSK_Data::ExtendedDenseMatrix::init()
+{
+   // Populate four corners
+   SW_ = 3.0 * ((*this)(0,0) - (*this)(1,1)) + (*this)(2,2);
+   SE_ = 3.0 * ((*this)(m_-1,0) - (*this)(m_-2,1)) + (*this)(m_-3,2);
+   NW_ = 3.0 * ((*this)(0,n_-1) - (*this)(1,n_-2)) + (*this)(2,n_-3);
+   NE_ = 3.0 * ((*this)(m_-1,n_-1) - (*this)(m_-2,n_-2))
+         + (*this)(m_-3,n_-3);
+
+   // Populate lowest rows
+   for (int j=0; j<n_; j++)
+   {
+      S_(1,j) = 3.0 * ((*this)(0,j) - (*this)(1,j)) + (*this)(2,j);
+      S_(0,j) = 3.0 * (2.0 * (*this)(0,j) + (*this)(2,j)) - 8.0 * (*this)(1,j);
+   }
+
+   // Populate highest rows
+   for (int j=0; j<n_; j++)
+   {
+      N_(1,j) = 3.0 * (2.0 * (*this)(m_-1,j) + (*this)(m_-3,j))
+                - 8.0 * (*this)(m_-2,j);
+      N_(0,j) = 3.0 * ((*this)(m_-1,j) - (*this)(m_-2,j)) + (*this)(m_-3,j);
+   }
+
+   // Populate lowest columns
+   for (int i=0; i<m_; i++)
+   {
+      W_(i,0) = 3.0 * (2.0 * (*this)(i,0) + (*this)(i,2)) - 8.0 * (*this)(i,1);
+      W_(i,1) = 3.0 * ((*this)(i,0) - (*this)(i,1)) + (*this)(i,2);
+   }
+
+   // Populate highest columns
+   for (int i=0; i<m_; i++)
+   {
+      E_(i,0) = 3.0 * ((*this)(i,n_-1) - (*this)(i,n_-2)) + (*this)(i,n_-3);
+      E_(i,1) = 3.0 * (2.0 * (*this)(i,n_-1) + (*this)(i,n_-3))
+                - 8.0 * (*this)(i,n_-2);
+   }
+}
+
+} // namespace plasma
+
+} // namespace mfem

miniapps/plasma/g_eqdsk_data.hpp

@@ -0,0 +1,506 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_G_EQDSK_DATA_HPP
+#define MFEM_G_EQDSK_DATA_HPP
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+#include "mfem.hpp"
+#include "../../general/text.hpp"
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+class G_EQDSK_Data
+{
+public:
+   G_EQDSK_Data(std::istream &is);
+
+   int GetNumPtsR() const { return NW_; }
+   int GetNumPtsZ() const { return NH_; }
+
+   double GetRExtent() const { return RDIM_; }
+   double GetZExtent() const { return ZDIM_; }
+
+   double GetRMin() const { return RLEFT_; }
+   double GetZMid() const { return ZMID_; }
+
+   double GetPsiCenter() const {return SIMAG_; }
+   double GetPsiBdry() const {return SIBRY_; }
+
+   std::vector<double> & GetPsi() { return PSIRZ_ ;}
+   // std::vector<double> & GetBTor() { return BTOR_; }
+
+   void PrintInfo(std::ostream &out = std::cout) const;
+   void DumpGnuPlotData(const std::string &file) const;
+
+   // double InterpFPol(double r);
+   // double InterpPres(double r);
+   // double InterpFFPrime(double r);
+   // double InterpPPrime(double r);
+   // double InterpQPsi(double r);
+   double InterpFPolRZ(const Vector &rz);
+   double InterpPresRZ(const Vector &rz);
+   double InterpFFPrimeRZ(const Vector &rz);
+   double InterpPPrimeRZ(const Vector &rz);
+   double InterpPsiRZ(const Vector &rz);
+   double InterpQRZ(const Vector &rz);
+   double InterpBTorRZ(const Vector &rz);
+   double InterpJTorRZ(const Vector &rz);
+
+   void InterpNxGradPsiRZ(const Vector &rz, Vector &nxdp);
+   void InterpBPolRZ(const Vector &rz, Vector &b);
+   // double InterpBTor(double r);
+
+   int GetNumBoundaryPts() const { return NBBBS_; }
+   const std::vector<double> & GetBoundaryRVals() const { return RBBBS_; }
+   const std::vector<double> & GetBoundaryZVals() const { return ZBBBS_; }
+
+   int GetNumLimiterPts() const { return LIMITR_; }
+   const std::vector<double> & GetLimiterRVals() const { return RLIM_; }
+   const std::vector<double> & GetLimiterZVals() const { return ZLIM_; }
+
+private:
+   class ShiftedVector;
+   class ShiftedDenseMatrix;
+   class ExtendedDenseMatrix;
+
+   enum FieldType {FPOL, PRES, FFPRIM, PPRIME, PSIRZ, QPSI/*, BTOR*/};
+
+   int init_flag_;
+   inline bool checkFlag(int flag) { return (init_flag_ >> flag) & 1; }
+   inline void   setFlag(int flag) { init_flag_ |= (1 << flag); }
+   inline void clearFlag(int flag) { init_flag_ &= ~(1 << flag); }
+
+   double checkPsiBoundary();
+
+   void initInterpR(const std::vector<double> &v,
+                    std::vector<double> &t);
+   void initInterpPsi(const std::vector<double> &v,
+                      std::vector<double> &t);
+   void initInterpRZ(const std::vector<double> &v,
+                     ShiftedDenseMatrix &c,
+                     ShiftedDenseMatrix &d,
+                     ShiftedDenseMatrix &e);
+
+   double interpR(double r, const std::vector<double> &v,
+                  const std::vector<double> &t);
+   double interpRZ(const Vector &rz,
+                   const std::vector<double> &v,
+                   const ShiftedDenseMatrix &c,
+                   const ShiftedDenseMatrix &d,
+                   const ShiftedDenseMatrix &e);
+   void interpNxGradRZ(const Vector &rz,
+                       const std::vector<double> &v,
+                       const ShiftedDenseMatrix &c,
+                       const ShiftedDenseMatrix &d,
+                       const ShiftedDenseMatrix &e,
+                       Vector &b);
+   double interpPsi(double psi, const std::vector<double> &v,
+                    const std::vector<double> &t);
+
+   std::vector<std::string> CASE_; // Identification character string
+
+   int NW_; // Number of horizontal R grid points
+   int NH_; // Number of vertical Z grid points
+
+   double RDIM_;    // Horizontal dimension in meter of computational box
+   double ZDIM_;    // Vertical dimension in meter of computational box
+   double RLEFT_;   // Minimum R in meter of rectangular computational box
+   double ZMID_;    // Z of center of computational box in meter
+   double RMAXIS_;  // R of magnetic axis in meter
+   double ZMAXIS_;  // Z of magnetic axis in meter
+   double SIMAG_;   // poloidal flux at magnetic axis in Weber /rad
+   double SIBRY_;   // poloidal flux at the plasma boundary in Weber /rad
+   double RCENTR_;  // R in meter of vacuum toroidal magnetic field BCENTR
+   double BCENTR_;  // Vacuum toroidal magnetic field in Tesla at RCENTR
+   double CURRENT_; // Plasma current in Ampere
+
+   // Poloidal current function in m-T, F = RBT on flux grid
+   std::vector<double> FPOL_;
+
+   // Plasma pressure in nt / m^2 on uniform flux grid
+   std::vector<double> PRES_;
+
+   // FF’(ψ) in (mT)2 / (Weber /rad) on uniform flux grid
+   std::vector<double> FFPRIM_;
+
+   // P’(ψ) in (nt /m2) / (Weber /rad) on uniform flux grid
+   std::vector<double> PPRIME_;
+
+   // Poloidal flux in Weber / rad on the rectangular grid points
+   std::vector<double> PSIRZ_;
+
+   // q values on uniform flux grid from axis to boundary
+   std::vector<double> QPSI_;
+
+   // Toroidal B field dervided from FPOL_
+   // std::vector<double> BTOR_;
+
+   int                 NBBBS_;  // Number of boundary points
+   std::vector<double> RBBBS_;  // R of boundary points in meter
+   std::vector<double> ZBBBS_;  // Z of boundary points in meter
+
+   int                 LIMITR_; // Number of limiter points
+   std::vector<double> RLIM_;   // R of surrounding limiter contour in meter
+   std::vector<double> ZLIM_;   // Z of surrounding limiter contour in meter
+
+   class ShiftedVector : public Vector
+   {
+   private:
+      int si_;
+   public:
+      ShiftedVector()
+         : si_(0) {}
+
+      ShiftedVector(int s, int si)
+         : Vector(s+2*si), si_(si) {}
+
+      void SetShift(int si) { si_ = si; }
+
+      ShiftedVector &operator=(double c)
+      { Vector::operator=(c); return *this; }
+
+      inline double &operator()(int i)
+      { return Vector::operator()(i + si_); }
+
+      inline const double &operator()(int i) const
+      { return Vector::operator()(i + si_); }
+   };
+
+   class ShiftedDenseMatrix : public DenseMatrix
+   {
+   private:
+      int si_, sj_;
+   public:
+      ShiftedDenseMatrix()
+         : si_(0), sj_(0) {}
+
+      ShiftedDenseMatrix(int m, int n, int si, int sj)
+         : DenseMatrix(m+2*si, n+2*sj), si_(si), sj_(sj) {}
+
+      void SetShifts(int si, int sj) { si_ = si; sj_ = sj; }
+
+      ShiftedDenseMatrix &operator=(double c)
+      { DenseMatrix::operator=(c); return *this; }
+
+      inline double &operator()(int i, int j)
+      { return DenseMatrix::operator()(i + si_, j + sj_); }
+
+      inline const double &operator()(int i, int j) const
+      { return DenseMatrix::operator()(i + si_, j + sj_); }
+   };
+
+   class ExtendedDenseMatrix
+   {
+   private:
+      int m_, n_;
+      const double *C_;
+      DenseMatrix N_;
+      DenseMatrix S_;
+      DenseMatrix E_;
+      DenseMatrix W_;
+      double SW_, SE_, NW_, NE_, DUMMY_;
+
+      void init();
+
+   public:
+      ExtendedDenseMatrix(const double *C, int m, int n)
+         : m_(m), n_(n), C_(C),
+           N_(2, n), S_(2, n),
+           E_(m, 2), W_(m, 2),
+           SW_(0.0), SE_(0.0), NW_(0.0), NE_(0.0), DUMMY_(0.0)
+      { N_ = 0.0; S_ = 0.0; E_ = 0.0; W_ = 0.0; init(); }
+
+      const double &operator()(int i, int j) const
+      {
+         if (i >= 0 && i < m_ && j >= 0 && j < n_)
+         {
+            return C_[n_ * i + j];
+         }
+         else if (i >= 0 && i < m_)
+         {
+            if (j < 0)
+            {
+               return W_(i, j + 2);
+            }
+            else
+            {
+               return E_(i, j - n_);
+            }
+         }
+         else if (j >= 0 && j < n_)
+         {
+            if (i < 0)
+            {
+               return S_(i + 2, j);
+            }
+            else
+            {
+               return N_(i - m_, j);
+            }
+         }
+         else if (i == -1 && j == -1)
+         {
+            return SW_;
+         }
+         else if (i == -1 && j == n_)
+         {
+            return SE_;
+         }
+         else if (i == m_ && j == -1)
+         {
+            return NW_;
+         }
+         else if (i == m_ && j == n_)
+         {
+            return NE_;
+         }
+         return DUMMY_;
+      }
+   };
+
+   // Divided differences for Akima's interpolation method
+   double dr_, dz_, dpsi_;
+
+   std::vector<double> FPOL_t_;
+   std::vector<double> PRES_t_;
+   std::vector<double> FFPRIM_t_;
+   std::vector<double> PPRIME_t_;
+   ShiftedDenseMatrix  PSIRZ_c_;
+   ShiftedDenseMatrix  PSIRZ_d_;
+   ShiftedDenseMatrix  PSIRZ_e_;
+   std::vector<double> QPSI_t_;
+
+   // std::vector<double> BTOR_t_;
+};
+
+class G_EQDSK_Psi_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_Psi_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpPsiRZ(transip);
+   }
+};
+
+class G_EQDSK_FPol_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_FPol_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpFPolRZ(transip);
+   }
+};
+
+class G_EQDSK_Pres_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_Pres_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpPresRZ(transip);
+   }
+};
+
+class G_EQDSK_Q_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_Q_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpQRZ(transip);
+   }
+};
+
+class G_EQDSK_BTor_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_BTor_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpBTorRZ(transip);
+   }
+};
+
+class G_EQDSK_JTor_Coefficient : public Coefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_JTor_Coefficient(G_EQDSK_Data &g_eqdsk) : eqdsk(g_eqdsk) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return eqdsk.InterpJTorRZ(transip);
+   }
+};
+
+class G_EQDSK_NxGradPsi_Coefficient : public VectorCoefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_NxGradPsi_Coefficient(G_EQDSK_Data &g_eqdsk)
+      : VectorCoefficient(2), eqdsk(g_eqdsk) {}
+
+   void Eval(Vector &b, ElementTransformation & T,
+             const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      eqdsk.InterpNxGradPsiRZ(transip, b);
+   }
+
+};
+
+class G_EQDSK_BPol_Coefficient : public VectorCoefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+
+public:
+
+   G_EQDSK_BPol_Coefficient(G_EQDSK_Data &g_eqdsk)
+      : VectorCoefficient(2), eqdsk(g_eqdsk) {}
+
+   void Eval(Vector &b, ElementTransformation & T,
+             const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      eqdsk.InterpBPolRZ(transip, b);
+   }
+
+};
+
+class G_EQDSK_BField_VecCoefficient : public VectorCoefficient
+{
+private:
+   G_EQDSK_Data &eqdsk;
+   bool unit_;
+
+public:
+
+   G_EQDSK_BField_VecCoefficient(G_EQDSK_Data &g_eqdsk, bool unit)
+      : VectorCoefficient(3), eqdsk(g_eqdsk), unit_(unit) {}
+
+   void Eval(Vector &V, ElementTransformation & T,
+             const IntegrationPoint & ip)
+   {
+      V.SetSize(3);
+      Vector b;
+      b.SetSize(2);
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      eqdsk.InterpBPolRZ(transip, b);
+      double btor = eqdsk.InterpBTorRZ(transip);
+
+      V[0] = b[0];
+      V[1] = b[1];
+      V[2] = btor;
+
+      if ( unit_ )
+      {
+         double bmag = sqrt(V * V);
+         V /= bmag;
+      }
+   }
+};
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_G_EQDSK_DATA_HPP

miniapps/plasma/interp_data.cpp

@@ -0,0 +1,300 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "interp_data.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+Interp_Data::Interp_Data(istream &is)
+   : init_flag_(0)
+{
+   double XDUM = 0.0;
+
+   const int buflen = 1024;
+   char buf[buflen];
+   is.getline(buf, buflen);
+   istringstream iss(buf);
+   string word;
+   iss >> std::ws;
+   while (!iss.eof())
+   {
+      iss >> word;
+      CASE_.push_back(word);
+      iss >> std::ws;
+   }
+
+   NW_ = to_int(CASE_[CASE_.size()-2]);
+   NH_ = to_int(CASE_[CASE_.size()-1]);
+
+   is >> RDIM_ >> ZDIM_ >> RLEFT_ >> ZMID_;
+
+   FIELD_.resize(NW_ * NH_);
+
+   for (int j=0; j<NH_; j++)
+   {
+      for (int i=0; i<NW_; i++)
+      {
+         is >> FIELD_[NH_ * i + j];
+      }
+   }
+
+   dr_ = RDIM_ / (NW_ - 1);
+   dz_ = ZDIM_ / (NH_ - 1);
+}
+
+void Interp_Data::PrintInfo(ostream & out) const
+{
+   out << endl << "Outside Plasma Field Info:" << endl;
+   out << "Size of grid: " << NW_ << " x " << NH_ << endl;
+   out << "Range of R: " << RLEFT_ << " -> " << RLEFT_ + RDIM_ << endl;
+   out << "Range of Z: " << ZMID_ - 0.5 * ZDIM_
+       << " -> " << ZMID_ + 0.5 * ZDIM_ << endl;
+}
+
+double Interp_Data::InterpDataRZ(const Vector &rz)
+{
+   initInterpRZ(FIELD_, DATA_c_, DATA_d_, DATA_e_);
+   return interpRZ(rz, FIELD_, DATA_c_, DATA_d_, DATA_e_);
+}
+
+void Interp_Data::initInterpRZ(const std::vector<double> &v,
+                                ShiftedDenseMatrix &c,
+                                ShiftedDenseMatrix &d,
+                                ShiftedDenseMatrix &e)
+{
+   ExtendedDenseMatrix ve(&v[0], NW_, NH_);
+
+   c.SetSize(NW_ + 3, NH_ + 2); c.SetShifts(2, 1); c = 0.0;
+   d.SetSize(NW_ + 2, NH_ + 3); d.SetShifts(1, 2); d = 0.0;
+   e.SetSize(NW_ + 1, NH_ + 1); e.SetShifts(1, 1); e = 0.0;
+
+   // x-directed divided differences
+   for (int i=-1; i<NW_; i++)
+   {
+      c(i,-1) = (ve(i+1,-1) - ve(i,-1)) / dr_;
+   }
+   for (int j=0; j<NH_; j++)
+   {
+      for (int i=-2; i<=NW_; i++)
+      {
+         c(i,j) = (ve(i+1,j) - ve(i,j)) / dr_;
+      }
+   }
+   for (int i=-1; i<NW_; i++)
+   {
+      c(i,NH_) = (ve(i+1,NH_) - ve(i,NH_)) / dr_;
+   }
+
+   // y-directed divided differences
+   for (int j=-1; j<NH_; j++)
+   {
+      d(-1,j) = (ve(-1,j+1) - ve(-1,j)) / dz_;
+   }
+   for (int i=0; i<NW_; i++)
+   {
+      for (int j=-2; j<=NH_; j++)
+      {
+         d(i,j) = (ve(i,j+1) - ve(i,j)) / dz_;
+      }
+   }
+   for (int j=-1; j<NH_; j++)
+   {
+      d(NW_,j) = (ve(NW_,j+1) - ve(NW_,j)) / dz_;
+   }
+
+   // Second order divided differences
+   for (int i=-1; i<NW_; i++)
+   {
+      for (int j=-1; j<NH_; j++)
+      {
+         e(i,j) = (c(i,j+1) - c(i,j)) / dz_;
+      }
+   }
+}
+
+double Interp_Data::interpRZ(const Vector &rz,
+                              const std::vector<double> &v,
+                              const ShiftedDenseMatrix &c,
+                              const ShiftedDenseMatrix &d,
+                              const ShiftedDenseMatrix &e)
+{
+   double r = rz[0];
+   double z = rz[1];
+
+   double rs = (r - RLEFT_) / RDIM_;
+   double zs = (z - ZMID_ + 0.5 * ZDIM_) / ZDIM_;
+
+   int i = std::max(0, std::min((int)floor(double(NW_-1) * rs), NW_-2));
+   int j = std::max(0, std::min((int)floor(double(NH_-1) * zs), NH_-2));
+
+   // Compute corners of local patch
+   double r0 = RLEFT_ + RDIM_ * i / (NW_ - 1);
+   double r1 = r0 + RDIM_ / (NW_ - 1);
+   double z0 = ZMID_ - 0.5 * ZDIM_ + ZDIM_ * j / (NH_ - 1);
+   double z1 = z0 + ZDIM_ / (NH_ - 1);
+
+   // Prepare position dependent factors
+   double wra = (r1 - r) / dr_;
+   double wrb = (r - r0) / dr_;
+   double wrc = (1.0 + 2.0 * wra);
+   double wrd = (1.0 + 2.0 * wrb);
+   double wra2 = wra * wra;
+   double wrb2 = wrb * wrb;
+
+   double wza = (z1 - z) / dz_;
+   double wzb = (z - z0) / dz_;
+   double wzc = (1.0 + 2.0 * wza);
+   double wzd = (1.0 + 2.0 * wzb);
+   double wza2 = wza * wza;
+   double wzb2 = wzb * wzb;
+
+   // Extract variable values at corners of local patch
+   double p00 = v[NH_ * i + j];
+   double p10 = v[NH_ * (i + 1) + j];
+   double p01 = v[NH_ * i + j + 1];
+   double p11 = v[NH_ * (i + 1) + j + 1];
+
+   double var = p00 * wra2 * wrd * wza2 * wzd
+                + p10 * wrb2 * wrc * wza2 * wzd
+                + p01 * wra2 * wrd * wzb2 * wzc
+                + p11 * wrb2 * wrc * wzb2 * wzc;
+
+   // Compute dvar/dx at corners of local patch
+   double wx00a = fabs(c(i-1,j) - c(i-2,j));
+   double wx00b = fabs(c(i+1,j) - c(i,j));
+
+   double wx10a = fabs(c(i,j) - c(i-1,j));
+   double wx10b = fabs(c(i+2,j) - c(i+1,j));
+
+   double wx01a = fabs(c(i-1,j+1) - c(i-2,j+1));
+   double wx01b = fabs(c(i+1,j+1) - c(i,j+1));
+
+   double wx11a = fabs(c(i,j+1) - c(i-1,j+1));
+   double wx11b = fabs(c(i+2,j+1) - c(i+1,j+1));
+
+   if (wx00a == 0.0 && wx00b == 0.0) { wx00a = 1.0; wx00b = 1.0; }
+   if (wx10a == 0.0 && wx10b == 0.0) { wx10a = 1.0; wx10b = 1.0; }
+   if (wx01a == 0.0 && wx01b == 0.0) { wx01a = 1.0; wx01b = 1.0; }
+   if (wx11a == 0.0 && wx11b == 0.0) { wx11a = 1.0; wx11b = 1.0; }
+
+   double px00 = (wx00b * c(i-1,j) + wx00a * c(i,j)) / (wx00b + wx00a);
+   double px10 = (wx10b * c(i,j) + wx10a * c(i+1,j)) / (wx10b + wx10a);
+   double px01 = (wx01b * c(i-1,j+1) + wx01a * c(i,j+1)) / (wx01b + wx01a);
+   double px11 = (wx11b * c(i,j+1) + wx11a * c(i+1,j+1)) / (wx11b + wx11a);
+
+   double varx = px00 * wra2 * wrb * wza2 * wzd
+                 - px10 * wrb2 * wra * wza2 * wzd
+                 + px01 * wrb * wra2 * wzb2 * wzc
+                 - px11 * wra * wrb2 * wzb2 * wzc;
+   var += varx * dr_;
+
+   // Compute dvar/dy at corners of local patch
+   double wy00a = fabs(d(i,j-1) - d(i,j-2));
+   double wy00b = fabs(d(i,j+1) - d(i,j));
+
+   double wy10a = fabs(d(i+1,j-1) - d(i+1,j-2));
+   double wy10b = fabs(d(i+1,j+1) - d(i+1,j));
+
+   double wy01a = fabs(d(i,j) - d(i,j-1));
+   double wy01b = fabs(d(i,j+2) - d(i,j+1));
+
+   double wy11a = fabs(d(i+1,j) - d(i+1,j-1));
+   double wy11b = fabs(d(i+1,j+2) - d(i+1,j+1));
+
+   if (wy00a == 0.0 && wy00b == 0.0) { wy00a = 1.0; wy00b = 1.0; }
+   if (wy10a == 0.0 && wy10b == 0.0) { wy10a = 1.0; wy10b = 1.0; }
+   if (wy01a == 0.0 && wy01b == 0.0) { wy01a = 1.0; wy01b = 1.0; }
+   if (wy11a == 0.0 && wy11b == 0.0) { wy11a = 1.0; wy11b = 1.0; }
+
+   double py00 = (wy00b * d(i,j-1) + wy00a * d(i,j)) / (wy00b + wy00a);
+   double py10 = (wy10b * d(i+1,j-1) + wy10a * d(i+1,j)) / (wy10b + wy10a);
+   double py01 = (wy01b * d(i,j) + wy01a * d(i,j+1)) / (wy01b + wy01a);
+   double py11 = (wy11b * d(i+1,j) + wy11a * d(i+1,j)) / (wy11b + wy11a);
+
+   double vary = py00 * wra2 * wrd * wza2 * wzb
+                 + py10 * wrb2 * wrc * wza2 * wzb
+                 - py01 * wra2 * wrd * wza * wzb2
+                 - py11 * wrb2 * wrc * wza * wzb2;
+   var += vary * dz_;
+
+   // Compute d^2var/dxdy at corners of local patch
+   double pxy00 = (wx00b * (wy00b * e(i-1,j-1) + wy00a * e(i-1,j)) +
+                   wx00a * (wy00b * e(i,j-1) + wy00a * e(i,j))) /
+                  ((wx00b + wx00a) * (wy00b + wy00a));
+   double pxy10 = (wx10b * (wy10b * e(i,j-1) + wy10a * e(i,j)) +
+                   wx10a * (wy10b * e(i+1,j-1) + wy10a * e(i+1,j))) /
+                  ((wx10b + wx10a) * (wy10b + wy10a));
+   double pxy01 = (wx01b * (wy01b * e(i-1,j) + wy01a * e(i-1,j+1)) +
+                   wx01a * (wy01b * e(i,j) + wy01a * e(i,j+1))) /
+                  ((wx01b + wx01a) * (wy01b + wy01a));
+   double pxy11 = (wx11b * (wy11b * e(i,j) + wy11a * e(i,j+1)) +
+                   wx11a * (wy11b * e(i+1,j) + wy11a * e(i+1,j+1))) /
+                  ((wx11b + wx11a) * (wy11b + wy11a));
+
+   double varxy = pxy00 * wra2 * wrb * wza2 * wzb
+                  - pxy10 * wra * wrb2 * wza2 * wzb
+                  - pxy01 * wra2 * wrb * wza * wzb2
+                  + pxy11 * wra * wrb2 * wza * wzb2;
+
+   var += dr_ * dz_ * varxy;
+
+   return var;
+}
+
+void Interp_Data::ExtendedDenseMatrix::init()
+{
+   // Populate four corners
+   SW_ = 3.0 * ((*this)(0,0) - (*this)(1,1)) + (*this)(2,2);
+   SE_ = 3.0 * ((*this)(m_-1,0) - (*this)(m_-2,1)) + (*this)(m_-3,2);
+   NW_ = 3.0 * ((*this)(0,n_-1) - (*this)(1,n_-2)) + (*this)(2,n_-3);
+   NE_ = 3.0 * ((*this)(m_-1,n_-1) - (*this)(m_-2,n_-2))
+         + (*this)(m_-3,n_-3);
+
+   // Populate lowest rows
+   for (int j=0; j<n_; j++)
+   {
+      S_(1,j) = 3.0 * ((*this)(0,j) - (*this)(1,j)) + (*this)(2,j);
+      S_(0,j) = 3.0 * (2.0 * (*this)(0,j) + (*this)(2,j)) - 8.0 * (*this)(1,j);
+   }
+
+   // Populate highest rows
+   for (int j=0; j<n_; j++)
+   {
+      N_(1,j) = 3.0 * (2.0 * (*this)(m_-1,j) + (*this)(m_-3,j))
+                - 8.0 * (*this)(m_-2,j);
+      N_(0,j) = 3.0 * ((*this)(m_-1,j) - (*this)(m_-2,j)) + (*this)(m_-3,j);
+   }
+
+   // Populate lowest columns
+   for (int i=0; i<m_; i++)
+   {
+      W_(i,0) = 3.0 * (2.0 * (*this)(i,0) + (*this)(i,2)) - 8.0 * (*this)(i,1);
+      W_(i,1) = 3.0 * ((*this)(i,0) - (*this)(i,1)) + (*this)(i,2);
+   }
+
+   // Populate highest columns
+   for (int i=0; i<m_; i++)
+   {
+      E_(i,0) = 3.0 * ((*this)(i,n_-1) - (*this)(i,n_-2)) + (*this)(i,n_-3);
+      E_(i,1) = 3.0 * (2.0 * (*this)(i,n_-1) + (*this)(i,n_-3))
+                - 8.0 * (*this)(i,n_-2);
+   }
+}
+
+} // namespace plasma
+
+} // namespace mfem

miniapps/plasma/interp_data.hpp

@@ -0,0 +1,226 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_INTERP_DATA_HPP
+#define MFEM_INTERP_DATA_HPP
+
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+
+#include "mfem.hpp"
+#include "../../general/text.hpp"
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+class Interp_Data
+{
+public:
+   Interp_Data(std::istream &is);
+
+   int GetNumPtsR() const { return NW_; }
+   int GetNumPtsZ() const { return NH_; }
+
+   double GetRExtent() const { return RDIM_; }
+   double GetZExtent() const { return ZDIM_; }
+
+   double GetRMin() const { return RLEFT_; }
+   double GetZMid() const { return ZMID_; }
+
+   void PrintInfo(std::ostream &out = std::cout) const;
+
+   double InterpDataRZ(const Vector &rz);
+
+private:
+   class ShiftedVector;
+   class ShiftedDenseMatrix;
+   class ExtendedDenseMatrix;
+
+   int init_flag_;
+
+   void initInterpRZ(const std::vector<double> &v,
+                     ShiftedDenseMatrix &c,
+                     ShiftedDenseMatrix &d,
+                     ShiftedDenseMatrix &e);
+
+   double interpRZ(const Vector &rz,
+                   const std::vector<double> &v,
+                   const ShiftedDenseMatrix &c,
+                   const ShiftedDenseMatrix &d,
+                   const ShiftedDenseMatrix &e);
+
+   std::vector<std::string> CASE_; // Identification character string
+
+   int NW_; // Number of horizontal R grid points
+   int NH_; // Number of vertical Z grid points
+
+   double RDIM_;    // Horizontal dimension in meter of computational box
+   double ZDIM_;    // Vertical dimension in meter of computational box
+   double RLEFT_;   // Minimum R in meter of rectangular computational box
+   double ZMID_;    // Z of center of computational box in meter
+
+   // Field on grid
+   std::vector<double> FIELD_;
+
+   class ShiftedVector : public Vector
+   {
+   private:
+      int si_;
+   public:
+      ShiftedVector()
+         : si_(0) {}
+
+      ShiftedVector(int s, int si)
+         : Vector(s+2*si), si_(si) {}
+
+      void SetShift(int si) { si_ = si; }
+
+      ShiftedVector &operator=(double c)
+      { Vector::operator=(c); return *this; }
+
+      inline double &operator()(int i)
+      { return Vector::operator()(i + si_); }
+
+      inline const double &operator()(int i) const
+      { return Vector::operator()(i + si_); }
+   };
+
+   class ShiftedDenseMatrix : public DenseMatrix
+   {
+   private:
+      int si_, sj_;
+   public:
+      ShiftedDenseMatrix()
+         : si_(0), sj_(0) {}
+
+      ShiftedDenseMatrix(int m, int n, int si, int sj)
+         : DenseMatrix(m+2*si, n+2*sj), si_(si), sj_(sj) {}
+
+      void SetShifts(int si, int sj) { si_ = si; sj_ = sj; }
+
+      ShiftedDenseMatrix &operator=(double c)
+      { DenseMatrix::operator=(c); return *this; }
+
+      inline double &operator()(int i, int j)
+      { return DenseMatrix::operator()(i + si_, j + sj_); }
+
+      inline const double &operator()(int i, int j) const
+      { return DenseMatrix::operator()(i + si_, j + sj_); }
+   };
+
+   class ExtendedDenseMatrix
+   {
+   private:
+      int m_, n_;
+      const double *C_;
+      DenseMatrix N_;
+      DenseMatrix S_;
+      DenseMatrix E_;
+      DenseMatrix W_;
+      double SW_, SE_, NW_, NE_, DUMMY_;
+
+      void init();
+
+   public:
+      ExtendedDenseMatrix(const double *C, int m, int n)
+         : m_(m), n_(n), C_(C),
+           N_(2, n), S_(2, n),
+           E_(m, 2), W_(m, 2),
+           SW_(0.0), SE_(0.0), NW_(0.0), NE_(0.0), DUMMY_(0.0)
+      { N_ = 0.0; S_ = 0.0; E_ = 0.0; W_ = 0.0; init(); }
+
+      const double &operator()(int i, int j) const
+      {
+         if (i >= 0 && i < m_ && j >= 0 && j < n_)
+         {
+            return C_[n_ * i + j];
+         }
+         else if (i >= 0 && i < m_)
+         {
+            if (j < 0)
+            {
+               return W_(i, j + 2);
+            }
+            else
+            {
+               return E_(i, j - n_);
+            }
+         }
+         else if (j >= 0 && j < n_)
+         {
+            if (i < 0)
+            {
+               return S_(i + 2, j);
+            }
+            else
+            {
+               return N_(i - m_, j);
+            }
+         }
+         else if (i == -1 && j == -1)
+         {
+            return SW_;
+         }
+         else if (i == -1 && j == n_)
+         {
+            return SE_;
+         }
+         else if (i == m_ && j == -1)
+         {
+            return NW_;
+         }
+         else if (i == m_ && j == n_)
+         {
+            return NE_;
+         }
+         return DUMMY_;
+      }
+   };
+
+   // Divided differences for Akima's interpolation method
+   double dr_, dz_;
+   ShiftedDenseMatrix  DATA_c_;
+   ShiftedDenseMatrix  DATA_d_;
+   ShiftedDenseMatrix  DATA_e_;
+};
+
+class Interp_Data_Coefficient : public Coefficient
+{
+private:
+   Interp_Data &interp_data;
+
+public:
+
+   Interp_Data_Coefficient(Interp_Data &i_data) : interp_data(i_data) {}
+
+   double Eval(ElementTransformation & T,
+               const IntegrationPoint & ip)
+   {
+      double x[3];
+      Vector transip(x, 3);
+
+      T.Transform(ip, transip);
+
+      return interp_data.InterpDataRZ(transip);
+   }
+};
+
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_INTERP_DATA_HPP

miniapps/plasma/makefile

@@ -0,0 +1,84 @@
+# 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/plasma/,)
+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 = stix1d stix2d stix3d stix1d_dh stix2d_dh
+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
+.PRECIOUS: %.o
+
+COMMON_O=../common/fem_extras.o ../common/pfem_extras.o \
+	../common/mesh_extras.o \
+	cold_plasma_dielectric_solver.o cold_plasma_dielectric_coefs.o \
+	cold_plasma_dielectric_dh_solver.o g_eqdsk_data.o interp_data.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(MINIAPPS)
+
+# Rules for building the miniapps
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -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: Specific execution options
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+stix1d-test-par:
+	@true
+stix2d-test-par:
+	@true
+stix3d-test-par:
+	@true
+
+# 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 *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf STIX1D-AMR-Parallel* STIX2D-AMR-Parallel* STIX3D-AMR-Parallel*
+	@rm -rf STIX1D-DH-AMR-Parallel* STIX2D-DH-AMR-Parallel*
+

miniapps/plasma/plasma.hpp

@@ -0,0 +1,63 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PLASMA_HPP
+#define MFEM_PLASMA_HPP
+
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "mfem.hpp"
+#include "../../general/text.hpp"
+
+namespace mfem
+{
+
+namespace plasma
+{
+
+// 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;
+
+// Speed of light in Free Space (units m/s)
+static const double c0_ = 1.0 / sqrt(epsilon0_ * mu0_);
+
+static const double q_     = 1.602176634e-19; // Elementary charge in coulombs
+static const double eV_    = 1.602176634e-19; // 1 eV in Joules
+static const double amu_   = 1.660539040e-27; // Atomic mass unit in kilograms
+static const double me_kg_ = 9.10938356e-31;  // Mass of electron in kilograms
+static const double me_u_  = 5.4857990907e-4; // Mass of electron in a.m.u
+
+/**
+   Returns the cyclotron frequency in radians/second
+   m is the mass in a.m.u
+   q is the charge in units of elementary electric charge
+   B is the magnetic field magnitude in tesla
+ */
+inline double cyclotronFrequency(double B, double m, double q)
+{
+   return fabs(q * q_ * B / (m * amu_));
+}
+
+} // namespace plasma
+
+} // namespace mfem
+
+#endif // MFEM_PLASMA_HPP

miniapps/plasma/simple_antenna.cpp

@@ -0,0 +1,332 @@
+#include "mfem.hpp"
+#include "../common/mesh_extras.hpp"
+
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::common;
+
+static double s[] =
+{
+   1.0,
+   0.6180339887498949,
+   0.5436890126920764,
+   0.5187900636758842,
+   0.5086603916420042,
+   0.5041382583616554,
+   0.5020170551781655,
+   0.5009941779228898,
+   0.5004931182865523,
+   0.5002454622667946,
+   0.5001224294760432,
+   0.5000611322390582,
+   0.5000305436878334,
+   0.500015265778675,
+   0.5000076312578446,
+   0.5000038151921251
+};
+
+int main(int argc, char ** argv)
+{
+   int mfb, mf, mb, na, nb, nt;
+   double af, ab, ba, bb, bt;
+   bool per_y = false;
+   bool visualization = true;
+
+   mf = mb = na = nb = nt = -1;
+   af = ab = ba = bb = bt = -1.0;
+   mfb = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mf, "-mf", "--num-front",
+                  "Number of elements in front of antenna (>= 1).");
+   args.AddOption(&mb, "-mb", "--num-back",
+                  "Number of elements behind antenna (>= 1).");
+   args.AddOption(&nb, "-nb", "--num-bottom",
+                  "Number of elements below antenna (>= 1).");
+   args.AddOption(&nt, "-nt", "--num-top",
+                  "Number of elements above antenna (>= 1).");
+   args.AddOption(&na, "-na", "--num-across",
+                  "Number of elements across antenna (>= 2).");
+   args.AddOption(&af, "-af", "--size-front",
+                  "Distance in front of antenna (> 0).");
+   args.AddOption(&ab, "-ab", "--size-back",
+                  "Distance behind antenna (> 0).");
+   args.AddOption(&bb, "-bb", "--size-bottom",
+                  "Distance below antenna (> 0).");
+   args.AddOption(&bt, "-bt", "--size-top",
+                  "Distance above antenna (> 0).");
+   args.AddOption(&ba, "-ba", "--size-across",
+                  "Distance across antenna (> 0).");
+   args.AddOption(&mfb, "-mfb", "--num-bdr-front",
+                  "Number of elements in boundry layer "
+                  "in front of antenna (>= 1).");
+   args.AddOption(&per_y, "-per-y", "--periodic-y", "-no-per-y",
+                  "--no-periodic-y",
+                  "Make the mesh periodic in the y direction.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+
+   if (mf < 0) { mf = 1; }
+   if (mb < 0) { mb = 1; }
+   if (na < 0) { na = 2; }
+   if (nb < 0) { nb = 1; }
+   if (nt < 0) { nt = 1; }
+   if (mfb < 1) { mfb = 1; }
+
+   if (af < 0) { af = 0.75; }
+   if (ab < 0) { ab = 0.25; }
+   if (ba < 0) { ba = 0.5; }
+   if (bb < 0) { bb = 0.25; }
+   if (bt < 0) { bt = 0.25; }
+
+   args.PrintOptions(cout);
+
+   MFEM_VERIFY(na >= 2,
+               "There must be at least two elements across "
+               "the face of the antenna");
+   MFEM_VERIFY(mf > 0 && na > 0 && nb > 0 && nt > 0 && mfb > 0,
+               "Numbers of elements must be greater than zero.");
+   MFEM_VERIFY(mfb <= 16, "Number of elements in boundary layer is too large.");
+   MFEM_VERIFY(af > 0.0 && ab > 0.0 && ba > 0.0 && bb > 0.0 && bt > 0.0,
+               "Distances must be greater than zero.");
+
+   int mx = mf + mb + mfb - 1;
+   int ny = nb + na + nt;
+
+   double ax = af + ab;
+   double by = bb + ba + bt;
+
+   int nelem = mx * ny;
+   int nnode = (mx + 1) * (ny + 1) + na - 1;
+   int nbdr  = 2 * mx + 2 * ny + 2 * na;
+
+   Mesh *mesh = new Mesh(2, nnode, nelem, nbdr);
+
+   // Create vertices
+   double c[2];
+   for (int j=0; j<=ny; j++)
+   {
+      double y0 = by * j / ny;
+      double ya = (j<=nb) ? (bb * j / nb) :
+                  ((j<=nb+na)? (bb + ba * (j - nb) / na) :
+                   (bb + ba + bt * (j - nb - na) / nt));
+
+      double dxf = af / mf;
+      double dxb = ab / mb;
+      double prev_cx = 0.0;
+
+      for (int i=0; i<=mx; i++)
+      {
+         if (i == 0)
+         {
+            c[0] = 0.0;
+            prev_cx = 0.0;
+         }
+         else if (mfb > 1 && i < mfb)
+         {
+            int p = mfb - i + 1;
+            double dc = dxf * pow(s[mfb-1], p);
+            c[0] = prev_cx + dc;
+            prev_cx = c[0];
+         }
+         else if (i <= mf + mfb - 1)
+         {
+            c[0] = dxf * (i - mfb + 1);
+         }
+         else
+         {
+            c[0] = af + dxb * (i - mf - mfb + 1);
+         }
+
+         if (i <= mf + mfb - 1)
+         {
+            c[1] = y0 + (ya - y0) * c[0] / af;
+         }
+         else
+         {
+            c[1] = y0 * (c[0] - af) / ab + ya * (ax - c[0]) / ab;
+         }
+
+         mesh->AddVertex(c);
+      }
+   }
+   for (int j=1; j < na; j++)
+   {
+      c[0] = (1.0 + 1.0e-4) * af;
+      c[1] = bb + ba * j / na;
+      mesh->AddVertex(c);
+   }
+
+   // Create elements
+   int v[4];
+   for (int j=0; j<nb; j++)
+   {
+      for (int i=0; i<mx; i++)
+      {
+         v[0] = j * (mx + 1) + i;
+         v[1] = j * (mx + 1) + i + 1;
+         v[2] = (j + 1) * (mx + 1) + i + 1;
+         v[3] = (j + 1) * (mx + 1) + i;
+
+         mesh->AddQuad(v);
+      }
+   }
+   for (int j=nb; j<nb + na; j++)
+   {
+      for (int i=0; i<mf+mfb-1; i++)
+      {
+         v[0] = j * (mx + 1) + i;
+         v[1] = j * (mx + 1) + i + 1;
+         v[2] = (j + 1) * (mx + 1) + i + 1;
+         v[3] = (j + 1) * (mx + 1) + i;
+
+         mesh->AddQuad(v);
+      }
+      for (int i=mf+mfb-1; i<mx; i++)
+      {
+         if (i == mf+mfb-1)
+         {
+            if (j == nb)
+            {
+               v[0] = j * (mx + 1) + i;
+               v[1] = j * (mx + 1) + i + 1;
+               v[2] = (j + 1) * (mx + 1) + i + 1;
+               v[3] = (mx + 1) * (ny + 1);
+               // v[3] = (j + 1) * (mx + 1) + i;
+            }
+            else if (j == nb + na -1)
+            {
+               v[0] = nnode - 1;
+               v[1] = j * (mx + 1) + i + 1;
+               v[2] = (j + 1) * (mx + 1) + i + 1;
+               v[3] = (j + 1) * (mx + 1) + i;
+            }
+            else
+            {
+               // v[0] = j * (mx + 1) + i;
+               v[0] = nnode - na + j - nb;
+               v[1] = j * (mx + 1) + i + 1;
+               v[2] = (j + 1) * (mx + 1) + i + 1;
+               v[3] = nnode - na + j - nb + 1;
+               // v[3] = (j + 1) * (mx + 1) + i;
+            }
+         }
+         else
+         {
+            v[0] = j * (mx + 1) + i;
+            v[1] = j * (mx + 1) + i + 1;
+            v[2] = (j + 1) * (mx + 1) + i + 1;
+            v[3] = (j + 1) * (mx + 1) + i;
+         }
+
+         mesh->AddQuad(v);
+      }
+   }
+   for (int j=nb + na; j<ny; j++)
+   {
+      for (int i=0; i<mx; i++)
+      {
+         v[0] = j * (mx + 1) + i;
+         v[1] = j * (mx + 1) + i + 1;
+         v[2] = (j + 1) * (mx + 1) + i + 1;
+         v[3] = (j + 1) * (mx + 1) + i;
+
+         mesh->AddQuad(v);
+      }
+   }
+
+   // Create boundary elements
+   for (int i=0; i<mx; i++)
+   {
+      v[0] = i;
+      v[1] = i + 1;
+      mesh->AddBdrSegment(v, 1);
+   }
+   for (int j=0; j<ny; j++)
+   {
+      v[0] = (mx + 1) * j + mx;
+      v[1] = (mx + 1) * (j + 1) + mx;
+      mesh->AddBdrSegment(v, 2);
+   }
+   for (int i=mx; i>0; i--)
+   {
+      v[0] = (mx + 1) * ny + i;
+      v[1] = (mx + 1) * ny + i - 1;
+      mesh->AddBdrSegment(v, 3);
+   }
+   for (int j=ny; j>0; j--)
+   {
+      v[0] = j * (mx + 1);
+      v[1] = (j - 1) * (mx + 1);
+      mesh->AddBdrSegment(v, 4);
+   }
+   for (int j=nb; j<na + nb; j++)
+   {
+      v[0] = (mx + 1) * j + mf + mfb - 1;
+      v[1] = (mx + 1) * (j + 1) + mf + mfb - 1;
+      mesh->AddBdrSegment(v, 5);
+   }
+   for (int j=nb+na; j>nb; j--)
+   {
+      if (j == nb + na)
+      {
+         v[0] = (mx + 1) * j + mf + mfb - 1;
+      }
+      else
+      {
+         v[0] = nnode - (nb + na - j);
+      }
+      if (j == nb + 1)
+      {
+         v[1] = (mx + 1) * (j - 1) + mf + mfb - 1;
+      }
+      else
+      {
+         v[1] = nnode - (nb + na - j + 1);
+      }
+      mesh->AddBdrSegment(v, 6);
+   }
+   mesh->FinalizeTopology();
+
+   if (per_y)
+   {
+      Array<int> v2v(mesh->GetNV());
+      for (int i=0; i<v2v.Size(); i++) { v2v[i] = i; }
+
+      for (int i=0; i<=mx; i++)
+      {
+         v2v[(mx + 1) * ny + i] = i;
+      }
+
+      Mesh * per_mesh = MakePeriodicMesh(mesh, v2v);
+      delete mesh;
+      mesh = per_mesh;
+   }
+
+   ofstream mesh_ofs("simple_antenna.mesh");
+   mesh_ofs.precision(8);
+   mesh->Print(mesh_ofs);
+
+   // Output the resulting mesh to GLVis
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "mesh\n" << *mesh << flush;
+   }
+
+   // Clean up and exit
+   delete mesh;
+}